MGMR: leveraging RNA-Seq population data to optimize expression estimation

Abstract

Background

RNA-Seq is a technique that uses Next Generation Sequencing to identify transcripts and estimate transcription levels. When applying this technique for quantification, one must contend with reads that align to multiple positions in the genome (multireads). Previous efforts to resolve multireads have shown that RNA-Seq expression estimation can be improved using probabilistic allocation of reads to genes. These methods use a probabilistic generative model for data generation and resolve ambiguity using likelihood-based approaches. In many instances, RNA-seq experiments are performed in the context of a population. The generative models of current methods do not take into account such population information, and it is an open question whether this information can improve quantification of the individual samples

Results

In order to explore the contribution of population level information in RNA-seq quantification, we apply a hierarchical probabilistic generative model, which assumes that expression levels of different individuals are sampled from a Dirichlet distribution with parameters specific to the population, and reads are sampled from the distribution of expression levels. We introduce an optimization procedure for the estimation of the model parameters, and use HapMap data and simulated data to demonstrate that the model yields a significant improvement in the accuracy of expression levels of paralogous genes.

Conclusions

We provide a proof of principal of the benefit of drawing on population commonalities to estimate expression. The results of our experiments demonstrate this approach can be beneficial, primarily for estimation at the gene level.

Introduction

With the rapid decline in the cost of sequencing, RNA-Seq has emerged as a legitimate competitor to mi-croarrays as a means of assessing global gene expression. Even as arrays currently enjoy a cost advantage, many new applications of information accessible only through sequencing further strengthen the case that sequencing may soon supplant arrays as the technology of choice for transcription analysis. One such application is fine-grained assessment of variation in expression and the sources for such variation, as exemplified by recent large-scale RNA-Seq studies [1, 2] of two different HapMap [3] populations. Such studies complement genomic DNA sequencing by elucidating the link between SNPs and expression.

Unfortunately, with any new technology come its limitations, and several studies have pointed out that RNA-Seq is not immune to bias [4–6]. Perhaps the most widely discussed hurdle to accurate estimation in the case of RNA-Seq is the problem of reads mapped to multiple locations in the target genome (or in the target transcript sequences). These reads, which are called multireads, can stem from either paralogous gene sequences or from different isoforms of the same gene that share exons.

Several methods have emerged to address the multiread problem for paralog and isoform estimation [7–10]. These methods are all based on probabilistic modeling that is optimized by an expectation maximization procedure. It has been repeatedly shown that by using such methods one can get better quantification of the expression levels compared to quantification based on naive approaches of read assignment.

In many applications, a set of samples is studied. For instance, one may be interested in comparing the expression levels in cases verses controls, or in tissues originating from different organs. In such cases, it is plausible that the commonality of expression patterns within each of the defined groups of studied samples may be used to improve the quantification results in each of the samples. We demonstrate that by analyzing expression profiles of a population together, one gets expression estimates more accurate than those obtained by estimating individual sample expression levels independently. We describe and implement a probabilistic model of the sequencing process that incorporates population commonalities, and demonstrate its advantages over existing methods in the population setting.

Methods

RNA-Seq multiread expression estimation

As we seek to extend the prevalent generative model of RNA-Seq [7–11], we begin by reviewing the basic elements of that model. Let G = (G₁, ...,G _M ) be the set of M transcribed regions considered and P = (P₁, ...,P _M ) be the proportions of RNA bases attributed to each transcript out of the total number of transcribed bases in a sequenced sample. Regions may be either genes or transcripts, depending on the level of transcription being investigated. We require P to satisfy ∑ _gϵG P _g = 1 and ∀gϵG, 0 ≤ P _g ≤ 1.

The model describes an RNA sequencing experiment where regions in G are randomly chosen according to the distribution P, start positions in these regions are chosen uniformly, and reads of length ℓ are generated by copying ℓ consecutive bases from each chosen region to produce a set of reads R = (r₁,..., r _ρ ). Sequencing is assumed to be error prone, leading to a certain probability of error for each read base. Based on the repetitions present in the set of regions and errors in alignment, reads may fail to map to the region from which they originate or may map to additional locations. Thus, we assign a probability of obtaining read r _j given that it originated from region G_k,P\left(r_j|G_k\right)\equiv \frac{{\left(1-\epsilon \right)}^{\left(\ell -erro{r}_{jk}\right)}{\epsilon }^{erro{r}_{jk}}}{{\ell }_k}. In this case we rely on the alignment of r _j to G _k to afford us the best match position instead of summing over all possible starting positions. ℓ _k is the effective length of G _k (i.e., the number of start positions from which a full length read can be derived) as defined in [11], ϵ is taken to be a constant per-base error rate, errors are assumed to be independent, and error _jk is the number of mismatches in the best alignment of r _j to G _k .

This formulation leads to the likelihood of observing the data:

L\left(P;R\right)=\prod _{j=1}^{\rho }P\left(r_j|G,P\right)=\prod _{j=1}^{\rho }\sum _{k=1}^{M}P\left(G_k\right)P\left(r_j|G_k\right)

(1)

This likelihood function is used to estimate P given the read alignments. Typically, one will use expectation maximization to find the P for which the likelihood is maximized. It is assumed that P(r _j |G _k ) is zero for all regions to which r _j is not aligned.

Common population extension

To estimate expression levels in N individuals from a defined population, we modify the above model by assuming that samples are drawn from a common population. This is imposed by having P = [(P₁₁, ...,P_1M), .., (P_{N 1}...,P _{N M} )] be probability densities drawn from a common Dirichlet distribution, defined by a set of hyper-parameters specific to the population: ∀iϵ[1, N], p _i = (P_{i 1},..., P _iM ) ~ Dirichlet (α₁,..., α _M ).

For sample i, we denote the set of reads as R_i=\left(r_{i1},\dots ,r_{i{\rho }_i}\right), where each r _ij is mapped to one or more regions in G. The output of a read alignment program defines the set of accepted regions for the read (in practice only alignments with up to 2 errors are accepted) and provides the number of errors in alignment for each read-region pair. This allows us to calculate P(r _ij |G _k ) as done above for one sample. For convenience we denote P(r _ij |G _k ) = q _ijk (taken to be zero for all regions not mapped to), which is independent of α and P.

As before, our objective is to estimate P, but in this case we must optimize by estimating P and α together. We begin by writing the likelihood function:

L\left({\mathbf{p}}_{\mathbf{1}},\dots ,{\mathbf{p}}_{\mathbf{N}},\alpha ;R\right)=Pr\left({\mathbf{p}}_1,\dots ,{\mathbf{p}}_{\mathbf{N}}|\alpha \right)Pr\left(R|{\mathbf{p}}_1,\dots ,{\text{p}}_{\mathbf{N}}\right)

(2)

Since expression values are sampled from the Dirichlet distribution,

Pr\left({\mathbf{p}}_{\mathbf{1}},\dots ,{\mathbf{p}}_{\mathbf{N}}|\alpha \right)=\prod _{i=1}^{N}P\left({\mathbf{p}}_i|\alpha \right)=\prod _{i=1}^{N}C\left(\alpha \right)\prod _{k=1}^M{P}_{ik}^{{\alpha }_k-1}

(3)

Where

C\left(\alpha \right)\equiv \frac{\Gamma \left(\sum _k{\alpha }_k\right)}{\prod _k\Gamma \left({\alpha }_k\right)}

(4)

and similar to (1) above,

Pr\left(R|{\mathbf{p}}_{\mathbf{1}},\dots ,{\mathbf{p}}_{\mathbf{N}}\right)=\prod _{i=1}^N\prod _{j=1}^{{\rho }_i}\sum _{k=1}^M{P}_{ik}{q}_{ijk}

(5)

This leads to

L\left({\mathbf{p}}_1,\dots ,{\mathbf{p}}_{\mathbf{N}},\alpha ;R\right)=\left[\prod _{i=1}^{N}C\left(\alpha \right)\prod _{k=1}^M{P}_{ik}^{{\alpha }_k-1}\right]\left[\prod _{i=1}^N\prod _{j=1}^{{\rho }_i}\sum _{k=1}^M{P}_{ik}{q}_{ijk}\right]

(6)

Taking the log, we get

\begin{array}{ll}\hfill log\left[L\left(p_1,\dots ,p_N,\alpha ;R\right)\right]& =NlogC\left(\alpha \right)+\sum _{k=1}^M\left({\alpha }_k-1\right)\sum _{i=1}^{N}log{P}_{ik}\\ +\sum _{i=1}^N\sum _{j=1}^{p_i}\left(log\left(\sum _{k=1}^M{P}_{ik}{q}_{ijk}\right)\right)\end{array}

(7)

Multi-Genome Multi-Read (MGMR) algorithm

We wish to estimate α and p ₁ ,...,p _n to maximize equation (7) above. For this purpose, we adopt an alternating iterative procedure of estimating α given the current estimate of p ₁ ,...,p _n and vice-versa until the total change in α becomes sufficiently small (or until a pre-set number of iterations have been executed).

Although for EM-based estimation methods convexity guarantees an optimal solution will be obtained, here (as shall be seen below) we have no such guarantee. Thus, we confine updates to be local by performing only one update for P and one for α. By one MGMR iteration, we refer to one EM-based P update followed by one α update.

Estimating P given α

If we assume α is given, we can write the EM steps to find p ₁ ,...,p _N :

E step Letting Match signify a matching between reads and regions, and Match(j) be the region from which read j originates, we get:

\begin{array}{ll}\hfill log\left[L\left(\mathbf{P},\alpha ;\mathbf{R},\mathbf{Match}\right)\right]& =NlogC\left(\alpha \right)+\sum _{k=1}^M\left({\alpha }_k-1\right)\sum _{i=1}^{N}log{P}_{ik}\\ +\sum _{i=1}^N\sum _{j=1}^{{\rho }_i}\left(log\left(P_{iMatch\left(j\right)}{q_{ij}}_{Match\left(j\right)}\right)\right)\end{array}

(8)

which leads to

Q\left(\mathbf{P},\alpha |{\mathbf{P}}^{\left(t\right)},{\alpha }^{\left(t\right)}\right)=E_{Match|R,{\mathbf{P}}^{\left(t\right)},{\alpha }^{\left(t\right)}}\left[log\left(L\right)\right]

(9)

\begin{array}{ll}\hfill =NlogC\left({\alpha }^{\left(t\right)}\right)& +\sum _{k=1}^M\left({\alpha }_k^{\left(t\right)}-1\right)\sum _{i=1}^{N}log{P}_{ik}\\ +\sum _{i=1}^N\sum _{j=1}^{{\rho }_i}\sum _{k=1}^M\left(log{P}_{ik}+log{q}_{ijk}\right)*\frac{P_{ik}^{\left(t\right)}{q}_{ijk}}{{\sum }_{k=1}^M{P}_{ik}^{\left(t\right)}{q}_{ijk}}\end{array}

(10)

\begin{array}{ll}\hfill =NlogC\left({\alpha }^{\left(t\right)}\right)+\sum _{k=1}^M\left({\alpha }_k^{\left(t\right)}-1\right)\sum _{i=1}^{N}log{P}_{ik}& +\sum _{i=1}^N\sum _{j=1}^{{\rho }_i}\sum _{k=1}^M{a}_{ijk}log{P}_{ik}\\ +\sum _{i=1}^N\sum _{j=1}^{{\rho }_i}\sum _{k=1}^M{a}_{ijk}log{q}_{ijk}\end{array}

(11)

where

a_{ijk}\equiv \frac{P_{ik}^{\left(t\right)}{q}_{ijk}}{{\sum }_{k=1}^M{P}_{ik}^{\left(t\right)}{q}_{ijk}}

(12)

M step Given that each q _ijk is fixed, the above reduces to maximizing

\begin{array}{ll}\hfill NlogC\left({\alpha }^{\left(t\right)}\right)& +\sum _{k=1}^M\left({\alpha }_k^{\left(t\right)}-1\right)\sum _{i=1}^{N}log{P}_{ik}\\ +\sum _{i=1}^N\sum _{j=1}^{{\rho }_i}\sum _{k=1}^M{a}_{ijk}*log{P}_{ik}\end{array}

(13)

Using Lagrange multipliers and differentiating, we see that this is maximized with

P_{ik}^{\left(t+1\right)}=\frac{{\alpha }_k^{\left(t\right)}-1+{\sum }_j{a}_{ijk}}{{\sum }_k\left({\alpha }_k^{\left(t\right)}-1+{\sum }_j{a}_{ijk}\right)}

(14)

Estimating α given P

Given a new estimate for P^(t), we can use a fixed point iteration [12] to get a new estimate of α

\begin{array}{ll}\hfill F\left(\alpha \right)& =N\left[log\Gamma \left(\sum _k{\alpha }_k\right)-\sum _{k}log\Gamma \left({\alpha }_k\right)\right]\\ +\sum _{k=1}^M\left({\alpha }_k-1\right)\sum _{i=1}^{N}log{P}_{ik}^{\left(t\right)}+Const\left(P^{\left(t\right)}\right)\end{array}

(15)

By using the known bound \Gamma \left(x\right)\ge \Gamma \left(\widehat{x}\right)exp\left(\left(x-\widehat{x}\right)\Psi \left(x\right)\right)\left(\text{having}\Psi \left(x\right)=\frac{dlog\Gamma \left(x\right)}{dx}\right), we can get a lower bound on F(α):

\begin{array}{l}F(\alpha )\ge N[({\displaystyle \sum _k{\alpha }_k)}\Psi ({\displaystyle \sum _k{\alpha }_k^{(t)}})-{\displaystyle \sum _{k=1}^{M}log\Gamma ({\alpha }_k)}\\ +{\displaystyle \sum _{k=1}^M{\alpha }_{k}log{\overline{P}}_k^{(t)}}+Const(P^{(t)})\end{array}

(16)

where log{\bar{P}}_k=\frac{1}{N}*{\sum }_{i}log{P}_{ik}.

We maximize this bound with a fixed point iteration similar to EM, noting that for fixed values of P convergence is guaranteed, and that for the Dirichlet distribution the maximum is the only stationary point [12]. This leads to the update

{\alpha }_k^{\left(t+\right)}={\Psi }^{-1}\left[\Psi \left(\sum \underset{k}{\overset{\left(t\right)}{\alpha }}\right)+log\underset{k}{\overset{\left(t\right)}{\bar{P}}}\right]

(17)

Heuristics/Implementation

As we have found F(α) presented in equation (15) is non-convex even in 2 dimensions (Figure 1), we confine updates to be local by allowing only one update for both the α and P estimation steps at each MGMR iteration. For genes with EM expression estimates equaling zero in all samples we substitute 10^-20 for their values in MGMR to avoid taking the log of zero. For P updates (e.g., equation 14), we avoid potentially negative P values by adding one to each alpha (thus ignoring -1 in the numerator and denominator). We implemented the algorithm in MATLAB, where the inputs required are read-gene map files for each sample as in SEQEM [7], and an initial P estimate matrix. Alphas are initialized as an M-length vector of ones.

Figure 1

A mesh representation of F(α) [equation (15)] showing non-convex behavior. P is a 10 × 2 constant matrix and α is varied on [0:50,0:50]. The case shown is for N = 10, M = 2 (ten samples, two α parameters). Non-convex behavior is demonstrated by the values on the plane defined by α₁ = .06 on the range [0,50] on the right.

Results

Experimental setup

As in [7, 9, 10], we examined MGMR's accuracy by comparing its estimates of known expression levels with those of existing methods, namely SEQEM [7] and RSEM [9, 10]. The initial "known" expression levels were estimates obtained from RNA-Seq samples; how these estimates were obtained is described below. In our case, we had to simulate a population sharing similar expression levels and the same set of gene regions. Our experiment differed in that we sought to use additional information to improve on the estimates of these existing methods. These methods were designed to estimate expression of single samples, and each had specific advantages which we disregarded in our comparison. For example, we ignored both SEQEM's ability to utilize SNP information and RSEM's ability to allow estimation on assembled transcripts by using only reference sequences.

Simulating data

To derive artificial reads, we first estimated expression on real biological samples using one method and then used the resulting distribution of expression values to generate simulated datasets for testing. Real samples were drawn from lymphoblastoid cells of the Yoruba in Ibadan (YRI) population [2, 3]. As MGMR requires initial expression estimates, the estimate derived from the method it was being compared with in each case was input to MGMR. Thus, the four initial estimates used were from SEQEM, MGMR(SEQEM) (namely, MGMR initialized by SEQEM's result), RSEM, and MGMR(RSEM) (namely, MGMR initialized by RSEM estimates). We denote these four estimates SEQEM-A, SEQEM-B, MGMR-A and MGMR-B, respectively. We simulated reads by randomly selecting start sites falling within gene boundaries and extracting sequences from those positions. Read lengths were defined for each experiment, and simulations were repeated multiple times to account for randomness in the sampling process.

To derive the sequence set for the SEQEM comparison, we expanded upon the procedure used in [7]. There, SEQEM was shown to improve estimation of paralogous gene expression on a set of exon sequences from 51 Homo Sapiens chromosome 1 paralogs from the HomoloGene [13] database. We extracted a larger set of sequences containing all HomoloGene paralogs in Homo Sapiens having at least one exon longer than twice the read length used that do not overlap in genomic coordinates. We required this minimal length because sequences were sampled from exons, and we needed to ensure enough positions existed for full length reads to be sampled from these exons. 285 such genes remained (for reads of length 35 bp), and these were the genes on which expression was tested and from which read sequences were derived. The SEQEM-A and SEQEM-B read sets were generated based on randomly selected exons from these genes and the expression levels from the SEQEM-A and SEQEM-B estimates taken on 20 YRI individuals. The read length of 35 bp corresponded to that of the YRI samples, and a coverage level of 20 was chosen, as this was the level at which SEQEM was shown to perform best in [7]. We performed a total of 30 repetitions of read simulations, where each repetition consisted of 20 samples (corresponding to the original 20 YRI samples used).

For the RSEM-A and RSEM-B read sets, the transcript set used was also obtained by filtering the HomoloGene database to avoid gene overlaps, but no length filtering was required: reads were now sampled directly from transcripts which all had effective lengths greater than the read length used. 524 transcripts corresponding to 265 genes survived this filtering. For these read sets, we produced 30 repetitions of 74 samples, where each consisted of 100 bp reads at a coverage level of 20. In all other respects the sampling process and read generation steps were identical to those performed for the SEQEM-A and SEQEM-B read sets.

Error measures

Accuracy was assessed by three error measures, the first two of which were applied in [7]: error rate, computed as \frac{1}{n}{\sum }_i\frac{\left|P_i-Q_i\right|}{Q_i},{\chi }^2 difference, computed as {\sum }_i\frac{{\left(P_i-Q_i\right)}^2}{Q_i}, and Kullback-Leibler divergence, computed as {\sum }_i{P}_i\text{log}\frac{P_i}{Q_i}. Here P is the estimated distribution generated by the tested algorithm and Q is the true distribution. Error measures were averaged over all repetitions per sample, and then over all samples.

Simulated data with priors based on real estimates - estimating paralogous gene expression

To test performance on paralogous gene estimates, we set out to compare independent sample SEQEM estimates with MGMR's common population estimates. Before applying SEQEM, we looked to address one criticism of it from [11], where it was suggested that SEQEM's estimation could be improved by incorporation of transcript length correction. Upon examination, we found the effect of this correction was an increase in accuracy, and thus we maintained it for subsequent tests. This improvement is documented in the appendix.

With this correction in place, we estimated expression levels on the SEQEM-A and SEQEM-B read sets, applying SEQEM and MGMR to each. Outputs were recorded at 1-10, 20, 30, 40, 50 and 100 iterations for MGMR and at 100 iterations for SEQEM. The results are shown in Figure 2. We observed that both error and variance levels dropped sharply within just a few iterations for MGMR, and converged to significantly better estimates on average than SEQEM. These trends were consistent across all three error measures [Table 1]. Variance seemed to diminish more with MGMR over time, as might be expected for a method that shares information across samples. Notably, MGMR outperformed SEQEM on estimates for samples based on SEQEM-A, where sample estimates were obtained independently (and thus we expect the variation inherent in the real samples to be maintained).

Figure 2

Relative error measured on SEQEM-A and SEQEM-B data sets. MGMR outputs on SEQEM-A and SEQEM-B initializations were compared with SEQEM up to 100 iterations. MGMR outputs were recorded at 1-10, 20, 30, 40, 50 and 100 iterations. The first few iterations have been trimmed to allow a compact presentation.

Table 1 MGMR vs. SEQEM error at 100 iterations on SEQEM-A and SEQEM-B data sets

Simulated data with priors based on real samples - estimating transcript level expression

We also sought to examine whether MGMR can improve results in the more challenging setting of estimating transcript level expression. Here, we expect estimates to be noisier due to low expression values in the real samples, and we must contend with multiread mappings due to paralogous genes as well as to isoforms of particular genes sharing subsequences as a result of alternative splicing. In anticipation of this challenge, we used a larger set consisting of 74 sample of single-end YRI samples as the real data source and simulated 100 bp reads instead of 35 bp. This was expected to be a difficult case for estimation, as all genes in the set are paralogs and many have multiple isoforms, as described in the section "Simulating data."

Once expression estimation was performed on the YRI samples and read sets RSEM-A and RSEM-B were generated, we again performed expression estimates with RSEM and MGMR on each set. In this case, unfortunately, we found the results did not exhibit a consistent trend as before and overall appeared inconclusive. These results are summarized in Table 2. It remains to be seen why the error results differ according to the level of estimation (gene vs. transcript) performed.

Table 2 MGMR vs. RSEM error at 100 iterations on RSEM-A and RSEM-B data sets

Conclusion

As shown by the 1000 Genomes and HapMap projects, one of the drives of modern genetics and bioinformatics research is to characterize variation in populations. Because of cost and time constraints, such projects have only recently become feasible. In addition to such studies assessing genomic variation and its relation to disease phenotypes based on DNA, it is anticipated that RNA-Seq population studies will also grow in popularity to more directly assign functional significance to variant loci by means of transcription measures. Thus, it becomes essential to accurately measure the expression levels from each individual to characterize such variation. Here, we have shown that for one common study design an unexpected benefit can arise. When individuals in these studies are drawn from the same population, the estimates made on each can be made more accurate because of the commonalities among population members.

A shortcoming of the MGMR approach is that since it assumes commonality among the samples, outlier samples will be attracted towards the common denominator, and thus appear more similar to the group profile than they really are. In particular, if the data are subject to differential expression analysis, MGMR may reduce the number of differentially expressed genes.

We have investigated the efficacy of MGMR in tackling two typical experimental settings - inferring expression levels of paralogs at the gene level, and of isoforms (also drawn from a difficult set of paralogs). Although substantial gains were obtained in the first, more inquiry is required to demonstrate a benefit in the latter. It is worth noting that in each case at least a quarter of the regions considered showed improvement, as shown in Table 3. With these results, we submit a proof of concept that population structure can aid in estimation of expression levels for RNA-Seq samples.

Table 3 Proportion of genes for which MGMR improves estimates on different data sets