RNA-seqlopedia

When designing an RNA-seq experiment researchers are faced with choosing between many experimental options, and decisions must be made at each step of the process. In some cases the choice one makes will have little impact on the quality of the experimental data. However, in other cases inappropriate decisions can result in spending a great deal of time and money only to obtain nearly useless data. It is therefore imperative that one carefully considers these choices before launching into a time consuming and costly RNA-seq experiment.

This stage of the experiment is arguably the most important but often receives the least amount of attention. Before one even begins an RNA-seq experiment they should have an understanding of each step and should put together a carefully devised experimental plan.

The types of information that can be gained from RNA-seq can be divided into two broad categories: qualitative and quantitative.

• Qualitative data includes identifying expressed transcripts, and identifying exon/intron boundaries, transcriptional start sites (TSS), and poly-A sites. Here, we will refer to this type of information as "annotation".
• Quantitative data includes measuring differences in expression, alternative splicing, alternative TSS, and alternative polyadenylation between two or more treatments or groups. Here we focus specifically on experiments to measure differential gene expression (DGE).

These two types of data are related and to some degree are inseparable. For instance, high quality annotation can lead to greater accuracy in mapping which can lead to higher quality measures of gene expression. Similarly, accurate annotation requires some degree of quantification. For instance, the simple existence of one or two reads that map to an intergenic region is not strong evidence for the existence of a gene at that locus.

It is important to realize that experiments that are designed to measure quantitative changes have requirements that differ from those that are designed for annotation, and it would be difficult and expensive to design a single experiment that satisfies the requirements of both. It is therefore important to first identify the primary goal of an RNA-seq experiment so that the experiment can be designed to maximize the objective. Below, we highlight some of the important distinctions between these two objectives with regards to RNA-seq. Further guidelines can be found in Table 1.1.

The goal of annotation is to identify genes and genic architecture by analyzing short-sequence reads derived from expressed RNA. Thus, the most important parameter is that the sequence reads evenly cover each transcript, including both ends.

Coverage is primarily a function of the method used to prepare the library. Each library prep methods suffer from specific types of biases. For instance, oligo-dT priming for first-strand synthesis can be used to accurately annotate 3′-ends but often fails to evenly cover the 5′-end. Random priming suffers from uneven coverage due to sequence/structure priming biases, and can result in poor coverage of either of the ends. Some library preparation methods have been specifically designed to provide even transcript coverage, such as the SMARTer system from Clontech. Currently few publications have been devoted to evaluating the coverage and biases of different library prep methods; however, in one study (Levin et al 2010) it was observed that the dUTP-stranded method gave the best overall performance.

Regardless of the method used to prepare the libraries, each method still results in uneven coverage of individual transcripts. In order to get sequence reads that span the entire transcript more reads are required (i.e. deeper sequencing). Since the primary goal of annotation is to characterize what is present rather than how much of it is present annotation experiments can greatly benefit from enrichment procedures such as DSN normalization which reduces the proportion of reads that come from highly expressed RNAs, effectively increasing the reads that come from rarer transcripts without requiring deeper sequencing.

Two of the biggest challenges related to using short-read sequencing for annotation are that the sequence reads are much shorter than the biological transcripts and that in many cases a reference genome is unavailable. Both of these issues create problems for mapping and de novo transcript assembly. These problems can be partially alleviated by using sequencing strategies that produce longer sequence reads, paired-end reads, and stranded library preparation methods. Third generation sequencing platforms, such as PacBio, have much longer reads, which also can overcome these challenges, albeit at a greater cost than short-read sequencing.

The primary goal of DGE is to quantitatively measure differences in the levels of transcripts between two or more treatments or groups. At one level this is as simple as comparing the counts for reads that come from each transcript. However, simple read counts alone are neither sufficient for identifying differentially expressed genes nor for quantifying the degree of differences. This is largely because there is no way to know how accurately such counts reflect what was actually present in the original sample, and there is no indication of how widely expression normally varies in a given group or without treatment. Therefore in order to appropriately interpret differences in read counts one must also obtain information regarding the variances that are associated with these numbers. Several types of variances contribute to RNA-seq data including:

1. Sampling variance: Even though a sequencing run is capable of producing millions of reads, these represent only a small fraction of the nucleic acid that is actually present in the library. There is therefore a built in sampling variance.
2. Technical variance: Library preparation and sequencing procedures involve a series of complex chemical reactions which all contribute to between sample variance.
3. Biological variance: Ultimately we are interested in measuring differences between different biological systems. Biological systems are inherently complex and very sensitive to perturbations. Thus, even in the absence of sampling and technical variance biological variance will always exist. Specifically the variance we are referring to is the nascent variance that is present within a treatment or control group.

You have probably noticed that these key considerations (especially target transcript attributes and effect sizes) seem likely to vary widely on a case-by-case basis. Factors like the nature of the treatment(s), the type of sample, or the ability to control the environment of subjects, for example, really differ from experiment to experiment. This means that published data from some other experimental system, or even naïvely simulated data, may not be informative enough for critical decisions up front. In light of this conundrum, our single best recommendation is to first sequence a subset of replicates from a larger experiment, get a handle on coverage levels across genes and effect sizes of treatments, and then use this information to estimate the level of replication necessary to reach a given level of statistical power for a given fraction of the transcriptome. One strategy, provided the experimental setup is not too costly, might be to run the desired experiment with a substantial level of replication (e.g. 10-20 individuals per treatment or combination of treatments), but only actually sequence 3 or 4 samples at moderate to high coverage (e.g. 10 million reads per individual) for each treatment level. This way you would have enough preliminary information to conduct a reliable power analysis, but still have plenty of remaining samples "in the bank" for additional sequencing.

Fortunately, the Marth Lab at Boston College has produced a handy, flexible tool for estimating power of RNA-seq study designs, given a suite of user-supplied parameters. This tool, called "Scotty" (Busby et al. 2013), is especially effective when fed user-supplied pilot data, and will output a range of sample size/coverage configurations acceptable for specified power and cost constraints. The tool may be found here. Although this tool is currently designed to accommodate a simple two-treatment experimental design, the framework is in principle extendable to more complicated designs.

Since the goal of RNA-seq is to characterize the transcriptome the first step naturally involves isolating and purifying cellular RNAs. Isolation and purification of RNA typically involves disrupting cells in the presence of detergents and chaotropic agents. Depending upon the starting material mechanical disruption is also recommended. After homogenization, RNA can be recovered and purified from the total cell lysate using either liquid-liquid partitioning or solid-phase extraction. Typically the total RNA is then enriched for messenger RNA (mRNA). This can be done by either directly selecting mRNA or by selectively removing ribosomal RNA (rRNA). To make the RNA suitable for RNA-seq it is typically fragmented and then the quality and fragmentation are assessed.

The success of RNA-seq experiments is highly dependent upon recovering pure and intact RNA. Because RNA is more labile than DNA and RNases are very stable enzymes extra care should be taken when purifying and working with RNA. Some general considerations are:

Other information about working with RNA can be found at: Working with RNA (Ambion).

The overall goal of RNA-seq is to characterize the transcriptome at the moment the cells were harvested; thus care should be taken to limit degradation of RNAs.

Two of the primary challenges of the isolation procedure are (1) to separate the RNA from other cellular materials and (2) to preserve the integrity of the RNA that is extracted. This is complicated by the fact that DNA and RNA are chemically very similar to each other and that RNases are notoriously stable proteins that can still be active under conditions that denature other proteins.

It is often necessary to enrich specific classes of RNA in the sample to be sequenced. Total RNA recovered using the procedures described above consists of >80% ribosomal RNA (rRNA) (Raz et al. 2011), so if rRNA were not removed, the majority of the final sequence reads would be from rRNA.

The four methods that are commonly used to enrich specific classes of RNAs are:

1. Selection of target RNAs via hybridization.
2. Removal of non-target RNAs via hybridization.
3. Copy-number normalization via duplex-specific nuclease digestion.
4. Target enrichment via size-selection

Most current sequencing platforms are capable of providing only relatively short sequence reads (~40-400bp depending upon the platform). Therefore, most protocols incorporate a fragmentation step to improve sequence coverage over the transcriptome. However, protocols differ as to when the fragmentation is performed. Most of the original protocols fragmented cDNA; however, fragmentation of the RNA (before converting it to cDNA) is becoming more popular. Four different methods are commonly used to fragment RNA: enzymatic, metal ion, heat, and sonication.

Although all of these methods are commonly used to fragment RNA, little information has been published comparing their relative performance. The goal is to produce a population of RNA fragments that are on average about 200 bp. The effectiveness of the fragmentation reaction should be assessed by evaluating the RNA on a gel or by using an Agilent Bioanalyzer. An example of fragmented RNA can be seen in Figure 2.3. This example and guidelines for fragmenting RNA can be found at the U. of Texas at Austin FGRS site.

Figure 2.3

After obtaining an RNA preparation that is suitable for RNA-seq the RNA must be converted to double-stranded complementary DNA (cDNA). Currently available sequencing technologies require a DNA template with platform-specific "adaptor" sequences at either end of each molecule. Generating the cDNA, adding the adaptors, and amplifying the DNA for sequencing round out the process of "library preparation". Competing RNA-seq procedures diverge during library preparation, and the protocol details are highly dependent upon which strategy and sequencing platform is used.

In order to convert RNA to DNA the RNA must be used as a template for DNA polymerase. Most DNA polymerases cannot use RNA as a template. However, retroviruses encode a unique type of polymerase known as reverse transcriptases, which are able to synthesize DNA using an RNA template.

All of the current protocols utilize the ability of reverse transcriptase (RT) to synthesize a DNA strand using RNA as a template. RT, like other polymerases, requires a primer annealed to either DNA or RNA to initiate polymerization. Several first-strand priming options are used. These include the following:

1. Using oligo-dT to prime off of the poly-A tail of mature mRNA.
2. Using random primers to prime at 'random' positions along the RNA molecule.
3. Priming off of oligos that are ligated onto the ends of the RNA.

Each of these priming options has strengths and weaknesses which are discussed more fully below.

More on first-strand synthesis

RTs, like other polymerases, can only initiate synthesis from a primer with a 3′-OH that is annealed to the template strand. Several first-strand priming options are used for RNA-seq. Priming options can be divided into two major categories (1) priming off of an oligo ligated to the 3′-end of the RNA template, or (2) using primers with random sequences that can direct synthesis at random locations throughout the RNA template. The most commonly used priming methods are discussed more fully below.

Figure 3.1

3.2.1 First-strand priming using oligo-dT

Copy this link to clipboard

One of the oldest methods for first-strand priming involves using oligo-dT to prime synthesis off of the poly-A tail that is found on most mature eukaryotic mRNAs. This method has the advantage that since the priming sequence is the same for all mRNAs then they should all be equally primed regardless of their coding sequence. However, this approach has several distinct disadvantages. The first is that non-polyadenylated RNAs are lost. If one is only interested in characterizing the expression of coding RNAs then this isn't really an issue. The other problem with this approach is related to the fact that RTs are not highly processive polymerases and can prematurely terminate. This phenomenon can lead to biased enrichment of 3′-ends relative to 5′-ends, which can interfere with downstream quantification and analysis.

Considerations before using oligo-dT for 1st strand priming

1. Oligo-dT priming is not compatible with RNA that has been fragmented since only the 3′-terminal fragments would be recovered. If using oligo-dT for 1st strand priming then fragmentation should be performed on the double stranded cDNA.
2. Oligo-dT priming cannot be used with bacterial mRNAs since they do not possess a poly-A tail.

3.2.2 First-strand priming using random oligos

Copy this link to clipboard

First-strand synthesis can also be primed using primers with random sequences. This is probably the most commonly used technique for first-strand synthesis for RNA-seq. This approach has several advantages. The first is that non-polyadenylated RNAs will be recovered, making it possible to recover ncRNAs and use RNA fragmentation. Secondly, since random primers can anneal throughout the length of the RNA the 3′-bias that is observed with oligo-dT priming is eliminated resulting in a more uniform transcript coverage. Several studies have demonstrated, however, that random priming is not random, and instead priming hot-spots are commonly observed (see for instance: Roberts et al. 2011 and Hansen 2010).

3.2.3 First-strand priming using a pre-ligated oligo

Copy this link to clipboard

Another option for first-strand synthesis is to first ligate an adapter with a known sequence to the 3′-end of the RNA using T4-RNA ligase. This sequence can subsequently be used to prime synthesis of the first strand. This has the advantage of reducing priming bias since all RNAs are primed using the same sequence, and, as discussed below, can be used to retain strand-specific information. This approach is used in both the Small RNA kit from Illumina and in the SOLiD RNA kits from Life Technologies.

The second cDNA strand is synthesized by a DNA polymerase using the RT-synthesized DNA-strand as a template.

Second-strand synthesis also requires a primer that is annealed to the first strand, and again several options exist. These include:

1. Synthesis by RNA nicking and displacement.
2. Using an oligo that is complementary to an adapter pre-ligated to the 5′-end of the RNA template.
3. Using a primer containing a 3′-oligo-dG (this method, referred to as SMART (Zhu et al. 2001) takes advantage of the phenomenon that the MMLV reverse-transcriptase leaves a terminal non-template poly-dC 3′-overhang).

Currently most RNA-seq cDNA libraries are constructed using RNA that has been fragmented as the initial template. However, there are situations where it is preferable to construct cDNA libraries using intact (i.e. unfragmented) RNA.

To sequence the cDNAs, specific adapter sequences must be present at the ends of the fragments. The roles and composition of the adapter sequences vary depending upon the sequencing platform. Adapters contain several different functional elements that are needed for sequencing and may contain one or more optional elements.

Most of the current sequencing platforms use clonal amplification to create clusters of identical molecules that are tethered next to each other on a solid support. For the Illumina platform the clusters are attached to the surface of a flow-cell, while for the 454, IonTorrent, and SOLiD platforms the clusters are generated on beads using emulsion PCR. Regardless of the platform, two types of sequence elements are required: (1) Terminal platform-dependent sequences that are required for clonal amplification and attachment to the sequencing support. (2) Sequences for priming the sequencing reaction. In addition several optional elements may be present. These include sequence tags to allow for multiplexing (known as barcodes or indices), and a second sequence-priming site to allow for sequencing of the insert from the other side (known as paired-end sequencing). Figure 3.2 and Table 3.1 provide additional information regarding these adapter elements.

Figure 3.2

As discussed above, platform specific sequences must be present at the end of the fragments before they can be sequenced (as shown in Figure 3.2). Several strategies have been developed for adding adapter sequences to create the asymmetrical molecule.

More on addition of adapters

3.6.1 Addition of adapters via Y-adapter PCR

Copy this link to clipboard

This method is one of the most commonly used and is essentially the same method that was developed for high-throughput sequencing of genomic DNA. Once the RNA is converted to double stranded cDNA the ends are blunted and adenosine overhangs are added. The adapter sequences can be added using a technique commonly referred to as Y-adapter PCR (see Figure 3.3).

Figure 3.3

3.6.3 Addition of adapters via ligation

Copy this link to clipboard

Although the technical details differ, this approach is used in the Illumina TruSeq Small RNA kit, the NEBNext Small RNA prep kit, and in the SOLiD RNA kits from Life Technologies. These kits use ligation procedures that allow two different adapters to be ligated onto each end of the target RNA. These adapters are then used to prime the first- and second-strand synthesis reactions resulting in cDNAs terminated by the appropriate adapter sequences.

Figure 3.4 shows an outline of the procedure used in both the Illumina TruSeq and in the NEBNext small RNA prep kits. The first step involves ligating an oligo to the 3′-OH of the RNA and takes advantage of the enzymatic activity of T4 RNA ligase II-truncated (T4 Rnl2-trunc) which is a modified form of T4 RNA ligase II. This enzyme can specifically ligate an oligo possessing a 5′-pre-adenylation group to RNA with a 3′-OH. Unlike the full-length ligase, T4 Rnl2-trunc is unable to adenylate the 5′-end of the RNA, preventing circularization of the RNA molecule.

The 5′-adapter is then admixed with T4 RNA ligase and ATP, all of which promotes ligation of the 5′-adapter onto the 5′-end of the RNA. One of the biggest challenges with this procedure is the formation of 5′-adapter-3′-adapter products that can form during the second reaction. In the NEBNext Small RNA prep kit this reaction is quenched by the addition of the RT-primer before the second ligation step. The RT-primer is a DNA oligo that contains sequence complementary to the 5′-end of the 3′-adapter oligo. Free 3′-adapters therefore form a double-stranded molecule that is unable to act as a substrate for the ligation reaction.

Figure 3.4

Conventional RNA-seq library protocols involve converting the RNA to double-stranded cDNA followed by the ligation of sequencing adapters. Although Y-adapter PCR results in an asymmetrical molecule, in reality the final libraries contain two populations of molecules with respect to the original RNA template. In one population the strand that is sequenced corresponds to the sense strand while in the other the strand that is sequenced corresponds to the anti-sense strand. Thus the 'strandedness' of the original RNA template is lost.

Why does this matter? Knowledge of which strand a read comes from may not be critical for measuring differential gene expression since most read aligners (see "Aligning RNA-seq reads to a reference") have settings that allow them to evaluate both the initial read and its reverse-complement. However, for genomes in which transcribed regions overlap (not most eukaryotic organisms), knowledge of strandedness may help assign reads to genes adjacent to one another but on opposite strands. In addition, if one is primarily interested in annotation strandedness can be used to discover features such as the existence of overlapping anti-sense transcripts.

More on stranded protocols

Various 'strand-specific' library preparation protocols have been developed and several of these were reviewed in Levin et.al. 2009. The method that has become the most popular (and performed well according to Levin et.al. 2009) takes advantage of the fact that E. coli DNA polymerase can synthesize DNA using either dUTP or dTTP. In this protocol dUTP is used for synthesis of the second strand, which essentially marks this strand. Y-adapters are then ligated onto both ends in the same manner as in the non-strand-specific protocol. Two approaches have been developed to selectively amplify the non-dUTP strand. One approach (see Figure 3.5 below and outlined in Borodina et al. 2011) involves adding uracil-DNA-glycosylase (UDG) which cleaves the second strand where dUTP has been incorporated, effectively preventing this strand from being amplified during the subsequent PCR amplification. The second approach (which is currently used in the Illumina TruSeq Stranded RNA-seq library preparation kit) involves using a DNA polymerase that is incapable of extending past dUTP bases.

Figure 3.5

To achieve the highest quality of data, it is important to validate the quality and accurately quantify the cDNA libraries before sequencing. As discussed in the RNA preparation section, each of the commonly used methods has certain drawbacks. Absorbance-based methods can lead to erroneous measurements due to residual primers and dNTPs and may not be sensitive enough to accurately measure the library concentrations if the library yield is low. Fluorescence-based techniques generally are more sensitive to RNA concentration but less sensitive to contaminants. If multiplexing multiple libraries in a single sequencing reaction, make sure to adjust the relative contributions from each library to suit your experimental design needs. Other robust options for quantification include quantitative PCR (qPCR) and estimation of fragment concentration using electrophoretic instruments like the ABI Bioanalyzer.

Assessing the fragment size distribution of final RNA-seq libraries is also important. The sizes of the molecules should fall into the expected size range (depending on protocol), and the libraries should be free of ligation/PCR artifacts (usually small relative to the desired molecules). Fragment sizes can be evaluated via electrophoresis, preferably using a sensitive instrument such as an ABI Bioanalyzer.

Sequencing technologies are rapidly expanding and new sequencing chemistries are being developed at a rapid pace. A number of sequencing platforms are widely available, including Illumina, Ion Torrent, and PacBio. The first nanopore sequencing devices are now being marketed by Oxford Nanopore while several legacy technologies can still be found, including Life’s 454 and SOLiD. Each is based upon different proprietary chemistries and technologies and each has unique strengths and weaknesses. A good overview of several platforms is provided by Metzker (2009).

One of the first steps that must be done when designing an RNA-seq experiment is choosing a sequencing platform. It is critical to understand that short-read sequencing requires the addition of specific adapter sequences to the ends of the cDNA molecules of interest. (This is covered in more detail under "Library Preparation"). Since each platform uses unique proprietary adapters the correct adapters must be added to match the sequencing platform.

Although many sequencing platforms exist some are better suited for RNA-seq than others. The current leading platform for RNA-seq is Illumina. This platform enables deep sequencing which is generally important for RNA-seq, and provides long enough, low-error reads that are suitable for mapping to reference genomes and transcriptome assembly. For detailed comparisons of the available platforms see Liu et al. (2012) and Glenn (2011). A useful online comparison can be found at bluseq. The PacBio platform is also beginning to be used for transcriptome assembly, although its use is still in the early stages. Each read from a PacBio is often long enough to recapitulate a full length cDNA transcript.

Before initiating an RNA-seq experiment users should first consult with the sequencing facility they intend to use to make sure they understand the sample preparation and submission requirements. Some general issues that need to be considered are:

1. That the samples are clean and free of major contaminants.
2. The primary DNA molecules contain inserts of the correct size.
3. The primary DNA molecules have adapters on each end.
4. The sample concentration is appropriate.
5. The samples are suspended in appropriate buffers.

Furthermore, the facility needs to know:

1. The sequence of the sequence-priming site.
2. The length of the read you desire.
3. Whether you want single-end or paired-end reads.
4. Whether there is a barcode or index sequence.
5. If using Illumina sequencing the facility also needs to be notified if the inserts contain a region of low sequence complexity immediately after the sequence-priming site (i.e. a barcode).

The analysis and interpretation of genomic data generated by sequencing are arguably among the most complex problems modern scientists face. A staggering degree of effort by biologists, computer scientists, mathematicians, and statisticians is currently leveled at curating, manipulating, and interpreting this information. RNA-seq is subject to many of the same issues other sequencing applications face, in addition to some nuances. Please keep in mind that the nature of questions people may address using RNA-seq data is effectively limitless, so there are even more facets to the analysis of transcriptomes than there are to generating the data itself. As a disclaimer, we emphasize that it is virtually impossible to stay abreast of all developments in approaches to analyzing RNA-seq data. Our objective rather is to provide in the following outline a rough guide to commonly encountered steps and questions one faces on the path from raw RNA-seq data to biological conclusion.

To generalize, however, embarking on an RNA-seq study usually requires initial filtering of sequencing reads, assembling those reads into transcripts or aligning them to reference sequences, annotating the putative transcripts, quantifying the number of reads per transcript, and statistical comparison of transcript abundance across samples or treatments.

Before either de novo assembly or alignment of RNA-seq reads to a reference, they must be processed for a number of reasons described in this section.

Outline:

1. Demultiplex by index or barcode
2. Remove adapter sequences
3. Trim reads by quality
4. Discard reads by quality/ambiguity
5. Filter reads by k-mer coverage
6. Normalize k-mer coverage

If you can't align RNA-seq data to a well-assembled reference genome from a relatively recently diverged organism, building a reference transcriptome out of RNA-seq reads is an alternative strategy. The primary goal in assembling a transcriptome de novo is to reconstruct a set of contiguous sequences (contigs) presumed to reflect accurately a large portion of the RNAs actually transcribed in the cells of your organism. This turns out to be no trivial task, because a limited amount of information about the original gene transcripts is retained in the short reads produced by a sequencer. The assembly process must place reads from the same transcript together in the face of variants introduced by polymorphism and sequencing errors, and the process must assemble reads from different but often similar, paralogous transcripts as separate contigs. Furthermore, repetitive regions that exist within genomes make assembly especially difficult, although this is probably less of a problem for transcriptome assembly than genome assembly.

Other important considerations to make when generating data for a reference transcriptome assembly are sequencing platform and read length. While there is arguably no single best type of data for a de novo transcriptome assembly, most have found that a combination of longer reads and high sequencing depth yields the best results. For example, generating 100 million or more, longer (e.g. 100 nt) Illumina paired-end reads will strike this balance. Alternatively, a mixture of shorter, high-throughput SOLiD and longer, lower-throughput 454 reads may converge on the same optimum. Read length and coverage do matter when it comes time to assemble the transcriptome, however, as different de novo assemblers employ different approaches that depend heavily on these parameters. A newer method involves using the PacBio platform to generate long sequencing reads (typicaly 1-10kb in length). The PacBio sequencer is generally able to generate a full, novel transcript with each read produced on the platform.

Here we do not attempt to present a detailed, comprehensive review of how de novo assembly works, nor do we dissect the nuances of different assembly approaches. Other published reviews have tackled these topics. Solutions to sequence assembly arose from the field of mathematics known as graph theory. For general reviews of how genomic sequencing data is assembled using graphs, see (Miller 2010; Compeau 2011; Li 2011). For a specific review of transcriptome assembly, see (Martin and Wang 2011), and for comparisons of how choice of assembly software and sequencing depth affect de novo transcriptome assemblies, see (Brautigam 2011; Francis 2013). Below we summarize at a superficial level the general approaches currently used for de novo transcriptome assembly, along with examples of commonly used software. These approaches were designed with genome assembly in mind but have been adapted for transcriptome assembly as necessary. Finally, we provide a basic comparison of actual de novo transcriptome assemblies for a non-model teleost (pipefish), specifically addressing assembler differences as a function of contig length.

Short reads generated by RNA-seq experiments must ultimately be aligned, or "mapped" to a reference genome or transcriptome assembly. Read alignment to a reference provides biological information in two basic ways. First, it generates a dictionary of the genomic features represented in each RNA-seq sample. That is, aligned reads become annotated, and the highly fragmented data is thus connected to the gene families, individual transcripts, small RNAs, or individual exons encompassed by the original tissue sample. Second, the number of reads aligned to each feature approximates abundances of those features in the original sample. Such measures of digital gene expression are then subject to comparison among samples or treatments in a statistical framework.

The general objective when mapping or aligning a collection of sequencing reads to a reference is to discover the true location (origin) of each read with respect to that reference. Obviously features of the reference such as repetitive regions, assembly errors, and missing information will render this objective impossible for a subset of the reads. Likewise, sequencing errors, polymorphisms, and limited complexity within the short reads will make them harder to align. In light of these factors, aligners must be flexible when applying mapping criteria. That is, they must allow for approximate matches. Otherwise, a large proportion of sequencing reads will not be assigned to a particular feature in the reference. Just how "approximate" the matches are depend on the features of the aligner, usually manifested as user-specified parameters such as the proportion of bases in each read required to match the reference identically, and penalties incurred by insertions, deletions, and gap sizes. There is no easy answer to which aligner one should choose or which alignment parameter values one should select for a particular aligner. Like many other steps in RNA-seq analysis, there really is no substitute for a little exploration of parameter space. Obviously alignment parameter settings that are extremely stringent will result in a small subset of your reads being mapped. Very liberal settings, on the other hand, will result in lost specificity, and many reads will map to multiple features of the reference. It therefore takes a little bit of experimentation to get the optimal balance of sensitivity and specificity for a given data set.

There is no short supply of short read aligners. In fact, dozens have been developed within the last several years, and some (but not all) were designed with specific utilities in mind. Several effective and popular aligners include Bowtie2, BWA-mem, and GSnap. Fonseca et al. 2012 maintain an updated, near-comprehensive database of aligners here. In the following section we outline several broad categories into which short read aligners fall, specifically with respect to RNA-seq. Keep in mind that using alignment software generally requires two basic steps. First, the reference genome or transcriptome (generally in ".fasta" format) is converted to an indexed "reference" by the software in order to allow fast read mapping, which is the second step. Both building the reference and aligning sequencing reads are usually quite straightforward given specific instructions in software manuals/websites. Depending on the speed of the aligner (some are much slower than others) and the number of reads you have to align, consider breaking up reads into multiple files and running alignments in parallel, across multiple nodes on a computer cluster, for example.

Most people are interested in the identities of the transcripts they assemble or identify as differentially expressed or alternatively spliced. Annotations connect novel RNA-seq data with biological information about the same (or homolgous) transcripts derived from other studies. For example, it is important to identify gene families or known pathways that respond transcriptionally to specific experimental treatments. Also, RNA-seq data derived from eukaryotic organisms may contain reads from other species (parasites, for example), and one would want to know about this.

In the event that a reference genome is available, annotation is relatively straightforward. Genomic coordinates from the reference genome are normally associated with various forms of annotation information through databases such as genome browsers. Given mapping coordinates from an alignment, for example, it is possible to extract the corresponding information from such databases.

A transcriptome assembled de novo, on the other hand, is often annotated from scratch. The most common way to annotate transcripts in this situation is to use the NCBI-supported BLAST (Basic Local Alignment Search Tool) (Altschul et al. 1990), the recent, widely used suite of tools being "BLAST+" (Camacho et al. 2009). A description of the BLAST+ command line applications can be found at www.ncbi.nlm.nih.gov/books/NBK1763/. In essence, the BLAST model allows the user to "match" query sequences to one or more databases of curated, annotated sequences, using an efficient local sequence alignment approach. In some situations it may be adequate to blast contigs against a database of known or predicted transcripts from the reference genome of a closely-related organism. In other situations it may be desirable to blast contigs against all nucleotide sequences in an inclusive database, such as the NCBI nucleotide database ("nt"), which itself is composed of entries from smaller databases. If the annotation emphasis is on protein-coding transcripts, BLASTx, which translates each query sequence (in all six reading frames) to amino acid sequences and uses these to query a protein database, is an appropriate tool. All NCBI databases are available at ftp.ncbi.nlm.nih.gov/blast/db/.

Once an automated BLAST search has been performed for a collection of sequences (contigs from a reference transcriptome, for example), one of several specifiable output formats may be parsed to extract information about the number of "hits", the probability that the alignment of a query to each hit occurred by chance (the "e-value"), the descriptions and accession numbers for each hit, and representations of the sequence alignments themselves. This information may then be used in a variety of ways, including basic annotation, evaluation of redundancy in the assembly, or functional enrichment analyses.

One generally useful way to compare groups of transcripts or contigs that are differentially expressed or otherwise is through Gene Ontology (GO) term analysis (www.geneontology.org/GO.doc.shtml). Entries in sequence databases such as those from NCBI are associated with a number of biologically relevant terms. The terms belong to one of three basic "ontologies:" cellular component, biological process, and molecular function. These terms are either associated with the database sequences directly via molecular genetic experimental evidence, or indirectly via homology with other sequence entries in the database. It is possible to take the database entry identifiers for a set of blast hits, extract GO term information, and quantify/compare GO term frequencies among different groups of original sequences. Several software platforms accomplish this type of workflow, including:

• Blast2GO (Conesa et al. 2005), available at: blast2go.de/b2ghome
• DAVID (Huang da, Sherman, Lempicki 2009), available at: david.abcc.ncifcrf.gov/
• ermineJ (Gillis, Mistry, Pavlidis 2010), available at: www.chibi.ubc.ca/ermineJ/

The ultimate objective of many RNA-seq studies is to identify features of the transcriptome that are differentially expressed between or among groups of individuals or tissues that are different with respect to some condition. It is important to keep in mind the nature of those "features" to which reads were mapped for counting. If, for example, reads are mapped to a set of references that include different splice variants of the same gene, those splice variants will each be analyzed separately unless treated as a unit. The nature of the features encompassed by read count data depends on what the mapping reference is (e.g. a genome or transcriptome assembly), and what type of aligner was used. Some pieces of analysis software are designed to handle isoform differences (Leng et al. 2013; Trapnell et al. 2013), and others analyze generic features of the transcriptome, so please consider this when selecting both read alignment and differential expression software.

A basic example of differential expression analysis would be the comparison of liver transcriptomes between a group of mice treated with a pharmacological compound and an untreated "control" group. For more details about experimental design in RNA-seq studies, please see the section on Experimental Design. In the past, microarrays were the most commonly used tool for transcriptome-wide interrogation of differential expression. In many contexts microarrays are still the most appropriate choice today, but RNA-seq is a generally more flexible approach. For example, as we have discussed, no previous genomic sequence information from a given organism is necessary for an RNA-seq study, whereas an existing library of expressed sequence tags is required for the construction of expression arrays. Like signals from RNAs hybridized to microarray probes, RNA-sequencing read counts correspond well with actual transcript levels in a sample (Mortazavi et al. 2008). In terms of differential expression analysis, RNA-seq data (counts of mapped reads) are fundamentally different from microarray data (computed scores of image intensity). Unnormalized RNA-seq data, therefore, must be modeled statistically using discrete distributions. In some differential expression analysis methods, however, RNA-seq data are first normalized to account for a number of library- and/or gene-specific biases (explained below), treated as a continuous variable of transcript abundance, and therefore modeled using continuous distributions for statistical inference.

The goal of this section is to provide a brief introduction to the general statistical approaches used to test hypotheses of differential expression, touch on normalization of RNA-seq data and how the statistical approaches implement normalization in a general sense, and provide an annotated list of currently available software for differential gene expression analysis. In lieu of repeating what other authors have accomplished, we refer our audience to several excellent reviews, which go into greater detail regarding the statistical models and test drive many analysis programs with real and simulated data (Bullard et al. 2010; Kvam, Liu, Si 2012; Soneson, Delorenzi 2013).

Anders, S, W Huber. (2010). Differential expression analysis for sequence count data. Genome Biol 11:R106

Benjamini, Y, Y Hochberg. (1995). Controlling the False Discovery Rate - a Practical and Powerful Approach to Multiple Testing. Journal of the Royal Statistical Society Series B-Methodological 57:289-300

Bullard, JH, E Purdom, KD Hansen, S Dudoit. (2010). Evaluation of statistical methods for normalization and differential expression in mRNA-Seq experiments. BMC Bioinformatics 11:94

Busby, MA, C Stewart, C Miller, K Grzeda, G Marth. (2013). Scotty: A Web Tool For Designing RNA-Seq Experiments to Measure Differential Gene Expression. Bioinformatics

Levin JZ, M. Yassour, X. Adiconis, et al. (2010) Comprehensive comparative analysis of strand-specific RNA sequencing methods. Nature methods. 7, 709-715

McIntyre, LM, KK Lopiano, AM Morse, V Amin, AL Oberg, LJ Young, SV Nuzhdin. (2011). RNA-seq: technical variability and sampling. BMC Genomics 12:293

Noble, WS. (2009). How does multiple testing correction work?. Nat Biotech 27:1135-1137

Sokal, R.R., and F.J. Rohlf. (1995). Biometry: The principles and practice of statistics in biological research. 3rd edition. W.H. Freeman, New York.

Storey, JD. (2002). A direct approach to false discovery rates. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 64:479-498

Wang, Y, N Ghaffari, CD Johnson, UM Braga-Neto, H Wang, R Chen, H Zhou. (2011). Evaluation of the coverage and depth of transcriptome by RNA-Seq in chickens. BMC Bioinformatics 12 Suppl 10:S5

Anders, S, W Huber. (2010). Differential expression analysis for sequence count data. Genome Biol 11:R106

Aranda R, Dineen SM, Craig RL, Guerrieri RA, Robertson JM. (2009). Comparison and evaluation of RNA quantification methods using viral, prokaryotic, and eukaryotic RNA over a 10⁴concentration range. Anal Biochem 387: 122-127.

Chomczynski P, Sacchi N. (1987). Single-step method of RNA isolation by acid guanidinium thiocyanate-phenol-chloroform extraction. Anal Biochem 162: 156-159.

Chomczynski P, Sacchi N. (2006) The single-step method of RNA isolation by acid guanidinium thiocyanate-phenol-chloroform extraction: twenty-something years on. Nat Protoc 1: 581-585.

Ekblom R., Slate J., Horsburgh, et al. (2012) Comparison between Normalised and Unnormalised 454-Sequencing Libraries for Small-Scale RNA-Seq Studies, Comparative and functional genomics. 2012, 281693.

Forconi M, Herschlag D. (2009) Metal ion-based RNA cleavage as a structural probe. Meth Enzymol 468: 91-106.

Kirby KS. (1956) A new method for the isolation of ribonucleic acids from mammalian tissues. Biochemical Journal 64: 405.

Matsumura H., Molina C., Kruger D.H., Terauchi R., and Kahl G. (2012) DeepSuperSAGE: High-Throughput Transcriptome Sequencing with Now- and Next-Generation Sequencing Technologies. Tag-based Next Generation Sequencing , First Edition. eds. M. Harbers & G. Kahl. Wiley-VCH Verlag GmbH & Co. KGaA.

Raz T, Kapranov P, Lipson D, Letovsky S, Milos PM, Thompson JF. (2011) Protocol dependence of sequencing-based gene expression measurements. PLoS ONE 6.

Volkin E, Carter CE. (1951) The Preparation and Properties of Mammalian Ribonucleic Acids. J Am Chem Soc 73. 1516-1519.

Zhulidov PA, Bogdanova EA, Shcheglov AS, Vagner LL, Khaspekov GL, Kozhemyako VB, Matz MV, Meleshkevitch E, Moroz LL, Lukyanov SA, Shagin DA. (2004) Simple cDNA normalization using kamchatka crab duplex-specific nuclease. Nucleic Acids Res. 32 (3):e37.

Zhulidov PA, Bogdanova EA, Shcheglov AS, Shagina IA, Wagner LL, Khazpekov GL, Kozhemyako | VV, Lukyanov SA, Shagin DA. (2005) A method for the preparation of normalized cDNA libraries enriched with full-length sequences. Bioorg Khim. 31 (2):186-94.

Borodina T, Adjaye J., and Sultan M. (2011) A strand-specific library preparation protocol for RNA sequencing, Methods in enzymology. 500, 79-98.

Hansen et al. (2010) Biases in Illumina transcriptome sequencing caused by random hexamer priming. Nucleic acids research. 38, e131.

Levesque-Sergerie J-P, Duquette M, Thibault C, Delbecchi L, Bissonnette N. (2007). Detection limits of several commercial reverse transcriptase enzymes: impact on the low- and high-abundance transcript levels assessed by quantitative RT-PCR. BMC Mol Biol 8: 93.

Levin J.Z., Yassour M. , Adiconis X., et al. (2010) Comprehensive comparative analysis of strand-specific RNA sequencing methods, Nature methods. 7, 709-715.

Okello et al. (2010). Quantitative assessment of the sensitivity of various commercial reverse transcriptases based on armored HIV RNA. PLoS ONE 5: e13931.

Roberts et al. (2011) Improving RNA-Seq expression estimates by correcting for fragment bias. Genome biology. 12,R22.

Ståhlberg A, Kubista M, Pfaffl M. (2004). Comparison of reverse transcriptases in gene expression analysis. Clin Chem 50:1678-1680.

Zhu YY, Machleder EM, Chenchik A, Li R, Siebert PD. (2001). Reverse transcriptase template switching: a SMART approach for full-length cDNA library construction. BioTechniques 30: 892-897.

Glenn T.C. (2011) Field guide to next-generation DNA sequencers, Molecular ecology resources. 11, 759-769.

Liu L., Li Y., Li S., et al. (2012) Comparison of next-generation sequencing systems, Journal of biomedicine & biotechnology. 2012, 251364.

M.L. Metzker M.L. (2009) Sequencing technologies - the next generation, Nature reviews Genetics. 11, 31-46.

Altschul, SF, W Gish, W Miller, EW Myers, DJ Lipman. (1990) Basic local alignment search tool. J Mol Biol 215:403-410

Anders, S, W Huber. (2010). Differential expression analysis for sequence count data. Genome Biol 11:R106.

Auer, PL, RW Doerge. (2011) A two-stage poisson model for testing RNA-Seq data.

Blanca, JM, L Pascual, P Ziarsolo, F Nuez, J Canizares. (2011). ngs_backbone: a pipeline for read cleaning, mapping and SNP calling using Next Generation Sequence. BMC Genomics 12.

Brautigam, A, T Mullick, S Schliesky, AP Weber. (2011) Critical assessment of assembly strategies for non-model species mRNA-Seq data and application of next-generation sequencing to the comparison of C(3) and C(4) species. J Exp Bot 62:3093-3102.

Brown, CT, A Howe, Q Zhang, AB Pyrkosz, TH Brom (2012) A Reference-Free Algorithm for Computational Normalization of Shotgun Sequencing Data arXiv:1203.4802v2 [q-bio.GN]

Bullard, JH, E Purdom, KD Hansen, S Dudoit. (2010) Evaluation of statistical methods for normalization and differential expression in mRNA-Seq experiments. BMC Bioinformatics 11:94

Busby, MA, C Stewart, C Miller, K Grzeda, G Marth (2013). Scotty: A Web Tool For Designing RNA-Seq Experiments to Measure Differential Gene Expression. Bioinformatics

Camacho, C, G Coulouris, V Avagyan, N Ma, J Papadopoulos, K Bealer, TL Madden. (2009). BLAST+: architecture and applications. BMC Bioinformatics 10:421

Catchen, JM, A Amores, P Hohenlohe, W Cresko, JH Postlethwait (2011). Stacks: Building and Genotyping Loci De Novo From Short-Read Sequences. G3: Genes, Genomes, Genetics 1:171-182

Chevreux, B, T Pfisterer, B Drescher, AJ Driesel, WE Muller, T Wetter, S Suhai (2004). Using the miraEST assembler for reliable and automated mRNA transcript assembly and SNP detection in sequenced ESTs. Genome Res 14:1147-1159

Chung, LM, JP Ferguson, W Zheng, F Qian, V Bruno, RR Montgomery, H Zhao (2013). Differential expression analysis for paired RNA-seq data. BMC Bioinformatics 14:110

Cock, PJ, CJ Fields, N Goto, ML Heuer, PM Rice (2010). The Sanger FASTQ file format for sequences with quality scores, and the Solexa/Illumina FASTQ variants. Nucleic Acids Res 38:1767-1771

Compeau, PE, PA Pevzner, G Tesler (2011). How to apply de Bruijn graphs to genome assembly. Nat Biotechnol 29:987-991

Conesa, A, S Gotz, JM Garcia-Gomez, J Terol, M Talon, M Robles (2005). Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics 21:3674-3676

Di, Y, DW Schafer, JS Cumbie, JH Chang (2011). The NBP negative binomial model for assessing differential gene expression from RNA-Seq. Stat Appl Genet Mol Biol 10:24

Dillies, MA, A Rau, J Aubert, et al (2012). A comprehensive evaluation of normalization methods for Illumina high-throughput RNA sequencing data analysis. Brief Bioinform

Fonseca, NA, J Rung, A Brazma, JC Marioni (2012). Tools for mapping high-throughput sequencing data. Bioinformatics 28:3169-3177

Francis, WR, LM Christianson, R Kiko, ML Powers, NC Shaner, SH Haddock (2013). A comparison across non-model animals suggests an optimal sequencing depth for de novo transcriptome assembly. BMC Genomics 14:167

Gillis, J, M Mistry, P Pavlidis (2010). Gene function analysis in complex data sets using ErmineJ. Nat. Protocols 5:1148-1159

Grabherr, MG, BJ Haas, M Yassour, et al (2011). Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol 29:644-652

Hardcastle, TJ, KA Kelly (2010). baySeq: empirical Bayesian methods for identifying differential expression in sequence count data. BMC Bioinformatics 11:422

Huang da, W, BT Sherman, RA Lempicki (2009). Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc 4:44-57

Huang, X, A Madan (1999). CAP3: A DNA sequence assembly program. Genome Res 9:868-877

Kvam, VM, P Liu, Y Si (2012). A comparison of statistical methods for detecting differentially expressed genes from RNA-seq data. Am J Bot 99:248-256

Leng, N, JA Dawson, JA Thomson, V Ruotti, AI Rissman, BMG Smits, JD Haag, MN Gould, RM Stewart, C Kendziorski (2013). EBSeq: An empirical Bayes hierarchical model for inference in RNA-seq experiments. Bioinformatics

Li, H, B Handsaker, A Wysoker, T Fennell, J Ruan, N Homer, G Marth, G Abecasis, R Durbin, S Genome Project Data Processing (2009). The Sequence Alignment/Map format and SAMtools. Bioinformatics 25:2078-2079

Li, J, R Tibshirani (2011). Finding consistent patterns: A nonparametric approach for identifying differential expression in RNA-Seq data. Stat Methods Med Res

Li, J, DM Witten, IM Johnstone, R Tibshirani (2012). Normalization, testing, and false discovery rate estimation for RNA-sequencing data. Biostatistics 13:523-538

Li, Z, Y Chen, D Mu, et al (2011). Comparison of the two major classes of assembly algorithms: overlap–layout–consensus and de-bruijn-graph. Brief Funct Genomics

Lohse, M, AM Bolger, A Nagel, AR Fernie, JE Lunn, M Stitt, B Usadel (2012). RobiNA: a user-friendly, integrated software solution for RNA-Seq-based transcriptomics. Nucleic Acids Res 40:W622-W627

Love, MI, Huber, W, S Anders (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology 15:550

Lund, SP, D Nettleton, DJ McCarthy, GK Smyth (2012). Detecting Differential Expression in RNA-sequence Data Using Quasi-likelihood with Shrunken Dispersion Estimates. Stat Appl Genet Mol Biol 11

Martin, JA, Z Wang (2011). Next-generation transcriptome assembly. Nat Rev Genet 12:671-682

Miller, JR, S Koren, G Sutton (2010). Assembly algorithms for next-generation sequencing data. Genomics 95:315-327

Mortazavi, A, BA Williams, K Mccue, L Schaeffer, B Wold (2008). Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nature Methods 5:621-628

Oshlack, A, MD Robinson, MD Young (2010). From RNA-seq reads to differential expression results. Genome Biol 11:220

Pevzner, PA, H Tang, MS Waterman (2001). An Eulerian path approach to DNA fragment assembly. PNAS 98:9748-9753

Robertson, G, J Schein, R Chiu, et al (2010). De novo assembly and analysis of RNA-seq data. Nat Meth 7:909-912

Robinson, MD, DJ McCarthy, GK Smyth (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26:139-140

Robinson, MD, A Oshlack (2010). A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol 11:R25

Schulz, MH, DR Zerbino, M Vingron, E Birney (2012). Oases: Robust de novo RNA-seq assembly across the dynamic range of expression levels. Bioinformatics

Simpson, JT, K Wong, SD Jackman, JE Schein, SJ Jones, I Birol (2009). ABySS: a parallel assembler for short read sequence data. Genome Res 19:1117-1123

Smyth, GK (2004). Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Stat Appl Genet Mol Biol 3:Article3

Soneson, C, M Delorenzi (2013). A comparison of methods for differential expression analysis of RNA-seq data. BMC Bioinformatics 14:91

Tarazona, S, F Garcia-Alcalde, J Dopazo, A Ferrer, A Conesa (2011). Differential expression in RNA-seq: a matter of depth. Genome Res 21:2213-2223

Trapnell, C, DG Hendrickson, M Sauvageau, L Goff, JL Rinn, L Pachter (2013). Differential analysis of gene regulation at transcript resolution with RNA-seq. Nat Biotech 31:46-53

Van De Wiel, MA, GGR Leday, L Pardo, H Rue, AW Van Der Vaart, WN Van Wieringen (2012). Bayesian analysis of RNA sequencing data by estimating multiple shrinkage priors. Biostatistics

Wang, LK, ZX Feng, X Wang, XW Wang, XG Zhang (2010). DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics 26:136-138

Yu, D, W Huber, O Vitek (2013). Shrinkage estimation of dispersion in Negative Binomial models for RNA-seq experiments with small sample size. Bioinformatics

Zerbino, DR, E Birney (2008). Velvet: algorithms for de novo short read assembly using de Bruijn graphs. Genome Res 18:821-829

Zhou, Y-H, K Xia, FA Wright (2011). A powerful and flexible approach to the analysis of RNA sequence count data. Bioinformatics 27:2672-2678.