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

# 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())
data(tmpwfc20132014_2treatments)

### for conviencence, rename info and counts
info<-info1314
counts<-counts1314

#################################
# Part 2.1: Full model
#################################
stats<-pipeLIMMA(counts=counts, info=info, block=info$sb_unique, formula="~ Treatment * Location + Year + order")
v<-stats$voom[["E"]]
stats.fullmodel<-stats$simpleStats
stats.allests<-stats$stats

#################################
# Part 1.2: All pairwise contrasts
#################################
## Generate a design matrix that specifies all relavant contrasts

design <- model.matrix(~ 0+factor(Year_Location_Treatment.mat) + order , data=info)
colnames(design)[1:8]<-c("TMP_13_high","TMP_13_low","WFC_13_high","WFC_13_low",
                         "TMP_14_high","TMP_14_low","WFC_14_high","WFC_14_low")

contrast.matrix <- makeContrasts(TMP_13_high-TMP_13_low, WFC_13_high-WFC_13_low, TMP_14_high-TMP_14_low, WFC_14_high-WFC_14_low,
                                 levels=design)
lim.contrasts<-anovaLIMMA(counts=counts, design=design, block=info$Sub_Block, contrast.matrix=contrast.matrix)

#################################
# Part 3: Run PCA
#################################
pca<-voom2PCA(v=v, info=info, ids=info$ID)
ggplot(pca, aes(x=PC1, y=PC2, col=Location, shape=Treatment))+ facet_wrap(~Year, nrow=2)
  theme_bw()+geom_point(size=4)+scale_shape_manual(values=c(2,19))

#################################
# Part 4: Run model with MDWP as the predictor
#################################
stats<-pipeLIMMA(counts=counts, info=info, block=info$sb_unique, formula="~ MDWP * Location + Year + order")
stats.fullmodel.mdwp<-stats$simpleStats
stats.allests.mdwp<-stats$stats

#################################
# Part 5: Run model with MDWP as the predictor, without controlling for location or year
#################################
stats<-pipeLIMMA(counts=counts, info=info, block=info$sb_unique, formula="~ MDWP + order")
stats.fullmodel.mdwponly<-stats$simpleStats
stats.allests.mdwponly<-stats$stats

jtlovell/physGenomicsPVFinal documentation built on May 20, 2019, 3:14 a.m.