RNAseq in Bioconductor exercises (part 2)

In todays session we will work with some of the RNAseq data of adult mouse tissues from Bing Ren’s lab, Liver and Heart.

• More information on liver data can be found here

• More information on heart data can be found here

• Precounted RNAseq reads in genes for these tissues can be found as an R data object in data/geneCounts_Tissue.RData.

• Salmon quantification for these tissues can be found in directory data/Salmon_Tissue/.

Exercises

1. Get counts into DESeq2

Load the data/geneCounts_Tissue.RData and create a DESeq2 object. Add 0.25 to every counts and plot a boxplot of the log2 of these updated counts.

library(DESeq2)
library(tximport)
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
library(RColorBrewer)

load("data/geneCounts_Tissue.RData")

colData(geneCounts_Tissue)$Tissue <- factor(c("Heart","Heart","Liver","Liver"))

ddsT <- DESeqDataSet(geneCounts_Tissue,design=~Tissue)

## renaming the first element in assays to 'counts'

boxplot(log2(counts(ddsT,normalized=FALSE)+0.25), names=c("Heart_1","Heart_2","Liver_1","Liver_2"), col=brewer.pal(4, "Paired"))

2. Normalize counts

Run the DEseq workflow function and retrieve normalized counts. Add 0.25 to normalized counts and plot a boxplot of the log2 of these updated counts.

ddsT <- DESeq(ddsT)
boxplot(log2(counts(ddsT,normalized=TRUE)+0.25), names=c("Heart_1","Heart_2","Liver_1","Liver_2"), col=brewer.pal(4, "Paired"))

3. Dispersion

Plot the dispersion fit and estimatesn for our DEseq2 object.

plotDispEsts(ddsT)

4. Multiple testing

Add a new padj value for all genes. Make a barplot of the number of significantly (padj < 0.05) up and down regulated genes by the original padj values and our new padj values.

library(ggplot2)
newRes <- results(ddsT,contrast = c("Tissue","Heart","Liver"))
newResDF <- as.data.frame(newRes)
newResDF$newPadj <- p.adjust(newResDF$pvalue,method = "BH")
newPadj <- ifelse(newResDF$newPadj < 0.05 & newResDF$log2FoldChange > 0,"HeartUp",
              ifelse(newResDF$newPadj < 0.05 & newResDF$log2FoldChange < 0,"HeartDown","NoChange"))
originalPadj <- ifelse(newResDF$padj < 0.05 & newResDF$log2FoldChange > 0,"HeartUp",
              ifelse(newResDF$padj < 0.05 & newResDF$log2FoldChange < 0,"HeartDown","NoChange"))
newPadjFrame <- data.frame(Method="All",Padj=newPadj)
originalPadjFrame <- data.frame(Method="Filtered",Padj=originalPadj)
toPlot <- rbind(newPadjFrame,originalPadjFrame)
ggplot(toPlot,aes(x=Padj,fill=Method))+geom_bar(position="dodge")+theme_bw()

5. MA plots

Produce an MA plot before and after logFC shrinkage.

newRes <- results(ddsT,contrast = c("Tissue","Liver","Heart"))
DESeq2::plotMA(newRes)

newResLfc <- lfcShrink(dds=ddsT,coef = 2)
DESeq2::plotMA(newResLfc)