Mixture models reveal multiple positional bias types in RNA-Seq data and lead to accurate transcript concentration estimates
[['RNA-Seq has established itself as a popular alternative to microarrays for the quantification of RNA transcripts .', 'In contrast to microarrays , which measure the quantity of an RNA transcript by hybridization to a transcript specific oligonucleotide , RNA-Seq generates cDNA for fragments of the RNA transcript , which are sequenced by a next generation ( NGS ) sequencer .', 'One advantage of RNA-Seq over microarrays is that it does not require prior knowledge of the nucleotide sequence of the RNA transcript , which is needed to produce a transcript specific hybridization probe , and that it can therefore detect and quantify novel RNA transcripts .', 'In addition , quantification by RNA-Seq covers a wider dynamic range since microarrays suffer from signal saturation resulting in the truncation of abundance estimates for highly abundant transcripts [1] .', 'Despite these advantages , obtaining accurate transcript quantification measurements from RNA-Seq has proven difficult .', 'One of the main reasons for the inaccuracy is the failure of the statistical models used in the derivation of the measurements to properly represent biases inherent in RNA-Seq data .', 'The statistical model of the original version of Cufflinks [2] , for instance , assumes that the cDNA fragments generated by RNA-Seq are uniformly distributed along the transcripts .', 'In reality , however , this assumption is rarely fulfilled and quantification measurements by this version of Cufflinks are therefore often inaccurate .', 'One type of bias affecting transcript quantification from RNA-Seq data is the result of a preference of the fragmentation , i . e . the process that generates cDNA fragments from RNA transcripts , to produce fragments at certain positions within the transcript , e . g . at the start and/or at the end of the transcript [3] .', 'Hence , this type of bias is referred to as positional bias [4] .', 'Positional bias can also be caused by a bias in the RNA itself , for instance , due to RNA degradation which results in a shortening of the RNA .', 'Another kind of bias in RNA-Seq is introduced during ligation , amplification and NGS sequencing [5] .', 'This bias is correlated to the RNA sequence of a transcript and is therefore called sequence specific bias [4] .', 'The present article focuses on the first type of bias , i . e . the positional bias , and develops a model , Mix2 ( rd . “mixquare” ) , which learns the positional bias in RNA-Seq data .', 'In our experiments we compare Mix2 to Cufflinks [2 , 4] , eXpress [6] , RSEM [7] and PennSeq [8] both on synthetic data and on real RNA-Seq data [9] generated from the Universal Human Reference ( UHR ) and Human Brain ( HBR ) samples of the Microarray Quality Control ( MAQC ) experiment [10] .', 'The inclusion of bias models into the statistical models of RNA-Seq data has been investigated before .', 'In [11] a model is proposed to account for the variability in read counts depending on the sequence surrounding the start of a fragment .', 'The intention is similar to that of the fragment specific bias model [4] , which has been implemented as an extension to Cufflinks [2] .', 'In addition , [4] includes a non-parametric positional bias model , which can theoretically be trained with the EM algorithm .', 'However , due to the large number of variables , this is only feasible for few transcript length dependent classes , for which statistics are collected in a small number of positional bins .', 'As a result , [4] implements a positional bias model depending exclusively on the length of a transcript .', 'Similar to [4] the generative model of RSEM [7 , 12] uses a hidden variable for the positional bias , where the latter is estimated from the global bias observed in the complete RNA-Seq data set .', 'Also in RSEM , therefore , does the positional bias model depend exclusively on the transcript length .', 'The generative model of eXpress [6] differs from Cufflinks mainly in the order of fragment-length selection and implements an online rather than a batch EM algorithm .', 'The implementation described in [6] is further restricted to a sequence specific bias , with a uniform positional bias similar to Cufflinks .', 'The model developed in PennSeq [8] is again non-parametric and the large number of variables makes its training computationally prohibitive .', 'For this reason , the bias model of PennSeq [8] is not included in the parameter update but is approximated by the overall bias in a gene locus and by the transcript specific reads .', 'The method described in [13] is a model for gene read counts , which models bias by exon specific weights , which are estimated both for the complete data set and for individual genes .', 'In [14] the authors focus on RNA-Seq data with 5’ bias which is the result of RNA degradation and use an exponential model for the fragment distributions .', 'The model proposed in [15] is , again , a model for the read counts of a gene .', 'Here the read counts are modelled by a quasi-multinomial distribution with a parameter that can be adapted to account for over and under dispersion .', 'Mix2 is , similar to [2 , 4 , 6–8 , 12] , a generative model for the probability of a fragment in an RNA-Seq data set .', 'In comparison , however , the model for the positional fragment bias in Mix2 is parametric .', 'This considerably simplifies its implementation removing the need for any restrictions of the non-parametric methods .', 'At the same time , the model of the positional fragment bias in Mix2 is very versatile since mixtures of probability distributions can approximate distributions of arbitrary complexity .', 'Section ( Materials and methods ) develops the theory of Mix2 in greater detail and shows how its parameters can be updated with the EM algorithm leading to simultaneous estimates for transcript abundances and transcript specific positional fragment biases .', 'Section ( Experiments on artificial data ) optimizes the number of mixture components of the Mix2 model and compares it with Cufflinks , RSEM and eXpress on artificial data sets .', 'Sections ( Experiments on the Microarray Quality Control ( MAQC ) data ) and ( Experiments on the Sequencing Quality Control ( SEQC ) data ) on the other hand , discuss experiments on two publicly available real RNA-Seq data sets with Mix2 , RSEM , eXpress , Cufflinks and PennSeq .', 'These experiments show that , in comparison to the other methods , Mix2 leads to better correlation between estimated and measured transcript concentrations , correct recovery of mixing ratios and yields consistent titration orders .', 'In addition , we show that the Mix2 concentration estimates are repeatable across lanes and laboratory sites and lead to a more accurate detection of differentially expressed transcripts .', 'In addition , Section ( Types of bias in the MAQC data ) shows that Mix2 can be used as an explorative tool to detect positional biases present in RNA-Seq data .', 'Mix2 has been implemented as an Octave script with readable code and as a closed source C++ implementation .', 'We used the latter for the majority of our experiments .', 'Both versions can be downloaded from https://www . lexogen . com/mix-square-scientific-license .', 'While the C++ version is considerably faster , we show at the end of Section ( Experiments on artificial data ) and in Fig E in S2 Appendix that quantification results for both implementations are virtually identical .', 'Hence , either version can be used to evaluate the accuracy of Mix2 .'], ['An essential part of Next Generation Sequencing ( NGS ) is the library preparation .', 'This process takes an RNA sample and produces a library of short cDNA fragments , each corresponding to a section of an RNA transcript .', 'The cDNA fragments are sequenced by an NGS sequencer resulting in single or paired end reads which are mapped to a reference genome .', 'Hence , the probability p ( r ) of a fragment r can be interpreted as the probability of its genomic coordinates .', 'In a genomic locus the probability p ( r ) is the superposition of the fragment distributions p ( r|t =', 'i ) for the N transcripts in the locus , i . e .', 'p', '(', 'r', ')', '=', '∑', 'i', '=', '1', 'N', 'α', 'i', 'p', '(', 'r', '|', 't', '=', 'i', ')', '( 1 )', 'where αi is the relative abundance of transcript t = i , i . e . the probability that transcript t = i generates any fragment , and p ( r|t =', 'i ) is the probability that transcript t = i generates fragment r .', 'Hence p ( r|t =', 'i ) models the transcript specific fragment bias .', 'An estimate for the concentration of transcript t = i is obtained by normalizing the relative abundance αi , yielding the RPKM [16] or FPKM values [2] .', 'Mix2 uses a mixture model for p ( r|t =', 'i ) , i . e .', 'p ( r | t = i ) = ∑ j = 1 M β i j p ( r | t = i , b = j ) ( 2 )', 'where p ( r|t = i , b =', 'j ) are the M components of the mixture and the βij are the non-negative component weights .', 'Hence , p ( r|t = i , b =', 'j ) is a probability distribution over r and', '∑ j = 1 M β i j = 1 .', '( 3 )', 'Since p ( r ) is itself a mixture of the p ( r|t =', 'i ) with weights αi , this implies that p ( r ) is a mixture of mixtures motivating the name of Mix2 .', 'For the p ( r|t = i , b =', 'j ) we use Gaussians placed equidistantly along the transcript t = i ( Materials and methods', ') .', 'An example for such a Mix2 model can be found in Fig 1 . Fig 1', '( a ) shows the mixture weights βij whereas Fig 1', '( b ) shows the weighted Gaussians , βijp ( r|t = i , b =', 'j ) , and the sum of the weighted Gaussians , p ( r|t =', 'i ) .', 'The distributions in Fig 1', '( b ) are given in transcript coordinates for a transcript of 2000 bp length , while the longer dashed curve in Fig 1', '( c ) shows p ( r|t =', 'i ) in genome coordinates .', 'The locus in Fig 1', '( c ) contains two transcripts which share a common junction .', 'The shorter of the transcripts in Fig 1', '( c ) has the same set of βij as in Fig 1', '( a ) but is only 1000 bp long .', 'The relative abundances of the long and short transcript are 0 . 7 and 0 . 3 , respectively , which results in the overall distribution p ( r ) given by the solid curve in Fig 1', '( c ) .', 'In comparison , Cufflinks [2] can , for this locus , only model fragment start distributions p ( r|t =', 'i ) as visualized by the dashed curves in Fig 1', '( d ) and is therefore inappropriate for 5’ biases as the one in Fig 1', '( c ) .', 'We use the Expectation Maximization ( EM ) algorithm ( Materials and methods ) to learn the parameters of Mix2 from RNA-Seq data .', 'This results in simultaneous estimates of the transcript abundances αi and the mixture weights βij , hence the transcript specific fragment distributions p ( r|t =', 'i ) .', 'Mix2 is identifiable in most cases and can otherwise be easily made identifiable ( see Section 1 . 2 in S1 Appendix ) .', 'It can therefore always be ensured that the EM algorithm converges to the unique maximum likelihood solution .', 'It should be pointed out that we did not check for identifiability in our implementation and , like Cufflinks , returned the maximum likelihood solution produced by our method .', 'The mixture weights βij determine the shape of the fragment distribution of transcript t = i .', 'Thus , if transcripts have a similar distribution they should share the same βij .', 'This results in their fragment distributions being identical and reduces the number of parameters in Mix2 making it less prone to over-fitting .', 'Consider , for instance , the fragment start distribution of the Cufflinks model in Fig 2', '( a ) .', 'Here the distributions are similar for transcripts with 2000 bp and 3000 bp length and for transcripts with 700 bp and 1000 bp length .', 'In this situation therefore , these four transcripts can be separated into two groups where the transcripts within each group share the same mixture weights βij .', 'In general , this leads to the scenario where each transcript t = i has an associated group g = k and the distributions p ( r|t =', 'i ) of transcripts within this group share the same βij .', 'Multiple factors might influence the similarity of fragment start distributions .', 'Here , we investigate gene membership and , as the bias correction methods in [4 , 12] , transcript length .', 'The rationale for choosing these two properties is that even if fragments are uniformly distributed immediately after fragmentation , fragment size selection introduces the transcript length dependent bias in Fig 2', '( a ) .', 'On the other hand , transcripts belonging to the same gene can share a substantial part of their sequence and exhibit therefore potentially similar fragmentation properties .', 'In the following , we refer to tying between all transcripts within a gene as global tying and to tying between all transcripts within a gene and the same length range as group tying ( Materials and methods ) .', 'Hence , Fig 1', '( c ) shows an example for a Mix2 model with global tying since both transcripts share the same set of weights βij .', 'In our experiments we do not treat multi-mapping and uniquely mapping reads differently and instead consider each read mapping to be associated with a unique fragment .', 'Multi-mapping reads and sequence specific bias can , however , be integrated into Mix2 as shown in Section 1 . 1 in S1 Appendix .', 'This section pursues two goals .', 'First , we find a sensible number of mixture components for each of the variants of Mix2 .', 'Second , we investigate the performance of Mix2 under conditions favoring other quantification methods .', 'The number of mixture components derived in this section is used as a default by Mix2 .', 'While this might be suboptimal in some cases it is a necessary compromise since reference measurements of isoform concentrations on which to optimize the number of mixture components are usually not available in RNA-Seq data .', 'In our first experiments we studied the 4 transcript length dependent biases in Fig 2 . These resemble the biases we detected in our experiments on real RNA-Seq datasets MAQC and SEQC .', 'The most dominant biases in the MAQC data are visualized in Fig 3 . The 5’ biased fragment distributions in Fig 3', '( a ) and 3', '( b ) , for instance , resemble the biases in Fig 2', '( b ) and the biases for short transcripts in Fig 2', '( a ) .', 'Similarly , the fragment distributions concentrated on the 3’ side in Fig 3', '( e ) resemble the biases in Fig 2', '( c ) .', 'While the biases in Fig 3 do not depend on transcript length as strongly as the biases in our artificial data , we see with the exception of Fig 3', '( d ) an increase in the average transcript length with increasing 3’ bias .', 'Overall , around 20 . 16% of transcripts have a 5’ bias in the MAQC data , 26 . 34% have a 3’ bias and 26 . 92% a uniform fragment distribution .', 'Hence , the biases in our artificial data represent a sensible starting point for evaluating the accuracy of quantification methods under real-life conditions .', 'Fig 2 also illustrates how an incorrect choice of bias model can affect the accuracy of quantification estimates .', 'If the 5’ bias in Fig 2', '( b ) is observed for a transcript of length 3000 bp but the bias model in Fig 2', '( a ) is used for quantification then the incorrect model will try to explain the bias in Fig 2', '( b ) with a shorter transcript , thus potentially leading to an overestimate of the concentration of the short transcript .', 'In our experiments we used a set of 7 test genes from the GRCh37/hg19 Ensembl annotation v75 containing between 4 and 15 transcripts as well as the main variants of differential splicing , see Table A in S2 Appendix .', 'While this set of genes might seem small , considering all parameter combinations resulted in 538k experiments .', 'The difference between true and estimated relative isoform abundances is measured with the L1 distance , which is the sum of the absolute differences between true and estimated relative isoform abundances .', 'This value lies between 0 and 2 . In our experiments in Fig 4', '( a ) we addressed the question if the optimal number of mixture components for Mix2 depends on the gene , the number of fragments in the gene or the fragment bias .', 'In the first three bars of Fig 4', '( a ) we optimized the number of mixture components for each combination of bias , gene and sample size separately and plotted the average L1 distance over all experiments for the three types of Mix2 .', 'In the second group of bars we optimized the number of mixture components for each combination of gene and sample size ignoring the bias , while in the third group we optimized for each gene ignoring bias and sample size .', 'Finally , for the last three bars in Fig 4', '( a ) we selected a single number of mixture components for all experiments independent of gene , bias and sample size .', 'Overall , Fig 4', '( a ) shows that little is gained by optimizing for each gene , sample size and bias separately .', 'The number of mixture components of Mix2 can be chosen independent of these factors .', 'Fig 4', '( b ) shows the influence of the number of mixture components on the quantification accuracy of Mix2 .', 'For each number of mixture components on the x-axis the average L1 distance is given for different sample sizes .', 'Each group of 4 blocks contains results for 500 , 1000 , 5000 and 10000 fragments , where the different groups show results for the different types of Mix2 .', 'Fig 4', '( b ) implies that the minimal L1 distance for Mix2 without tying is obtained for 3 mixture components , while for Mix2 with group and global tying the minimum is obtained with any number between 4 and 10 .', 'In all our subsequent experiments we therefore chose 3 mixture components for Mix2 without tying and 4 mixture components for Mix2 with global and group tying .', 'We further studied the convergence of the fragment start distributions p ( r|t =', 'i ) estimated by Mix2 to the correct start distributions from which the data was sampled .', 'For this purpose , we calculated the average L1 distance between estimated and true distribution for the transcripts within a gene weighted by the estimated relative abundance of the transcript .', 'Fig 4', '( c ) and 4', '( d ) visualize the difference between the L1 distance before ( tall bars ) and after estimation of the fragment start distributions .', 'The initial fragment start distribution in Mix2 is close to uniform and independent of tying structure and read depth ( Materials and methods ) .', 'Fig 4', '( c ) shows that estimates of the fragment start distributions in Mix2 converge to the correct solution since the L1 distance for the 7 test genes decreases during the course of the parameter estimation by between 50% and 65% .', 'Mix2 with group tying and a read depth of 10000 fragments achieves in many cases the smallest final L1 distance .', 'In Fig 4', '( d ) the relative decrease in L1 distance is considerably smaller for the Cufflinks bias as the initial fragment start distribution in Mix2 is already close to the uniform Cufflinks distribution .', 'Overall , the final L1 distance is similar for all 4 biases .', 'Next , we compared the accuracy of Mix2 , Cufflinks , RSEM , and eXpress on the full transcriptome .', 'We generated an exponential transcript expression profile from an exponential gene expression profile by FluxSimulator [17] , as shown in Fig F and Fig G in S2 Appendix , and divided the gene reads between the gene transcripts according to relative abundances drawn from a uniform Dirichlet distribution .', 'In total we generated around 50 mio 100 bps read pairs for each bias in Fig 2 , which we aligned with Tophat2 .', 'Data of this kind correspond to the positional bias models in Cufflinks [4] and RSEM which are derived by scaling a single prototype distribution to the transcript length .', 'Our experiments therefore show whether Mix2 can learn data structures hard-coded into the statistical models of Cufflinks and RSEM .', 'Fig 4', '( e ) shows that on 5’ and 5’+3’ bias RSEM slightly outperforms Mix2 with group tying , while for 3’ bias Mix2 yields the best results .', 'On Cufflinks bias all methods apart from eXpress perform similarly .', 'Fig 4', '( e ) further contains results for Mix2 with initial fragment start distributions before estimation of mixture weights βij ( Mix2 init ) .', 'This shows that adapting the mixture weights improves relative abundance estimates .', 'For Cufflinks bias this improvement is minor since the initial fragment start distributions in Mix2 are already close to the correct solution .', 'The bad performance of eXpress is likely due to the fact that it models only sequence-specific and not positional bias .', 'Fig 4', '( f ) shows that Mix2 without any tying and with group tying perform consistently across biases and slightly better than RSEM with mean L1 distance of 0 . 36 versus 0 . 39 .', 'For Mix2 without parameter tying this is remarkable , as its statistical model is completely assumption free with regards to the nature of fragment bias .', 'We also performed experiments for a smaller read-depth of 5 mio read pairs .', 'These experiments which are summarized in Fig H in S2 Appendix show an increase of L1 distance for all methods but otherwise a similar trend as the experiments with 50 mio read pairs .', 'Finally , we compared the accuracy of the Octave and C++ version of our code , since on larger data sets we used , for efficiency reasons , the closed source C++ version .', 'We found on our 7 test genes that the median L1 distance between the relative abundances estimated by the two versions was between 0 . 01 and 0 . 004 for 500 and 10k reads per gene .', 'The boxplots for these experiments can be found in Fig E in S2 Appendix .', 'Due to the small difference , either version of our code can be used to evaluate the accuracy of Mix2 .', 'This section compares the accuracy of Mix2 , Cufflinks [4 , 18] , PennSeq [8] , RSEM [7] and eXpress [6] on two publicly available real RNA-Seq data sets generated from the Universal Human Reference ( UHR ) RNA and human brain ( HBR ) RNA of the Microarray Quality Control ( MAQC ) data [10] ( Materials and methods ) .', 'Both Cufflinks and RSEM were run with bias correction .', 'In accordance with our experiments on artificial data we used three mixture components for Mix2 with no tying and 4 mixture components for Mix2 with global and group tying .', 'We treat the qPCR measurements as a reference and compare them to the FPKM values generated by the quantification methods .', 'The qPCR concentrations in the MAQC data are unevenly distributed spanning several orders of magnitude .', 'It is therefore customary to compress the concentrations by taking the logarithm .', 'This achieves a more even distribution but leads , on the other hand , to outliers for FPKM values close to zero whose logarithm approaches minus infinity .', 'In order to reduce the influence of these outliers on the quality metrics for the quantification methods it is necessary to either remove transcripts with small FPKM values or to truncate the logarithm of their FPKM values .', 'We chose the latter strategy since the former ignores the fact that some quantification methods fail to detect highly abundant transcripts and furthermore leads to test sets varying significantly in size reducing the comparability of the quality metrics .', 'In our experiments , the logarithm of an FPKM value was truncated if it was below the first quartile minus 1 . 5 times the interquartile range of the logarithms of the method’s FPKM values on the test set .', 'This threshold , which was between 0 . 01 and 0 . 001 in our experiments , is the point below which values are often considered to be outliers in boxplots .', 'We use the logarithm to the basis of 10 in our experiments .', 'In order to evaluate the quantification methods on more recent data than MAQC , we also performed experiments on the SEQC data set [25] .', 'The latter contains 100 bps paired-end RNA-Seq reads generated with Illumina HiSeq 2000 at 6 laboratory sites and from four different RNA samples .', 'Samples A and B correspond to UHR and HBR from the MAQC data set .', 'Sample C and D were created by mixing A and B in 3:1 and 1:3 ratios , respectively .', 'This allows tests for titration order consistency and the correct recovery of mixing ratios .', 'These tests are independent of a “gold standard” such as qPCR which is biased by its own technical limitations .', 'In this section , therefore , we evaluate the accuracy of quantification with respect to the built-in ground truths of SEQC rather than with respect to qPCR measurements .', 'For this purpose , we downloaded RNA-Seq data for samples A to D for laboratory site BGI .', 'To study repeatability across sites we further downloaded RNA-Seq data for sample A and sites AGR , CNL and COH .', 'We compare the three variants of Mix2 with Cufflinks , RSEM and eXpress .', 'PennSeq had to be excluded as it failed to produce any output on the SEQC data .', 'As before , we used GRCh37/hg19 and Ensembl annotation version 75 in our experiments and 3 mixture components for Mix2 without tying and 4 mixture components for Mix2 with global and group tying .'], ['This article introduced Mix2 which uses a mixture of probability distributions to model transcript specific positional fragment distributions in RNA-Seq data .', 'Due to the flexibility of mixture models , Mix2 can adapt to multiple positional fragment biases of arbitrary complexity .', 'The parameters of Mix2 are efficiently trained with the EM algorithm resulting in simultaneous estimates for fragment distributions and relative abundances .', 'In addition , parameters of Mix2 can be tied between transcripts with similar fragment distribution leading to improved estimates of the relative abundances .', 'Even though Mix2 accommodates a sequence specific bias , we currently implement only a model for positional bias .', 'Sequence specific bias can , however , be a prominent feature in RNA-Seq data .', 'In [11] , for instance , it was found that a linear model on the sequences of 3 RNA-Seq data sets accounted for between 40% and 50% of the variance in sequence frequencies .', 'Implementing a model for the sequence specific bias to Mix2 might therefore result in further improvements in the accuracy of transcript quantification .', 'Experiments were conducted on artificial data to determine the optimal number of mixture components of Mix2 .', 'These experiments showed that optimization can be performed independent of gene , bias and sample size and that the optimal number of mixture components is 3 for Mix2 without tying and any number between 4 and 10 for Mix2 with tying .', 'These numbers are fairly small given the biases in Fig 2 .', 'In particular , it seems implausible that Mix2 with only 3 mixture components should be able to accurately model the 5’+3’ bias in Fig 2 ( d ) .', 'This suggests that the potential of Mix2 has yet to be fully exploited .', 'Our experiments , further , showed that the estimate of Mix2 for the transcript specific fragment bias converges to the correct distribution and that these estimates can therefore be used to detect positional bias present in RNA-Seq data .', 'Experiments were also performed on RNA-Seq data generated from Universal Human Reference ( UHR ) RNA and Human Brain ( HBR ) RNA for the Microarray and Sequencing Quality Control ( MAQC , SEQC ) data sets .', 'On MAQC , we obtained improved correlation between qPCR measurements and quantification estimates with Mix2 , while on the SEQC data set , Mix2 produced improved titration order consistency and recovery of mixing ratios .', 'In addition , correlation and standard deviation of Mix2 quantification estimates were superior across lanes in MAQC and across laboratory sites in SEQC , implying reduced sensitivity to technical variance .', 'Furthermore , correlation of qPCR and FPKM fold changes between UHR and HBR on MAQC were noticeably higher for the Mix2 model than for the other methods .', 'We showed in a classification experiment that this leads to higher accuracy in the detection of differential expression .', 'In general , quantification accuracy affects sensitivity and specificity of statistical tests for differential expression .', 'Inconsistent quantification estimates , for instance , increase technical variability and lead to an increase in the variance of the data models underlying statistical tests .', 'This , in turn , leads to a decrease of the chi-square distributed test statistic of the Wald test in DSS [27] , DESeq2 [23] , edgeR [22] and the score test in PoissonSeq [28] resulting in a loss of sensitivity .', 'The same holds for the Wilcoxon test in SamSeq [32] where the test statistic moves closer to the mean of the test distribution under the 0-hypothesis .', 'Likewise , conditional model probabilities in baySeq [29] and conditional probabilities in the exact tests of DESeq [30] and edgeR [31] become more uniform and therefore less distinctive between different conditions .', 'The noise distributions in NOISeq [33 , 34] and DEGSeq [36] become wider with increasing technical variability and fold change thresholds for the detection of differential expression increase , again reducing the sensitivity of these methods .', 'For DEXUS [35] the variance of the major and minor conditions increases resulting in greater overlap .', 'Since the prior of the model probabilities favors a single condition this makes it less likely that minor conditions and therefore differential expression will be detected .', 'Overall , the aforementioned statistical tests benefit both in terms of sensitivity and selectivity from more accurate quantification estimates and we expect therefore to see improved differential expression calls for transcripts when using these tests in conjunction with Mix2 .', 'In terms of resource usage , both RSEM and Mix2 take about the same time to process the 7 lanes of UHR and HBR in MAQC and are both faster than eXpress and Cufflinks .', 'Memory usage of RSEM is slightly smaller than that of Mix2 and Cufflinks but memory consumption of all 4 methods is low given the specifications of current computing environments .', 'In contrast to the experiments on artificial data , experiments on MAQC and SEQC showed a degradation in the performance of Mix2 when tying parameters .', 'This cannot be attributed to a suboptimal choice of parameters in the transcript clustering procedure .', 'Instead , it seems that positional fragment bias does not exclusively depend on gene membership and transcript length .', 'This fact was also highlighted in our experiments on bias types in MAQC and SEQC .', 'In these experiments we found dominant bias types by clustering fragment start distributions estimated by Mix2 .', 'On MAQC we obtained 6 clusters , containing 73 . 43% of the distributions , of which 5 clusters exhibited non-uniform distributions .', 'The cluster with uniform distributions contained only 26 . 92% .', 'On SEQC we also see the majority of distributions located in non-uniform clusters .', 'In addition , clusters are similar across laboratory sites .', 'Contrary to our experiments on artitifical data there is no obvious relationship between bias and transcript length , although correlations between the two do exist .', 'For instance , transcripts whose fragment start distributions are 3’ biased or uniform tend to be longer , whereas 5’ biased transcripts tend to be shorter .', 'A more detailed analysis of biases might reveal relations between positional bias and RNA sequences that will lead to a better tying strategy for the Mix2 model on real RNA-Seq data .', 'In summary , Mix2 can be used as an explorative tool to investigate the positional biases present in RNA-Seq data and thereby study the influence of library preparation , sequencing and data processing on the accuracy of transcript concentration estimates .', 'In addition , and more importantly , our results show that Mix2 yields improved transcript concentration estimates for RNA-Seq data with higher repeatability for technical replicates and leads , furthermore , to improved accuracy in the detection of differential expression .'], ['We factorize p ( r|t = i ) as follows ,', 'p', '(', 'r', '|', 't', '=', 'i', ')', '=', 'p', '(', 'l', '(', 'r', ')', '|', 't', '=', 'i', ',', 's', '(', 'r', ')', ')', '∑', 'j', 'β', 'i', 'j', 'p', '(', 's', '(', 'r', ')', '|', 't', '=', 'i', ',', 'b', '=', 'j', ')', '( 4 )', 'where s ( r ) and l ( r ) are the start and length of fragment r and 
['Introduction', 'Results', 'Discussion', 'Materials and methods']
['Accuracy of transcript quantification with RNA-Seq is negatively affected by positional fragment bias .', 'This article introduces Mix2 ( rd . “mixquare” ) , a transcript quantification method which uses a mixture of probability distributions to model and thereby neutralize the effects of positional fragment bias .', 'The parameters of Mix2 are trained by Expectation Maximization resulting in simultaneous transcript abundance and bias estimates .', 'We compare Mix2 to Cufflinks , RSEM , eXpress and PennSeq; state-of-the-art quantification methods implementing some form of bias correction .', 'On four synthetic biases we show that the accuracy of Mix2 overall exceeds the accuracy of the other methods and that its bias estimates converge to the correct solution .', 'We further evaluate Mix2 on real RNA-Seq data from the Microarray and Sequencing Quality Control ( MAQC , SEQC ) Consortia .', 'On MAQC data , Mix2 achieves improved correlation to qPCR measurements with a relative increase in R2 between 4% and 50% .', 'Mix2 also yields repeatable concentration estimates across technical replicates with a relative increase in R2 between 8% and 47% and reduced standard deviation across the full concentration range .', 'We further observe more accurate detection of differential expression with a relative increase in true positives between 74% and 378% for 5% false positives .', 'In addition , Mix2 reveals 5 dominant biases in MAQC data deviating from the common assumption of a uniform fragment distribution .', 'On SEQC data , Mix2 yields higher consistency between measured and predicted concentration ratios .', 'A relative error of 20% or less is obtained for 51% of transcripts by Mix2 , 40% of transcripts by Cufflinks and RSEM and 30% by eXpress .', 'Titration order consistency is correct for 47% of transcripts for Mix2 , 41% for Cufflinks and RSEM and 34% for eXpress .', 'We , further , observe improved repeatability across laboratory sites with a relative increase in R2 between 8% and 44% and reduced standard deviation .']
['RNA-Seq is a powerful tool for detecting and quantifying genes and gene isoforms .', 'However , accurate quantification in genomic loci with multiple isoforms has proven difficult .', 'This is due to the fact that the transcript generating an RNA-Seq fragment cannot be identified if multiple transcripts share the fragment sequence .', 'Due to this ambiguity , transcript concentration is usually determined in a statistical framework by calculating the probability that a transcript generates an RNA-Seq fragment .', 'Accurate estimation of this probability requires an accurate model of the transcript specific distributions of RNA-Seq fragments .', 'However , fragment distributions in statistical models of RNA-Seq data are usually over-simplified .', 'This article introduces the Mix2 ( rd . “mixquare” ) model which uses mixtures of probability distributions to model the transcript specific positional fragment distributions .', 'Mix2 learns the mixture weights and approximates therefore the fragment bias in RNA-Seq data .', 'We compare Mix2 on artificial and real RNA-Seq data to four state-of-the-art quantification methods .', 'Our experiments show that Mix2 yields more accurate and repeatable quantification estimates and that it leads to more accurate detection of differential expression .', 'We further show that the biases detected by Mix2 contradict the common assumption of a uniform fragment distribution .']
['sequencing techniques', 'engineering and technology', 'molecular probe techniques', 'probability distribution', 'mathematics', 'statistics (mathematics)', 'industrial engineering', 'quality control', 'forms of dna', 'bioassays and physiological analysis', 'molecular biology techniques', 'dna', 'rna sequencing', 'research and analysis methods', 'probe hybridization', 'rna hybridization', 'molecular biology', 'probability theory', 'complementary dna', 'microarrays', 'biochemistry', 'statistical models', 'nucleic acids', 'genetics', 'biology and life sciences', 'physical sciences', 'statistical data']