scripts/de42012_plotting.R

In jtlovell/physGenomicsPVFinal: Differential Gene Expression for several Switchgrass Projects

library(physGenomicsPVFinal)
rm(list=ls())

load("/Users/John/Desktop/dropbox/Switchgrass_PlantPhys/stats_output/tempe2012_allstats.RData")
data(temple2012_6treatments)

### for conviencence, rename info and counts
info<-info12
counts<-counts12
library(ggplot2)
#################################
### Part 1: Statistical presentation of main model
# stats<-pipeLIMMA(counts=counts, info=info, block=info$Sub_Block, formula="~ Treatment")
# v<-stats$voom[["E"]]
# stats.fullmodel<-stats$simpleStats
# stats.allests<-stats$stats
#################################
ggplot(info[!is.na(info$PDWP),], aes(x=PDWP, y=MDWP, col=Treatment))+geom_point(size=4)+ theme_bw() +
  stat_smooth(method="lm", se=F, lty=2, alpha=.5, lwd=.5)+
  scale_color_manual(values=c("darkred","forestgreen","skyblue"))+
  scale_x_continuous("Pre-dawn Leaf Water Potential (MPa)")+
  scale_y_continuous("Mid Day Leaf Water Potential (MPa)")+
  ggtitle("Effect of sampling order on measured leaf water potential")

pqHists(stats.fullmodel, what.p="pvalue", what.q="qvalue", main="main effect treatment 2012", breaks=100)

#################################
### Part 2: Comparison of contrasts
# contrast.matrix<-makeContrasts(f25th-flow, fmean-flow, fambient-flow, f75th-flow, fhigh-flow,
#                                fmean-f25th, fambient-f25th, f75th-f25th, fhigh-f25th,
#                                fambient-fmean, f75th-fmean, fhigh-fmean,
#                                f75th-fambient, fhigh-fambient,
#                                fhigh-f75th,
#                                levels=design)
# lim.contrasts<-anovaLIMMA(counts=counts, design=design, block=info$Sub_Block, contrast.matrix=contrast.matrix)
#################################
sig.q.05<-makeBinarySig(lim.contrasts, what="q.value", alpha=0.05)
counts2Venn(x=sig.q.05, cols=c(6, 10, 13, 15), names=c("l.v.h","25th.v.h","h.v.mean","amb.v.h"), colors=c("darkblue","green","cyan","red"),
            main="comparison of genes affected by \n ~ Treatment Contrasts in 2012")

counts2Venn(x=sig.q.05, cols=c(1,6,15), names=c("F-stat","l.v.h","amb.v.high"), colors=c("darkblue","green","red"),
            main="comparison of  genes affected by \n ~ Contrasts and F-tests in 2012")
v.means<-voom2MeanHeatMaps(v=v[sig.q.05[,"q.value_Ftest"]==1,], grps=info$Treatment,rowids=info$ID,thresh=8)

#################################
# Part 3: Plotting of PCAs
#################################
library(ggplot2)
ggplot(pca, aes(x=PC1, y=PC2))+geom_point(aes(col=Treatment), size=4) + theme_bw() +
  scale_color_manual(values=c("darkred","darkorange","yellow","black","green","lightblue"))+
  scale_y_continuous("PCA #2, 15% variance explained") + scale_x_continuous("PCA #1, 34% variance explained") +  ggtitle("DESeq2 PCA vs. Leaf Water Potential")

ggplot(pca, aes(x=MDWP, y=PC1))+geom_point(aes(col=Treatment), size=4) + theme_bw() +
  scale_color_manual(values=c("darkred","darkorange","yellow","black","green","lightblue"))+
  scale_y_continuous("PCA #1, 34% variance explained") + scale_x_continuous("Mid Day Leaf Water Potential (MPa)") +  ggtitle("DESeq2 PCA vs. Leaf Water Potential")

ggplot(pca, aes(x=MDWP, y=PC2))+geom_point(aes(col=Treatment), size=4) + theme_bw() +
  scale_color_manual(values=c("darkred","darkorange","yellow","black","green","lightblue"))+
  scale_y_continuous("PCA #2, 15% variance explained") + scale_x_continuous("Mid Day Leaf Water Potential (MPa)") +  ggtitle("DESeq2 PCA vs. Leaf Water Potential")

#################################
# Part 4: Effect of midday water potential
# stats<-pipeLIMMA(counts=counts, info=info, block=info$Sub_Block, formula="~ MDWP")
# stats.fullmodel.mdwp<-stats$simpleStats
# stats.allests.mdwp<-stats$stats

#################################
mdwp.q.05<-makeBinarySig(stats.allests.mdwp, what="mdwp_q.value", alpha=0.05)
trt.q.05<-makeBinarySig(stats.fullmodel.subtrt, what="Fqvalue", alpha=0.05)
mdwp.q.1<-makeBinarySig(stats.allests.mdwp, what="mdwp_q.value", alpha=0.1)
trt.q.1<-makeBinarySig(stats.fullmodel.subtrt, what="Fqvalue", alpha=0.1)
hl.q.05<-makeBinarySig(lim.contrasts, what="q.value_fhigh...flow", alpha=0.05)
hl.q.1<-makeBinarySig(lim.contrasts, what="q.value_fhigh...flow", alpha=0.1)
trt.mdwp.qs<-data.frame(trt.05=trt.q.05,mdwp.05=mdwp.q.05, hl.05=hl.q.05,  trt.1=trt.q.1, mdwp.1=mdwp.q.1, hl.1=hl.q.1)
colSums(trt.mdwp.qs)
par(mfrow=c(1,1))
counts2Venn(x=trt.mdwp.qs, cols=c(1:3), names=c("trt","MDWP","hl"), colors=c("darkblue","red","cyan"),
            main="comparison of  genes affected by \n ~ Treatment vs. Midday Water Potential F-tests alpha = 0.05")
counts2Venn(x=trt.mdwp.qs, cols=c(4:6), names=c("trt","MDWP","hl"), colors=c("darkblue","red","cyan"),
            main="comparison of  genes affected by \n ~ Treatment vs. Midday Water Potential F-tests alpha = 0.1")

#################################
# Part 4: Effect of sampling order
# stats<-pipeLIMMA(counts=counts, info=info, block=info$Sub_Block, formula="~ order")
# stats.fullmodel.order<-stats$simpleStats
# stats.allests.order<-stats$stats
#################################
ggplot(info[!is.na(info$PDWP),], aes(x=order, y=MDWP, col=Treatment))+geom_point(size=4)+ theme_bw() +
  stat_smooth(method="lm", se=F, lty=2, alpha=.5, lwd=.5)+
  scale_color_manual(values=c("darkred","forestgreen","cornflowerblue"))+
  scale_x_continuous("Sampling order")+
  scale_y_continuous("Mid Day Leaf Water Potential (MPa)")+
  ggtitle("Effect of sampling order on measured leaf water potential")

sum(stats.allests.order$ebayes_order_q.value<=0.05)
order.q.05<-makeBinarySig(stats.allests.order, what="ebayes_order_q.value", alpha=0.05)

vorder<-v
stats.order<-stats.allests.order
maxord<-as.character(stats.order[order(stats.order$ebayes_order_q.value)[1:20],"gene"])
vorder.pos<-data.frame(ID=colnames(vorder), t(vorder[maxord,]))
vorder2<-merge(info, vorder.pos, by="ID")
vtp<-melt(vorder2, id.vars=colnames(vorder2)[grep("Pavir", colnames(vorder2), invert=T)])

ggplot(vtp[!is.na(vtp$order),], aes(x=order, y=value, col=Treatment))+ geom_point()+
  facet_wrap(~variable, scales="free_y", nrow=5, ncol=4)+
  scale_color_manual(values=c("darkred","forestgreen","cornflowerblue"))+
  stat_smooth(span = 200,se=F, lty=2, alpha=.2, lwd=.5)+
  theme_bw()+
  theme(strip.text.x = element_text(size = 8))+
  scale_y_continuous("variance stabilized / normalized counts")+
  scale_x_continuous("order of RNA extraction, first plants were measured at 11:00, last plants ~ 13:00")+
  ggtitle("Effect of sampling time on 20 of the genes with strongest effects")