BACKGROUND
next-generation sequencing based technology for rna profiling  has become the predominant method to quantify the transcript abundance in cells.
compared to microarray technology, rna-seq offers broader quantification range and enables the detection of novel transcripts  <cit> .
however, due to the fragmentation of sequencing material, there is greater complexity in quantification and analysis of rna-seq data  <cit> .
current state-of-the-art quantification tools for rna-seq data can be divided into two major categories  <cit> : alignment-based and alignment-free.
alignment-based quantification methods will first map each sequenced reads to a reference genome or transcriptome and then estimate the abundance of transcripts based on the alignment.
alignment-free quantification methods rely on light-weight pseudo-alignment in k-mer space to quantify the transcript abundance.
an analytic challenge raised from these quantification methods is that different method generates abundance measurements in different units, including counts, fpkm , rpkm , and tpm   <cit> .
furthermore, various transformation strategies can be applied to quantification values in purpose of specific downstream analysis like differential gene expression analyses  <cit>  or novel splicing site detection  <cit> .
although several studies have provided assessment of analysis tools for rna-seq data, little consensus on the optimal analysis pipeline is obtained .
deconvolution is a mathematical process used to extract constituent elements from a mixture of multiple signals  <cit> .
in the field of biomedical research, deconvolution is widely applied to retrieve cell-type or tissue specific gene expression profiles from heterogeneous tissue samples.
most deconvolution algorithms in the literature assume a linear model , in which the expression signal of the mixture is a weighted sum of the expression for its constitutive cell types.
previous analysis has shown the necessity of using anti-log expression microarray data to avoid unwanted bias introduced by non-linear transformation  <cit> .
however, no study has assessed the linearity of transcript abundance in rna-seq data.
therefore, in this study, we conducted a comprehensive comparison of seven rna-seq quantification methods on the linearity of the estimated abundance using a deep sequencing dataset where rna samples were mixed at known proportions.
our results will provide a good recommendation to researchers considering deconvolution on rna-seq data.
RESULTS
data
we employed the benchmark dataset used to assess rna-seq measurement performance in different application sites and platforms from the sequencing quality control  project  <cit> .
in order to have minimal inter-sample variability in the linearity evaluation analyses, we included samples from the same platform  and same sequencing center .
specifically, raw sequenced reads for four biological replicates of four types of samples  were obtained; where sample a is derived from universal human reference rna, sample b is derived from human brain reference rna, sample c is obtained by mixing a and b in ratio 3: <dig>  and sample d is obtained by mixing a and b in ratio 1: <dig>  out of  <dig> samples from a, b and c, nine samples have about eighty million pairs of raw reads and three samples have double the depth.
overall, the mappability of all the samples is around 70–80%.
a brief summary about the samples is given in additional file 6: table s <dig>
quantification methods
we performed a literature survey and selected seven prevalent quantification methods for comparison.
to increase the comparability of the estimated transcript abundance, all the alignment-based quantification methods were applied on mappings processed with tophat <dig>  <cit> .
htseq-count  <cit>  provides the number of reads/fragments mapped unambiguously to a single feature, referred as count.
cufflinks  <cit> , which is also the most popular quantification method, uses comparative algorithm assembly to produce minimal set of transcript supported by the transcript alignment.
the resulting transcript abundance is measured in fpkm.
edger  <cit>  models count data based on an overdisposed poisson model and uses an empirical bayes procedure to moderate the degree of overdispersion across genes, which is intended for downstream gene expression analysis.
rsem  <cit>  uses a generative model of rna-seq reads and the expectation-maximization  algorithm to estimate abundances of transcript feature.
the two alignment-free based quantification methods, kallisto  <cit>  and salmon  <cit>  apply pseudo-alignment to find potential transcript origins of rna-seq reads.
data distribution
to our surprise, the distribution of the estimated abundance from each quantification method results in distinct distribution at both the gene level  and isoform level .
except cufflinks fpkm, all distributions contain two sharp peaks for abundance at gene level .
while count values are normalized to tpm or fpkm, the second peak is weakened and results in a smoother curve.
although the scaling/normalization factors based on library size and gene length used in tpm and fpkm will explain the reduced range of the quantifications, it cannot explain the reduced height and smoothened second peak in the distribution.
for example, logarithm transformation reduces the range of ht-count quantifications in scale, without weakening the second peak .
the sharpness of the second peak observed at the gene level diminished at isoform level for all quantification methods .
nevertheless, the distribution pattern at both gene and isoform levels remains consistent for the same quantification method.
fig.
<dig> distribution of quantifications at gene level  and isoform level
assessment of linearity
to evaluate the linearity of the rna-seq data produced by different quantification methods, we used a multiple linear regression model: c∼m×a+n×b+ε.
since sample c is derived by mixing sample a and b at ratio of 3: <dig>  the expected values for the parameters m, n, and ε should be  <dig> ,  <dig> , and zero respectively.
we reason that if the rna-seq data are linear, the fitted multiple linear regression model will provide precise estimation for the parameters m, n, and ε. since there are four biological replicates for each sample types , there is a total of  <dig> possible models to be fitted.
we fitted all  <dig> models and studied the performance based on the average of the parameters and predicted values.
as shown in fig.
<dig>  the true value of c is linearly correlated to the fitted values from models with known mixture proportion  <dig>  and  <dig>  as parameters in all quantification methods.
the linear relationships are especially pronounced at the gene level .
although there are more data points away from the diagonal at the isoform level , the high density along the diagonal represents strong linear relationships.
results on the estimated coefficients and intercepts from the  <dig> models show that count data produce large intercept values at both the gene level and isoform level .
the intercept values are also sensitive to the counts data from different samples fitted, resulted in highly variable results across models of different samples.
fpkm values from cufflinks give the largest variation in estimated m and n among all non-count quantifications.
overall, tpm and fpkm quantifications reported by salmon, kallisto and rsem resulted in the best estimate of coefficients with small 95% confidence intervals.
similar analyses on isoform level abundance give the same results as observed in gene level data .
fig.
<dig> concordant analysis between rank of quantifications of  <dig> ×Ā+ <dig> ×b¯ and c¯  at gene level  and isoform level .
rankes were normalized by the number of quantifications in each plot
fig.
<dig> jitter boxplot of estimated coefficients and intercepts from linear model c∼m×a+n×b+ε at gene level  and isoform level .
red line indicates expected estimates if c, a and b satisfy linear assumptiotn
we observed that the estimation of m resulted in two clusters from count data : one close to true value of  <dig>  and another center around  <dig> .
we found that this phenomenon is due to the library size difference between samples.
specifically, results from cluster close to  <dig>  are from  <dig> combinations with the second replicate of sample a. from additional file 6: table s <dig>  we can find that the second replicate of a is the only sample among the four replicates to have similar library size with other samples of c and b. from this result, we could conclude that normalization of quantified abundance is essential to eliminate the inter-sample variability to meet the linear assumption of deconvolution analysis.
to further assess the linearity of the rna-seq quantifications, we evaluated the fitted model’s prediction through three analyses: 1) a concordance analysis between measured c and fitted value Ĉ, 2) a receiver operating characteristic  curve-like analysis on the absolute residual values, and 3) residual analysis for rescaled model.
similar to parameter estimation analysis mentioned above, linear regression is performed on  <dig> models, which is constructed upon combinations of  <dig> replicates in each sample, then estimated value Ĉ and residuals are averaged for final analysis.
in the concordance analysis, we evaluated if the rank of genes from the fitted value c is consistent with the rank from the true observed c. plots of the rank of c against the rank of Ĉ in fig.
4
a demonstrate that count data provide the best concordance while the cufflink’s fpkm result in the worst concordance.
in addition, the log transformation of count data induces a slight underestimation.
furthermore, the concordance of abundance estimated from kallisto and salmon enhanced tremendously as the expression level increases.
although good concordance is observed in abundances quantified by some methods at gene level , concordance of abundance at isoform levels is poor in all methods .
the poor concordance of isoform level data might due to the bias introduced by isoform abundance quantification methods.
moreover, the concordance analysis  on measured c against constructed c based on ground truth represents consistent result to the concordance analysis between c against Ĉ .
in the roc-like analysis of absolute residual values, we evaluated the proportions of genes or isoforms with residual at a given threshold level t. results on the abundance at the gene level show that all studied methods, except for log transformation of count, result in reasonable performance.
among which, tpm from salmon and kallisto perform the best .
similarly, tpm from salmon and kallisto gives the best performance at isoform level and count data perform poorly in general
fig.
<dig> concordant analysis between rank of estimated quantifications and rank of measured abundance value at gene level  and isoform level .
the fitted value in the y-axis is estimated from model c∼m×a+n×b+ε.
ranks were normalized by the number of quantifications in each plot
fig.
<dig> roc-like curve evaluating linearity of quantified abundance at gene level  and isoform level  based on residuals from model c∼m×a+n×b+ε.
proportion of variables with residuals smaller than a threshold is computed
in the residual analysis, in order to make the residual comparable across quantification methods in different ranges, we rescaled quantifications and use the model c−μcσc∼m×a−μaσa+n×b−μbσb+ε to conduct multiple linear regression.
from the residual plots, we could observe that all the methods except for cufflinks and log transformation resulted in small residuals across fitted values .
fig.
<dig> residual plot for rescaled model c−μcσc∼m×a−μaσa+n×b−μbσb+ε at gene level  and isoform level
discussion
from the results in fig.
<dig>  we observed that the performance of linear model on count data is largely affected by the library size of samples fitted.
therefore, we suggest normalization as a compulsory preprocess step prior the deconvolution analysis.
the poor estimation of isoform based data  might due to the challenge of isoform abundance quantification.
in addition to the analyses presented in the “results” section above, we also conducted the same set of analyses on sample d, using the model d∼m×a+n×b+ε.
since d is made up of a and b with proportion 1: <dig>  the expected value for m and n is  <dig>  and  <dig>  respectively in this d-based model.
all of the results from the d-based model are consistent with the findings observed from the c-based model discussed .
therefore, we can conclude that the evaluation of linearity of quantification data is not affected by the specific proportion of the original mixture.
in summary, the main objective of this study is to provide reference to the researcher considering deconvolution as downstream analysis of rna-seq data.
although linear assumption is the priority of deconvolution analysis, it does not guarantee the performance.
for future study, assessment of deconvolution performance is still required.
CONCLUSIONS
we conducted a comprehensive study to assess the linearity of gene and isoform abundance reported by different rna-seq quantification methods based on the performance how these quantifications fitted in a multiple regression linear model.
from our analysis, we observed that abundance at both the gene and isoform level from different quantification methods exhibit distinct distribution patterns and thus give diverse results.
abundance reported in the units of counts gave poor linearity, demonstrated from the worst estimated parameters and intercept values of the model.
this indicates the necessity of normalizing the abundance data prior the deconvolution analysis.
in total, tpm reported by salmon and kallisto is the best abundance data for linear models as the estimated parameters are close to true mixture proportions and the fitted values are tightly linearly correlated to the sample’s measured abundance.
moreover, when comparing within the same quantification method at both gene and isoform levels, the correlation between the measured and models’ fitted values is lower at isoform level compared to gene level while the performance on estimating parameters of linear models are similar in both levels.

 next-generation sequencing based technology for rna profiling has become the predominant method to quantify the transcript abundance in cells. compared to microarray technology, rna-seq offers broader quantification range and enables the detection of novel transcripts. however, due to the fragmentation of sequencing material, there is greater complexity in quantification and analysis of rna-seq data. current state-of-the-art quantification tools for rna-seq data can be divided into two major categories: alignment-based and alignment-free. alignment-based quantification methods will first map each sequenced reads to a reference genome or transcriptome and then estimate the abundance of transcripts based on the alignment. alignment-free quantification methods rely on light-weight pseudo-alignment in k-mer space to quantify the transcript abundance. an analytic challenge raised from these quantification methods is that different method generates abundance measurements in different units, including counts, fpkm, rpkm, and tpm. furthermore, various transformation strategies can be applied to quantification values in purpose of specific downstream analysis like differential gene expression analyses or novel splicing site detection. although several studies have provided assessment of analysis tools for rna-seq data, little consensus on the optimal analysis pipeline is obtained. deconvolution is a mathematical process used to extract constituent elements from a mixture of multiple signals. in the field of biomedical research, deconvolution is widely applied to retrieve cell-type or tissue specific gene expression profiles from heterogeneous tissue samples. most deconvolution algorithms in the literature assume a linear model, in which the expression signal of the mixture is a weighted sum of the expression for its constitutive cell types. previous analysis has shown the necessity of using anti-log expression microarray data to avoid unwanted bias introduced by non-linear transformation. however, no study has assessed the linearity of transcript abundance in rna-seq data. therefore, in this study, we conducted a comprehensive comparison of seven rna-seq quantification methods on the linearity of the estimated abundance using a deep sequencing dataset where rna samples were mixed at known proportions. our results will provide a good recommendation to researchers considering deconvolution on rna-seq data.  data we employed the benchmark dataset used to assess rna-seq measurement performance in different application sites and platforms from the sequencing quality control project. in order to have minimal inter-sample variability in the linearity evaluation analyses, we included samples from the same platform and same sequencing center. specifically, raw sequenced reads for four biological replicates of four types of samples were obtained; where sample a is derived from universal human reference rna, sample b is derived from human brain reference rna, sample c is obtained by mixing a and b in ratio 3: and sample d is obtained by mixing a and b in ratio 1: out of samples from a, b and c, nine samples have about eighty million pairs of raw reads and three samples have double the depth. overall, the mappability of all the samples is around 70–80%. a brief summary about the samples is given in additional file 6: table s quantification methods we performed a literature survey and selected seven prevalent quantification methods for comparison. to increase the comparability of the estimated transcript abundance, all the alignment-based quantification methods were applied on mappings processed with tophat. htseq-count provides the number of reads/fragments mapped unambiguously to a single feature, referred as count. cufflinks, which is also the most popular quantification method, uses comparative algorithm assembly to produce minimal set of transcript supported by the transcript alignment. the resulting transcript abundance is measured in fpkm. edger models count data based on an overdisposed poisson model and uses an empirical bayes procedure to moderate the degree of overdispersion across genes, which is intended for downstream gene expression analysis. rsem uses a generative model of rna-seq reads and the expectation-maximization algorithm to estimate abundances of transcript feature. the two alignment-free based quantification methods, kallisto and salmon apply pseudo-alignment to find potential transcript origins of rna-seq reads. data distribution to our surprise, the distribution of the estimated abundance from each quantification method results in distinct distribution at both the gene level and isoform level. except cufflinks fpkm, all distributions contain two sharp peaks for abundance at gene level. while count values are normalized to tpm or fpkm, the second peak is weakened and results in a smoother curve. although the scaling/normalization factors based on library size and gene length used in tpm and fpkm will explain the reduced range of the quantifications, it cannot explain the reduced height and smoothened second peak in the distribution. for example, logarithm transformation reduces the range of ht-count quantifications in scale, without weakening the second peak. the sharpness of the second peak observed at the gene level diminished at isoform level for all quantification methods. nevertheless, the distribution pattern at both gene and isoform levels remains consistent for the same quantification method. fig. distribution of quantifications at gene level and isoform level assessment of linearity to evaluate the linearity of the rna-seq data produced by different quantification methods, we used a multiple linear regression model: c∼m×a+n×b+ε. since sample c is derived by mixing sample a and b at ratio of 3: the expected values for the parameters m, n, and ε should be,, and zero respectively. we reason that if the rna-seq data are linear, the fitted multiple linear regression model will provide precise estimation for the parameters m, n, and ε. since there are four biological replicates for each sample types, there is a total of possible models to be fitted. we fitted all models and studied the performance based on the average of the parameters and predicted values. as shown in fig. the true value of c is linearly correlated to the fitted values from models with known mixture proportion and as parameters in all quantification methods. the linear relationships are especially pronounced at the gene level. although there are more data points away from the diagonal at the isoform level, the high density along the diagonal represents strong linear relationships. results on the estimated coefficients and intercepts from the models show that count data produce large intercept values at both the gene level and isoform level. the intercept values are also sensitive to the counts data from different samples fitted, resulted in highly variable results across models of different samples. fpkm values from cufflinks give the largest variation in estimated m and n among all non-count quantifications. overall, tpm and fpkm quantifications reported by salmon, kallisto and rsem resulted in the best estimate of coefficients with small 95% confidence intervals. similar analyses on isoform level abundance give the same results as observed in gene level data. fig. concordant analysis between rank of quantifications of ×Ā+ ×b¯ and c¯ at gene level and isoform level. rankes were normalized by the number of quantifications in each plot fig. jitter boxplot of estimated coefficients and intercepts from linear model c∼m×a+n×b+ε at gene level and isoform level. red line indicates expected estimates if c, a and b satisfy linear assumptiotn we observed that the estimation of m resulted in two clusters from count data: one close to true value of and another center around. we found that this phenomenon is due to the library size difference between samples. specifically, results from cluster close to are from combinations with the second replicate of sample a. from additional file 6: table s we can find that the second replicate of a is the only sample among the four replicates to have similar library size with other samples of c and b. from this result, we could conclude that normalization of quantified abundance is essential to eliminate the inter-sample variability to meet the linear assumption of deconvolution analysis. to further assess the linearity of the rna-seq quantifications, we evaluated the fitted model’s prediction through three analyses: 1) a concordance analysis between measured c and fitted value Ĉ, 2) a receiver operating characteristic curve-like analysis on the absolute residual values, and 3) residual analysis for rescaled model. similar to parameter estimation analysis mentioned above, linear regression is performed on models, which is constructed upon combinations of replicates in each sample, then estimated value Ĉ and residuals are averaged for final analysis. in the concordance analysis, we evaluated if the rank of genes from the fitted value c is consistent with the rank from the true observed c. plots of the rank of c against the rank of Ĉ in fig. 4 a demonstrate that count data provide the best concordance while the cufflink’s fpkm result in the worst concordance. in addition, the log transformation of count data induces a slight underestimation. furthermore, the concordance of abundance estimated from kallisto and salmon enhanced tremendously as the expression level increases. although good concordance is observed in abundances quantified by some methods at gene level, concordance of abundance at isoform levels is poor in all methods. the poor concordance of isoform level data might due to the bias introduced by isoform abundance quantification methods. moreover, the concordance analysis on measured c against constructed c based on ground truth represents consistent result to the concordance analysis between c against Ĉ. in the roc-like analysis of absolute residual values, we evaluated the proportions of genes or isoforms with residual at a given threshold level t. results on the abundance at the gene level show that all studied methods, except for log transformation of count, result in reasonable performance. among which, tpm from salmon and kallisto perform the best. similarly, tpm from salmon and kallisto gives the best performance at isoform level and count data perform poorly in general fig. concordant analysis between rank of estimated quantifications and rank of measured abundance value at gene level and isoform level. the fitted value in the y-axis is estimated from model c∼m×a+n×b+ε. ranks were normalized by the number of quantifications in each plot fig. roc-like curve evaluating linearity of quantified abundance at gene level and isoform level based on residuals from model c∼m×a+n×b+ε. proportion of variables with residuals smaller than a threshold is computed in the residual analysis, in order to make the residual comparable across quantification methods in different ranges, we rescaled quantifications and use the model c−μcσc∼m×a−μaσa+n×b−μbσb+ε to conduct multiple linear regression. from the residual plots, we could observe that all the methods except for cufflinks and log transformation resulted in small residuals across fitted values. fig. residual plot for rescaled model c−μcσc∼m×a−μaσa+n×b−μbσb+ε at gene level and isoform level discussion from the results in fig. we observed that the performance of linear model on count data is largely affected by the library size of samples fitted. therefore, we suggest normalization as a compulsory preprocess step prior the deconvolution analysis. the poor estimation of isoform based data might due to the challenge of isoform abundance quantification. in addition to the analyses presented in the “results” section above, we also conducted the same set of analyses on sample d, using the model d∼m×a+n×b+ε. since d is made up of a and b with proportion 1: the expected value for m and n is and respectively in this d-based model. all of the results from the d-based model are consistent with the findings observed from the c-based model discussed. therefore, we can conclude that the evaluation of linearity of quantification data is not affected by the specific proportion of the original mixture. in summary, the main objective of this study is to provide reference to the researcher considering deconvolution as downstream analysis of rna-seq data. although linear assumption is the priority of deconvolution analysis, it does not guarantee the performance. for future study, assessment of deconvolution performance is still required.  we conducted a comprehensive study to assess the linearity of gene and isoform abundance reported by different rna-seq quantification methods based on the performance how these quantifications fitted in a multiple regression linear model. from our analysis, we observed that abundance at both the gene and isoform level from different quantification methods exhibit distinct distribution patterns and thus give diverse results. abundance reported in the units of counts gave poor linearity, demonstrated from the worst estimated parameters and intercept values of the model. this indicates the necessity of normalizing the abundance data prior the deconvolution analysis. in total, tpm reported by salmon and kallisto is the best abundance data for linear models as the estimated parameters are close to true mixture proportions and the fitted values are tightly linearly correlated to the sample’s measured abundance. moreover, when comparing within the same quantification method at both gene and isoform levels, the correlation between the measured and models’ fitted values is lower at isoform level compared to gene level while the performance on estimating parameters of linear models are similar in both levels. 
using a benchmark rna-seq dataset, we investigated the linearity of abundance estimated from seven most popular rna-seq quantification methods both at the gene and isoform levels. linearity is evaluated through parameter estimation, concordance analysis and residual analysis based on a multiple linear regression model. salmon and kallisto tpm data gives the best fit to the linear model studied. in the field of biomedical research, deconvolution analysis is applied to obtain single cell-type or tissue specific signatures from a mixed signal and most of them follow the linearity assumption. although recent development of next generation sequencing technology suggests rna-seq as a fast and accurate method for obtaining transcriptomic profiles, few studies have been conducted to investigate best rna-seq quantification methods that yield the optimum linear space for deconvolution analysis. results show that count data gives poor parameter estimations, large intercepts and high inter-sample variability; while tpm value from kallisto and salmon shows high linearity in all analyses. deconvolution is a mathematical process of resolving an observed function into its constituent elements.
 deconvolution is a mathematical process of resolving an observed function into its constituent elements. in the field of biomedical research, deconvolution analysis is applied to obtain single cell-type or tissue specific signatures from a mixed signal and most of them follow the linearity assumption. although recent development of next generation sequencing technology suggests rna-seq as a fast and accurate method for obtaining transcriptomic profiles, few studies have been conducted to investigate best rna-seq quantification methods that yield the optimum linear space for deconvolution analysis.  using a benchmark rna-seq dataset, we investigated the linearity of abundance estimated from seven most popular rna-seq quantification methods both at the gene and isoform levels. linearity is evaluated through parameter estimation, concordance analysis and residual analysis based on a multiple linear regression model. results show that count data gives poor parameter estimations, large intercepts and high inter-sample variability; while tpm value from kallisto and salmon shows high linearity in all analyses.  salmon and kallisto tpm data gives the best fit to the linear model studied. this suggests that tpm values estimated from salmon and kallisto are the ideal rna-seq measurements for deconvolution studies. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords rna-seqdeconvolutionlinearitythird international workshop on computational network biology: modeling, analysis, and control cnb-mac 2016seattle, wa, usa 02-oct- issue-copyright-statement© the author 2017 