Contents

1 Introduction

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)

2 Reading in unaligned and aligned read data

2.1 Unaligned 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.2 Aligned 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

3 Read-level EDA

3.1 Numbers 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 package leeBamViews. 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 signature c(x="BamFileList", y="FastqFileList"). See the manual pages for details.

3.2 Read 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.3 Individual 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.4 Read 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]])

4 Gene-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.1 Classes 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=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. As in the CountDataSet class, 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.2 Between-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))

4.3 Over-dispersion

Although the Poisson distribution is a natural and simple way to model count data, it has the limitation of assuming equality of the mean and variance. For this reason, the negative binomial distribution has been proposed as an alternative when the data show over-dispersion. The function meanVarPlot can be used to check whether the count data are over-dispersed (for the Poisson distribution, one would expect the points in the following Figures to be evenly scattered around the black line).

meanVarPlot(data[,1:2], log=TRUE, ylim=c(0,16))

meanVarPlot(data, log=TRUE, ylim=c(0,16))

Note that the mean-variance relationship should be examined within replicate lanes only (i.e., conditional on variables expected to contribute to differential expression). For the yeast dataset, it is not surprising to see no evidence of over-dispersion for the two mutant technical replicate lanes; likewise for the two wild-type lanes. However, one expects over-dispersion in the presence of biological variability, when considering at once all four mutant and wild-type lanes (Anders and Huber 2010,Bullard et al. (2010),Robinson, McCarthy, and Smyth (2010)).

4.4 Gene-specific effects on read counts

Several authors have reported selection biases related to sequence features such as gene length, GC-content, and mappability (Bullard et al. 2010,(???),Oshlack and Wakefield (2009),Risso et al. (2011)).

In the following figure, obtained using biasPlot, one can see the dependence of gene-level counts on GC-content. The same plot could be created for gene length or mappability instead of GC-content.

biasPlot(data, "gc", log=TRUE, ylim=c(1,5))

To show that GC-content dependence can bias differential expression analysis, one can produce stratified boxplots of the log-fold-change of read counts from two lanes using the biasBoxplot method. Again, the same type of plots can be created for gene length or mappability.

lfc <- log(counts(data)[,3]+0.1) - log(counts(data)[,1]+0.1)
biasBoxplot(lfc, fData(data)$gc)

5 Normalization

Following (Risso et al. 2011), we consider two main types of effects on gene-level counts: (1) within-lane gene-specific (and possibly lane-specific) effects, e.g., related to gene length or GC-content, and (2) effects related to between-lane distributional differences, e.g., sequencing depth. Accordingly, withinLaneNormalization and betweenLaneNormalization adjust for the first and second type of effects, respectively. We recommend to normalize for within-lane effects prior to between-lane normalization.

We implemented four within-lane normalization methods, namely: loess robust local regression of read counts (log) on a gene feature such as GC-content (loess), global-scaling between feature strata using the median (median), global-scaling between feature strata using the upper-quartile (upper), and full-quantile normalization between feature strata (full). For a discussion of these methods in context of GC-content normalization see (Risso et al. 2011).

dataWithin <- withinLaneNormalization(data,"gc", which="full")
dataNorm <- betweenLaneNormalization(dataWithin, which="full")

Regarding between-lane normalization, the package implements three of the methods introduced in (Bullard et al. 2010): global-scaling using the median (median), global-scaling using the upper-quartile (upper), and full-quantile normalization (full).

The next figure shows how after full-quantile within- and between-lane normalization, the GC-content bias is reduced and the distribution of the counts is the same in each lane.

biasPlot(dataNorm, "gc", log=TRUE, ylim=c(1,5))

boxplot(dataNorm, col=colors)

5.1 Offset

Some authors have argued that it is better to leave the count data unchanged to preserve their sampling properties and instead use an offset for normalization purposes in the statistical model for read counts (Anders and Huber 2010,(???),Robinson, McCarthy, and Smyth (2010)). This can be achieved easily using the argument offset in both normalization functions.

dataOffset <- withinLaneNormalization(data,"gc",
                                      which="full",offset=TRUE)
dataOffset <- betweenLaneNormalization(dataOffset,
                                       which="full",offset=TRUE)

Note that the dataOffset object will have both normalized counts and offset stored in their respective slots.

6 Differential expression analysis

One of the main applications of RNA-Seq is differential expression analysis. The normalized counts (or the original counts and the offset) obtained using the EDASeq package can be supplied to packages such as edgeR (Robinson, McCarthy, and Smyth 2010) or DESeq (Anders and Huber 2010) to find differentially expressed genes. This section should be considered only as an illustration of the compatibility of the results of EDASeq with two of the most widely used packages for differential expression; our aim is not to compare differential expression strategies (e.g., normalized counts vs. offset).

6.1 edgeR

We can perform a differential expression analysis with edgeR based on the original counts by passing an offset to the generalized linear model. For simplicity, we estimate a common dispersion parameter for all genes. See the edgeR vignette for details about how to perform a differential expression analysis using a gene-specific dispersion or more complex designs.

library(edgeR)
design <- model.matrix(~conditions, data=pData(dataOffset))
disp <- estimateGLMCommonDisp(counts(dataOffset),
                              design, offset=-offst(dataOffset))

fit <- glmFit(counts(dataOffset), design, disp, offset=-offst(dataOffset))

lrt <- glmLRT(fit, coef=2)
topTags(lrt)
## Coefficient:  conditionswt 
##               logFC   logCPM       LR        PValue
## YGL088W   -5.553631 27.34850 923.9754 6.024802e-203
## YPL198W   -7.425858 26.51938 899.4352 1.302000e-197
## YMR013W-A -6.236487 26.16157 775.4217 1.191073e-170
## YHR156C   -8.132766 25.10965 647.5159 7.750965e-143
## YLR106C   -4.162281 25.93346 612.7063 2.884815e-135
## YGL076C   -4.206669 26.10030 585.3729 2.543201e-129
## YOR040W   -4.974040 25.46016 567.4563 2.007391e-125
## YBL004W   -5.261041 25.03250 556.0820 5.983023e-123
## YAL003W   -2.991170 27.49022 517.4098 1.549064e-114
## YJL047C-A -4.346547 24.96907 416.7858  1.221958e-92
##                     FDR
## YGL088W   3.008786e-199
## YPL198W   3.251093e-194
## YMR013W-A 1.982739e-167
## YHR156C   9.677080e-140
## YLR106C   2.881353e-132
## YGL076C   2.116791e-126
## YOR040W   1.432130e-122
## YBL004W   3.734902e-120
## YAL003W   8.595581e-112
## YJL047C-A  6.102457e-90

6.2 DESeq

We can perform a differential expression analysis with DESeq based on the normalized counts by using the coerce method from the SeqExpressionSet class to the CountDataSet class of DESeq. When working with data that have been normalized for both within- and between-lane effects, we force the size factors to be one, since differences in lane sequencing depths have already been accounted for in our between-lane normalization. One could also consider only within-lane normalization and account for differences in sequencing depth by estimating the size factors using DESeq.

library(DESeq)
counts <- as(dataNorm, "CountDataSet")
sizeFactors(counts) <- rep(1,4)
counts <- estimateDispersions(counts)
res <- nbinomTest(counts, "wt", "mut")
head(res)
##        id baseMean baseMeanA baseMeanB foldChange
## 1 YAL001C    92.00      58.0     126.0  2.1724138
## 2 YAL002W    73.75     111.0      36.5  0.3288288
## 3 YAL003W   565.00     126.0    1004.0  7.9682540
## 4 YAL004W   100.00     133.0      67.0  0.5037594
## 5 YAL005C   259.25     330.5     188.0  0.5688351
## 6 YAL007C    32.25      37.0      27.5  0.7432432
##   log2FoldChange         pval         padj
## 1      1.1192989 3.415271e-04 3.732136e-03
## 2     -1.6045913 2.106666e-06 4.826006e-05
## 3      2.9942636 3.882861e-31 1.615917e-28
## 4     -0.9891932 1.181444e-03 1.022553e-02
## 5     -0.8139176 1.418425e-03 1.178638e-02
## 6     -0.4280937 3.346353e-01 5.776740e-01

7 Definitions and conventions

7.1 Rounding

After either within-lane or between-lane normalization, the expression values are not counts anymore. However, their distribution still shows some typical features of counts distribution (e.g., the variance depends on the mean). Hence, for most applications, it is useful to round the normalized values to recover count-like values, which we refer to as “pseudo-counts”.

By default, both withinLaneNormalization and betweenLaneNormalization round the normalized values to the closest integer. This behavior can be changed by specifying round=FALSE. This gives the user more flexibility and assures that rounding approximations do not affect subsequent computations (e.g., recovering the offset from the normalized counts).

7.2 Zero counts

To avoid problems in the computation of logarithms (e.g. in log-fold-changes), we add a small positive constant (namely \(0.1\)) to the counts. For instance, the log-fold-change between \(y_1\) and \(y_2\) is defined as \[\begin{equation*} \frac{\log(y_1 + 0.1)}{\log(y_2 + 0.1)}. \end{equation*}\]

7.3 Offset

We define an offset in the normalization as \[\begin{equation*} o = \log(y_{norm} + 0.1) - \log(y_{raw} + 0.1), \end{equation*}\] where \(y_{norm}\) and \(y_{raw}\) are the normalized and raw counts, respectively.

One can easily recover the normalized data from the raw counts and offset, as shown here:

dataNorm <- betweenLaneNormalization(data, round=FALSE, offset=TRUE)

norm1 <- normCounts(dataNorm)
norm2 <- exp(log(counts(dataNorm) + 0.1 ) + offst(dataNorm)) - 0.1

head(norm1 - norm2)
##                 mut_1         mut_2          wt_1
## YAL001C -1.421085e-14 -1.421085e-14  0.000000e+00
## YAL002W  3.552714e-15 -7.105427e-15 -2.131628e-14
## YAL003W  4.547474e-13 -2.273737e-13  1.136868e-13
## YAL004W -2.842171e-14 -1.421085e-14  5.684342e-14
## YAL005C -1.136868e-13  5.684342e-14 -3.410605e-13
## YAL007C  3.552714e-15  3.552714e-15 -3.552714e-15
##                  wt_2
## YAL001C  0.000000e+00
## YAL002W  0.000000e+00
## YAL003W -5.684342e-14
## YAL004W -5.684342e-14
## YAL005C -1.136868e-13
## YAL007C  3.552714e-15

Note that the small constant added in the definition of offset does not matter when pseudo-counts are considered, i.e.,

head(round(normCounts(dataNorm)) - round(counts(dataNorm) * exp(offst(dataNorm))))
##         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

We defined the offset as the log-ratio between normalized and raw counts. However, the edgeR functions expect as offset argument the log-ratio between raw and normalized counts. One must use -offst(offsetData) as the offset argument of edgeR.

8 Retrieving gene length and GC-content

Two essential features the gene-level EDA normalizes for are gene length and GC-content. As users might wish to automatically retrieve this information, we provide the function getGeneLengthAndGCContent. Given selected ENTREZ or ENSEMBL gene IDs and the organism under investigation, this can be done either based on BioMart (default) or using BioC annotation utilities.

getGeneLengthAndGCContent(id=c("ENSG00000012048", "ENSG00000139618"), org="hsa")

Accordingly, we can retrieve the precalculated yeast data that has been used throughout the vignette via

fData(data) <- getGeneLengthAndGCContent(featureNames(data),
                                              org="sacCer3", mode="org.db")

9 SessionInfo

sessionInfo()
## R version 3.6.0 (2019-04-26)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows Server 2012 R2 x64 (build 9600)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=C                          
## [2] LC_CTYPE=English_United States.1252   
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.1252    
## 
## attached base packages:
## [1] stats4    parallel  stats     graphics  grDevices
## [6] utils     datasets  methods   base     
## 
## other attached packages:
##  [1] DESeq_1.36.0                lattice_0.20-38            
##  [3] locfit_1.5-9.1              edgeR_3.26.0               
##  [5] limma_3.40.0                leeBamViews_1.19.0         
##  [7] BSgenome_1.52.0             rtracklayer_1.44.0         
##  [9] yeastRNASeq_0.21.0          EDASeq_2.18.0              
## [11] ShortRead_1.42.0            GenomicAlignments_1.20.0   
## [13] SummarizedExperiment_1.14.0 DelayedArray_0.10.0        
## [15] matrixStats_0.54.0          Rsamtools_2.0.0            
## [17] GenomicRanges_1.36.0        GenomeInfoDb_1.20.0        
## [19] Biostrings_2.52.0           XVector_0.24.0             
## [21] IRanges_2.18.0              S4Vectors_0.22.0           
## [23] BiocParallel_1.18.0         Biobase_2.44.0             
## [25] BiocGenerics_0.30.0         knitr_1.22                 
## [27] BiocStyle_2.12.0           
## 
## loaded via a namespace (and not attached):
##  [1] httr_1.4.0             bit64_0.9-7           
##  [3] splines_3.6.0          R.utils_2.8.0         
##  [5] assertthat_0.2.1       BiocManager_1.30.4    
##  [7] aroma.light_3.14.0     latticeExtra_0.6-28   
##  [9] blob_1.1.1             GenomeInfoDbData_1.2.1
## [11] yaml_2.2.0             progress_1.2.0        
## [13] RSQLite_2.1.1          digest_0.6.18         
## [15] RColorBrewer_1.1-2     htmltools_0.3.6       
## [17] Matrix_1.2-17          R.oo_1.22.0           
## [19] XML_3.98-1.19          pkgconfig_2.0.2       
## [21] biomaRt_2.40.0         genefilter_1.66.0     
## [23] bookdown_0.9           zlibbioc_1.30.0       
## [25] xtable_1.8-4           annotate_1.62.0       
## [27] GenomicFeatures_1.36.0 survival_2.44-1.1     
## [29] magrittr_1.5           crayon_1.3.4          
## [31] memoise_1.1.0          evaluate_0.13         
## [33] R.methodsS3_1.7.1      hwriter_1.3.2         
## [35] tools_3.6.0            prettyunits_1.0.2     
## [37] hms_0.4.2              stringr_1.4.0         
## [39] AnnotationDbi_1.46.0   compiler_3.6.0        
## [41] rlang_0.3.4            grid_3.6.0            
## [43] RCurl_1.95-4.12        bitops_1.0-6          
## [45] rmarkdown_1.12         DBI_1.0.0             
## [47] R6_2.4.0               bit_1.1-14            
## [49] KernSmooth_2.23-15     stringi_1.4.3         
## [51] Rcpp_1.0.1             geneplotter_1.62.0    
## [53] xfun_0.6

References

Anders, Simon, and Wolfgang Huber. 2010. “Differential Expression Analysis for Sequence Count Data.” Genome Biology 11 (10). BioMed Central:R106.

Bullard, James H, Elizabeth Purdom, Kasper D Hansen, and Sandrine Dudoit. 2010. “Evaluation of Statistical Methods for Normalization and Differential Expression in mRNA-Seq Experiments.” BMC Bioinformatics 11 (1). BioMed Central:94.

Langmead, Ben, Cole Trapnell, Mihai Pop, and Steven L Salzberg. 2009. “Ultrafast and Memory-Efficient Alignment of Short Dna Sequences to the Human Genome.” Genome Biology 10 (3). BioMed Central:R25.

Lee, Albert, Kasper Daniel Hansen, James Bullard, Sandrine Dudoit, and Gavin Sherlock. 2008. “Novel Low Abundance and Transient Rnas in Yeast Revealed by Tiling Microarrays and Ultra High–Throughput Sequencing Are Not Conserved Across Closely Related Yeast Species.” PLoS Genetics 4 (12). Public Library of Science:e1000299.

Oshlack, Alicia, and Matthew J Wakefield. 2009. “Transcript Length Bias in Rna-Seq Data Confounds Systems Biology.” Biology Direct 4 (1). BioMed Central:14.

Risso, Davide, Katja Schwartz, Gavin Sherlock, and Sandrine Dudoit. 2011. “GC-Content Normalization for Rna-Seq Data.” BMC Bioinformatics 12 (1). BioMed Central:480.

Robinson, Mark D, Davis J McCarthy, and Gordon K Smyth. 2010. “EdgeR: A Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data.” Bioinformatics 26 (1). Oxford University Press:139–40.

“Saccharomyces Genome Database.” n.d. http://www.yeastgenome.org.