BACKGROUND
tag-based transcriptome sequencing libraries consist of a collection of short sequences of dna called tags along with tabulated counts of the number of times each tag is observed in a sample.
these observed tag counts represent a sample from a much larger pool of mrna tags in a tissue or organism.
in the past, serial analysis of gene expression  was the most commonly used technology for generating libraries of tag counts.
sage libraries have been used to address a number of biological questions including: estimating transcriptome size, and estimating the density of relative expression levels.
most frequently, sage was used to assess differential expression across cells from different tissues or strains, or cells grown under different experimental conditions.
next generation methods such as digital gene expression  tag profiling  <cit>  now provide a more efficient method to generate tag libraries and are growing in popularity.
libraries based on dge are used to address the same questions and provide much larger numbers of tags leading to increased statistical power.
however, the close similarities between dge and sage, particularly the use of restriction enzymes, lead both techniques to share the same inherent biases.
as has been repeatedly noted in both sage and dge studies e.g.
<cit> , the mrna from a single gene can lead to the production of more than one type of tag within the same library.
observing multiple tags from a single gene complicates data analysis and has a number of potential causes including alternative mrna splicing and incomplete cdna digestion.
incomplete digestion typically results in a number of different tags being formed from the same mrna transcript.
alternative splicing can lead to multiple mrna isoforms within a sample.
when combined with incomplete digestion, each of these isoforms may give rise to a different tag or group of tags.
while it is intuitively sensible to analyze the prevalence of an mrna by aggregating multiple tags derived from it, such tags are often analyzed separately  <cit> .
at best, this approach diminishes the statistical power that is available when testing differences in mrna abundances and reduces the accuracy of inferences based on asymptotic statistical methods .
more seriously, such comparisons can result in dramatic biases leading to erroneous conclusions in testing differential expression due to varying levels of incomplete digestion across libraries.
the current work extends an earlier analysis  <cit>  based on sage libraries generated using the yeast saccharomyces cerevisiae  <cit> .
alternative splicing is generally rare in microorganisms such as s. cerevisiae so that the tag multiplicity observed in these libraries is primarily attributed to incomplete cdna digestion.
the first objective of the current work is to provide a methodology that allows multiple tags, which arise due to incomplete digestion, to be combined and used to infer expression levels of the underlying mrna transcripts.
aggregation is complicated by the fact that certain tags are more likely to form than others due to processing methods  <cit> .
unless the tag formation probabilities are accounted for, inferences based on aggregate tag counts may be highly biased.
previously, we have addressed the problem of estimation bias by directly modeling the variation in tag formation and deriving bias corrected maximum likelihood and posterior mode estimators of the composition of the mrna population  <cit> .
however, a number of assumptions and difficulties limit the application of these techniques.
we address these in the present study.
three new hierarchical models for sage data are introduced that directly correct for the bias by embedding a model of tag formation within the probability distribution describing the data.
corresponding gibbs sampling algorithms suitable for model fitting and inference are also derived.
these methods are computationally efficient and allow a more general class of prior distributions than earlier analytic approaches.
the results are also applicable for both aggregate and standard tag counts.
we use simulation studies and analysis of publicly available data from s. cerevisiae  <cit>  in order to contrast the accuracy of our sampling methods and analyze the effects of priors.
the approach we present should serve as a building block for a more comprehensive method that takes into account the effects of both incomplete digestion and alternative splicing.
a second major objective of the current work is to explore the general effects of tag formation bias on differential expression.
we carefully demonstrate that, in the absence of proper adjustments, variation in tag formation probability across experiments can lead to both false positive and false negative conclusions regarding the differential expression of a gene.
a brief description of the implementation of proposed algorithms in differential expression analysis is also provided.
while the approaches presented here require more work in processing and preparation than studies that focus on individual tags, they may offer a significant advantage with respect to statistical power in studies of differential expression.
a number of statistical methodologies have been developed to analyze sage libraries.
for studies that focus on relative expression levels, individual tags are viewed as possible outcomes from a multinomial sampling model  <cit> .
<cit>  offered an improvement on this by directly applying a bayesian multinomial-mixture dirichlet model to the observed vector of tag counts which improves accuracy in estimation of tag frequency.
marginal analysis may be more useful in studies of transcriptome size and was developed by  <cit>  and  <cit> .
statistical methods to assess differential expression across a number of samples were discussed by  <cit> .
these analyses typically test the differential expression of each tag across libraries independently; see  <cit> .
the sage methodology
the goal of a sage experiment  <cit>  is to sample the mrna population within a group of cells.
broadly speaking, the sage method generates a set of short sequence-based cdna tags from the mrna population.
initially, a pool of mrna is extracted from a group of cells.
the unstable single stranded mrna is reverse transcribed into a double stranded dna copy  using a modified primer that allows the cdna to be bound to a streptavidin bead.
using restriction enzymes, the bead-bound portion of cdna is cut into small 'tags' that are concatenated together into longer cdna molecules referred to as concatemers.
these concatemers are amplified using pcr and then sequenced.
a cleavage motif of the anchoring enzyme allows the start and stop points of tags to be identified within a concatemer.
each time a tag from a specific gene is observed in the sequence data, it contributes a count to the library.
the dge technology follows a similar procedure in the first several steps.
however, tags are neither concatenated nor cloned but sequenced immediately.
in both technologies the restriction enzymes used to create the tags cut at very specific sites within the cdna.
for example, the restriction enzyme used by  <cit>  cleaves only at the four nucleotide motif catg.
thus, tags can only be created at specific points within a cdna.
because the site where the cdna is cut serves as an 'anchor' between tags, the restriction enzyme used is often referred to as the anchoring enzyme .
we refer to the specific points in the cdna cleaved by the anchoring enzyme as ae sites.
while some genes will have no ae sites, many genes have multiple sites at which the ae could cleave the cdna.
given that the ae is expected to act in a site-independent manner, a single cdna molecule can be cut by the ae in multiple places.
however, because only the fragment of cdna that is attached to a streptavidin bead is retained during the experimental process, the site closest to the bead  that is actually cleaved is the only site that can lead to an observable tag .
if the ae worked with 100% efficiency, then each mrna could only lead to one observable tag.
however, the cutting efficacy of the ae is always less than 100%, sometimes referred to as incomplete digestion, and as a result multiple tags may be generated from the cdna of a single gene.
the critical point is that differences in the number of ae sites between genes result in different probabilities of tag formation.
the most extreme cases are genes whose mrna transcript lacks any ae site.
such genes are impossible to observe using the sage methodology.
for genes with ae sites, differences in digestion efficiencies lead to differences in tag formation probabilities between experiments.
a further complication is that many of these sites lead to non-unique or 'ambiguous' tags that cannot be readily assigned to a particular gene.
ambiguous tags are 'uninformative' using current technologies; tag-to-gene mapping is discussed at length in  <cit> .
the tags generated in the experiment analyzed below were  <dig> bp long leading to a number of ambiguous tags.
experimental advances, for example, 'supersage' techniques can lead to tags up to  <dig> bp long  <cit>  significantly diminishing this problem.
modeling tag formation probabilities
the previous discussion leads us to postulate that genes with larger numbers of ae sites have a increased probability of tag formation and further that there is large variation in the probability of tag formation within a gene.
below we present a quantitative model for tag formation probabilities within a gene that is key to later bias correction methods.
because only cdna that is attached to a streptavidin bead is retained during the experimental process, tags are created from the 3' most ae site that is actually cleaved .
let ki be the number of restriction enzyme sites which may be cleaved by the ae on mrna generated from gene i, let p be the probability the ae will cleave a site, and assume that cleavage probability is independent between sites and does not vary between genes.
if we label sites  <dig> through ki starting at the 3' most site , the probability of generating a tag through ae cleavage at site j ∈ { <dig>  ..., ki} follows a geometric distribution, j-1p.
this corresponds to the probability of no cleaving in sites  <dig> to j -  <dig>  and a cleaving at site j. note that this probability is independent of what happens at the ae sites 5' from the jth site.
the fact that the expected distribution of tags varies with ae cleavage probability p can be used to estimate p from the library of tag counts; see  <cit> .
it follows that the probability of generating any of the possible tags from the ith category of mrna follows a geometric distribution
where the sum is over all non-ambiguous tagging sites.
we refer to this quantity as the tag formation probability.
for purposes of our analysis, we will assume that, because each of the thousands of genes provides an independent estimate, the value of p is known to a very high precision.
the variation in the tag formation probability between genes stems from variation in the number of ae sites as well as the number of ambiguous tagging sites.
given p, the only uncertainty in the estimate of ϕi results from ambiguous tags; this problem may be partially eliminated through techniques that generate longer tags; see  <cit> .
hence, for methodological purposes we reasonably assume the ϕi to be known constants .
as a result, the tag pool represents a biased sample of the mrna population of interest being a function of both the known distribution of tag formation probabilities across the genome, and the distribution of these genes within the mrna population itself.
traditionally, a gene's frequency in the tag pool was assumed to be an unbiased estimate of its frequency in the mrna population, the population of interest.
this equivalence only holds when all genes have the same tag formation probability, a condition that is never met.
consequently, genes with greater than average tag formation probabilities will be over-represented while those with lower than the average probability will be under-represented.
equation  allows us to accurately estimate this sampling bias and correct the inferences made from tag libraries.
before the tag counts are considered, we assume that the data is processed to retain only informative tags, i.e.
all ambiguous or orphan tags are removed, as is standard practice.
let the vector t =  represent the aggregated observed tag counts for individual genes.
here l is the number of genes that contain at least one unambiguous ae site and ti is the total number of observed tags that can uniquely be attributed to an ae site in gene i.
.)
it is natural to view this vector as a sample from a multinomial distribution with l categories, t ~multi  where θ = , θi represents the frequency of tags from gene i in the tag pool, and  is the total count of all observed tags.
because the variation in tag formation probabilities ϕi distorts the abundances of observed tags, it is clear that θi represents a biased estimate of the true proportion of mrna from gene i in the overall pool.
unfortunately, this true proportion of mrna, denoted mi, is clearly the quantity of interest.
<cit>  relates the biased proportion θi to mi with the following simple expression:
where ϕi is the known tag formation probability based on eq.
and the proportions mi are positive and sum to one.
the bias corrected likelihood model is then,
direct maximization of the likelihood with respect to the parameters mi is straightforward due to the fact that the likelihood function is maximized at the same location irrespective of the way the model is parameterized.
consider the observed sample proportion  = ti/ttot.
the fact that the mi are positive and must sum to one along with eq.
force the maximum likelihood estimates of  to satisfy the following equality,
while computation of the maximum likelihood estimates  is relatively straightforward, constructing confidence intervals for m is difficult.
the existence of the normalization term ∑jmjϕj in the denominator of the left hand side of eq.
makes computation of fisher's information matrix taxing, particularly for the large dimensional vectors encountered when working with sage data .
in addition, because the number of categories is generally within an order of magnitude of, or possibly even greater than, the number of tags sampled, most categories have either zero or only a few observations.
these relatively small counts call into question inferences such as chi-squared tests which are based on asymptotic approximations.
as an alternative, we consider a bayesian approach to the problem.
the constraints lead naturally to the assumption of a dirichlet prior on m,
where αi is the parameter describing any prior information we have on the value of mi.
combining this prior distribution with the likelihood function given in eq.
leads to a posterior distribution proportional to the product
<cit>  introduces methods for the direct maximization of this quantity and discusses the choice of a prior and its consequences on the marginal inferences for individual values of mi.
however, a number of difficulties severely limit the range of prior parameters that can be used; i.e.
those priors with appreciable weight relative to the observed sample sizes.
given the bias in the observed data, the gibbs sampling based approaches for inference based on the posterior of the mrna proportions introduced here are more flexible, and require fewer assumptions than analytic techniques.
the methods are also highly computationally efficient.
while the assumption of independent cleaving resulted in a geometric model, it is important to note that alternative models of tag formation can be substituted into the bayesian algorithms without making any further adjustments to the model.
shrinkage estimation in multivariate models
joint estimators of the relative proportions of individual tags in sage experiments based on posterior means from a conjugate dirichlet-multinomial model were explored by  <cit>  and have similar properties to the methods proposed here.
in the case of flat priors , they note that while the posterior mean provide improved average accuracy with respect to sums of squared errors across all categories, estimates of categories where large counts are observed were shrunk excessively in order to boost the estimated probabilities of cells with few or no counts.
due to the massive number of categories and the extreme imbalance in observed frequencies , the estimates for frequently observed genes tend to consistently underestimate the true proportions while the reverse is true of genes with very low expression proportions.
this observation highlights the main weakness of shrinkage estimators such as the posterior mean.
while they perform better on average across categories, they may perform poorly on particular categories that may be of primary interest.
<cit>  addresses this issue by proposing a mixture of dirichlet-multinomial's but does not discuss the issue of sampling bias.
note that after aggregation, the large reduction in the number of categories significantly reduces the side effects of shrinkage leading to improved inferences.
these observations will be useful in analyzing the bias corrected models proposed below.
RESULTS
hierarchical models and gibbs sampling
three hierarchical approaches are proposed to model aggregate tag counts while taking into account the tag formation probability.
the dirichlet-poisson-binomial model  extends the method of  <cit> , who modeled tag counts as independent random variables from a poisson-gamma mixture to a multivariate setting that models the joint distribution of tags.
the dirichlet-multinomial-binomial model  views the unobserved true counts of mrna in the population as being derived from a multinomial distribution.
the missing-data model  assumes that if the number of mrna transcripts that were not converted into tags due to incomplete digestion by the anchoring enzyme, denoted r, were known then the data could be modeled as a standard dirichlet-multinomial model.
by proposing a poisson distribution for r, posterior inferences can be made.
detailed descriptions for these algorithms are given in the methods section.
gibbs samplers based on these algorithms all produce posterior samples of m, the vector of relative proportions of each gene in the mrna population.
in addition, dpb produces samples representing n, the hypothetical total number of mrna transcripts in the population while md produces posterior samples of r, the number of unconverted transcripts.
these quantities may be relevant in studies of the number of unique transcripts.
a dirchlet prior distribution was assumed for the proportion vector m in all three algorithms.
two parameter vectors  =  were tested.
the flat prior sets αi =  <dig> for all i while the tub prior sets αi = 1/l.
the flat prior effectively assumes that a hypothetical previous experiment produced  <dig> tag for each category.
as the name suggests, the flat prior represents a uniform density over the parameter space.
the tub prior assumes that the prior information is equivalent to that of a single tag.
the density of this prior takes a bowl or tub shape.
the preponderance of mass along the edges forces the estimates of mi to be small unless significant counts are observed for that gene.
to analyze performance, we applied each algorithm to data collected from the log-phase growth of s. cerevisiae and analyzed by  <cit> .
here, l =  <dig>   <dig> is the number of genes in the sample.
for each of the three methods, we compared posterior means from the gibbs sampler to unadjusted multinomial mle's, corrected mle's based on eq.
, and the analytically derived joint posterior mode given in  <cit> .
in addition, 95% marginal posterior intervals are also computed for each gene from the gibbs sampler output.
figure  <dig> presents results from dpb analysis.
both estimates and intervals generated using the flat prior are plotted versus rank for the  <dig> genes with largest tag counts.
equations  and  dictate that the corrected mle and the posterior mode coincide perfectly in this case.
both posterior means and posterior modes based on the dpb bias correction deviate as much as 40% from the unadjusted mle for many of the most frequently tagged genes indicating the importance of the correction.
posterior modes and means follow nearly identical trends but with the modal estimates being uniformly larger.
figure  <dig> summarizes results when using the tub prior in the dpb case.
here one can see the ithcorrected mle and posterior mean now coincide up to simulation precision .
as with the flat prior, all estimators generated using the tub prior deviate significantly from the unadjusted mle's.
both figures  <dig> and  <dig> show that modal estimates typically lie slightly above the  <dig> % bound of the simulated marginal posterior distribution.
the reader might wonder why the multivariate mode and mean are so far apart.
the explanation lies in the fact that posterior distribution is spread across a a large parameter space consisting of many thousands of dimensions.
while the posterior distribution has a skewed density with a very slight slope, it is essentially flat over most of the parameter space.
while the mode is located at the most likely point, an extreme point where only a fraction of the categories have non-zero tags, the mass spread across the rest of the space pulls the mean away from this point leading to the observed separation.
this is similar to the effect of skewness in a one dimensional distribution which pulls the mean away from the mode.
graphical displays based on the results of the dmb and md algorithms had very similar appearances and so are not shown.
however, meaningful differences in posterior means generated by the methods and across the two priors do exist and are discussed in the following section.
in order to ensure that autocorrelation did not adversely effect parameter estimates, convergence analysis was performed.
posterior samples for the proportion m <dig>  which corresponds to the open reading frame yal003w, were investigated for all approaches.
this gene was observed in  <dig> of the  <dig>  tags and was selected randomly among the set of genes with medium to large tag counts.
autocorrelations of both sequences were tested at various lags.
results given in table  <dig> suggest that the flat prior may generate slightly larger correlations than the tub but in both cases autocorrelations are very low between samples with lags as small as  <dig>
autocorrelation estimates for three proposed algorithms across posterior samples of the mrna proportion at the the open reading frame yal003w in s. cerevisiae log-phase data.
autocorrelation estimates for mrna population size n in the dpb algorithm and number of unconverted transcripts r in the md algorithm.
in terms of computational efficiency, with l =  <dig>   <dig>   <dig> million samples from the dpb algorithm required  <dig>  and  <dig>  minutes with the flat and tub priors respectively.
the dmb was slightly slower taking  <dig>  and  <dig>  minutes to run.
finally, the md method completed  <dig> million samples in  <dig>  and  <dig>  minutes, respectively.
in order to avoid any effects due to autocorrelation our real data experiments stored only the last of every hundred samples for analysis.
inference for means was based upon the final  <dig> of these stored samples drawn after an extensive burn-in period of  <dig>  simulations.
because the vector m is usually the quantity of interest, far fewer samples would be needed in practice.
this would significantly reduce the computation times for the algorithm.
a simulation study
although all methods produce similar trends in their estimates, quantitative differences do exist.
it is also important to quantify the effects of the prior distributions on the analysis.
because the inferential quality is often more relevant in scientific investigations of relative and differential expression than point estimates, we attempt to contrast inferential quality of the three methods by examining marginal 95% posterior intervals constructed for mi from simulated libraries.
the proportion of times the computed interval covers the true known simulation proportion mi is then recorded.
accuracy is assessed by how close the observed coverage percentage is to the desired level of 95%.
the average interval length of the 95% posterior intervals was considered as a second criteria of inferential accuracy but no meaningful differences in average interval length were observed among the three methods.
coverage percentages for 95% posterior intervals for mcmc methods.
percentages represent the average number of intervals out of  <dig> which covered the true proportion mi over the given range of mi values.
standard errors of proportions reported are below  <dig> % and typically close to 1%.
not surprisingly, shrinkage effects become more pronounced as the number of genes l increases further degrading the quality of the estimates.
this shrinkage is most acute for dmb and md methods applied to larger library sizes, which perform poorly for the segments with the.
as noted above, the flat prior works well across a wide range of values but performs poorly in the extremes.
figure  <dig> displays this effect for a subset of  <dig> genes.
alternatively the tub prior works well for the majority of genes where the flat prior fails.
for example, the tub works well in large mi categories because it shrinks the estimates very little.
its weakness; it adds almost no mass to categories with small positive values of mi and, if no counts are observed, intervals are produced with both endpoints near  <dig>  figures  <dig> and  <dig> provide detailed views from a sample of  <dig> simulated libraries for a large and small mi gene.
focusing on the left panel of figure  <dig>  the shrinkage effect is evidenced by the fact that posterior mean estimates of mi systematically underestimate its true value.
this leads to upper interval endpoints being systematically too small.
the tub prior on the right corrects this phenomenon well.
on the left of figure  <dig> we see posterior mean estimates of mi being drawn well above their actual values by the reversing of the shrinkage for small mi categories.
this leads to reasonably accurate coverage.
in the right panel we see that tub prior fails miserably in this case.
intervals for genes with  <dig> counts collapse close to zero systematically underestimating the true mi values.
because these methods are very efficient to compute, both methods can be applied to analyze a single library with the tub prior being relied upon for larger count categories while the flat prior is used for the remaining categories.
differential expression in sage between libraries
tag based transcriptome sequencing such as sage is most commonly used to identify differential expression of genes across groups of libraries.
it is important to point out the consequences of differences in tag formation probabilities in such studies along with the advantages of the proposed methods in this context.
because the anchoring enzyme ae cutting probability p is experiment dependent  <cit> , it is sensible for the tag formation mechanism to be taken into account when evaluating differential expression across experiments.
indeed, this should lead to increased power and accuracy in such studies.
to briefly explore the consequences of variation in cutting probability, consider two experiments a and b. letting pa and mia represent the cutting probability and relative proportion of mrna for the ith tag in experiment a and similarly for b, the theoretical odds ratio for the tag proportions across libraries a and b is
where  is the tag formation probability for tag j in experiment a, Ω is the true odds ratio, and kj is the number of ae positions upstream from which the tag was derived.
the approximation holds when both pa and pb are large  and the ith expression level mi is small in an absolute sense for both experiments, conditions that are almost universally met in sage studies.
if only data derived from the 3' most tagging site of each gene is analyzed so that k =  <dig>  and ϕja = pa then the odds simply reduces to )/), the exact odds ratio desired.
however, suppose all tags are allowed and the experimenter tests differential expression for the ith tag, which derives from an ae site  <dig> position upstream from the 3' end.
the tag formation probability now takes the form p. assuming cutting probabilities pa = .
<dig> and pb = .
<dig> we see that our estimate of the odds ratio is overestimated by a factor of ~ <dig>  potentially leading to a false positive diagnosis.
in fact, this phenomenon could result in a sizeable proportion of both false positive and false negative conclusions.
if sage analysis were restricted to 3' most ae sites and all data from upstream tags were eliminated from the library, ϕi would be constant and sampling bias would be eliminated by analysis using odds ratios .
once again, the caveat to this approach is that a significant proportion of the observed tags may be discarded.
the need for the adjustments we propose stems from an effort to accurately combine multiple tags arising from the same mrna and offers the ability to utilize a much larger fraction of the data collected leading to more power for differential testing.
monte-carlo sampling mechanisms such as dpb can easily be applied to evaluate differential expression across groups of libraries.
in the case of comparison for gene i across two libraries, individual monte carlo samples from each library can be differenced, mi <dig> - mi <dig>  and the distribution of these differences used for inference; see .
as in a paired t-test, if the bulk of differences fall far from  <dig> we may conclude that differential expression exists.
alternatively, it may be more appropriate to form odds ratios based on the sampled quantities and see if the bulk of these are far from  <dig>
the monte-carlo sampling nature of the algorithms presented here and their computational efficiency leads to several advantages over recent methods  <cit> .
first, and foremost the ability to use a larger fraction of the data collected while simultaneously correcting bias is critical to improving power.
in addition, mcmc samples allow both marginal and joint testing for differential expression across all genes using a single sampling run for each library involved.
joint testing of differential expression is possible because posterior samples represent the joint posterior distribution of the proportions m. importantly, our method requires no asymptotic conditions to be met in order to guarantee the validity of the results.
CONCLUSIONS
this work focuses on a general biasing mechanism that may have a significant impact on the interpretation of libraries generated by sage and closely related methods.
beyond drawing attention to the consequences of incomplete digestion, two main contributions are presented.
first, we derive and analyze three bayesian algorithms that provide corrected posterior inference for relative proportions of genes or tags in the overall population.
second, through calculations we deduce the consequences of tag formation bias on testing of differential expression and discuss how the algorithms derived here can be used to correctly assess differential expression.
this work compliments the earlier work in  <cit>  by providing a flexible and efficient methodology which extends the range of the prior distributions that are available.
it also allows researchers access to interval estimates in order to assess uncertainty in the estimated quantities.
the proposed methods are appropriate when incomplete digestion is the main source of multiple tag generation.
they may be applied to both tag level data and to inference at the gene level based on aggregated counts.
if possible, aggregating tags will provide better inferences due to larger available counts.
of the three approaches, the dpb method provided the most accurate confidence intervals over the widest range of relative proportions and we recommend its exclusive use.
in addition, the choice of prior also plays an important role in the effectiveness of the algorithm.
of the two priors tested, the flat prior was effective over a wide range of categories but resulted in excessive shrinkage for the most abundant categories.
for a fixed number of sampled tags, this shrinkage effect becomes more severe as one increases the number of categories.
hence, aggregating tags provides a second advantage since it greatly reduces the number of categories used in an analysis.
the proposed tub prior provides little shrinkage and posterior intervals based on this prior are very accurate for abundant categories.
in practice, both priors should be applied and inference should be based on the samples with the larger average.
we note that the sampling algorithms presented are independent of the model for tag formation probability and models other than the geometric model discussed can easily be substituted.
it should also be possible to extend the proposed methods using a mixture methodology as suggested by  <cit> .
correction of tag formation bias may play a significant role in improving power in differential expression studies based on sage or related tag libraries.
to ensure valid testing, one must either include bias correction or eliminate non-3' tag counts from the analysis.
furthermore, the ability to effectively test for changes in groups of genes  is now possible due to the multivariate nature of the posterior distribution.
this advantage may become significantly more important as more elaborate studies and techniques are developed in the future.
due to the overall simplicity of genome architecture, alternative splicing is generally not a problem when analyzing gene expression in s. cerevisiae and other microorganisms.
however, alternative splicing is widespread in most multicellular organisms  <cit>  and greatly complicates transcriptome analysis.
while the adjustments we propose do not solve the interpretation problems caused by alternative splicing, properly identifying multiple tags from common mrna is a necessary first step in doing so.
a logical next step would be to expand the tag formation model underlying this work by allowing a single gene to produce multiple mrna isoforms instead of a single type.
one approach would be to extend our model in a manner similar to the approach developed for analyzing mrna sequence data developed by  <cit> .
extending our model in such a way should allow us to make inferences about the distribution of mrna isoforms produced by a gene based on the distribution of tags experimentally observed.
conceptually, the problems caused by alternative splicing are similar to those caused by tag ambiguity between different genes.
for example, in the case of alternative splicing, the challenge is in determining the mrna isoform from which each tag originates.
similarly, in the case of ambiguous tags, the challenge is in determining the gene from which each tag originates.
thus it is plausible that data interpretation problems caused by both alternative splicing and ambiguous tags could be overcome by extending our model to include multiple sources of a single tag.

 tag-based transcriptome sequencing libraries consist of a collection of short sequences of dna called tags along with tabulated counts of the number of times each tag is observed in a sample. these observed tag counts represent a sample from a much larger pool of mrna tags in a tissue or organism. in the past, serial analysis of gene expression was the most commonly used technology for generating libraries of tag counts. sage libraries have been used to address a number of biological questions including: estimating transcriptome size, and estimating the density of relative expression levels. most frequently, sage was used to assess differential expression across cells from different tissues or strains, or cells grown under different experimental conditions. next generation methods such as digital gene expression tag profiling now provide a more efficient method to generate tag libraries and are growing in popularity. libraries based on dge are used to address the same questions and provide much larger numbers of tags leading to increased statistical power. however, the close similarities between dge and sage, particularly the use of restriction enzymes, lead both techniques to share the same inherent biases. as has been repeatedly noted in both sage and dge studies e.g., the mrna from a single gene can lead to the production of more than one type of tag within the same library. observing multiple tags from a single gene complicates data analysis and has a number of potential causes including alternative mrna splicing and incomplete cdna digestion. incomplete digestion typically results in a number of different tags being formed from the same mrna transcript. alternative splicing can lead to multiple mrna isoforms within a sample. when combined with incomplete digestion, each of these isoforms may give rise to a different tag or group of tags. while it is intuitively sensible to analyze the prevalence of an mrna by aggregating multiple tags derived from it, such tags are often analyzed separately. at best, this approach diminishes the statistical power that is available when testing differences in mrna abundances and reduces the accuracy of inferences based on asymptotic statistical methods. more seriously, such comparisons can result in dramatic biases leading to erroneous conclusions in testing differential expression due to varying levels of incomplete digestion across libraries. the current work extends an earlier analysis based on sage libraries generated using the yeast saccharomyces cerevisiae. alternative splicing is generally rare in microorganisms such as s. cerevisiae so that the tag multiplicity observed in these libraries is primarily attributed to incomplete cdna digestion. the first objective of the current work is to provide a methodology that allows multiple tags, which arise due to incomplete digestion, to be combined and used to infer expression levels of the underlying mrna transcripts. aggregation is complicated by the fact that certain tags are more likely to form than others due to processing methods. unless the tag formation probabilities are accounted for, inferences based on aggregate tag counts may be highly biased. previously, we have addressed the problem of estimation bias by directly modeling the variation in tag formation and deriving bias corrected maximum likelihood and posterior mode estimators of the composition of the mrna population. however, a number of assumptions and difficulties limit the application of these techniques. we address these in the present study. three new hierarchical models for sage data are introduced that directly correct for the bias by embedding a model of tag formation within the probability distribution describing the data. corresponding gibbs sampling algorithms suitable for model fitting and inference are also derived. these methods are computationally efficient and allow a more general class of prior distributions than earlier analytic approaches. the results are also applicable for both aggregate and standard tag counts. we use simulation studies and analysis of publicly available data from s. cerevisiae in order to contrast the accuracy of our sampling methods and analyze the effects of priors. the approach we present should serve as a building block for a more comprehensive method that takes into account the effects of both incomplete digestion and alternative splicing. a second major objective of the current work is to explore the general effects of tag formation bias on differential expression. we carefully demonstrate that, in the absence of proper adjustments, variation in tag formation probability across experiments can lead to both false positive and false negative conclusions regarding the differential expression of a gene. a brief description of the implementation of proposed algorithms in differential expression analysis is also provided. while the approaches presented here require more work in processing and preparation than studies that focus on individual tags, they may offer a significant advantage with respect to statistical power in studies of differential expression. a number of statistical methodologies have been developed to analyze sage libraries. for studies that focus on relative expression levels, individual tags are viewed as possible outcomes from a multinomial sampling model. offered an improvement on this by directly applying a bayesian multinomial-mixture dirichlet model to the observed vector of tag counts which improves accuracy in estimation of tag frequency. marginal analysis may be more useful in studies of transcriptome size and was developed by and. statistical methods to assess differential expression across a number of samples were discussed by. these analyses typically test the differential expression of each tag across libraries independently; see. the sage methodology the goal of a sage experiment is to sample the mrna population within a group of cells. broadly speaking, the sage method generates a set of short sequence-based cdna tags from the mrna population. initially, a pool of mrna is extracted from a group of cells. the unstable single stranded mrna is reverse transcribed into a double stranded dna copy using a modified primer that allows the cdna to be bound to a streptavidin bead. using restriction enzymes, the bead-bound portion of cdna is cut into small'tags' that are concatenated together into longer cdna molecules referred to as concatemers. these concatemers are amplified using pcr and then sequenced. a cleavage motif of the anchoring enzyme allows the start and stop points of tags to be identified within a concatemer. each time a tag from a specific gene is observed in the sequence data, it contributes a count to the library. the dge technology follows a similar procedure in the first several steps. however, tags are neither concatenated nor cloned but sequenced immediately. in both technologies the restriction enzymes used to create the tags cut at very specific sites within the cdna. for example, the restriction enzyme used by cleaves only at the four nucleotide motif catg. thus, tags can only be created at specific points within a cdna. because the site where the cdna is cut serves as an'anchor' between tags, the restriction enzyme used is often referred to as the anchoring enzyme. we refer to the specific points in the cdna cleaved by the anchoring enzyme as ae sites. while some genes will have no ae sites, many genes have multiple sites at which the ae could cleave the cdna. given that the ae is expected to act in a site-independent manner, a single cdna molecule can be cut by the ae in multiple places. however, because only the fragment of cdna that is attached to a streptavidin bead is retained during the experimental process, the site closest to the bead that is actually cleaved is the only site that can lead to an observable tag. if the ae worked with 100% efficiency, then each mrna could only lead to one observable tag. however, the cutting efficacy of the ae is always less than 100%, sometimes referred to as incomplete digestion, and as a result multiple tags may be generated from the cdna of a single gene. the critical point is that differences in the number of ae sites between genes result in different probabilities of tag formation. the most extreme cases are genes whose mrna transcript lacks any ae site. such genes are impossible to observe using the sage methodology. for genes with ae sites, differences in digestion efficiencies lead to differences in tag formation probabilities between experiments. a further complication is that many of these sites lead to non-unique or'ambiguous' tags that cannot be readily assigned to a particular gene. ambiguous tags are'uninformative' using current technologies; tag-to-gene mapping is discussed at length in. the tags generated in the experiment analyzed below were bp long leading to a number of ambiguous tags. experimental advances, for example,'supersage' techniques can lead to tags up to bp long significantly diminishing this problem. modeling tag formation probabilities the previous discussion leads us to postulate that genes with larger numbers of ae sites have a increased probability of tag formation and further that there is large variation in the probability of tag formation within a gene. below we present a quantitative model for tag formation probabilities within a gene that is key to later bias correction methods. because only cdna that is attached to a streptavidin bead is retained during the experimental process, tags are created from the 3' most ae site that is actually cleaved. let ki be the number of restriction enzyme sites which may be cleaved by the ae on mrna generated from gene i, let p be the probability the ae will cleave a site, and assume that cleavage probability is independent between sites and does not vary between genes. if we label sites through ki starting at the 3' most site, the probability of generating a tag through ae cleavage at site j ∈ {..., ki} follows a geometric distribution, j-1p. this corresponds to the probability of no cleaving in sites to j - and a cleaving at site j. note that this probability is independent of what happens at the ae sites 5' from the jth site. the fact that the expected distribution of tags varies with ae cleavage probability p can be used to estimate p from the library of tag counts; see. it follows that the probability of generating any of the possible tags from the ith category of mrna follows a geometric distribution where the sum is over all non-ambiguous tagging sites. we refer to this quantity as the tag formation probability. for purposes of our analysis, we will assume that, because each of the thousands of genes provides an independent estimate, the value of p is known to a very high precision. the variation in the tag formation probability between genes stems from variation in the number of ae sites as well as the number of ambiguous tagging sites. given p, the only uncertainty in the estimate of ϕi results from ambiguous tags; this problem may be partially eliminated through techniques that generate longer tags; see. hence, for methodological purposes we reasonably assume the ϕi to be known constants. as a result, the tag pool represents a biased sample of the mrna population of interest being a function of both the known distribution of tag formation probabilities across the genome, and the distribution of these genes within the mrna population itself. traditionally, a gene's frequency in the tag pool was assumed to be an unbiased estimate of its frequency in the mrna population, the population of interest. this equivalence only holds when all genes have the same tag formation probability, a condition that is never met. consequently, genes with greater than average tag formation probabilities will be over-represented while those with lower than the average probability will be under-represented. equation allows us to accurately estimate this sampling bias and correct the inferences made from tag libraries. before the tag counts are considered, we assume that the data is processed to retain only informative tags, i.e. all ambiguous or orphan tags are removed, as is standard practice. let the vector t = represent the aggregated observed tag counts for individual genes. here l is the number of genes that contain at least one unambiguous ae site and ti is the total number of observed tags that can uniquely be attributed to an ae site in gene i..) it is natural to view this vector as a sample from a multinomial distribution with l categories, t ~multi where θ =, θi represents the frequency of tags from gene i in the tag pool, and is the total count of all observed tags. because the variation in tag formation probabilities ϕi distorts the abundances of observed tags, it is clear that θi represents a biased estimate of the true proportion of mrna from gene i in the overall pool. unfortunately, this true proportion of mrna, denoted mi, is clearly the quantity of interest. relates the biased proportion θi to mi with the following simple expression: where ϕi is the known tag formation probability based on eq. and the proportions mi are positive and sum to one. the bias corrected likelihood model is then, direct maximization of the likelihood with respect to the parameters mi is straightforward due to the fact that the likelihood function is maximized at the same location irrespective of the way the model is parameterized. consider the observed sample proportion = ti/ttot. the fact that the mi are positive and must sum to one along with eq. force the maximum likelihood estimates of to satisfy the following equality, while computation of the maximum likelihood estimates is relatively straightforward, constructing confidence intervals for m is difficult. the existence of the normalization term ∑jmjϕj in the denominator of the left hand side of eq. makes computation of fisher's information matrix taxing, particularly for the large dimensional vectors encountered when working with sage data. in addition, because the number of categories is generally within an order of magnitude of, or possibly even greater than, the number of tags sampled, most categories have either zero or only a few observations. these relatively small counts call into question inferences such as chi-squared tests which are based on asymptotic approximations. as an alternative, we consider a bayesian approach to the problem. the constraints lead naturally to the assumption of a dirichlet prior on m, where αi is the parameter describing any prior information we have on the value of mi. combining this prior distribution with the likelihood function given in eq. leads to a posterior distribution proportional to the product introduces methods for the direct maximization of this quantity and discusses the choice of a prior and its consequences on the marginal inferences for individual values of mi. however, a number of difficulties severely limit the range of prior parameters that can be used; i.e. those priors with appreciable weight relative to the observed sample sizes. given the bias in the observed data, the gibbs sampling based approaches for inference based on the posterior of the mrna proportions introduced here are more flexible, and require fewer assumptions than analytic techniques. the methods are also highly computationally efficient. while the assumption of independent cleaving resulted in a geometric model, it is important to note that alternative models of tag formation can be substituted into the bayesian algorithms without making any further adjustments to the model. shrinkage estimation in multivariate models joint estimators of the relative proportions of individual tags in sage experiments based on posterior means from a conjugate dirichlet-multinomial model were explored by and have similar properties to the methods proposed here. in the case of flat priors, they note that while the posterior mean provide improved average accuracy with respect to sums of squared errors across all categories, estimates of categories where large counts are observed were shrunk excessively in order to boost the estimated probabilities of cells with few or no counts. due to the massive number of categories and the extreme imbalance in observed frequencies, the estimates for frequently observed genes tend to consistently underestimate the true proportions while the reverse is true of genes with very low expression proportions. this observation highlights the main weakness of shrinkage estimators such as the posterior mean. while they perform better on average across categories, they may perform poorly on particular categories that may be of primary interest. addresses this issue by proposing a mixture of dirichlet-multinomial's but does not discuss the issue of sampling bias. note that after aggregation, the large reduction in the number of categories significantly reduces the side effects of shrinkage leading to improved inferences. these observations will be useful in analyzing the bias corrected models proposed below.  hierarchical models and gibbs sampling three hierarchical approaches are proposed to model aggregate tag counts while taking into account the tag formation probability. the dirichlet-poisson-binomial model extends the method of, who modeled tag counts as independent random variables from a poisson-gamma mixture to a multivariate setting that models the joint distribution of tags. the dirichlet-multinomial-binomial model views the unobserved true counts of mrna in the population as being derived from a multinomial distribution. the missing-data model assumes that if the number of mrna transcripts that were not converted into tags due to incomplete digestion by the anchoring enzyme, denoted r, were known then the data could be modeled as a standard dirichlet-multinomial model. by proposing a poisson distribution for r, posterior inferences can be made. detailed descriptions for these algorithms are given in the methods section. gibbs samplers based on these algorithms all produce posterior samples of m, the vector of relative proportions of each gene in the mrna population. in addition, dpb produces samples representing n, the hypothetical total number of mrna transcripts in the population while md produces posterior samples of r, the number of unconverted transcripts. these quantities may be relevant in studies of the number of unique transcripts. a dirchlet prior distribution was assumed for the proportion vector m in all three algorithms. two parameter vectors = were tested. the flat prior sets αi = for all i while the tub prior sets αi = 1/l. the flat prior effectively assumes that a hypothetical previous experiment produced tag for each category. as the name suggests, the flat prior represents a uniform density over the parameter space. the tub prior assumes that the prior information is equivalent to that of a single tag. the density of this prior takes a bowl or tub shape. the preponderance of mass along the edges forces the estimates of mi to be small unless significant counts are observed for that gene. to analyze performance, we applied each algorithm to data collected from the log-phase growth of s. cerevisiae and analyzed by. here, l = is the number of genes in the sample. for each of the three methods, we compared posterior means from the gibbs sampler to unadjusted multinomial mle's, corrected mle's based on eq., and the analytically derived joint posterior mode given in. in addition, 95% marginal posterior intervals are also computed for each gene from the gibbs sampler output. figure presents results from dpb analysis. both estimates and intervals generated using the flat prior are plotted versus rank for the genes with largest tag counts. equations and dictate that the corrected mle and the posterior mode coincide perfectly in this case. both posterior means and posterior modes based on the dpb bias correction deviate as much as 40% from the unadjusted mle for many of the most frequently tagged genes indicating the importance of the correction. posterior modes and means follow nearly identical trends but with the modal estimates being uniformly larger. figure summarizes results when using the tub prior in the dpb case. here one can see the ithcorrected mle and posterior mean now coincide up to simulation precision. as with the flat prior, all estimators generated using the tub prior deviate significantly from the unadjusted mle's. both figures and show that modal estimates typically lie slightly above the % bound of the simulated marginal posterior distribution. the reader might wonder why the multivariate mode and mean are so far apart. the explanation lies in the fact that posterior distribution is spread across a a large parameter space consisting of many thousands of dimensions. while the posterior distribution has a skewed density with a very slight slope, it is essentially flat over most of the parameter space. while the mode is located at the most likely point, an extreme point where only a fraction of the categories have non-zero tags, the mass spread across the rest of the space pulls the mean away from this point leading to the observed separation. this is similar to the effect of skewness in a one dimensional distribution which pulls the mean away from the mode. graphical displays based on the results of the dmb and md algorithms had very similar appearances and so are not shown. however, meaningful differences in posterior means generated by the methods and across the two priors do exist and are discussed in the following section. in order to ensure that autocorrelation did not adversely effect parameter estimates, convergence analysis was performed. posterior samples for the proportion m which corresponds to the open reading frame yal003w, were investigated for all approaches. this gene was observed in of the tags and was selected randomly among the set of genes with medium to large tag counts. autocorrelations of both sequences were tested at various lags. results given in table suggest that the flat prior may generate slightly larger correlations than the tub but in both cases autocorrelations are very low between samples with lags as small as autocorrelation estimates for three proposed algorithms across posterior samples of the mrna proportion at the the open reading frame yal003w in s. cerevisiae log-phase data. autocorrelation estimates for mrna population size n in the dpb algorithm and number of unconverted transcripts r in the md algorithm. in terms of computational efficiency, with l = million samples from the dpb algorithm required and minutes with the flat and tub priors respectively. the dmb was slightly slower taking and minutes to run. finally, the md method completed million samples in and minutes, respectively. in order to avoid any effects due to autocorrelation our real data experiments stored only the last of every hundred samples for analysis. inference for means was based upon the final of these stored samples drawn after an extensive burn-in period of simulations. because the vector m is usually the quantity of interest, far fewer samples would be needed in practice. this would significantly reduce the computation times for the algorithm. a simulation study although all methods produce similar trends in their estimates, quantitative differences do exist. it is also important to quantify the effects of the prior distributions on the analysis. because the inferential quality is often more relevant in scientific investigations of relative and differential expression than point estimates, we attempt to contrast inferential quality of the three methods by examining marginal 95% posterior intervals constructed for mi from simulated libraries. the proportion of times the computed interval covers the true known simulation proportion mi is then recorded. accuracy is assessed by how close the observed coverage percentage is to the desired level of 95%. the average interval length of the 95% posterior intervals was considered as a second criteria of inferential accuracy but no meaningful differences in average interval length were observed among the three methods. coverage percentages for 95% posterior intervals for mcmc methods. percentages represent the average number of intervals out of which covered the true proportion mi over the given range of mi values. standard errors of proportions reported are below % and typically close to 1%. not surprisingly, shrinkage effects become more pronounced as the number of genes l increases further degrading the quality of the estimates. this shrinkage is most acute for dmb and md methods applied to larger library sizes, which perform poorly for the segments with the. as noted above, the flat prior works well across a wide range of values but performs poorly in the extremes. figure displays this effect for a subset of genes. alternatively the tub prior works well for the majority of genes where the flat prior fails. for example, the tub works well in large mi categories because it shrinks the estimates very little. its weakness; it adds almost no mass to categories with small positive values of mi and, if no counts are observed, intervals are produced with both endpoints near figures and provide detailed views from a sample of simulated libraries for a large and small mi gene. focusing on the left panel of figure the shrinkage effect is evidenced by the fact that posterior mean estimates of mi systematically underestimate its true value. this leads to upper interval endpoints being systematically too small. the tub prior on the right corrects this phenomenon well. on the left of figure we see posterior mean estimates of mi being drawn well above their actual values by the reversing of the shrinkage for small mi categories. this leads to reasonably accurate coverage. in the right panel we see that tub prior fails miserably in this case. intervals for genes with counts collapse close to zero systematically underestimating the true mi values. because these methods are very efficient to compute, both methods can be applied to analyze a single library with the tub prior being relied upon for larger count categories while the flat prior is used for the remaining categories. differential expression in sage between libraries tag based transcriptome sequencing such as sage is most commonly used to identify differential expression of genes across groups of libraries. it is important to point out the consequences of differences in tag formation probabilities in such studies along with the advantages of the proposed methods in this context. because the anchoring enzyme ae cutting probability p is experiment dependent, it is sensible for the tag formation mechanism to be taken into account when evaluating differential expression across experiments. indeed, this should lead to increased power and accuracy in such studies. to briefly explore the consequences of variation in cutting probability, consider two experiments a and b. letting pa and mia represent the cutting probability and relative proportion of mrna for the ith tag in experiment a and similarly for b, the theoretical odds ratio for the tag proportions across libraries a and b is where is the tag formation probability for tag j in experiment a, Ω is the true odds ratio, and kj is the number of ae positions upstream from which the tag was derived. the approximation holds when both pa and pb are large and the ith expression level mi is small in an absolute sense for both experiments, conditions that are almost universally met in sage studies. if only data derived from the 3' most tagging site of each gene is analyzed so that k = and ϕja = pa then the odds simply reduces to )/), the exact odds ratio desired. however, suppose all tags are allowed and the experimenter tests differential expression for the ith tag, which derives from an ae site position upstream from the 3' end. the tag formation probability now takes the form p. assuming cutting probabilities pa =. and pb =. we see that our estimate of the odds ratio is overestimated by a factor of ~ potentially leading to a false positive diagnosis. in fact, this phenomenon could result in a sizeable proportion of both false positive and false negative conclusions. if sage analysis were restricted to 3' most ae sites and all data from upstream tags were eliminated from the library, ϕi would be constant and sampling bias would be eliminated by analysis using odds ratios. once again, the caveat to this approach is that a significant proportion of the observed tags may be discarded. the need for the adjustments we propose stems from an effort to accurately combine multiple tags arising from the same mrna and offers the ability to utilize a much larger fraction of the data collected leading to more power for differential testing. monte-carlo sampling mechanisms such as dpb can easily be applied to evaluate differential expression across groups of libraries. in the case of comparison for gene i across two libraries, individual monte carlo samples from each library can be differenced, mi - mi and the distribution of these differences used for inference; see. as in a paired t-test, if the bulk of differences fall far from we may conclude that differential expression exists. alternatively, it may be more appropriate to form odds ratios based on the sampled quantities and see if the bulk of these are far from the monte-carlo sampling nature of the algorithms presented here and their computational efficiency leads to several advantages over recent methods. first, and foremost the ability to use a larger fraction of the data collected while simultaneously correcting bias is critical to improving power. in addition, mcmc samples allow both marginal and joint testing for differential expression across all genes using a single sampling run for each library involved. joint testing of differential expression is possible because posterior samples represent the joint posterior distribution of the proportions m. importantly, our method requires no asymptotic conditions to be met in order to guarantee the validity of the results.  this work focuses on a general biasing mechanism that may have a significant impact on the interpretation of libraries generated by sage and closely related methods. beyond drawing attention to the consequences of incomplete digestion, two main contributions are presented. first, we derive and analyze three bayesian algorithms that provide corrected posterior inference for relative proportions of genes or tags in the overall population. second, through calculations we deduce the consequences of tag formation bias on testing of differential expression and discuss how the algorithms derived here can be used to correctly assess differential expression. this work compliments the earlier work in by providing a flexible and efficient methodology which extends the range of the prior distributions that are available. it also allows researchers access to interval estimates in order to assess uncertainty in the estimated quantities. the proposed methods are appropriate when incomplete digestion is the main source of multiple tag generation. they may be applied to both tag level data and to inference at the gene level based on aggregated counts. if possible, aggregating tags will provide better inferences due to larger available counts. of the three approaches, the dpb method provided the most accurate confidence intervals over the widest range of relative proportions and we recommend its exclusive use. in addition, the choice of prior also plays an important role in the effectiveness of the algorithm. of the two priors tested, the flat prior was effective over a wide range of categories but resulted in excessive shrinkage for the most abundant categories. for a fixed number of sampled tags, this shrinkage effect becomes more severe as one increases the number of categories. hence, aggregating tags provides a second advantage since it greatly reduces the number of categories used in an analysis. the proposed tub prior provides little shrinkage and posterior intervals based on this prior are very accurate for abundant categories. in practice, both priors should be applied and inference should be based on the samples with the larger average. we note that the sampling algorithms presented are independent of the model for tag formation probability and models other than the geometric model discussed can easily be substituted. it should also be possible to extend the proposed methods using a mixture methodology as suggested by. correction of tag formation bias may play a significant role in improving power in differential expression studies based on sage or related tag libraries. to ensure valid testing, one must either include bias correction or eliminate non-3' tag counts from the analysis. furthermore, the ability to effectively test for changes in groups of genes is now possible due to the multivariate nature of the posterior distribution. this advantage may become significantly more important as more elaborate studies and techniques are developed in the future. due to the overall simplicity of genome architecture, alternative splicing is generally not a problem when analyzing gene expression in s. cerevisiae and other microorganisms. however, alternative splicing is widespread in most multicellular organisms and greatly complicates transcriptome analysis. while the adjustments we propose do not solve the interpretation problems caused by alternative splicing, properly identifying multiple tags from common mrna is a necessary first step in doing so. a logical next step would be to expand the tag formation model underlying this work by allowing a single gene to produce multiple mrna isoforms instead of a single type. one approach would be to extend our model in a manner similar to the approach developed for analyzing mrna sequence data developed by. extending our model in such a way should allow us to make inferences about the distribution of mrna isoforms produced by a gene based on the distribution of tags experimentally observed. conceptually, the problems caused by alternative splicing are similar to those caused by tag ambiguity between different genes. for example, in the case of alternative splicing, the challenge is in determining the mrna isoform from which each tag originates. similarly, in the case of ambiguous tags, the challenge is in determining the gene from which each tag originates. thus it is plausible that data interpretation problems caused by both alternative splicing and ambiguous tags could be overcome by extending our model to include multiple sources of a single tag. 
incomplete digestion during the tag formation process may allow for multiple tags to be generated from a given mrna transcript. three new hierarchical models, which directly embed a model for the variation in tag formation probability, are proposed and their associated bayesian inference algorithms are developed. the consequences of tag formation bias are discussed in the context of testing differential expression. several bayesian inference algorithms that account for tag formation effects are compared with the dpb algorithm providing clear evidence of superior performance. calculations demonstrate the potential for false positive and negative findings due to variation in tag formation probabilities across samples when testing for differential expression. taking this bias into account should allow more of the available data to be used leading to increased statistical power. as a result, the observed tag counts represent a biased sample of the actual transcript pool. the multivariate nature of the approach easily allows both univariate and joint tests of differential expression.
 tag-based techniques, such as sage, are commonly used to sample the mrna pool of an organism's transcriptome. incomplete digestion during the tag formation process may allow for multiple tags to be generated from a given mrna transcript. the probability of forming a tag varies with its relative location. as a result, the observed tag counts represent a biased sample of the actual transcript pool. in sage this bias can be avoided by ignoring all but the 3' most tag but will discard a large fraction of the observed data. taking this bias into account should allow more of the available data to be used leading to increased statistical power.  three new hierarchical models, which directly embed a model for the variation in tag formation probability, are proposed and their associated bayesian inference algorithms are developed. these models may be applied to libraries at both the tag and aggregate level. simulation experiments and analysis of real data are used to contrast the accuracy of the various methods. the consequences of tag formation bias are discussed in the context of testing differential expression. a description is given as to how these algorithms can be applied in that context.  several bayesian inference algorithms that account for tag formation effects are compared with the dpb algorithm providing clear evidence of superior performance. the accuracy of inferences when using a particular non-informative prior is found to depend on the expression level of a given gene. the multivariate nature of the approach easily allows both univariate and joint tests of differential expression. calculations demonstrate the potential for false positive and negative findings due to variation in tag formation probabilities across samples when testing for differential expression. 