BACKGROUND
transcriptome analysis is integral in understanding the role of genetic and genomic variants in disease progression, classification, and diagnosis.
in recent years, next generation sequencing  technologies have catapulted transcriptome profiling to new heights by providing greater precision of cellular rna content.
as a result, rna-seq is more sensitive to subtle changes in expression levels and typically includes many transcripts with low read count integers.
researchers have taken interest in the low regions of a sample for several reasons, all with the primary intent of quantifying the level of expression that reflects important biological meaning.
one of the major problems in rna-seq data analysis is that there is no consensus on what is considered low expression.
low counts have been referenced as less than  <dig> reads when summed across treatments
<cit> , less than  <dig> reads on average
<cit> , less than  <dig> reads on average
<cit> , and less than  <dig> reads
<cit> .
additionally, count-based pre-filtering methods have been adopted to exclude genes with minimal expression from differential testing.
for example, risso et al.
<cit>  filtered out genes with an average count below  <dig>  robinson et al.
<cit>  recommended removing genes that did not have at least  <dig> counts per million reads in at least two samples in the edger user’s guide.
in the deseq user’s guide, anders
<cit>  took a more aggressive approach and suggested removing up to 40% of genes that ranked lowest in regard to total count across all experimental samples.
thus, at present, whether low expressed transcripts are simply identified or deliberately omitted from analysis, methods of differentiating spurious rna-seq data from meaningful information have not been explicitly defined.
furthermore, it is not clear whether it is justifiable to extrapolate the aforementioned cutoffs to other rna-seq data.
a number of factors affect the distribution of read counts for a given study.
specifically, the manufacturer, library preparation and construction, alignment algorithm, gene length, sequencing depth, and experimental design all play a role in determining the number of reads that are mapped to a gene.
since each of these factors varies from study to study, it may be misleading to ignore properties of the sample distribution by applying an arbitrary cutoff to classify the low expressed region of a sample.
to address the question of what should be considered the lower bound of functional gene expression, it is useful to consider the premise that genes can be categorized into a group of high expressed , meaningful genes and a group of low expressed , non-informative genes.
in previous studies, the empirical distribution of a sample was used to identify classes of mrna abundance levels
<cit> .
this concept has also been utilized to differentiate functional expression states of microarray data
<cit> .
in order to optimally separate expression classes, a priori knowledge of the expression distribution is useful.
in many cases, the precise distribution of noise is unknown; however, characterizations of a global bimodal distribution and a normally distributed he component have been reported in previous studies
<cit> .
hebenstreit et al.
<cit>  studied the global distribution of rna sequences in mice th <dig> cells.
their work demonstrated that log2-transformed reads per kilobase per million  could be separated into two classes of mrna abundance.
the researchers used the likelihood ratio test
<cit>  to determine whether the data was best modeled as a mixture of n or n +  <dig> components, where  <dig> < n <  <dig>  mclust
<cit> , an r package that uses an expectation-maximization  algorithm, was used to compute maximum likelihood estimates of all parameters of the mixture distribution.
for one-dimensional data, mclust evaluates a fitted model with equal variance terms and a fitted model with unequal variance and selects the model with higher bayesian information criterion
<cit> .
however, the use of mixture modeling to identify mixture components of rna-seq data may not be appropriate for two reasons.
first, it is difficult to capture the true number of mixture components.
it is known that the em algorithm performs well at estimating the parameters of a finite mixture model; however, when mixture distributions are unimodal and there is no clear separation between components, the em algorithm commonly returns more components than what seems logical based on visual inspection of the data
<cit> .
second, when a two-component mixture model is forced to fit data that does appear to be bimodally distributed, the fitted model does not always approximate very well the observed empirical distribution.
hebenstreit et al.
<cit>  were able to characterize the mixture distribution of the le and he components and determine the peak of the le region by using a poisson distribution to estimate the proportion of undetected genes at each expression level.
however, they did not discuss or provide an expression cutoff for separating the two overlapping regions.
in this study, we propose dafs, a data-adaptive method for identifying and subsequently flagging expressions in the le region by estimating the lower bound of high expression in a given rna-seq sample.
in light of the drawbacks of the mixture modeling approach, dafs was constructed without imposing a finite mixture model on the data.
several real rna-seq datasets and simulated data are used to present our findings and demonstrate the robustness of our methodology.
RESULTS
when the le and he regions of rna-seq data are well separated , dafs and mclust separate the components at similar cutoffs .
cutoff values determined by dafs and the theoretical intersection of the two mixture components obtained by mclust were - <dig>  and - <dig> , respectively.
in their study of th <dig> cells, hebenstreit et al.
<cit>  identified a sample of known expressed and unexpressed genes in th <dig> cells.
based on their classification of low and high expression, they mapped all the known expressed genes to the he component and mapped most of the known unexpressed genes to the le component .
we were able to reproduce their findings using the cutoff estimated by dafs.
analysis of bulk rna positive controls generated from cultured hct- <dig> cells
<cit>  allowed for an additional assessment of the ability of dafs to separate known expressed genes into one mixture component and known unexpressed genes into a separate mixture component.
in their analysis of single-cell whole transcriptome amplification, the authors identified  <dig> known expressed and unexpressed genes in hct- <dig> cells .
fragments transformed by fragments per kilobase of exon per million mapped reads
<cit>  were downloaded from the geo database , and as specified in wu et al.
<cit> , were normalized according to median expression across all transcripts from a single cell and log <dig> transformed.
the estimated dafs cutoff of log <dig> median-adjusted fpkm data from rna bulk was  <dig> .
at this cutoff all of the known expressed genes in the subset of  <dig> genes were mapped to the he component and all of the known unexpressed genes were mapped to the le component.
the gene closest to the boundary of our cutoff for separating low and high expression was mettl <dig>  with a log <dig> median-adjusted fpkm value of  <dig> .
to assess the expression level of mettl <dig> in hct- <dig> cells, we examined two microarray datasets from the geo database  and confirmed that mettl <dig> is expressed in hct- <dig> cells with low expression levels.
thus, to classify mettl <dig> in the he component appears to be the appropriate decision based on biological evidence.
to demonstrate the robustness of dafs on distinct empirical distributions, four rna-seq samples with unique expression profiles of log2-transformed raw counts were obtained from four different experimental studies.
similar to expression patterns observed in hebenstreit et al.
<cit> , every density exhibited a normally distributed component on the right  and a cluster of the remaining data on the left .
however, each dataset differed in the degree of separation between the le and he component and in the distinct characterization of the le region.
despite variation in the four different patterns of expression, dafs performed consistently well at classifying data as le and he .
in each sample, we were able to select a cutoff that clearly defined the he component.
from a visual perspective, the lower bound of the he component, as identified by our approach, separated the data at a value that preserves as much of the normal he component as possible.
in other words, one would question normality of the he component with a cutoff placed to the left of the estimated optimal quantile cutoff, q
c
, and one would forego valuable normally distributed data with a higher q
c
value.
densities of the two-component mixture distribution obtained by mclust are also presented in figure
<dig>  when the normal he component is well characterized by mclust, the theoretical cutoff and dafs are close.
the gap between estimated cutoffs widens as the model-data misfit increases.
dafs estimates of q
c
for each sample in figure
<dig> demonstrate two additional important findings.
not only does the estimated quantile cutoff differ in each dataset, but the raw rna-seq count corresponding to the computed quantile also differs.
we identified quantile cutoff values of  <dig>  ,  <dig>  ,  <dig>  , and  <dig>  .
these percentiles mapped to raw counts of  <dig>   <dig>   <dig>  and  <dig>  respectively.
this indicates the necessity of a data-adaptive feature.
as evidenced by the results, it would not be appropriate to apply an arbitrary cutoff based on percentiles of the data nor on raw counts without consideration of the individual data structure of the sample.
additionally, we analyzed two different rna-seq datasets with technical replicates to assess the consistency of dafs.
overall, dafs demonstrated great consistency in estimating cutoffs from technical replicates.
for seqc, a q
c
value of  <dig>  was computed for each replicate of sample a. values of q
c
for replicates of sample b fluctuated between  <dig>  and  <dig> .
in hammer et al.
<cit> , mrna-seq in rats was measured 2 weeks and 2 months after l <dig> spinal nerve ligation .
the study included two technical replicates for each treatment condition and also included a control for each time point.
dafs returned q
c
values of  <dig>  and  <dig>  ;  <dig>  and  <dig>  ;  <dig>  and  <dig>  ;  <dig>  and  <dig>  .
overall, dafs showed small variability in estimating quantile cutoffs for technical replicated sequencing data.
these results further demonstrate the impracticality of an arbitrary cutoff even within the same study.
as evidenced by the data, it is reasonable to suspect that the separation of le and he genes is more homogeneous within replicates of the same experimental treatment.
however, we may presume less agreement across biological replicates, treatments, and experiments.
to explore further the spectrum of possible expression profiles, we measured the performance of dafs on sequencing reads from cdna fragments of cultured human b-cells
<cit>  in order to evaluate the functionality of dafs on multi-modal data and data measured at various sequencing depths.
as indicated in toung et al.
<cit> , we pooled all  <dig> unrelated samples obtained from the center d’etude du polymorphisme humain to create an  <dig> million read dataset.
in addition, we randomly selected  <dig> sample and randomly pooled  <dig> and  <dig> samples to simulate sequencing depths of 30– <dig> million reads, 100+ million reads, and 500+ million reads, respectively.
densities of each distribution of log <dig> transformed fpkm data are provided in figure
<dig>  dafs cutoffs for data composed of 30– <dig> million reads, 100+ million reads, 500+ million reads, and  <dig> million reads are  <dig> ,  <dig> ,  <dig> , and  <dig> , respectively.
all are higher than the value of  <dig>  that toung et al.
<cit>  used to classify low expression.
since the authors suggest that  <dig> million reads are necessary to accurately measure transcript expression, our findings indicate that cutoffs estimated by dafs stabilize as the quality of the expression improves.
most importantly, although the log <dig> fpkm values are clearly multi-model, dafs adapts to preserve most of the regions with high expression.
simulation
to evaluate the performance of dafs on simulated log <dig> count data and log <dig> rpkm data, sensitivity, specificity, positive predictive value , and negative predictive value  were evaluated for dafs and mclust .
the mean,  <dig> % percentile, median, and  <dig> % percentile of estimated cutoffs for dafs, the point of intersection between theoretical distributions of the two-component normal mixture fitted by mclust, and for empirical bounds associated with achieving sensitivity values of  <dig>  ,  <dig>  , and  <dig>   are reported in table
<dig>  for better visualization of the results, average cutoff estimates  and 95% confidence intervals for q
c
are provided for both simulation scenarios.
, sen
<dig>
, and sen
<dig>
on simulated data
p
25
p
<dig>
for log <dig> transformed count data, which presented a clear bimodal distribution, the estimated cutoffs from dafs and mclust were similar.
average sensitivity for both methods exceeded 95%.
on average, the theoretical cutoff estimated by mclust was closer to the 95% sensitivity estimate.
although mclust shows an improvement in specificity and dafs does a better job of optimizing the sensitivity and npv, the differences in performance measures were small.
analysis of the 95% confidence intervals for cutoff values demonstrated that both dafs and mclust were not significantly different from estimates obtained from achieving 95% sensitivity.
in evaluating cutoffs for log <dig> transformed rpkm data, dafs returned an average sensitivity estimate that exceeded 98%.
however, mclust estimates fell between values obtained for achieving 90% and 95% sensitivity.
analysis of the 95% confidence intervals for cutoff values presented similar findings.
the results indicated that dafs captured a higher proportion of true high expression.
as such, dafs demonstrated superior performance in sensitivity and npv, whereas ppv was comparable between mclust and dafs.
the low specificity of dafs is explained by the proportion of low expression.
since the low expressed region represented roughly 10% of the total sample size, a small number of misclassified observations will decrease specificity rapidly.
based on the two different scenarios, the results suggest that dafs is comparable to mclust at maximizing sensitivity when the mixture components are well separated; however, dafs is superior at providing a cutoff that maximizes sensitivity when the le and he components are not distinguishable.
the simulation results clearly demonstrate how dafs adapts to different data patterns.
the pattern of log <dig> rpkm data is relatively closer to one normal distribution than the log <dig> raw counts.
therefore, dafs adjusts and retains as many genes as possible when analyzing log <dig> rpkm data.
consequently, relative to the results obtained from simulating log <dig> raw counts, sensitivity and ppv are increased but specificity and npv are decreased.
discussion
this study was motivated by the need to provide a data-adaptive algorithm for separating rna-seq data with low and high expression, particularly when the distributions of expression abundance are not distinctly separated.
in order to compute the optimal cutoff between low and high expression, our method relied heavily on the assumption that high-expressed data are normally distributed
<cit> .
an advantage of using the kolmogorov-smirnov distance as a measure of agreement is that the method is easily modifiable if a different distributional assumption is required to characterize high expression.
for any continuous distribution, the k-s statistic tests the null hypothesis h0 : f = f <dig>  for all x versus the alternative h1 : f ≠ f <dig>  for some x. thus, the reference distribution f <dig>  e.g.
normal, lognormal, student’s t, etc., is completely specifiable by the researcher.
to our advantage, dafs performed consistently well when the he component departed from normality .
careful consideration should be taken with regard to the number of observations used to model the distribution of k-s statistics.
the benefit of adding more data to characterize the underlying density should be balanced with the disadvantage of modeling added noise.
if the predictor space is segmented too finely, then it is possible for multiple percentiles to map to the same k-s statistic.
multiple many-to-one mappings would make it difficult for the mars algorithm to differentiate true variation from random noise.
in the present study, incrementing p by step sizes of  <dig>  or  <dig>  provided a good balance between parsimonious and over-saturated input.
nevertheless, the increment selection is a variable in the proposed methodology that must be investigated by the researcher.
we were not remiss to consider the presence of more than two components of expression abundance.
since the gaussian mixture model can well approximate the shape of any density
<cit> , the number of gaussian mixture components was estimated for multiple datasets.
when mclust was allowed to estimate the number of gaussian mixture components, the algorithm often returned multiple mixture components.
a similar finding was presented in the supplemental material presented of hebenstreit et al.
<cit> .
in their analysis of real datasets, values of aic and bic indicated that the data would be better fit by a k >  <dig> component gaussian mixture.
in our own analysis, many of the identified components were not separated enough to be heterogeneous populations.
we employed a number of methods/packages to merge gaussian mixture components  with no success.
nearly every method struggled by either distinguishing no separation or overly characterizing the distribution of abundance levels.
the latter scenario was more pronounced in le regions, where it seemed apparent that a number of mixture components were necessary to estimate non-gaussian density regions.
it became clear that variability in the expression data made it difficult to ascertain whether homogeneous sub-mixtures could be interpreted as a single component.
in our web search of the literature, the question of a cutoff for low expression in rna-seq was frequently asked.
some questions were motivated by a desire to quantify what is considered expressed in rna-seq.
others were motivated by a need to classify the level of measurement that could be trusted in assessing the significance of differential expression in low-expressed regions, particularly since research shows that the precision of rna-seq data analysis improves as genes are more highly expressed.
whether transcripts with low expression are simply flagged or removed prior to testing via independent filtering, the work presented here provides a data-driven methodology for separating rna-seq expression into meaningful components.
providing an accurate separation of rna-seq data that is not based on ad hoc techniques or methodology that may be prone to model-data misfit will facilitate interpreting the quality of sequencing reads and lead to improved power for differential detection of high expressed, reliable data.
CONCLUSIONS
in this study, we presented a method for classifying transcripts with low and high expression that promises wide-range application.
the robustness of dafs was demonstrated by applying dafs to a number of rna-seq data samples  that varied by sequencing depth, species, normalization, and density shape.

 transcriptome analysis is integral in understanding the role of genetic and genomic variants in disease progression, classification, and diagnosis. in recent years, next generation sequencing technologies have catapulted transcriptome profiling to new heights by providing greater precision of cellular rna content. as a result, rna-seq is more sensitive to subtle changes in expression levels and typically includes many transcripts with low read count integers. researchers have taken interest in the low regions of a sample for several reasons, all with the primary intent of quantifying the level of expression that reflects important biological meaning. one of the major problems in rna-seq data analysis is that there is no consensus on what is considered low expression. low counts have been referenced as less than reads when summed across treatments, less than reads on average, less than reads on average, and less than reads. additionally, count-based pre-filtering methods have been adopted to exclude genes with minimal expression from differential testing. for example, risso et al. filtered out genes with an average count below robinson et al. recommended removing genes that did not have at least counts per million reads in at least two samples in the edger user’s guide. in the deseq user’s guide, anders took a more aggressive approach and suggested removing up to 40% of genes that ranked lowest in regard to total count across all experimental samples. thus, at present, whether low expressed transcripts are simply identified or deliberately omitted from analysis, methods of differentiating spurious rna-seq data from meaningful information have not been explicitly defined. furthermore, it is not clear whether it is justifiable to extrapolate the aforementioned cutoffs to other rna-seq data. a number of factors affect the distribution of read counts for a given study. specifically, the manufacturer, library preparation and construction, alignment algorithm, gene length, sequencing depth, and experimental design all play a role in determining the number of reads that are mapped to a gene. since each of these factors varies from study to study, it may be misleading to ignore properties of the sample distribution by applying an arbitrary cutoff to classify the low expressed region of a sample. to address the question of what should be considered the lower bound of functional gene expression, it is useful to consider the premise that genes can be categorized into a group of high expressed, meaningful genes and a group of low expressed, non-informative genes. in previous studies, the empirical distribution of a sample was used to identify classes of mrna abundance levels. this concept has also been utilized to differentiate functional expression states of microarray data. in order to optimally separate expression classes, a priori knowledge of the expression distribution is useful. in many cases, the precise distribution of noise is unknown; however, characterizations of a global bimodal distribution and a normally distributed he component have been reported in previous studies. hebenstreit et al. studied the global distribution of rna sequences in mice th cells. their work demonstrated that log2-transformed reads per kilobase per million could be separated into two classes of mrna abundance. the researchers used the likelihood ratio test to determine whether the data was best modeled as a mixture of n or n +  components, where < n <  mclust, an r package that uses an expectation-maximization algorithm, was used to compute maximum likelihood estimates of all parameters of the mixture distribution. for one-dimensional data, mclust evaluates a fitted model with equal variance terms and a fitted model with unequal variance and selects the model with higher bayesian information criterion. however, the use of mixture modeling to identify mixture components of rna-seq data may not be appropriate for two reasons. first, it is difficult to capture the true number of mixture components. it is known that the em algorithm performs well at estimating the parameters of a finite mixture model; however, when mixture distributions are unimodal and there is no clear separation between components, the em algorithm commonly returns more components than what seems logical based on visual inspection of the data. second, when a two-component mixture model is forced to fit data that does appear to be bimodally distributed, the fitted model does not always approximate very well the observed empirical distribution. hebenstreit et al. were able to characterize the mixture distribution of the le and he components and determine the peak of the le region by using a poisson distribution to estimate the proportion of undetected genes at each expression level. however, they did not discuss or provide an expression cutoff for separating the two overlapping regions. in this study, we propose dafs, a data-adaptive method for identifying and subsequently flagging expressions in the le region by estimating the lower bound of high expression in a given rna-seq sample. in light of the drawbacks of the mixture modeling approach, dafs was constructed without imposing a finite mixture model on the data. several real rna-seq datasets and simulated data are used to present our findings and demonstrate the robustness of our methodology.  when the le and he regions of rna-seq data are well separated, dafs and mclust separate the components at similar cutoffs. cutoff values determined by dafs and the theoretical intersection of the two mixture components obtained by mclust were - and -, respectively. in their study of th cells, hebenstreit et al. identified a sample of known expressed and unexpressed genes in th cells. based on their classification of low and high expression, they mapped all the known expressed genes to the he component and mapped most of the known unexpressed genes to the le component. we were able to reproduce their findings using the cutoff estimated by dafs. analysis of bulk rna positive controls generated from cultured hct- cells allowed for an additional assessment of the ability of dafs to separate known expressed genes into one mixture component and known unexpressed genes into a separate mixture component. in their analysis of single-cell whole transcriptome amplification, the authors identified known expressed and unexpressed genes in hct- cells. fragments transformed by fragments per kilobase of exon per million mapped reads were downloaded from the geo database, and as specified in wu et al., were normalized according to median expression across all transcripts from a single cell and log transformed. the estimated dafs cutoff of log median-adjusted fpkm data from rna bulk was. at this cutoff all of the known expressed genes in the subset of genes were mapped to the he component and all of the known unexpressed genes were mapped to the le component. the gene closest to the boundary of our cutoff for separating low and high expression was mettl with a log median-adjusted fpkm value of. to assess the expression level of mettl in hct- cells, we examined two microarray datasets from the geo database and confirmed that mettl is expressed in hct- cells with low expression levels. thus, to classify mettl in the he component appears to be the appropriate decision based on biological evidence. to demonstrate the robustness of dafs on distinct empirical distributions, four rna-seq samples with unique expression profiles of log2-transformed raw counts were obtained from four different experimental studies. similar to expression patterns observed in hebenstreit et al., every density exhibited a normally distributed component on the right and a cluster of the remaining data on the left. however, each dataset differed in the degree of separation between the le and he component and in the distinct characterization of the le region. despite variation in the four different patterns of expression, dafs performed consistently well at classifying data as le and he. in each sample, we were able to select a cutoff that clearly defined the he component. from a visual perspective, the lower bound of the he component, as identified by our approach, separated the data at a value that preserves as much of the normal he component as possible. in other words, one would question normality of the he component with a cutoff placed to the left of the estimated optimal quantile cutoff, q c, and one would forego valuable normally distributed data with a higher q c value. densities of the two-component mixture distribution obtained by mclust are also presented in figure  when the normal he component is well characterized by mclust, the theoretical cutoff and dafs are close. the gap between estimated cutoffs widens as the model-data misfit increases. dafs estimates of q c for each sample in figure  demonstrate two additional important findings. not only does the estimated quantile cutoff differ in each dataset, but the raw rna-seq count corresponding to the computed quantile also differs. we identified quantile cutoff values of,,, and. these percentiles mapped to raw counts of and respectively. this indicates the necessity of a data-adaptive feature. as evidenced by the results, it would not be appropriate to apply an arbitrary cutoff based on percentiles of the data nor on raw counts without consideration of the individual data structure of the sample. additionally, we analyzed two different rna-seq datasets with technical replicates to assess the consistency of dafs. overall, dafs demonstrated great consistency in estimating cutoffs from technical replicates. for seqc, a q c value of was computed for each replicate of sample a. values of q c for replicates of sample b fluctuated between and. in hammer et al., mrna-seq in rats was measured 2 weeks and 2 months after l spinal nerve ligation. the study included two technical replicates for each treatment condition and also included a control for each time point. dafs returned q c values of and; and; and; and. overall, dafs showed small variability in estimating quantile cutoffs for technical replicated sequencing data. these results further demonstrate the impracticality of an arbitrary cutoff even within the same study. as evidenced by the data, it is reasonable to suspect that the separation of le and he genes is more homogeneous within replicates of the same experimental treatment. however, we may presume less agreement across biological replicates, treatments, and experiments. to explore further the spectrum of possible expression profiles, we measured the performance of dafs on sequencing reads from cdna fragments of cultured human b-cells in order to evaluate the functionality of dafs on multi-modal data and data measured at various sequencing depths. as indicated in toung et al., we pooled all unrelated samples obtained from the center d’etude du polymorphisme humain to create an million read dataset. in addition, we randomly selected sample and randomly pooled and samples to simulate sequencing depths of 30– million reads, 100+ million reads, and 500+ million reads, respectively. densities of each distribution of log transformed fpkm data are provided in figure  dafs cutoffs for data composed of 30– million reads, 100+ million reads, 500+ million reads, and million reads are,,, and, respectively. all are higher than the value of that toung et al. used to classify low expression. since the authors suggest that million reads are necessary to accurately measure transcript expression, our findings indicate that cutoffs estimated by dafs stabilize as the quality of the expression improves. most importantly, although the log fpkm values are clearly multi-model, dafs adapts to preserve most of the regions with high expression. simulation to evaluate the performance of dafs on simulated log count data and log rpkm data, sensitivity, specificity, positive predictive value, and negative predictive value were evaluated for dafs and mclust. the mean, % percentile, median, and % percentile of estimated cutoffs for dafs, the point of intersection between theoretical distributions of the two-component normal mixture fitted by mclust, and for empirical bounds associated with achieving sensitivity values of,, and are reported in table  for better visualization of the results, average cutoff estimates and 95% confidence intervals for q c are provided for both simulation scenarios., sen, and sen on simulated data p 25 p for log transformed count data, which presented a clear bimodal distribution, the estimated cutoffs from dafs and mclust were similar. average sensitivity for both methods exceeded 95%. on average, the theoretical cutoff estimated by mclust was closer to the 95% sensitivity estimate. although mclust shows an improvement in specificity and dafs does a better job of optimizing the sensitivity and npv, the differences in performance measures were small. analysis of the 95% confidence intervals for cutoff values demonstrated that both dafs and mclust were not significantly different from estimates obtained from achieving 95% sensitivity. in evaluating cutoffs for log transformed rpkm data, dafs returned an average sensitivity estimate that exceeded 98%. however, mclust estimates fell between values obtained for achieving 90% and 95% sensitivity. analysis of the 95% confidence intervals for cutoff values presented similar findings. the results indicated that dafs captured a higher proportion of true high expression. as such, dafs demonstrated superior performance in sensitivity and npv, whereas ppv was comparable between mclust and dafs. the low specificity of dafs is explained by the proportion of low expression. since the low expressed region represented roughly 10% of the total sample size, a small number of misclassified observations will decrease specificity rapidly. based on the two different scenarios, the results suggest that dafs is comparable to mclust at maximizing sensitivity when the mixture components are well separated; however, dafs is superior at providing a cutoff that maximizes sensitivity when the le and he components are not distinguishable. the simulation results clearly demonstrate how dafs adapts to different data patterns. the pattern of log rpkm data is relatively closer to one normal distribution than the log raw counts. therefore, dafs adjusts and retains as many genes as possible when analyzing log rpkm data. consequently, relative to the results obtained from simulating log raw counts, sensitivity and ppv are increased but specificity and npv are decreased. discussion this study was motivated by the need to provide a data-adaptive algorithm for separating rna-seq data with low and high expression, particularly when the distributions of expression abundance are not distinctly separated. in order to compute the optimal cutoff between low and high expression, our method relied heavily on the assumption that high-expressed data are normally distributed. an advantage of using the kolmogorov-smirnov distance as a measure of agreement is that the method is easily modifiable if a different distributional assumption is required to characterize high expression. for any continuous distribution, the k-s statistic tests the null hypothesis h0: f = f for all x versus the alternative h1: f ≠ f for some x. thus, the reference distribution f e.g. normal, lognormal, student’s t, etc., is completely specifiable by the researcher. to our advantage, dafs performed consistently well when the he component departed from normality. careful consideration should be taken with regard to the number of observations used to model the distribution of k-s statistics. the benefit of adding more data to characterize the underlying density should be balanced with the disadvantage of modeling added noise. if the predictor space is segmented too finely, then it is possible for multiple percentiles to map to the same k-s statistic. multiple many-to-one mappings would make it difficult for the mars algorithm to differentiate true variation from random noise. in the present study, incrementing p by step sizes of or provided a good balance between parsimonious and over-saturated input. nevertheless, the increment selection is a variable in the proposed methodology that must be investigated by the researcher. we were not remiss to consider the presence of more than two components of expression abundance. since the gaussian mixture model can well approximate the shape of any density, the number of gaussian mixture components was estimated for multiple datasets. when mclust was allowed to estimate the number of gaussian mixture components, the algorithm often returned multiple mixture components. a similar finding was presented in the supplemental material presented of hebenstreit et al.. in their analysis of real datasets, values of aic and bic indicated that the data would be better fit by a k >  component gaussian mixture. in our own analysis, many of the identified components were not separated enough to be heterogeneous populations. we employed a number of methods/packages to merge gaussian mixture components with no success. nearly every method struggled by either distinguishing no separation or overly characterizing the distribution of abundance levels. the latter scenario was more pronounced in le regions, where it seemed apparent that a number of mixture components were necessary to estimate non-gaussian density regions. it became clear that variability in the expression data made it difficult to ascertain whether homogeneous sub-mixtures could be interpreted as a single component. in our web search of the literature, the question of a cutoff for low expression in rna-seq was frequently asked. some questions were motivated by a desire to quantify what is considered expressed in rna-seq. others were motivated by a need to classify the level of measurement that could be trusted in assessing the significance of differential expression in low-expressed regions, particularly since research shows that the precision of rna-seq data analysis improves as genes are more highly expressed. whether transcripts with low expression are simply flagged or removed prior to testing via independent filtering, the work presented here provides a data-driven methodology for separating rna-seq expression into meaningful components. providing an accurate separation of rna-seq data that is not based on ad hoc techniques or methodology that may be prone to model-data misfit will facilitate interpreting the quality of sequencing reads and lead to improved power for differential detection of high expressed, reliable data.  in this study, we presented a method for classifying transcripts with low and high expression that promises wide-range application. the robustness of dafs was demonstrated by applying dafs to a number of rna-seq data samples that varied by sequencing depth, species, normalization, and density shape. 
due to the large dynamic range of expression levels, rna-seq is more prone to detect transcripts with low expression. to date, there is no consensus on the definition of low expression. since random variation is high in regions with low expression and distributions of transcript expression are affected by numerous experimental factors, methods to differentiate low and high expressed data in a sample are critical to interpreting classes of abundance levels in rna-seq data. the kolmgorov-smirnov statistic and multivariate adaptive regression splines were used to determine the optimal cutoff value for separating transcripts with high and low expression. the robustness of the proposed method was demonstrated by analyzing different rna-seq datasets that varied by sequencing depth, species, scale of measurement, and empirical density shape. a data-adaptive approach was developed to estimate the lower bound of high expression for rna-seq data. the analysis of real and simulated data presented here illustrates the need to employ data-adaptive methodology in lieu of arbitrary cutoffs to distinguish low expressed rna-seq data from high expression. the ability to ascertain stably expressed rna-seq data is essential in the filtering process of data analysis, and methodologies that consider the underlying data structure demonstrate superior performance in preserving most of the interpretable and meaningful data. results from the proposed method were compared to results obtained by estimating the theoretical cutoff of a fitted two-component mixture distribution. our results also present the drawbacks of characterizing the data by a two-component mixture distribution when classes of gene expression are not well separated. the proposed algorithm for classifying low and high regions of transcript abundance promises wide-range application in the continuing development of rna-seq analysis.
 next-generation sequencing has advanced the application of high-throughput sequencing technologies in genetic and genomic variation analysis. due to the large dynamic range of expression levels, rna-seq is more prone to detect transcripts with low expression. it is clear that genes with no mapped reads are not expressed; however, there is ongoing debate about the level of abundance that constitutes biologically meaningful expression. to date, there is no consensus on the definition of low expression. since random variation is high in regions with low expression and distributions of transcript expression are affected by numerous experimental factors, methods to differentiate low and high expressed data in a sample are critical to interpreting classes of abundance levels in rna-seq data.  a data-adaptive approach was developed to estimate the lower bound of high expression for rna-seq data. the kolmgorov-smirnov statistic and multivariate adaptive regression splines were used to determine the optimal cutoff value for separating transcripts with high and low expression. results from the proposed method were compared to results obtained by estimating the theoretical cutoff of a fitted two-component mixture distribution. the robustness of the proposed method was demonstrated by analyzing different rna-seq datasets that varied by sequencing depth, species, scale of measurement, and empirical density shape.  the analysis of real and simulated data presented here illustrates the need to employ data-adaptive methodology in lieu of arbitrary cutoffs to distinguish low expressed rna-seq data from high expression. our results also present the drawbacks of characterizing the data by a two-component mixture distribution when classes of gene expression are not well separated. the ability to ascertain stably expressed rna-seq data is essential in the filtering process of data analysis, and methodologies that consider the underlying data structure demonstrate superior performance in preserving most of the interpretable and meaningful data. the proposed algorithm for classifying low and high regions of transcript abundance promises wide-range application in the continuing development of rna-seq analysis. rna-sequencinglow expressiondata-adaptiveflagmixture distribution 