dsi_heat_maturation.R

labelHeat <- function(rsvars,dsivars,Title,tifname) {
 # rsvars="Ramyg_cm2_16"
  #dsivars="Ramyg_cm2_16_final_qa0"
  
#location of completed analyses
rootpath <- "/Users/mariaj/Dropbox/Pitt/Luna/amyg_con/01_data/Z_scores_4_Fig3"
#where I want to put my images
imgpath <- "/Users/mariaj/Dropbox/Pitt/Luna/amyg_con/01_data/Z_scores_4_Fig3"

#make plot in variable that show the age effects identified in spline analyses
#rsvars <- c("fornix") 

#the number of the directory from spline analysis you want to pull (2 is for age )
m <- 2
#1.5 negative and positive SD for a normal distribution
percentiles <- 100 * pnorm(c(-2, 2))
#rsvars="Ramyg_cm2_16"
path <- file.path(rootpath,m)
#this is only if we want to load in the actual data into the plot (20150812-not currently in there)
rsydata<-read.delim(file.path(path,rsvars, "data.txt"),sep=" ")
dsiydata<-read.delim(file.path(path,dsivars, "data.txt"),sep=" ")


#for each variable (rsvars defined above) read in the stages_pred_data.txt and combine into one big data frame
rsdataFit <- ldply(rsvars, function(i) data.frame(var = i, read.table(file.path(path, i, "stages_pred_data.txt"), head = T)))
dsidataFit <- ldply(dsivars, function(i) data.frame(var = i, read.table(file.path(path, i, "stages_pred_data.txt"), head = T)))
#create the 1.5 sd lines (the light grey ones) for each variable
rsdataPerc <- ldply(rsvars, function(v) ldply(percentiles, function(i) {
  data.frame(var = v, perc = i, age = rsdataFit$age[which(rsdataFit$var==v)], 
             pred = with(rsdataFit[which(rsdataFit$var==v),], pred + qnorm(i/100) * sim.sd))}))

dsidataPerc <- ldply(dsivars, function(v) ldply(percentiles, function(i) {
  data.frame(var = v, perc = i, age = dsidataFit$age[which(dsidataFit$var==v)], 
             pred = with(dsidataFit[which(dsidataFit$var==v),], pred + qnorm(i/100) * sim.sd))}))

#for each variable, make a plot of your fit line (rsdataFit)
#also plot your 1.5 sd away from the rsdataFit (rsdataPerc)
#facet.grid is a ggplot command that allows you to lay out your plots on multiple panels on 1 pg
#need empty geom_path so that it plots your rsdataFit line
xlab <- " "
ylab <- "Z-Score"
ymin<-rsdataPerc[rsdataPerc$perc=="2.27501319481792",]
ymax<-rsdataPerc[rsdataPerc$perc=="97.7249868051821",]
        
pFit <- ggplot(data=rsdataFit, aes(x =age, y =pred)) +
        geom_line(colour="blue",size=2) +
        geom_ribbon(aes(ymin=ymin$pred,ymax=ymax$pred),alpha=0.2,fill="blue")+
        labs(list(x = xlab, y = ylab)) +
        theme_bw()+
        ggtitle(Title)

dsiymin<-dsidataPerc[dsidataPerc$perc=="2.27501319481792",]
dsiymax<-dsidataPerc[dsidataPerc$perc=="97.7249868051821",]
pFit2<-
  pFit+ 
  geom_path(data=dsidataFit, aes(x =age, y =pred),colour="dark cyan",size=2)+
  geom_ribbon(aes(ymin=dsiymin$pred,ymax=dsiymax$pred),alpha=0.2,fill="dark cyan")+
  theme_bw()+
  theme(

    # panel.border = element_blank(), 
    # panel.grid.major = element_blank(),
    # panel.grid.minor = element_blank(),
   
    text=element_text(size=20),
    axis.title=element_text(vjust=0,size=20),
    axis.line = element_line(colour = "black"),
   
    axis.title.x=element_blank(),
    axis.text.x=element_blank()
  )

# for each variable read in the stages_predDeriv_data.txt
#tells you how many SD it is away from??
dataHeat <- ldply(rsvars, function(i) data.frame(var = i, read.table(file.path(path, i, "stages_predDeriv_data.txt"), head = T)))
dataHeat$z <- with(dataHeat, pred / sim.sd)
#dataHeat$z[which(dataHeat$p > 0.05)] <- 0
levels(dataHeat$var)="      "   
pHeat <- ggplot(dataHeat, aes(x = age, y = var, fill = z)) + 
         geom_raster() + 
         scale_fill_gradient2(
           " ", 
           low = "blue", 
           high = "red",
           limits=c(-4,4)
         ) + 
         labs(list(x = xlab, y=" \n")) + 
         theme_bw() +
         theme(
          panel.grid.minor = element_blank(),
          panel.grid.major = element_blank(),
          legend.position = c(-.057, .5),
          axis.ticks.y=element_line(colour="white"),
          legend.text=element_text(size=6),
          legend.title=element_text(size=8),
          axis.text=element_text(size=20),
          axis.title.x=element_blank(),
          axis.text.x=element_blank()
         )



dsidataHeat <- ldply(dsivars, function(i) data.frame(var = i, read.table(file.path(path, i, "stages_predDeriv_data.txt"), head = T)))
dsidataHeat$z <- with(dsidataHeat, pred / sim.sd)
#dsidataHeat$z[which(dsidataHeat$p > 0.05)] <- 0
levels(dsidataHeat$var)="      "   
 
pHeatdsi <- ggplot(dsidataHeat, aes(x = age, y = var, fill = z)) + 
  geom_raster() + 
  scale_fill_gradient2(
           "", 
           low = "dark cyan", 
           high = "red",
           limits=c(-4,4)
  ) + 
  labs(list(x = xlab, y=" \n")) + 
  theme_bw() +
  theme(panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        legend.position = c(-.057, .5),
        axis.ticks.y=element_line(colour="white"),
        legend.text=element_text(size=6),
        legend.title=element_text(size=8),
        axis.text=element_text(size=20))

 
  test<-paste("/Users/mariaj/Dropbox/Pitt/Luna/amyg_con/01_data/Z_scores_4_Fig3/2/Final_Fig3/",tifname,sep="") 
              pdf(test,width = 8, height = 8)
  grid.arrange( pFit2, pHeat, pHeatdsi, nrow=3,heights = c(1,.31,.4))
  dev.off()


}