Mapping the non-standardized biases of ribosome profiling

Harvesting the cells and antibiotic pretreatment

The most delicate step in the sample preparation is the isolation of intact ribosome-mRNA complexes. Ideally, the isolation procedure should faithfully freeze the translating ribosomes and avoid conditions that stimulate ribosomal drop-off and, most importantly, ribosome relocation on the mRNA during the sample processing.

Early in the development of the ribosome profiling approach, cells were pre-incubated with elongation inhibitors (mainly chloramphenicol for bacteria and cychloheximide for eukaryotes) to inhibit further movement of the elongating ribosomes along the mRNA (Ingolia et al., 2009). The antibiotic treatment markedly affects the coverage profiles and introduces some bias in the results; the elongation inhibitors do not uniformly stall elongating ribosomes but rather show a codon-dependent mode of action (Orelle et al., 2013). Cycloheximide also allows one complete translocation cycle before blocking the ribosome (Pestova and Hellen, 2003; Schneider-Poetsch et al., 2010) and thus diffuses the read-out when determining codon-dependent stalling (Nedialkova and Leidel, 2015), while non-antibiotic treated cells deliver much sharper pause sites corresponding to rare codons (Pelechano et al., 2015).

In addition, a broad cumulative peak downstream of the start codon has been seen in the earlier profiling papers that use elongation inhibitors and interpreted as slow initiation (Ingolia et al., 2009). The initial peak, albeit still present, significantly decreases when cell mass is flash-frozen and elongation inhibiters are omitted (Guydosh and Green, 2014; Lareau et al., 2014).

The disproportionately high accumulation of reads at initiation is rather an artifact of the antibiotic pretreatment (Becker et al., 2013) and results from inhibition of translation elongation with ongoing initiation (Ingolia et al., 2011). The antibiotic does not immediately reach the threshold of complete inhibition of elongation; instead its concentration increases gradually in the cell (Gerashchenko and Gladyshev, 2014). Hence, upon treatment, some initiating ribosomes continue into the elongation cycle until they encounter the drug, which results in an excess of ribosomal footprints over the first five to ten codons from the coding sequence (Gerashchenko and Gladyshev, 2014). Additionally, an 80S ribosome stalled in the proximity of the start codon will prevent any subsequent scanning ribosome from reaching the initiation codon, which may result in an apparent stalling at an upstream open-reading frame (uORF). Thus, an initiation site with mediocre context in uORF will be occupied because of the highly efficient but blocked downstream start site (Jackson and Standart, 2015), which may lead to an erroneous interpretation of alternative uORF-induced initiation.

Careful consideration of the effect of antibiotics on ribosome coverage offers little support that the large number of genes with uORFs is involved in shaping the resistance to oxidative stress (Gerashchenko and Gladyshev, 2014). Ribosome profiling without antibiotics prior to cell harvesting revealed that translation of only a small fraction of uORF-bearing mRNA was refractory to oxidative stress (Andreev et al., 2015). Elongation inhibitors added prior to harvesting the ribosome-mRNA complexes alter the distribution of reads in the cumulative ribosome profiles (namely, the aligned and averaged profiles of many genes). For example, emetine-stalled elongating ribosomes give slightly longer fragments than those isolated from cycloheximide-treated mammalian cells, suggesting that various antibiotics stabilize different ribosome conformation (Ingolia et al., 2011). Conversely, drug pretreatment may eliminate some features of biological importance in the cumulative ribosome profiles. For example, antibiotic pretreatment in mammalian cells eliminates the ribosomal peak at the end of the open-reading frames, which is observed in untreated cells (Ingolia et al., 2011).

The most widely applied cell harvesting procedure involves rapid cooling of the cell suspension and centrifugation (Becker et al., 2013) (Figure 1). Bacteria are cooled by pouring the cell suspension over crushed ice, while eukaryotic (mammalian) cells cultured in monolayer are re-suspended in ice-cold PBS supplemented with elongation inhibitor and immediately pelleted by centrifugation (Guo et al., 2010). Tissues are usually flash-frozen and grinded in the lysis buffer supplemented with elongation inhibitor (Gonzalez et al., 2014). An alternative approach for harvesting of cells growing in suspension is a rapid filtration of the cells in a pre-warmed glass nitrocellulose filtration system and flash-freezing the membrane with the cells (Figure 1). So far, this filtration approach has been mainly used in unicellular organisms (yeast and Escherichia coli, for example) (Ingolia et al., 2009; Oh et al., 2011; Li et al., 2012). Both harvesting protocols show good reproducibility between biological replicates (r=0.99, Pearson correlation coefficient) (Becker et al., 2013). Importantly, however, the RPF accumulation at native stalling sites, e.g. SecM and TnaC, is higher using the filtration harvesting (Becker et al., 2013). Most likely, the filtration approach compared to the centrifugation is less susceptible to variations and faithfully halts the translating ribosomes. Still, harvesting by centrifugation might be the only option for cells that cannot be rapidly filtered. However, it is important to perform it as quickly as possible using pre-chilled devices.

In summary, the procedure for isolation of ribosome-mRNA complexes is of crucial importance. While drug pretreatment may not influence differential expression analysis, as the expression of each gene is compared under two different conditions with an otherwise uniform protocol, the use of elongation inhibitors or the harvesting procedure may alter the interpretation of position-specific information.

Cell lysis

Similar to the cell harvesting procedure, the aim at this step is to recover the ribosome-mRNA complexes with minimal losses from ribosomal dissociation (or drop-off) and mRNA degradation. The composition of the lysis buffer is optimized to stabilize the ribosome-mRNA complexes with high concentration of magnesium (between 5 and 20 mm) and an additional salt, such as KCl or NaCl and NH4Cl.

The isolation of intact polysomes is a procedure established in early ribosome research and is still applied today almost unchanged (Wettstein et al., 1963; Dresden and Hoagland, 1965). The composition of the lysis buffer underwent several variations. However, some components of the lysis buffer, if overdosed, may distort the ribosome profiles. For example, high NaCl concentration decreases the monosome peak and enhances the fraction of dissociated ribosomal subunits (Becker et al., 2013); high salt concentration increases the fraction of vacant ribosomes that are not engaged in translation (Blobel and Sabatini, 1971) and consequently decreases the number of RPF. Magnesium stabilizes the translating ribosomes (Ron et al., 1968) and at high concentrations freezes the conformational changes in the bacterial ribosome (Blanchard et al., 2004). Moreover, high magnesium concentration induces folding of the mRNA, which hinders the subsequent nucleolytic digestion (Andreev et al., 2015). Lowering the magnesium concentration from 15 mm to 5 mm greatly improves the codon positioning of the footprints and the resolution of the ribosome profiling (Ingolia et al., 2012). Also, low magnesium conditions permit conformational flexibility of the ribosome and create heterogeneity in the length ribosomal footprints (Lareau et al., 2014); the variant ribosomal footprints are informative on distinct stages of the translating ribosome during the elongation cycle.

The lysis buffer also contains an elongation inhibitor to additionally stabilize the ribosome-mRNA complexes during sample processing. The binding kinetics of the antibiotic when present in the cell lysis is rapid compared to the diffusion-driven process of antibiotic enrichment in intact cells during the pretreatment procedure. Generation of cell extracts from Saccharomyces cells in the cycloheximide-containing lysis buffer faithfully halted the ribosomes along the mRNA with no distortion (Guydosh and Green, 2014). Cycloheximide should be preferred over alternative substances that stabilize eukaryotic ribosome-mRNA complexes, e.g. the non-hydrolyzable GTP analog GMP-PNP, as they slightly increase the size of the ribosome footprints (Guydosh and Green, 2014). Although such studies with bacterial elongation inhibitors are missing, it can be expected that their mode of action will be similar to that of the cycloheximide when added to the lysis buffer.

Along with variations in the composition of the lysis buffer, the lysis procedure also varies. In general, despite the presence of components stabilizing the ribosome-mRNA complexes (e.g. elongation inhibitors, magnesium) to avoid ribosomal reallocation or dissociation, lysis is usually carried out at low temperatures by either adding frozen drops of lysis buffer to a frozen cell powder or flash-freezing with the cell mass. When this is not applicable, i.e. by ribosome profiling of tissues, the lysis buffer is generally added to the sample ice-cold (Gonzalez et al., 2014).

Eukaryotic cells are lysed on ice by repeated micropipetting or homogenization (Guo et al., 2010; Becker et al., 2013; Chew et al., 2013). Pulverized bacteria or monocellular eukaryotes are homogenized in a mill with liquid nitrogen (Oh et al., 2011; Guydosh and Green, 2014; Woolstenhulme et al., 2015). This method is transferrable to any cell type and frozen tissue and should be the preferred lysis approach as it allows treatment of the sample at very low temperatures. During the homogenization, local temperature fluctuations in the sample should be avoided by careful choice of the conditions, i.e. short homogenization pulses and pre-cooling the grinder jar before and after each homogenization cycle (Oh et al., 2011; Guydosh and Green, 2014; Woolstenhulme et al., 2015).

Nucleolytic generation of ribosomal footprints

The clarified lysate is then digested with a nuclease to generate monosomes (Figure 1). RNase I has been exclusively used in eukaryotic ribosome profiling (Ingolia et al., 2012) and micrococcal nuclease (MNase) from Staphylococcus aureus in bacteria; RNase I is inactive in bacteria (Datta and Burma, 1972). MNase can also be used in eukaryotic lysates (Reid and Nicchitta, 2012; Dunn et al., 2013), and, in fact, it leads to a reduced amount of ribosomal RNA (rRNA) contamination compared to RNase I treatment (Oh et al., 2011; Miettinen and Bjorklund, 2015). The activity of the MNase is modulated by calcium ions. A disadvantage of MNase is its preferential cleavage at A or T nucleotides (Dingwall et al., 1981) and consequently, the MNase-generated ribosome footprints might be enriched in A or T nucleotides at their 5′ ends. Compared to fragments derived from yeast lysates treated with RNase I, the MNase-generated footprints are more heterogeneous in length (Becker et al., 2013) due to steric effects and less precise 5′ cleavage (Woolstenhulme et al., 2015). In contrast, MNase cleaves precisely at the 3′ end contour of the ribosome, thus the calibration of the reads in bacterial system should be preferably done using the 3′ ends of the reads (Woolstenhulme et al., 2015) (see section ‘Analysis of the sequencing data’). RNase I cleavages are precise at both 5′ and 3′ ends, enabling calibration using both termini. Conversely, RNase I-treated samples show a slight bias towards enrichment of short genes (Miettinen and Bjorklund, 2015), although the reason for this remains unclear.

Contamination with rRNA fragments released by the nucleolytic digestion substantially decreases the amount of informative sequencing data. Importantly, the rRNA fragments generated during the nucleolysis of the polysomes are species-specific, but are limited to only few fragments and can be efficiently removed to near completeness by using few complementary oligonucleotides. Thus, in setting up a protocol for ribosome profiling in a new cell line or species, it is recommendable by to perform a pioneer sequencing run to identify the contaminant rRNA species and design specific oligonuclotides for the depletion of rRNA-derived fragements.

Finally, the amount of each nuclease needs a careful determination; enhanced nuclease activity (caused either by large amounts of enzyme, pH variations or long digestion times) leads primarily to an increased contamination of the ribosome footprint libraries with rRNA fragments. By contrast, insufficient amount of MNase causes less stringent cleavage of the mRNA and results in longer fragments which migrate outside of the range selected for ribosomal fragments during the gel purification procedure. Consequently, it will yield lower depth and coverage of the mRNAs and it will decrease the accuracy in determining ribosome positions along mRNAs (Becker et al., 2013).

Generation of the deep-sequencing library

The preparation of libraries for deep sequencing involves fusion of adapters to the generated small DNA or RNA fragments. This process also contains biases and a detailed knowledge is of crucial importance to avoid erroneous interpretation of the data. A recent review summarizes the critical caveats in each step of library preparation (van Dijk et al., 2014). Here, we only compare various methods for adaptor ligations to the ribosomal footprints, which are unique to the ribosome profiling procedure. In principle, after nucleolytic digestion the ribosome profiling follows the typical steps of library preparation in the micro RNA-Seq methodology (Guo et al., 2010), including sequential adaptor ligation, reverse transcription of the RNA fragments and PCR amplification of the transcribed DNA. The earliest approach uses circularization of the fragments to fuse adaptors at both ends (Ingolia et al., 2009). Prior to this, each fragment is polyadenylated at its 3′ ends with poly(A)-polymerase (Ingolia et al., 2009), which serves as a priming site for the reverse transcription. Polyadenylation was also introduced to produce uniform 3′ ends of all fragments and to reduce the bias in the ligation (Ingolia, 2010), however the sequenced fragments are enriched in adenines at their 3′ termini (Artieri and Fraser, 2014). Furthermore, in the circularization procedure, an additional preference for adenine at the first 5′-position is observed (Lamm et al., 2011; Artieri and Fraser, 2014): it does not depend on the polyA-tails of the fragments and the origin of this bias is unknown.

Later developments in the library preparation of ribosomal footprints use ligation approaches established in the sequencing of miRNAs, in which 3′ and 5′ adaptors are ligated sequentially to the fragments without circularization (Guo et al., 2010). This allowed capture of low-abundance fragments and omitted the sequence bias (i.e. the preference for adenines at 5′ and 3′ positions). Note that direct ligation of a 3′ adapter might be applied also as an alternative to polyA-tailing, preceding the circularization approach. However, some sequences in the libraries generated with sequential adaptor ligation were overrepresented compared to a sequencing in which the adaptors were ligated using the circularization protocol. The overrepresented fragments are a consequence of local secondary structure preferences (Hafner et al., 2011; Zhuang et al., 2012) and their propensity to co-fold with the adaptor sequences (Jackson et al., 2014). Using truncated T4 RNA Ligase 2 instead of the previously used full-length, non-truncated version decreased the amount of those fragments by a half (Jackson et al., 2014). Introducing short (2–4 nt) randomized sequences at the 5′ and 3′ ends of the adaptors also reduced the adaptor ligation bias (Jayaprakash et al., 2011; Sorefan et al., 2012; Zhang et al., 2013).

Read mapping

Read mapping is the most crucial procedure. Although principally the ribosomal footprints are in their core an RNA-Seq data set, there is no standardized pipeline with recommended mapping parameters. Mapping can be performed to genomes or transcriptomes, but the short single-end reads generated in the ribosome profiling experiment cannot be used for de novo assembly of genomes or transcriptomes (Simpson and Pop, 2015). Mapping to the genome should be preferred as it is unbiased towards known exon and intron annotations and allows for discovery of previously undescribed ORFs (Andreev et al., 2015). Genome mapping usually gives greater coverage than mapping to transcriptomes (the loss of reads on exon junctions is minor) (Oshlack et al., 2010). Furthermore, genomes are better defined than transcriptomes, which are constructed in several different ways (reviewed in Garber et al., 2011). Also, mapping to genomes is less computationally intense and thus faster.

A prerequisite to good results is complete genome annotation, i.e. the availability of the gene coordinates. Genome annotation is a subject of intensive and constant improvement. For example the E.coli genome hosted on the NCBI server (Freddolino et al., 2012) is updated daily and the number of genes constantly changes. Although this fast adjustment makes new findings immediately available, it creates a gap with the hand-curated databases, some of which may offer more precise annotation of additional features. For example, RegulonDB (Salgado et al., 2013) offers more information on additional features than the NCBI annotations, including genes organized in operons, 5′ and 3′ UTRs. For eukaryotes the development is equally fast with frequently updated versions of genomes and their annotations. Three important webservers host various eukaryotic genomes: NCBI reference sequences, RefSeq (Pruitt et al., 2007), ensembl (Cunningham et al., 2015) and UCSC (Kent et al., 2002). The genome annotation choice may significantly influence the downstream quantification of expression and differential analysis (Zhao and Zhang, 2015), although a simple advice on which database to use is not possible and should be driven by the purpose of the analysis. For research aiming at reproducible and robust gene expression estimates, RefSeq might be preferred (Wu et al., 2013). More exploratory questions may rely on more complex annotations, e.g. ensembl.

The mapping tools can be classified into two major groups: hash-table based (Li et al., 2008; Homer et al., 2009); or Burrows-Wheeler Transform (BWT) algorithms (Langmead et al., 2009; Li and Durbin, 2009). While BWT-based approaches are faster and less computationally demanding, the hash-table-based algorithms are more flexible in aligning reads with non-perfect matches. Also the efficiency of BWT-based mapping approaches inversely correlates with the number of mismatches [reviewed in (Li and Homer, 2010; Garber et al., 2011)]. Comparison of the tools is not trivial and differs depending on the data set, thus only few objective investigations have been performed so far (Giannoulatou et al., 2014). In the majority of ribosome profiling experiments (Ingolia et al., 2009, 2011; Guo et al., 2010; Gerashchenko et al., 2012; Li et al., 2012; Chew et al., 2013; Guttman et al., 2013; Aspden et al., 2014; Baudin-Baillieu et al., 2014; Bazzini et al., 2014; Subramaniam et al., 2014), Bowtie (Langmead et al., 2009) is used as a BWT-based mapping program. Bowtie offers two ways of mapping a read to a reference sequence: seed- (parameter n) and mismatch-based approach (parameter v, Table 1). The seed approach aligns first a seed (or core) of a read and then extends the alignment further along the read length. Thereby, the mismatches in the seed count stronger than those in the extensions. Mostly, default Bowtie parameters (parameter n for the seed-based approach) are used (Guo et al., 2010; Li et al., 2012; Baudin-Baillieu et al., 2014; Subramaniam et al., 2014). Some studies apply the mismatch approach (Ingolia et al., 2011; Gerashchenko et al., 2012) which scores every base of each read equally. As the default seed length of 28 nt remains unchanged when using the default parameter settings, the seed-based strategy effectively works as a mismatch approach.

A general drawback of Bowtie is its inability to map splice junctions. One commonly used tool to align short reads across junctions is TopHat (Trapnell et al., 2009; Kim et al., 2013) which can also find junctions de novo. First, the TopHat pipeline maps to all reads to a reference genome using Bowtie and allows reporting more than one alignment of a read (i.e. m=inf k=20 [translated to Bowtie parameters]). TopHat then assembles the mapped reads using the assembly module in Maq (Li et al., 2008) in contiguous sequences inferring them to be putative exons, then uses seed and extended alignment to match reads to possible splice sites (Trapnell et al., 2009). The pipeline of TopHat is more structured, with fewer possibilities for changing the mapping parameters (Table 1), whereas Bowtie allows flexible adjustment of the mapping parameters. A new version of Bowtie, Bowtie2, has been launched (Langmead and Salzberg, 2012) which however differs conceptually from Bowtie and can find gapped alignments of reads resulting from insertions or deletions or sequencing errors. Note that Bowtie2 is suitable for reads longer than 50 nt.

In general, mapping can be defined as a procedure to find the unique position of each read in the reference genome (Oshlack et al., 2010). The logical consequence of this is to discard all reads with more than one best position. As ribosomal footprints are very short, the proportion of reads mapping at more than one position increases with the size of the genome. In the ribosome profiling datasets a large fraction of the reads map at multiple positions to the genome (in the range of 30% and more), however there is no uniform strategy on how to handle them. Strategies range from considering only reads uniquely mapping to the genome (Guo et al., 2010; Baudin-Baillieu et al., 2014), to allowing more than 200 or unlimited number of positions (Ingolia et al., 2011). One strategy to estimate the true location of reads that map to many (multimappable) locations involves proportional assignment of reads based on the read density of the neighboring positions (Trapnell et al., 2010). We varied the number of the mapping positions when aligning reads from ribosome profiling data (Figure 3A) and observed a clear difference between unique mapping (Figure 3B,C) and allowing multiple positions (Figure 3D,E). The number of the mapped reads increases when multiple mapping positions are allowed (Figure 3A). Such scenarios are relevant mostly to genes with duplications or highly homologous isoforms, nevertheless the fraction of the discarded non-uniquely mappable reads can be in the range of 30%. Stringent criterion leads to sparse and incomplete coverage of each gene (Figure 3B–C), while allowing mapping to multiple positions in the genome improves significantly the coverage of a single gene (Figure 3D–E). The assignment of the reads mapping to multiple positions is also of crucial importance. Multimappable reads can be assigned randomly to one of the possible positions they map (Figure 3D) or to all possible positions (Figure 3E). However, choosing the mapping parameter in such a way that the first hit position is reported (Figure 3D) bears some caveats as the origin of the reads is unclear, i.e. whether they are from the same gene or originate from another position in the genome. Thus, it might artificially increase the total number of reads on a gene. In this context, mapping to multiple positions with equal weighting of all positions (Figure 3E) might be a better choice as it does not prefer between positions and maps uniformly to all best mappable positions. For some analysis, to avoid overinterpretation of the data (for example by differential analysis), the most conservative mapping with uniquely mappable reads (Figure 3C) might be the best choice.

Figure 3:

Effects of different mapping strategies on gene coverage.

(A) Total reads mapped to the genome allowing different number of maximal best mappable positions per read (at m=1, a read with one best position (uniquely mappable) is allowed, otherwise discarded; at m=100 a read can be mapped to the 100 best positions and all positions are recorded). Mapping of ribosome profiling data of mouse brain was performed with Bowtie using the mouse genome (assembly GRCm38) allowing two mismatches per read (-v2 –strata –best). (B–E) Different mapping strategies result in variations in the read coverage of a human rplA gene mapped from ribosome profiling data with Bowtie using the human genome (assembly GRCm38). Mapping with (B) a single hit (uniquely mappable) fulfilling the restrictions of maximum two mismatches (-v2 –m1), (C) with a single best hit (uniquely mappable) (-v2 –m1 –strata –best), (D) with default parameter and restrictions to maximum 100 positions with two mismatches, but with only one listed in the output (-v2 –m100), (E) with multimapping restricted to 100 multimappable best positions (i.e. lowest number of mismatches) and best positions listed in the output (-v2 –m100 –k100). Best, the parameters strata best are given to ensure that a multimappable position is counted as such only by the same minimum number of mismatches; no best, default mode with no strata best parameters chosen; m, maximum number of multiple positions per read; k, maximum number of reported alignment. Note settings as in B should be avoided. Parameters as in C show the most reliable data, although the coverage is incomplete. The loss of reads in D compared to E is most likely due to the reporting of only one of the valid alignment positions.

The majority of the ribosome profiling datasets mapped with Bowtie do not set parameters to evaluate the quality of the alignments for a read (e.g. strata best), that compares, for example, whether a zero-mismatch mapping is better than an alignment with two mismatches. Usually, the first encountered alignment of a read is assigned to it (Gerashchenko et al., 2012; Li et al., 2012; Subramaniam et al., 2014). Thus, when multimapping is allowed, a read with zero mismatches in a certain position may also be mapped to a different position with two mismatches. Consequently, it creates a bias because the best alignment would not be satisfied, but a read is randomly assigned to one of the two positions independent of the number of the mismatches. The choice of the parameters for the mapping are of crucial importance as they can result in significant variations in the mapping and gene coverage profiles (Figure 3B–E). For reproducibility of the results it is advisable to clearly state the mapping parameter in each publication.

As both nucleases (RNase I for eukaryotic and MNase for bacterial systems) that are used to produce ribosomal footprints, cleave also rRNA, the rRNA reads comprise a large fraction of the sequencing reads, despite their removal in the experimental procedure. rRNA mapping and subtraction of those reads can be done in an extra round before or after mapping to the genome. Thereby, the mapping of the rRNA reads should be strict, i.e. allowing only a single mismatch.

In summary, mapping defines the shape of dataset to be used for further analysis and hence is a crucial step for which the parameters should be chosen carefully. In studies aiming at reproducible and robust gene expression estimates, uniquely mappable reads aligned to a reference genome (i.e. m=1 strata best) should be selected as they bear the lowest bias. However, for some genes (e.g. isoforms, duplicates), using only uniquely mapped reads may result in a partial coverage of a gene (Figure 3B,C); an incomplete coverage cannot be used to extract specific positions on which the ribosomes may pause or enrichment of reads over specific codons. For such analysis a multiple alignment of the multimappable reads (i.e. m=10 a strata best) might be chosen to ensure a maximal gene coverage. This parameter set bears drawbacks in analyzing the coverage of simultaneous expression of genes sharing large sequence identity (i.e. isoforms and duplicates). Such genes should be carefully assessed and might be separately compared only with their uniquely mappable reads or their expression should be confirmed with alternative methods (e.g. qRT-PCR).

Normalization of the read counts

Following mapping, read counts, also called gene counts, are collected and assigned to each gene or non-coding RNAs. Overlapping genes can be an issue here. As the ribosome profiling protocol is strand-specific, overlapping genes on different strands are well resolved. For genes overlapping on the same strand, as commonly observed in E. coli in which the coding sequence of the one gene falls into the end of the coding sequence of another, some read counting tools correct for this by randomly distributing the reads to the two overlapping genes (Anders et al., 2015), while other tools do not recognize overlapping features (Quinlan, 2014).

A commonly applied approach for normalization of the read counts is reads per kilobase of exon per million mapped reads, rpkM, (Mortazavi et al., 2008), which accounts for the differences in the sequencing depth (i.e. total number of the mapped reads) between sequencing libraries and for the length variation of each gene (i.e. per kilobase). Note that for short genes this normalization can give quite high rpkM values despite the presence of only few raw counts. Thus the detection limit should be set up using the raw counts (Ingolia et al., 2011). Other normalization approaches frequently applied in the RNA sequencing (RNA-Seq) might be applied too (Anders and Huber, 2010; Robinson et al., 2010; Dillies et al., 2013) but the statistical behavior of ribosome profiling data with those normalization procedures has not yet been tested (Olshen et al., 2013).

Further downstream analysis and post-processing

In the RNA-Seq datasets, several tools are used to identify differentially expressed (DE) genes (Guo et al., 2013), some of which (e.g. DESeq tool) have been applied in a few ribosome profiling studies (Baudin-Baillieu et al., 2014; Sidrauski et al., 2015). Still, they require a test that the ribosome profiling read counts follow the underlying distributions required by many tools designed for DE analysis of RNA-Seq, for example DESeq (Anders and Huber, 2010), EdgeR (Robinson et al., 2010), and baySeq (Hardcastle and Kelly, 2010). In all cases, a careful and conservative interpretation of the data is needed because, unlike RNA-Seq (Dillies et al., 2013), no uniform pipeline exists for ribosome profiling data. So far, only one tool has been developed specifically for ribosome profiling data (Olshen et al., 2013). Instead of performing DE analysis, a simple fold-change analysis can be carried out (Dunn et al., 2013) with the assumption that most of the genes are unchanged.

Still, a fascinating issue of ribosome profiling is the ability to record the position of ribosomes with single nucleotide resolution (Ingolia et al., 2009; Woolstenhulme et al., 2015), which enables detecting ribosomal pausing (i.e. specific positions at which ribosomes pause) or encoding events (e.g. readthrough or frameshifting) (Li et al., 2012; Michel et al., 2012; O’Connor et al., 2013). The alignment of the ribosomal reads to the open-reading frame is called calibration in which the start codon is assigned to the ribosomal P-site (Ingolia et al., 2009) or the stop codon is assigned to the A-site (Woolstenhulme et al., 2015). If the ribosomes are not completely halted during the isolation procedure, it will compromise the calibration and would not allow for codon resolution (Ingolia et al., 2009). While ribosomal footprints of eukaryotic ribosomes can be calibrated using both stop and start codons, i.e. both 5′ and 3′ of the reads, reads from bacterial systems give only codon resolution when calibrated using their 3′ ends, most likely because of the sharp cleavage of the MNase at the 3′ of the reads but not at the 5′ ends (Woolstenhulme et al., 2015). Another approach to gain positional information of the translating ribosomes is center-weighted or center-assigned approach (Li et al., 2012). A defined number of nucleotides are excluded from both 5′ and 3′ sides of a read and the remaining centrally positioned nucleotides are weighted equally. This approach delivers less sharp resolution and defines the position of the ribosomal A- or P-sites with a subcodon resolution. Thus, it has limited applications and cannot be used for determining the reading frame (Woolstenhulme et al., 2015). Both, calibrated and center-weighted ribosomal reads can be used to assess ribosomal enrichment over specific codons (Li et al., 2012; Ishimura et al., 2014) or to determine sequences over which ribosomes transiently pause (Li et al., 2012; Woolstenhulme et al., 2015).

In the library preparation, RNA fragments over a typical length range of 25–35 nt, tightly distributed around a peak of ~28 nt, are selected from the gel upon ribonucleolytic digestion (Ingolia et al., 2009; Guydosh and Green, 2014). It should be noted that reads outside this range may also bear some biological information and, dependent on the specific question, might also be included in the library preparation. Reads shorter than the average length of ~28 nt represent different conformational states of the elongating ribosome (Lareau et al., 2014) or report on ribosomes stalled over 3′ truncated mRNAs (Guydosh and Green, 2014). In turn, longer reads may be informative on frameshifting events (O’Connor et al., 2013). When comparing expression level on a gene basis in the DE analysis, all reads independent of their length might be considered under the assumption that each ribosome read produces one protein. For more specific analysis, including ribosomal stalling at specific positions, the reads should be separated by their length and each length group should be treated separately.

Computational demand and infrastructure

Raw data from one sequencing lane of Illuminas HiSeq machine (Illumina, San Diego, USA) can reach a size of more than 20 GB (uncompressed). Preprocessing and mapping of these raw files easily exceeds another 20 GB; discarding the intermediate preprocessing file and keeping only compressed raw files requires hard disk space for one lane of about 20 GB. The demand of RAM varies dependent on the type of analysis and programming languages. For example, using a simple Perl hash index build on each of the ~4 million nucleotides of the relatively small E. coli genome requires more than 4 GB of RAM. Mapping with BWT-based algorithms demands relatively low memory (Langmead et al., 2009; Li and Durbin, 2009). For example, the human genome can be mapped with <8 GB of RAM (Langmead et al., 2009; Li and Durbin, 2009). The mapping programs offer an option to use more than one CPU in parallel to increase speed (Langmead et al., 2009). Many of pre- and post-processing steps are not implemented as full programs but as a collection of scripts or even in-house scripts (Anders et al., 2015).