Quantification of tumour evolution and heterogeneity via Bayesian epiallele detection

TL;DR

This paper introduces a Bayesian model to analyze DNA methylation epialleles across tumor regions, revealing tumor evolution, heterogeneity, and epigenetic disorder with improved statistical power and uncertainty quantification.

Contribution

A novel Bayesian approach for inferring and analyzing epialleles across multiple tumor regions, enhancing understanding of tumor heterogeneity and evolution.

Findings

Identified tumor-specific epialleles and their distribution.

Reconstructed tumor phylogenetic history based on epiallele patterns.

Quantified epigenetic disorder within tumor samples.

Abstract

Motivation: Epigenetic heterogeneity within a tumour can play an important role in tumour evolution and the emergence of resistance to treatment. It is increasingly recognised that the study of DNA methylation (DNAm) patterns along the genome -- so-called `epialleles' -- offers greater insight into epigenetic dynamics than conventional analyses which examine DNAm marks individually. Results: We have developed a Bayesian model to infer which epialleles are present in multiple regions of the same tumour. We apply our method to reduced representation bisulfite sequencing (RRBS) data from multiple regions of one lung cancer tumour and a matched normal sample. The model borrows information from all tumour regions to leverage greater statistical power. The total number of epialleles, the epiallele DNAm patterns, and a noise hyperparameter are all automatically inferred from the data. Uncertainty as to which epiallele an observed sequencing read originated from is explicitly incorporated by marginalising over the appropriate posterior densities. The degree to which tumour samples are contaminated with normal tissue can be estimated and corrected for. By tracing the distribution of epialleles throughout the tumour we can infer the phylogenetic history of the tumour, identify epialleles that differ between normal and cancer tissue, and define a measure of global epigenetic disorder.

Tables

Table 1. Table 1: In the middle column are estimates of tumour purity based on a comparison of epiallele distributions between normal tissue and tumour tissue. The third column contains estimates obtained from a separate study of exome data from the same tumour samples.

Tumour sample	Epiallele purity estimate	Exome purity estimate
R1	35%	32%
R2	54%	51%
R3	75%	73%
R4	53%	67%
R5	25%	28%
R6	<20%	13%
R7	30%	36%

Equations

$$p(y_{i\mu}|x_{q\mu},\epsilon,Q)=\left\{\begin{array}[]{ll}\epsilon&\quad\text{if $y_{i\mu}\neq x_{q\mu}$}\\ 1-\epsilon&\quad\text{if $y_{i\mu}=x_{q\mu}$.}\end{array}\right.$$

$$p(\mathbf{X},\mathbf{w},ϵ|\mathbf{Y},Q)=\frac{p(\mathbf{Y}|\mathbf{X},\mathbf{w},ϵ,Q)p(\mathbf{X}|Q)p(\mathbf{w}|Q)}{p(\mathbf{Y}|Q)}$$

$$p(\mathbf{Y}|\mathbf{X},\mathbf{w},ϵ,Q)$$

$$w_i^{\ast }={\text{argmax}}_{q}p({\mathbf{y}}_i|{\mathbf{x}}_q,ϵ,Q).$$

$$x_{q\mu }^{\ast }$$

$${ϵ}^{\ast }=\frac{{\alpha }_0}{{\alpha }_0+{\alpha }_1}$$

$$\text{AIC}(Q)=-2\log p(\mathbf{Y}|\hat{\mathbf{X}},\hat{\mathbf{w}},\widehat{ϵ},Q)+2Qd$$

$$p(w_i|{\hat{\mathbf{w}}}_{-i},\hat{\mathbf{X}},\widehat{ϵ},\hat{Q})=\frac{p(\mathbf{Y}|\hat{\mathbf{X}},{\hat{\mathbf{w}}}_{-i},w_i,\widehat{ϵ},\hat{Q})p(\hat{\mathbf{X}}|\hat{Q})p(\hat{\mathbf{w}}|\hat{Q})}{p(\mathbf{Y}|\hat{Q})}$$

$${\varphi }_q=\frac{1}{N}\sum _{i=1}^{N}p(w_i=q|{\hat{\mathbf{w}}}_{-i},\hat{\mathbf{X}},\widehat{ϵ},\hat{Q}).$$

$${\varphi }_q^s=\frac{1}{N_s}\sum _{i\in I_s}p(w_i=q|{\hat{\mathbf{w}}}_{-i},\hat{\mathbf{X}},\widehat{ϵ},\hat{Q})$$

$$\mathbf{\varphi }=\rho \mathbf{t}+(1-\rho )\mathbf{n}$$

$$\xi =\frac{1}{2}\sum _{\mu =1}^{\hat{Q}}\text{abs}({\varphi }_{\mu }-n_{\mu })$$

$${\hat{t}}_{\mu }=\frac{{\varphi }_{\mu }-(1-\rho )n_{\mu }}{\rho }\text{for μ=1,…,Q^.}$$

$$\mathcal{L}(ϵ)$$

$$=\log \prod _{i=1}^N\prod _{\mu =1}^{d}p(y_{i\mu }|x_{w_i\mu },ϵ,Q)$$

$$={\alpha }_0\log ϵ+{\alpha }_1\log (1-ϵ).$$

Full text

Results: We have developed a Bayesian model to infer which epialleles are present in multiple regions of the same tumour. We apply our method to reduced representation bisulfite sequencing (RRBS) data from multiple regions of one lung cancer tumour and a matched normal sample. The model borrows information from all tumour regions to leverage greater statistical power. The total number of epialleles, the epiallele DNAm patterns, and a noise hyperparameter are all automatically inferred from the data. Uncertainty as to which epiallele an observed sequencing read originated from is explicitly incorporated by marginalising over the appropriate posterior densities. The degree to which tumour samples are contaminated with normal tissue can be estimated and corrected for. By tracing the distribution of epialleles throughout the tumour we can infer the phylogenetic history of the tumour, identify epialleles that differ between normal and cancer tissue, and define a measure of global epigenetic disorder.

Availability: R code is available at github.com/james-e-barrett.

Contact: regmjeb@ucl.ac.uk

1 Introduction

Epigenetic variability allows greater phenotypic diversity and plasticity within a population of genetically similar cells. Epigenetic diversity within a tumour provides a mechanism for clonal evolution and the emergence of resistance to therapy (Mazor et al., 2016). Persistence of treatment-resistant subclonal populations may explain the failure of some therapies, and higher levels of heterogeneity have been associated with poorer clinical outcomes (Landau et al., 2014). Analysing multiple tissue samples from different tumour regions facilitates quantification of tumour heterogeneity and phylogenetic analyses. Pan et al. (2015) showed that intra-tumour DNAm heterogeneity is predictive of time-to-relapse in diffuse B-cell lymphomas. Brocks et al. (2014) have shown that both epigenetic and genetic alterations reflect the evolutioary history of prostate cancers. A recent study of Ewing sarcoma found substantial levels of epigenetic heterogeneity within tumours (Sheffield et al., 2017).

Epigenetic modifications play an important role in the regulation of gene expression. One of the most common types is DNA methylation (DNAm) — where a methyl group is added to cytosine. We will focus on DNAm in the canonical CpG context where cytosine (C) is followed by guanine (G). High levels of DNAm in promoter regions are associated with suppressed gene expression whereas increased methylation in gene body regions tends to have the opposite effect (Suzuki and Bird, 2008).

Reduced representation bisulfite sequencing (RRBS) is a sequencing technique that measures DNAm (Gu et al., 2011). The experimental protocol consists of treating DNA with bisulphite which converts unmethylated cytosines into uracils. During the amplification process uracils are converted into thymines. After sequencing and comparison to a reference genome, unconverted CpGs are identified as unmethylated and vice versa. The RRBS technique does not sequence the entire genome, but rather regions of the genome that are enriched for CpGs. This naturally splits the genome into distinct loci which can be analysed separately.

Conventional analyses of DNAm have focused on the average DNAm level per CpG site. This is obtained by examining all of the sequencing reads which contain a given CpG and simply counting how many times it is methylated. This type of analysis, however, fails to take into account the full methylation pattern at a given locus which can be observed by looking at all contiguous CpGs along a sequencing read. If there are $d$ CpG sites on one read then there are $2^{d}$ possible methylation patterns, which are called epialleles (Richards, 2006). Sequencing reads that cover the same $d$ CpG sites can be compared, and the frequency of distinct epialleles that are present can be calculated. Since each DNA fragment comes from a different cell (more precisely a different allele) this provides a snapshot of how many distinct cellular subpopulations are present within the sample. The additional information acquired from contiguous CpG sites on sequencing reads is not present using array-based platforms. It is becoming clear that leveraging this extra information offers potential insights into the epigenetic landscape that would otherwise be missed (Li et al., 2014; Lin et al., 2015; He et al., 2013).

If multiple samples are taken from the same tumour then each sample can be analysed to see which epialleles are present, and in what proportion, at a given loci. By tracing the presence and absence of different epialleles across different regions of the tumour and matched normal tissue it is possible to reconstruct the evolutionary history of the tumour regions, and to probe for significant differences between normal and tumour tissue. Moreover, the diversity of epialleles within the tumour provides a measure of overall epigenetic heterogeneity.

The acquisition of tumour samples may result in a mixture of both tumour and normal tissue. By comparing the expression of epialleles within the tumour samples and matched normal tissue it is possible to estimate the sample purity — the proportion of the sample which is tumour tissue. Furthermore, it is possible to decontaminate the tumour samples by effectively ‘subtracting’ that component of the epiallele profile which can be attributed to the contaminating normal tissue. An analysis of differential epiallele expression and phylogenetics can be conducted after decontamination.

We present a Bayesian statistical model to infer which epialleles are present at a given locus. The model infers the epialleles that are present and which epiallele each observed sequencing read corresponds to. One hyperparameter controls the level of noise in the model (which represents errors due to bisulfite conversion, PCR amplification, and sequencing) and this is also inferred from the data. Finally, the total number of distinct epialleles is inferred. This final step is a model selection problem and we use the Akaike Information Criterion to avoid overfitting the model. The Bayesian approach allows the quantification of uncertainty regarding the model parameters. In particular, there may be some ambiguity as to which epiallele a certain observed read corresponds to (if some epialleles are very similar to each other for instance). This uncertainty is incorporated into the epiallele distribution by averaging over the appropriate model parameters with respect to the corresponding posterior density.

1.1 Related work

The additional information garnered from adjacent CpGs can be used to define a measure of variability or heterogeneity within a biological sample. The concept of ‘epipolymorphism’, for instance, has been proposed by Landan et al. (2012). Xie et al. (2011) define a measure of ‘methylation entropy’ based on the Shannon entropy and Landau et al. (2014) developed the concept of ‘proportion of discordant reads’.

The term allele-specific methylation has also been used to refer to epialleles. Statistical models have been developed by Peng and Ecker (2012); Fang et al. (2012); Wu et al. (2015) to identify epialleles at a given locus and which epiallele each observed read originated from. These models can infer multiple epialleles but in applications only two epialleles have been assumed. An algorithm to estimate tumour purity and deconvolve the epigenomes of tumour and normal tissue uses a very similar statistical model (Zheng et al., 2014).

Li et al. (2014) compare the epiallele distribution at two disease stages using a ‘composition entropy difference calculation’. They identify loci with substantial shifts in epiallele composition. They confine their analysis to epialleles defined by four CpG sites. Lee et al. (2015) used multinomial logistic regression to test for differences in the epiallele distribution between normal and cancer cells. They report performance that is very similar to the method of Li et al. (2014), but do not constrain their approach to four CpGs. In both of these approaches the epialleles are identified from the raw sequencing data, without any inference step to account for experimental noise.

Lin et al. (2015) develop a statistical model that explicitly takes into account measurement noise due to bisulfite conversion and sequencing errors. The model allows identification of ‘spurious’ epialleles that are due to measurement error (spurious epialleles will tend to have low counts and be very similar to a dominant epiallele). Noise parameters are manually estimated from experimental data, and missing data are not facilitated by their model.

In summary, an adequate epiallele analysis of DNAm sequencing data should have the following features. It should answer the basic research question of whether there is a difference in the epiallele composition between two or more groups of samples — and identify the loci at which there are significant differences. Ideally, some measures need to be taken to avoid spurious epiallele detection due to experimental noise. In addition, an analysis method will generally need to accommodate variable sequencing depth per loci, a variable number of contiguous CpGs per sequencing read, and missing data. Missing data can arise from partially overlapping reads or gaps in a read due to non-overlapping paired-end sequencing protocols.

In addition to the above features, our Bayesian approach automatically infers all model parameters (including the total number of epialleles) from the observed data. Ambiguity in model parameters is explicitly incorporated in our analysis by averaging over the appropriate Bayesian posterior density (descried in detail below). We have applied our method to data from multiple tumour regions and matched normal tissue. We have developed a protocol for estimating the tumour sample purity and consequently decontaminating the inferred epiallele profiles. Although we have focused on multi-region tumour sampling our method could be applied to a single sample also.

2 Methods

Sequencing reads are aligned to the reference genome and organised into different genomic loci. A locus is a region of the genome containing $d$ CpG sites ($d$ can take different values to each locus). Due to the nature of RRBS data the sequencing reads naturally tend to form non-overlapping loci. Some additional steps were taken to modify loci in order to control the amount of missing data per locus. See Supplementary Material A for full details.

Let $N$ denote the number of sequencing reads at a given locus. To keep our notation compact we will avoid indexing each locus and what follows here is applicable to any locus of the genome. A sequencing read is represented by a $d$-dimensional vector $\boldsymbol{\mathbf{y}}_{i}\in\{0,1\}^{d}$ where $i=1,\ldots,N$ and 0 and 1 correspond to unmethylated and methylated CpG sites respectively. An example is plotted in Figure 1 (a). It is assumed that each observed read can be attributed to one of $Q$ epialleles $\boldsymbol{\mathbf{x}}_{q}$ with $q=1,\ldots,Q$ and $Q\leq N$. The parameter $w_{i}\in(1,\ldots,Q)$ specifies which epiallele read $\boldsymbol{\mathbf{y}}_{i}$ originated from. The observed methylation status of each CpG may differ from the corresponding epiallele status with probability $\epsilon\in[0,1/2]$. Supposing $w_{i}=q$ we can therefore write $p(\boldsymbol{\mathbf{y}}_{i}|\boldsymbol{\mathbf{x}}_{q},\epsilon,Q)=\prod_{\mu=1}^{d}p(y_{i\mu}|x_{q\mu},\epsilon,Q)$ where

[TABLE]

The epialleles are analogous to latent variables in a latent variable model. Our goal is to infer the quantities $\boldsymbol{\mathbf{X}}=(\boldsymbol{\mathbf{x}}_{1},\ldots,\boldsymbol{\mathbf{x}}_{Q})$ and $\boldsymbol{\mathbf{w}}=(w_{1},\ldots,w_{N})$ as well as the hyperparameter $\epsilon$ and the number of epialleles $Q$ from the observed data $\boldsymbol{\mathbf{Y}}=(\boldsymbol{\mathbf{y}}_{1},\ldots,\boldsymbol{\mathbf{y}}_{N})$. Using Bayes’ theorem the posterior over the unknown quantities is

[TABLE]

where the likelihood is

[TABLE]

The delta function is defined by $\delta_{xy}=1$ if $x=y$ and $\delta_{xy}=0$ otherwise. The marginal density $p(\boldsymbol{\mathbf{Y}}|Q)=\sum_{\boldsymbol{\mathbf{X}}^{\prime}}\sum_{\boldsymbol{\mathbf{w}}^{\prime}}\int\text{d}\epsilon^{\prime}\,p(\boldsymbol{\mathbf{Y}}|\boldsymbol{\mathbf{X}},\boldsymbol{\mathbf{w}},\epsilon,Q)p(\boldsymbol{\mathbf{X}}|Q)p(\boldsymbol{\mathbf{w}}|Q)$ serves to normalise the posterior density where the summation is over all possible values of $\boldsymbol{\mathbf{X}}$ and $\boldsymbol{\mathbf{w}}$. We will use maximum entropy priors which are uniform densities over the $2^{Qd}$ possible epiallele configurations $\boldsymbol{\mathbf{X}}$ and $Q^{N}$ possible values of $\boldsymbol{\mathbf{w}}$.

2.1 Bayesian inference

For fixed $\boldsymbol{\mathbf{X}}$, $\epsilon$, and $Q$, the maximum a posteriori (MAP) estimate for $\boldsymbol{\mathbf{w}}$ is given by attributing each read $\boldsymbol{\mathbf{y}}_{i}$ to the epiallele that is most similar to it. That is,

[TABLE]

Next we wish to obtain the MAP estimate for $x_{q\mu}$ for fixed $\boldsymbol{\mathbf{w}}$, $\epsilon$ and $Q$. Let $N_{1}$ denote the total number of methylated CpGs at site $\mu$ in observed reads that have been attributed to epiallele $q$. That is, $N_{1}=\sum_{i}y_{i\mu}$ where the sum is restricted to indices for which $w_{i}=q$. Similarly, $N_{0}$ is the total number of unmethylated CpGs at site $\mu$ in reads stemming from epiallele $q$. It is straightforward to show that the MAP estimate is

[TABLE]

We now define the total matches at a given locus as $\alpha_{1}=\sum_{i,\mu}\delta_{y_{i\mu},x_{w_{i}\mu}}$ and mismatches as $\alpha_{0}=\sum_{i,\mu}1-\delta_{y_{i\mu},x_{w_{i}\mu}}$. It can be shown (see Supplementary Material) that the MAP estimate for $\epsilon$ is

[TABLE]

which is simply the proportion of observed CpGs at that locus that differ from the underlying epialleles. Some values of $y_{i\mu}$ may be missing and these are handled by simply omitting them from sums and products over $i$ and $\mu$.

2.1.1 Algorithm

Note that the MAP estimates $\boldsymbol{\mathbf{w}}^{*}$ and $\boldsymbol{\mathbf{X}}^{*}$ are independent of $\epsilon$. Given a set of observed data $\boldsymbol{\mathbf{Y}}$ the first task is to determine optimal values for $\boldsymbol{\mathbf{w}}$ and $\boldsymbol{\mathbf{X}}$. This is done according to the following algorithm:

Initialise $\boldsymbol{\mathbf{w}}$ by using hierarchical clustering to group the observed reads $\boldsymbol{\mathbf{Y}}$ into $Q$ groups. The hamming distance (the proportion of CpGs that differ between two sequencing reads) is used as a distance measure. 2. 2.

Compute $\boldsymbol{\mathbf{X}}$ according to (5) using the current estimate of $\boldsymbol{\mathbf{w}}$. 3. 3.

Compute $\boldsymbol{\mathbf{w}}$ according to (4) using the current estimate of $\boldsymbol{\mathbf{X}}$. 4. 4.

Repeat steps 2 and 3 until $\boldsymbol{\mathbf{w}}$ and $\boldsymbol{\mathbf{X}}$ converge to a steady solution (typically two or three iterations).

Denote the final parameter values as $\hat{\boldsymbol{\mathbf{w}}}$ and $\hat{\boldsymbol{\mathbf{X}}}$. The value for $\hat{\epsilon}$ is then given by (6).

2.1.2 Model selection

In principle, the marginal density $p(\boldsymbol{\mathbf{Y}}|Q)$ could be used to compare models with different values of $Q$. In practice, however, $p(\boldsymbol{\mathbf{Y}}|Q)$ is analytically intractable. Instead we use the Akaike information criterion (AIC) (Akaike, 1998) in order to select the optimal number of epialleles

[TABLE]

where $\hat{Q}=\text{argmin}_{Q}\text{AIC}(Q)$. For a model with $Q$ epialleles the $Qd$ parameters that make up the matrix $\boldsymbol{\mathbf{X}}$ are regarded as free parameters. The term $2Qd$ penalises more complex models (i.e. models with larger $Q$). A more complex model will only be selected if the evidence from the data is sufficiently strong to overcome the penalty term.

2.1.3 Marginalisation of $\boldsymbol{\mathbf{w}}$

Finally, it may not be completely clear which epiallele an observed read should be attributed to (there could be several epialleles an equal edit distance away). This ambiguity manifests itself as the uncertainty surrounding the parameter $w_{i}$. The Bayesian approach allows this uncertainty to be incorporated into our analysis. The marginal density over $w_{i}$ is given by fixing all other parameters to their MAP values

[TABLE]

where $\hat{\boldsymbol{\mathbf{w}}}_{-i}$ is a $(d-1)$-dimensional vector obtained from $\hat{\boldsymbol{\mathbf{w}}}$ by removing element $i$. At the given locus in question the proportion of observed reads originating from epiallele $q$ is given by

[TABLE]

The quantity $\boldsymbol{\phi}=(\phi_{1},\ldots,\phi_{\hat{Q}})$ specifies the distribution of epialleles within that locus. An example is given in Figure 1 (d).

2.2 Application to multi-region tumour sampling

We will now describe our analysis protocol. In our application we are considering sequencing data from multiple regions of the same tumour. The number of distinct epialleles present at a particular locus is determined by pooling sequencing reads from all tissue samples (tumour and normal) in order to boost statistical power. Suppose there are $s=1,\ldots,S$ tumour samples with $N_{s}$ reads per sample (at a given locus). The total number of reads in the pool is now $N=\sum_{s}N_{s}$. Using the pooled reads a model is fitted as described above. The vector $\hat{\boldsymbol{\mathbf{w}}}$ defines which epiallele each sequencing read originated from. The distribution of epialleles within region $s$ is given by

[TABLE]

where $I_{s}$ is the set of indices of reads belonging to sample $s$. The vectors $\boldsymbol{\phi}^{s}$ serve to characterise each sample in terms of their epiallele distributions.

2.2.1 Estimation of sample purity

Suppose $\hat{Q}$ epialleles are inferred at particular locus of a particular tumour sample (for the sake of compactness we will not index the loci or samples). The locus is characterised by $\boldsymbol{\phi}$, the inferred probability distribution over the $\hat{Q}$ epialleles. If the tumour sample is contaminated with normal tissue then we can write

[TABLE]

where $\rho\in[0,1]$ is the proportion of observed tissue that comes from the tumour (the sample ‘purity’), and $\boldsymbol{\mathbf{t}}$ and $\boldsymbol{\mathbf{n}}$ are the epiallele distributions in the tumour and normal tissues respectively (at the particular locus in question). We can estimate $\boldsymbol{\phi}$ and $\boldsymbol{\mathbf{n}}$ from the observed data at a particular locus. Estimation of both $\rho$ and $\boldsymbol{\mathbf{t}}$ requires solving the $\hat{Q}$ equations in (11) for $\hat{Q}+1$ variables which generally is not possible. However, the quantity

[TABLE]

can be computed at each locus of the observed tissue sample. If we substitute (11) into (12) we can see that $\xi$ takes a minimum value of 0 when $\boldsymbol{\mathbf{t}}=\boldsymbol{\mathbf{n}}$. At a locus in which the tumour and normal tissues have a completely different epiallele composition then we say that if $t_{\mu}>0$ then $n_{\mu}=0$ and if $n_{\mu}>0$ then $t_{\mu}=0$ for $\mu=1,\ldots,\hat{Q}$. It is straightforward to show that if this is the case then $\xi=\rho$ and that this is the maximum value $\xi$ can take.

We therefore expect that $\xi$ will take values in the range $[0,\rho]$ when computed across all loci of the observed sample. If we plot the empirical density of $\xi$ values the parameter $\rho$ can be estimated from the maximum value of $\xi$. Since $\boldsymbol{\phi}$ and $\boldsymbol{\mathbf{n}}$ are estimated from finite data samples we expect the distribution of $\xi$ to be ‘smoothed’ by sampling noise. This is precisely what we observe in practice. An example of the empirical density of $\xi$ is plotted in Figure 11.

2.2.2 Decontamination of normal tissue

Finally, we note that once estimates of $\rho$ have been obtained we can calculate the ‘decontaminated’ tumour epiallele profiles at each locus according to

[TABLE]

We have used the notation $\hat{t}_{\mu}$ to emphasise that this is an estimate of the tumuor epiallele distribution. Due to the fact that $\boldsymbol{\phi}$, $\boldsymbol{\mathbf{n}}$ and $\rho$ are estimated from finite data samples it is possible that $\hat{t}_{\mu}$ can take values outside $[0,1]$. Any cases where $\hat{t}_{\mu}<0$ are set to 0 and any cases where $\hat{t}_{\mu}>1$ are set to 1.

A conventional anlaysis of DNAm sequencing data will typically ‘call’ a methylation level at each CpG site by computing the proportion of reads on which a CpG is observed in a methylated state. Using our method a methylation level for each CpG site can readily be computed after decontamination of normal tissue and used in existing analysis pipelines.

2.2.3 Construction of a phylogenetic tree

Using the decontaminated representation of a sample $\hat{\boldsymbol{\mathbf{t}}}_{s}$ the euclidean distance between $\hat{\boldsymbol{\mathbf{t}}}_{s}$ and $\hat{\boldsymbol{\mathbf{t}}}_{s^{\prime}}$ can be used as a distance measure between samples $s$ and $s^{\prime}$. Each locus provides a distance matrix that depends on the distribution of epialleles at that particular locus. To obtain an overall distance matrix we average over distance matrices from all loci. Any distance based phylogenetic inference method can subsequently be used to construct a phylogenetic tree. We used the ‘fastme.bal’ function as part of the ‘ape’ R package (Paradis et al., 2004).

3 Results

3.1 Simulations

Simulations of a single locus were performed to study what effect the number of CpGs, $d$, the number of sequencing reads, $N$, and the noise level, $\epsilon$, have on our ability to correctly detect the underlying epialleles. The simulated reads were noise corrupted versions of three distinct randomly generated epialleles, and on average each epiallele corresponded to one third of the observed reads. To assess model performance we counted the proportion of observed reads that were attributed to their correct underlying epiallele (which requires both inference of the correct epialleles and attribution to the correct epiallele). For every value of the parameters results were averaged over 100 simulations.

We found that $N=100$ and $d=6$ gave a success rate of approximately 95% at a 5% noise level. These values were used to guide the selection of viable loci in subsequent analyses of experimental data. Dropping to $N=50$ gave a performance of just over 90% (Supplementary Figure 3). Sequencing depth beyond $N=100$ did not yield any additional performance gain. The performance saturates at 100% for $d>15$ (Supplementary Figure 4). Since the number of possible epialleles is $2^{d}$ a larger $d$ will typically make it easier to resolve distinct epialleles. Additionally, since the underlying epialleles are randomly generated it is possible that some may be within one edit distance from each other, making it difficult for the model to distinguish between very similar epialleles and noise when $d$ is small. Performance was observed to decrease sharply for increasing noise levels (Supplementary Figure 5).

3.2 Cell line data: detection of low frequency epialleles

In order to test whether our statistical methods could detect low frequency epialleles in practice we mixed a fully unmethylated and fully methylated cell line in a 9:1 ratio prior to sequencing. Loci with six or more CpGs and 50 or more reads were identified. Within these loci 6.3% of observed CpGs were methylated overall. The two cell lines were sequenced separately and we found that the fully methylated and unmethylated cells were in fact $97.3\%$ and $3.8\%$ methylated respectively.

The Bayesian model was used to detect the presence of epialleles at each loci. We found that 5.2% of methylated CpGs were attributed to methylated epialleles (defined as epialleles with $\geq 50\%$ methylation). The mean noise level was inferred as 1.1%. This suggests that the majority of methylation is correctly identified as corresponding to a methylated profile and therefore our method is capable of resolving a distinct low frequency cellular subpopulation.

3.3 Multi-region tumour sampling case study

Our case study data consisted of seven tissue samples from a single lung tumour (CRUK00620) along with one matched normal tissue sample. These tissue samples were acquired as part of the larger TRACERx study (Jamal-Hanjani et al., 2017). The raw sequencing data were trimmed and aligned to a reference genome. Sequencing reads were subsequently organised into distinct genomic loci as described in the Supplementary Material. We demanded that no more than 25% of data were missing per locus (due to partially overlapping paired-end reads or reads not covering the whole locus). Any data from chromosomes X and Y were discarded. At each locus $\hat{Q}$ epialleles are inferred and any epialleles that accounted for less than 5% of observed reads were discarded prior to the computation of $\boldsymbol{\phi}_{s}$ for $s=1,\ldots,S$. This was done in order to focus on the dominant shifts in epiallele profiles and to minimise the risk of inferring spurious epialleles.

In order to compare the distribution of epialleles within different tumour samples it was necessary to identify all of the loci which occurred in two or more samples. That is, the loci themselves must ‘match up’ between tumour samples in order for a comparison to be made (partially overlapping loci were permitted provided they met the minimum number of non-missing CpG requirements). Only loci with a median read depth $\geq 100$ across normal and tumour tissue samples and six or more CpGs were considered. A total of 39,940 loci were analysed out of which 73% were found to contain a single epiallele, 13% contained two, 7% contained three, 4% contained four, and 3% had five our more (up to a maximum of thirteen).

3.3.1 Comparison of epiallele distribution throughout the tumour

At each locus the Bayesian model is used to infer the epialleles present, the total number of epialleles, and which epialleles each observed sequence came from. An example locus with seven CpGs from chromosome one is presented in Figure 2. At this locus five distinct epialleles were detected. Both the observed and decontaminated profiles are shown. The normal tissue is predominantly composed of methylated epialleles whereas the tumour samples have a greater proportion of less methylated epialleles. This suggests that within the tumour there exist cellular subpopulations that are undergoing a transition from a methylated state to an unmethylated one.

In order to understand shifts in epiallele frequency at a global level we plotted a heatmap of the top 200 most variable epialleles in Figure 4 (a) and (c). Both the observed and decontaminated epiallele profiles were used. Tumour samples are characterised by both a loss and gain of numerous epialleles when compared to the normal tissue sample. The variability in epiallele expression throughout different parts of the tumours suggests that a substantial level of tumour heterogeneity exists at the epigenetic level. Note that in the contaminated samples 71 out of the 200 epialleles were located on CpG islands, and 54 were located on a CpG shore (defined as 2 kilobases either side of an island). In the decontaminated version 124 epialleles were located on an island and 38 on a shore.

3.3.2 Estimation of sample purity

The sample purities were estimated as described in Section 2.2.1. An example of the empirical density of $\xi$ within tumour region 2 is plotted in Figure 11. From the location of the rightmost maximum we estimate $\rho=0.535$. Plots for all tumour regions are given in Supplementary Figure 6. Estimates of purity for the seven tumour samples are given in Table 1. For tumour region 6 the rightmost maxima was not visible presumably due to very low tumour purity. The purity estimates are compared to estimates obtained from an analysis of exome data from the same tissue samples performed independently in Jamal-Hanjani et al. (2017).

3.3.3 Inference of a phylogenetic tree

Phylogenetic trees were generated as described in Section 2.2.3. The trees for both contaminated and decontaminated samples are plotted in Figure 4 (b) and (d). The structure of the contaminated tree is dominated by the sample purities, with low purity samples clustering together. The decontaminated tree has a totally different structure and this is broadly similar to a phylogenetic tree obtained from from a separate genetic analysis of the same patient and shown in Supplementary Figure 7.

3.3.4 Quantification of epigenetic disorder

The Shannon entropy provides a measure of how disordered a random variable is. In particular, the entropy of the epiallele distribution $\boldsymbol{\phi}_{s}$ quantifies how disordered or heterogeneous each loci is in sample $s$. In Figure 5 box plots summarise the distribution of entropies across tumour and normal tissues (without decontamination). The tumour tissue samples have a substantially elevated entropy in comparison to the normal tissue.

4 Discussion

Analysis of epialleles allows for a deeper interrogation of the underlying biology than a pointwise examination of CpG methylation states. Tracing the patterns of DNA methylation along epialleles allows one to tease apart different cellular subpopulations and acquire a richer quantification of heterogeneity and disorder that would not be possible by looking at individual CpG sites. In particular, the distribution of epialleles throughout a tumour can shed light on the evolutionary history of the tumour.

Our analysis protocol specifically pools sequencing reads from multiple tissue samples in order to leverage greater statistical power in epiallele detection. Our Bayesian approach will automatically detect the number of epialleles present, and infer what the methylation pattern of those epialleles are. One strength of the Bayesian approach is that it provides a framework for averaging over uncertainty in model parameters. If there is uncertainty as to which epiallele an observed sequencing read may have originated from, then a natural solution is to average over that uncertainty by marginalising over the appropriate posterior distribution. In addition to the above features our model can easily accommodate missing data and can handle an arbitrary sequencing depth and number of CpG sites per locus. Furthermore, by comparing the distribution of epialleles within normal and tumour tissue samples it is possible to estimate the purity of each sample and to subsequently decontaminate them. Methylation levels at each CpG site can be extracted from the decontaminated samples and subsequently used in standard analysis pipelines.

Tracking the presence or absence of epialleles throughout the tumour opens up an additional layer of complexity beyond that of conventional methylation analyses. Pointwise methylation analysis protocols typically average over sequencing reads – to ‘call’ the methylation status at single CpGs – that potentially come from a diverse and heterogenous population of cells. Detecting which epialleles are present allows one to distinguish between these cellular subpopulations and identify tumour subclones that are defined by distinct epialleles. One can then probe changes between normal and cancerous tissue at a finer resolution. As we have demonstrated here, studying epiallele frequencies in different parts of the tumour reveals the evolutionary history of the tumour and allows a phylogenetic tree to be constructed. A measure of disorder or heterogeneity inside the tumour can be obtained through measures such as Shannon’s entropy.

5 Conclusion

Understanding tumour heterogeneity is an important step towards understanding why certain therapies fail and why resistance to treatment can emerge. Subclonal populations of treatment-resistant cells can persist after treatment even if they only account for a small fraction of the original tumour. Epigenetic diversity within the tumour may play an important role in tumour evolution alongside genetic variability. It is increasingly recognised that for DNA methylation sequencing data studying the patterns of methylation along the genome – ‘epialleles’ – can provide greater insight into the underlying dynamics of epigenetic regulation than a conventional pointwise analysis.

We have exploited this opportunity to study the distribution of epialleles throughout a tumour by performing reduced representation bisulfite sequencing on seven regions of the same tumour and one matched normal tissue sample. Our new Bayesian approach infers which epialleles are present at a given locus. A comparison of the frequency of different epialleles across the tumour and normal tissue highlights changes between normal and cancerous tissue and allows the extraction of a phylogenetic history. The concept of entropy can be used as a measure of global disorder within the tumour. Our method can be applied more generally to any type of DNAm sequencing data.

Future work will focus on larger scale studies of multiple patients with multi-region tumour sampling in order to probe for systematic alterations in epiallele expression between normal and cancerous tissue. Previously, measures of epigenetic disorder were found to be associated with clinical outcome and it will be interesting to see if quantification of disorder at the level of epialleles will provide a more refined measure of tumour aggressiveness. Ultimately, it is hoped that a clearer elucidation of epigenetic dynamics will complement our genetic knowledge of cancer and provide a more comprehensive understanding of the disease.

Acknowledgments

The authors would like to thank Pawan Dhami (UCL Cancer Institute Genomics Core Facility) for sequencing support.

Funding

JB was supported by the CRUK & EPSRC Comprehensive Cancer Imaging Centre at King’s College London and University College London jointly funded by Cancer Research UK and the EPSRC; AF by the MRC (MR/M025411/1); JH by the UCL Cancer Institute Research Trust; MT by the People Programme (Marie Curie Actions) of the EU Seventh Framework Programme (FP7/2007-2013/608765) and the Danish Council for Strategic Research (1309-00006B); GAW is funded by Cancer Research UK (grant number C11496/A17786); SB by NIHR-BRC (BRC275/CN/SB/101330) and the Wellcome Trust (99148); CS is Royal Society Napier Research Professor; This work was supported by the Francis Crick Institute which receives its core funding from Cancer Research UK (FCI01), the UK Medical Research Council (FC001169), and the Wellcome Trust (FC001169); by the UK Medical Research Council (MR/FC001169/1); CS is funded by Cancer Research UK (TRACERx), the CRUK Lung Cancer Centre of Excellence, Stand Up 2 Cancer (SU2C), the Rosetrees Trust, NovoNordisk Foundation (ID 16584), the Prostate Cancer Foundation, the Breast Cancer Research Foundation (BCRF), the European Research Council (THESEUS) and Marie Curie Network PloidyNet. Support was also provided to CS by the National Institute for Health Research, the University College London Hospitals Biomedical Research Centre, and the Cancer Research UK University College London Experimental Cancer Medicine Centre.

Appendix A Extraction of viable loci

In Figure 6 is an example of an observed locus before any preprocessing steps have been taken. The locus is defined naturally since it does not overlap with any other observed sequencing reads. Note that several sequencing reads only partially overlap and many do not overlap at all. In addition, several CpG measurements are missing from the middle of some reads due to the paired-end sequencing protocol that was used to generate the data (sometimes the paired ends may not span the full length of the DNA fragment).

Intuitively this locus should be split into two loci as the two blocks of reads overlap by a single CpG which is not enough to phase the inferred epialleles. We implemented the following algorithm to split the observed sequencing reads into sensible loci.

Specify the minimum number of contiguous CpGs $d_{min}$ and the minimum number of reads $N_{min}$ required in order for a locus to be admissible (we used $d_{min}=6$ and $N_{min}=100$ in practice). 2. 2.

Specify the maximum proportion of missing data that is allowed (we used 25% in practice). 3. 3.

If the observed locus contains more than 25% missing data move to the next step, otherwise skip to step 4.

• 3(a).

Represent the observed reads as a matrix with all non-missing measurements equal to 1 and all missing values represented with 0. Use hierarchical clustering to split the reads into two groups using the hamming distance.

• 3(b).

If the two groups contain less than 25% of the missing data then proceed to the next step. Otherwise repeat the clustering above with three groups and so forth. 4. 4.

Discard any loci that fail to meet the minimum values of $d$ and $N$.

We demanded that, in addition to the loci as a whole, each CpG site should not contain more than 25% missing data. Any CpGs that failed to meet this constraint were discarded. This helped to trim low-coverage CpGs from at the edges of observed loci (for example, the very leftmost CpG in Figure 12). The resulting loci after preprocessing are shown in Figure 7.

Appendix B MAP estimate for $\epsilon$

Given $\boldsymbol{\mathbf{X}}$, $\boldsymbol{\mathbf{w}}$, $Q$ and uniform priors $p(\boldsymbol{\mathbf{X}}|Q)$ and $p(\boldsymbol{\mathbf{w}}|Q)$ the MAP estimate for the hyperparameter $\epsilon$ is given by maximising the log of the posterior (2):

[TABLE]

Note that on the second line that if read $\boldsymbol{\mathbf{y}}_{i}$ originates from epiallele $q$ then $w_{i}=q$. Recall that $\alpha_{1}=\sum_{i,\mu}\delta_{y_{i\mu},x_{w_{i}\mu}}$ and $\alpha_{0}=\sum_{i,\mu}1-\delta_{y_{i\mu},x_{w_{i}\mu}}$ denote the total number of matches and mismatches between the observed reads $\boldsymbol{\mathbf{y}}$ and the corresponding epialleles $\boldsymbol{\mathbf{x}}$ at this particular loci. Solving $\text{d}\mathcal{L}/\text{d}\epsilon=0$ yields $\epsilon=\alpha_{0}/(\alpha_{0}+\alpha_{1})$.

Appendix C Simulation results

Appendix D Purity estimation

Appendix E Supplementary figures

Appendix F Experimental datasets

The RRBS data are available at the European Nucleotide Archive under accession numbers ERS1546024, ERS1546025 and ERS1546026.

Appendix G TRACERx consortium members

The TRACERx study (Clinicaltrials.gov no: NCT01888601) is sponsored by University College London (UCL/12/0279) and has been approved by an independent Research Ethics Committee (13/LO/1546). TRACER is funded by Cancer Research UK (grant number C11496/A17786) and coordinated through the Cancer Research UK & UCL Cancer Trials Centre.

Consortium members

Charles Swanton1,2,5, Mariam Jamal-Hanjani1, Selvaraju Veeriah1, Seema Shafi1, Justyna Czyzewska-Khan1, Diana Johnson1, Joanne Laycock1, Leticia Bosshard-Carter1, Gerald Goh1, Rachel Rosenthal1, Pat Gorman1, Nirupa Murugaesu1, Robert E Hynds1,3, Gareth Wilson1,2, Nicolai J Birkbak1,2, Thomas B K Watkins2, Nicholas McGranahan1,2, Stuart Horswell2, Richard Mitter2, Mickael Escudero2, Aengus Stewart2, Peter Van Loo2, Andrew Rowan2, Hang Xu2, Samra Turajlic2,4, Crispin Hiley2, Christopher Abbosh1, Jacki Goldman2, Richard Kevin Stone2, Tamara Denner2, Nik Matthews2, Greg Elgar2, Sophia Ward2, Jennifer Biggs2, Marta Costa2, Sharmin Begum2, Ben Phillimore2, Tim Chambers2, Emma Nye2, Sofia Graca2, Maise Al Bakir2, Kroopa Joshi1, Andrew Furness1, Assma Ben Aissa1, Yien Ning Sophia Wong1, Andy Georgiou1,Sergio Quezada1, John A Hartley1, Helen L Lowe1, Javier Herrero1, David Lawrence5, Martin Hayward5, Nikolaos Panagiotopoulos5, Shyam Kolvekar5, Mary Falzon5, Elaine Borg5, Teresa Marafioti5, Celia Simeon5, Gemma Hector5, Amy Smith5, Marie Aranda5, Marco Novelli5, Dahmane Oukrif5, Sam M Janes5, Ricky Thakrar5, Martin Forster5, Tanya Ahmad5, Siow Ming Lee5, Dionysis Papadatos-Pastos5, Dawn Carnell5, Ruheena Mendes5, Jeremy George5, Neal Navani5, Asia Ahmed5, Magali Taylor5, Junaid Choudhary5, Yvonne Summers6, Raffaele Califano6, Paul Taylor6, Rajesh Shah6, Piotr Krysiak6, Kendadai Rammohan6, Eustace Fontaine6, Richard Booton6, Matthew Evison6, Phil Crosbie6, Stuart Moss6, Faiza Idries6, Leena Joseph6, Paul Bishop6, Anshuman Chaturved6, Anne Marie Quinn6, Helen Doran6, Angela leek7, Phil Harrison7, Katrina Moore7, Rachael Waddington7, Juliette Novasio7, Fiona Blackhall8, Jane Rogan7, Elaine Smith6, Caroline Dive9, Jonathan Tugwood9, Ged Brady9, Dominic G Rothwell9, Francesca Chemi9, Jackie Pierce9, Sakshi Gulati9, Babu Naidu10, Gerald Langman10, Simon Trotter10, Mary Bellamy10, Hollie Bancroft10, Amy Kerr10, Salma Kadiri, 10, Joanne Webb10, Gary Middleton10, Madava Djearaman10, Dean Fennell11, Jacqui A Shaw11, John Le Quesne11, David Moore11, Apostolos Nakas12, Sridhar Rathinam12, William Monteiro13, Hilary Marshall13, Louise Nelson12, Jonathan Bennett12, Joan Riley12, Lindsay Primrose12, Luke Martinson12, Girija Anand14, Sajid Khan15, Anita Amadi16, Marianne Nicolson17, Keith Kerr17, Shirley Palmer17, Hardy Remmen17, Joy Miller17, Keith Buchan17, Mahendran Chetty17, Lesley Gomersall17, Jason Lester18, Alison Edwards18, Fiona Morgan19, Haydn Adams19, Helen Davies19, Malgorzata Kornaszewska20, Richard Attanoos21, Sara Lock22, Azmina Verjee22, Mairead MacKenzie23, Maggie Wilcox23, Harriet Bell24, Natasha Iles24, Allan Hackshaw24, Yenting Ngai24, Sean Smith24, Nicole Gower24, Christian Ottensmeier25, Serena Chee25, Benjamin Johnson25, Aiman Alzetani25, Emily Shaw25, Eric Lim26, Paulo De Sousa26, Monica Tavares Barbosa26, Alex Bowman26, Simon Jordan26, Alexandra Rice26, Hilgardt Raubenheimer26, Chiara Proli26, Maria Elena Cufari26, John Carlo Ronquillo26, Angela Kwayie26, Harshil Bhayani26, Morag Hamilton26, Yusura Bakar26, Natalie Mensah26, Lyn Ambrose26, Anand Devaraj26, Silviu Buderi26, Jonathan Finch26, Leire Azcarate26, Hema Chavan26, Sophie Green26, Hillaria Mashinga26, Andrew G Nicholson26, 27, Kelvin Lau28, Michael Sheaff28, Peter Schmid28, John Conibear28, Veni Ezhil29, Babikir Ismail29, Melanie Irvin-sellers29, Vineet Prakash29, Peter Russell30, Teresa Light30, Tracey Horey30, Sarah Danson31, Jonathan Bury31, John Edwards31, Jennifer Hill31, Sue Matthews31, Yota Kitsanta31, Kim Suvarna31, Patricia Fisher31, Allah Dino Keerio31, Michael Shackcloth32, John Gosney32, Pieter Postmus32, Sarah Feeney32, Julius Asante-Siaw32, Tudor Constatin33, Raheleh Salari33, Nicole Sponer33, Ashwini Naik33, Bernhard Zimmermann33, Hugo J.W.L. Aerts34, Stefan Dentro35, Christophe Dessimoz36,37,38.

Affiliations

1. Cancer Research UK Lung Cancer Centre of Excellence, University College London Cancer Institute, United Kingdom

2. The Francis Crick Institute, United Kingdom

3. Lungs for Living, UCL Respiratory, University College London, United Kingdom

4. The Royal Marsden Hospital, United Kingdom

5. University College London Hospitals NHS Foundation Trust, United Kingdom

6. University Hospital of South Manchester, United Kingdom

7. Manchester Cancer Research Centre Biobank, United Kingdom

8. Christie NHS Foundation Trust, Manchester, United Kingdom

9. Cancer Research UK Manchester Institute, United Kingdom

10. Heart of England NHS Foundation Trust, Birmingham, United Kingdom

11. Cancer Studies and Molecular Medicine, University of Leicester, United kingdom

12. Leicester University Hospitals, United Kingdom

13. National Institute for Health Research Leicester Respiratory Biomedical, Research Unit, United Kingdom

14. North Middlesex Hospital, United Kingdom

15. Royal Free Hospital, United Kingdom

16. Barnet Hospital, United Kingdom

17. Aberdeen Royal Infirmary, United Kingdom

18. Velindre Cancer Centre, Cardiff, Wales, United Kingdom

19. Cardiff & Vale University Health Board, Cardiff, Wales, United Kingdom

20. University Hospital Of Wales Heath Park, Cardiff, Wales, United Kingdom

21. Department of Pathology, University Hospital of Wales and Cardiff University, Heath Park, Cardiff, Wales, United Kingdom

22. The Whittington Hospital NHS Trust, United Kingdom

23. Independent Cancer Patients Voice, United Kingdom

24. Cancer Research UK & UCL Cancer Trials Centre, United Kingdom

25. University Hospital Southampton NHS Foundation Trust, United Kingdom

26. Royal Brompton and Harefield NHS Foundation Trust, United Kingdom

27. National Heart and Lung Institute, Imperial College, United Kingdom

28. Barts Health NHS Trust, United Kingdom

29. Ashford and St. Peter’s Hospitals NHS Foundation Trust, United Kingdom

30. The Princess Alexandra Hospital NHS Trust, United Kingdom

31. Sheffield Teaching Hospitals NHS Foundation Trust, United Kingdom

32. Liverpool Heart and Chest Hospital NHS Foundation Trust, United Kingdom

33. Natera Inc., 201 Industrial Road, Suite 410, San Carlos, CA 94070

34. Dana-Farber Cancer Institute, Brigham & Women’s Hospital, Harvard Medical School, 450 Brookline Ave, JF518, Boston, MA 02115-5450, USA

35. Wellcome Trust Sanger Institute, Hinxton, CB10 1SA, United Kingdom

36. Bioinformatics Group, Department of Computer Science, University College London

37. University of Lausanne

38. Swiss Institute of Bioinformatics