R包:EDASeq–RNA-Seq 的数据分析和标准化

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 BamFileListfrom 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)

R包:EDASeq--RNA-Seq 的数据分析和标准化

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))

R包:EDASeq--RNA-Seq 的数据分析和标准化

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)

R包:EDASeq--RNA-Seq 的数据分析和标准化

barplot(bfs[[1]],las=2)

R包:EDASeq--RNA-Seq 的数据分析和标准化

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]])

R包:EDASeq--RNA-Seq 的数据分析和标准化

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])

R包:EDASeq--RNA-Seq 的数据分析和标准化

The MDPlot method produces a mean-difference plot (MD-plot) of read counts for two lanes.

MDPlot(data,c(1,3))