Modeling Ribosome Footprinting Data

In sequencing-based ribosome footprinting, the RF read count is naturally confounded by mRNA abundance (Fig. 1A). We seek a strategy to compare RF measurements taking mRNA abundance into account in order to accurately discern the translation effect in case–control experiments. We model the vector of RNA-Seq and RF read counts $y_{\text{mRNA}}^i$ and $y_{\text{RF}}^i$ , respectively, for gene i with Negative Binomial (NB) distributions, as described before (for instance, Love et al., 2014 (link); Drewe et al., 2013 (link); Robinson et al., 2010 (link)): $y^i\hspace{0.17em}\sim \hspace{0.17em}NB({\mu }^i,{\kappa }^i),$ where μⁱ is the expected count and κⁱ is the estimated dispersion across biological replicates. Here yⁱ denotes the observed counts normalized by the library size factor (Supplementary Section A). Formulating the problem as a generalized linear model (GLM) with the logarithm as link function, we can express expectations on read counts as a function of latent quantities related to mRNA abundance β_C in the two conditions ( $C=\{0,1\}$ ), a quantity ${\beta }_{\text{RNA}}$ that relates mRNA abundance to RNA-Seq read counts, a quantity ${\beta }_{\text{RF}}$ that relates mRNA abundance to RF read counts and a quantity ${\beta }_{\Delta ,C}$ that captures the effect of the treatment on translation. In particular, the expected RNA-Seq read count ${\mu }_{\text{mRNA},C}^i$ is given by the equation $\text{log}({\mu }_{\text{mRNA},C}^i)={\beta }_C^i+{\beta }_{\text{RNA}}^i$ .
We assume that transcription and translation are successive cellular processing steps and that abundances are linearly related. The expected RF read count, ${\mu }_{\text{RF},C}^i$ , is given by $\text{log}({\mu }_{\text{RF},C}^i)={\beta }_C^i+{\beta }_{RF}^i+{\beta }_{\Delta ,C}^i$ . A key point to note is that ${\beta }_C^i$ is revealed to be a shared parameter between the expressions governing the expected RNA-Seq and RF counts. It can be considered to be a proxy for shared transcriptional/translation activity under condition C in this context. Then, ${\beta }_{\Delta ,C}^i$ indicates the deviation from that activity under condition C, with ${\beta }_{\Delta ,C}^i=0$ for C = 0 and free otherwise (See Supplementary Section B for more details).
Fitting the GLM consists of learning the parameters βⁱ and dispersions κⁱ given mRNA and RF counts for the two conditions $C=\{0,1\}$ . We perform alternating optimization of the parameters βⁱ given dispersions κⁱ and the dispersion parameters κⁱ given βⁱ, similar to the EM algorithm (Supplementary Sections B and C):
${\beta }^i=\underset{{\beta }^i}{\text{arg}\hspace{0.17em}\text{max}}{\ell }_{glm}({\beta }^i|y^i,{\kappa }^i)\hspace{1em}\text{and}\hspace{1em}{\kappa }^i=\underset{{\kappa }^i}{\text{arg}\hspace{0.17em}\text{max}}{\ell }_{NB}({\kappa }^i|y^i,{\mu }^i).$
As experimental procedures for measuring mRNA counts and RF counts differ, we enable the estimating of separate dispersion parameters for the data sources of RNA-Seq and RF profiling to account for different characteristics (Supplementary Section E).
As in Anders et al. (2012) (link), with raw dispersions estimated from previous steps, we regress all κⁱ given the mean counts to obtain a mean-dispersion relationship $f(\mu )={\lambda }_1/\mu +{\lambda }_0$ . We perform empirical Bayes shrinkage (Love et al., 2014 (link)) to shrink κⁱ towards $f(\mu )$ to stabilize estimates (see Supplementary Section D). The proposed model in RiboDiff with a joint dispersion estimate is conceptually identical to using the following GLM design matrix $\text{protocol}+\text{condition}+\text{condition}:\text{protocol}$ (for instance, in conjunction with edgeR or DESeq1/2).
In a treatment/control setting, we can then evaluate whether a treatment (C = 1) has a significant differential effect on translation efficiency compared to the control (C = 0). This is equivalent to determining whether the parameter ${\beta }_{\Delta ,1}$ differs significantly from 0 and whether the relationship denoted by the dashed arrow in Figure 1A is needed or not. We can compute significance levels based on the ${\chi }^2$ distribution by analyzing $\text{log}$ -likelihood ratios of the Null model ( ${\beta }_{\Delta ,1}^i=0$ ) and the alternative model ( ${\beta }_{\Delta ,1}^i\cancel{=}0$ ).