--- title: "RNAseq DE Analysis" output: word_document: fig_height: 8 fig_width: 8 pdf_document: default html_document: default --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) knitr::opts_chunk$set(cache=TRUE) setwd("USER_PATH") ``` ### Treatment name legend R1 = T1.30°(±0) R2 = T1.30°(±5) R3 = T1.30°(±10) R1S1 = T2.30°(±0) R1S2 = T2.35°(±0) R1S3 = T2.40°(±0) R2S2 = T2.35°(±5) R3S3 = T2.40°(±10) ### Load counts data ```{r, message = FALSE} library(DESeq2) MS_table <- read.table("MS_sample_table.txt", sep="\t", header=TRUE, colClasses=c(rep("character",2),rep("factor",8),rep("character",2),"factor")) directory = "RNAseq_counts" dds_full <- DESeqDataSetFromHTSeqCount(sampleTable=MS_table, directory=directory, design=~group) raw_counts = as.data.frame(counts(dds_full)) barplot(colSums(raw_counts), las=2, col=c(rep("tomato",3),rep("tomato4",3),rep("steelblue2",3),rep("steelblue4",3))) title("raw read counts per sample") ``` ### Preliminary analysis ```{r, message = FALSE} library(vsn) library(ggplot2) dds <- collapseReplicates(dds_full, dds_full$sample_rep) cds <- estimateSizeFactors(dds) cds <- estimateDispersions(cds) vsd = varianceStabilizingTransformation(cds, blind=TRUE) theme_set(theme_bw()) meanSdPlot(assay(vsd)) plotDispEsts(cds) library(ggthemes) nudge <- position_nudge(y = 10) PCA_data <- plotPCA(vsd, intgroup = c("treatment", "pop"), returnData = TRUE) percentVar <- round(100 * attr(PCA_data, "percentVar")) ggplot(PCA_data, aes(x = PC1, y = PC2, color = treatment, shape = pop, position_nudge(y=10))) + geom_point(size =3, position = position_jitter(w=0.05, h=0.05)) + xlab(paste0("PC1: ", percentVar[1], "% variance")) + ylab(paste0("PC2: ", percentVar[2], "% variance")) + coord_fixed()+ theme(panel.background = element_rect(fill = "white", colour = "black"))+ theme(plot.title = element_text(hjust = 0.5)) ``` ### Normalize read counts and run analysis ```{r, message=FALSE} norm_counts = as.data.frame(counts(cds, normalized=TRUE)) barplot(colSums(norm_counts), las=2, col=c(rep("tomato",3),rep("tomato4",3),rep("steelblue2",3),rep("steelblue4",3))) title("normalized read counts per sample, collapsed TR") data <- nbinomWaldTest(cds, maxit=100000) ``` ### Comparison 1: Effect of Rearing Temperature ```{r, message=FALSE} #lab RTa.L = results(data, contrast=c("group", "R2L", "R1L"), alpha=0.01) RTa.L = lfcShrink(data, contrast=c("group", "R2L", "R1L"), res=RTa.L) RTa.L.filter = RTa.L[which(abs(RTa.L$log2FoldChange) > 1 & RTa.L$padj < 0.01),] summary(RTa.L.filter) RTb.L = results(data, contrast=c("group", "R3L", "R1L"), alpha=0.01) RTb.L = lfcShrink(data, contrast=c("group", "R3L", "R1L"), res=RTb.L) RTb.L.filter = RTb.L[which(abs(RTb.L$log2FoldChange) > 1 & RTb.L$padj < 0.01),] summary(RTb.L.filter) #field RTa.F = results(data, contrast=c("group", "R2F", "R1F"), alpha=0.01) RTa.F = lfcShrink(data, contrast=c("group", "R2F", "R1F"), res=RTa.F) RTa.F.filter = RTa.F[which(abs(RTa.F$log2FoldChange) > 1 & RTa.F$padj < 0.01),] summary(RTa.F.filter) RTb.F = results(data, contrast=c("group", "R3F", "R1F"), alpha=0.01) RTb.F = lfcShrink(data, contrast=c("group", "R3F", "R1F"), res=RTb.F) RTb.F.filter = RTb.F[which(abs(RTb.F$log2FoldChange) > 1 & RTb.F$padj < 0.01),] summary(RTb.F.filter) par(mfrow=c(2,2)) plot(RTa.L$log2FoldChange ~ RTa.L$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray70", main="Lab: R2 vs R1", ylim=c(-5,8)) abline(h=0, col='red') points(RTa.L.filter$log2FoldChange ~ RTa.L.filter$baseMean, col = "red", pch = 16, cex=0.7) plot(RTb.L$log2FoldChange ~ RTb.L$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray70", main="Lab: R3 vs R1", ylim=c(-5,8)) abline(h=0, col='red') points(RTb.L.filter$log2FoldChange ~ RTb.L.filter$baseMean, col = "red", pch = 16, cex=0.7) plot(RTa.F$log2FoldChange ~ RTa.F$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray70", main="Field: R2 vs R1", ylim=c(-5,8)) abline(h=0, col='red') points(RTa.F.filter$log2FoldChange ~ RTa.F.filter$baseMean, col = "red", pch = 16, cex=0.7) plot(RTb.F$log2FoldChange ~ RTb.F$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray70", main="Field: R3 vs R1", ylim=c(-5,8)) abline(h=0, col='red') points(RTb.F.filter$log2FoldChange ~ RTb.F.filter$baseMean, col = "red", pch = 16, cex=0.7) ``` ### Comparison 2: Effect of Sample Time ```{r, message=FALSE} STa.L = results(data, contrast=c("group", "R1S1L", "R1L"), alpha=0.01) STa.L = lfcShrink(data, contrast=c("group", "R1S1L", "R1L"), res=STa.L) STa.L.filter = STa.L[which(abs(STa.L$log2FoldChange) > 1 & STa.L$padj < 0.01),] summary(STa.L.filter) STa.F = results(data, contrast=c("group", "R1S1F", "R1F"), alpha=0.01) STa.F = lfcShrink(data, contrast=c("group", "R1S1F", "R1F"), res=STa.F) STa.F.filter = STa.F[which(abs(STa.F$log2FoldChange) > 1 & STa.F$padj < 0.01),] summary(STa.F.filter) par(mfrow=c(2,2)) plot(STa.L$log2FoldChange ~ STa.L$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray32", main="Lab: R1S1 vs R1", ylim=c(-5,5)) abline(h=0, col='red') points(STa.L.filter$log2FoldChange ~ STa.L.filter$baseMean, col = "red", pch = 16, cex=0.7) plot(STa.F$log2FoldChange ~ STa.F$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray32", main="Field: R1S1 vs R1", ylim=c(-5,5)) abline(h=0, col='red') points(STa.F.filter$log2FoldChange ~ STa.F.filter$baseMean, col = "red", pch = 16, cex=0.7) ``` ### Comparison 3: Novel heat shock vs control ```{r, message=FALSE} Sa.L = results(data, contrast=c("group", "R1S2L", "R1S1L"), alpha=0.01) Sa.L = lfcShrink(data, contrast=c("group", "R1S2L", "R1S1L"), res=Sa.L) Sa.L.filter = Sa.L[which(abs(Sa.L$log2FoldChange) > 1 & Sa.L$padj < 0.01),] summary(Sa.L.filter) Sb.L = results(data, contrast=c("group", "R1S3L", "R1S1L"), alpha=0.01) Sb.L = lfcShrink(data, contrast=c("group", "R1S3L", "R1S1L"), res=Sb.L) Sb.L.filter = Sb.L[which(abs(Sb.L$log2FoldChange) > 1 & Sb.L$padj < 0.01),] summary(Sb.L.filter) Sa.F = results(data, contrast=c("group", "R1S2F", "R1S1F"), alpha=0.01) Sa.F = lfcShrink(data, contrast=c("group", "R1S2F", "R1S1F"), res=Sa.F) Sa.F.filter = Sa.F[which(abs(Sa.F$log2FoldChange) > 1 & Sa.F$padj < 0.01),] summary(Sa.F.filter) Sb.F = results(data, contrast=c("group", "R1S3F", "R1S1F"), alpha=0.01) Sb.F = lfcShrink(data, contrast=c("group", "R1S3F", "R1S1F"), res=Sb.F) Sb.F.filter = Sb.F[which(abs(Sb.F$log2FoldChange) > 1 & Sb.F$padj < 0.01),] summary(Sb.F.filter) par(mfrow=c(2,2)) plot(Sa.L$log2FoldChange ~ Sa.L$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray32", main="Lab: R1S2 vs R1S1", ylim=c(-5,8)) abline(h=0, col='red') points(Sa.L.filter$log2FoldChange ~ Sa.L.filter$baseMean, col = "red", pch = 16, cex=0.7) plot(Sb.L$log2FoldChange ~ Sb.L$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray32", main="Lab: R1S3 vs R1S1", ylim=c(-5,8)) abline(h=0, col='red') points(Sb.L.filter$log2FoldChange ~ Sb.L.filter$baseMean, col = "red", pch = 16, cex=0.7) plot(Sa.F$log2FoldChange ~ Sa.F$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray32", main="Field: R1S2 vs R1S1", ylim=c(-5,8)) abline(h=0, col='red') points(Sa.F.filter$log2FoldChange ~ Sa.F.filter$baseMean, col = "red", pch = 16, cex=0.7) plot(Sb.F$log2FoldChange ~ Sb.F$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray32", main="Field: R1S3 vs R1S1", ylim=c(-5,8)) abline(h=0, col='red') points(Sb.F.filter$log2FoldChange ~ Sb.F.filter$baseMean, col = "red", pch = 16, cex=0.7) ``` ### Comparison 4: Novel heat shock vs thermal history ```{r, message=FALSE} SvAa.L = results(data, contrast=c("group", "R1S2L", "R2S2L"), alpha=0.01) SvAa.L = lfcShrink(data, contrast=c("group", "R1S2L", "R2S2L"), res=SvAa.L) SvAa.L.filter = SvAa.L[which(abs(SvAa.L$log2FoldChange) > 1 & SvAa.L$padj < 0.01),] summary(SvAa.L.filter) SvAb.L = results(data, contrast=c("group", "R1S3L", "R3S3L"), alpha=0.01) SvAb.L = lfcShrink(data, contrast=c("group", "R1S3L", "R3S3L"), res=SvAb.L) SvAb.L.filter = SvAb.L[which(abs(SvAb.L$log2FoldChange) > 1 & SvAb.L$padj < 0.01),] summary(SvAb.L.filter) SvAa.F = results(data, contrast=c("group", "R1S2F", "R2S2F"), alpha=0.01) SvAa.F = lfcShrink(data, contrast=c("group", "R1S2F", "R2S2F"), res=SvAa.F) SvAa.F.filter = SvAa.F[which(abs(SvAa.F$log2FoldChange) > 1 & SvAa.F$padj < 0.01),] summary(SvAa.F.filter) SvAb.F = results(data, contrast=c("group", "R1S3F", "R3S3F"), alpha=0.01) SvAb.F = lfcShrink(data, contrast=c("group", "R1S3F", "R3S3F"), res=SvAb.F) SvAb.F.filter = SvAb.F[which(abs(SvAb.F$log2FoldChange) > 1 & SvAb.F$padj < 0.01),] summary(SvAb.F.filter) par(mfrow=c(2,2)) plot(SvAa.L$log2FoldChange ~ SvAa.L$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray32", main="Lab: R1S2 vs R2S2", ylim=c(-5,8)) abline(h=0, col='red') points(SvAa.L.filter$log2FoldChange ~ SvAa.L.filter$baseMean, col = "red", pch = 16, cex=0.7) plot(SvAb.L$log2FoldChange ~ SvAb.L$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray32", main="Lab: R1S3 vs R3S3", ylim=c(-5,8)) abline(h=0, col='red') points(SvAb.L.filter$log2FoldChange ~ SvAb.L.filter$baseMean, col = "red", pch = 16, cex=0.7) plot(SvAa.F$log2FoldChange ~ SvAa.F$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray32", main="Field: R1S2 vs R2S2", ylim=c(-5,8)) abline(h=0, col='red') points(SvAa.F.filter$log2FoldChange ~ SvAa.F.filter$baseMean, col = "red", pch = 16, cex=0.7) plot(SvAb.F$log2FoldChange ~ SvAb.F$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray32", main="Field: R1S3 vs R3S3", ylim=c(-5,8)) abline(h=0, col='red') points(SvAb.F.filter$log2FoldChange ~ SvAb.F.filter$baseMean, col = "red", pch = 16, cex=0.7) ``` ### Comparison 5: Thermal history vs control ```{r, message=FALSE} HSa.L = results(data, contrast=c("group", "R2S2L", "R1S1L"), alpha=0.01) HSa.L = lfcShrink(data, contrast=c("group", "R2S2L", "R1S1L"), res=HSa.L) HSa.L.filter = HSa.L[which(abs(HSa.L$log2FoldChange) > 1 & HSa.L$padj < 0.01),] summary(HSa.L.filter) HSb.L = results(data, contrast=c("group", "R3S3L", "R1S1L"), alpha=0.01) HSb.L = lfcShrink(data, contrast=c("group", "R3S3L", "R1S1L"), res=HSb.L) HSb.L.filter = HSb.L[which(abs(HSb.L$log2FoldChange) > 1 & HSb.L$padj < 0.01),] summary(HSb.L.filter) HSa.F = results(data, contrast=c("group", "R2S2F", "R1S1F"), alpha=0.01) HSa.F = lfcShrink(data, contrast=c("group", "R2S2F", "R1S1F"), res=HSa.F) HSa.F.filter = HSa.F[which(abs(HSa.F$log2FoldChange) > 1 & HSa.F$padj < 0.01),] summary(HSa.F.filter) HSb.F = results(data, contrast=c("group", "R3S3F", "R1S1F"), alpha=0.01) HSb.F = lfcShrink(data, contrast=c("group", "R3S3F", "R1S1F"), res=HSb.F) HSb.F.filter = HSb.F[which(abs(HSb.F$log2FoldChange) > 1 & HSb.F$padj < 0.01),] summary(HSb.F.filter) par(mfrow=c(2,2)) plot(HSa.L$log2FoldChange ~ HSa.L$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray32", main="Lab: R2S2 vs R1S1", ylim=c(-5,8)) abline(h=0, col='red') points(HSa.L.filter$log2FoldChange ~ HSa.L.filter$baseMean, col = "red", pch = 16, cex=0.7) plot(HSb.L$log2FoldChange ~ HSb.L$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray32", main="Lab: R3S3 vs R1S1", ylim=c(-5,8)) abline(h=0, col='red') points(HSb.L.filter$log2FoldChange ~ HSb.L.filter$baseMean, col = "red", pch = 16, cex=0.7) plot(HSa.F$log2FoldChange ~ HSa.F$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray32", main="Field: R2S2 vs R1S1", ylim=c(-5,8)) abline(h=0, col='red') points(HSa.F.filter$log2FoldChange ~ HSa.F.filter$baseMean, col = "red", pch = 16, cex=0.7) plot(HSb.F$log2FoldChange ~ HSb.F$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray32", main="Field: R3S3 vs R1S1", ylim=c(-5,8)) abline(h=0, col='red') points(HSb.F.filter$log2FoldChange ~ HSb.F.filter$baseMean, col = "red", pch = 16, cex=0.7) ``` ### Field vs lab population comparisons by treatment ```{r, message=FALSE} R1.P = results(data, contrast=c("group", "R1F", "R1L"), alpha=0.01) R1.P = lfcShrink(data, contrast=c("group", "R1F", "R1L"), res=R1.P) R1.P.filter = R1.P[which(abs(R1.P$log2FoldChange) > 1 & R1.P$padj < 0.01),] summary(R1.P.filter) R2.P = results(data, contrast=c("group", "R2F", "R2L"), alpha=0.01) R2.P = lfcShrink(data, contrast=c("group", "R2F", "R2L"), res=R2.P) R2.P.filter = R2.P[which(abs(R2.P$log2FoldChange) > 1 & R2.P$padj < 0.01),] summary(R2.P.filter) R3.P = results(data, contrast=c("group", "R3F", "R3L"), alpha=0.01) R3.P = lfcShrink(data, contrast=c("group", "R3F", "R3L"), res=R3.P) R3.P.filter = R3.P[which(abs(R3.P$log2FoldChange) > 2 & R3.P$padj < 0.01),] summary(R3.P.filter) R1S1.P = results(data, contrast=c("group", "R1S1F", "R1S1L"), alpha=0.01) R1S1.P = lfcShrink(data, contrast=c("group", "R1S1F", "R1S1L"), res=R1S1.P) R1S1.P.filter = R1S1.P[which(abs(R1S1.P$log2FoldChange) > 1 & R1S1.P$padj < 0.01),] summary(R1S1.P.filter) R2S2.P = results(data, contrast=c("group", "R2S2F", "R2S2L"), alpha=0.01) R2S2.P = lfcShrink(data, contrast=c("group", "R2S2F", "R2S2L"), res=R2S2.P) R2S2.P.filter = R2S2.P[which(abs(R2S2.P$log2FoldChange) > 1 & R2S2.P$padj < 0.01),] summary(R2S2.P.filter) R3S3.P = results(data, contrast=c("group", "R3S3F", "R3S3L"), alpha=0.01) R3S3.P = lfcShrink(data, contrast=c("group", "R3S3F", "R3S3L"), res=R3S3.P) R3S3.P.filter = R3S3.P[which(abs(R3S3.P$log2FoldChange) > 1 & R3S3.P$padj < 0.01),] summary(R3S3.P.filter) R1S2.P= results(data, contrast=c("group", "R1S2F", "R1S2L"), alpha=0.01) R1S2.P = lfcShrink(data, contrast=c("group", "R1S2F", "R1S2L"), res=R1S2.P) R1S2.P.filter = R1S2.P[which(abs(R1S2.P$log2FoldChange) > 1 & R1S2.P$padj < 0.01),] summary(R1S2.P.filter) R1S3.P = results(data, contrast=c("group", "R1S3F", "R1S3L"), alpha=0.01) R1S3.P = lfcShrink(data, contrast=c("group", "R1S3F", "R1S3L"), res=R1S3.P) R1S3.P.filter = R1S3.P[which(abs(R1S3.P$log2FoldChange) > 1 & R1S3.P$padj < 0.01),] summary(R1S3.P.filter) par(mfrow=c(3,3)) plot(R1.P$log2FoldChange ~ R1.P$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray32", main="R1: F vs L", ylim=c(-10,10)) abline(h=0, col='red') points(R1.P.filter$log2FoldChange ~ R1.P.filter$baseMean, col = "red", pch = 16, cex=0.7) plot(R2.P$log2FoldChange ~ R2.P$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray32", main="R2: F vs L", ylim=c(-10,10)) abline(h=0, col='red') points(R2.P.filter$log2FoldChange ~ R2.P.filter$baseMean, col = "red", pch = 16, cex=0.7) plot(R3.P$log2FoldChange ~ R3.P$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray32", main="R3: F vs L", ylim=c(-10,10)) abline(h=0, col='red') points(R3.P.filter$log2FoldChange ~ R3.P.filter$baseMean, col = "red", pch = 16, cex=0.7) plot(R1S1.P$log2FoldChange ~ R1S1.P$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray32", main="R1S1: F vs L", ylim=c(-10,10)) abline(h=0, col='red') points(R1S1.P.filter$log2FoldChange ~ R1S1.P.filter$baseMean, col = "red", pch = 16, cex=0.7) plot(R2S2.P$log2FoldChange ~ R2S2.P$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray32", main="R2S2: F vs L", ylim=c(-10,10)) abline(h=0, col='red') points(R2S2.P.filter$log2FoldChange ~ R2S2.P.filter$baseMean, col = "red", pch = 16, cex=0.7) plot(R3S3.P$log2FoldChange ~ R3S3.P$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray32", main="R3S3: F vs L", ylim=c(-10,10)) abline(h=0, col='red') points(R3S3.P.filter$log2FoldChange ~ R3S3.P.filter$baseMean, col = "red", pch = 16, cex=0.7) plot(R1S2.P$log2FoldChange ~ R1S2.P$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray32", main="R1S2: F vs L", ylim=c(-10,10)) abline(h=0, col='red') points(R1S2.P.filter$log2FoldChange ~ R1S2.P.filter$baseMean, col = "red", pch = 16, cex=0.7) plot(R1S3.P$log2FoldChange ~ R1S3.P$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray32", main="R1S3: F vs L", ylim=c(-10,10)) abline(h=0, col='red') points(R1S3.P.filter$log2FoldChange ~ R1S3.P.filter$baseMean, col = "red", pch = 16, cex=0.7) ``` ### Field vs lab population comparison all ```{r, message=FALSE} MS_table$pop <- relevel(MS_table$pop, "lab") dds2 <- DESeqDataSetFromHTSeqCount(sampleTable=MS_table, directory=directory, design=~pop) dds2 <- collapseReplicates(dds2, dds2$sample_rep) data2 <- estimateSizeFactors(dds2) data2 <- estimateDispersions(data2) data2 <- nbinomWaldTest(data2, maxit=100000) Pop = results(data2, contrast=c("pop", "field", "lab"), alpha=0.01, cooksCutoff=50) Pop = lfcShrink(data2, contrast=c("pop", "field", "lab"), res=Pop) Pop.filter = Pop[which(abs(Pop$log2FoldChange) > 1 & Pop$padj < 0.01),] summary(Pop.filter) plot(Pop$log2FoldChange ~ Pop$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray32", main="Field vs Lab", ylim=c(-10,10)) abline(h=0, col='red') points(Pop.filter$log2FoldChange ~ Pop.filter$baseMean, col = "red", pch = 16, cex=0.7) ``` ### Heat shock proteins ```{r, message=FALSE} HSP <- scan("HSP_list.txt", what="", sep="\n") RTa.L.HSP <- RTa.L[HSP,] RTa.L.HSP_filter = RTa.L.HSP[which(abs(RTa.L.HSP$log2FoldChange) > 1 & RTa.L.HSP$padj < 0.01),] summary(RTa.L.HSP_filter) RTb.L.HSP <- RTb.L[HSP,] RTb.L.HSP_filter = RTb.L.HSP[which(abs(RTb.L.HSP$log2FoldChange) > 1 & RTb.L.HSP$padj < 0.01),] summary(RTb.L.HSP_filter) RTa.F.HSP <- RTa.F[HSP,] RTa.F.HSP_filter = RTa.F.HSP[which(abs(RTa.F.HSP$log2FoldChange) > 1 & RTa.F.HSP$padj < 0.01),] summary(RTa.F.HSP_filter) RTb.F.HSP <- RTb.F[HSP,] RTb.F.HSP_filter = RTb.F.HSP[which(abs(RTb.F.HSP$log2FoldChange) > 1 & RTb.F.HSP$padj < 0.01),] summary(RTb.F.HSP_filter) Sa.L.HSP <- Sa.L[HSP,] Sa.L.HSP_filter = Sa.L.HSP[which(abs(Sa.L.HSP$log2FoldChange) > 1 & Sa.L.HSP$padj < 0.01),] summary(Sa.L.HSP_filter) Sb.L.HSP <- Sb.L[HSP,] Sb.L.HSP_filter = Sb.L.HSP[which(abs(Sb.L.HSP$log2FoldChange) > 1 & Sb.L.HSP$padj < 0.01),] summary(Sb.L.HSP_filter) Sa.F.HSP <- Sa.F[HSP,] Sa.F.HSP_filter = Sa.F.HSP[which(abs(Sa.F.HSP$log2FoldChange) > 1 & Sa.F.HSP$padj < 0.01),] summary(Sa.F.HSP_filter) Sb.F.HSP <- Sb.F[HSP,] Sb.F.HSP_filter = Sb.F.HSP[which(abs(Sb.F.HSP$log2FoldChange) > 1 & Sb.F.HSP$padj < 0.01),] summary(Sb.F.HSP_filter) SvAa.L.HSP <- SvAa.L[HSP,] SvAa.L.HSP_filter = SvAa.L.HSP[which(abs(SvAa.L.HSP$log2FoldChange) > 1 & SvAa.L.HSP$padj < 0.01),] summary(SvAa.L.HSP_filter) SvAb.L.HSP <- SvAb.L[HSP,] SvAb.L.HSP_filter = SvAb.L.HSP[which(abs(SvAb.L.HSP$log2FoldChange) > 1 & SvAb.L.HSP$padj < 0.01),] summary(SvAb.L.HSP_filter) SvAa.F.HSP <- SvAa.F[HSP,] SvAa.F.HSP_filter = SvAa.F.HSP[which(abs(SvAa.F.HSP$log2FoldChange) > 1 & SvAa.F.HSP$padj < 0.01),] summary(SvAa.F.HSP_filter) SvAb.F.HSP <- SvAb.F[HSP,] SvAb.F.HSP_filter = SvAb.F.HSP[which(abs(SvAb.F.HSP$log2FoldChange) > 1 & SvAb.F.HSP$padj < 0.01),] summary(SvAb.F.HSP_filter) par(mfrow=c(2,2)) plot(RTa.L$log2FoldChange ~ RTa.L$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray60", main="Lab: R2 vs R1", ylim=c(-5,8)) abline(h=0, col='red') points(RTa.L.filter$log2FoldChange ~ RTa.L.filter$baseMean, col = "red", pch = 16, cex=0.7) points(RTa.L.HSP_filter$log2FoldChange ~ RTa.L.HSP_filter$baseMean, col="royalblue1", pch=16, cex=0.7) plot(RTb.L$log2FoldChange ~ RTb.L$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray60", main="Lab: R3 vs R1", ylim=c(-5,8)) abline(h=0, col='red') points(RTb.L.filter$log2FoldChange ~ RTb.L.filter$baseMean, col = "red", pch = 16, cex=0.7) points(RTb.L.HSP_filter$log2FoldChange ~ RTb.L.HSP_filter$baseMean, col="royalblue1", pch=16, cex=0.7) plot(RTa.F$log2FoldChange ~ RTa.F$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray60", main="Field: R2 vs R1", ylim=c(-5,8)) abline(h=0, col='red') points(RTa.F.filter$log2FoldChange ~ RTa.F.filter$baseMean, col = "red", pch = 16, cex=0.7) points(RTa.F.HSP_filter$log2FoldChange ~ RTa.F.HSP_filter$baseMean, col="royalblue1", pch=16, cex=0.7) plot(RTb.F$log2FoldChange ~ RTb.F$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray60", main="Field: R3 vs R1", ylim=c(-5,8)) abline(h=0, col='red') points(RTb.F.filter$log2FoldChange ~ RTb.F.filter$baseMean, col = "red", pch = 16, cex=0.7) points(RTb.F.HSP_filter$log2FoldChange ~ RTb.F.HSP_filter$baseMean, col="royalblue1", pch=16, cex=0.7) par(mfrow=c(2,2)) plot(Sa.L$log2FoldChange ~ Sa.L$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray60", main="Lab: R1S2 vs R1S1", ylim=c(-5,8)) abline(h=0, col='red') points(Sa.L.filter$log2FoldChange ~ Sa.L.filter$baseMean, col = "red", pch = 16, cex=0.7) points(Sa.L.HSP_filter$log2FoldChange ~ Sa.L.HSP_filter$baseMean, col="royalblue1", pch=16, cex=0.7) plot(Sb.L$log2FoldChange ~ Sb.L$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray60", main="Lab: R1S3 vs R1S1", ylim=c(-5,8)) abline(h=0, col='red') points(Sb.L.filter$log2FoldChange ~ Sb.L.filter$baseMean, col = "red", pch = 16, cex=0.7) points(Sb.L.HSP_filter$log2FoldChange ~ Sb.L.HSP_filter$baseMean, col="royalblue1", pch=16, cex=0.7) plot(Sa.F$log2FoldChange ~ Sa.F$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray60", main="Field: R1S2 vs R1S1", ylim=c(-5,8)) abline(h=0, col='red') points(Sa.F.filter$log2FoldChange ~ Sa.F.filter$baseMean, col = "red", pch = 16, cex=0.7) points(Sa.F.HSP_filter$log2FoldChange ~ Sa.F.HSP_filter$baseMean, col="royalblue1", pch=16, cex=0.7) plot(Sb.F$log2FoldChange ~ Sb.F$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray60", main="Field: R1S3 vs R1S1", ylim=c(-5,8)) abline(h=0, col='red') points(Sb.F.filter$log2FoldChange ~ Sb.F.filter$baseMean, col = "red", pch = 16, cex=0.7) points(Sb.F.HSP_filter$log2FoldChange ~ Sb.F.HSP_filter$baseMean, col="royalblue1", pch=16, cex=0.7) par(mfrow=c(2,2)) plot(SvAa.L$log2FoldChange ~ SvAa.L$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray60", main="Lab: R1S2 vs R2S2", ylim=c(-5,8)) abline(h=0, col='red') points(SvAa.L.filter$log2FoldChange ~ SvAa.L.filter$baseMean, col = "red", pch = 16, cex=0.7) points(SvAa.L.HSP_filter$log2FoldChange ~ SvAa.L.HSP_filter$baseMean, col="royalblue1", pch=16, cex=0.7) plot(SvAb.L$log2FoldChange ~ SvAb.L$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray60", main="Lab: R1S3 vs R3S3", ylim=c(-5,8)) abline(h=0, col='red') points(SvAb.L.filter$log2FoldChange ~ SvAb.L.filter$baseMean, col = "red", pch = 16, cex=0.7) points(SvAb.L.HSP_filter$log2FoldChange ~ SvAb.L.HSP_filter$baseMean, col="royalblue1", pch=16, cex=0.7) plot(SvAa.F$log2FoldChange ~ SvAa.F$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray60", main="Field: R1S2 vs R2S2", ylim=c(-5,8)) abline(h=0, col='red') points(SvAa.F.filter$log2FoldChange ~ SvAa.F.filter$baseMean, col = "red", pch = 16, cex=0.7) points(SvAa.F.HSP_filter$log2FoldChange ~ SvAa.F.HSP_filter$baseMean, col="royalblue1", pch=16, cex=0.7) plot(SvAb.F$log2FoldChange ~ SvAb.F$baseMean,log='x', ylab="log2 Fold Change", xlab="Mean Expression", pch=16, cex=0.6, col="gray60", main="Field: R1S3 vs R3S3", ylim=c(-5,8)) abline(h=0, col='red') points(SvAb.F.filter$log2FoldChange ~ SvAb.F.filter$baseMean, col = "red", pch = 16, cex=0.7) points(SvAb.F.HSP_filter$log2FoldChange ~ SvAb.F.HSP_filter$baseMean, col="royalblue1", pch=16, cex=0.7) R1.P.HSP <- R1.P[HSP,] R1.P.HSP_filter = R1.P.HSP[which(abs(R1.P.HSP$log2FoldChange) > 1 & R1.P.HSP$padj < 0.01),] summary(R1.P.HSP_filter) R2.P.HSP <- R2.P[HSP,] R2.P.HSP_filter = R2.P.HSP[which(abs(R2.P.HSP$log2FoldChange) > 1 & R2.P.HSP$padj < 0.01),] R3.P.HSP <- R3.P[HSP,] R3.P.HSP_filter = R3.P.HSP[which(abs(R3.P.HSP$log2FoldChange) > 1 & R3.P.HSP$padj < 0.01),] R1S1.P.HSP <- R1S1.P[HSP,] R1S1.P.HSP_filter = R1S1.P.HSP[which(abs(R1S1.P.HSP$log2FoldChange) > 1 & R1S1.P.HSP$padj < 0.01),] R2S2.P.HSP <- R2S2.P[HSP,] R2S2.P.HSP_filter = R2S2.P.HSP[which(abs(R2S2.P.HSP$log2FoldChange) > 1 & R2S2.P.HSP$padj < 0.01),] R3S3.P.HSP <- R3S3.P[HSP,] R3S3.P.HSP_filter = R3S3.P.HSP[which(abs(R3S3.P.HSP$log2FoldChange) > 1 & R3S3.P.HSP$padj < 0.01),] R1S2.P.HSP <- R1S2.P[HSP,] R1S2.P.HSP_filter = R1S2.P.HSP[which(abs(R1S2.P.HSP$log2FoldChange) > 1 & R1S2.P.HSP$padj < 0.01),] R1S3.P.HSP <- R1S3.P[HSP,] R1S3.P.HSP_filter = R1S3.P.HSP[which(abs(R1S3.P.HSP$log2FoldChange) > 1 & R1S3.P.HSP$padj < 0.01),] ```