1Introduction
In this document, we show how to conduct Exploratory Data Analysis (EDA) and normalization for a typical RNA-Seq experiment using the package EDASeq
.
One can think of EDA for RNA-Seq as a two-step process: “read-level” EDA helps in discovering lanes with low sequencing depths, quality issues, and unusual nucleotide frequencies, while “gene-level’’ EDA can capture mislabeled lanes, issues with distributional assumptions (e.g., over-dispersion), and GC-content bias.
The package also implements both “within-lane” and “between-lane” normalization procedures, to account, respectively, for within-lane gene-specific (and possibly lane-specific) effects on read counts (e.g., related to gene length or GC-content) and for between-lane distributional differences in read counts (e.g., sequencing depths).
To illustrate the functionality of the EDASeq
package, we make use of the Saccharomyces cerevisiae RNA-Seq data from (Lee et al. 2008). Briefly, a wild-type strain and three mutant strains were sequenced using the Solexa 1G Genome Analyzer. For each strain, there are four technical replicate lanes from the same library preparation. The reads were aligned using Bowtie
(Langmead et al. 2009), with unique mapping and allowing up to two mismatches.
The leeBamViews
package provides a subset of the aligned reads in BAM format. In particular, only the reads mapped between bases 800,000 and 900,000 of chromosome XIII are considered. We use these reads to illustrate read-level EDA.
The yeastRNASeq
package contains gene-level read counts for four lanes: two replicates of the wild-type strain (“wt”) and two replicates of one of the mutant strains (“mut”). We use these data to illustrate gene-level EDA.
library(EDASeq)
library(yeastRNASeq)
library(leeBamViews)
2Reading in unaligned and aligned read data
2.1Unaligned reads
Unaligned (unmapped) reads stored in FASTQ format may be managed via the class FastqFileList
imported from ShortRead
. Information related to the libraries sequenced in each lane can be stored in the elementMetadata
slot of the FastqFileList
object.
files <- list.files(file.path(system.file(package = "yeastRNASeq"),
"reads"), pattern = "fastq", full.names = TRUE)
names(files) <- gsub("\\.fastq.*", "", basename(files))
met <- DataFrame(conditions=c(rep("mut",2), rep("wt",2)),
row.names=names(files))
fastq <- FastqFileList(files)
elementMetadata(fastq) <- met
fastq
## FastqFileList of length 4
## names(4): mut_1_f mut_2_f wt_1_f wt_2_f
2.2Aligned reads
The package can deal with aligned (mapped) reads in BAM format, using the class BamFileList
from Rsamtools
. Again, the elementMetadata
slot can be used to store lane-level sample information.
files <- list.files(file.path(system.file(package = "leeBamViews"), "bam"),
pattern = "bam$", full.names = TRUE)
names(files) <- gsub("\\.bam", "", basename(files))
gt <- gsub(".*/", "", files)
gt <- gsub("_.*", "", gt)
lane <- gsub(".*(.)$", "\\1", gt)
geno <- gsub(".$", "", gt)
pd <- DataFrame(geno=geno, lane=lane,
row.names=paste(geno,lane,sep="."))
bfs <- BamFileList(files)
elementMetadata(bfs) <- pd
bfs
## BamFileList of length 8
## names(8): isowt5_13e isowt6_13e ... xrn1_13e xrn2_13e
3Read-level EDA
3.1Numbers of unaligned and aligned reads
One important check for quality control is to look at the total number of reads produced in each lane, the number and the percentage of reads mapped to a reference genome. A low total number of reads might be a symptom of low quality of the input RNA, while a low mapping percentage might indicate poor quality of the reads (low complexity), problems with the reference genome, or mislabeled lanes.
colors <- c(rep(rgb(1,0,0,alpha=0.7),2),
rep(rgb(0,0,1,alpha=0.7),2),
rep(rgb(0,1,0,alpha=0.7),2),
rep(rgb(0,1,1,alpha=0.7),2))
barplot(bfs,las=2,col=colors)
The figure, produced using the barplot
method for the BamFileList
class, displays the number of mapped reads for the subset of the yeast dataset included in the packageleeBamViews
. Unfortunately, leeBamViews
does not provide unaligned reads, but barplots of the total number of reads can be obtained using the barplot
method for the FastqFileList
class. Analogously, one can plot the percentage of mapped reads with the plot
method with signaturec(x="BamFileList", y="FastqFileList")
. See the manual pages for details.
3.2Read quality scores
As an additional quality check, one can plot the mean per-base (i.e., per-cycle) quality of the unmapped or mapped reads in every lane.
plotQuality(bfs,col=colors,lty=1)
legend("topright",unique(elementMetadata(bfs)[,1]), fill=unique(colors))
3.3Individual lane summaries
If one is interested in looking more thoroughly at one lane, it is possible to display the per-base distribution of quality scores for each lane and the number of mapped reads stratified by chromosome or strand. As expected, all the reads are mapped to chromosome XIII.
plotQuality(bfs[[1]],cex.axis=.8)
barplot(bfs[[1]],las=2)
3.4Read nucleotide distributions
A potential source of bias is related to the sequence composition of the reads. The function plotNtFrequency
plots the per-base nucleotide frequencies for all the reads in a given lane.
plotNtFrequency(bfs[[1]])
4Gene-level EDA
Examining statistics and quality metrics at a read level can help in discovering problematic libraries or systematic biases in one or more lanes. Nevertheless, some biases can be difficult to detect at this scale and gene-level EDA is equally important.
4.1Classes and methods for gene-level counts
There are several Bioconductor packages for aggregating reads over genes (or other genomic regions, such as, transcripts and exons) given a particular genome annotation, e.g., IRanges
,ShortRead
, Genominator
, Rsubread
. See their respective vignettes for details.
Here, we consider this step done and load the object geneLevelData
from yeastRNASeq
, which provides gene-level counts for 2 wild-type and 2 mutant lanes from the yeast dataset of lee2008novel
(see the Genominator
vignette for an example on the same dataset).
data(geneLevelData)
head(geneLevelData)
## mut_1 mut_2 wt_1 wt_2
## YHR055C 0 0 0 0
## YPR161C 38 39 35 34
## YOL138C 31 33 40 26
## YDR395W 55 52 47 47
## YGR129W 29 26 5 5
## YPR165W 189 180 151 180
Since it is useful to explore biases related to length and GC-content, the EDASeq
package provides, for illustration purposes, length and GC-content for S. cerevisiae genes (based on SGD annotation, version r64 (“Saccharomyces Genome Database,” n.d.)).
Functionality for automated retrieval of gene length and GC-content is introduced in the last section of the vignette.
data(yeastGC)
head(yeastGC)
## YAL001C YAL002W YAL003W YAL004W YAL005C YAL007C
## 0.3712317 0.3717647 0.4460548 0.4490741 0.4406428 0.3703704
data(yeastLength)
head(yeastLength)
## YAL001C YAL002W YAL003W YAL004W YAL005C YAL007C
## 3483 3825 621 648 1929 648
First, we filter the non-expressed genes, i.e., we consider only the genes with an average read count greater than 10 across the four lanes and for which we have length and GC-content information.
filter <- apply(geneLevelData,1,function(x) mean(x)>10)
table(filter)
## filter
## FALSE TRUE
## 1988 5077
common <- intersect(names(yeastGC),
rownames(geneLevelData[filter,]))
length(common)
## [1] 4994
This leaves us with 4994 genes.
The EDASeq
package provides the SeqExpressionSet
class to store gene counts, (lane-level) information on the sequenced libraries, and (gene-level) feature information. We use the data frame met
created in Section secRead
for the lane-level data. As for the feature data, we use gene length and GC-content.
feature <- data.frame(gc=yeastGC,length=yeastLength)
data <- newSeqExpressionSet(counts=as.matrix(geneLevelData[common,]),
featureData=feature[common,],
phenoData=data.frame(
conditions=factor(c(rep("mut",2),rep("wt",2))),
row.names=colnames(geneLevelData)))
data
## SeqExpressionSet (storageMode: lockedEnvironment)
## assayData: 4994 features, 4 samples
## element names: counts, normalizedCounts, offset
## protocolData: none
## phenoData
## sampleNames: mut_1 mut_2 wt_1 wt_2
## varLabels: conditions
## varMetadata: labelDescription
## featureData
## featureNames: YAL001C YAL002W ... YPR201W (4994
## total)
## fvarLabels: gc length
## fvarMetadata: labelDescription
## experimentData: use 'experimentData(object)'
## Annotation:
Note that the row names of counts
and featureData
must be the same; likewise for the row names of phenoData
and the column names of counts
. The expression values can be accessed with counts
, the lane information with pData
, and the feature information with fData
.
head(counts(data))
## mut_1 mut_2 wt_1 wt_2
## YAL001C 80 83 27 40
## YAL002W 33 38 53 66
## YAL003W 1887 1912 270 270
## YAL004W 90 110 276 295
## YAL005C 325 316 874 935
## YAL007C 27 30 19 24
pData(data)
## conditions
## mut_1 mut
## mut_2 mut
## wt_1 wt
## wt_2 wt
head(fData(data))
## gc length
## YAL001C 0.3712317 3483
## YAL002W 0.3717647 3825
## YAL003W 0.4460548 621
## YAL004W 0.4490741 648
## YAL005C 0.4406428 1929
## YAL007C 0.3703704 648
The SeqExpressionSet
class has two additional slots: normalizedCounts
and offset
(matrices of the same dimension as counts
), which may be used to store a matrix of normalized counts and of normalization offsets, respectively, to be used for subsequent analyses (see Section and the edgeR
vignette for details on the role of offsets). If not specified, the offset is initialized as a matrix of zeros.
head(offst(data))
## mut_1 mut_2 wt_1 wt_2
## YAL001C 0 0 0 0
## YAL002W 0 0 0 0
## YAL003W 0 0 0 0
## YAL004W 0 0 0 0
## YAL005C 0 0 0 0
## YAL007C 0 0 0 0
4.2Between-lane distribution of gene-level counts
One of the main considerations when dealing with gene-level counts is the difference in count distributions between lanes. The boxplot
method provides an easy way to produce boxplots of the logarithms of the gene counts in each lane.
boxplot(data,col=colors[1:4])
The MDPlot
method produces a mean-difference plot (MD-plot) of read counts for two lanes.
MDPlot(data,c(1,3))