HiLDA: a statistical approach to investigate differences in mutational signatures

Hierarchical Bayesian Mixture Model

We introduce a hierarchical latent Dirichlet allocation model (HiLDA) using the following notation, also summarized in Table 1. Let i index the mutational catalog and j the mutation. The nucleotide substitutions are reduced to six possible types (C >A, C >T, C>G, T >A, T>C, T>G) to eliminate redundancy introduced by the complementary strands. Each observed mutation is characterized by a vector, X_i,j describing the nucleotide substitution (e.g., C >T) and a set of genomic features in the neighborhood. Example features include the base(s) 3′ and 5′ of the nucleotide substitution (C, G, A, T), and the transcription strand (+, −). Each observed feature characteristic, x_i,j,l for mutation feature l, takes values in the set {1, 2, …, M_l} (where M_l = 6 for the nucleotide substitution, or 4 for a flanking base, and 2 for the transcription strand).

We assume each mutation belongs to one of K distinct signatures. A specific mutational signature k is defined by an l-tuple of probability vectors, F_k, denoting the relative frequencies of the M_l discrete values for the l features, i.e., a vector f_k,l for the M_l values corresponding to feature l. We let z_i,j denote the unique latent assignment of mutation X_i,j to a particular signature. Then, given the signature to which a mutation belongs, the probability of observing a mutational pattern is calculated as the product of the mutation feature probabilities for that signature. Thus, for signature k we write Pr(X_i,j|z_i,j) = ∏_lf_k,l(x_i,j,l|z_i,j). This assumes independent contributions of each feature to the signature. To model each multinomial distribution of f_k,l, we use a non-informative Dirichlet prior distribution with all concentration parameters equal to one.

The unique personal exposure history of each individual leads to them having a particular (latent) vector, q_i, indicating the resulting contribution of each of the K signatures to that individual’s mutational catalog. These qs are modeled using a Dirichlet distribution with concentration parameters α, i.e., q_i ∼ Dir(α). Extending this model to the two-group setting, we allow the Dirichlet parameters to depend on group, Dir(α^(g_i)), with g_i indexing the group corresponding to the ith catalog (g_i = 1 or 2). The mean mutational exposures, E(q_i), denoted by μ^(g_i), are represented by using the concentration parameters, i.e., μ^(g_i) = α^(g_i)∕∑α^(g_i).

With this extension, we can infer differences in mutational processes between groups of catalogs by testing whether the mean mutational exposures differ between the two sets, i.e., at least one μ_k⁽¹⁾ ≠ μ_k⁽²⁾. The likelihood and prior of the multi-level model is specified as follows,

${\displaystyle x_{i,j,l}|z_{i,j}\sim Multinomial\left({\mathit{f}}_{z_{i,j},l}\right)}$ ${\displaystyle z_{i,j}\sim Multinomial\left({\mathit{q}}_i|g\right)}$ ${\displaystyle {\mathit{q}}_{\mathit{i}}|g_i\sim Dir\left({\alpha }^{\left(g_i\right)}\right)}$

For full details see See Text S1 and Fig. S1.

Testing for differences in signature exposures

To characterize the signature contributions for different sets of tumor catalogs, we wish to conduct a hypothesis test that there is no difference in mean exposures versus the alternative that the mean exposure of at least one signature differs between the two groups, i.e., H₀:μ⁽¹⁾ = μ⁽²⁾ vs. H₁: at least one μ_k⁽¹⁾ ≠ μ_k⁽²⁾. We propose both local and global tests, implemented in a Bayesian framework. The former provides signature-level evaluations to determine where the differences in mean mutational exposures occur, while the latter provides an overall conclusion about any difference in mean mutational exposures. The details of our implementation are given in our Just Another Gibbs Sampler (JAGS) scripts and Source code is freely available in GitHub at https://github.com/USCbiostats/HiLDA (Plummer, 2003).

Two-stage inference methods using the point estimates of mutational exposures

An alternative approach is to perform hypothesis testing using point estimates of the mutational exposures, ${\hat{\mathit{q}}}_i$, in a two-stage analysis, which we refer to as the “two-stage” method (TS). We used the R package pmsignature to estimate $\hat{\mathit{q}}$ (Shiraishi et al., 2015). Other methods are also available, but we selected pmsignature for the purpose of comparisons to the results from HiLDA since it assumes the same model for estimating signatures under independence of features. We summarize the steps of the TS method as follows:

1. Jointly estimate the vectors of mutational signature exposures, q_i, for each mutational catalog.

2. Test for differential mutational exposures for signature k by performing the Wilcoxon rank-sum test on the $\widehat{{\mathit{q}}_k}$.

However, we note that the Wilcoxon rank-sum test in stage 2 is also sensitive to changes in variance across the two groups, which might lead to significant results even when there has been no change in mean exposures (Kasuya, 2001; Ruxton, 2006). We implemented the two-stage method using R version 3.5.0 (R Core Team, 2017). A two-sided P value of less than 0.05 was considered statistically significant.

Choosing the number of signatures

The number of signatures, K, needs to be determined prior to any of the above analyses. We adopted the method of Shiraishi et al. (2015) to determine K. Their method is based on the following criteria:

1. The optimal value of K is selected over a range of K values such that the likelihood remains relatively high while simultaneously having relatively low standard errors for the parameters.

2. Pairwise correlations between any two signatures (the kth signature and the k′th signature, say) are measured by calculating the Pearson correlation between their estimated mutational exposures across all samples, (i.e., the correlation between $\left({\hat{q}}_{1,k},\dots ,{\hat{q}}_{I,k}\right)$ and $\left({\hat{q}}_{1,k^{\prime }},\dots ,{\hat{q}}_{I,k^{\prime }}\right)$). K is chosen such that no strong correlation (i.e., >0.6) exists between any pair.

For full details see Shiraishi et al. (2015).

Application

Tumor evolution in USC colon cancer data

A total of 12,554 somatic single-nucleotide substitutions were identified, with a median of 277 per sample (range: 82–1,762) (see Table S1). One tumor with microsatelite instability has more than double the number of somatic mutations (1,751 side A, 1,762 side B) than any of the remaining 30 catalogs (all <750 mutations). In our first analysis, we compared the mutational exposures in side A to those in side B. If the tumors represent a single clonal expansion, we would expect similar mutational exposure frequencies in the two catalogs from the same tumor. Indeed, this is what we found (Table 2).

We identified a median of 174 trunk and 186 branch mutations per tumor. The numbers ranged from 49 to 1,578 trunk mutations and from 66 to 503 branch mutations (Fig. 1A). Interestingly, the microsatellite instable tumor had the most trunk mutations, but not the most branch mutations, suggesting that during tumor growth the mutation frequency is similar in microsatellite stable and instable tumors. Figs. 1B and 1C show that the C >T substitution is most common in all trunk catalogs, and most branch catalogs. The spontaneous deamination of methylated Cs in CpGs is known to contribute to hotspots of C >T mutation in the genome.

We identified three mutational signatures in our data (see Fig. S2). Those three signatures, and their corresponding exposures, are depicted in Figs. 2A, 2B and 2C. For each mutational signature, we compute the probabilities for the 1536 possible five-base signature patterns by taking the product of the feature component probabilities. We use these multinomial vectors to calculate the cosine similarity between pairs of signatures (Yang et al., 2019). The signature shown in the yellow box in Fig. 2D, involving C >T mutations at CpG sites, resembles signature 7 in Shiraishi et al. (2015) (cosine similarity 0.95), where it was identified in 25 out of 30 cancer types and likely relates to the deamination of 5-methylcytosine (‘aging’); the signature in the orange box in Fig. 2E, involving T>G mutations at GpGpTpG sites, is novel; the third signature, in the red box in Fig. 2F, is most similar to signature 23 in Shiraishi et al. (2015) (cosine similarity 0.85), where it was identified in four other cancer types. The pairwise cosine similarities between pairs of our yellow, orange and red signatures are 0.12, 0.01, and 0.02 which are rather dissimilar from each other given the [0, 1] range for cosine similarity. Using HiLDA, we test whether the three signatures differ in mean exposure between trunk and branch mutations.

Our global test strongly suggests that, in our data, the signature exposures differ between trunk and branch catalogs (Bayes Factor = 1265.0). A Bayes Factor greater than 10 is considered strong evidence for model H₁ (Jeffreys, 1998). Each of the individual signatures (depicted in Figs. 2D, 2E and 2F) is found to differ in exposure between the two sample groups, a conclusion supported by both HiLDA and the two-stage method (Table 3). From Fig. 2C, it is evident that the exposures of the first (‘aging’) signature in trunk mutations is almost always greater than that for the matching catalog of branch mutations, which is intuitively consistent with the fact that trunk mutations may well reflect an accumulation of mutations over the life of the subject, whereas branch mutations are accumulated only after tumor initiation. For the previously unseen signature, the higher exposures in branch catalogs might suggest that this signature’s underlying mechanism for generating mutations might be associated with the processes occurring during tumor evolution as opposed to normal development. From Fig. 2G, we observed that the distributional ranges of the two groups of mutational exposures have some overlaps, but that the centers of each group, i.e., the means of mutational exposures, are clearly deviated from each other. However, the distributional radii, indicating the variances of mutational exposures, do not substantially differ between the groups.

We sought to validate the discovery of the previously unseen signature by repurposing targeted sequencing data from the same tumor set (Siegmund & Shibata, 2016) and using publicly available data from the Cancer Genome Atlas. Four T>G substitutions that we assigned to the previously unseen signature were part of a much larger independent validation set of mutations subjected to targeted, high-coverage Ampliseq technology (Siegmund & Shibata, 2016); all four of these T>G substitutions failed to validate. Further, a systematic analysis of data from the Cancer Genome Atlas also did not find evidence for this signature (Williams et al., 2016). Therefore, we cannot rule out that the signature is the result of sequencing error.

Esophageal adenocarcinoma

We reanalyzed the 146 EAC previously studied by Shiraishi et al. (2015) and recovered the same four mutational signatures, C >T at CpG (S7), C >T or A at TpC (S14), T >G or C at Cp(T >G/C)pT (S21), and a signature capturing the remaining mutations, i.e., those that do not fall into the previous three signatures. We tested for differences in mutational exposures by sex, age, smoking status, and tumor site. Only tumor site showed some evidence of differences in mutational exposure by patient subgroup (Fig. 3). The TS approach showed the mutational exposure for signature S21 was lower in the cardia/GEJ compared to the esophagus (p = 0.019) (Fig. 3A). HiLDA only identified a significant deficit in the mutational exposure for S21 in the cardia/GEJ location (−7.3% with 95% credible interval: [−11.8%, −2.7%]; HiLDA-Wald p = 0.002) (Fig. 3B). However, the HiLDA global test showed no strong evidence for associations between mutational exposures and any of the four clinical variables age, sex, smoking status or tumor site (all Bayes Factors < 1), suggesting the differences with tumor site may not be real. Still, both HiLDA-CI and HiLDA-Wald tests return significant results even when using the Bonferroni method to adjust for multiple comparisons (−7.3% with 98.75% credible interval: [−12.9%, −1.5%]; adjusted p = 0.019). See Fig. S3 for more details. We now go on to assess the reliability of results using a simulation study.

Simulation study

We conducted a simulation study to assess the performance of both HiLDA and the two-stage approach in terms of the false-positive rate (FPR) and true-positive rate (TPR), in local, univariate tests of the difference in mean exposure between two groups of mutational catalogs. We assess the functionality of the methods in a setting similar to that of the USC data, simulating somatic mutations directly using the estimated signatures (f_k) from Figs. 2D, 2E and 2F for the same number of mutational catalogs (two groups of 16 catalogs each) and somatic mutations per catalog (J_i in Table S1). The mutational exposures (q_i) were indirectly used to derive the concentration parameters of the Dirichlet distributions. The scenarios are as follows:

1. The two groups of mutational catalogs are from separate Dirichlet distributions with parameters α⁽¹⁾ = (9.2, 0.2, 7.5) and α⁽²⁾ = (4.2, 0.6, 7.3). Here, the αs corresponds to the maximum-likelihood estimated parameters from the three exposure distributions in the trunk and branch mutational catalogs. This gives mean exposures of μ⁽¹⁾ = (0.54, 0.01, 0.44) and μ⁽²⁾ = (0.35, 0.05, 0.60) in trunk and branch catalogs, respectively, for the aging signature, new signature, and random signature.

2. The two groups of mutational catalogs are from the same Dirichlet distribution, Dir(4.2, 0.6, 7.3), (so here we use the concentration parameters estimated from the branch mutational catalogs).

For each tumor, mutational exposures q_i, are drawn from the Dirichlet distribution. Each set of probabilities parameterize a multinomial distribution later used to probabilistically choose the underlying mutational signature for a mutation (See Fig. S4). Then, every mutation feature in the mutational pattern of the mutation is simulated independently from a corresponding multinomial distribution of the chosen signature. To estimate the FPRs, 1,000 sets of data were simulated for scenario 2, when there is no difference in the exposure distribution between two groups of mutational catalogs. The two-stage method is slightly conservative for 1st and 3rd signatures (resulting FPRs of 4.3%, 5.2%, and 4.3%) when testing at the 5% significant level (Table 4). In comparison, HiLDA showed better control of the FPR by using the 95% credible interval of the posterior distributions (4.8%, 5.0%, and 5.1%). The Wald test also showed control of the FPR, except in the case of the rare signature when it was noticeably lower (3.7%), presumably due to the asymmetric posterior distribution.

We then moved to scenario 1, where we simulated 200 data sets with a difference in mean exposures between the two groups of catalogs. Here, the statistical powers of both HiLDA and the two-stage method are high when detecting the difference in exposures for the 1st and 3rd signatures (Table 4). In contrast, for the 2nd signature, which has the lowest mean mutational exposure, the TPRs of all methods are lower (77.5%–85.5%). By using the 95% credible interval of posterior distributions, HiLDA is able to distinguish a difference more often than the two-stage method (99.5% vs. 99.0%, 85.5% vs. 77.5%, and 91.5% vs. 88.0%). At the same time, using the credible interval resulted in higher TPRs compared to performing a Wald test (85.5% vs. 80.5% for the 2nd signature). In summary, across tests involving these three mutational signatures, HiLDA provides higher statistical power to the TS method with a tendency of better improvement for signatures with lower mutational exposures, i.e., the power difference between HiLDA and the TS method is the highest (8%) for signature 2 with the lowest mean mutational exposures. The improvements in the power to detect the mean exposure difference is presumably due to the fact that HiLDA accounts for the uncertainty in the estimated mutational exposures and provides better model fit of the posterior distributions. All data were simulated in R 3.5.0 using the hierarchical Bayesian mixture model described in the methods section. All replicates reached convergence with an Rhat value less than 1.05 for each of the scenarios shown in Tables 2–4.