MGMR: leveraging RNA-Seq population data to optimize expression estimation
© Rozov et al; licensee BioMed Central Ltd. 2012
Published: 19 April 2012
Skip to main content
© Rozov et al; licensee BioMed Central Ltd. 2012
Published: 19 April 2012
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
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.
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.
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  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.
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 1, ...,G M ) be the set of M transcribed regions considered and P = (P 1, ...,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 1,..., 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 . 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 , ϵ 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 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.
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 11, ...,P 1M ), .., (P N1 ...,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 i1,..., P iM ) ~ Dirichlet (α 1,..., α M ).
For sample i, we denote the set of reads as , 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.
We wish to estimate α and p 1 ,...,p n to maximize equation (7) above. For this purpose, we adopt an alternating iterative procedure of estimating α given the current estimate of p 1 ,...,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.
If we assume α is given, we can write the EM steps to find p 1 ,...,p N :
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  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.
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 . 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  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 . 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.
Accuracy was assessed by three error measures, the first two of which were applied in : error rate, computed as difference, computed as , and Kullback-Leibler divergence, computed as . 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.
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 , 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.
MGMR vs. SEQEM error at 100 iterations on SEQEM-A and SEQEM-B data sets
1 * 10-2
6 * 10-3
2 * 10-3
4 * 10-3
1 * 10-4
7 * 10-4
1 * 10-4
2 * 10-4
1 * 10-4
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."
MGMR vs. RSEM error at 100 iterations on RSEM-A and RSEM-B data sets
1 * 10-3
1 * 10-3
1 * 10-4
1 * 10-3
6 * 10-4
9 * 10-4
3 * 10-4
1 * 10-3
6 * 10-4
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.
Proportion of genes for which MGMR improves estimates on different data sets
relative error rate
E.H. is a faculty fellow of the Edmond J. Safra Bioinformatics program at Tel-Aviv University. R.S. was supported in part by the European Community's Seventh Framework Programme (grant HEALTH-F4-2009-223575 for the TRIREME project) and by the Israel Science Foundation (grant 802/08).
This article has been published as part of BMC Bioinformatics Volume 13 Supplement 6, 2012: Proceedings of the Second Annual RECOMB Satellite Workshop on Massively Parallel Sequencing (RECOMB-seq 2012).