Supporting files:
Objective: To review transcriptional responses to influenza virus infection using mRNA-seq. The looked at lineage-/MHCII+/CD11b-/CD11c+ cDC sorted from dLN of WT and Ifnlr-/- mice on day 4 following IAV infection, corresponding to the time point of DC accumulation from the lung to the dLN to prime the CD8+ T cell response.
## Loading required package: limma
We will use the limma bioconductor package to perform linear modeling.
source("http://bioconductor.org/biocLite.R")
## Bioconductor version 3.2 (BiocInstaller 1.20.3), ?biocLite for help
## A new version of Bioconductor is available after installing the most
## recent version of R; see http://bioconductor.org/install
# load limma package
library(limma)
library(stringi)
counts <- read.csv(file="//Users//greener//Downloads//counts.csv",row.names = 1, header=T)
targets <- read.csv(file="//Users//greener//Downloads//targets.csv", header=T)
dim(counts)
## [1] 24060 19
dim(targets)
## [1] 19 6
means <- rowMeans(counts)
filter <- rowSums(counts) >= 75
table(filter)
## filter
## FALSE TRUE
## 10804 13256
counts <- counts[filter,]
dim(counts)
## [1] 13256 19
Plot data with log2 tranformation.
boxplot(log2(counts+1), las=2)
Normalize with voom and create model matrix
f <- factor(targets$Treatment, levels = unique(targets$Treatment))
design <- model.matrix(~0 + f)
colnames(design) <- levels(f)
eset_voom <- voom(counts, design,plot=TRUE)
plotMDS(eset_voom, labels = targets$Treatment, cex = .8)
norm_matrix <- eset_voom$E
write.csv(norm_matrix, file="norm_matrix.csv")
Lets look at the distribution of counts after normalization
boxplot(eset_voom$E, las=2, names= colnames(eset_voom$E))
Build Linear model.
fit <- lmFit(eset_voom, design)
cont.matrix2 <- makeContrasts(IFNLR_IAV - IFNLR_Mock,
WT_IAV - WT_Mock,
IFNAR_IAV - IFNLR_IAV,levels=design)
fit2 <- contrasts.fit(fit,cont.matrix2)
fit2 <- eBayes(fit2)
results <- decideTests(fit2, lfc=(.58), method="separate", adjust.method="BH", p.value=0.05);
a <- vennCounts(results)
print(a)
## IFNLR_IAV - IFNLR_Mock WT_IAV - WT_Mock IFNAR_IAV - IFNLR_IAV Counts
## 1 0 0 0 11140
## 2 0 0 1 800
## 3 0 1 0 184
## 4 0 1 1 43
## 5 1 0 0 593
## 6 1 0 1 126
## 7 1 1 0 239
## 8 1 1 1 131
## attr(,"class")
## [1] "VennCounts"
vennDiagram(a)
vennDiagram(results,include=c("up","down"),counts.col=c("red","green"), cex = .6)
write.fit(fit2, file="emily_DE_updated.txt", digits=3, method="separate", adjust="BH");
Run MDS plots but different variables.
plotMDS(eset_voom$E, labels=targets$FileName)
plotMDS(eset_voom$E, labels=targets$Treatment)
MDS_Hulls(eset_voom$E,shape_var=targets$virus, color_var=targets$KO)
## Loading required package: ggplot2
## initial value 23.437676
## final value 23.422391
## converged
Run Hierarchical Clustering on the voom normalized logCPMs.
dat.dist<-dist(t(eset_voom$E))
plot(hclust(dat.dist), cex=.75)
dat.dist<-dist(t(eset_voom$E))
plot(hclust(dat.dist), labels = targets$Treatment, cex=.75)
Print Session Info
sessionInfo()
## R version 3.2.0 (2015-04-16)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: OS X 10.10.5 (Yosemite)
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggplot2_2.2.1 stringi_1.1.1 BiocInstaller_1.20.3
## [4] plyr_1.8.4 limma_3.26.9 MASS_7.3-45
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.7 knitr_1.14 magrittr_1.5 munsell_0.4.3
## [5] colorspace_1.2-6 stringr_1.2.0 tools_3.2.0 grid_3.2.0
## [9] gtable_0.2.0 htmltools_0.3.5 assertthat_0.1 yaml_2.1.13
## [13] lazyeval_0.2.0 rprojroot_1.2 digest_0.6.10 tibble_1.2
## [17] formatR_1.4 evaluate_0.9 rmarkdown_1.6 labeling_0.3
## [21] scales_0.5.0 backports_1.1.0