# Don't delete this chunk if you are using the mosaic package # This loads the mosaic and dplyr packages require(mosaic)
# Some customization. You can alter or delete as desired (if you know what you are doing). # This changes the default colors in lattice plots. trellis.par.set(theme=theme.mosaic()) # knitr settings to control how R chunks work. require(knitr) opts_chunk$set( tidy=FALSE, # display code as typed size="small" # slightly smaller font for code )
# Load additional packages here. Uncomment the line below to use Project MOSAIC data sets. # require(mosaicData)
library(devtools) install_github("jtlovell/physGenomicsPVFinal") library(physGenomicsPVFinal)
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, 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, 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")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.