--- title: "Emily Hemann IFN_IAV" author: "Richard Green" date: "Aug 27, 2018" --- 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. ```{r , echo=FALSE, include=TRUE} library(MASS) if (!require("limma")) { install.packages("limma", repos="http://cran.rstudio.com/") library("limma") } ``` ```{r MDS function, echo=FALSE} library(plyr) MDS_Hulls <- function(data, color_var, shape_var, is_shape_numeric = FALSE){ require(MASS); require(ggplot2); require(plyr) d <- dist(t(data)) d[d==0] = 0.0000000000000000000000000000001 fit <- isoMDS(d) fit # view results, stress must < 30 x <- fit$points[,1] y <- fit$points[,2] min_lim = min(min(x),min(y)) max_lim = max(max(x),max(y)) datai = data.frame(x =x, y =y) datai$id = names(x) datai$cols = as.factor(color_var) datai$shape = as.factor(shape_var) datai$hull = as.factor(paste0(color_var,shape_var)) find_hull <- function(df) df[chull(df$x, df$y), ] # Computes the subset of points which lie on the convex hull of the set of points specified. hulls <- ddply(datai, "hull", find_hull) p <- ggplot(datai, aes(x= x, y=y,colour=cols,shape=shape,label = id)) png <- p + geom_point(alpha=1,size = 3) + scale_shape_manual(values=(0:(length(levels(datai$shape))-1))) + xlim(min_lim, max_lim)+ylim (min_lim, max_lim) + #xlim(-100,100)+ylim(-100,100)+ geom_polygon(data = hulls,aes(group = hull, fill=cols), alpha = 0.4, size = 0.5) + annotate("text", x=min_lim, y=min_lim, label=paste("stress=",round(fit$stress,2)), size = 2)+ theme_bw() print(png) } ``` We will use the limma bioconductor package to perform linear modeling. ```{r, echo=TRUE} source("http://bioconductor.org/biocLite.R") # 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) dim(targets) means <- rowMeans(counts) filter <- rowSums(counts) >= 75 table(filter) counts <- counts[filter,] dim(counts) ``` Plot data with log2 tranformation. ```{r boxplot_before_normalization, fig.width=10, fig.height=10,echo=TRUE} boxplot(log2(counts+1), las=2) ``` Normalize with voom and create model matrix ```{r} 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 ```{r boxplot_after_normalization, fig.width=10, fig.height=10,echo=TRUE} boxplot(eset_voom$E, las=2, names= colnames(eset_voom$E)) ``` Build Linear model. ```{r, eval=TRUE} 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) 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. ```{r MDS, fig.width=10, fig.height=10,eval=TRUE} 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) ``` Run Hierarchical Clustering on the voom normalized logCPMs. ```{r} dat.dist<-dist(t(eset_voom$E)) plot(hclust(dat.dist), cex=.75) ``` ```{r heatmap1, fig.width=10, fig.height=8,echo=TRUE} dat.dist<-dist(t(eset_voom$E)) plot(hclust(dat.dist), labels = targets$Treatment, cex=.75) ``` Print Session Info ```{r echo=TRUE} sessionInfo() ```