Skip to main content
Science
View as Markdown Suggest changes

Mutational signature estimation with Hierarchical Dirichlet Process

· Reading time: 11 min
Mutational signature estimation with Hierarchical Dirichlet Process

The Hierarchical Dirichlet Process (HDP) is used in several fields, including cancer genomics. While researching methods for mutational signature estimation, I found several studies that applied HDPs but gave little detail on how the model maps to this problem.

Nicola Roberts’ dissertation describes parts of the mathematical model in its appendix. This post fills in the derivation, starting with a finite mixture model for mutational signatures, then moving to the Dirichlet Process and the Hierarchical Dirichlet Process.

Introduction

Cancer cells have elevated rates of somatic mutation. The resulting mutation patterns help identify mechanisms that drive the disease (Fischer et al., 2013), which can inform treatment and prevention.

One type of somatic mutation is the single nucleotide variant (SNV), also referred to as a single base substitution. This involves the replacement of one base pair in the genome with another. There are six possible base pair mutations: C:G>A:T, C:G>G:C, C:G>T:A, T:A>A:T, T:A>C:G, and T:A>G:C. When considering the context of the immediate 5’ and 3’ neighboring bases, the total number of mutations increases to 96, calculated as 6×4×4=966 \times 4 \times 4 = 96. These are known as trinucleotide mutation channels. Examples of such mutations include ACA:TGT>AAA:TTT and ACT:TGA>AGT:TCA. For simplicity, we will abbreviate these mutations as ACA>AAA and ACT>AGT respectively.

Mutational signatures

A mutational signature θ\theta with respect to trinucleotide mutation channels is a discrete probability distribution over the 96 mutational channels i.e. θ∈Δ96\theta\in \Delta^{96} where

Δn={(x1,x2,…,xn)∈Rn∣xi≥0 for all i and ∑i=1nxi=1}.\Delta^{n} = \left\{ (x_1, x_2, \ldots, x_{n}) \in \mathbb{R}^{n} \mid x_i \geq 0 \text{ for all } i \text{ and } \sum_{i=1}^{n} x_i = 1 \right\}.

A mutational signature θ\theta can be visualised as in the following.

500

The framework for studying somatic mutations through mutational signatures was introduced in a landmark study by Alexandrov et al. (2013), where over 7,000 bulk-sequenced cancer samples were analyzed. Conceptually, a mutational signature represents a biological process acting on the genome, leaving a distinct imprint captured in a probability vector, θ\theta. The COSMIC project maintains a list of mutational signatures and their proposed underlying biological processes. Some of the signatures and their probability vectors θ\theta are consistently identifiable across most examined cohorts and signature identification algorithms (cite).

To estimate mutational signatures and their activities, Alexandrov et al. employed Non-Negative Matrix Factorization (NMF). NMF and its variations remain the most widely used methods for estimating mutational signatures and mutational signature activities.

Dirichlet Process

The Dirichlet Process (DP), denoted DP(α,H)\mathrm{DP}(\alpha, H), is a distribution over probability measures. Sampling G∼DP(α,H)G \sim \mathrm{DP}(\alpha, H) gives a probability measure GG. Its two parameters are the concentration α>0\alpha>0 and the base measure HH; the support of GG is a subset of the support of HH.

To sample from the Dirichlet process, we begin by drawing an infinite sequence of i.i.d. samples θ~1,θ~2,…∼H\tilde{\theta}_1, \tilde{\theta}_2, \ldots \sim H. Then, we draw a vector e={e1,e2,…}∼Stick⁡(α)\boldsymbol{e} = \{e_1, e_2, \ldots\} \sim \operatorname{Stick}(\alpha), where Stick⁡(α)\operatorname{Stick}(\alpha) refers to the stick-breaking process. Intuitively, the stick-breaking process works by repeatedly breaking off and discarding a random fraction of a “stick” that is initially of length 11. Each broken off part of the stick has a length eke_k. The stick-breaking process generates a sequence e∈Δ∞\boldsymbol{e} \in \Delta^\infty, where the total length of all the broken pieces is 1.

Finally, we construct the random measure GG as:

G=∑k=1∞ekδθ~kG = \sum_{k=1}^{\infty} e_k \delta_{\tilde{\theta}_k}

where δθ~k\delta_{\tilde{\theta}_k} is a Dirac delta function centered at θ~k\tilde{\theta}_k. The resulting measure GG is discrete, with countable support given by {θ~1,θ~2,…}\{\tilde{\theta}_1, \tilde{\theta}_2, \ldots\}.

Mathematical properties

The base measure HH is the mean of the random measure GG. Specifically, for any measurable set AA, the expectation of G(A)G(A) under the Dirichlet process is given by [@millerDirichletProcessModels]:

EG∼DP(α,H)[G(A)]=H(A)\mathrm{E}_{G \sim \mathrm{DP}(\alpha, H)}[G(A)]=H(A)

As α→∞\alpha \rightarrow \infty, the distribution of GG converges to HH in the weak topology, meaning that larger values of α\alpha make GG more closely resemble the base measure HH.

GG is centered around the prior distribution HH, with the concentration parameter α\alpha controlling how tightly GG clusters around HH. A higher α\alpha reduces variability, causing GG to approximate HH more closely.

Dirichlet Process resources

Recommended material for Dirichlet Process

Notebook implementations of the Dirichlet Process

Modelling mutational signatures with a Dirichlet Process

Previous studies have used Dirichlet processes to model trinucleotide mutations associated with mutational signatures (cite) and (cite). Appendix B (p. 234) and Chapter 4 (p. 132) in (cite) provide an overview but omit parts of the derivation. This section works through the model, beginning with a finite mixture. Consider the trinucleotide mutations within a single sample, encoded as x∈{1,…,96}x \in \{1, \ldots, 96\} for the 96 possible mutation types.

Step 1: Known number of mutational signatures KK

First assume we know the number of active mutational signatures, denoted by KK. Each signature is a discrete probability distribution θ~k∈Δ96\tilde{\theta}_k \in \Delta^{96}, where k=1,…,Kk=1, \ldots, K. With no prior knowledge of the specific signatures, we model these distributions as:

θ~k∼Dir⁡(196⋅196) for k=1,…,K,\begin{aligned} \tilde{\theta}_k & \sim \operatorname{Dir}\left(\frac{1}{96} \cdot \mathbf{1}_{96}\right) \text{ for } k=1,\ldots, K\\ \end{aligned},

where Dir⁡\operatorname{Dir} denotes the Dirichlet distribution, serving as a symmetric prior. Similarly, we assume no prior information about the activity of these mutational signatures. We model the activity of the KK mutational signatures e∈ΔK\boldsymbol{e} \in \Delta^K in a sample as:

e∣α,K∼Dir⁡(αK⋅1K).\begin{aligned} \boldsymbol{e} \mid \alpha,K & \sim \operatorname{Dir}\left(\frac{\alpha}{K} \cdot \mathbf{1}_K\right). \end{aligned}

where α\alpha is a hyperparameter controlling the concentration of the distribution. Next, let zi=kz_i=k represent the event that the ii-th mutation was generated by the kk-th mutational signature. If we observe a total of MM trinucleotide mutations in our sample, we draw each zi∈{{1,…,K}}z_i\in \{\{1,\ldots, K\}\} according to the mutational signature activity e\boldsymbol{e}, as follows:

zi∣e∼Categorical⁡(e) for i=1,…,M,\begin{aligned} z_i \mid \boldsymbol{e} & \sim \operatorname{Categorical} \left(\boldsymbol{e}\right) \text{ for } i=1,\ldots, M, \end{aligned}

indicating that the probability P(zi=k)=ek\mathrm{P}\left(z_i=k\right)=e_k, where eke_k is the activity level of signature kk. Finally, the observed trinucleotide mutation xi∈{1,…,96}x_i \in\{1, \ldots, 96\} is drawn from a categorical distribution based on the corresponding mutational signature zi=kz_i=k and its distribution θ~k\tilde{\theta}_k :

{1,…,96}∋xi∣zi=k,θ~k∼Categorical⁡(θ~k) for i=1,…,M\begin{aligned} \{{1,\ldots, 96}\} \ni x_i \mid z_i=k, \tilde{\theta}_k & \sim \operatorname{Categorical} \left(\tilde{\theta}_{k}\right) \text{ for } i=1,\ldots, M \end{aligned}

This model describes how the mutational signatures {θ~1,…,θ~K}\left\{\tilde{\theta}_1, \ldots, \tilde{\theta}_K\right\} generate the observed trinucleotide mutations {x1,…,xM}\left\{x_1, \ldots, x_M\right\}.

Graphical model for known number of components KK

We can visualize the process that generates trinucleotide mutations in a sample using a graphical model:

300

Here the grey background color indicates that xix_i is an observed random variable.

Step 2: Equivalent model as mixture

Next, we integrate out the indicator variable zi=kz_i=k to obtain an equivalent mixture representation. As before, we draw the mutational signature activities as follows:

e∣α,K∼Dir⁡(αK⋅1K).\begin{aligned} \boldsymbol{e} \mid \alpha,K & \sim \operatorname{Dir}\left(\frac{\alpha}{K} \cdot \mathbf{1}_K\right). \end{aligned}

Let H=Dir⁡(196⋅196)H=\operatorname{Dir}\left(\frac{1}{96} \cdot \mathbf{1}_{96}\right). Similar to the last step, we generate the mutational signatures as

θ~k∼H for k=1,…,K.\begin{aligned} \tilde{\theta}_k & \sim H \text{ for } k=1,\ldots, K. \end{aligned}

Now, instead of assuming that each observation i=1,…,Mi=1, \ldots, M is first assigned a mutational signature zi∈{1,…,K}z_i \in\{1, \ldots, K\} and then drawn from the corresponding signature θ~zi\tilde{\theta}_{z_i} we directly draw θi\theta_i from {θ~1,…,θ~K}\left\{\tilde{\theta}_1, \ldots, \tilde{\theta}_K\right\}, with probabilities {e1,…,eK}\left\{e_1, \ldots, e_K\right\} for each observation:

θi∣e,θ~1,…θ~K∼G=∑k=1Kekδθ~k(θi) for i=1,…,M.\begin{aligned} \theta_i \mid \boldsymbol{e}, \tilde{\theta}_1, \ldots \tilde{\theta}_K & \sim G=\sum_{k=1}^K e_k \delta_{\tilde{\theta}_k}\left(\theta_i\right) \text{ for } i=1,\ldots, M. \end{aligned}

Here, GG represents a finite mixture distribution with support on the acting mutational signatures. Finally, we draw the trinucleotide mutation as follows:

xi∣θi∼Categorical⁡(θi) for i=1,…,M.\begin{aligned} x_i \mid \theta_i & \sim \operatorname{Categorical} \left(\theta_i\right) \text{ for } i=1,\ldots, M. \end{aligned}

Step 3: Unknown, unbounded number of components KK

Let

Δ∞={(θ1,θ2,…,)∈R∞∣θi≥0 for all i and ∑i=1∞θi=1}.\Delta^{\infty} = \left\{ (\theta_1, \theta_2, \ldots, ) \in \mathbb{R}^{\infty} \mid \theta_i \geq 0 \text{ for all } i \text{ and } \sum_{i=1}^{\infty} \theta_i = 1 \right\}.

We now extend the model from the previous step to allow for an unknown number of components KK, with any K∈NK \in \mathbb{N} being possible. As in the previous step, we generate the underlying mutational signatures θ~k∈Δ96\tilde{\theta}_k\in \Delta^{96} from HH as follows:

θ~k∼H for k=1,2,…\begin{aligned} \tilde{\theta}_k & \sim H \text{ for } k=1,2,\ldots \end{aligned}

Previously, we sampled a finite set of KK mutational signatures θ~k\tilde{\theta}_k, but now we generate an infinite sequence of signatures. Accordingly, the mutational signature activities e\boldsymbol{e} now lie in Δ∞\Delta^{\infty}. To model these activities, we use the Stick-breaking process, which defines a distribution over Δ∞\Delta^{\infty}. Thus, we sample the mutational signature activities as:

e∣α∼Stick⁡(α).\begin{aligned} \boldsymbol{e} \mid \alpha &\sim \operatorname{Stick}(\alpha). \end{aligned}

Intuitively, the stick-breaking process works by repeatedly breaking off and discarding a random fraction of a “stick” that is initially of length 11. Each broken off part of the stick has a length eke_k. The stick-breaking process generates a sequence e∈Δ∞\boldsymbol{e} \in \Delta^\infty, where the total length of all the broken pieces is 11. Consequently, the mixture distribution GG, which previously had finite support, now has countably infinite support {θ~1,θ~2,…}\left\{\tilde{\theta}_1, \tilde{\theta}_2, \ldots\right\} and takes the form:

P(Δ96)∋G=∑k=1∞ekδθ~k(θi),\begin{aligned} \mathcal{P}(\Delta^{96}) \ni G&=\sum_{k=1}^{\infty} e_k \delta_{\tilde{\theta}_k}\left(\theta_i\right) , \end{aligned}

where P(Δ96)\mathcal{P}(\Delta^{96}) denotes the space of probability distributions on Δ96\Delta^{96}. We draw θi\theta_i from the support of GG, which is {θ~1,θ~2,…}\left\{\tilde{\theta}_1, \tilde{\theta}_2, \ldots\right\}, as follows:

θi∣G∼G for i=1,…,M\begin{aligned} \theta_i \mid G & \sim G \text{ for } i=1,\ldots, M \\ \end{aligned}

Finally, we observe the mutations as

{1,…,96}∋xi∣θi∼Categorical⁡(θi) for i=1,…,M\begin{aligned} \{{1,\ldots, 96}\} \ni x_i \mid \theta_i & \sim \operatorname{Categorical} \left(\theta_i\right) \text{ for } i=1,\ldots, M \end{aligned}

Step 4: Shorthand notation

The Dirichlet Process notation summarizes the construction of GG from α\alpha and HH:

G∼DP(α,H)Δ96∋θi∣G∼G for i=1,…,M{1,…,96}∋xi∣θi∼Categorical⁡(θi) for i=1,…,M\begin{aligned} G &\sim \mathrm{DP}(\alpha, H)\\ \Delta^{96}\ni \theta_i \mid G &\sim G \text{ for } i=1,\ldots, M \\ \{{1,\ldots, 96}\} \ni x_i \mid \theta_i & \sim \operatorname{Categorical} \left(\theta_i\right) \text{ for } i=1,\ldots, M \end{aligned}

This completes the derivation of a Dirichlet Process model for trinucleotide mutations, as used in (cite). It provides an alternative to the classical Non-negative Matrix Factorisation approach used in (cite) for estimating mutational signatures and their activities.

Graphical model of a Dirichlet Process

We can visualise the Dirichlet process as a graphical model:

ThegraphicalmodelofaDirichletProcess.

Hierarchical Dirichlet Process

A Hierarchical Dirichlet Process (HDP) composes multiple Dirichlet Processes (cite). For trinucleotide mutations, it gives each sample a separate mixture while sharing signatures across samples. We start with the prior distribution HH and draw a common base distribution G0G_0:

G0∼DP(α0,H).\begin{aligned} G_0 &\sim \mathrm{DP}(\alpha_0, H).\\ \end{aligned}

Assume we have NN samples. For each sample, we draw a separate Dirichlet Process GjG_j from the common base distribution G0G_0 as

Gj∣G0∼DP(αj,G0) for j=1,2,…N.\begin{aligned} G_j \mid G_0 &\sim \mathrm{DP}(\alpha_j, G_0) \text{ for } j = 1,2,\ldots N. \end{aligned}

The Dirichlet Process G0G_0 sits at the top, with the processes GjG_j below it. Each GjG_j varies around G0G_0, allowing samples to share information through the common base distribution. If we observe MjM_j mutations in sample jj, we draw them as before:

θji∣Gj∼Gj for i=1,2,…Mj,xji∣θji∼Categorical⁡(θji).\begin{aligned} \theta_{j i} \mid G_j &\sim G_j \text{ for } i = 1,2,\ldots M_j , \\ x_{ji} \mid \theta_{j i} & \sim \operatorname{Categorical} \left(\theta_{j i}\right). \end{aligned}

The corresponding graphical model is:

500

Estimating signatures and their activities

Assume the mutational signatures (θk)k∈N⊂Δ96\left(\theta_k\right)_{k\in \mathbb{N}}\subset \Delta^{96} act on the genome. Borrowing notation from the finite model, let zji∈Nz_{ji}\in \mathbb{N} denote the signature that generated mutation xji∈{1,…,96}x_{ji}\in \{1,\ldots,96\}. Fitting an HDP gives an estimate z^ji\hat{z}_{ji} of zjiz_{ji}. The estimate θ^k\hat{\theta}_k for the probability vector θk∈Δ96\theta_k \in \Delta^{96} of signature kk is

θ^k=∑l=196∑j=1N∑i=1MjI(z^ji=k)I(xji=l)δl∑j=1N∑i=1MjI(z^ji=k),\hat{\theta}_k=\frac{\sum_{l=1}^{96}\sum_{j=1}^N\sum_{i=1}^{M_j} \mathbb{I}(\hat{z}_{ji}=k)\mathbb{I}(x_{ji}=l)\delta_l}{\sum_{j=1}^N\sum_{i=1}^{M_j} \mathbb{I}(\hat{z}_{ji}=k)},

where δl∈Δ96\delta_l\in \Delta^{96} has mass 11 at position ll. Further, we estimate the mutational signature activity ekje_{kj} of mutational signature kk in sample jj by

e^kj=1/Mj∑i=1MjI(z^ji=k).\hat{e}_{kj}=1/M_j\sum_{i=1}^{M_j} \mathbb{I}(\hat{z}_{ji}=k).

Comparison with NMF

Compared with classical Non-negative Matrix Factorization (NMF), the HDP can incorporate prior knowledge and group structure (cite). It estimates signature activity while discovering signatures de novo, and can match observations to an existing signature library while identifying new signatures by pseudo-counting that library as observational data. It also learns the number of signatures from the data. Some NMF variants can quantify uncertainty, but the HDP provides this more directly.

References

Python implementations 1

Python implementations 2

Recommended material for HDP

AI Chat

Messages you send are processed by the Google Gemini API to generate responses. Do not share sensitive personal data. See the privacy policy for details.