 Method
 Open access
 Published:
Haplotype and isoform specific expression estimation using multimapping RNAseq reads
Genome Biology volume 12, Article number: R13 (2011)
Abstract
We present a novel pipeline and methodology for simultaneously estimating isoform expression and allelic imbalance in diploid organisms using RNAseq data. We achieve this by modeling the expression of haplotypespecific isoforms. If unknown, the two parental isoform sequences can be individually reconstructed. A new statistical method, MMSEQ, deconvolves the mapping of reads to multiple transcripts (isoforms or haplotypespecific isoforms). Our software can take into account nonuniform read generation and works with pairedend reads.
Background
Highthroughput sequencing of RNA, known as RNAseq, is a promising new approach to transcriptome profiling. RNAseq has a greater dynamic range than microarrays, which suffer from nonspecific hybridization and saturation biases. Transcriptional subsequences spanning multiple exons can be directly observed, allowing more precise estimation of the expression levels of splice variants. Moreover, unlike traditional expression arrays, RNAseq produces sequence information that can be used for genotyping and phasing of haplotypes, thus permitting inferences to be made about the expression of each of the two parental haplotypes of a transcript in a diploid organism.
The first step in RNAseq experiments is the preparation of cDNA libraries, whereby RNA is isolated, fragmented and synthesized to cDNA. Sequencing of one or both ends of the fragments then takes place to produce millions of short reads and an associated base call uncertainty measure for each position in each read. The reads are then aligned, usually allowing for sequencing errors and polymorphisms, to a set of reference chromosomes or transcripts. The alignments of the reads are the fundamental data used to study biological phenomena such as isoform expression levels and allelic imbalance. Methods have recently been developed to estimate these two quantities separately but no approaches exist to make inferences about them simultaneously to estimate expression at the haplotype and isoform ('haploisoform') level. In diploid organisms, this level of analysis can contribute to our understanding of cis vs. trans regulation [1] and epigenetic effects such as genomic imprinting [2]. We first set out the problems of isoform level expression, allelic mapping biases and allelic imbalance, and then propose a pipeline and statistical model to deal with them.
Isoform level expression
Multiple isoforms of the same gene and multiple genes within paralogos gene families often exhibit exonic sequence similarity or identity. Therefore, given the short length of reads relative to isoforms, many reads map to multiple transcripts (Table 1). Discarding multimapping reads leads to a significant loss of information as well as a systematic underestimation of expression estimates. For reads that map to multiple locations, one solution is to distribute the multimapping reads according to the coverage ratios at each location using only singlemapping reads [3]. However, this does not address the problem of inferring expression levels at the isoform level.
Essentially, the estimation of isoform level expression can be done by constructing a matrix of indicator functions M_{ it } = 1 if region i belongs to transcript t. The 'regions' may for now be thought of as exons or part exons, though we later define them more generally. Using this construction it is natural to define a model:
where X_{ it } are the (unobserved) counts of reads from region i of transcript t, b is a normalization constant used when comparing experiments, μ_{ t } is a parameter representing the expression of transcript t and s_{ i } is the effective length of region i (that is the number of possible start positions for reads in the region). This model can be fit using an expectation maximization (EM) algorithm, since the X_{ it } are unobserved but their sums across transcripts {k}_{i}\equiv {\displaystyle {\sum}_{t}{X}_{it}} are observed.
This model has been used by [4] in their POEM software, with i representing exons. Their method does not use reads that span multiple exons or reads that map to multiple genes. The same model has been used in [5], with i representing exons or part exons, or regions spanning exon junctions, enabling good estimation of isoform expression within genes. They do not, however, include reads mapping to multiple genes. The RSEM method [6] employs a similar model, but models the probability of each read individually, rather than read counts. This method allows reads to come from multiple genes as well as multiple isoforms of the same gene. The modeling of individual reads allows RSEM to accommodate general positionspecific biases in the generation of reads. However, two recent papers [7, 8] have shown that deviations from uniformity in the generation of reads are in great part sequence rather than positiondependent for a given experimental protocol and sequencing platform. Furthermore, the computational requirements of modeling individual reads increasing proportionately with read depth, which, in the case of RSEM, is exacerbated further by the use of computationally intensive bootstrapping procedures to estimate standard errors. None of the above methods are compatible with pairedend data. A recently published method, Cufflinks [9], focuses on transcript assembly as well as expression estimation using an extension of the [5] model that is compatible with pairedend data. However, this method does not model sequencespecific uniformity biases and uses a fixed downweighting scheme to account for reads mapping to more than one transcription locus, meaning that the abundances of transcripts in different regions are estimated independently.
Allelic imbalance
Studies of imbalances between the expression of two parental haplotypes have mostly been restricted to testing the null hypothesis of equal expression between two alleles at a single heterozygous base, typically with a binomial test [1, 2, 10]. However, as transcripts may contain multiple heterozygotes, a more powerful approach is to assess the presence of a consistent imbalance across all the heterozygotes in a gene together. This has been done on a casebycase basis using read pairs that overlap two heterozygous SNPs [11] while [12] propose an extension to the binomial test for detecting allelic imbalance that takes into account all SNPs and their positions in a gene. However, this approach, which is a statistical test rather than a method of quantifying haplotypespecific expression, assumes imbalances to be homogeneous along genes and thus does not take into account the possibility of asymmetric imbalances between isoforms of the same gene.
Allelic mapping biases
Aligners usually have a maximum tolerance threshold for mismatches between reads and the reference. Reads containing nonreference alleles are less likely to align than reads matching the reference exactly, so genes with a high frequency of nonreference alleles may be underestimated. Ideally, aligners would accept ambiguity codes for alleles that segregate in the species (cf. Novoalign [13]), but no free software is currently able to do this. A possible workaround is to change the nucleotide at each SNP to an allele that does not segregate in the species, as has been proposed to remove biases when estimating allelic imbalance [10]. However, in the context of gene expression analysis, this leads to even greater underestimation of genes with many nonreference alleles and an increase in incorrect alignments to homologous regions. Instead, we propose aligning to a samplespecific transcriptome reference, constructed from (potentially phased) genotype calls.
MMSEQ
In this paper we present a new pipeline, including a novel statistical method called MMSEQ, for estimating haplotype, isoform and gene specific expression. The MMSEQ software is straightforward to use, fully documented and freely available online [14] and as part of ArrayExpressHTS [15]. Our pipeline exploits all reads that can be mapped to at least one annotated transcript sequence and reduces the number of alignments missed due to the presence of nonreference alleles. It is compatible with pairedend data and makes use of inferred insert size information to choose the best alignments. Our method permits estimating the expression of the two versions of each heterozygotecontaining isoform ('haploisoform') individually and thus it can detect asymmetric imbalances between isoforms of the same gene. Our software further takes into account sequencespecific deviations from uniform sampling of reads using the model described in [8] but can flexibly accommodate other models. We validate our method at the isoform level with a simulation study, comparing our results to RSEM's, and applying it to a published Illumina dataset consisting of lymphoblastoid cell lines from 61 HapMap individuals [16]. We validate our method at the haploisoform level by showing we can deconvolve the expression estimates of haploisoforms on the nonpseudoautosomal (nonPAR) region of the X chromosome using a pooled dataset of two HapMap males. We further apply our method to a published dataset of F_{1} initial and reciprocal crosses of CAST/EiJ (CAST) and C57BL/6J (C57) inbred mice [2] and demonstrate that MMSEQ is able to detect parental imbalance between the two haplotypes of each isoform.
Results
Overview of the pipeline
The pipeline can be depicted as a flow chart with two different start positions (Figure 1):

(a)
Expression estimation using alignments to a predefined transcriptome reference,

(b)
Expression estimation using alignments to a transcriptome reference that is obtained from the RNAseq data.
In case (a), the level of estimation (haploisoform or isoform) depends on whether the reference includes two copies of heterozygous transcripts. In case (b), it depends on whether the genotypes are phased. The most exhaustive use of the pipeline proceeds as follows. First, the reads are aligned to the standard genome reference using TopHat [17]. Then, genotypes are called with SAMtools pileup [18]. Genotypes are then phased with polyHap [19] using population genotype data to produce a pair of haplotypes for all gene regions on the genome. The standard transcriptome reference is then edited for each individual to match the inferred haplotypes. The reads are realigned to the individualized haplotype specific transcriptome reference with Bowtie [20], finding alignments for reads that originally failed to align due to having too many mismatches with the standard reference (approximately 0:3% more reads recovered, with some transcripts receiving up to 13% more hits, in the HapMap dataset [16]). Finally, our new method, MMSEQ, is used to disaggregate the expression level of each haploisoform.
MMSEQ
Poisson model
We use the model in Equation 1 as a starting point for modeling gene isoforms and extend it to apply to haploisoforms. First, we employ a more general definition of 'region': each read maps to one set of transcripts, which may belong to the same gene or to various different genes, and which can have two versions, one containing the paternal and the other the maternal haplotype. These sets are labeled by i. Many reads will map to the exact same set, hence we can model reads counts (k_{ i } ) for the set. The M_{ it } are defined very straightforwardly as the indicator functions for transcript t belonging to set i. The region length s_{ i } is the effective length of the sequence shared between the whole set. If the set of transcripts all belong to the same gene and haplotype, then s_{ i } may be the effective length of an exon or part exon. However, aligned reads often map to multiple genes equally well (Table 1) so the region need not correspond to an actual region on the genome. Using our definition of a region, the s_{ i } would be difficult to calculate given the sheer number of overlaps and regions, but in fact the s_{ i } are not needed in the calculation of the model (see Materials and Methods). Hence we have a model for read counts in which the data and fixed quantities (k_{ i } and M_{ it } ) are calculated in a straightforward way, and which allows for reads mapping to multiple isoforms of the same or different genes in exons or exon junctions and to paternal and maternal haplotypes separately.
Without loss of generality, Figure 2a illustrates our formulation for a gene with an alternatively spliced cassette exon and Figure 2b illustrates it for a gene with a single heterozygous base. The heterozygote casts a 'shadow' upstream of length equal to the read length, which acts like an alternative middle exon. This is because reads with starting positions within the shadow cover the heterozygote and contain one of the two alleles, thus mapping to only one of the two haplotypes.
We now formulate a Poisson model for read counts from transcript sets:
where b is a normalization constant, ∑ _{ t } M_{ it }μ_{ t } is the total expression from the transcript set i and s_{ i } is the effective length of the region of shared sequence between transcripts in set i. Figure 2a shows how the s_{ i } can be calculated for the gene with a cassette exon. Note that the sum of lengths of all the regions shared by transcript t add up to its effective length (transcript length minus read length plus one for uniformly generated reads): ∑ _{ i } s_{ i }M_{ it } = l_{ t } , so the transcriptset model is consistent with the usual Poisson model. Setting l_{ t } to the transcript length minus read length plus one is appropriate if a constant Poisson rate is assumed along all positions in a transcript: {r}_{t}\sim Pois(b{\displaystyle {\sum}_{p=1}^{{l}_{t}}{\mu}_{t}})\sim Pois(b{l}_{t}{\mu}_{t}), where r_{ t } is the number of reads originating from transcript t and the sum is over all possible starting read positions p. The nonuniformity of read generation demonstrated in [8], however, suggests a variablerate Poisson model:
where {\tilde{l}}_{t} is an adjusted effective length, referred to as the sum of sequence preferences (SSP) in [8]. We use their Poisson regression model to adjust the length of each transcript based on its sequence, but other adjustment procedures may be used instead. Briefly, the logarithm of the sequencing preference of each possible start position in a transcript is calculated as the sum of an intercept term plus a set of coefficients determined by the sequence immediately upstream and downstream of the start position. It would also be possible to integrate the method described in [7], which uses a weighting for reads based on the first seven nucleotides of their sequences, by applying this weighting in our calculation of k_{ i } . However, this approach does not incorporate the effects of the sequence composition on the reference upstream of the read start positions or further downstream than seven bases, and we thus prefer to use the [8] method instead. The normalization constant b is used to make lanes with different read depths comparable. We set b to the total number of reads (in millions) and measure transcript lengths in kilobases, which means the scale of the expression parameter μ_{ t } is equivalent to RPKM (reads per kilobase per million mapped reads) described in [3]. In downstream analysis, a more robust measure can be used, such as the library size parameter suggested by [21].
The only unknown parameters in the model are the μ_{ t } . The observed data are the k_{ i } and the matrix M and effective transcript lengths l_{ t } are known. In principle the effective lengths of the transcript sets s_{ i } can be calculated, but in fact, they are not needed (see Materials and Methods).
Inference
The maximum likelihood (ML) estimate of μ_{ t } cannot be obtained analytically, so instead we use an expectation maximization (EM) algorithm to compute it, an approach also taken by [4, 6] for isoforms. After convergence of the algorithm, we output the estimates of μ_{ t } and refer to them as MMSEQ EM estimates.
The usual approach to estimating statistical standard errors of ML estimators requires inversion of the observed information matrix. When analyzing the expression of thousands of transcripts, the high dimensionality of the observed information matrix and the possibility of identical columns due to gene homology make this approach impracticable. Bootstrapping may also be used to estimate errors, as in [6], but it is a computationally intensive method requiring repeated runs of the EM algorithm. Instead we use a simple Bayesian model with a vague prior on μ_{ t } . As before, we use the augmented data reads per region and transcript, X_{ it } . The full model is:
Again, the only lengths needed are the l_{ t } . The conjugacy of the PoissonGamma model makes the sampling fast and straightforward as the full conditionals are in closed form (see Materials and Methods). We use the final EM estimate of the μ_{ t } as the initial values for the Gibbs sampling. We then produce samples from the whole posterior distributions of the μ_{ t } and calculate the sample means and their respective Monte Carlo standard errors (MCSE), which take into account the autocorrelations of the samples [22]. We set the hyperprior parameters to α = 1.2 and β = 0.001, producing a vague prior on the μ_{ t } that captures the wellknown broad and skewed distribution of gene expression values. We output the means of the Gibbs samples of μ_{ t } , which we refer to as MMSEQ GS estimates. As we shall show, the regularization afforded by the Bayesian algorithm produces estimates with a lower error than the MMSEQ EM estimates. Moreover, it can readily be shown that for transcript with low coverage, the ML estimate is often zero, even though this is likely to be an underestimate of the expression. For example, suppose there exist two equallyexpressed haploisoforms differing by only one heterozygote. Under the assumption of uniform sampling of 0.01 reads per nucleotide for both haploisoforms, if the read length is 35, then the probability of observing a read containing one allele but no reads containing the other allele is fairly high (2(1e^{0.35})e^{0.35} ≃ 0.42). The ML estimate of the haploisoform with the unsampled allele under this scenario is zero while the ML estimate of the haploisoform with the sampled allele is overestimated. With Gibbs sampling, on the other hand, this effect is tempered by the Gamma prior. The MMSEQ GS estimates are thus our preferred expression measures.
Best mismatch stratum filter
While a read may align to multiple transcripts, not all alignments may be equally reliable. We therefore filter out all alignments that do not have the minimal number of mismatches for a given read or read pair (similar to the strata switch in Bowtie, but compatible with paired as well as single end data). In the case of pairedend data, the number of mismatches from both ends is added up to determine the 'mismatch stratum' of a read pair. This filter is crucial in order to correctly discriminate between the two versions of an isoform at a heterozygous position, since reads from one haplotype also match the alternative haplotype with an additional mismatch. The stratum filter thus ensures that reads are properly assigned to the correct haplotype.
Insert size filter for pairedend data
For pairedend data, both reads in a pair must align to a transcript for the mapping to be considered. If the fragments are sufficiently large, the alignments may span three exons and align to transcripts that both retain and skip the middle exon. However, the alignment with an inferred fragment size (also called insert size) that is nearer to the expected insert size from the fragmentation protocol, is more likely to be correct. We exploit this information by applying an insert size filter to alignments in the best mismatch stratum for each read. If an alignment's insert size is nearer than x bp (for example equivalent to one standard deviation) away from the expected insert size, then all other alignments for that read with an insert size greater than x bp away from the expected insert size are removed. This filter can be thought of as an extension of mismatchbased filtering for reporting only alignments with moderately high probability of being true. Although full probabilistic modeling is more principled, filtering is a commonplace approach to reducing alignment candidates for each read to a set that can be dealt with pragmatically. For the HapMap dataset, mistakes in the protocol resulted in two distributions of insert sizes within some samples, so we omitted this filter.
MMSEQ output
The mmseq program produces three files each containing EM and GS expression estimates with associated MCSEs. The first file provides estimates at the transcript/haploisoform level, the second file provides aggregate estimates for sets of transcripts that have been amalgamated due to having identical sequences (and therefore indistinguishable expression levels), and the third file aggregates transcript estimates into genes, thus providing genelevel estimates. Homozygous transcripts are aggregated together, whereas heterozygous transcripts are aggregated separately to produce 'haplogene' level estimates. With respect to transcripts that have identical sequences and hence indistinguishable and unidentifiable expression levels, the posterior samples exhibit high variance and strong anticorrelation but the sum of their expression can be precisely estimated (Additional file 1). We therefore recommend use of the amalgamated estimates.
Performance and scalability
The performance of the EM and Gibbs algorithms is determined principally by the size of the M matrix, which is bounded by the total number of known transcripts and the total number of combinations of transcripts that share sequence. Marginal increases in the total number of observed reads do not result in commensurate increases in the size of M, because additional reads tend to map to transcript sets that have been mapped to by previous reads (Table 2). Consequently, the mmseq program exhibits economies of scale which allow it to cope with future increases in throughput. This contrasts with the RSEM method, which represents each read separately in their indicator matrix that maps reads to isoforms [6].
Correction for nonuniform read sampling
We have assessed the effect of applying the Poisson regression [8] correction for nonuniform sampling using read data from three Illumina Genome Analyzer II (GAII) lanes from the HapMap dataset [16] (described below). Two of the samples were from the same run (ID 3125) and a third from a separate run (ID 3122). We obtained Poisson regression coefficients for 20 bases upstream and downstream of each possible start position using the first 10 million alignments for each lane. The regression model was fitted using only the most highly expressed transcripts, as these have the best signaltonoise ratio [8]. Specifically, from the 500 transcripts with the highest average number of nucleotides per position, we selected a subset containing only one transcript per gene so as to avoid doublecounting of sequence preferences. As shown in Additional file 2, the coefficients are highly stable across both lanes and runs. The timeconsuming task of calculating adjusted transcript lengths separately for each lane is therefore unnecessary. Instead, our software can reuse the adjusted transcript lengths calculated from one sample when analyzing other samples. Variations in the Poisson rate from base to base tend to average out over the length of each transcript, and thus the adjustments to the lengths are generally slight (Additional file 3). As expected from the Poisson model (Equation 3), changes in the expression estimates (estimates of μ_{ t } ) tend to be inversely proportional to adjustments to the lengths. Nevertheless, as transcripts sharing reads may be adjusted in opposite directions, for some transcripts even a small change in the length has a significant impact on the expression estimate (Figure 3).
Simulation study of isoform expression estimation
We simulated reads from human and mouse Ensembl cDNA files under the assumption of uniform sampling of reads and ran the MMSEQ workflow. We found good correlation between simulated and estimated expression values and between dispersion around the true values and estimated MCSEs. We did however observe a small upward bias in our estimates of transcripts with low expression levels, attributable to our use of the mean to summarize highly skewed distributions. We evaluated our genelevel estimates by summing over the isoform components within each gene. As anticipated, we obtained more precise estimates for genes than for transcripts (Figure 4).
We also observed better estimates for mouse, which has 45,452 annotated transcripts, than for human, which has higher splicing complexity manifested in 122,636 annotated transcripts (Figure 5). Transcripts may be connected to other transcripts via reads that align to regions shared by isoforms of the same gene or to different genes with sequence homology. The complexity of the graph that connects transcripts with each other reflects the ambiguity in the assignment of reads to transcripts and thus the errors in our estimates. A bar plot of the number of transcripts that each transcript is connected to in human and mouse demonstrates a significant difference in complexity between the annotated transcriptomes of the two species (Additional file 4).
Comparison of isoform expression estimation between MMSEQ and RSEM
Like MMSEQ, the RSEM method [6] makes use of all classes of reads to estimate isoform expression. The authors have shown an improvement of their method for genelevel estimation over strategies that discard multiply aligned reads or allocate them to mapped transcripts according to the coverage by singlemapping reads (as in [3]). However, isoformlevel results for their method have not been assessed. We obtained RSEM estimates for Ensembl transcripts using our simulated human sequence dataset for the purposes of comparison.
We scaled our simulated and estimated expression values to add up to one in order to make them comparable to RSEM's fractional expression estimates. We found that RSEM and MMSEQ EM are comparable but, unlike the MMSEQ EM algorithm, RSEM tended to overestimate some mediumexpression transcripts. Both the RSEM and MMSEQ EM algorithms tended to underestimate some lowexpression transcripts, pushing them very close to zero and thus producing very large errors on the log scale. This was avoided by the regularization of the Gibbs algorithm, which produced tighter estimates and only overestimated slightly some very lowly expressed transcripts (Figure 5 and Additional file 5), showing the benefits of using the whole posterior distribution of μ_{ t } to estimate expression rather than a maximization strategy.
Isoformlevel application to the HapMap dataset
The HapMap pairedend Illumina GAII dataset [16] consists of 73 lanes: 7 lanes for the same Yoruban individual, another 7 lanes for the same CEU individual and the remaining 59 lanes each for different CEU individuals. The authors assessed exoncount correlations between the lanes. Here we look at transcript and genelevel correlations. We analyzed the data using the MMSEQ pipeline, aligning approximately 75% of reads to Ensembl human reference transcripts. The average rank correlation was 0.92 and 0.84 respectively at the gene and transcript level (Figure 6). When comparing identical samples at the gene level the rank correlation ranged from 0.96 to 0.97 for the Yoruban individual and from 0.92 to 0.97 for the CEU individual. At the transcript level, the ranges were 0.91 to 0.92 and 0.90 to 0.91 for the Yoruban and CEU individuals respectively. The transcriptlevel values are comparable to exoncount correlations found by [16]. Both are lower than the genelevel correlation, as might be expected due to the inclusion of withingene variance.
Although the ordering of transcripts and genes was broadly maintained even between lanes belonging to different individuals and runs, we found a striking contrast in the distribution of expression values between lanes of the same individual and lanes of different individuals (Additional file 6). The consistency of expression values for lanes of the same individual indicates that the technical replicability of the Illumina GAII sequencer is extremely high and therefore that the variation observed between lanes from different individuals is mostly a reflection of biological variability. This is in line with previous research showing that sequence count data follow a negative binomial distribution in biological replicates and a Poisson distribution in technical replicates [21]. As such, we expect the variance of our estimates to be proportional and greater than proportional to the expression values for technical and biological replicates respectively. This is indeed borne out both at the gene and transcript level (Additional file 7) and corroborates the need to take into account extra variability for highlyexpressed transcripts in differential expression analysis with biological replication (see Discussion).
Validation of haploisoform deconvolution
The nonpseudoautosomal region (nonPAR) of the X chromosome in human males is haploid, and thus the alleles in that region can be called directly without the need for phasing. We validated our method for deconvolving expression between two haplotypes of the same isoform as follows. We used the RNAseq data of two males from the HapMap data (NA12045 and NA12872) to call their haplotypes. We identified 117 isoforms on the nonPAR of the X chromosome that differed between the two individuals. We created custom transcriptome references for each of the two males, containing their individual versions of the 117 isoforms. We then created a third hybrid reference containing two copies of the 117 isoforms, one matching the haplotype of one male and the second matching the haplotype of the other. This hybrid reference mimics the case of a female with two X chromosomes with unknown expression of the two parental copies of each isoform. We obtained individual expression estimates of the 117 isoforms using the separate transcriptome references in each male and compared them with estimates obtained by aligning a dataset pooled from the data of both males to the hybrid reference. Although the original correlation between the two males was 0.85, the correlation between the individual estimates and the deconvolved estimates was 0.96 and 0.98, showing MMSEQ is capable of disaggregating the expression from paternal and maternal isoforms (Additional file 8).
To test whether MMSEQ is able to recover greater imbalances than found naturally between the two male individuals, we divided the genes of the 117 isoforms that are heterozygous in the hybrid reference into three equalsized groups. For one group, we artificially removed 90% of the reads hitting one male and, for another group, we artificially removed 90% of the reads hitting the other male. This reduction of reads mimics what would be observed if more extreme imbalances existed. We thus reduced the correlation between the log expression of the two males from 0.85 to 0.48. Despite this large imbalance, there was a correlation of 0.91 and 0.95 between the individual and the deconvolved estimates obtained from the pooled dataset (Figure 7), showing that MMSEQ is able to accurately disaggregate haplotypespecific expression in the presence of large imbalances.
Demonstration of haploisoform expression estimation using an F_{1}hybrid mouse brain dataset
We have applied MMSEQ to a published murine embryonic day 15 RNAseq dataset of CAST/C57 initial (F_{1}i) and reciprocal (F_{1}r) crosses [2]. Each RNA sample was a pool from four individuals. The C57 reference transcriptome used by the authors is available from the UCSC Genome Browser [23]. The authors called SNPs by aligning reads from the CAST samples to the C57 reference. We created a CAST reference transcriptome by changing alleles in the C57 reference sequences according to those SNP calls. The two references were combined in a hybrid reference containing two entries for isoforms that differed in sequence between C57 and CAST. Thus there is a onetoone mapping between SNPs called in the parents and heterozygotes in the hybrids. The data consist of 152 and 159 million 36 bp Illumina GAII reads for F_{1}i and F_{1}r respectively.
A scatterplot of the CAST/C57 differential expression between F_{1}i and F_{1}r crosses reveal a clear clustering of points into three groups (Figure 8). Firstly, the points on the upperleft to lowerright diagonal correspond to transcripts which show imbalance towards the parent of origin, suggesting they are imprinted. Those on the upperleft quadrant and bottomright quadrant correspond to maternally and paternally imprinted transcripts respectively. Transcripts termed 'consensus imprinted' by [2] are highlighted in color. These were defined arbitrarily by the authors as transcripts with more than two heterozygotes exhibiting imbalance in favor of the same parental sex, at least one of which was significant in a χ^{2} goodnessof fit test with a Pvalue threshold of 0.05.
We also identified a clustering of transcripts that exhibited CAST overexpression in the F_{1}i hybrids but approximately balanced expression in the F_{1}r's. We identified the cluster as consisting wholly of transcripts on the X chromosome (Additional file 9), which suggests that the initial crosses were male and the reciprocal crosses female. The sexes of the hybrid mice are in fact unknown. There was a slight skew in favor of the CAST strand in the reciprocal crosses. We think it is unlikely that this was due to mapping biases, since the CAST reference was produced from SNP calls against the C57 reference and was thus of lower quality, so any mapping bias would be expected to be in favor of C57. Moreover, [24] found a similar skew in adult samples of the same crosses. It is possible that the skew is the result of a selective bias in favor of C57 Xinactivated cells [25], possibly caused by one or more of the three mutations on the Xinactivation transcript (UCSC ID uc009tzp.1) or mutations in its promoter region.
The third grouping in the plot is on the lowerleft to upperright diagonal. These transcripts demonstrate consistent CAST/C57 differential expression regardless of the sexstrain combination of the parents, and are thus indicative of cis regulation.
One advantage of MMSEQ is that imbalances are assessed at the transcript level rather than for individual SNPs. Thus it is not necessary to set arbitrary thresholds on the numbers of heterozygotes or the magnitude and significance of the imbalances to make claims about transcriptlevel imbalances. Indeed, some of the transcripts that contain one or more heterozygote with a significant Pvalue but were not classified as 'consensus imprinted' by [2] are clearly shown to be imprinted by our results. Note however that 27 transcripts had significant heterozygotes with imbalances in opposing directions, demonstrating that it is not always appropriate to generalize from a single locus to make claims of imbalance at the transcript level (Figure 8 and Additional file 10).
For genes containing heterozygotes with opposing imbalances, one approach is to scan the transcript annotations to identify isoform structures consistent with the observed SNP positions and imbalances. This approach was taken by [2], who defined these genes as 'complex' as long as at least one SNP was significant. An example of a complex gene is H13, which has two short isoforms and three longer isoforms with several additional exons towards the 3' end (Figure 9). The short isoforms contained heterozygotes with a paternal bias in their 3' exons while the heterozygotes on the 3' and intermediary exons of the longer isoforms had a maternal bias (cf. Figure S9 of [2] for a SNPbySNP visualization of the results of their preoptic area F_{1} samples). Using MMSEQ, we were able to discern this effect by direct quantification of haploisoforms. The two short isoforms were clearly imbalanced towards the paternally inherited haplotype while two of the long isoforms were clearly imbalanced towards the maternally inherited haplotype. An additional gene within the boundaries of H13, Mcts2, was also found to be paternally overexpressed (Table 3). By exploiting the data and annotation simultaneously, MMSEQ can be used to detect opposing imbalances between isoforms of the same gene directly.
Discussion
We have presented a pipeline and statistical method that can disaggregate expression between isoforms and even between the two haplotypes of each isoform within an individual. MMSEQ produces improved isoform estimates compared to RSEM for medium to low expression transcripts, is more scalable, and estimates standard errors more efficiently. Furthermore, our principled approach to haploisoform quantification obviates the need for adhoc interpretations of SNPbySNP imbalances in terms of transcripts. Two aspects of our method, however, deserve further discussion.
Transcript discovery
MMSEQ aims to quantify the abundance of known transcripts, and as such relies on the comprehensiveness of the transcriptome's annotation. It is usually possible to align a very large proportion of the reads to Ensembl transcripts (approximately 75% in the HapMap study using Ensembl version 56). However, samples may contain previously unobserved genes or isoforms. MMSEQ can in such cases work in tandem with transcript discovery methods by adding newly predicted isoform sequences to the reference transcript FASTA file and using it in the alignment and mapping steps of the MMSEQ workflow.
Modeling biological variability
The Poisson distribution captures technical variability arising in repeated sequencing experiments with the same biological sample. The true expression value is, in effect, fixed by the experiment, and the only source of variability arises from measurement error and mapping uncertainty. However, between biological replicates such as different individuals in the HapMap study, there is, additionally, variability of a biological origin. As has been previously reported, this results in expression values between replicates that show overdispersion, captured, for example, by a negative binomial distribution [21].
Here we have focused on the problem of estimating the posterior distribution of expression values independently per sample. Nevertheless, it would be possible to add a further level to our Bayesian model to capture overdispersion across samples flexibly. For example, if exchangeable Gamma priors are set on the μ_{ t } , a suitable negative binomial model can be induced.
Phasing with pairedend data
In this work, we have phased genotype calls obtained from SAMtools pileups  an approach that works well with both single and pairedend data. However, in the case of pairedend data, the haplotypes observed directly at multiple SNPs spanned by overlapping read pairs could be used to increase confidence in the phasing calls. Although incorporating this information would benefit phasing estimates only for some sets of SNPs, we believe it is a worthwhile area of future research. As phasing is a distinct step in our pipeline, improved methodologies can be integrated flexibly as they become available.
Conclusions
RNAseq is a promising and rapidly developing technology that provides sequence and expression intensity information of a sample in a single experiment. We have presented a novel pipeline and fast, scalable methodology to estimate expression of diploid organisms at the haplotype, isoform and gene levels. This allows researchers to go beyond allelespecific expression analysis and assess imbalance between paternal and maternal copies of isoforms, which in turn may be compared to differential isoform expression between individuals. We have shown that our method is able to deconvolve the expression of transcripts on each of two X chromosomes from human males in a pooled dataset, and that it can be successfully applied to detect genomic imprinting and cisregulated transcripts in mouse hybrids. Our method retains reads that emanate from junctions as well as wholly within exons, models alignments to multiple transcripts, potentially across genes, exploits insert size information in pairedend data to choose the best alignments and flexibly incorporates corrective models for nonuniform read sampling. The pipeline, the MMSEQ software and related documentation are freely available online [14].
Materials and methods
Expectation maximization
We augment the data with the reads per region and transcript, X_{ it } , where Σ _{ t } X_{ it } = k_{ i } and use the Poisson approximation for the augmented data likelihood:
The distribution of the augmented data conditional on the observed data and the parameters is multinomial:
The derivative of the expected Poisson log likelihood over X given k and μ^{(p) }with respect to μ_{ t } is linear in X, and hence
The EM algorithm can be thus be expressed as repeatedly updating the {\mu}_{t}^{(p)} at each iteration p using a form of the Poisson ML estimator in which the X_{ it } have been substituted with \mathbb{E}({X}_{it}{k}_{it},{\mu}_{t}^{(p)}):
which converges to the ML estimate of μ_{ t } . To initialize the algorithm, we set μ^{(0)} equal to {\scriptscriptstyle \frac{1}{b{l}_{t}}}{\displaystyle {\sum}_{i}{\scriptscriptstyle \frac{{M}_{it}{k}_{i}}{{\Sigma}_{{t}^{\prime}}{M}_{i{t}^{\prime}}}}}, which is equivalent to distributing k_{ i } evenly between cells of X_{ i } where M_{ it } is one. For a given region i, the probability of reads being allocated to a given transcript depends only on the μ_{ t } and not on s_{ i } (as the region is the same length on all transcripts). Hence, the s_{ i } do not appear in the update steps.
Bayesian model and Gibbs sampling
As before, we use the augmented data reads per region and transcript, X_{ it } . The full model is:
The full conditionals are:
Again, the s_{ i } are not needed as they are absent from the full conditionals.
TopHat settings
Gapped alignment to the genome is performed with TopHat. We use a GFF file (specified with G) based on the Ensembl annotation. We set nonoveljuncs, minisoformfraction 0.0 and minanchorlength 3. The expected inner distance between mate pairs is specified with the r switch.
SAMtools pileup settings
Genotypes output by SAMtools pileup were filtered using samtools.pl varFilter with default options and setting a minimum Phredscaled probability of the genotype being identical to the reference ('SNP quality') threshold of 20.
Bowtie settings
Alignment to the transcriptome with Bowtie is performed with the a best switches, which ensure all the best alignments in terms of mismatches are produced. Additionally, we recommend using strata to output only alignments with the minimum number of mismatches, although it currently has no effect on pairedend data. The minimum and maximum insert sizes should be set appropriately with the I and X switches respectively, as should norc/nofw for stranded protocols.
Abbreviations
 CEU:

Utah residents with ancestry from northern and western Europe
 EM:

expectation maximization
 GAII:

Genome Analyzer II
 GS:

Gibbs sampling
 Haploisoform:

haplotypespecific isoform
 MCSE:

Monte Carlo standard errors
 ML:

maximum likelihood
 PAR:

pseudoautosomal region
 RPKM:

reads per kilobase per million mapped reads
 SNP:

single nucleotide polymorphism
 UCSC:

University of California: Santa Cruz.
References
McManus CJ, Coolon JD, Duff MO, EipperMains J, Graveley BR, Wittkopp PJ: Regulatory divergence in Drosophila revealed by mRNAseq. Genome Res. 2010, 20: 816825. 10.1101/gr.102491.109.
Gregg C, Zhang J, Weissbourd B, Luo S, Schroth GP, Haig D, Dulac C: Highresolution analysis of parentoforigin allelic expression in the mouse brain. Science. 2010, 329: 643648. 10.1126/science.1190830.
Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B: Mapping and quantifying mammalian transcriptomes by RNASeq. Nat Methods. 2008, 5: 621628. 10.1038/nmeth.1226.
Richard H, Schulz MH, Sultan M, Nürnberger A, Schrinner S, Balzereit D, Dagand E, Rasche A, Lehrach H, Vingron M, Haas SA, Yaspo ML: Prediction of alternative isoforms from exon expression levels in RNASeq experiments. Nucleic Acids Res. 2010
Jiang H, Wong WH: Statistical inferences for isoform expression in RNASeq. Bioinformatics. 2009, 25: 10261032. 10.1093/bioinformatics/btp113.
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.
Hansen KD, Brenner SE, Dudoit S: Biases in Illumina transcriptome sequencing caused by random hexamer priming. Nucleic Acids Res. 2010, 38: e13110.1093/nar/gkq224.
Li J, Jiang H, Wong WH: Modeling nonuniformity in shortread rates in RNASeq data. Genome Biol. 2010, 11: R5010.1186/gb2010115r50.
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.
Degner JF, Marioni JC, Pai AA, Pickrell JK, Nkadori E, Gilad Y, Pritchard JK: Effect of readmapping biases on detecting allelespecific expression from RNAsequencing data. Bioinformatics. 2009, 25: 32073212. 10.1093/bioinformatics/btp579.
Heap GA, Yang JHM, Downes K, Healy BC, Hunt KA, Bockett N, Franke L, Dubois PC, Mein CA, Dobson RJ, Albert TJ, Rodesch MJ, Clayton DG, Todd JA, van Heel DA, Plagnol V: Genomewide analysis of allelic expression imbalance in human primary cells by highthroughput transcriptome resequencing. Hum Mol Genet. 2010, 19: 122134. 10.1093/hmg/ddp473.
Fontanillas P, Landry CR, Wittkopp PJ, Russ C, Gruber JD, Nusbaum C, Hartl DL: Key considerations for measuring allelic expression on a genomic scale using highthroughput sequencing. Mol Ecol. 2010, 19 (Suppl 1): 212227. 10.1111/j.1365294X.2010.04472.x.
Novocraft. [http://novocraft.com]
Bayesian Gene eXpression. [http://bgx.org.uk]
Goncalves A, Tikhonov A, Brazma A, Kapushesky M: A pipeline for RNAseq data processing and quality assessment. Bioinformatics. 2011
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.
Trapnell C, Pachter L, Salzberg SL: TopHat: discovering splice junctions with RNASeq. Bioinformatics. 2009, 25: 11051111. 10.1093/bioinformatics/btp120.
Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R, 1000 Genome Project Data Processing Subgroup: The sequence alignment/map format and SAMtools. Bioinformatics. 2009, 25: 20782079. 10.1093/bioinformatics/btp352.
Su SY, Asher JE, Jarvelin MR, Froguel P, Blakemore AIF, Balding DJ, Coin LJM: Inferring combined CNV/SNP haplotypes from genotype data. Bioinformatics. 2010, 26: 14371445. 10.1093/bioinformatics/btq157.
Langmead B, Trapnell C, Pop M, Salzberg SL: Ultrafast and memoryefficient alignment of short DNA sequences to the human genome. Genome Biol. 2009, 10: R2510.1186/gb2009103r25.
Anders S, Huber W: Differential expression analysis for sequence count data. Genome Biol. 2010, 11: R10610.1186/gb20101110r106.
Law AM: Confidence intervals in discrete event simulation: a comparison of replication and batch means. Naval Res Logist Q. 1977, 23: 667678. 10.1002/nav.3800240414.
UCSC Genome Browser. [http://genome.ucsc.edu]
Gregg C, Zhang J, Butler JE, Haig D, Dulac C: Sexspecific parentoforigin allelic expression in the mouse brain. Science. 2010, 329: 682685. 10.1126/science.1190831.
Puck JM, Willard HF: X inactivation in females with Xlinked disease. N Engl J Med. 1998, 338: 325328.
Acknowledgements
This research was supported by BBSRC grant BBG0003521. We thank William Astle for constructive comments on the statistical model.
Author information
Authors and Affiliations
Corresponding author
Additional information
Authors' contributions
ET and AL developed the statistical model with the advice of SR and drafted the manuscript. ET implemented the MMSEQ software and ran validation experiments. SS and ET implemented the haploisoform pipeline and ran validation experiments. LC proposed and helped develop the EM algorithm and supervised the haploisoform validation experiment. AG proposed the application to mouse crosses and provided guidance on RNAseq analysis. AG and ET applied the method to the mouse brain dataset. AL and SR supervised the project. AG, LC, SR and SS reviewed and revised the manuscript. All authors have read and approved the final manuscript.
Electronic supplementary material
13059_2010_2479_MOESM1_ESM.PDF
Additional file 1: Gibbs traces of identical transcripts. Gibbs traces for two transcripts that have identical sequences, ENST00000436491 and ENST00000415119, and their sums. The individual transcript estimates exhibit high variability and anticorrelation, but the total expression level of the two transcripts can be well estimated. (PDF 13 KB)
13059_2010_2479_MOESM2_ESM.PDF
Additional file 2: Poisson regression coefficients for three lanes in the HapMap dataset. Plots of the Poisson regression coefficients obtained using the method described in [8] from three lanes in the HapMap dataset. The first two plots are for two lanes of the same Illumina GAII run (3125_2 and 3125_7), while the last plot is for a lane in a separate run (3122_7). The coefficients are highly stable across both lanes and runs. (PDF 44 KB)
13059_2010_2479_MOESM3_ESM.PNG
Additional file 3: Plots of adjusted transcript lengths. Scatterplot of log10 true vs. adjusted transcript lengths (top) and histogram of the log10 fold change in transcript length after adjustment (bottom). The adjustments are in general very slight. (PNG 194 KB)
13059_2010_2479_MOESM4_ESM.PDF
Additional file 4: Transcript connectivity bar plot. Bar plot of the number of transcripts that each transcript is connected to via shared reads for human and mouse. (PDF 4 KB)
13059_2010_2479_MOESM5_ESM.PNG
Additional file 5: MMSEQ vs. RSEM scatterplots. Normalized simulated expression vs. log ratio between simulated and estimated normalized expression for RSEM (left) and MMSEQ GS (right) (note the difference in the scales of the yaxes). The RSEM estimates tend to underestimate some lowtomedium expression values and set them very close to zero, which translates to large negative log ratios. This also applies to MMSEQ EM estimates. The posterior means estimated using MMSEQ Gibbs sampling are less biased except for a slight upwards bias for very lowly expressed transcripts. (PNG 539 KB)
13059_2010_2479_MOESM6_ESM.PDF
Additional file 6: Quantilequantile plots between pairs of lanes of the same individual and between pairs of lanes of different individuals. Quantilequantile plots of transcript expression estimates between pairs of lanes in the HapMap dataset. The lane IDs are shown along the diagonal. The bottomleft triangle shows pairwise comparisons for a single individual sequenced in seven lanes of the same run. The upperright triangle shows pairwise comparisons between different individuals all sequenced in different lanes. There is a striking contrast in the consistency of the distribution of high values between pairs in the two triangles. (PDF 60 KB)
13059_2010_2479_MOESM7_ESM.PNG
Additional file 7: Logbase meanvariance correlation between technical and biological replicates. Scatterplots of log mean expression values against the log of the variance across technical and biological replicates at the transcript and gene levels. Each scatterplot has a line with a gradient of one if it shows technical replicates and two if it shows biological replicates. The variance is approximately proportional to the mean for technical replicates and the square of the mean for biological replicates. (PNG 679 KB)
13059_2010_2479_MOESM8_ESM.PDF
Additional file 8: Scatterplots of log expression estimates from individual and pooled data. Left: scatterplot of log expression estimates of male NA12045 vs. NA12872 obtained from individual datasets. Center: scatterplot of log expression estimates of male NA12045 obtained from the individual vs. pooled data. Right: scatterplot of log expression estimates of male NA12872 obtained from the individual vs. pooled data. (PDF 28 KB)
13059_2010_2479_MOESM9_ESM.PNG
Additional file 9: Reciprocal vs. initial cross, omitting transcripts on the X chromosome. Scatterplot of log fold changes between haploisoforms in the reciprocal (F_{1}r) and the initial (F_{1}i) cross, omitting transcripts on the X chromosome. (PNG 411 KB)
13059_2010_2479_MOESM10_ESM.PNG
Additional file 10: Reciprocal vs. initial cross, highlighting isoforms containing at least one significant SNP. Scatterplot of log fold changes between haploisoforms in the reciprocal (F_{1}r) and the initial (F_{1}i) cross, highlighting in green circles and red crosses isoforms containing at least one significant SNP imbalanced towards the paternal and maternal strain respectively. SNPs were called significant using a χ^{2} goodnessoffit test with a Pvalue threshold of 0.05 and are listed in [2]. Some transcripts contain significant SNPs with opposing imbalances, one example of which is clearly visible in the bottomright quadrant. (PNG 643 KB)
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
Rights and permissions
About this article
Cite this article
Turro, E., Su, SY., Gonçalves, Â. et al. Haplotype and isoform specific expression estimation using multimapping RNAseq reads. Genome Biol 12, R13 (2011). https://doi.org/10.1186/gb2011122r13
Received:
Revised:
Accepted:
Published:
DOI: https://doi.org/10.1186/gb2011122r13