---
title: "Mutational signature estimation with Hierarchical Dirichlet Process"
pubDate: 2024-09-15T00:00:00.000Z
updatedDate: 2024-09-15T00:00:00.000Z
tags:
- Science
math: true
image: /images/blog/mutational-signature-with-hierarchical-dirichlet-process/featured.png
---
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](https://doi.org/10.17863/CAM.22674) 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 \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. $\theta\in \Delta^{96}$ where
$$
\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.
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](https://cancer.sanger.ac.uk/signatures/sbs/) 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)](https://pubmed.ncbi.nlm.nih.gov/34316057/).
To estimate mutational signatures and their activities, [Alexandrov et al.](https://pubmed.ncbi.nlm.nih.gov/23318258/) 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 $\mathrm{DP}(\alpha, H)$, is a distribution over probability measures. Sampling $G \sim \mathrm{DP}(\alpha, H)$ gives a probability measure $G$. Its two parameters are the concentration $\alpha>0$ and the base measure $H$; the support of $G$ is a subset of the support of $H$.
To sample from the Dirichlet process, we begin by drawing an infinite sequence of i.i.d. samples $\tilde{\theta}_1, \tilde{\theta}_2, \ldots \sim H$. Then, we draw a vector $\boldsymbol{e} = \{e_1, e_2, \ldots\} \sim \operatorname{Stick}(\alpha)$, where $\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 $1$. Each broken off part of the stick has a length $e_k$. The stick-breaking process generates a sequence $\boldsymbol{e} \in \Delta^\infty$, where the total length of all the broken pieces is 1.
Finally, we construct the random measure $G$ as:
$$G = \sum_{k=1}^{\infty} e_k \delta_{\tilde{\theta}_k}$$
where $\delta_{\tilde{\theta}_k}$ is a Dirac delta function centered at $\tilde{\theta}_k$. The resulting measure $G$ is discrete, with countable support given by $\{\tilde{\theta}_1, \tilde{\theta}_2, \ldots\}$.
### Mathematical properties
The base measure $H$ is the mean of the random measure $G$. Specifically, for any measurable set $A$, the expectation of $G(A)$ under the Dirichlet process is given by \[@millerDirichletProcessModels\]:
$$
\mathrm{E}_{G \sim \mathrm{DP}(\alpha, H)}[G(A)]=H(A)
$$
As $\alpha \rightarrow \infty$, the distribution of $G$ converges to $H$ in the weak topology, meaning that larger values of $\alpha$ make $G$ more closely resemble the base measure $H$.
$G$ is centered around the prior distribution $H$, with the concentration parameter $\alpha$ controlling how tightly $G$ clusters around $H$. A higher $\alpha$ reduces variability, causing $G$ to approximate $H$ more closely.
### Dirichlet Process resources
Recommended material for Dirichlet Process
- https://en.wikipedia.org/wiki/Dirichlet_process
- https://www.geeksforgeeks.org/dirichlet-process-mixture-models-dpmms/
- https://towardsdatascience.com/dont-be-afraid-of-nonparametric-topic-models-d259c237a840
- https://pawel-czyz.github.io/posts/dirichlet-process.html
Notebook implementations of the Dirichlet Process
- https://datamicroscopes.github.io/ncluster.html
- https://github.com/tdhopper/notes-on-dirichlet-processes/tree/master
- https://github.com/tdhopper/notes-on-dirichlet-processes/blob/master/pages/2015-07-28-dirichlet-distribution-dirichlet-process.ipynb
## Modelling mutational signatures with a Dirichlet Process
Previous studies have used Dirichlet processes to model trinucleotide mutations associated with mutational signatures [(cite)](https://doi.org/10.17863/CAM.22674) and [(cite)](https://www.nature.com/articles/s41586-019-1913-9). Appendix B (p. 234) and Chapter 4 (p. 132) in [(cite)](https://doi.org/10.17863/CAM.22674) 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 \in \{1, \ldots, 96\}$ for the 96 possible mutation types.
#### Step 1: Known number of mutational signatures $K$
First assume we know the number of active mutational signatures, denoted by $K$. Each signature is a discrete probability distribution $\tilde{\theta}_k \in \Delta^{96}$, where $k=1, \ldots, K$. With no prior knowledge of the specific signatures, we model these distributions as:
$$
\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 $\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 $K$ mutational signatures $\boldsymbol{e} \in \Delta^K$ in a sample as:
$$
\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 $z_i=k$ represent the event that the $i$-th mutation was generated by the $k$-th mutational signature. If we observe a total of $M$ trinucleotide mutations in our sample, we draw each $z_i\in \{\{1,\ldots, K\}\}$ according to the mutational signature activity $\boldsymbol{e}$, as follows:
$$
\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 $\mathrm{P}\left(z_i=k\right)=e_k$, where $e_k$ is the activity level of signature $k$. Finally, the observed trinucleotide mutation $x_i \in\{1, \ldots, 96\}$ is drawn from a categorical distribution based on the corresponding mutational signature $z_i=k$ and its distribution $\tilde{\theta}_k$ :
$$
\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 $\left\{\tilde{\theta}_1, \ldots, \tilde{\theta}_K\right\}$ generate the observed trinucleotide mutations $\left\{x_1, \ldots, x_M\right\}$.
##### Graphical model for known number of components $K$
We can visualize the process that generates trinucleotide mutations in a sample using a graphical model:
Here the grey background color indicates that $x_i$ is an observed random variable.
#### Step 2: Equivalent model as mixture
Next, we integrate out the indicator variable $z_i=k$ to obtain an equivalent mixture representation. As before, we draw the mutational signature activities as follows:
$$
\begin{aligned}
\boldsymbol{e} \mid \alpha,K & \sim \operatorname{Dir}\left(\frac{\alpha}{K} \cdot \mathbf{1}_K\right).
\end{aligned}
$$
Let $H=\operatorname{Dir}\left(\frac{1}{96} \cdot \mathbf{1}_{96}\right)$. Similar to the last step, we generate the mutational signatures as
$$
\begin{aligned}
\tilde{\theta}_k & \sim H \text{ for } k=1,\ldots, K.
\end{aligned}
$$
Now, instead of assuming that each observation $i=1, \ldots, M$ is first assigned a mutational signature $z_i \in\{1, \ldots, K\}$ and then drawn from the corresponding signature $\tilde{\theta}_{z_i}$ we directly draw $\theta_i$ from $\left\{\tilde{\theta}_1, \ldots, \tilde{\theta}_K\right\}$, with probabilities $\left\{e_1, \ldots, e_K\right\}$ for each observation:
$$
\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, $G$ represents a finite mixture distribution with support on the acting mutational signatures. Finally, we draw the trinucleotide mutation as follows:
$$
\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 $K$
Let
$$
\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 $K$, with any $K \in \mathbb{N}$ being possible. As in the previous step, we generate the underlying mutational signatures $\tilde{\theta}_k\in \Delta^{96}$ from $H$ as follows:
$$
\begin{aligned}
\tilde{\theta}_k & \sim H \text{ for } k=1,2,\ldots
\end{aligned}
$$
Previously, we sampled a finite set of $K$ mutational signatures $\tilde{\theta}_k$, but now we generate an infinite sequence of signatures. Accordingly, the mutational signature activities $\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:
$$
\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 $1$. Each broken off part of the stick has a length $e_k$. The stick-breaking process generates a sequence $\boldsymbol{e} \in \Delta^\infty$, where the total length of all the broken pieces is $1$. Consequently, the mixture distribution $G$, which previously had finite support, now has countably infinite support $\left\{\tilde{\theta}_1, \tilde{\theta}_2, \ldots\right\}$ and takes the form:
$$
\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 $\mathcal{P}(\Delta^{96})$ denotes the space of probability distributions on $\Delta^{96}$. We draw $\theta_i$ from the support of $G$, which is $\left\{\tilde{\theta}_1, \tilde{\theta}_2, \ldots\right\}$, as follows:
$$
\begin{aligned}
\theta_i \mid G & \sim G \text{ for } i=1,\ldots, M \\
\end{aligned}
$$
Finally, we observe the mutations as
$$
\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 $G$ from $\alpha$ and $H$:
$$
\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)](https://www.nature.com/articles/s41586-019-1913-9). It provides an alternative to the classical Non-negative Matrix Factorisation approach used in [(cite)](https://www.nature.com/articles/nature12477) 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)](https://www.tandfonline.com/doi/abs/10.1198/016214506000000302). For trinucleotide mutations, it gives each sample a separate mixture while sharing signatures across samples. We start with the prior distribution $H$ and draw a common base distribution $G_0$:
$$
\begin{aligned}
G_0 &\sim \mathrm{DP}(\alpha_0, H).\\
\end{aligned}
$$
Assume we have $N$ samples. For each sample, we draw a separate Dirichlet Process $G_j$ from the common base distribution $G_0$ as
$$
\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 $G_0$ sits at the top, with the processes $G_j$ below it. Each $G_j$ varies around $G_0$, allowing samples to share information through the common base distribution. If we observe $M_j$ mutations in sample $j$, we draw them as before:
$$
\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:
### Estimating signatures and their activities
Assume the mutational signatures $\left(\theta_k\right)_{k\in \mathbb{N}}\subset \Delta^{96}$ act on the genome. Borrowing notation from the finite model, let $z_{ji}\in \mathbb{N}$ denote the signature that generated mutation $x_{ji}\in \{1,\ldots,96\}$. Fitting an HDP gives an estimate $\hat{z}_{ji}$ of $z_{ji}$. The estimate $\hat{\theta}_k$ for the probability vector $\theta_k \in \Delta^{96}$ of signature $k$ is
$$\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 $\delta_l\in \Delta^{96}$ has mass $1$ at position $l$. Further, we estimate the mutational signature activity $e_{kj}$ of mutational signature $k$ in sample $j$ by
$$\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)](https://doi.org/10.17863/CAM.22674). 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
- https://radimrehurek.com/gensim/models/hdpmodel.html
- Gensim is a Python library for *topic modelling*, *document indexing* and *similarity retrieval* with large corpora. Target audience is the *natural language processing* (NLP) and *information retrieval* (IR) community.
- https://github.com/piskvorky/gensim/blob/develop/gensim/models/hdpmodel.py
- based on https://github.com/blei-lab/online-hdp
- https://towardsdatascience.com/dont-be-afraid-of-nonparametric-topic-models-d259c237a840
- https://towardsdatascience.com/dont-be-afraid-of-nonparametric-topic-models-part-2-python-e5666db347a
- https://github.com/ecoronado92/hdp
- https://github.com/morrisgreenberg/hdp-py
Python implementations 2
- https://radimrehurek.com/gensim/models/hdpmodel.html
- Gensim is a Python library for *topic modelling*, *document indexing* and *similarity retrieval* with large corpora. Target audience is the *natural language processing* (NLP) and *information retrieval* (IR) community.
- https://github.com/piskvorky/gensim/blob/develop/gensim/models/hdpmodel.py
- based on https://github.com/blei-lab/online-hdp
- https://towardsdatascience.com/dont-be-afraid-of-nonparametric-topic-models-d259c237a840
- https://towardsdatascience.com/dont-be-afraid-of-nonparametric-topic-models-part-2-python-e5666db347a
- https://github.com/ecoronado92/hdp
- https://github.com/morrisgreenberg/hdp-py
- Simple Jupyter notebook: https://github.com/tdhopper/notes-on-dirichlet-processes/blob/master/pages/2015-07-30-sampling-from-a-hierarchical-dirichlet-process.ipynb (just for sampling?)
Recommended material for HDP
- https://en.wikipedia.org/wiki/Hierarchical_Dirichlet_process
- https://mlg.eng.cam.ac.uk/zoubin/tut06/ywt.pdf
