Extending HDP for Mutational Signatures: Incorporating Evolutionary Trees

My previous blog post derived a Hierarchical Dirichlet Process (HDP) for discovering mutational signatures, estimating their activities across samples, and learning the number of signatures from the data. That model shares information across samples but does not encode their evolutionary relationships.

Here I extend it with a phylogenetic tree. The prior then favors similar mutational signature activities in evolutionarily related cells.

Biological motivation

Single-cell sequencing lets us model how mutational signature activities vary within a tumor. Evolutionarily related cells should have more similar activities than distantly related ones.

In cancer, a subclone represents a set of tumor cells descending from the same ancestor and hence sharing mutations. For many measurement modalities, there are established methods to infer relationships between subclones, often represented as tree structures. Examples include:

• SNV: SCITE OncoNEM
• CN: SCICoNE  
• SNV+CN: COMPASS
• SCI Phi
• SIEVE

The model needs a reliable tree topology but does not depend on a particular construction method.

Consider the following example tree, where edges represent evolutionary relatedness:

The nodes can have different numbers of cells, including none. Dropping the cell counts and enumerating the nodes gives:

Let $\mathcal{T}$ represent the topology of the tree. Without loss of generality, we assume a single tree for our input data. If multiple trees exist, such as when separate trees are constructed for different tumors, we can combine them by introducing a new root node and connecting the root nodes of the individual trees to this new root.

This tree consists of nodes $V(\mathcal{T})=\{\color{#BEBEBE} v_1, \color{#BB6438} v_2, \color{#B25870} v_3, \color{#7630A9} v_4, \color{#6672C9} v_5\}$. Each node can represent a subclone or serve as a connecting node. For each node $v_j$ we observe trinucleotide mutations $x_{j,1},\ldots, x_{j,M_j}\in \{1,\ldots,96\}$. Note that $M_j=0$ is allowed. Let $e_{\cdot j}\in \Delta^{\infty}$ denote the activity catalogue in node $j$.

Subclones with a closer evolutionary history should have more similar mutational signature activity. Let $d$ be a measure of dissimilarity on $\Delta^{\infty}$. In the above example we expect that $d(\color{#7630A9} e_{\cdot 4}, \color{#6672C9} e_{\cdot 5}) \leq d(\color{#BB6438} e_{\cdot 2}, \color{#6672C9} e_{\cdot 5})$. When we have access to tree topology $\mathcal{T}$, we can incorporate $\mathcal{T}$ into our model as prior knowledge when inferring $\hat{z}_{ji}$ and consequently $\hat{\theta}_k$ and $\hat{e}_{kj}$.

The tree-structured HDP model

Instead of making every sample depend directly on a common base distribution $G_0$, each tree node gets a Dirichlet Process based on its parent.

Sticking to the color coding from above, the graphical model of the proposed model looks as follows:

Each variable $G$ is a Dirichlet Process. For simplicity, we have omitted the scaling parameters $\alpha_j$ for each Dirichlet process from the graph. In this example, we observe trinucleotide mutations in nodes $v_2, v_4, v_5$ but not in $v_1, v_3$.

Because $G_4$ and $G_5$ vary around $G_3$, while $G_2$ and $G_3$ vary around $G_1$, the prior makes $G_5$ more similar to $G_4$ than to $G_2$. This induces the same relationship among $\color{#6672C9} e_{\cdot 5}$, $\color{#7630A9} e_{\cdot 4}$, and $\color{#BB6438} e_{\cdot 2}$ that motivated the extension.

In the above example, $G_2$ is the base probability distribution of $G_1$. Therefore, the support of $G_2$ is a subset of the support of $G_1$, i.e., $\operatorname{supp}(G_2) \subset \operatorname{supp}(G_1) = {\tilde{\theta}_1^{G_1}, \tilde{\theta}_2^{G_1}, \dots}$. A major shift in signature activity from node $v_1$ to $v_2$ could present a challenge for the model. It would be valuable to investigate whether such a shift is biologically plausible and, if so, to test the model in that scenario.

Mathematical formulation

The sampling of $G_0$, $\theta_{j i}$, and $x_{ji}$ works exactly as in the classical HDP model from my previous post:

$$\begin{aligned} G_0 &\sim \mathrm{DP}(\alpha_0, H),\\ \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 change is in the sampling of $G_j$. Consider a tree topology $\mathcal{T}$. We write $(v_i,v_s)\in \mathcal{T}$ if and only if the $i$-th node is a child of the $s$-th node. Then we model the Dirichlet Process at each node $v_j$ of the tree $\mathcal{T}$ as:

$$\begin{aligned} G_j \mid G_s &\sim \mathrm{DP}(\alpha_j, G_s) \text{ for } (v_i,v_s)\in \mathcal{T}. \end{aligned}$$

In words, at each node of the tree we draw $G_j$ from a Dirichlet Process with an individual scaling parameter $\alpha_j$ and the base probability distribution $G_s$ from the parent node. By the mathematical properties of the Dirichlet Process, $G_j$ varies around the probability distribution $G_s$ of its parent node. Hence, distribution properties of $G_s$ get propagated to $G_j$ and further down the tree, ensuring that evolutionarily related nodes have similar mutational signature activities.

Data requirements

The model requires sequencing at the single-cell or subclonal level to construct a tree. The tumor sample must also contain enough trinucleotide mutations, typically seen in late-stage tumors, and the sequencing must cover a large enough genomic region to measure them.

Implementation and evaluation

Implementation

The current prototype is a fork of the HDP implementation by Nicola Roberts in R with some helper functions. A Python port would make it accessible to more users.

Benchmarking against standard methods

The tree-structured HDP also needs to be benchmarked against standard approaches. The evaluation would involve:

1. Dataset Selection: Choose a dataset suitable for the proposed method with available tree topology
2. Method Comparison: Estimate mutational signatures using:
   • Classical NMF approach: $\Theta^{\text{NMF}}$ and $E^{\text{NMF}}$
   • Tree-structured HDP: $\Theta^{\text{Tree-HDP}}$ and $E^{\text{Tree-HDP}}$
3. Evaluation: Compare how well the estimated signature activities $D^{\text{NMF}}$ and $D^{\text{Tree-HDP}}$ respect the known tree topology $\mathcal{T}$

Alternative approaches

Instead of the Hierarchical Dirichlet Process approach, we could use variations of Non-negative Matrix Factorization. Given observed mutation catalogue $X$ and a constructed tree topology $\mathcal{T}$, we could introduce a penalization term $R_{\mathcal{T}}(E)$ to ensure that the tree topology is respected:

$$\min_{S, E \geq 0} D(X, S E) + \lambda R_{\mathcal{T}}(E).$$

A desirable property of the function $R_{\mathcal{T}}$ would be that for samples $i,j\in \{1,...,l \}$ attached to nodes in the graph $\mathcal{T}$ which are close to each other, we have similar $E_{\cdot i}$ and $E_{\cdot j}$.

Alternatively, we could explore hierarchical NMF approaches. There is abundant literature on NMF and its variations, including Sugahara et al., Ding et al., Ferreira et al., and Schmidt & Raphael.

Extension to other mutational events

Instead of focusing solely on trinucleotide mutational events, we could study other mutational events with the proposed model. For example, we could study chromosomal instability events using the approach described in Drews et al.. However, this would introduce additional requirements to the input data, such as high read depth, and such datasets are not yet widely available.

Recommended material for chromosomal instability events

• https://markowetz.cruk.cam.ac.uk/cincompendium/
• https://github.com/markowetzlab/CINSignatureQuantification
  • calculateFeatures function
  • R package to quantify signatures of chromosomal instability on absolute copy number profiles as described in Drews et al.

Comparison with existing literature

The proposed approach should also be compared with existing work such as Alam et al. to identify its theoretical and practical differences.

Conclusion

The tree-structured HDP encodes the expectation that evolutionarily related cells should have more similar mutational signature activities.

Each node has its own Dirichlet Process based on its parent, so distribution properties propagate down the tree instead of treating samples independently.

Making the model practical still requires suitable data, an implementation, and benchmarks against the alternatives described above.