BACKGROUND
the possibility of genomic phenotype classification arose with the inception of gene-expression microarrays.
from the outset, two fundamental problems have frustrated the endeavor:  the inaccuracy of microarray measurements, and  small samples.
our particular application of interest is classification using rna-seq data.
modern rna-seq technologies sequence small rna fragments  to measure gene expression, where the number of reads mapped to a gene on the reference genome defines the count data.
given that rna-seq data has advantages over microarray data, in particular, more accurate measurement, we still confront the second fundamental problem, which is statistical, not technological: small samples cause re-sampling-based classifier error estimators to be very inaccurate due to excessive variance and lack of regression with the true error  <cit> .
since the error rate of a classifier quantifies its predictive accuracy, it is the salient epistemological attribute of any classifier.
the inability to satisfactorily estimate the error with model-free methods with small samples implies that genomic classifier error estimation is virtually impossible without the use of prior information, so that the whole small-sample classification problem becomes unapproachable in a model-free framework  <cit> .
the situation has been addressed by utilizing prior knowledge via a bayesian approach that considers a prior distribution on an uncertainty class of feature-label distributions  <cit> .
for expression-based classification, prior distributions have been constructed using expression data not employed in classifier design  <cit>  and known regulatory pathways  <cit> .
given that a prior model must be assumed to achieve satisfactory error estimation, an obvious course of action is to derive an optimal classifier based on the prior knowledge and the sample data, the result being an optimal bayesian classifier  that is guaranteed to have the best average performance of any classifier relative to the posterior distribution derived from the prior distribution and data  <cit> .
while bayesian classification does not depend on particular distributional forms, closed-form solutions have been derived for the multinomial model and gaussian models using linear classifiers for the minimum mean squared error  error estimate  <cit> , the mse of the error estimate  <cit> , and an optimal bayesian classifier  relative to the prior distribution  <cit> , the latter being expressed in terms of effective class conditional distributions, which are expectations relative to the posterior distribution of the class-conditional distributions.
the closed-form solutions depend on particular models  and the existence of conjugate priors, which can be too constraining for practical applications such as rna-seq classification.
much of the statistical literature concerning classification of rna-seq data attempts to address differential expression testing, that is, univariate statistical testing on an individual gene basis.
these attempts typically model rna-seq data via negative binomial  <cit>  and poisson distributions  <cit> .
in addition, network inference has been attempted using a hierarchical poisson log-normal model  <cit> , and clustering of rna-seq data points has utilized various approaches  <cit> .
however, in clinical settings one is often interested in sample classification: the problem of classifying the rna-seq data from unlabeled patients using a set of labeled training data.
one of the few rna-seq-specific attempts towards this goal uses a poisson modeling assumption with independent features  <cit> .
the poisson model is completely parameterized by its mean and thus is known to exhibit problems in fitting rna-seq data due to the overdispersion typically observed in such datasets.
in this paper, we focus on modeling the pipeline that starts with extracting the gene concentrations from the biological samples and their subsequent processing by the sequencing instrument  <cit> .
this is accomplished using a hierarchical, multivariate poisson model .
specifically, gene concentration levels are modeled by a log-normal distribution and the sequencing instrument sampling of those is modeled via a poisson process.
this allows us to accurately model the rna-seq data overdispersion as demonstrated by marginal variance calculations and posterior predictive model diagnostics in section ‘overdispersion’.
in addition, this hierarchical model allows for inferring any covariance structure observed between the features.
whereas dalton and dougherty have presented a computational method for nonlinear classifiers in the gaussian model  <cit> , this still depends upon conjugate priors.
in this work, we remove the constraints imposed by the requirement of a closed-form solution by developing the optimal bayesian classifier using a markov-chain-monte-carlo  methodology.
this provides a computational framework for calculating the obc for any parameterized class conditional-density and any prior distribution.
most notably, this allows us to use distributions designed to closely model particular datasets and a prior distribution of any form to improve classification performance in small-sample settings, in particular, for rna-seq-based classification.
methods
notation
throughout, we use capital letters to indicate random variables, lower case letters to indicate individual realizations of random variables or indices, bold latin characters for observed vectors, and greek letters for latent features and parameters.
we write p as the probability measure over the random variable x. p may be a probability mass function, probability density function, or arbitrary probability measure.
p denotes the conditional probability p. similarly, following bayesian convention, we write parameterized distributions by conditioning on the parameter, for instance, p, and posterior expectations by conditioning on the sample, such as e, where sn and all other values are defined in section ‘review of optimal bayesian classification’.
if it is unclear which density an expectation is taken with respect to, then we denote it in subscript notation, such as eθ|sn·, where the expectation is taken with respect to the density p.
review of optimal bayesian classification
binary classification considers a set of n labeled training data points, sn=1n, where y∈{ <dig> } is the class label and x∈x is the feature vector over a feature space .
an example of binary classification in a clinical setting might include class  <dig> and  <dig> being two types of cancers, or normal and cancerous tissues.
available features would then be the gene or genes that will eventually be used in the designed classifier to assign this label.
the feature space  would be the set of possible gene expression measurements for all genes in the feature vector.
the labeled training data sn would be the set of gene expression measurements from samples which had undergone further testing  to identify the type or malignancy of the tissue.
using sn, we design a classifier ψ that hopefully performs well on the unknown joint feature-label distribution p. in the same clinical example, the classifier ψ could then identify the type of cancer using gene expression measurements alone.
by parameterizing this unknown joint distribution in a model-based bayesian framework one can derive an optimal bayesian classifier  that minimizes the expected error over the space of all classifiers under assumed forms of the class-conditional densities.
specifically, under gaussian and multinomial class-conditional densities and their corresponding conjugate prior distributions, closed-form solutions for the obc  <cit>  and the first two moments of the error estimate conditioned on the sample  <cit>  have been obtained.
the parameterization of the feature-label distribution consists of the marginal class probability c and the class-conditional densities p, where a particular value θy∈Θy specifies a single class-conditional density contained in the class of densities defined over the space Θy, which will be a cartesian product as described in section ‘the multivariate poisson model’.
therefore, for a two-class problem, we specify a parameterized joint feature-label distribution as θ=∈Θ= <cit> ×Θ0×Θ <dig>  in the bayesian classification framework, these values are then treated as random variables, so that we may consider quantities such as the expectation of c, or another random variable conditioned on the value of the parameter vector θ.
figure  <dig>
bayesian classification derivation tree.
a tree summarizing the relationships between several important quantities in the general theoretical framework of bayesian classification.
a directed edge between a parent and its child indicates that the child can be derived from the parent by the equations indicated in the edge label.
the root of the tree p is the posterior distribution of the feature label parameters and by taking expectations with respect to this distribution, we can derive the effective class conditional densities p and the distribution of the classifier error p. then these quantities give rise to the obc, and mmse and mse estimates for the error as described in the text.
quantities highlighted in grey are given in closed form for gaussian and multinomial distributions in  <cit> .
we now examine the tree in more detail.
starting at the far left of the tree, p is the posterior distribution of the parameterized feature-label distribution – posterior to the labeled samples in sn.
typically, error estimates and the optimal classifier are our primary interest, so that this posterior distribution is traditionally used as a means to compute other quantities and is not of interest by itself.
the effective class-conditional density is the marginal predictive posterior of the feature vector x conditioned sn and the class variable y,
px|y,sn=∫Θypx|y,θypθy|sndθy.
it gives the distribution of the feature vector using a weighted average over all the parameterized class-conditional densities in Θy given a class y. the weights in this expectation are the posterior, p, evaluated at each θy.
the true error of classifier ψ is ε=p≠y).
given the sample data sn, ε is a random unknown quantity in the bayesian framework.
the mmse estimate given in  <cit>  can be written as
eε|sn=pψ≠y|sn=eθ|snpψ≠y|θ,sn=ĉε0θ <dig> ψ+ε1θ <dig> ψ=∫xĉpx| <dig> snix∈r1+px| <dig> snix∈r0dx,
where ia is the indicator function for event a, ĉ=ec|sn is the posterior expectation of c, ry is the region of the feature space the classifier predicts to be class y,  is the feature space, and εy is the error of classifier ψ contributed by class y on the fixed distribution θy.
we can also obtain the full posterior distribution of the error,
pε|sn=∫Θpε|θpθ|sndθ=eθ|snpε|θ,
where p is the true error for a fixed feature-label distribution and fixed classifier.
we denote this deterministic function by ε. as shown in figure  <dig>  the mmse estimate and the sample conditioned mse for this error can also be calculated using the first two moments of the error distribution.
with the mmse estimator defined, the optimal bayesian classifier  is the classifier minimizing the expected error by pointwise minimization of the integral   <cit> :
ψobc=0ifĉp≥p,1otherwise..
conditional error estimator
if the true feature-label distribution were known, then we could compute the true error of a classifier exactly as an expectation over the conditional error  <cit> :
ε=p≠y)=∫xp≠y|x)pdx.
treating ε as a random variable, one can similarly derive its posterior distribution by conditioning on the feature vector:
pε|sn=∫xpε,x|sndx=∫x∫Θpε,θ,x|sndθdx=∫Θpθ|sn∫xpε|x,θpx|sndxdθ,
which is different than the derivation of the same quantity in .
this introduces the idea of the conditional error estimator, which we define as the mmse estimate of the classification error conditioned on the feature vector x,
ε^=eθ|snε|x,sn=pψ≠y|x,sn=px|y≠ψ,snpy≠ψ|snpx|sn=z−1px|y≠ψ,snpy≠ψ|sn,
as expanded through application of bayes’ theorem, where z is a normalizing constant given by
z=px|sn=∑y∈ <dig> px|y,snpy|sn.
in addition to being useful in the above alternative derivation of the classifier’s error posterior, the conditional error estimate has other practical applications.
when classifying an unlabeled data point, we would like to estimate the error of the classifier output for that particular data point, as opposed to the overall error estimate for the classifier.
for the obc, from  the conditional error estimator can be written as
ε^ψobc,x=z−1miny∈ <dig> px|y,snpy|sn.
in sum, using the effective class-conditional densities and the posterior marginal probabilities one can calculate conditional error estimates for points in the feature space in addition to the earlier quantities described.
the multivariate poisson model
with the widespread use of next-generation sequencing techniques, classification approaches must be developed to account for the discrete nature of the mapped sequence data and to accommodate the various types of prior information available regarding these experiments.
gene concentration levels can be modeled using a log-normal distribution  <cit> .
as discussed in the introduction, we assume that the sequencing instrument samples this mrna concentration through a poisson process and obtains xi,j reads for sample point i and gene j. we model this as
pxi,j|λi,j∼poissondiexpλi,j,
where λi,j is the location parameter of the log-normal distribution for sample i and gene j, and di is a variable accounting for the sequencing depth as determined by the sequencing process  <cit> .
for each i, we model the location parameter vector λi with a multivariate gaussian distribution, λi∼normal.
we then consider the mean μ and covariance Σ of the gene concentrations as independent quantities for each class y.
the entire mp model is represented in figure  <dig> as a plate diagram.
the distribution of a single class y is parameterized by θy=, where d= and λ=,i= <dig> ,…n,j= <dig> ,…,d, for n sample points and d total genes.
therefore, θy∈Θy=ℝd×ℝd×d×ℝn×ℝd×n.
the feature-label distribution parameterization for the two-class problem is then given by θ=, where c=p, the prior probability for class  <dig> multivariate poisson model plate diagram.
a plate diagram for the multivariate poisson model.
the outermost plate represents the classes that we are interested in classifying against, where i is the index of the sample in class y, and j are the genes being modeled.
to ensure a proper posterior with unit integral, we place weakly informative priors over the latent variables in the mp model.
in choosing these values, we have aimed to avoid the complications that can occur with overly diffuse priors, such as lindley’s paradox  <cit> .
we choose:
μy∼normalηy,ν2idΣy∼inverse-wishartκy,syc∼beta,  where each element of μy is distributed according to a univariate gaussian.
unless otherwise stated, η is the d dimensional zero vector, ν2= <dig>  κ= <dig>  and s=id.
for computational and identifiability reasons, d is fixed to be a vector of normalization constants in order to match the different sequencing depths across all the samples.
in practice, d can be approximated by an upper quartile normalization, which has been shown to be effective  <cit> .
in any bayesian approach the choice of prior affects the results, especially when only a few data points are given.
in the case of mmse classifier error estimation in the bayesian framework, robustness to incorrect modeling assumptions has been extensively studied in  <cit>  and in those studies performance held up well for various kinds of incorrect modeling assumptions.
robustness of optimal bayesian classifiers to false modeling assumptions was extensively studied in  <cit> .
again, good robustness was exhibited.
of course, one can get bad small-sample results by intentionally selecting an inaccurate prior.
in general, if one is confident in his knowledge, then a tight prior is called for because tighter priors require less data for good performance; on the other hand, when one is not confident, then prudence calls for a less informative prior.
as proven in  <cit> , obc classification is consistent under very general conditions; however, a prior whose mass is concentrated far away from the true parameters will perform worse than one that is non-informative.
these issues have been extensively discussed in the bayesian literature  <cit> .
in the end, performance is the measure of worth and our results with synthetic and real data indicate solid performance for the modeling approach used herein.
overdispersion
the mp model uses the poisson distribution in a hierarchical scheme.
it is important to note that, while the read counts are modeled as conditionally poisson in equation  <dig>  the observed read counts are not marginally poisson distributed.
to demonstrate this, consider a one-dimensional simplification of the mp model in which x is the number of reads observed, λ is the log of the rna concentration, and
λ∼normalμ,σ2x∼poisson).
then for the marginal variance of x,
var=evar+var=e+e≥eμ=varpoissoneμ  where μ and σ <dig> are the mean and variance of the log of the concentration.
therefore, when σ2> <dig>  the marginal variance of x is always greater than that of a poisson random variable with the same effective rate.
in addition, by carrying out a posterior predictive model check  <cit> , p.  <dig>  by computing marginal posterior p-values against real rna-seq data, we can quantitatively assess the ability of the mp model to fit the dispersion of the tcga data.
for a test statistic t, we compute the p-value by comparing the test statistic on the true data t and the value of the statistic averaged across the posterior predictive distribution t, where xrep∼p:
pt=pr≥t|sn)=pr>t|sn)+pr=t|sn)≈1m∑i=0mi{t)>t}+ <dig> i{t)=t},  where xrep are monte carlo samples taken from the posterior predictive distribution p using the m monte carlo samples from the posterior distribution of θ as described in section ‘computation’.
the term pr=t|sn) is necessary due to the discrete nature of rna-seq data.
p-values away from  <dig> and  <dig> indicate that the model posterior produces test statistics both above and below that measured on the real data.
we also consider where the real test statistic falls in relation to credible intervals of the test statistic to consider the magnitude of any differences.
we apply the inter-quartile distance test statistic to provide a measure of the mp model’s ability to fit the dispersion of rna-seq data.
we also consider several other test quantities in the additional file 1: table s1-s <dig>
prior calibration using discarded features
since designed classifiers typically use very few of the totality of observed genes, only a small fraction of the data is used for classifier design.
similarly to  <cit> , we can use the discarded features to calibrate the inverse-wishart prior for our mp obc.
our goal is to obtain hyperparameters s,m,κ, and ν <dig> for each class from our training data sn.
in general, we do not expect the discarded features to give us information about any particular genes and the specific covariances between genes, so we make the simplifying assumptions that we learn information from the discarded genes in an aggregate sense.
thus, we consider the following structure on the hyperparameters: m=mt and
where m∈ℝ, σ2> <dig>  and −1≤ρ≤ <dig>  for each class, we need to determine values for five scalar quantities: m,ν <dig> σ <dig> ρ, and κ.
due to the hierarchical design of the mp model, we cannot apply the method of moments in a direct fashion, as did  <cit> .
instead, we utilize a sampling based approach to the method of moments.
this mcmc sampling approach has been examined in  <cit>  as an extension to the generalized method of moments  <cit> .
the sampling approach uses the discarded features in an additional mcmc run evaluated prior to the primary classification mcmc procedure as discussed in section ‘computation’ – and then proceeds to the method of moments.
in this calibration mcmc, we initialize all prior distributions with flat priors and use the discarded features to obtain samples from the posterior distribution of μ and Σ. typically, the number of discarded features f is much larger than the dimensionality d of the classification problem.
therefore, due to computation time, we uniformly sample fs pairs of features from f and average the resulting runs rather than using all or large groups of discarded features in a single mcmc run.
we use the following procedure :
for each randomly chosen discarded feature pair :
obtain mcmc samples using the feature pair as data and flat priors.
record posterior averages of μ and Σ.
average over these posterior averages as given by equations - .
using the resulting five hyperparameter estimates, run the final mcmc for classification.
following  <cit> , we use the moments of the posterior samples to determine the hyperparameters through the following relations: the mean of an inverse-wishart distribution is
eΣ=sκ−d− <dig>
which together with our simplified covariance structure implies
σ2=κ−d−1eΣ <dig>
ρ=eΣ12eΣ <dig>
the variance of the first diagonal of an inverse-wishart matrix can be used to solve for κ via
κ=2eΣ112varΣ11+d+ <dig>
as we have samples of μ directly from our posterior, we obtain
m=eμ <dig>
ν=varμ <dig>
in order to use equations -, we obtain estimates of the moments from mcmc performed over the fs discarded feature pairs.
for the i-th feature pair we obtain the posterior means μ^ <dig> Σ^ <dig>  and Σ^ <dig> and then average:
e^μ1=1fs∑i=1fsμ^ <dig>
var^μ1=1fs−1∑i=1fse^μ1−μ^ <dig>
e^Σ11=1fs∑i=1fsΣ^11+Σ^ <dig>
e^Σ12=1fs∑i=1fsΣ^ <dig>
var^Σ11=1fs−1∑i=1fse^Σ11−Σ^ <dig>
we substitute the estimates from equations - back into equations - to obtain the final hyperparameter estimates.
one must keep in mind that the calibration procedure explicitly assumes the mp model.
hence, one can only expect an improvement in the classification performance if the data follow the mp model.
computation
to obtain the mp obc, we approximate the effective class conditional densities in order to minimize the expected error in a pointwise fashion:
px|y,sn=∫Θypx|y,θypθy|sndθy≈1m∑i=1mpx|y,θy,
where θy are m samples of θy from the model posterior distributions.
for clarity of presentation, we do not consider the class variable y, and we assume a single class.
we do this because the computation can be performed per-class due to the assumed independence between the classes and the marginal probability, p=ppp.
to obtain posterior samples of θ using the metropolis hastings mcmc algorithm we define a proposal distribution p to obtain a new value for the class parameters θ′ from the old values θ. we then calculate the acceptance ratio
r=min <dig> pθ′|snpθ′|θpθ|snpθ|θ′=min <dig> psn|θ′pθ′psn|θpθ,  under the assumption of a symmetric proposal distribution =p).
the process of proposing and accepting samples from this distribution with the probability r induces a markov chain.
positivity of the proposal distribution > <dig> for any θ) is a sufficient condition for ergodicity of this markov chain.
furthermore, this markov chain admits a steady-state distribution equal to our desired posterior distribution p  <cit> .
from the definition of the likelihood,
psn|θ=∏i=1npxi|θ=∏i=1npxi|λi=∏i=1n∏d=1dpxi,d|λi,d,  where p=p owing to conditional independence.
from the definition of the prior,
p=p=ppp=∏i=1nppp.
the posterior predictive distribution in  is approximated by
px|sn≈1m∑i=1mpx|θ=1m∑i=1m∫Λpλ|θpx|λdλ=1m∑i=1m∫Λpλ|θ∏k=1dpxk|λkdλ≈1mt∑i=1m∑g=1t∏k=1dpxk|λk,  where, p∼poisson), λ∼ normal , Λ=ℝn×d, and the λ are t vector-valued samples drawn from the appropriate class’s posterior distribution used to approximate the inner intractable integral.
in addition, we use this approximation of the effective class-conditional density to calculate the conditional error estimates of  in a pointwise fashion.
finally, because we have assumed a conjugate prior distribution for the marginal class probability c, the posterior expectation takes the closed form
eθ|snc=n0+α0n0+n1+α0+α <dig>   where the ny are the number of training samples obtained from class y and the αy are hyperparameters set to  <dig> for an uninformative prior.
conjugacy was used for this one parameter because the increased flexibility of the full sampling approach was deemed not necessary due to the constrained, univariate nature of the parameter.
if more complex relationships between c and other parameters were desired, then a sampling approach using non-conjugate priors would be straightforward to implement.
synthetic data
to evaluate obc performance in the setting of the mp model, we generate synthetic data using the method proposed in  <cit>  to simulate gene expression/mrna concentrations .
these gene expression values are then statistically sampled to emulate modern sequencing machines as described in  <cit> .
parameter values are drawn from the following distributions to examine a wide variety of classification problems:
μy∼normal,σy∼inverse-gamma,ρ=uniform,dlow= <dig> dhigh= <dig> blocksize= <dig>
with these parameters, ten global, twenty heterogeneous, and ten non-marker features are generated.
then four features are randomly chosen to represent a mixture of features of various classification quality.
following  <cit> , the features in the data are zero mean and unit standard deviation normalized except for the mp obc.
the exception occurs because the mp model expects features to be positive integers and normalization is not necessary.
the discarded features are used for calibration of the mp obc priors, and  <dig> samples are generated from each class to estimate the true classification rate for each classifier.
we use four features in this synthetic data classification study owing to limited computational resources as discussed in section ‘computational limitations’.
the synthetic data generation method proposed in  <cit>  imposes the strong assumption of a homogeneous covariance  structure between the two classes of data.
this assumption does not hold for biological situations where interactions between features are not necessarily preserved between classes, and this occurs frequently in biology when considering the possible effects of canalizing genes, nonlinear gene regulation, and mutations in the case of cancer  <cit> .
specifically, if the canalizing gene is not observed, and differs in activity between the two classes, then the measured correlation between two downstream genes could potentially be negligible in one class while strong in the other class.
similarly, for highly nonlinear gene regulation, if a gene in one class is in the saturation region of its response curve from a master gene, then the correlation will be low, while a lower expression level in the other class would allow for a large measured correlation with the same canalizing gene.
and finally, if one class represents normal gene expression and the other tumor-related expression, then a correlation might exist from a functioning pathway in the normal tissue, but a mutation could result in a lack of correlation effects in the tumor.
hence, we modify the synthetic data generation procedure in an attempt to produce synthetic datasets more representative of such nonlinear phenomena in biology.
in this modified procedure, we allow independent covariance  matrices for the features of the two classes.
to generate these covariance matrices, Σy, we utilize independent draws from inverse-wishart distributions with parameters κy= <dig> d= <dig>  and scale matrix s=σy <dig> id.
to examine the effects of feature correlation in ic datasets, we can also generate low-correlation covariance matrices by zeroing the off-diagonal terms.
once the covariance matrix for class y is obtained, location parameters for gene-expression values for each sample point are drawn from the respective multivariate normal distribution λy∼n.
each sample point is then assumed to be normalized through an upper quartile or other suitable method, but in practice any sample-based normalization is imperfect.
we reflect this variation by drawing the sequencing depth di from a uniform distribution, giving the rate of the poisson process as di expλi.
the number of reads for a single gene from a single sample is then drawn from this poisson distribution.
see additional file  <dig> for more detail.
the obc is optimal on average across the space of distributions determined by its prior distributions.
to avoid biasing the performance comparison, we draw the classification problem datasets using different distributions than those of the obc priors.
see additional file  <dig> for more detail.
real data
we consider a real rna-seq dataset composed of level  <dig>  rnaseqv <dig> data from the cancer genome atlas  project.
it contains  <dig> and  <dig> specimens from lung adenocarcinoma and lung squamous cell carcinoma tumor biopsies, respectively.
the samples are mapped read counts against  <dig> known human rna transcripts as generated by the university of north carolina at chapel hill, one of the genome sequencing centers for the tcga.
the data for each cancer type is the result of processing approximately  <dig> billion reads and the read count files for each are one gigabyte apiece.
the problem is to classify the tumor types.
because the class- <dig>  and class- <dig>  sample sizes,  <dig> and  <dig>  are not chosen randomly, we are confronted with the problem of separate sampling.
this means that there is no way to obtain a posterior distribution for c and therefore c must be known in advance.
based upon records from 2006- <dig>  we have a very accurate estimate,  <dig> / <dig> ≈ <dig>   <cit> .
whereas we can use the value of c directly, along with all of the data, in designing the obc, for classification rules that do not use c explicitly, the separately sampled data must be maximally subsampled to the proper sampling ratio c before applying the classification rule  <cit> .
this means that for ntrn desired samples, the sample subsets will contain roundntrn) and round for class  <dig> and  <dig>  respectively.
moreover, holdout error estimation, which we use here, must be properly adapted for separate sampling for all design methods, including the obc.
the holdout estimate is given by
ε^c=cε^0+ε^ <dig>   where ε^ <dig> and ε^ <dig> are the ordinary holdout estimators  for the class- <dig> and class- <dig> errors, respectively  <cit> .
we note that many studies have made the mistake of using classification rules designed for random sampling when sampling is separate.
this can have devastating effects on classifier performance  <cit> .
while averaging over sample subsets for holdout error estimation, we also average over uniformly, randomly selected gene subsets of size  <dig>  this sampling occurs from low  expression genes.
we sample from these lower expression genes because we are ultimately interested in classification problems where the delineation between phenotypes is determined by genes with low expression.
we used  <dig>  for averaging in order to obtain a large enough sample over this feature and sample subset space to achieve repeatable results .
computational runtime for each sample and gene subset was similar to the synthetic data.
RESULTS
the additional file  <dig> contains a simple two-class, two-feature demonstration of the overall procedure to allow for easy visualization and interpretation.
here we discuss the results for the synthetic and real data.
synthetic data
to evaluate classification performance, classifiers were trained using 3nn, lda, and c-support vector machines with a radial basis function kernel  <cit> .
starting with the homogeneous-covariance model, figure 3a shows that the performance of the multivariate poisson obc is better than nonlinear svm when more than  <dig> samples are available and is significantly better than any other classifier when using calibrated features.
equivalently, by using discarded features, we can obtain the same classification performance while requiring fewer training samples.synthetic rna-seq classification.
synthetic data classification results with  homogeneous-covariance,  high correlation independent-covariance,  low correlation independent-covariance, and  high correlation independent-covariance data at several problem difficulties.
in the case of independent-covariance data with highly correlated features, figure 3b shows superior classification performance of the mp obc at nearly all sample sizes considered.
in addition, for calibrated prior distributions, the performance of the mp obc improves.
this improvement is greater when the sample sizes are small, which demonstrates the importance of additional knowledge  when data are expensive to obtain or not readily available.
the superior performance of the obc relative to lda, 3nn, and svm in figure 3b is on account of classification optimization relative to the model, which characterizes prior information.
to further investigate obc improvement, we again considered heterogeneous covariance matrices but with independent features to determine if there is any difference in the relative performance between the classifiers.
in fact, the results provided in figure 3c show identical relative performance to the error curves in figure 3b, thereby indicating that both the standard classifiers and the obc, relative performance  is not affected by whether or not the features are correlated.
indeed, comparing figure 3a with figures 3b and 3c, we see that the relative performance of svm, mp obc, and calibrated mp obc is the same in both the homogeneous and heterogeneous models.
the switch in relative performance between lda and 3nn between figure 3a and figures 3b and 3c is not surprising because lda is optimal for a fixed  homogeneous gaussian model but not for a heterogeneous gaussian model.
the larger overall classification errors in figure 3a as compared to figures 3b and 3c are due to the different covariance matrices generated by the hc and ic models.
each model required different generating distributions for {σy,ρ} and {s,κ} for the hc and ic cases, respectively, and the particular choices made in section ‘synthetic data’ resulted in larger dispersions and higher errors in the hc models than the ic models.
to demonstrate this, we tested lda with  <dig> training and testing samples across  <dig> random generating distributions and found the average hc classification error to be  <dig>  and the ic error to be  <dig> .
this is despite lda being optimal for homogeneous, fixed, known gaussian cases and sub-optimal for heterogeneous, fixed, known gaussian cases, where the former is similar to the hc case.
still using independent-covariance data, we fixed the mean of class  <dig> at μ0= <dig>  in figure 3d, and varied μ <dig> from  <dig>  to  <dig>  to make the classification problem harder and easier, respectively.
across this range of classification problems, the mp obc had better classification performance than the other classification methods.
in addition, calibrated priors improved performance further, especially for harder classification problems.
real data
in table  <dig>  we chose ten genes at random from adenocarcinoma tumor tcga samples and performed model diagnostics  <cit> , p.  <dig>  by calculating posterior predictive p-values for interquartile distance  as a measure of dispersion.
in the additional file  <dig>  we provide additional test statistics and graphical predictive posterior model diagnostics.
these results indicate that rna-seq overdispersion is modeled sufficiently with the mp model.posterior predictive model diagnostics are given for  <dig> randomly selected genes from adenocarcinoma tcga samples
gene id
iqr
95% int.
for iqr
p-value
inter-quartile distance  is used as a robust measure of dispersion.
in the table, iqr is the training data’s iqr, followed by the 95-th credible interval, and the posterior predictive p-value.
in cases where the p-value is close to  <dig> or  <dig>  the true test statistic’s distance from the 95-th credible interval can be used to determine the magnitude of the mis-fit.
in figure  <dig>  we see mean holdout errors averaged over  <dig>  training sets and testing sets of tcga data as described in section ‘real data’.
here the mp obc performs better than all other classifiers across most training sample sizes considered, but calibration does not improve performance for this particular dataset.
recall that improvement owing to calibration depends on the extent to which the data satisfy the mp model.
if the aim of this paper were to build an operational classifier based on the tcga data, then we would have to go back and extensively study the data set to examine deviations from the model – for instance, outliers; however, here our aim is to show the functionality of the obc with non-gaussian data based on mcmc and apply it to the mp model.
the fact that the mp obc performs well on the real data satisfies this aim.
calibration is a tricky business and it would be a major separate study to characterize the manner in which model variation affects calibration, even if we were to perform an intensive study of this particular data set.
performance on the synthetic data demonstrates the effectiveness of the calibration when the model is satisfied.tcga rna-seq classification.
average holdout errors were computed over  <dig>  training sets and feature subsets using two types of lung cancer rna-seq data from tcga.
mp obc with and without calibrated priors demonstrates superior performs across a range of training sample sizes.
in addition, providing the mp obc with calibrated priors does not appear to improve performance in this particular dataset.
computational limitations
the results in figure  <dig> and figure  <dig> required tens of thousands of mcmc runs.
owing to limited available computational resources, we could only allocate around  <dig> seconds on a single cpu core for each mcmc run.
this necessitated using only four genes for these classification results as each iteration of the mcmc procedure has time complexity of o, where d is the number of features.
in practice, one would have a small number of data sets and could use parallel computing to devote more time and computing effort for the classification.
for example, in timescales on the order of hours on a typical workstation, we have successfully performed classification using  <dig> genes.
the other classification methods compared in this study have smaller computational requirements and can correspondingly handle larger numbers of features given the same available resources.
however, for the small sample sizes often available in biology,  <dig> genes is typically beyond the “peaking” point where most classifiers decrease in classification performance as more features are added   <cit> .
incidentally, the obc does not suffer this “peaking phenomenon” as shown in  <cit> .
in addition, the computational time requirements of classification is typically not a bottleneck in translational medicine given the timescales used in collecting biological data.
in these settings, the accuracy of classification is much more valuable than rapid runtimes, and this is the primary advantage of the computational obc framework proposed in this paper.
CONCLUSIONS
we have demonstrated that bayesian classification can be applied to specific problem domains such as rna-seq through statistical modeling and mcmc computation.
the resulting classifier provides superior classification performance compared to state-of-the-art classifiers such as svm with a radial basis kernel.
although we have not discussed error estimation – our interest in the present paper being classification, ipso facto, the mcmc approach to optimal bayesian classification can be applied, via  <cit>  and  <cit> , to obtain optimal mmse error estimators for any classification rule and sample-conditioned evaluation of the mse for error estimation.
future work includes examining the normalization parameter d and determining if additional performance improvements can be made by considering the distribution over d rather than transforming the original data through the process of data normalization.
additionally, more efficient computational techniques could be used to allow for larger feature sizes, including program optimization and utilizing structure in the feature covariance to reduce the size of the parameter space.
additional file
additional file  <dig>
supplementary materials.
algorithms and model diagnostics.
supporting details including in-depth algorithms and model diagnostic plots and figures are given in a single multi-page pdf.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
jk developed the mcmc implementation of the optimal bayesian classifier, worked on the mp modeling, and wrote the draft of the manuscript.
ii worked on the mp modeling, rna-seq modeling, and finalizing the manuscript.
ed worked on the application of the optimal bayesian classifier and finalizing the manuscript.
all authors read and approved the final manuscript.
funding
this work was supported in part by the national institute of health grant u01ca <dig>  the niehs center for translational environmental health research  p30es <dig>  and the center for nonlinear studies at the los alamos national laboratory.

 the possibility of genomic phenotype classification arose with the inception of gene-expression microarrays. from the outset, two fundamental problems have frustrated the endeavor: the inaccuracy of microarray measurements, and small samples. our particular application of interest is classification using rna-seq data. modern rna-seq technologies sequence small rna fragments to measure gene expression, where the number of reads mapped to a gene on the reference genome defines the count data. given that rna-seq data has advantages over microarray data, in particular, more accurate measurement, we still confront the second fundamental problem, which is statistical, not technological: small samples cause re-sampling-based classifier error estimators to be very inaccurate due to excessive variance and lack of regression with the true error. since the error rate of a classifier quantifies its predictive accuracy, it is the salient epistemological attribute of any classifier. the inability to satisfactorily estimate the error with model-free methods with small samples implies that genomic classifier error estimation is virtually impossible without the use of prior information, so that the whole small-sample classification problem becomes unapproachable in a model-free framework. the situation has been addressed by utilizing prior knowledge via a bayesian approach that considers a prior distribution on an uncertainty class of feature-label distributions. for expression-based classification, prior distributions have been constructed using expression data not employed in classifier design and known regulatory pathways. given that a prior model must be assumed to achieve satisfactory error estimation, an obvious course of action is to derive an optimal classifier based on the prior knowledge and the sample data, the result being an optimal bayesian classifier that is guaranteed to have the best average performance of any classifier relative to the posterior distribution derived from the prior distribution and data. while bayesian classification does not depend on particular distributional forms, closed-form solutions have been derived for the multinomial model and gaussian models using linear classifiers for the minimum mean squared error error estimate, the mse of the error estimate, and an optimal bayesian classifier relative to the prior distribution, the latter being expressed in terms of effective class conditional distributions, which are expectations relative to the posterior distribution of the class-conditional distributions. the closed-form solutions depend on particular models and the existence of conjugate priors, which can be too constraining for practical applications such as rna-seq classification. much of the statistical literature concerning classification of rna-seq data attempts to address differential expression testing, that is, univariate statistical testing on an individual gene basis. these attempts typically model rna-seq data via negative binomial and poisson distributions. in addition, network inference has been attempted using a hierarchical poisson log-normal model, and clustering of rna-seq data points has utilized various approaches. however, in clinical settings one is often interested in sample classification: the problem of classifying the rna-seq data from unlabeled patients using a set of labeled training data. one of the few rna-seq-specific attempts towards this goal uses a poisson modeling assumption with independent features. the poisson model is completely parameterized by its mean and thus is known to exhibit problems in fitting rna-seq data due to the overdispersion typically observed in such datasets. in this paper, we focus on modeling the pipeline that starts with extracting the gene concentrations from the biological samples and their subsequent processing by the sequencing instrument. this is accomplished using a hierarchical, multivariate poisson model. specifically, gene concentration levels are modeled by a log-normal distribution and the sequencing instrument sampling of those is modeled via a poisson process. this allows us to accurately model the rna-seq data overdispersion as demonstrated by marginal variance calculations and posterior predictive model diagnostics in section ‘overdispersion’. in addition, this hierarchical model allows for inferring any covariance structure observed between the features. whereas dalton and dougherty have presented a computational method for nonlinear classifiers in the gaussian model, this still depends upon conjugate priors. in this work, we remove the constraints imposed by the requirement of a closed-form solution by developing the optimal bayesian classifier using a markov-chain-monte-carlo methodology. this provides a computational framework for calculating the obc for any parameterized class conditional-density and any prior distribution. most notably, this allows us to use distributions designed to closely model particular datasets and a prior distribution of any form to improve classification performance in small-sample settings, in particular, for rna-seq-based classification. methods notation throughout, we use capital letters to indicate random variables, lower case letters to indicate individual realizations of random variables or indices, bold latin characters for observed vectors, and greek letters for latent features and parameters. we write p as the probability measure over the random variable x. p may be a probability mass function, probability density function, or arbitrary probability measure. p denotes the conditional probability p. similarly, following bayesian convention, we write parameterized distributions by conditioning on the parameter, for instance, p, and posterior expectations by conditioning on the sample, such as e, where sn and all other values are defined in section ‘review of optimal bayesian classification’. if it is unclear which density an expectation is taken with respect to, then we denote it in subscript notation, such as eθ|sn·, where the expectation is taken with respect to the density p. review of optimal bayesian classification binary classification considers a set of n labeled training data points, sn=1n, where y∈{ } is the class label and x∈x is the feature vector over a feature space. an example of binary classification in a clinical setting might include class and being two types of cancers, or normal and cancerous tissues. available features would then be the gene or genes that will eventually be used in the designed classifier to assign this label. the feature space would be the set of possible gene expression measurements for all genes in the feature vector. the labeled training data sn would be the set of gene expression measurements from samples which had undergone further testing to identify the type or malignancy of the tissue. using sn, we design a classifier ψ that hopefully performs well on the unknown joint feature-label distribution p. in the same clinical example, the classifier ψ could then identify the type of cancer using gene expression measurements alone. by parameterizing this unknown joint distribution in a model-based bayesian framework one can derive an optimal bayesian classifier that minimizes the expected error over the space of all classifiers under assumed forms of the class-conditional densities. specifically, under gaussian and multinomial class-conditional densities and their corresponding conjugate prior distributions, closed-form solutions for the obc and the first two moments of the error estimate conditioned on the sample have been obtained. the parameterization of the feature-label distribution consists of the marginal class probability c and the class-conditional densities p, where a particular value θy∈Θy specifies a single class-conditional density contained in the class of densities defined over the space Θy, which will be a cartesian product as described in section ‘the multivariate poisson model’. therefore, for a two-class problem, we specify a parameterized joint feature-label distribution as θ=∈Θ= ×Θ0×Θ in the bayesian classification framework, these values are then treated as random variables, so that we may consider quantities such as the expectation of c, or another random variable conditioned on the value of the parameter vector θ. figure bayesian classification derivation tree. a tree summarizing the relationships between several important quantities in the general theoretical framework of bayesian classification. a directed edge between a parent and its child indicates that the child can be derived from the parent by the equations indicated in the edge label. the root of the tree p is the posterior distribution of the feature label parameters and by taking expectations with respect to this distribution, we can derive the effective class conditional densities p and the distribution of the classifier error p. then these quantities give rise to the obc, and mmse and mse estimates for the error as described in the text. quantities highlighted in grey are given in closed form for gaussian and multinomial distributions in. we now examine the tree in more detail. starting at the far left of the tree, p is the posterior distribution of the parameterized feature-label distribution – posterior to the labeled samples in sn. typically, error estimates and the optimal classifier are our primary interest, so that this posterior distribution is traditionally used as a means to compute other quantities and is not of interest by itself. the effective class-conditional density is the marginal predictive posterior of the feature vector x conditioned sn and the class variable y, px|y,sn=∫Θypx|y,θypθy|sndθy. it gives the distribution of the feature vector using a weighted average over all the parameterized class-conditional densities in Θy given a class y. the weights in this expectation are the posterior, p, evaluated at each θy. the true error of classifier ψ is ε=p≠y). given the sample data sn, ε is a random unknown quantity in the bayesian framework. the mmse estimate given in can be written as eε|sn=pψ≠y|sn=eθ|snpψ≠y|θ,sn=ĉε0θ ψ+ε1θ ψ=∫xĉpx| snix∈r1+px| snix∈r0dx, where ia is the indicator function for event a, ĉ=ec|sn is the posterior expectation of c, ry is the region of the feature space the classifier predicts to be class y, is the feature space, and εy is the error of classifier ψ contributed by class y on the fixed distribution θy. we can also obtain the full posterior distribution of the error, pε|sn=∫Θpε|θpθ|sndθ=eθ|snpε|θ, where p is the true error for a fixed feature-label distribution and fixed classifier. we denote this deterministic function by ε. as shown in figure the mmse estimate and the sample conditioned mse for this error can also be calculated using the first two moments of the error distribution. with the mmse estimator defined, the optimal bayesian classifier is the classifier minimizing the expected error by pointwise minimization of the integral: ψobc=0ifĉp≥p,1otherwise.. conditional error estimator if the true feature-label distribution were known, then we could compute the true error of a classifier exactly as an expectation over the conditional error: ε=p≠y)=∫xp≠y|x)pdx. treating ε as a random variable, one can similarly derive its posterior distribution by conditioning on the feature vector: pε|sn=∫xpε,x|sndx=∫x∫Θpε,θ,x|sndθdx=∫Θpθ|sn∫xpε|x,θpx|sndxdθ, which is different than the derivation of the same quantity in. this introduces the idea of the conditional error estimator, which we define as the mmse estimate of the classification error conditioned on the feature vector x, ε^=eθ|snε|x,sn=pψ≠y|x,sn=px|y≠ψ,snpy≠ψ|snpx|sn=z−1px|y≠ψ,snpy≠ψ|sn, as expanded through application of bayes’ theorem, where z is a normalizing constant given by z=px|sn=∑y∈ px|y,snpy|sn. in addition to being useful in the above alternative derivation of the classifier’s error posterior, the conditional error estimate has other practical applications. when classifying an unlabeled data point, we would like to estimate the error of the classifier output for that particular data point, as opposed to the overall error estimate for the classifier. for the obc, from the conditional error estimator can be written as ε^ψobc,x=z−1miny∈ px|y,snpy|sn. in sum, using the effective class-conditional densities and the posterior marginal probabilities one can calculate conditional error estimates for points in the feature space in addition to the earlier quantities described. the multivariate poisson model with the widespread use of next-generation sequencing techniques, classification approaches must be developed to account for the discrete nature of the mapped sequence data and to accommodate the various types of prior information available regarding these experiments. gene concentration levels can be modeled using a log-normal distribution. as discussed in the introduction, we assume that the sequencing instrument samples this mrna concentration through a poisson process and obtains xi,j reads for sample point i and gene j. we model this as pxi,j|λi,j∼poissondiexpλi,j, where λi,j is the location parameter of the log-normal distribution for sample i and gene j, and di is a variable accounting for the sequencing depth as determined by the sequencing process. for each i, we model the location parameter vector λi with a multivariate gaussian distribution, λi∼normal. we then consider the mean μ and covariance Σ of the gene concentrations as independent quantities for each class y. the entire mp model is represented in figure as a plate diagram. the distribution of a single class y is parameterized by θy=, where d= and λ=,i=,…n,j=,…,d, for n sample points and d total genes. therefore, θy∈Θy=ℝd×ℝd×d×ℝn×ℝd×n. the feature-label distribution parameterization for the two-class problem is then given by θ=, where c=p, the prior probability for class multivariate poisson model plate diagram. a plate diagram for the multivariate poisson model. the outermost plate represents the classes that we are interested in classifying against, where i is the index of the sample in class y, and j are the genes being modeled. to ensure a proper posterior with unit integral, we place weakly informative priors over the latent variables in the mp model. in choosing these values, we have aimed to avoid the complications that can occur with overly diffuse priors, such as lindley’s paradox. we choose: μy∼normalηy,ν2idΣy∼inverse-wishartκy,syc∼beta, where each element of μy is distributed according to a univariate gaussian. unless otherwise stated, η is the d dimensional zero vector, ν2= κ= and s=id. for computational and identifiability reasons, d is fixed to be a vector of normalization constants in order to match the different sequencing depths across all the samples. in practice, d can be approximated by an upper quartile normalization, which has been shown to be effective. in any bayesian approach the choice of prior affects the results, especially when only a few data points are given. in the case of mmse classifier error estimation in the bayesian framework, robustness to incorrect modeling assumptions has been extensively studied in and in those studies performance held up well for various kinds of incorrect modeling assumptions. robustness of optimal bayesian classifiers to false modeling assumptions was extensively studied in. again, good robustness was exhibited. of course, one can get bad small-sample results by intentionally selecting an inaccurate prior. in general, if one is confident in his knowledge, then a tight prior is called for because tighter priors require less data for good performance; on the other hand, when one is not confident, then prudence calls for a less informative prior. as proven in, obc classification is consistent under very general conditions; however, a prior whose mass is concentrated far away from the true parameters will perform worse than one that is non-informative. these issues have been extensively discussed in the bayesian literature. in the end, performance is the measure of worth and our results with synthetic and real data indicate solid performance for the modeling approach used herein. overdispersion the mp model uses the poisson distribution in a hierarchical scheme. it is important to note that, while the read counts are modeled as conditionally poisson in equation the observed read counts are not marginally poisson distributed. to demonstrate this, consider a one-dimensional simplification of the mp model in which x is the number of reads observed, λ is the log of the rna concentration, and λ∼normalμ,σ2x∼poisson). then for the marginal variance of x, var=evar+var=e+e≥eμ=varpoissoneμ where μ and σ are the mean and variance of the log of the concentration. therefore, when σ2> the marginal variance of x is always greater than that of a poisson random variable with the same effective rate. in addition, by carrying out a posterior predictive model check, p. by computing marginal posterior p-values against real rna-seq data, we can quantitatively assess the ability of the mp model to fit the dispersion of the tcga data. for a test statistic t, we compute the p-value by comparing the test statistic on the true data t and the value of the statistic averaged across the posterior predictive distribution t, where xrep∼p: pt=pr≥t|sn)=pr>t|sn)+pr=t|sn)≈1m∑i=0mi{t)>t}+ i{t)=t}, where xrep are monte carlo samples taken from the posterior predictive distribution p using the m monte carlo samples from the posterior distribution of θ as described in section ‘computation’. the term pr=t|sn) is necessary due to the discrete nature of rna-seq data. p-values away from and indicate that the model posterior produces test statistics both above and below that measured on the real data. we also consider where the real test statistic falls in relation to credible intervals of the test statistic to consider the magnitude of any differences. we apply the inter-quartile distance test statistic to provide a measure of the mp model’s ability to fit the dispersion of rna-seq data. we also consider several other test quantities in the additional file 1: table s1-s prior calibration using discarded features since designed classifiers typically use very few of the totality of observed genes, only a small fraction of the data is used for classifier design. similarly to, we can use the discarded features to calibrate the inverse-wishart prior for our mp obc. our goal is to obtain hyperparameters s,m,κ, and ν for each class from our training data sn. in general, we do not expect the discarded features to give us information about any particular genes and the specific covariances between genes, so we make the simplifying assumptions that we learn information from the discarded genes in an aggregate sense. thus, we consider the following structure on the hyperparameters: m=mt and where m∈ℝ, σ2> and −1≤ρ≤ for each class, we need to determine values for five scalar quantities: m,ν σ ρ, and κ. due to the hierarchical design of the mp model, we cannot apply the method of moments in a direct fashion, as did. instead, we utilize a sampling based approach to the method of moments. this mcmc sampling approach has been examined in as an extension to the generalized method of moments. the sampling approach uses the discarded features in an additional mcmc run evaluated prior to the primary classification mcmc procedure as discussed in section ‘computation’ – and then proceeds to the method of moments. in this calibration mcmc, we initialize all prior distributions with flat priors and use the discarded features to obtain samples from the posterior distribution of μ and Σ. typically, the number of discarded features f is much larger than the dimensionality d of the classification problem. therefore, due to computation time, we uniformly sample fs pairs of features from f and average the resulting runs rather than using all or large groups of discarded features in a single mcmc run. we use the following procedure: for each randomly chosen discarded feature pair: obtain mcmc samples using the feature pair as data and flat priors. record posterior averages of μ and Σ. average over these posterior averages as given by equations -. using the resulting five hyperparameter estimates, run the final mcmc for classification. following, we use the moments of the posterior samples to determine the hyperparameters through the following relations: the mean of an inverse-wishart distribution is eΣ=sκ−d− which together with our simplified covariance structure implies σ2=κ−d−1eΣ ρ=eΣ12eΣ the variance of the first diagonal of an inverse-wishart matrix can be used to solve for κ via κ=2eΣ112varΣ11+d+ as we have samples of μ directly from our posterior, we obtain m=eμ ν=varμ in order to use equations -, we obtain estimates of the moments from mcmc performed over the fs discarded feature pairs. for the i-th feature pair we obtain the posterior means μ^ Σ^ and Σ^ and then average: e^μ1=1fs∑i=1fsμ^ var^μ1=1fs−1∑i=1fse^μ1−μ^ e^Σ11=1fs∑i=1fsΣ^11+Σ^ e^Σ12=1fs∑i=1fsΣ^ var^Σ11=1fs−1∑i=1fse^Σ11−Σ^ we substitute the estimates from equations - back into equations - to obtain the final hyperparameter estimates. one must keep in mind that the calibration procedure explicitly assumes the mp model. hence, one can only expect an improvement in the classification performance if the data follow the mp model. computation to obtain the mp obc, we approximate the effective class conditional densities in order to minimize the expected error in a pointwise fashion: px|y,sn=∫Θypx|y,θypθy|sndθy≈1m∑i=1mpx|y,θy, where θy are m samples of θy from the model posterior distributions. for clarity of presentation, we do not consider the class variable y, and we assume a single class. we do this because the computation can be performed per-class due to the assumed independence between the classes and the marginal probability, p=ppp. to obtain posterior samples of θ using the metropolis hastings mcmc algorithm we define a proposal distribution p to obtain a new value for the class parameters θ′ from the old values θ. we then calculate the acceptance ratio r=min pθ′|snpθ′|θpθ|snpθ|θ′=min psn|θ′pθ′psn|θpθ, under the assumption of a symmetric proposal distribution =p). the process of proposing and accepting samples from this distribution with the probability r induces a markov chain. positivity of the proposal distribution > for any θ) is a sufficient condition for ergodicity of this markov chain. furthermore, this markov chain admits a steady-state distribution equal to our desired posterior distribution p. from the definition of the likelihood, psn|θ=∏i=1npxi|θ=∏i=1npxi|λi=∏i=1n∏d=1dpxi,d|λi,d, where p=p owing to conditional independence. from the definition of the prior, p=p=ppp=∏i=1nppp. the posterior predictive distribution in is approximated by px|sn≈1m∑i=1mpx|θ=1m∑i=1m∫Λpλ|θpx|λdλ=1m∑i=1m∫Λpλ|θ∏k=1dpxk|λkdλ≈1mt∑i=1m∑g=1t∏k=1dpxk|λk, where, p∼poisson), λ∼ normal, Λ=ℝn×d, and the λ are t vector-valued samples drawn from the appropriate class’s posterior distribution used to approximate the inner intractable integral. in addition, we use this approximation of the effective class-conditional density to calculate the conditional error estimates of in a pointwise fashion. finally, because we have assumed a conjugate prior distribution for the marginal class probability c, the posterior expectation takes the closed form eθ|snc=n0+α0n0+n1+α0+α where the ny are the number of training samples obtained from class y and the αy are hyperparameters set to for an uninformative prior. conjugacy was used for this one parameter because the increased flexibility of the full sampling approach was deemed not necessary due to the constrained, univariate nature of the parameter. if more complex relationships between c and other parameters were desired, then a sampling approach using non-conjugate priors would be straightforward to implement. synthetic data to evaluate obc performance in the setting of the mp model, we generate synthetic data using the method proposed in to simulate gene expression/mrna concentrations. these gene expression values are then statistically sampled to emulate modern sequencing machines as described in. parameter values are drawn from the following distributions to examine a wide variety of classification problems: μy∼normal,σy∼inverse-gamma,ρ=uniform,dlow= dhigh= blocksize= with these parameters, ten global, twenty heterogeneous, and ten non-marker features are generated. then four features are randomly chosen to represent a mixture of features of various classification quality. following, the features in the data are zero mean and unit standard deviation normalized except for the mp obc. the exception occurs because the mp model expects features to be positive integers and normalization is not necessary. the discarded features are used for calibration of the mp obc priors, and samples are generated from each class to estimate the true classification rate for each classifier. we use four features in this synthetic data classification study owing to limited computational resources as discussed in section ‘computational limitations’. the synthetic data generation method proposed in imposes the strong assumption of a homogeneous covariance structure between the two classes of data. this assumption does not hold for biological situations where interactions between features are not necessarily preserved between classes, and this occurs frequently in biology when considering the possible effects of canalizing genes, nonlinear gene regulation, and mutations in the case of cancer. specifically, if the canalizing gene is not observed, and differs in activity between the two classes, then the measured correlation between two downstream genes could potentially be negligible in one class while strong in the other class. similarly, for highly nonlinear gene regulation, if a gene in one class is in the saturation region of its response curve from a master gene, then the correlation will be low, while a lower expression level in the other class would allow for a large measured correlation with the same canalizing gene. and finally, if one class represents normal gene expression and the other tumor-related expression, then a correlation might exist from a functioning pathway in the normal tissue, but a mutation could result in a lack of correlation effects in the tumor. hence, we modify the synthetic data generation procedure in an attempt to produce synthetic datasets more representative of such nonlinear phenomena in biology. in this modified procedure, we allow independent covariance matrices for the features of the two classes. to generate these covariance matrices, Σy, we utilize independent draws from inverse-wishart distributions with parameters κy= d= and scale matrix s=σy id. to examine the effects of feature correlation in ic datasets, we can also generate low-correlation covariance matrices by zeroing the off-diagonal terms. once the covariance matrix for class y is obtained, location parameters for gene-expression values for each sample point are drawn from the respective multivariate normal distribution λy∼n. each sample point is then assumed to be normalized through an upper quartile or other suitable method, but in practice any sample-based normalization is imperfect. we reflect this variation by drawing the sequencing depth di from a uniform distribution, giving the rate of the poisson process as di expλi. the number of reads for a single gene from a single sample is then drawn from this poisson distribution. see additional file for more detail. the obc is optimal on average across the space of distributions determined by its prior distributions. to avoid biasing the performance comparison, we draw the classification problem datasets using different distributions than those of the obc priors. see additional file for more detail. real data we consider a real rna-seq dataset composed of level rnaseqv data from the cancer genome atlas project. it contains and specimens from lung adenocarcinoma and lung squamous cell carcinoma tumor biopsies, respectively. the samples are mapped read counts against known human rna transcripts as generated by the university of north carolina at chapel hill, one of the genome sequencing centers for the tcga. the data for each cancer type is the result of processing approximately billion reads and the read count files for each are one gigabyte apiece. the problem is to classify the tumor types. because the class- and class- sample sizes, and are not chosen randomly, we are confronted with the problem of separate sampling. this means that there is no way to obtain a posterior distribution for c and therefore c must be known in advance. based upon records from 2006- we have a very accurate estimate, / ≈. whereas we can use the value of c directly, along with all of the data, in designing the obc, for classification rules that do not use c explicitly, the separately sampled data must be maximally subsampled to the proper sampling ratio c before applying the classification rule. this means that for ntrn desired samples, the sample subsets will contain roundntrn) and round for class and respectively. moreover, holdout error estimation, which we use here, must be properly adapted for separate sampling for all design methods, including the obc. the holdout estimate is given by ε^c=cε^0+ε^ where ε^ and ε^ are the ordinary holdout estimators for the class- and class- errors, respectively. we note that many studies have made the mistake of using classification rules designed for random sampling when sampling is separate. this can have devastating effects on classifier performance. while averaging over sample subsets for holdout error estimation, we also average over uniformly, randomly selected gene subsets of size this sampling occurs from low expression genes. we sample from these lower expression genes because we are ultimately interested in classification problems where the delineation between phenotypes is determined by genes with low expression. we used for averaging in order to obtain a large enough sample over this feature and sample subset space to achieve repeatable results. computational runtime for each sample and gene subset was similar to the synthetic data.  the additional file contains a simple two-class, two-feature demonstration of the overall procedure to allow for easy visualization and interpretation. here we discuss the results for the synthetic and real data. synthetic data to evaluate classification performance, classifiers were trained using 3nn, lda, and c-support vector machines with a radial basis function kernel. starting with the homogeneous-covariance model, figure 3a shows that the performance of the multivariate poisson obc is better than nonlinear svm when more than samples are available and is significantly better than any other classifier when using calibrated features. equivalently, by using discarded features, we can obtain the same classification performance while requiring fewer training samples.synthetic rna-seq classification. synthetic data classification results with homogeneous-covariance, high correlation independent-covariance, low correlation independent-covariance, and high correlation independent-covariance data at several problem difficulties. in the case of independent-covariance data with highly correlated features, figure 3b shows superior classification performance of the mp obc at nearly all sample sizes considered. in addition, for calibrated prior distributions, the performance of the mp obc improves. this improvement is greater when the sample sizes are small, which demonstrates the importance of additional knowledge when data are expensive to obtain or not readily available. the superior performance of the obc relative to lda, 3nn, and svm in figure 3b is on account of classification optimization relative to the model, which characterizes prior information. to further investigate obc improvement, we again considered heterogeneous covariance matrices but with independent features to determine if there is any difference in the relative performance between the classifiers. in fact, the results provided in figure 3c show identical relative performance to the error curves in figure 3b, thereby indicating that both the standard classifiers and the obc, relative performance is not affected by whether or not the features are correlated. indeed, comparing figure 3a with figures 3b and 3c, we see that the relative performance of svm, mp obc, and calibrated mp obc is the same in both the homogeneous and heterogeneous models. the switch in relative performance between lda and 3nn between figure 3a and figures 3b and 3c is not surprising because lda is optimal for a fixed homogeneous gaussian model but not for a heterogeneous gaussian model. the larger overall classification errors in figure 3a as compared to figures 3b and 3c are due to the different covariance matrices generated by the hc and ic models. each model required different generating distributions for {σy,ρ} and {s,κ} for the hc and ic cases, respectively, and the particular choices made in section ‘synthetic data’ resulted in larger dispersions and higher errors in the hc models than the ic models. to demonstrate this, we tested lda with training and testing samples across random generating distributions and found the average hc classification error to be and the ic error to be. this is despite lda being optimal for homogeneous, fixed, known gaussian cases and sub-optimal for heterogeneous, fixed, known gaussian cases, where the former is similar to the hc case. still using independent-covariance data, we fixed the mean of class at μ0= in figure 3d, and varied μ from to to make the classification problem harder and easier, respectively. across this range of classification problems, the mp obc had better classification performance than the other classification methods. in addition, calibrated priors improved performance further, especially for harder classification problems. real data in table we chose ten genes at random from adenocarcinoma tumor tcga samples and performed model diagnostics, p. by calculating posterior predictive p-values for interquartile distance as a measure of dispersion. in the additional file we provide additional test statistics and graphical predictive posterior model diagnostics. these results indicate that rna-seq overdispersion is modeled sufficiently with the mp model.posterior predictive model diagnostics are given for randomly selected genes from adenocarcinoma tcga samples gene id iqr 95% int. for iqr p-value inter-quartile distance is used as a robust measure of dispersion. in the table, iqr is the training data’s iqr, followed by the 95-th credible interval, and the posterior predictive p-value. in cases where the p-value is close to or the true test statistic’s distance from the 95-th credible interval can be used to determine the magnitude of the mis-fit. in figure we see mean holdout errors averaged over training sets and testing sets of tcga data as described in section ‘real data’. here the mp obc performs better than all other classifiers across most training sample sizes considered, but calibration does not improve performance for this particular dataset. recall that improvement owing to calibration depends on the extent to which the data satisfy the mp model. if the aim of this paper were to build an operational classifier based on the tcga data, then we would have to go back and extensively study the data set to examine deviations from the model – for instance, outliers; however, here our aim is to show the functionality of the obc with non-gaussian data based on mcmc and apply it to the mp model. the fact that the mp obc performs well on the real data satisfies this aim. calibration is a tricky business and it would be a major separate study to characterize the manner in which model variation affects calibration, even if we were to perform an intensive study of this particular data set. performance on the synthetic data demonstrates the effectiveness of the calibration when the model is satisfied.tcga rna-seq classification. average holdout errors were computed over training sets and feature subsets using two types of lung cancer rna-seq data from tcga. mp obc with and without calibrated priors demonstrates superior performs across a range of training sample sizes. in addition, providing the mp obc with calibrated priors does not appear to improve performance in this particular dataset. computational limitations the results in figure and figure required tens of thousands of mcmc runs. owing to limited available computational resources, we could only allocate around seconds on a single cpu core for each mcmc run. this necessitated using only four genes for these classification results as each iteration of the mcmc procedure has time complexity of o, where d is the number of features. in practice, one would have a small number of data sets and could use parallel computing to devote more time and computing effort for the classification. for example, in timescales on the order of hours on a typical workstation, we have successfully performed classification using genes. the other classification methods compared in this study have smaller computational requirements and can correspondingly handle larger numbers of features given the same available resources. however, for the small sample sizes often available in biology, genes is typically beyond the “peaking” point where most classifiers decrease in classification performance as more features are added. incidentally, the obc does not suffer this “peaking phenomenon” as shown in. in addition, the computational time requirements of classification is typically not a bottleneck in translational medicine given the timescales used in collecting biological data. in these settings, the accuracy of classification is much more valuable than rapid runtimes, and this is the primary advantage of the computational obc framework proposed in this paper.  we have demonstrated that bayesian classification can be applied to specific problem domains such as rna-seq through statistical modeling and mcmc computation. the resulting classifier provides superior classification performance compared to state-of-the-art classifiers such as svm with a radial basis kernel. although we have not discussed error estimation – our interest in the present paper being classification, ipso facto, the mcmc approach to optimal bayesian classification can be applied, via and, to obtain optimal mmse error estimators for any classification rule and sample-conditioned evaluation of the mse for error estimation. future work includes examining the normalization parameter d and determining if additional performance improvements can be made by considering the distribution over d rather than transforming the original data through the process of data normalization. additionally, more efficient computational techniques could be used to allow for larger feature sizes, including program optimization and utilizing structure in the feature covariance to reduce the size of the parameter space. additional file additional file supplementary materials. algorithms and model diagnostics. supporting details including in-depth algorithms and model diagnostic plots and figures are given in a single multi-page pdf. competing interests the authors declare that they have no competing interests. authors’ contributions jk developed the mcmc implementation of the optimal bayesian classifier, worked on the mp modeling, and wrote the draft of the manuscript. ii worked on the mp modeling, rna-seq modeling, and finalizing the manuscript. ed worked on the application of the optimal bayesian classifier and finalizing the manuscript. all authors read and approved the final manuscript. funding this work was supported in part by the national institute of health grant u01ca the niehs center for translational environmental health research p30es and the center for nonlinear studies at the los alamos national laboratory. 
we also introduce the bayesian minimum mean squared error conditional error estimator and demonstrate its computation over the feature space. through model-based, optimal bayesian classification, we demonstrate superior classification performance for both synthetic and real rna-seq datasets. thus, we introduce a hierarchical multivariate poisson model and the associated optimal bayesian classifier for classifying samples using sequencing data. we demonstrate superior or equivalent classification performance compared to typical classifiers for two synthetic datasets and over a range of classification problem difficulties. lacking closed-form solutions, we employ a monte carlo markov chain approach to perform classification.
 sequencing datasets consist of a finite number of reads which map to specific regions of a reference genome. most effort in modeling these datasets focuses on the detection of univariate differentially expressed genes. however, for classification, we must consider multiple genes and their interactions.  thus, we introduce a hierarchical multivariate poisson model and the associated optimal bayesian classifier for classifying samples using sequencing data. lacking closed-form solutions, we employ a monte carlo markov chain approach to perform classification. we demonstrate superior or equivalent classification performance compared to typical classifiers for two synthetic datasets and over a range of classification problem difficulties. we also introduce the bayesian minimum mean squared error conditional error estimator and demonstrate its computation over the feature space. in addition, we demonstrate superior or leading class performance over an rna-seq dataset containing two lung cancer tumor types from the cancer genome atlas.  through model-based, optimal bayesian classification, we demonstrate superior classification performance for both synthetic and real rna-seq datasets. a tutorial video and python source code is available under an open source license at http://bit.ly/1gimnss. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords classificationrna-seqmodel-basedbayesianissue-copyright-statement© the author 2014 