From RNA-seq reads to differential expression results.

Many methods and tools are available for preprocessing high-throughput RNA sequencing data and detecting differential expression.

High-throughput sequencing technologies are now in common use in biology. These technologies produce millions of short sequence reads and are routinely being applied to genomes, epigenomes and transcriptomes. Sequencing steady-state RNA in a sample, known as RNA-seq, is free from many of the limitations of previous technologies, such as the dependence on prior knowledge of the organism, as required for microarrays and PCR (see Box 1: Comparisons of microarrays and sequencing for gene expression analysis). In addition, RNA-seq promises to unravel previously inaccessible complexities in the transcriptome, such as allele-specific expression and novel promoters and isoforms [1-4]. However, the datasets produced are large and complex and interpretation is not straightforward. As with any high-throughput technology, analysis methodology is critical to interpreting the data, and RNA-seq analysis procedures are continuing to evolve. Therefore, it is timely to review currently available data analysis methods and comment on future research directions.

Making sense of RNA-seq data depends on the scientific question of interest. For example, determining differences in allele-specific expression requires accurate determination of the prevalence of transcribed single nucleotide polymorphisms (SNPs) [5]. Alternatively, fusion genes or aberrations in cancer samples can be detected by finding novel transcripts in RNA-seq data [6,7]. In the past year, several methods have emerged that use RNA-seq data for abundance estimation [8,9], detection of alternative splicing [10-12], RNA editing [13] and novel transcripts [11,14]. However, the primary objective of many biological studies is gene expression profiling between samples. Thus, in this review we focus on the methodologies available to detect differences in gene level expression between samples. This sort of analysis is particularly relevant for controlled experiments comparing expression in wild-type and mutant strains of the same tissue, comparing treated versus untreated cells, cancer versus normal, and so on. For example, comparison of expression changes between the cultured pathogen Acinetobacter baumannii and the pathogen grown in the presence of ethanol - which is known to increase virulence - revealed 49 differentially expressed genes belonging to a range of functional categories [15]. Here we outline the processing pipeline used for detecting differential expression (DE) in RNA-seq and examine the available methods and open-source software tools to perform the analysis. We also highlight several areas that require further research.

Most RNA-seq experiments take a sample of purified RNA, shear it, convert it to cDNA and sequence on a high-throughput platform, such as the Illumina GA/HiSeq, SOLiD or Roche 454 [16]. This process generates millions of short (25 to 300 bp) reads taken from one end of the cDNA fragments. A common variant on this process is to generate short reads from both ends of each cDNA fragment, known as 'paired-end' reads. The platforms differ substantially in their chemistry and processing steps, but regardless of the precise details, the raw data consist of a long list of short sequences with associated quality scores; these form the entry point for this review.

An overview of the typical RNA-seq pipeline for DE analysis is outlined in Figure 1. First, reads are mapped to the genome or transcriptome. Second, mapped reads for each sample are assembled into gene-level, exon-level or transcript-level expression summaries, depending on the aims of the experiment. Next, the summarized data are normalized in concert with the statistical testing of DE, leading to a ranked list of genes with associated P-values and fold changes. Finally, biological insight from these lists can be gained by performing systems biology approaches, similar to those performed on microarray experiments. We critique below the currently available methodologies for each of these steps for RNA-seq data analysis. Rather than providing a complete list of all available tools, we focus on examples of commonly used open-source software that illustrate the methodology (Table 1). For a complete list of RNA-seq analysis software, see [17,18].

Figure 1

Overview of the RNA-seq analysis pipeline for detecting differential expression. The steps in the pipeline are in red boxes; the methodological components of the pipeline are shown in blue boxes and bold text; software examples and methods for each step (a non-exhaustive list) are shown by regular text in blue boxes. References for the tools and methods shown are listed in Table 1. First, reads are mapped to the reference genome or transcriptome (using junction libraries to map reads that cross exon boundaries); mapped reads are assembled into expression summaries (tables of counts, showing how may reads are in coding region, exon, gene or junction); the data are normalized; statistical testing of differential expression (DE) is performed, producing and a list of genes with associated P-values and fold changes. Systems biology approaches can then be used to gain biological insights from these lists.

Table 1

Software methods and tools for differential expression analysis of RNA-seq

Analysis step	Method	Implementation	References
Mapping	General aligner	GMAP/GSNAP	[91]
		BFAST	[20]
		BOWTIE	[25]
		CloudBurst	[92]
		GNUmap	[93]
		MAQ/BWA	[23]
		PerM	[19]
		RazerS	[94]
		Mrfast/mrsfast	[22]
		SOAP/SOAP2	[24,95]
		SHRiMP	[21]
	De novo annotator	QPALMA/GenomeMapper/PALMapper	[37]
		SpliceMap	[96]
		SOAPals	[95]
		G-Mo.R-Se	[97]
		TopHat	[40]
		SplitSeek	[36]
	De novo transcript assembler	Oases	[98]
		MIRA	[99]
Summarization	Isoform-based	Cufflinks	[11]
		ALEXA-seq	[10]
	Gene-based	Count exons only	For example, [34,45]
		Exon junction libraries	[34,44]
Normalization	Library size		For example, [34]
	RPKM	ERANGE	[32]
	TMM	edgeR	[48]
	Upper quartile	Myrna	[45,47]
Differential expression	Poisson GLM	DEGseq	[100]
		Myrna	[47]
	Negative binomial	edgeR	[57]
		DESeq	[46]
		baySeq	[58]
Systems biology	Gene Ontology analysis	GOseq	[68]

Abbreviations: GLM, generalized linear model; RPKM, reads per kilobase of exon model per million mapped reads; TMM, trimmed mean of M-values.

Mapping

To use RNA-seq data to compare expression between samples, it is necessary to turn millions of short reads into a quantification of expression. The first step in this procedure is the read mapping or alignment. At its simplest, the task of mapping is to find the unique location where a short read is identical to the reference. However, in reality the reference is never a perfect representation of the actual biological source of RNA being sequenced. In addition to sample-specific attributes such as SNPs and indels (insertions or deletions), there is also the consideration that the reads arise from a spliced transcriptome rather than a genome. Furthermore, short reads can sometimes align perfectly to multiple locations and can contain sequencing errors that have to be accounted for. Therefore, the real task is to find the location where each short read best matches the reference, while allowing for errors and structural variation.

Although research into how best to align reads to a reference is ongoing, all solutions by necessity involve some compromise between the computational requirements of the algorithm and the fuzziness allowed in matching to the reference. Almost all short read aligners use a strategy of a first pass 'heuristic' match, which quickly finds a reduced list of possible locations, followed by thorough evaluation of all candidate alignments by a complex 'local alignment' algorithm. Without this initial heuristic search to reduce the number of potential alignment locations, performing local alignment of millions of short reads would be computationally impossible on current hardware.

Current aligners enable fast heuristic matching by using either hash tables [19-22] or the Burrows Wheeler transform (BWT) [23-25]. Hash-table aligners have the advantage of being easily extendable to detect complicated differences between read and reference, at the cost of ever increasing computational requirements. Alternatively, BWT-based aligners can map reads that closely match the reference very efficiently but are prohibitively slow once more complex misalignments are considered. A detailed explanation of these techniques is beyond the scope of this review, but can be found in [23,26-30].

Aligners also differ in how they handle 'multimaps' (reads that map equally well to several locations). Most aligners either discard multimaps [25], allocate them randomly [29] or allocate them on the basis of an estimate of local coverage [31,32], although a statistical method incorporating alignment scores has also been proposed [33]. Paired-end reads reduce the problem of multi-mapping, as both ends of the cDNA fragment from which the short reads were generated should map nearby on the transcriptome, allowing the ambiguity of multimaps to be resolved in most circumstances.

When considering reads from genomic DNA, mapping to a relevant reference genome is all that is needed. However, RNA-seq is sequencing fragments of the transcriptome. This difference is dealt with in several ways. Given that the transcriptome is 'built from' the genome, the most commonly used approach (at least initially) is to use the genome itself as the reference. This has the benefit of being easy and not biased towards any known annotation. However, reads that span exon boundaries will not map to this reference. Thus, using the genome as a reference will give greater coverage (at the same true expression level) to transcripts with fewer exons, as they will contain fewer exon junctions. Longer reads are more likely to cross exon boundaries, thus causing the fraction of junction reads to increase [2].

In order to account for junction reads, it is common practice to build exon junction libraries in which reference sequences are constructed using boundaries between annotated exons [2,32,34,35]. To map reads that cross exon boundaries without relying on existing annotations, it is possible to use the dataset itself to detect splice junctions de novo [36-41]. Another option is the de novo assembly of the transcriptome, for use as a reference, using genome assembly tools [42,43]. All de novo methods can identify novel transcripts and may be the only option for organisms for which no genomic reference or annotation is available. However, de novo methods are computationally intensive and may require long, paired-end reads and high levels of coverage to work reliably. For example, Trapnell et al. [11] used over 430 million paired-end reads for de novo assembly of the mouse myoblast transcriptome in order to quantify expression during cell differentiation.

A commonly used approach for transcriptome mapping is to progressively increase the complexity of the mapping strategy to handle the unaligned reads [44]. For example, in a large study investigating expression variation in 69 Nigerian HapMap samples, Pickrell et al. [35] found that for 46 bp Illumina reads, 87% mapped to the reference genome with two mismatches using MAQ (a hash-table-based aligner) [23]. An additional 7% could be mapped to an exon-exon junction library, constructed from all possible combinations of Ensembl exons. The remaining unmapped reads were examined for evidence of the sequencer having erroneously sequenced the poly(A) tail. If a read began or ended with at least four As or Ts, these bases were trimmed and the rest of the read was mapped to the reference, resulting in a further 0.005% of reads being mapped. This large dataset enabled the annotation of over 100 new exons and identified more than a thousand genes in which genetic variation influences overall expression levels or splicing. This would not have been possible without a method for handling reads that cross exon boundaries.

Summarizing mapped reads

Having obtained genomic locations for as many reads as possible, the next task is to summarize and aggregate reads over some biologically meaningful unit, such as exons, transcripts or genes. The simplest and most common approach counts the number of reads overlapping the exons in a gene (for example, [32,34,45]). However, a significant proportion of reads map to genomic regions outside annotated exons, even in well-annotated organisms, such as mouse and human. For example, Pickrell et al. [35] found that about 15% of mapped reads were located outside annotated exons for their Nigerian HapMap samples and these extra-exonic reads were more likely to be cell-type-specific exons. Similarly, Figure 2a shows an example of transcription occurring outside annotated exons in the RNA-binding protein 39 (RBM39) gene in LNCaP prostate cancer cells. Reads from other normal tissue cell types are more limited to known exons, but also show evidence for transcription outside of known exons.

Figure 2

Summarizing mapped reads into a gene level count. (a) Mapped reads from a small region of the RNA-binding protein 39 (RBM39) gene are shown for LNCaP prostate cancer cells [90], human liver and human testis from the UCSC track. The three rows of RNA-seq data (blue and black graphs) are shown as a 'pileup track', where the y-axis at each location measures the number of mapped reads that overlap that location. Also shown are the genomic coordinates, gene model (labeled RBM39; blue boxes indicate exons) and conservation score across vertebrates. It is clear that many reads originate from regions with no known exons. (b) A schematic of a genomic region and reads that might arise from it. Reads are color-coded by the genomic feature from which they originate. Different summarization strategies will result in the inclusion or exclusion of different sets of reads in the table of counts. For example, including only reads coming from known exons will exclude the intronic reads (green) from contributing to the results. Splice junctions are listed as a separate class to emphasize both the potential ambiguity in their assignment (such as which exon should a junction read be assigned to) and the possibility that many of these reads may not be mapped because they are harder to map than continuous reads. CDS, coding sequence.

One alternative summarization is to include reads along the whole length of the gene and thereby incorporate reads from 'introns'. This will include unannotated exons in the summary and account for poorly annotated or variable exon boundaries. However, including introns might also capture overlapping transcripts, which share a genomic location but originate from different genes. There are many other possible variations that could be used for summarization, such as including only reads that map to coding sequence or summarizing from de novo predicted exons [40]. Junction reads can also be added into the gene summary count or be used to model the abundance of splicing isoforms [11]. These different possibilities are illustrated schematically in Figure 2b. With these options, the choice of summarization has the potential to change the count for each gene as substantially as, or more substantially than, the choice of mapping strategy. Despite this, little research has been carried out on which summarization method is the most appropriate for DE detection.

Normalization

Normalization enables accurate comparisons of expression levels between and within samples [2,32,34]. It has been shown that normalization is an essential step in the analysis of DE from RNA-seq data [45-48]. Normalization methods differ for between- and within-library comparisons.

Within-library normalization allows quantification of expression levels of each gene relative to other genes in the sample. Because longer transcripts have higher read counts (at the same expression level), a common method for within-library normalization is to divide the summarized counts by the length of the gene [32,34]. The widely used RPKM (reads per kilobase of exon model per million mapped reads) accounts for both library size and gene length effects in within-sample comparisons. To validate this approach, Mortazavi et al. [32] introduced several Arabidopsis RNAs into their mouse tissue samples, across a range of gene lengths and expression levels. These non-native RNAs are known as 'spike-ins' and demonstrated that RPKM gives accurate comparisons of expression levels between genes. However, it has been shown that read coverage along expressed transcripts can be non-uniform because of sequence content [49] and RNA preparation methods, such as random hexamer priming [50]. Incorporating this understanding into the within-library normalization method may improve the ability to compare expression levels. Using RNA-seq data to estimate the absolute number of transcripts in a sample is possible, but it requires RNA standards and additional information, such as the total number of cells from which RNA is extracted and RNA preparation yields [32].

When testing individual genes for DE between samples, technical biases, such as gene length and nucleotide composition, will mainly cancel out because the underlying sequence used for summarization is the same between samples. However, between-sample normalization is still essential for comparing counts from different libraries relative to each other. The simplest and most commonly used normalization adjusts by the total number of reads in the library [34,51], accounting for the fact that more reads will be assigned to each gene if a sample is sequenced to a greater depth. However, it has been shown that more sophisticated normalization is required to account for composition effects [48], or for the fact that a small number of highly expressed genes can consume a significant amount of the total sequence [45]. To account for these features, scaling factors can be estimated from the data and used within the statistical models that test for DE [45,46,48]. Scaling factors have the advantage that the raw count data are preserved for subsequent analysis. Alternatively, quantile normalization and a method using matching power law distributions [52,53] have also been proposed for between-sample normalization of RNA-seq. The non-linearity of both of these transformations removes the count nature of the data, making it unclear how to appropriately test for DE. So far, quantile normalization does not seem to improve DE detection to the same extent as an appropriate scaling factor [45] and it is not clear that the power law distribution applies to all datasets [48].

Differential expression

The goal of a DE analysis is to highlight genes that have changed significantly in abundance across experimental conditions. In general, this means taking a table of summarized count data for each library and performing statistical testing between samples of interest.

Many methods have been developed for the analysis of differential expression using microarray data. However, RNA-seq gives a discrete measurement for each gene whereas microarray intensities have a continuous intensity distribution. Although microarray intensities are typically log-transformed and analyzed as normally distributed random variables, transformation of count data is not well approximated by continuous distributions, especially in the lower count range and for small samples. Therefore, statistical models appropriate for count data are vital to extracting the most information from RNA-seq data.

In general, the Poisson distribution forms the basis for modeling RNA-seq count data. In an early RNA-seq study using a single source of RNA, sequenced on multiple lanes of an Illumina GA sequencer, goodness-of-fit statistics suggested that the distribution of counts across lanes for the majority of genes was indeed Poisson distributed [34]. This has been independently confirmed using a technical experiment [45] and software tools are readily available to perform these analyses [54]. However, biological variability is not captured well by the Poisson assumption [47,51]. Hence, Poisson-based analyses for datasets with biological replicates will be prone to high false positive rates resulting from the underestimation of sampling error [46,47,55]. Despite the low background and high sensitivity of the RNA-seq platform, designing experiments with biological replication is still critical for identifying changes in RNA abundance that generalize to the population being sampled. Design of RNA-seq experiments in general, including the fundamental considerations of blocking, randomization and replication, has recently been discussed in depth [56].

In order to account for biological variability, methods that have been developed for serial analysis of gene expression (SAGE) data have recently been applied to RNA-seq data [57]. The major difference between SAGE and RNA-seq data is the scale of the datasets. To account for biological variability, the negative binomial distribution has been used as a natural extension of the Poisson distribution, requiring an additional dispersion parameter to be estimated. A few variations of negative-binomial-based DE analysis of count data have emerged, including common dispersion models [55], sharing information over all genes using weighted likelihood [51], empirical estimation of the mean-variance relationship [46] and an empirical Bayesian implementation using equivalence classes [58]. Extensions to the Poisson model to include overdispersion have also been proposed, through the generalized Poisson distribution [59] or a two-stage Poisson model, which tests for differential expression in two modes depending on the evidence for overdispersion in the data [60]. Several tools for either simultaneous transcript discovery and quantification [11] or alternative isoform expression analysis [10] also perform DE analysis. However, it is worth noting that these methods use either the Poisson distribution or Fisher's exact test, neither of which explicitly deal with the biological variation discussed above.

Many of the current strategies for DE analysis of count data are limited to simple experimental designs, such as pairwise or multiple group comparisons. To the best of our knowledge, no general methods have been proposed for the analysis of more complex designs, such as paired samples or time course experiments, in the context of RNA-seq data. In the absence of such methods, researchers have transformed their count data and used tools appropriate for continuous data [31,47,61]. Generalized linear models provide the logical extension to the count models presented above, and clever strategies to share information over all genes will need to be developed; software tools now provide these methods (such as edgeR [57]). Furthermore, the methods discussed above are predominantly aimed at summarizing expression levels at which annotation exists. Methods, such as the maximum mean discrepancy test [62], have recently been proposed to detect DE in an untargeted manner.

Systems biology: going beyond gene lists

In many cases, creating lists of DE genes is not the final step of the analysis; further biological insight into an experimental system can be gained by looking at the expression changes of sets of genes. Many tools focusing on gene set testing, network inference and knowledge databases have been designed for analyzing lists of DE genes from microarray datasets [63-65]. However, RNA-seq is affected by biases not present in microarray data. For example, gene length bias is an issue in RNA-seq data, in which longer genes have higher counts (at the same expression level) [66]. This results in greater statistical power to detect DE for long and highly expressed genes. These biases can dramatically affect the results of downstream analyses, such as testing Gene Ontology (GO) terms for enrichment among DE genes [66,67]. In order to enable gene set analyses, Bullard et al. [45] suggested modifying a DE t-statistic by dividing by the square root of gene length to minimize the effect of length bias on DE. Alternatively, GO-seq is an approach developed specifically for RNA-seq data that can incorporate length or total count bias into gene set tests [68]. As the understanding of biases in RNA-seq data grows, systems biology tools that incorporate this understanding will be critical to extracting biological insight.

There is wide scope for integrating the results of RNA-seq data with other sources of biological data to establish a more complete picture of gene regulation [69]. For example, RNA-seq has been used in conjunction with genotyping data to identify genetic loci responsible for variation in gene expression between individuals (expression quantitative trait loci or eQTLs) [35,70]. Furthermore, integration of expression data with transcription factor binding, RNA interference, histone modification and DNA methylation information has the potential for greater understanding of a variety of regulatory mechanisms. A few reports of these 'integrative' analyses have emerged recently [71-73]. For example, Lister and coauthors [71] highlighted a striking difference in the correlations of RNA-seq expression with CG and non-CG methylation levels in gene bodies. Similarly, combinations of sequencing-based datasets are beginning to provide insights into the mono-allelic associations between expression, histone modifications and DNA methylation [74].

Outlook

In this review, we have outlined the major steps in processing the millions of short reads produced by RNA-seq into an analysis of DE between samples. In brief, the process is to map and summarize short read sequences, then normalize between samples and perform a statistical test of DE. Further biological insight can be gained by looking for patterns of expression changes within sets of genes and integrating the RNA-seq data with data from other sources.

Although many parts of this pipeline have been the focus of extensive research, there are still areas that offer the possibility of further refinements. So far, there has been little work researching which summarization metric is best suited to finding DE between samples. There is also scope for expanding existing statistical methods for DE detection to enable the analysis of more complex experimental designs. Moreover, the relative merits of the many approaches now available deserve further study, in terms of their flexibility to analyze various study designs, their performance in small and large studies, dependence on sequencing depth and the accuracy of the assumptions (such as mean-variance relationships) that are imposed. Furthermore, although there are many examples of using RNA-seq for the detection of alternative splicing, there is scope to extend current methods to detect differences in gene isoform preference [10,11] when biological variability is prominent, perhaps using the count-based statistical methods mentioned above.

Given that there are substantial differences in the protocols that generate short reads, it will be important to formally compare RNA-seq platforms and the relative merits of the many data analysis methodologies. Such investigations may reveal benefits of platform-specific DE analysis methods and will also facilitate greater data integration. As the field is still relatively young, we expect many new methods and tools for the analysis of RNA-seq data to emerge in the near future.

Box 1: Comparisons of microarrays and sequencing for gene expression analysis

Several comparisons of RNA-seq and microarray data have now been made. These include proof-of-principle demonstrations of the sequencing platform [2,31,32], dedicated comparison studies [34,75-77] and analysis methodology development [10]. The results are unanimous: sequencing has higher sensitivity and dynamic range, coupled with lower technical variation. Furthermore, comparisons have highlighted strong concordance between microarrays and sequencing in measures of both absolute and differential expression. Nevertheless, microarrays have been, and continue to be, highly successful in interrogating the transcriptome in many biological settings. Examples include defining the cell of origin for breast cancer subtypes [78] and investigating the effect of evolution on gene expression in Drosophila [79].

Microarrays and sequencing each have their own specific biases that can affect the ability of a platform to measure DE. It is well known that cross-hybridization of microarray probes affects expression measures in a non-uniform way [80,81] and sequence content influences measured probe intensities [82]. Meanwhile, several studies have observed a GC bias in RNA-seq data [45] and RNA-seq can suffer from mapping ambiguity for paralogous sequences. Furthermore, there is a higher statistical power to detect changes at higher counts (for example, a twofold difference of 200 reads to 100 reads is more statistically significant than 20 reads to 10, under the null hypothesis of no difference); this bias typically manifests in RNA-seq as an association between DE and gene length, an effect not present in microarray data [66,68]. Other studies indicate that specific sequencing protocols produce biases in the generated reads, which can be related to the sequence composition and distance along the transcript [49,50,83,84]. For example, library preparation for small RNAs has been found to strongly affect the set of observed sequences [85]. Furthermore, transcriptome assembly approaches are necessarily biased by expression level because less information is available for genes expressed at a low level [11,14]. Many of these biases are still being explored and clever statistical methods that harness this knowledge may be able to provide improvements on existing methods.

In addition to the larger dynamic range and sensitivity of RNA-seq, several additional factors have contributed to the rapid uptake of sequencing for differential expression analysis. First, microarrays are simply not available for many non-model organisms (for example, Affymetrix offers microarrays for approximately 30 species [86]). By contrast, genomes and sequence information are readily available for thousands of species [87]. Moreover, even when genomes are not available, RNA-seq can still be performed and the transcriptome can still be interrogated (for instance, a recent study used RNA-seq to investigate the cell origin of the Tasmanian Devil facial tumor [88]). Second, sequencing gives unprecedented detail about transcriptional features that arrays cannot, such as novel transcribed regions, allele-specific expression, RNA editing and a comprehensive capability to capture alternative splicing. For example, a recent RNA-seq study [11] was able to show several examples of isoform switching during cell differentiation, and RNA-seq was used to show parent-of-origin expression in mouse brain [5].

Sequencing is not without its challenges, of course. The cost of the platform may be limiting for some studies. However, with the expansion in total sequencing capacity and the ability to multiplex, the cost per sample to generate sufficient sequence depth will soon be comparable to that of microarrays. However, the cost of informatics to house, process and analyze the data is substantial [89]. Researchers with limited access to computing staff and resources may elect to use microarrays because data analysis procedures are relatively mature. Finally, it is clear that data analysis methodologies for sequencing data will continue to evolve for some time yet.

Author information

All authors contributed equally to this review.