Volume 13 Supplement 6
Proceedings of the Second Annual RECOMB Satellite Workshop on Massively Parallel Sequencing (RECOMBseq 2012)
MGMR: leveraging RNASeq population data to optimize expression estimation
 Roye Rozov^{1},
 Eran Halperin^{1, 2, 3}Email author and
 Ron Shamir^{1}
DOI: 10.1186/1471210513S6S2
© Rozov et al.; licensee BioMed Central Ltd. 2012
Published: 19 April 2012
Abstract
Background
RNASeq 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 RNASeq 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 likelihoodbased approaches. In many instances, RNAseq 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 RNAseq 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, RNASeq has emerged as a legitimate competitor to microarrays 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 finegrained assessment of variation in expression and the sources for such variation, as exemplified by recent largescale RNASeq 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 RNASeq is not immune to bias [4–6]. Perhaps the most widely discussed hurdle to accurate estimation in the case of RNASeq 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
RNASeq multiread expression estimation
As we seek to extend the prevalent generative model of RNASeq [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 ${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 perbase 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.
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_{11}, ...,P_{1M}), .., (P_{N 1}...,P_{ N M })] be probability densities drawn from a common Dirichlet distribution, defined by a set of hyperparameters specific to the population: ∀iϵ[1, N], p_{ i }= (P_{i 1},..., P_{ iM }) ~ Dirichlet (α_{1},..., α_{ 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 readregion 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.
MultiGenome MultiRead (MGMR) algorithm
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 viceversa until the total change in α becomes sufficiently small (or until a preset number of iterations have been executed).
Although for EMbased 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 EMbased P update followed by one α update.
Estimating P given α
If we assume α is given, we can write the EM steps to find p_{ 1 },...,p_{ N }:
Estimating α given P
where $log{\stackrel{\u0304}{P}}_{k}=\frac{1}{N}*{\sum}_{i}log{P}_{ik}$.
Heuristics/Implementation
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 RNASeq 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 SEQEMA, SEQEMB, MGMRA and MGMRB, 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 SEQEMA and SEQEMB read sets were generated based on randomly selected exons from these genes and the expression levels from the SEQEMA and SEQEMB 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 RSEMA and RSEMB 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 SEQEMA and SEQEMB 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 KullbackLeibler 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.
MGMR vs. SEQEM error at 100 iterations on SEQEMA and SEQEMB data sets
SEQEMA sampling  SEQEMB sampling  

SEQEM  MGMR  SEQEM  MGMR  
Error  SD  Error  SD  Error  SD  Error  SD  
E  1.27  1 * 10^{2}  1.03  0.14  1.50  0.70  0.82  6 * 10^{3} 
χ^{2}  0.66  2 * 10^{3}  0.22  4 * 10^{3}  0.69  0.05  0.27  1 * 10^{4} 
KL  0.29  7 * 10^{4}  0.14  1 * 10^{4}  0.18  2 * 10^{4}  0.17  1 * 10^{4} 
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 singleend 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 RSEMA and RSEMB data sets
RSEMA Sampling  RSEMB Sampling  

RSEM  MGMR  RSEM  MGMR  
Error  SD  Error  SD  Error  SD  Error  SD  
E  0.1  1 * 10^{3}  0.69  1 * 10^{3}  1.0  1 * 10^{4}  0.61  1 * 10^{3} 
χ^{2}  0.02  6 * 10^{4}  1.25  0.01  0.02  9 * 10^{4}  0.58  3 * 10^{4} 
KL  1.5  0.22  0.6  1 * 10^{3}  0.8  0.11  0.38  6 * 10^{4} 
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 RNASeq 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
SEQEMA  SEQEMB  RSEMA  RSEMB  

Proportion  104/285  78/285  126/524  173/524 
%  36.5  27.3  24.0  33.0 
List of abbreviations
 E:

relative error rate
 χ ^{2} :

Chisquared error
 KL:

KullbackLiebler divergence
 SD:

standard deviation
 bp:

base pair
Declarations
Acknowledgements
E.H. is a faculty fellow of the Edmond J. Safra Bioinformatics program at TelAviv University. R.S. was supported in part by the European Community's Seventh Framework Programme (grant HEALTHF42009223575 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 (RECOMBseq 2012).
Authors’ Affiliations
References
 Montgomery SB, Sammeth M, GutierrezArcelus M, Lach RP, Ingle C, Nisbett J, Guigo R, Dermitzakis ET: Transcriptome genetics using second generation sequencing in a Caucasian population. Nature. 2010, 464: 773777. 10.1038/nature08903.View ArticlePubMed
 Pickrell JK, Marioni JC, Pai AA, Degner JF, Engelhardt BE, Nkadori E, Veyrieras JB, Stephens M, Gilad Y, Pritchard JK: Understanding mechanisms underlying human gene expression variation with RNA sequencing. Nature. 2010, 464: 768772. 10.1038/nature08872.PubMed CentralView ArticlePubMed
 Frazer KA et al.: A second generation human haplotype map of over 3.1 million SNPs. Nature. 2007, 449: 851861. 10.1038/nature06258.View ArticlePubMed
 Robinson MD, Oshlack A: A scaling normalization method for differential expression analysis of RNAseq data. Genome Biol. 2010, 11: R2510.1186/gb2010113r25.PubMed CentralView ArticlePubMed
 Bullard JH, Purdom E, Hansen KD, Dudoit S: Evaluation of statistical methods for normalization and differential expression in mRNASeq experiments. BMC Bioinformatics. 2010, 11: 9410.1186/147121051194.PubMed CentralView ArticlePubMed
 Oshlack A, Wakefield MJ: Transcript length bias in RNAseq data confounds systems biology. Biol Direct. 2009, 4: 1410.1186/17456150414.PubMed CentralView ArticlePubMed
 Pasaniuc B, Zaitlen N, Halperin E: Accurate estimation of expression levels of homologous genes in RNAseq experiments. J Comput Biol. 2011, 18: 459468. 10.1089/cmb.2010.0259.View ArticlePubMed
 Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, Salzberg SL, Wold BJ, Pachter L: Transcript assembly and quantification by RNASeq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010, 28: 511515. 10.1038/nbt.1621.PubMed CentralView ArticlePubMed
 Li B, Ruotti V, Stewart RM, Thomson JA, Dewey CN: RNASeq gene expression estimation with read mapping uncertainty. Bioinformatics. 2010, 26: 493500. 10.1093/bioinformatics/btp692.PubMed CentralView ArticlePubMed
 Li B, Dewey CN: RSEM: accurate transcript quantification from RNASeq data with or without a reference genome. BMC Bioinformatics. 2011, 12: 32310.1186/1471210512323.PubMed CentralView ArticlePubMed
 Pachter L: Models for transcript quantification from RNASeq. ArXiv eprints. 2011
 Minka TP: Estimating a Dirichlet distribution. 2003, [http://research.microsoft.com/\~minka]
 [http://www.ncbi.nlm.nih.gov/homologene/]
Copyright
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.