High-throughput sequencing technologies are now in common use in biology. ese 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] [2] [3] [4] . However, the datasets produced are large and complex and inter pretation is not straight forward. As with any high-throughput technology, analysis methodology is critical to interpreting the data, and RNA-seq analysis procedures are continuing to evolve. erefore, 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] [11] [12] , RNA editing [13] and novel transcripts [11, 14] . However, the primary objective of many biological studies is gene expression profiling between samples. us, in this review we focus on the methodologies available to detect differences in gene level expression between samples. is sort of analy sis 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 RNAseq 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 highthroughput platform, such as the Illumina GA/ HiSeq, SOLiD or Roche 454 [16] . is 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. e 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  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] [20] [21] [22] or the Burrows Wheeler transform (BWT) [23] [24] [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] [27] [28] [29] [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 multimapping, 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 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.
However, the Poisson distribution assumption on the number of reads cannot capture the variability of the RNA-Seq data well, leading to high-false positive rates in the detection of DEGs REF .
From RNA-seq reads to differential expression results
Genome Biology
Genome Biology
Biology Medicine