Topical Hidden Genome: Discovering Latent Cancer Mutational Topics using a Bayesian Multilevel Context-learning Approach

Consider a training dataset comprising $N$ tumors indexed $n=1,\dots,N$ with $c_{n}\in\{1,\dots,K\}$ denoting the known cancer site/class label for tumor $n$, where $K$ is the total number of cancer sites. Label the $J$ distinct mutations observed in the training dataset as $j=1,\dots,J$; let $x_{nj}$ be a binary indicator of the presence of the $j$-th mutation in the $n$-th tumor ($1\equiv$ presence; $0\equiv$ absence) and let $M_{n}=\sum_{j=1}^{J}x_{nj}$ be the total mutation burden in tumor $n$. Construct the $N\times J$ variant design matrix $X$ by stacking $x_{n}=(x_{n1},\dots,x_{nJ})^{T}$ as rows: $X=((x_{nj}))$. Thus, $e_{n;N}^{T}X=x_{n}^{T}$ where $e_{n;N}$ denotes the $N$-component unit vector whose $n$-th component is $1$ and the remaining components are $0$. We are interested in regressing the cancer site labels $\{c_{n}\}$ on the mutation occurrences $\{x_{nj}\}$, to (a) jointly infer the tissue-site specificities of the individual mutations $\{j\}$, and (b) predict the cancer label $c_{n^{\prime}}$ of a new (test set) tumor $n^{\prime}\notin\{1,\dots,K\}$ through its mutational footprint $x_{n}$ Note that the new tumor $n^{\prime}$ may feature hitherto unobserved mutations $\{j^{\prime}\}\notin\{1,\dots,J\}$, and our regression model needs to appropriately acknowledge this possibility. Below we review the hidden genome model (Chakraborty, Begg, and Shen 2020) that utilizes contextualization of variants for this purpose.

The hidden genome model leverages meta-features, higher level features about the individual variants themselves, embodying mutation contexts. These meta-features are available for all variants, including rare variants observed in the training data and new variants appearing only in the prediction (test) data. Genomic meta-features are commonly categorical, but they may also be continuous. Examples of categorical meta-features include the gene itself, single base substitution types (SBS; 96 possible types), and topological position on the chromosome (e.g., location among $\sim 3000$ 1-megabase length windows). Examples of continuous meta-features include epigenomic feature scores such as histone marks and chromatin accessibility. Here we focus primarily on categorical meta-features to aid topic model based dimension reduction. Continuous meta-feature variables, if present, can be directly incorporated without topic modeling into the hidden genome framework following Chakraborty, Martin, et al. (2021).

To aid distinction we use the term meta-feature source to denote a specific meta-feature variable (e.g., gene, chromosome window) and meta-feature categories to denote individual categories of a meta-feature source (e.g., gene KRAS, chromosome 1 window 1 etc.). Individual meta-feature categories from all meta-feature sources are combined to produce say $P$ total categories labeled as $1,2,\dots,P$. Subsequently we create dummy coded (“one-hot” encoded) binary indicators cataloging correspondence between individual variants and meta-feature categories. For variant $j$ let $u_{j}=(u_{j1},\dots,u_{jP})$ be meta-feature indicators with $u_{jp}$ being $1$ if the $j$-th variant is associated with meta-feature category $p$ and $0$ otherwise. Continuous meta-feature values for variant $j$, if present, are subsequently appended to $u_{j}$. Finally we construct the meta-design matrix $U$ by stacking $u_{j}$ as rows: $U=((u_{jp}))$. The key idea behind the hidden genome model is to quantify the tissue-specific effects of individual variants via a meta-feature regression based on the meta-feature design matrix $U$.

We briefly review below the multilevel hidden genome model and its point estimation as suggested in Chakraborty, Begg, and Shen (2020). The likelihood layer of the model considers a multi-logistic regression connecting the cancer type probability $\Pr(c_{n}=k)$ to the observed variant indicators $\{x_{nj}\}$. A subsequent hierarchical meta-feature regression layer connects variant $j$ and cancer site $k$-specific regression coefficient $\beta_{jk}$ to variant $j$ specific meta-features $u_{j}$ as $\beta_{j,k}=u_{j}^{T}\omega_{\bullet,k}+\beta^{0}_{j,k}$ where $\omega_{\bullet,k}=(\omega_{p,k})_{p=1,\dots,P}$ is the vector of meta-feature effects specific to cancer site $k$, and $\beta^{0}_{j,k}$ is the residual effect of variant $j$ in cancer site $k$ that is not explained by the context. Finally, a group lasso prior for $\{\beta_{j,k}^{0}\}$ and $\{\omega_{p,k}\}$ is considered for regularized estimation. The model is succinctly expressed as follows.

The key meta-feature-regression layer (3.3.2) regresses $\{\beta_{jk}\}$ on the contexts $\{u_{j}\}$ effectuating signal condensation in all, including rare, variants. The residual effects $\{\beta_{jk}^{0}\}$ are non-zero and indeed can be realistically be estimated for only a few strongly discriminative non-rare variants. Following computation of the $XU$ matrix arising after substituting (3.3.2) into (3.3.1)), one may therefore perform a feature screening to keep all but a few recurring tissue-specific variants for estimation of their $\beta^{0}_{j,k}$. For scalable point estimation of the model Chakraborty, Begg, and Shen (2020) marginalize out all but the key parameters of interest, viz., ${{\alpha}},{{\beta}}^{0},{{\omega}}$, given the sparsity level parameter $\lambda$ to produce the following log marginal posterior density.

This is a group lasso penalized multi-logistic log-likelihood with intercept ${{\alpha}}$, and regression coefficients ${{\beta}}$ and ${{\omega}}$ with associated predictors ${{x}}_{n}$ and ${{x}}_{n}^{T}{{U}}$ respectively, and penalty parameter $\lambda$. The corresponding (marginal) posterior mode can be efficiently computed using existing software (Friedman, Hastie, and Tibshirani 2010), and consistency of the resulting estimates follow from group-lasso theory. Note however that uncertainty quantification for statistical inference in this framework is infeasible due to voluminousness of the meta-feature categories.

At the outset, we deviate from Chakraborty, Begg, and Shen (2020) by using instead a normalized version of the predictor $\{\tilde{x}_{nj}=x_{nj}/M_{n}\}$ measuring proportions of total mutation burden $M_{n}$ attributable to individual variants $\{j\}$ and producing a normalized version $\tilde{X}$ of the variant design matrix. This normalization acknowledges mutation burden heterogeneity observed in large-scale genomic data. Furthermore, it aids a product matrix $\tilde{X}U$ whose columns catalog analogous proportions of $\{M_{n}\}$ attributable to individual meta-feature categories. Next, we embed into the hidden genome model a hierarchical topic model layer for interpretable dimension reduction of the mutation context categories. For notational simplicity below we restrict to a single, categorical meta-feature source (e.g., only SBS) for the model; a generalization with several categorical and continuous meta-feature sources is described in Section 3.6. The proposed topical hidden genome model is first succinctly expressed as follows.

The posterior distribution for the proposed model is highly intractable. For implementation we propose a doubly data augmented, blocked, independence Metropolis-within-Gibbs algorithm for MCMC sampling from the target posterior. We provide a summary of the algorithm below; a detailed description with notes on initialization is provided in Supplement A. We note that to aid computation, prior to using the MCMC sampler one needs to screen the columns of the training predictor matrix $\tilde{X}$ after computing the product matrices $XU$ and $\tilde{X}U$. This effectively amounts to setting the residual effects $\{\beta_{jk}^{0}\}$ to be zero apriori for all but a few (say $\leq 50$) highly discriminating variants $\{j\}$. Following Chakraborty, Begg, and Shen (2020) we use a computationally efficient mutual information based variable screening for this.

The proposed Gibbs-type MCMC algorithm iteratively draws (a) topic model parameters, and (b) group-lasso multi-logistic parameters in two blocks each containing a separate data-augmentation step. In block (a) latent Poisson data $\{Z_{nsp}\}$ (Liang and Hoffman 2014) are augmented to aid MCMC sampling of $(H,W)$. To account for the supervised contribution of the multi-logistic layer we propose an independence Metropolis-Hastings step that requires no manual tuning. In block (b) a Gibbs draw is performed for the hierarchical group-lasso (Kyung et al. 2010) multi-logistic parameters with a column-scaled $(\tilde{X},\tilde{X}UW^{-})$ used as predictor matrix (see Chakraborty, Begg, and Shen 2020). Conditional conjugacy for the multi-logistic parameters are achieved via Pólya-Gamma data augmentation (Polson, Scott, and Windle 2013).

To handle multiple categorical meta-feature sources, one may introduce into the model separate, independent topic model layers (3.4.7)-(3.4.8) and a topic-to-observed coefficient connection layer (3.4.3) for each source. Continuous meta-feature sources, if present, would directly produce observed meta-feature regression coefficients $\omega$. All observed meta-feature regression coefficients are then stacked for use in layer (3.4.2). The MCMC algorithm discussed in Section 3.5 requires small modifications to reflect these new layers. With several categorical meta-feature sources one would sequentially update the topic parameters for each source conditional on the “supervised” multi-logistic contributions of the other source topics. The hierarchical group-lasso multi-logistic regression is performed on an expanded column scaled predictor matrix column-stacking $X$, all $\tilde{X}UW^{-}$ matrices arising out of categorical meta-feature sources, and all $XU$ matrices arising out of continuous sources.

A full Bayesian “ensemble” prediction using posterior MCMC draws $\{\alpha^{(t)},\beta^{0,(t)},\omega^{(t)}:t=1,\dots,T\}$ for the multi-logistic parameters is described as follows. Given a “new”/test tumor with normalized variant indicators $\tilde{x}_{\text{new}}$, one first finds variants $\{j\}$ in $\tilde{x}_{\text{new}}$ with non-zero residual effects $\{\beta_{j,k}^{0,(t)}:k=1,\dots,K\}$ for each draw $t=1,\dots,T$. All remaining variants $\{j^{*}\}$, including “novel” variants in test data and screened out variants in training data are assigned $\{\beta_{j^{*},k}^{0,(t)}=0:k=1,\dots,K\}$. The meta-feature values $u_{j*}$ for a “novel” variant $j^{*}$ are cataloged as a new row of $U$ for computation of $\tilde{x}_{\text{new}}^{T}U$. Subsequently, for each draw $t\in\{1,\dots,T\}$ one computes first the linear predictions using the layers (3.4.1) and (3.4.2): $\{\eta_{k}^{(t)}=\alpha_{k}^{(t)}+\tilde{x}_{\text{new}}^{T}\beta^{0,(t)}+\tilde{x}_{\text{new}}^{T}U\omega^{0,(t)}:k=1,\dots,K\}$ for all $K$ cancer sites. These draw-specific linear predictions are then “softmax”-ed and subsequently averaged across draws to compute cancer-site specific posterior predictive probabilities:

Note that unlike a point estimate-based prediction (e.g., MAP prediction) the above acknowledges model estimation uncertainty, and is thus preferred in full Bayesian inference.

The topic related parameters $W$, $H$, and $\{\theta_{sk}\}$ may reveal important biological insights. However, their estimation is challenging due to their non-identifiability caused by the mixture structure of the topic model. Directly using MCMC outputs for these parameters for point/interval estimation may produce misleading results.⁵⁵5Note only direct inference on the (topic) parameters $H$, $W$, and $\{\theta_{sk}\}$ is challenging. Inference on posterior quantities that marginalizes over these parameters, e.g., prediction for new data and estimation of $\{\beta^{0}_{jk}\}$ and $\{\omega_{pk}\}$ is straightforward. Instead, using a point process representation of MCMC draws (Frühwirth-Schnatter 2011), we introduce a simple clustering approach to post-hoc identification and approximate inference for these parameters.

Briefly, we first create a list of all topics (the rows of $W$) drawn across all iterations of the MCMC run, and cluster these pooled rows using $k$-means with an appropriately chosen $k$ (via elbow method in our computations). The key idea is then to use the consequent cluster assignments to relabel MCMC draws for the topics (rows of $W$), topic exposures (columns of $H$) and topic regression coefficients $\{\theta_{sk}\}$. These cluster-relabeled draws are then treated as posterior MCMC draws for “post-hoc identified” parameters and used for MCMC output-analysis as usual. We make a few notes on the use of $k$-means as a clustering algorithm in the current context. First, a pairwise distance based clustering algorithm, e.g., $k$-medoids, might be conceptually preferable to $k$-means given the simplex nature of the topics. However they are not practicable when the number of topics $S$ and/or the number of MCMC iterations $T$ is even moderately large. Second, $k$-means results are sensitive to outliers. Here outliers correspond to posterior draws for the rows of $W$ that are far from any other rows of $W$ drawn in any iteration, and hence may represent configurations with low posterior probabilities. To avoid the influence of these outliers, we considered an approximate nearest neighbors-based filtering prior to $k$-means clustering. Finally, $k$ means algorithm suffer from high computational expense and reduced stability as the number of meta-feature categories (number of columns of $W$) grows. This is particularly challenging for gene and window meta-feature sources with hundreds to thousands of categories. To manage the computation load, prior to $k$-means (and outlier detection) we perform PCA and retain the first few (50 in our computations) principal components.

The above approach resembles the $k$-means finite mixture model identification approach of Frühwirth-Schnatter (2011). However, instead of only resolving label-switching/permutation modes as done in finite mixtures here we use clustering to also possibly combine duplicated topics/modes arising within the same MCMC iteration. Furthermore, while heuristic, Bayesian estimation via “cluster identified” MCMC draws may still be viewed as an approximate solution to a rigorously defined decision theoretic problem; see Stephens (2000).

On each dataset, we perform 10 fold cross-validations to assess predictive performance of the model. The folds are created using stratified (based on cancer sites) random partitions of the PCAWG tumors; these produce 10 different training/test set combinations obtained by pooling each 9 of the 10 folds for training and the remaining fold for testing. On each training set a separate model is fitted and then used for Bayesian (posterior) prediction of the corresponding test set. The predicted class probabilities for all tumors in all folds are subsequently combined. To aid comparison, alongside we also obtain analogous MAP predictions from the hidden genome model (without any topics) of Chakraborty, Begg, and Shen (2020) with similar normalized $\tilde{X}$ as predictors. For each model, we derive one-vs-rest classification probabilities for each of the $K=10$ cancer sites, and compute the associated area under the one-vs-rest precision-recall curve (PR AUC) as a measure of site-specific predictive performance (Saito and Rehmsmeier 2015). PR AUCs reflect class-size imbalances, and are more informative than ROC AUCs for one-vs-rest comparisons obtained from a multi-class classifier. The site-specific PR AUCs are then averaged to produce an overall metric for each model. We consider 10 random replications of this cross-validation, and obtain model specific average PR AUCs across replications. These average PR AUCs are displayed as barplots in Figure 1 with results from the three different subsets plotted along panels. The figure also shows the null baseline PR AUC value for each site, corresponding to a null classifier randomly assigning observations to classes. The null baseline PR AUCs are proportional to the sample size of the corresponding class (see Chakraborty, Martin, et al. (2021)).

We make the following observations from Figure 1. First, the full Bayesian prediction from the proposed model has a better overall-level performance than the MAP prediction from the original model. The improvement is most prominent in targeted gene panel data followed by exome level and whole genome level data. This overall improvement is likely due to the topic model-based dimension reduction which facilitates more stable estimation and prediction. Second, all except one cancer site demonstrate improvement in site-specific prediction in the proposed model. Certain cancer sites, such as lung and pancreas, show consistently better prediction in all three datasets; while for certain other sites, notably skin cancer in targeted gene sequencing, the improvement is noticeable only at specific subsets. Third, for only one site, namely prostate, is there a drop in PR AUC in the proposed approach. This indicates the lack/non-existence of discriminative topics for prostate cancer, which is known to have a rather sporadic mutational landscape. However, the nominal drop in prostate is more than compensated by improvements in the other sites, making the proposed approach superior to the hidden genome model of Chakraborty, Begg, and Shen (2020) as a predictive model.

This section illustrates how to make inference from the proposed model though its applications on the PCAWG datasets. For brevity here we focus on the full whole genome and the exome-subsetted datasets; analogous inferences can also be drawn on the targeted cancer gene panel subsets. Using the proposed MCMC algorithm we fit our model to each dataset similarly as discussed in Section 4.1, except now we use the full datasets instead of cross-validation partitions to train the models. Subsequently in each dataset we compute point (posterior mean and median) and 80% (highest posterior density) interval estimates of topic model parameters and one-vs-rest log odds ratio parameters obtained from regression coefficients for latent meta-feature topics, observed meta-feature categories, and residual variant effects, using the MCMC draws for model parameters.