RiboProfiling: a Bioconductor package for standard Ribo-seq pipeline processing

The ribosome profiling technique (Ribo-seq) allows the selective sequencing of translated RNA regions. Recently, the analysis of genomic sequences associated to Ribo-seq reads has been widely employed to assess their coding potential. These analyses led to the identification of differentially translated transcripts under different experimental conditions, and/or ribosome pausing on codon motifs. In the context of the ever-growing need for tools analyzing Ribo-seq reads, we have developed ‘RiboProfiling’, a new Bioconductor open-source package. ‘RiboProfiling’ provides a full pipeline to cover all key steps for the analysis of ribosome footprints. This pipeline has been implemented in a single R workflow. The package takes an alignment (BAM) file as input and performs ribosome footprint quantification at a transcript level. It also identifies footprint accumulation on particular amino acids or multi amino-acids motifs. Report summary graphs and data quantification are generated automatically. The package facilitates quality assessment and quantification of Ribo-seq experiments. Its implementation in Bioconductor enables the modeling and statistical analysis of its output through the vast choice of packages available in R. This article illustrates how to identify codon-motifs accumulating ribosome footprints, based on data from Escherichia coli.

Ribosome profiling (Ribo-seq) is a recently developed high throughput sequencing technique (Ingolia et al., 2009) that allows the identification of RNA fragments resistant to RNAse digestion. Fragments mainly correspond to coding sequences protected against RNAse digestion by translating ribosomes. Ribo-seq data have been widely used to assess the translational status of open reading frames (ORFs) (Bazzini et al., 2014;Fields et al., 2015;Ingolia et al., 2009;Popa et al., 2016), and to identify ORFs differentially translated between experimental conditions (Schafer et al., 2015).
Ribo-seq bioinformatics analyses comprise the selection of reads consistent with ribosome footprints, a recalibration of the read start or end (5' or 3' extremity of the read depending on the RNA digestion step) to the peptidyl site (P-site) position of the ribosome, and quantification of the reads on specific features of interest (i.e. transcript, codons, multi-codons motifs). Several tools for processing ribosome profiling data have previously been proposed. RiboTools ' (Legendre et al., 2015) and 'RUST' (O'Connor et al., 2015)) were developed in python, 'riboSeqR' (Hardcastle, 2014) corresponds to an R package. Each of the above software integrates some, but not all, of the functions necessary for a standard Ribo-seq workflow from reads to quality assessment, recalibration and quantification. A previous effort to group the different approaches for Ribo-seq analyses has been developed with the Galaxy instance RiboGalaxy (Michel et al., 2016).
These tools have been developed to answer specific questions related to ribosome occupancy: normalization of Ribo-seq reads ('RUST'), detection and characterization of reading frame usage ('riboSeqR'), or irregular translational behavior such as translational ambiguities (Ribo-seq footprints in different phases) and stop-codon read through ('RiboTools'). However, prior to any advanced Ribo-seq data processing, it is necessary to have standard pipelines for quality assessment of the experiments and specific ribosome footprint assignment to sequences.
We implemented the features for a standard Ribo-seq workflow in a pipeline entitled 'RiboProfiling'. The pipeline takes an alignment file (BAM) as input, performs identification of the read offset, generates transcript and (multi-) codon coverage quantification data, and performs statistical analyses as well as graphical representations. Our pipeline is, to our knowledge, the most complete integration of a ribosome profiling standard analysis pipeline in an unique R framework. This includes the crucial step prior to quantifying ribosome footprints that consists in identifying the offset between Ribo-seq reads and the P-site of the ribosome. We have given special attention to this step as it is essential to correctly associate ribosome footprints with codon resolution. Depending on the RNA digestion protocol, the assignment of the ribosome must be made specifically either to the 5', 3' or the center of the read. To our knowledge, there is only one implementation for the determination of ribosome offset that was published so far. It was proposed in 'riboSeqR' as a metagene plot and the determination of the offset was only possible from the 5' read end and for a particular read length. Our package offers several options, to compute the offset and recalibrate reads based either on the 5' or the 3' read ends. RiboProfiling also enables the graphical representation of ribosome density around the Translation Start Site (TSS) for multiple read lengths. This option allows to perform analysis for a single read length or on the merge of several read lengths. It enables to group all lengths sharing a same and unique offset value. The R/Bioconductor implementation provides an easy-to-use comprehensive set of functions that requires a minimal knowledge in R programming. The package contains a function entitled 'riboSeqFromBAM' that treats multiple Ribo-seq BAM files in parallel. The automated workflow generates report summary graphs and data quantification.
We illustrate the main features of the package using a Ribo-seq control sample from murine ES cells (GSM1655059), taken from a recently published ribosome profiling study using translation inhibitors (Popa et al., 2016). We then detail the analysis of ribosome accumulation on certain codons and tri-peptide motifs on a public dataset in Escherichia coli (GSE64488) (Woolstenhulme et al., 2015). The script for performing all these analyses is publicly available.

Requirements
The 'RiboProfiling' package, v.1.2.0 can be used with R 3.3.0 and Bioconductor version 3.3. The script needs a minimum memory limit in R of 3 Gb when analyzing tripeptide motifs. The package starts from alignment BAM files, from either Ribo-seq or RNA-seq experiments. We have validated BAM files from bowtie/ tophat, Hisat2, STAR, and Lifescope (Solid), both single-or paired-end (for RNA-seq reads). Reads from rRNA, tRNA, and PCR duplicates (if unique molecular identifiers are available) should be removed from the BAM files before starting the analysis (see package vignette for details).

RiboProfiling package
All analyses can be performed either through a call to a function called 'riboSeqFromBAM', or through a step by step approach. Figure 1 describes the workflow of the package starting from BAM files with reads mapped to the genome of interest. The first step in processing Ribo-seq reads is to select only those with match lengths compatible with standard ribosome footprints. The function 'histMatchLength' allows the visual inspection of read match sizes distribution across the BAM file ( Figure 2.A), which should be enriched in reads of sizes between 20 to 40 nucleotides (Ingolia et al., 2009;Popa et al., 2016).
A second particularity in the handling of Ribo-seq data comes from the shift existing between the extremities of the read (i.e. 5' or 3') and the P-site position of the ribosome. Our package allows the identification of an offset from the 5' end of the read, but also from the 3' end. The function 'readStartCov' computes the read frequency distribution centered on the translation start site (TSS) of the most expressed protein coding transcripts (by default the 3% most expressed). Based on this frequency distribution, the 'plotSummarizedCov' function enables the visual quantification of the offset between the reads and the ribosome P-site (Figure 2.B).  In our Ribo-seq example, the 5' read end is shifted 13 bp from the TSS. The innovation of this feature consists in the visualization of read lengths independently and as a summary figure. When computed, the offset can be applied on all reads based on the transcript referential with the function 'countShiftReads' and coverage on three different sequence features: 5'-UTR, coding sequences (CDS), and 3'-UTR (Figure 3.A). In Figure 2.C we observe that the majority of reads accumulate on the CDS of protein coding sequences and are practically lacking in the 3' non-coding UTR regions. Ribosome footprints are also detectable in the 5' UTR regions of protein coding genes, suggesting either the presence of coding upstream ORFs (Popa et al., 2016) or possible confounding information from missing annotations.
Finally, our package provides quantification of ribosome footprints at codon resolution (Figure 3.B). A PCA analysis of codon occupancy can be performed and several graphical functions are implemented. In Figure 2.D we employed the 'codonInfo' and 'codonPCA' functions to analyze the codons accumulating ribosome footprints. As expected, the codon ATG is the most discriminant codon in the PCA analyses, since ribosome accumulation peaks are observed at the start codon of coding regions. Detailed descriptions with examples of the pipeline from BAM files to Ribo-seq reads quantification and processing are available in the vignette of our package: https://www.bioconductor.org/packages/ release/bioc/vignettes/RiboProfiling/inst/doc/RiboProfiling.pdf.

Analysis of ribosome stalling on sequence motifs
'RiboProfiling' can also be useful for analyzing ribosome occupancy on multi-codons motifs. Codons accumulating ribosome footprints are indicative of slowed ribosome progression (stalling) during the translation elongation process. The 'RiboProfiling' package offers several features for quantifying footprint accumulation on sequence motifs (ranging from one to three consecutive codons), performs principal component analyses, and allows graphical representation of those data.
To illustrate how 'RiboProfiling' can be used to explore the influence of sequence motifs (in this case tri-amino-acid sequences) on ribosome pausing, we analyzed an Escherichia coli Ribo-seq dataset (Woolstenhulme et al., 2015). We downloaded, filtered and mapped the reads of an efp-knockout sample (ΔEfp2, GSE64488), the elongation factor EFP being essential for the translation of polyproline motifs. Uniquely mapped reads from the resulting BAM file were analyzed with our 'RiboProfiling' package. After quality assessment of the reads size distribution, we quantified the offset between the 3' end of the reads and the TSS for different alignment match sizes (Figure 4.A). We can clearly observe the 15 nucleotides offset that was reported by the authors for reads with alignment sizes >= 30 nucleotides. Smaller match lengths exhibited either a strong variation in the distribution of reads around the TSS (i.e. the 29 mers), or a different offset (i.e. offset of 20 for 21mers) (Figure 4.A). We selected the reads with alignment match sizes between 30 and 40 nucleotides and quantified codon coverage by positioning the ribosome P-site 15 nucleotides upstream of the read 3' extremity.
An important stalling has been reported after the incorporation of two consecutive prolines (Pro-Pro) in the peptide chain. This stalling is highly dependent on the nature of the codon that follows the aminoacyl-tRNA reacting at the A-site (Woolstenhulme et al., 2015). Following the article's analysis of stalling at PP(X) (Proline -Proline -3rd codon) motifs, we used the 'countShiftReads' and 'codonInfo' functions to quantify the ribosome footprints on these motifs, in the ΔEfp2 sample. We then where, Reads PPX|ORF is the ribosome density of motif PPX in a given ORF; Reads ORF is the ribosome density on the ORF; Nbr PPX|ORF is the number of time a given PPX motif is present in the ORF; Length ORF is the total length of the ORF. We averaged the ribosome occupancy for each PP(X) motif on all ORFs. Figure 4.B shows a strong stalling in sample ΔEfp2 when the ribosome encounters PPN, PPW, PPD, in agreement with their previous identification as pause sites in Escherichia coli (Woolstenhulme et al., 2015). A step by step R script implementing this entire analysis is provided at: http://genomique.info/data/public/RiboProfiling/scriptWoolsten-hulme_Defp2.R.

�ummary
Our 'RiboProfiling' Bioconductor package offers a collection of tools for Ribo-seq data analysis. It provides an unique, straightforward R implementation of a ribosome profiling pipeline from BAM, to P-site calibration, quantification of reads on sequence features, and codon coverage. The packages' graphical features offer quality assessment and result representation across the analyses. Following the overview of Ribo-seq experiments with 'RiboProfiling, the output tables can then be easily integrated into more specialized dowstream analyses, either using more specialisez riboseq tools such as (XXX, YYY) or directly within R.
We highlighted here the features of our package in characterizing ribosome stalling at sequence motifs along ORFs based on an example dataset from Woolstenhulme et al. (Woolstenhulme et al., 2015). The workflow we propose for the analysis of ribosome occupancy on codon motifs using the 'RiboProfiling' package will most surely prove an useful asset in the context of recent ribosome profiling applications such as the detection of tumor sensitivity to differential amino acid depletion (Loayza-Puch et al., 2016).
In the vignette page 5, I would add the statement library(GenomicAlignments) before the readGAlignments call in order to make the code easier to work.. Vignete page 6. The instruction txdb <-TxDb.Hsapiens.UCSC.hg19.knownGene does not work. I do not have a TxDb.Hsapiens.UCSC.hg19.knownGene object in my session. I couldn't go further in running the provided tutorial. Could you clarify on how to obtain that object ?
Apart from the vignette, the paper itself is well written and can be followed simply through. An analysis case of ribosome stalling study is presented through a complete example and R code (analysis of Woolstenhulme et al dataset). The code works without problem until the invocation of library(BSgenome.Ecoli.NCBI.K12.MG1655). I couldn't get this package on the Bioconductor Web Site, and some clarification might be helpful for this step.
Minor remark: The authors missed a couple of references in the last paragraph: (...) either using more specialisez riboseq tools such as (XXX, YYY) or directly within R.

I have read this submission. I believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.
No competing interests were disclosed. RiboProfiling allows a number of ribo-seq specific analyses to be carried out using a single Bioconductor package in R and is supported by instructions in the corresponding reference and 'package' manuals. The automatic assignment of the P-site from either the 5' or 3' end of the ribosome footprint read is a useful feature. In addition the tool allows the quantification of the distribution of footprints across coding regions (CDS), 5' leader and 3' trailer regions at the metagene and individual gene transcript level. The authors use the reported stalling events in the Woolstenhulme 2015 ribo-seq data to illustrate how et al., RiboProfiling can be used to determine the ribosome occupancy on multiple codon motifs.
Overall, the RiboProfiling package allows some ribo-seq specific analyses to be automatically carried out using the RiboProfiling 'quick command' or by using the step-by-step commands. However I would argue that RiboProfiling is not more complete than existing tools. There is no assessment of triplet periodicity for example. The quality assessment carried out by RiboProfilig is the distribution of footprint read lengths ( Figure 1). While determining the read lengths in a ribo-seq dataset is a useful assessment, it is only one 1 F1000Research data and this package generates a great summary of these. However this package is not as complete as suggested by the authors. Indeed an important point is the periodicity of footprints that is not addressed here. Moreover it could be interesting to add the ability to detect which reading phase is indeed read. This should be feasible since the P-site is well defined. I include below a few comments for the authors.
The algorithm of the readStartCov function is not sufficiently described. The only one sentence found, does not give enough information to understand how it works. For instance, "the read frequency distribution centered on the translation start site": "Centered" needs more explanation. This is important, because using our own datasets we were not able to obtain any convincing results about the position of the P-site. The lack of information prevented us from understanding why.
In the Package Manual, when the aroundPromoter function is defined, the TSS is described as the "Transcript Start Site", which does not match with the "Translational Start Site" found in the article.
The pipeline requires an annotation file with a specific format (TxDb object). Users have to be aware of that if they want to use their own genome.
We have read this submission. We believe that we have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.
No competing interests were disclosed. Competing Interests: