RNA-seq analysis in R

Pre-processsing RNA-seq data - Challenge Solutions

Data

library(DESeq2)
library(gplots)
library(RColorBrewer)
library(limma)
library(tidyverse)

# Read the sample information into R
sampleinfo <- read.delim("../data/SampleInfo.Corrected.txt", 
                         stringsAsFactors = F)
# Read the data into R
seqdata <- read.delim("../data/GSE60450_Lactation.featureCounts", 
                      stringsAsFactors = F,
                      comment = "#")
# Remove first two columns from seqdata
countdata <- as.data.frame(seqdata) %>% 
    column_to_rownames("Geneid") %>% # turn the geneid column into rownames
    rename_all(str_remove, ".bam") %>% # remove the ".bam" from the column names
    select(sampleinfo$Sample) %>% # keep sample columns using sampleinfo
    as.matrix()
# filter data
keep <- rowSums(countdata) > 5
countdata <- countdata[keep,]
# rlogcounts
rlogcounts <- rlog(countdata)
# We estimate the variance for each row in the logcounts matrix
countVar <- apply(rlogcounts, 1, var)
# DGE list
design <- as.formula(~ CellType)
# create the DESeqDataSet object
ddsObj <- DESeqDataSetFromMatrix(countData = countdata,
                              colData = sampleinfo,
                              design = design)

some variables in design formula are characters, converting to factors

ddsObj <- estimateSizeFactors(ddsObj)
# normalise dcounts
logcounts <- log2(countdata + 1)
normalizedCounts <- counts(ddsObj, normalized=TRUE) 
logNormalizedCounts <- log2(normalizedCounts + 1)

Challenge 1

1. Use the DESeq2 function rlog to transform the count data. This function also normalises for library size.
2. Plot the count distribution boxplots with this data How has this effected the count distributions?

rlogcounts <- rlog(countdata)
statusCol <- as.numeric(factor(sampleinfo$Status)) + 1
# Check distributions of samples using boxplots
boxplot(rlogcounts, 
        xlab="", 
        ylab="Log2(Counts)",
        las=2,
        col=statusCol)
# Let's add a blue horizontal line that corresponds to the median logCPM
abline(h=median(as.matrix(logcounts)), col="blue")

Challenge 2

Redo the heatmap using the top 500 LEAST variable genes. Change the colour scheme to “PiYG” and redo the heatmap. Try ?RColorBrewer and see what other colour schemes are available. Change the sample names to group using the labCol argument Remove the gene names from the righthand side of the plot using labRow

Solution

# Get the gene names for the top 500 least variable genes
lowVar <- order(countVar)[1:500]
# Subset logcounts matrix
hmData <- rlogcounts[lowVar,]
## Get some nicer colours
mypalette <- brewer.pal(11,"PiYG")
## http://colorbrewer2.org/#type=sequential&scheme=BuGn&n=3
morecols <- colorRampPalette(mypalette)
# Set up colour vector for celltype variable
col.cell <- c("purple","orange")[sampleinfo$CellType]
# Plot the heatmap
heatmap.2(hmData, 
          col=rev(morecols(50)),
          trace="none", 
          main="Top 500 most variable genes across samples",
          ColSideColors=col.cell,scale="row",
          labCol=sampleinfo$Group, 
          labRow = NA)

Challenge 3

Plot the biased and unbiased MA plots both samples side by side to see the before and after normalisation.

par(mfrow=c(2,2))
plotMA(logcounts, array = 7, main="MCL1.LA - Raw")
abline(h=0,col="grey")

plotMA(logNormalizedCounts, array = 7, main="MCL1.LA - Normalised")
abline(h=0,col="grey")

plotMA(logcounts, array = 11, main="MCL1.LE - Raw")
abline(h=0,col="grey")

plotMA(logNormalizedCounts, array = 11, main="MCL1.LE - Normalised")
abline(h=0,col="grey")