A hierarchical Bayesian mixture model for inferring the expression state of genes in transcriptomes

A Generative Model.

To introduce our model, it is helpful to imagine using it to generate data. We begin in the upper level of the hierarchical model, which reflects the true expression level of genes in the transcriptome; this is a mixture distribution composed of inactive (blue) and active (red) genes (Fig. 1A). For each gene, we randomly draw the expression state from this mixture distribution: specifically, a gene is inactive (active) with probability proportional to the area under the blue (red) distribution. If the selected expression state is inactive, it will either have zero transcripts (with probability proportional to the blue spike) or nonzero transcripts, in which case its expression level is drawn from the inactive (blue) normal distribution. Conversely, if the selected expression state is active, its expression level will be drawn from the active (red) normal distribution.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

A hierarchical Bayesian mixture model for inferring the expression state (active or inactive) of genes from replicate transcriptomic libraries. We introduce our model by describing how it could be used to simulate transcriptomic libraries. Panel A depicts the true expression state—inactive (blue) or active (red)—and expression level of all genes in the transcriptome. We simulate each gene by randomly sampling from this mixture distribution. Our first four draws include two inactive genes—one with zero expression (gene 1, in the “spike” at left) and one with nonzero expression (gene 2)—and two active genes (3 and 4). Panel B depicts the probability that a gene is not detected (Left) and—given detection—the observed expression level of each gene across libraries (Right). For each simulated gene, we first determine whether it is detected in each library; if a gene is not detected in a given library, it will have an observed expression level of zero (i.e., be assigned to the library-specific spikes at left). If a gene is detected in a given library, its observed expression level will be drawn from a normal distribution (the gene-specific distributions at right) that describes its variation across all libraries in which it is detected. These normal distributions have a mean equal to the true expression level of each gene and a gene-specific variance. Panel C depicts the observed expression level of all genes—with zero transcripts (Left) or nonzero transcripts (Right)—in two replicate libraries. For example, gene 1 was not detected in either library because its true expression is zero (panel A). The observed expression levels of genes 2 to 4 were drawn from their corresponding normal distributions (panel B), resulting in zero transcripts for gene 2 in library b and nonzero transcripts for the remaining genes in both libraries. To generate a complete library with n genes, we repeat the above procedure n times. Like real datasets, transcriptomes simulated under our model have bimodal expression levels, albeit the active and inactive distributions are obscured by library-specific factors (e.g., sequencing depth) and gene-specific factors (e.g., gene length, true expression level, and gene-specific variance). When used as a generative model, we assume the parameter values are known and the data are unknown; conversely, when used as an inference model, we assume the data are known (observed) and the parameter values are unknown (inferred).

Having simulated the true expression level for each gene, we now simulate their observed expression levels (Fig. 1B). For each gene, we first determine whether it is detected in each transcriptomic library. The probability that a gene is detected in a given library depends on its true expression level, its length, and library-specific factors (e.g., sequencing depth). For any library in which the gene is not detected, the observed expression level will be zero (Fig. 1 B, Left). For all libraries in which the gene is detected, we draw its observed expression level from a normal distribution, with a mean equal to its true expression level and a gene-specific variance (Fig. 1 B, Right).

Like empirical datasets, transcriptomic libraries simulated under our model have a characteristic bimodal distribution, with a dominant right mode and a left shoulder (Fig. 1C). We simulate a set of transcriptomic libraries by repeatedly drawing from the gene-specific distributions described above. Biological and technical sources of variation (in the lower level of our hierarchical model) largely obscure the distinct inactive and active distributions of true expression levels (in the upper level of our hierarchical model). Note that the number of genes with zero transcripts may differ among libraries owing to differences in their sequencing depth. Additionally, the rank order of the expression level of genes may vary among libraries owing to variation in the observed expression level of each gene; e.g., an inactive gene may have a higher observed expression level than an active gene in a given libraries (Fig. 1C, arrows).

Model Parameters and Inference.

Our goal is to infer the expression state of each gene from replicate (observed) transcriptomic libraries using our hierarchical Bayesian mixture model. When used as a generative model (as above), we assume that the parameter values are known and the data are unknown. To perform inference under our model, we treat the data as known (observed) and treat the parameter values as unknown. Here, we describe the parameters of the lower and upper levels of the hierarchical model and adopt a Bayesian approach to estimate those parameters from observed transcriptomic data.

Our hierarchical Bayesian model describes the processes that give rise to our observed dataset, which we denote X, that is composed of two or more replicate transcriptomic libraries. The lower level of our model describes the observed expression levels for each gene across all libraries. Specifically, we model the variation in observed expression levels for each gene across libraries with a gene-specific variance parameter. We assume that the variance parameter is inversely related to the (log) expression level, where genes of similar expression levels have similar levels of variation across replicate libraries. Additionally, the lower level includes library-specific parameters that impact the probability that a gene is detected in each library. We represent all of the parameters in the lower level of our model with the container parameter θ1.

The upper level of our hierarchical Bayesian model describes the distribution of true (unobserved) expression levels, which we denote Y. We assume that the true expression levels of genes can be divided into two components: those genes that are actively expressed and those that are not actively expressed. The assignment of each gene to these (in-)active expression-state components is described by parameter zga, where zga=1 indicates that the gene is assigned to the active component and zga=0 indicates that it is assigned to the inactive component. We refer to the assignments of all genes to the (in-)active expression states as za. We further assume that inactive genes can be subdivided into two subcomponents: one with zero expression and another with nonzero expression. Similarly, active genes may be subdivided into one or more subcomponents with distinctly different expression levels (e.g., housekeeping genes may collectively have higher true expression levels relative to other genes). A given model assumes a specific number of active subcomponents (e.g., the model in Fig. 1A has a single active subcomponent); a model with two active subcomponents would have two red distributions. We can specify a set of distinct models with different numbers of active subcomponents and compare their fit to a given dataset (see below). We represent all of the parameters in the upper level of our model—describing true active and inactive distributions—with the container parameter θ2.

We infer the joint posterior probability distribution of the hierarchical model parameters—including the set of parameters describing the expression state of all genes, za—given our observed transcriptomic data, X, by applying Bayes’ theorem:P(za,θ1,θ2,Y∣X)︷posterior distribution=P(X∣Y,θ1)P(θ1)︷lower levelP(Y∣za,θ2)P(za,θ2)︷upper levelP(X)︸marginal likelihood,where the first term in the numerator is the joint probability of the lower level of the hierarchical model given the local model parameters, θ1; the second term is the joint probability of the upper level of the hierarchical model given the local model parameters, θ2; and the denominator is the average probability of the data under the model (the marginal likelihood).

The posterior probability distribution, P(za,θ1,θ2,Y∣X), cannot be calculated analytically because the marginal likelihood, P(X), is impossible to evaluate. Accordingly, we use a numerical algorithm—Markov chain Monte Carlo (MCMC) (22⇓⇓–25)—to approximate the posterior probability distribution. The MCMC algorithm samples parameter values in proportion to their posterior probability. From these MCMC samples, we compute the posterior probability that a given gene is active as the fraction of MCMC samples where zga=1. We validated our MCMC implementation by running it under the prior and by measuring coverage probabilities using simulated data.

Model Checking.

The Bayesian approach for assessing model adequacy is called posterior-predictive assessment (26). This approach is based on the following premise: if our inference model provides an adequate description of the process that gave rise to our observed data, then we should be able to use that model to simulate datasets that resemble our original data. The resemblance between the observed and simulated datasets is quantified using a summary statistic. We use three summary statistics: 1) the upper-level Wasserstein statistic (which measures the discrepancy between the expected and realized true expression levels), 2) the lower-level Wasserstein statistic (which measures the discrepancy between the expected and realized observed expression levels), and 3) the Rumsfeld statistic (which measures the discrepancy between the observed and expected number of undetected genes).

Replicate Libraries Improve Our Ability to Correctly Infer Expression States.

We expect that our ability to correctly infer expression states will depend on the disparity between the true distributions of active and inactive expression levels and the number of replicate libraries. Specifically, we expect the power to increase as we 1) decrease the degree of overlap between the true (in-)active distributions and 2) increase the number of replicate libraries used to estimate the model parameters.

We simulated data with low, moderate, and high levels of overlap between the true active and inactive distributions. For each condition, we simulated datasets comprising two, four, and six replicate libraries. For each unique combination of overlap and library number, we simulated 100 datasets. We measured power by evaluating the posterior probability of the true expression state, averaged across all of the genes in the transcriptome. Our results reveal that our method generally has good power (on average, we inferred the correct expression state for ≈90% of the simulated datasets), which increases with the number of replicate libraries and the disparity between true active and inactive distributions (Fig. 2, Left).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Exploring the power and robustness of our hierarchical Bayesian mixture model to correctly infer the expression state of genes from simulated data. (Left) We explored the power of our method by simulating datasets with two, four, and six replicate libraries under varying degrees of overlap (low, moderate, and high) between the true active (red) and inactive (blue) distributions. The power to infer the true expression state is generally high and increases with the number of replicate libraries and/or the degree of separation between the true (in-)active distributions. (Right) We explored the robustness of our method to model misspecification by simulating datasets with one or two active (red) subcomponents. For simulated datasets with two active subcomponents, we varied their degree of overlap (low, moderate, and high). We analyzed each simulated dataset under two models: a model with one active subcomponent and a model with two active subcomponents. The power of the method to correctly infer expression states is robust to model overspecification: estimates from datasets with one active subcomponent are virtually identical under the correct and overspecified models (leftmost pair of box plots). Similarly, the power of the method is robust to moderate model underspecification: estimates from datasets with two active subcomponents are virtually identical under the correct and underspecified models (two middle pairs of box plots), except when the degree of disparity between the two active subcomponents is extreme (rightmost pair of box plots).

Estimates of Expression State Are Robust to Model Misspecification.

We expect that our ability to correctly infer expression states of genes will be adversely affected when the number of assumed active subcomponents in the hierarchical mixture model differs from the true number of active subcomponents (i.e., when the model is misspecified). The model may include either too many active subcomponents (overspecified) or too few active subcomponents (underspecified).

We simulated datasets with one or two active subcomponents. In the latter scenario, we varied the degree of overlap (low, moderate, and high) between the active subcomponents. For each scenario, we simulated 100 datasets, each with four replicate libraries. For each simulated dataset, we inferred expression states under a model with one or two active subcomponents. When the model was overspecified (i.e., with one true active subcomponent and two assumed active subcomponents), the expression-state estimates were virtually identical to those inferred under the correctly specified model. When the model was underspecified (i.e., with two true active subcomponents and one assumed active subcomponent), the accuracy of expression-state estimates decreased as the disparity between the two true active subcomponents increased. These results indicate that expression-state estimates are robust to overspecification and moderate underspecification but are sensitive to severe underspecification (Fig. 2, Right).

Human-Lung Transcriptomes.

Our first empirical benchmark is a human-lung dataset comprising 427 libraries with 19,154 protein-coding genes sourced from the Genotype-Tissue Expression (GTEx) RNA-seq database (27, 28). We inferred the expression state of each gene using the extensive epigenomic dataset for human-lung tissues from the Roadmap Epigenomics Consortium (27, 29, 30); this dataset includes 15 epigenetic marks that are strongly associated with expression state. Using this epigenetic evidence, we were able to confidently classify the expression state of 11,968 genes (SI Appendix has details); we identified 7,261 active and 4,707 inactive genes (Fig. 3, Left). These expression-state assignments (treated as known) provide an empirical benchmark to assess the power of our method.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Exploring the power of our hierarchical Bayesian mixture model to accurately infer known expression states of genes in the human-lung transcriptome. (Left) Average observed expression levels of genes in the human-lung transcriptome; active (red) and inactive (blue) genes are known from epigenetic marks, providing an empirical benchmark to assess the performance of our method (genes that could not be classified using epigenetic marks are shown in gray; here, genes with log 0 TPM are represented by the gray and blue bars on the left). Note that the expression levels of active (red) and inactive (blue) genes fall into two distinct but overlapping distributions. (Center) We used our model to infer the expression state of all genes from datasets consisting of 2, 4, 8, and 16 randomly selected libraries: we depict estimates for three example genes, where active (red) or inactive (blue) expression states were inferred from a dataset with 4 randomly selected replicate libraries (gray distributions; here, genes with log 0 TPM are not shown for clarity). (Right) We compared our inferred expression states with the known expression states; the power of our method to correctly infer the known expression states is generally high. The use of multiple replicate libraries improves power, and this benefit is realized with only a modest number (four) of replicate libraries. Box plots represent variation in estimates of power across the 10 sets of randomly selected datasets for each number of libraries.

Next, we used our hierarchical Bayesian mixture model to infer the expression state of all 19,154 protein-coding genes. We analyzed data subsets consisting of 2, 4, 8, and 16 randomly selected replicate libraries. For each number of libraries, we sampled 10 independent datasets (e.g., 10 sets of two libraries, 10 sets of four libraries, etc.). We measured power by evaluating the posterior probability of the true expression state for each gene, averaged across all of the genes in the transcriptome. The results of these empirical analyses confirm the findings of our simulation study; the power is generally high (>90% in all cases), and the method performs well with four libraries (Fig. 3, Right).

Drosophila-Testis Transcriptomes.

Our second empirical benchmark is a Drosophila-testis transcriptomic dataset (with four libraries) that we generated for this study. In this experiment, we assessed the power of our method to correctly identify expression states in a challenging empirical setting (i.e., where genes are known to be active in a small number of cells within a tissue, with correspondingly low tissue-wide expression levels). Germline stem cells and several types of somatic cells collectively comprise the stem-cell niche at the tip of the testis (restricted to 20 to 30 cells per testis); we used developmental-genetic evidence to identify 39 active genes in the stem-cell niche (SI Appendix, Table S4 shows a list of studies). Conversely, we identified 119 genes that are known to encode odorant and gustatory receptors that are unlikely to be active in the testis; we therefore classify these genes as inactive.

We used our method to infer the expression state of all genes in the Drosophila-testis dataset. We measured power by evaluating the posterior probability of the true expression state for each gene, averaged across all of the genes in the transcriptome. As previously, the power of our method is generally high, even for genes that are actively expressed in a tiny fraction of cells in the tissue (SI Appendix, Fig. S11). Among genes that are known to be actively expressed in the stem-cell niche, the median inferred posterior probability of being in the active expression state was 0.96, with 37 of the 39 genes inferred to be active (P[active]>0.5). Among olfactory- and gustatory-receptor genes, which we assume are inactive in the testis, the median inferred posterior probability of being in the active expression state was 0.005, with 111 of the 119 genes inferred to be inactive (P[active]<0.5).

Theoretical Power Analysis.

Our analyses of simulated and empirical-benchmark datasets demonstrate that our method generally has high power to infer true/known expression states. Here, we attempt to evaluate the absolute power of our method. To this end, we first establish an upper bound on the power to infer expression states (under a method that requires known expression states) and then compare the power of our method with this reference.

Specifically, we imagine a threshold-based method (i.e., where a gene is inferred to be active if its relative expression level exceeds a fixed threshold value). Unlike actual threshold-based methods (28, 31, 32), this “omniscient” threshold-based method knows the true expression state of each gene. Because this method is aware of the true expression states, it can choose the perfect threshold value that simultaneously maximizes the number of true active genes it infers to be active (the true-positive rate) and minimizes the number of true inactive genes it infers to be active (the false-positive rate).

We first characterize the power of the omniscient threshold method by applying it to the empirical-benchmark datasets (where the expression state of each gene is known). Specifically, we characterize its power by plotting receiver operating characteristic (ROC) curves: for each possible threshold value, we compute the true- and false-positive rates and plot the true-positive rate as a function of the false-positive rate (Fig. 4, orange curves). Note that a method with perfect power would exhibit an L-shaped ROC curve as it would simultaneously achieve a 100% true-positive rate and a 0% false-positive rate. Conveniently, we can compare the power of two methods by comparing their ROC curves.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

The power of our method to infer expression states approaches the practical limit for this inference problem. We used our method to infer expression states of all genes in the two empirical-benchmark datasets, human-lung (Left) and Drosophila-testis (Right) transcriptomes. A gene is assigned to the active expression state if its posterior probability of being active is greater than P. For all possible values of P, we plot the true-positive rate (the fraction of active genes correctly assigned to the active expression state) against the false-positive rate (the fraction of inactive genes incorrectly assigned to the active expression state; blue curves). The resulting ROC curves characterize the discriminatory power of a method (i.e., its ability to distinguish between active and inactive genes). The power of our method (which infers the unknown expression states) is equivalent to that of an omniscient threshold-based method that requires knowledge of the true expression states (orange curves; one for each library). Symbols indicate conventional P thresholds.

Next, we plot ROC curves for our method based on the same empirical-benchmark datasets. Our method infers the posterior probability that each gene is (in-)active. In principle, we could adopt any posterior probability threshold to classify the expression state of each gene. Accordingly, we plot ROC curves by computing the true- and false-positive rates for all posterior probability thresholds between zero and one (Fig. 4, blue curves).

Remarkably, the power of our method is virtually identical to that of the omniscient threshold-based method for both the human-lung and the Drosophila-testis datasets (Fig. 4). These results demonstrate that our method—under the typical inference scenario, where the true expression states of genes are unknown—is able to correctly infer expression states as well as a method that requires a priori knowledge of the true expression states.