Genomic evolution shapes prostate cancer disease type

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contactDavid C. WedgePh.D. ([email protected]).

Materials availability

There are no tangible materials produced by this study that are available for distribution.

Data and code availability

Sequencing data generated for this study have been deposited in the European Genome-phenome Archive with accession code EGAS00001000262. Processed data and code used in this manuscript is available at https://github.com/woodcockgrp/evotypes_p1/ and via https://doi.org/10.5281/zenodo.10214795.

DNA preparation and DNA sequencing

DNA from frozen tumor tissue and whole blood samples (matched controls) was extracted and quantified using a ds-DNA assay (UK-Quant-iT PicoGreen dsDNA Assay Kit for DNA) following the manufacturer’s instructions with a Fluorescence Microplate Reader (Biotek SynergyHTBiotek). Acceptable DNA had a concentration of at least 50 ng/μl in TE (10mM Tris/1mM EDTA)and displayed an optical density 260/280 (${\text{OD}}_{260}$ / ${\text{OD}}_{280}$) ratio between 1.8 and 2.0. Whole Genome Sequencing (WGS) was performed at IlluminaInc. (Illumina Sequencing FacilitySan DiegoCA USA) or the BGI (Beijing Genome InstituteHong Kong)as described previously,^38,39 to a target depth of 50X for the cancer samples and 30X for matched controls.³⁸ The Burrows-Wheeler Aligner³³ (BWA) was used to align the sequencing data to the GRCh37 reference human genome.

Generation of summary measurements

We generated 123 summary measurements from the WGS data using a number previously published algorithmsso we briefly outline those below. These are grouped into measurements that were generated with similar or related algorithms; default parameters were used unless otherwise stated. The processed data is given alongside the code at https://github.com/woodcockgrp/evotypes_p1/.

Statistics

Prior to the study we predetermined we would use Fisher’s Exact Test for 2x2 contingency tables and Chi-squared test for contingency tables of greater dimensionality and this is applied throughout. Associations between genetic alterations and ARBS clusters was identified using one-tailed Fisher Exact Test with $p<$ 0.05corrected for multiple testing using the False Discovery Rate. Relationships were determined dependent on the variable type: for Bernoulli variablessignificance was determined with Chi-squared test followed by a one-tailed Fisher exact test for each pairwise relationship; one-tailed tests were used as a two-tailed test would not have revealed the direction of the relationship. For continuous variables a Kruskal-Wallace test with Tukey’s HSD was used (adjusted $p<$ 0.05 for all tests). Significance of Depleted groups across countries clustering together was determined using the Approximately Unbiased Multiscale Bootstrap procedure. Associations between evotypes and individual genetic alterations was conducted with a two-tailed Fisher Exact Testcorrected for multiple testing using the False Discovery Rate. The associations with ARBS pairs were established with a one-sided Mann-Whitney U-test with $p<$ 0.05. Statistics associated with the Kaplan-Meier plot were calculated using log rank methodsand significance level was set at 0.05. Cramer’s V statistic was used to determine the strength of the associations in between the cluster assignments. As we only claim an association between patients assigned to MC-A (Metaclusters)the Depleted group (ARBS) and Ordering IIwe combined metaclusters MC-B1 and MC-B2 into one class and the Enriched and Depleted ARBS groups into one class for this comparison.

Unsupervised feature extraction

The summary measurements detailed above form the dataset for further analysis. Howeverit contains a number of different data types (binaryproportionscontinuousinteger counts)it is high dimensional relative to the number of patientsand it undoubtedly contains highly correlatedcooccurring or equivalent events that may confound the analysis. To address this we performed a feature extraction preprocessing step prior to the analysis. As our downstream analysis will be investigating genomic patterns that are indicative of evolutionary behaviorit is critical that the results of these analyses can be easily interpreted. This necessitates methodology where the links between input variables that correspond to the features are identifiable. We therefore opted for a latent feature approach as the basis of our feature extraction as these can provide an interpretable representation of the relationships between the inputs.⁴² Latent feature (or latent variable) analysis provides a way of reformulating the data into a reduced set of features that encapsulate the underlying relationships between the original inputs. The data can be recast in terms of these latent featureswhich is known as the latent feature representationand the downstream analysis performed directly on this.

There have been many latent feature models proposedeach with associated inference methods for the features (a process called feature learning). These included methods such as non-negative matrix factorization,⁴³ Bayesian non-parametric methods⁴⁴ and neural networks.⁴⁵ Howevernone of these were able to fulfill all of our requirements above. We therefore created a bespoke method for feature extraction on this dataset.

Neural networks for feature extraction

We utilized a Restricted Boltzmann Machine⁴⁶ (RBM) neural network as the basis of our feature learning method. We chose to use an RBM as it is extensible to multiple data types^47,48 and can provide interpretable hidden unitswith appropriate modifications.⁴⁹ An RBM is functionally similar to another type of neural network architecture called an autoencoder.⁴⁵ Autoencoders are a class of network types that compress (encode) the data into a transformed representation (the code)and then decompress (decode) in an attempt to reconstruct the original data. A measure of the error between the reconstruction and the original data is used to update the parameters through backpropagation. Typically the code layer contains fewer units than the input/output layers and this bottleneck means that the learning process attempts to compress the information in the dataset into a more a compact representation in the code.

In contrastthe basic RBM unit consists of only two layersknown as the visible and the hidden layers. The RBM is formulated as a probabilistic networkmeaning each unit represents a random variable rather than a fixed value. As suchthe hidden layer performs a similar function to the code layer in the autoencoderalbeit with a probabilistic representation. It has been shown that the RBM is equivalent to the graphical model of factor analysis⁵⁰ and so each hidden unit can be interpreted as a latent feature. Another distinction from the autoencoder formulation is that there is only one weight matrixwhich used to update both the visible and hidden layers. This means that the information on the transformation from visible units (input representation) to the hidden units (feature representation) is encapsulated in this matrix. Hence we also refer to it as the input-feature map.

The restricted Boltzmann machine

The standard RBM formulation⁴⁶ consists of Bernoulli random variables for all visible $\mathbf{v}=\left\{v_i\right\}$ and hidden units $\mathbf{h}=\left\{h_i\right\}$where $v_i,h_j\in \left\{0,1\right\}$with respective biases $\mathbf{a}=\left\{a_i\right\},\mathbf{b}=\left\{b_j\right\};a_i,b_j\in \left(-\infty ,\infty \right)$and a matrix of weights$W;w_{ij}\in \left(-\infty ,\infty \right)$. Training of an RBM is based on minimizing the free-energy of the visible unitsas a low free-energy corresponds to a state where the data is explained well through the model parameterization. Energy-based probability distributions take the form

where $E\left(\mathbf{v},\mathbf{h}\right)$ is the energy function and Z is a normalizing factor. This is the probability of observing the joint $\mathbf{v},\mathbf{h}$ pair. The energy function in an RBM is given as

In this formulation,

which is difficult to calculate due to the number of possible combinations of $\mathbf{v}$ and $\mathbf{h}$.

As we want training to be conducted with respect to the energy at the visible unitswe need to marginalize over $\mathbf{h}$ in Equation 1 to calculate the likelihood of observing the visible unit corresponding to a single data sample ${\mathbf{d}}_k$ from dataset $D=\left\{{\mathbf{d}}_k,k=1,2,\dots ,K\right\}$.

where $\theta \in \left\{\left\{a_i\right\},\left\{b_j\right\},\left\{w_{ij}\right\}\right\}$ is the full parameter set. To simplify notationwe write $\mathcal{L}\left(\theta |\mathbf{v}={\mathbf{d}}_k\right)$ as $\mathcal{L}\left({\mathbf{d}}_k\right)$ with no loss of generality. To perform training through gradient descentwe need to calculate the gradient of the negative log likelihood for each parameter we wish to update$\partial \left(-\log \mathcal{L}\left({\mathbf{d}}_k\right)\right)/\partial \theta$. The partial derivative of the logarithm of Equation 4 takes the form

We then calculate the expected values using the entire training set

which can be used to update the model parameters via gradient descent. The ${\mathbb{E}}_{P\left(\mathbf{h}|D\right)}$ term corresponds to the expected energy state invoked from observing the data samplesand the ${\mathbb{E}}_{P\left(\mathbf{v},\mathbf{h}\right)}$ is the expected energy state of the model configurationsboth contingent on the current model parameters. As suchthey are often called ${\mathbb{E}}_{data}$ and ${\mathbb{E}}_{model}$ respectively. Calculating the partial derivatives with respect to the parameters gives

which are used to construct the update equations

for learning rates ν and η. The ${\mathbb{E}}_{data}$ values can be estimated easily by taking the arithmetic mean.

The ${\mathbb{E}}_{model}$ terms are generally difficult to calculate as they involve summation over all possible configurations of $\mathbf{v}$ and $\mathbf{h}$. An alternative is to perform Gibbs sampling using the conditional probabilities as these are far easier to calculate due to the conditional independence between units in the same layer. We can estimate the conditional probability of values of the hidden layer from the visible layer and vice versa thus

The form of $P\left(h_j|\mathbf{v}\right)$ and $P\left(v_i|\mathbf{h}\right)$ depends on the activation function. This function that inputs the products of the units in one layer and their corresponding weightsand outputs a probability that a unit is active. In this studywe use a logistic sigmoid (or simply “sigmoid”) functionwhich is given by

where x is dependent on the layer we are samplingand so the individual hidden and visible probabilities can be written as

A sample is drawn by setting the corresponding unit to 1 with probability given by the value for $P\left(h_j|\mathbf{v}\right)$ or $P\left(v_i|\mathbf{h}\right)$ as appropriate. These can then be used to calculate estimates for $P\left(\mathbf{v}\right)$ and $P\left(\mathbf{h}\right)$ by marginalization over the conditional variable. In practice a full Gibbs sample every update iteration would be prohibitively slow and so we used an approximation called contrastive divergence,⁴⁶ in which the Gibbs sampler is initialized using the input data and a limited number of Gibbs steps are performed. In our implementation we use one contrastive divergence step (i.e.$CD\left(1\right)$)and so the data (or mini-batches of the data) is presented as a matrix and used to sample the hidden unit valueswhich are then used to update the values of the visible units. These values are used to update the network parameters using stochastic gradient descent (SGD).⁵¹

During trainingthe results of these updates are stored in three matrices that correspond to the weights as well as the network representation of the tumor data at the visible and hidden layers. These matrices correspond to the network reconstruction of the data (visible layer$\mathbf{V}$) the latent feature representation of the data (hidden layer$\mathbf{H}$)and the input-feature mapping (weights$\mathbf{W}$). When the network is trainedthese can be extracted and utilized in the analysis.

Modifications to the base RBM

We made a number of simple modifications to the base RBM described above to ensure the feature representation was interpretablegeneralizablestable and reproducible. These modifications are described below.

Amalgamation of feature matrices

Each individual network run provides a similarbut not identicalweight matrix. As suchweight matrices from each network run were amalgamated and filtered to form the final input-feature map. Numbers of featuresthe inputs they representtheir magnitude and order would not necessarily occur the same in each network and so we constructed an algorithm based on the cosine similarityis which depicted in Figure S10and outlined in Algorithm 1. Note that Low magnitude weights were those less than 50% of the maximum weight value for each hidden unit.

Synthetic data

To investigate whether our RBM network can identify true associations in data of multiple typeswe trained the network on a synthetic dataset with known associations and data generation methods. The values in the synthetic dataset were generated from function that encapsulated simple relationships when applied to binary latent variables; we also utilized various statistical distributions on top of these relationships to model types like proportions and counts that we might find in the real dataset. The synthetic data is constructed in seven ‘blocks’ to aid interpretation. In the first six blocks5 latent variables are mapped to 10 observed variables in exactly the same way. The difference between the blocks is the statistical distributions used to generate the values in the data. The final block consists of latent variables mapped to distributions from all the previous types. In total there were 72 ‘observed variables’ generated from 34 ‘latent variables’. To generate the data for a single synthetic samplethe 34 binary latent variables were sampled uniformly with probability of the latent variable being 1 set at 0.5. This was done for 200 synthetic samples to create the synthetic dataset.

Latent variable to observed variable mappings

We denote the ith simulated observed variables by $v_i$ and the jth latent variables as $l_j$. The first block consists of a simple binary mapping designed to determine if the network can extract the correct relationships with no sources of noise. These relationships are written as

From this we model several different types of unambiguous relationships as described on the rightincluding logical relationships AND $\left(v_2,v_3,v_4,v_5\right)$OR $\left(v_8\right)$ and NOT $\left(v_{10}\right)$.

We build the next five blocks on exactly the same relationships. Block 2 utilizes introduces noise into the mappingrequiring a uniformly sampled random value $\text{Unif}\left(\left[0,1\right]\right)$ to be greater than a threshold $t=0.25$ as well as the latent variable being equal to 1 for the relationship to be passed through to the observed variable(s). $\text{Unif}\left(\left[0,1\right]\right)>t$ returns 1 if the condition if fulfilled and 0 otherwise. Block 2 mappings can be written thus:

Note that where the latent variable maps to many observed variablesa random value is sampled separately for each observed value so they are correlated rather than identical.

Block 3 reflects binary to continuous value [0,1] relationships in which the observed value(s) are zero if the latent value is zero but take a uniformly distributed random value [0,1] if the latent variable is 1. The mappings therein can be written as:

Block 4 introduces sampling from a parametric probability distributionnamely the beta distribution $\text{Beta}\left(a,b\right)$to reflect proportions and other continuous values bounded by 0 and 1. Here we aim to model situations where there are generally lower or higher values depending on the status of the latent variable. Therefore the main difference between this block and previous ones is that the observed value is sampled from one of two differently parameterized beta distributions depending if the latent variable is 1 of 0. The mappings are:

Block 5 is similar to block 3where the observed value(s) are zero if the latent feature is zero. Howeverwhen the latent feature is 1 then the observed value is sampled from a Poisson distribution to simulate count data. All observed variables are rescaled by the maximum of ${\tilde{v}}_i=v_i/\mathit{max}\left({\mathbf{v}}_i\right)$but we leave this step out of the mapping equations below for simplicity and to aid comparison with the other blocks. We therefore write the mappings as:

In Block 6 we aim to model the situation where a perturbation to a cellular process elicits different levels of counts (such as gene expression in mRNA data for instance). The format is similar to block 4 but with Poisson distributions used instead of beta distributions. We adopt the same rescaling scheme as used in block 5.

Finally we include latent variables that are mapped to generating functions of more than one of the types in the blocks above. Simulated count data is again rescaled as before. The maps are given as

RBM results on synthetic data

We trained 2000 networks using 80% of the synthetic data as the training set (chosen uniformly at random). The remainder of the data was used as a validation set for early stopping using the procedure described above. To investigate if overfitting occurs and if our early stopping procedure could identify potential overfittingwe allowed training to continue if overfitting was suspected and monitored the behavior of the free energy. We found that overfitting was suspected in 123/2000 (6.15%) of the training runs and early stopping would have been invoked in these cases. We investigated what occurred by plotting the free energy of the training sets along with that of the validation set. We provide an example of a well-behaved profile (Figure S11A) and some examples of training runs that would have been stopped in Figures S11B–S11D. We observe that the free energy values oscillate with the periodicity of the cyclical learning rate described above and in the second half of training there are iterations where the free energy values decrease significantly across all setswhich corresponds to removal of dead hidden units in network pruning. In the samples where overfitting was not suspectedthe free energy of the validation set decreased at the same rate at the free energy values of the training setsremaining below the maximum value set by the training sets. This is exactly how we would expect the free energies to behave if there were no overfitting.

In training runs where overfitting was suspectedwe found that this was generally because the free energy of the validation set had increased to a slightly higher value than the maximum value of the training sets; this always occurred during the latter half of training when network pruning was taking place (as in Figures S11B and S11C). This could indicate the network is starting to overfit. Occasionally there would be runs where the free energy of the validation set was notably higher than the maximum of the training sets (e.g.Figure S11D)which is more concerning and could indicate a higher degree of overfitting. However in every case we noted that in the general trend of the free energy of the validation set was to decrease until training stopped at the maximum number of iterations; this actually indicates that the data is not being overfit as we would expect the free energy to increase if the model was losing generalizability by incorporating aspects only present in the training sets. Therefore it is inconclusive whether these runs are actually overfit. Nonethelessour early stopping procedure is very conservative and would have caught and removed these suspect runsleaving only indisputably non-overfit runsand so we are confident that overfit networks will not contribute to the final results.

We next performed training with early stopping enforced to investigate how well the network captures known relationships in the data. When training was completewe amalgamated the weight matrices using the procedure described in Figure S10 and the final input feature encoding is given in Figure S12.

We can use the direct binary mappings of Block 1 to investigate how the RBM attempts to encode the relationships provided in the hidden data. The algorithm unambiguously identifies the one-to-one mapping of feature 1 to input 1 of this blockas well as the one-to-many mapping from feature 2 to inputs 234 and 5. In features 3 and 4the algorithm can identify that input 6 and 8 arise from the same featureas do inputs 7 and 8but inputs 6 and 7 are not directly associated. The algorithm cannot identify the inverse relationship between inputs 9 and 10 and encodes them in two separate features (5 and 6). This is consistent with the way our approach is constructed as a positive weight matrix can only encode relationships in which a feature is active leading to inputs that are active rather than inactive feature giving rise to active inputs.

A similar pattern is seen in blocks 2–6where the algorithm is generally able to extract the correct (or logical) associations encapsulated in the features. There are two exceptions to thiswhich we have highlighted by annotations A and B. Annotation A shows a different encoding of the situation in which the original features both map to the same input (as in features 3 and 4 in block 1). Herethese are encoded as three inferred features where the first two encode a strong association with each input exclusive to each feature and a weak association with the shared input and the third feature encodes a strong association with the shared input with weak associations to both of the exclusive inputs. It is important to note that although this encoding does not precisely replicate the original input mappingthe network has still learned a logical way of encoding this relationship. Thereforewe do not consider this encoding incorrect. Annotation B shows inputs that are not assigned to any feature. Both of these are one-to-one mappings in block 6 (Poisson low/high) indicating that these relationships in this data type might be too subtle for the algorithm to distinguish. Note the one-to-one mapping of the 9th input in the block and the one-to-many mapping are identified correctly. In block 7the block consisting of one-to-many mappings of data of multiple typesthe algorithm identifies the correct number of features (4) and the correct inputs are assigned to each feature. These results provide empirical evidence that our RBM algorithm can extract relationships across a number of data types.

Consistency of results

There are several sources of variation between runs on the same dataset: the RBM is intrinsically a probabilistic networktraining it is a stochastic process and we are using different sets of data in each run. Howeveralthough variation between the features extracted in each run is inevitableit is important that they are consistent as an ensemble so we can be confident that the result is stable and the final amalgamated weights reflect the true relationships in the data.

To quantify the consistency between feature setswe use an approach based on the Hamming distance. If we describe a latent feature as a binary string that is equal to 1 with a non-zero weight is present and zero elsewherewe can then calculate the Hamming distance between individual featureswhich returns the number of input mappingsmthat are not shared between those features. We can use this to identify which feature in set ${\mathbf{F}}_j$ is equivalent to a given feature in set ${\mathbf{F}}_i$ (as that has the minimum Hamming distance) and then identify the feature in set $F_i$ that is most incongruous to their equivalent feature (as this has the maximum Hamming distance). We write the Hamming distance of the incongruous feature as ${\tilde{m}}_{i,j}$which can be expressed as

We can use this metric to assess how the input/feature mapping differs between runs. We randomly sampled $w=1000$ weight matrices (without replacement) from the 2000 matrices generated by networks trained on the synthetic data and amalgamated these as described above. We repeated this 100 times to give 100 feature setsand then we calculated ${\tilde{m}}_{i,j}$ for all pairwise combinations of i and j. Of the 10,000 resulting comparisonsthe greatest number of differences between equivalent features was 2which occurred 72 times (0.72% of the time). Figure S13 shows the histogram of the ${\tilde{m}}_{i,j}$ values for $w=1000$which reveals that the most common difference was 1which occurred 6001 times (60.01%). There was no difference in the feature maps in 39.27% of the runs. This is remarkably consistent given the aforementioned sources of variation.

Network training on real data

We used the real data to train various versions of the networks across a number of runs to obtain distinct outputs for use in the analysis. The goal of the first run was to extract the amalgamated weight matrix that describes the input to feature mapping (Figure S2). We used this to determine the latent feature representation used in the clustering (MP Figure 1) by training a new network run with the weight matrix initialized to the amalgamated weight matrix and setting the weight learning rate to zero. Learning of the biases was enabledas these may be different from the biases in the previous networks due to the removal of low magnitude weights. Once the remaining network parameters have converged during trainingtaking further iterations is equivalent to sampling the hidden units/feature representation for each patient. We therefore averaged the hidden unit values taken every 10 iterations during the final 1000 iterations to obtain the final feature representation.

The input-feature map was used to extract an informed subset of genetic alterations from the original inputs - these were used to determine associations in the ARBS analysis (MP Figure 2)as well as the inputs for the Ordering Analysis (MP Figure 3).

Two-stage clustering

The dimensionality of the feature representation is still quite large for conventional clustering techniques. Therefore we adopted a two-stage approach where we first clustered by those features that were most informative of clinical outcomecalculated the centroids of these first-stage clusters for all featuresand then clustered these in the second-stage of clustering to produce MP Figure 1. Here we provide more details on identification of informative features using a discrimination score and the clustering methods used.

Discrimination score

There have been several methods proposed for quantifying the relative importance of the units of a neural network.⁶² Howevermost of these are generally formulated to discover the inputs that are important in discerning the output.^63,64 In our applicationwe wish to quantify the discriminative capacity of each of the features (hidden layer) with respect to the clinical outcome. As we utilize non-negative weights to determine the relevance of the inputs to the hidden units in the feature extractionfor consistency we adopt a similar strategy to determine the relevance of the hidden units to adverse clinical outcome as determined by biochemical relapse.

To obtain the discrimination scores for each featurewe modified the architecture of the base RBM so that it was similar to ClassRBM.⁶⁵ This adds an extra classification layerwhich is fully connected to the hidden layerthe units of which contain the values of the classes. In ClassRBMthere is another set of weights that denote the strength of the connection between the hidden and classification layersand these are trained in the same bi-directional fashion as the input weights.

Howeverin our application we wish to uncover underlying relationships in the data (encapsulated by the features) in an unbiased way and then determine how relevant these features are to determining the clinical outcome. We therefore performed the learning of the latent feature representation and the discrimination scores separately to ensure that learning the classification weights to ensure that the latent representation remains unbiased by the knowledge of the clinical outcomeand the algorithm for feature learning described above can still be considered as unsupervised. This was done by fixing the input weights to the amalgamated weight matrix described above but then training the class weights using contrastive divergence as described above.

We also enforced a non-negative constraint on these class weightssimilar to the input weights. To get our discrimination scorewe take the absolute value of the weights corresponding to relapse minus the weights corresponding to no-relapse$\mathbf{s}$. This can be expressed mathematically as

were ${\mathbf{c}}_r$ are the class-weights associated with relapseand ${\mathbf{c}}_{r^{\prime }}$ are those associated with no relapse.

These $\mathbf{s}$ values can be considered as heuristic quantity relating importance of the corresponding feature to the clinical outputsimilar to how the component loadings quantify the explained variance of the corresponding principal component in principle component analysis (PCA). As there is no set rule for determining the number of featuresso we followed a similar approach to that conventionally used in PCA and selected the number of features using the cumulative distribution. We chose a cut off of 0.9 of the total cumulative discrimination scorewhich resulted in 14 out of 30 features being selected for the initial clustering phase.

Clustering

Clustering of tumors was performed on the latent feature representation in a two-stage process to facilitate the identification of clusters that were relevant to clinical outcome. As the feature representation for each patient can be considered as a vector containing the probabilities that the corresponding feature is activeit is appropriate to use a distance measure that quantifies the distance between probabilities. As suchwe calculated the mean Jensen-Shannon (J-S) divergence⁶⁶ between tumors in a pairwise fashion.

For a pair of patients A and Brepresented by the latent feature representation in hidden layers ${\mathbf{h}}_A$ and ${\mathbf{h}}_B$the mean J-S divergence can be written as

where $\mathbf{m}=\frac{1}{2}\left({\mathbf{h}}_A+{\mathbf{h}}_B\right)\text{,}$ is the midpoint of ${\mathbf{h}}_A$ and ${\mathbf{h}}_B$. The additive terms in the square brackets in Equation 73 represent the Kullback-Leibler divergence between each element of the latent feature representation for either patient and the corresponding element of the midpoint vector$\mathbf{m}$.

As we are not using a Euclidean distance metricclustering through k-means is not appropriate and so we used k-medoid clustering for the first stage; this is similar to k-means but selects a representative data point (medoid) as the centroid for each cluster instead of the mean. Using the silhouette method,⁶⁷ we determined that 11 clusters was optimal. For the second stage of clusteringwe used hierarchical clustering to cluster the medoids themselves (again using the J-S divergence)and this was used to generate and order clusters by the dendrogram MP Figure 1.

DNA breakpoint proximity to androgen receptor binding site

To examine the proximity of DNA breakpoints to androgen receptor binding sites (ARBS)we designed a permutation approach that quantifies the departure from a random distribution of the breakpoints across the genome. We downloaded processed ChIP-seq data targeting AR for 13 primary prostate cancer tumors from Gene Expression Omnibus (GEO: GSE70079)⁶⁸ and amalgamated them for use as the ARBS locations. For each of our 159 sampleswe simulated the scenario whereby breakpoints were randomly distributed across the genome. The simulation of breakpoints was performed chromosome-wise. For each chromosome we simulate 1000 sets of N breakpoint positionswhere N is the number of breakpoints we observed on that chromosome. These positions are randomly distributed (with a uniform distribution) across the full chromosome. Therefore the simulations intrinsically take into account the size of the chromosomes. We used the R package RegioneR⁶⁹ with genome assembly GRCh37masked for assembly gaps (AGAPS mask) and intra-contig ambiguities (AMB mask) to keep the possible chromosomal locations consistent with what could be observed in the real data.

To detect significant departure from a uniform random distributionwe calculated the proportion of breakpoints within 20,000 base pairs (bp) of an ARBS for the observed and permuted data ($B_{obs}$ and $B_{perm}$respectively). If $B_{obs}>p_{97.5\%}\left(B_{perm}\right)$the tumor was classified as Enrichedelse if $B_{obs}<p_{2.5\%}\left(B_{perm}\right)$the tumor was classified as Depleted. Otherwise the difference is not significant and the tumor was classified as Indeterminate. The level of enrichment or depletion of breakpoints in the proximity of ARBS used in MP Figure 2A was estimated according to the following formula:

To check the method was not inherently biasedwe performed the analysis on the breakpoints derived from the UK dataset (as in MP Figure 2) and compared these to the classes derived if the position of the AR-binding sites was distributed uniformly across the genome (Figure S14). As expectedwe find that almost all tumors in the randomized set were classed as Indeterminatein stark contrast to the real data.

The agglomerative hierarchical clustering of the ARBS groups across AustralianCanadian and UK datasets was generated using the R package pvclust⁷⁰ v2.0.0 using the ward.D2 clustering method with squared Euclidean distance (100,000 iterations). This package also enabled the estimation of the Approximately Unbiased Multiscale Bootstrap (AU) p values for the Depleted group. These clustering results were confirmed by a partitional clustering approach using the R packages cluster v2.1.0 and factoextra v1.0.5.

ARBS pairs required for DNA loop formation

In MP Figure 5 we investigated the role of AR in the formation of DNA loops that can precipitate DNA double-strand breaks (DSBs). For each ARBS that was previously identified as proximal to a DNA breakpoint in each patientwe determined the proportion that were also proximal to another ARBSas required by the mechanism for DNA loop formation described previously²¹ and depicted in MP Figure 5A. As it has been shown that the two ARBS involved in the DNA loop involved in the TMPRSS2/ERG fusion are separated by 19972 base pairs,²¹ we set the threshold for the second proximal ARBS as 30,000 base pairs in the direction away from the DNA breakpoint. p values were determined as before.

In all boxplots the red line denotes the medianthe blue box encapsulates the interquartile rangeand the black dashed lines denote the range of data not considered outliers; outliers (red dots) are as defined in the MATLAB boxplot() function default settings. The size of the angular ‘notch’ corresponds to a confidence interval around the mediansuch that if two notches are not overlapping then there is approximately 95% confidence that the median of the two groups differ. The folded notch in MP Figure 5D indicates that the notch extends past the interquartile range by the folded amount.

Ordering

We previously estimated consensus ordering of events by estimating phylogenetic trees from the cancer cell fraction (CCF) that contained each aberrationand applying the Bradley-Terry model to determine the most consistent order of events.¹¹ We recently released a study⁷¹ that improves this approach using a Plackett-Luce model,^72,73 which we also utilize in this study. We provide a complete description of the method for reproducibility.

There are a number of sources of uncertainty when attempting to determine the order of events from bulk DNA sequencing. In particularwe often cannot infer the true phylogenetic tree for each patientand furthermore it is impossible to determine the relative timing of events on parallel branches. Howeverwe can estimate the set of possible trees using the relative cancer cell fractions (CCFs) of the genomic aberrations involvedand from these we can estimate a set of possible orderings. Therefore we created an algorithm where we (sampled) a single possible tree from the dataand using this we sampled a viable order of events for each patient. This is repeated multiple times so that the uncertainty in these estimates is encapsulated in the output distributions. Algorithms of this type are called Monte Carlo simulations to emphasise the use of randomness in the procedure.

We adopted the Plackett-Luce model^72,73 to construct a probability distribution over the relative rankings of a finite set of itemsthe parameters of which can then be estimated from a number of individual rankings. This can be used to quantify the expected rank of each item relative to the others across the population. In our applicationan item corresponds to an eventnamely the emergence and fixation of a novel copy number alteration (CNA) identified in the extracted features. Ranking these events therefore relates to the order in which they would be expected to occur. We also utilized a Plackett-Luce mixture model,²⁰ which allow us to determine whether there are subpopulations in the data with different orderings.

The Plackett-Luce model

Given a set of CNA occurrences for each patient with associated subclonalitywe would like to infer the order which these events generally occur. To do this we used a Plackett-Luce modelwhich is formulated as a ranking methodand returns a value quantifying the ranking preference. We use a different interpretationnamely the orderingwhich is defined as the inverse of the ranking pref. ⁷⁴ Like the Bradley-Terry modelthe Plackett-Luce model does not return any temporal information outside the expected order of events.

We have a set of N copy number events we are interested in,

then we can apply Luce’s choice axiom,⁷² which states that the probability of selecting one event over another from a set of events is independent of the presence or absence of the other events in the set. We can therefore write the probability of observing event i as

where $\left\{{\alpha }_i\right\}$ are the coefficients that quantify the relative probability of observing the $i^{\text{th}}$ event. To reflect the ordering aspect of our application we refer to this value as the proclivity. Plackett⁷³ used this formalism to construct a generative model in which all N events are randomly sampled from C without replacement (i.e.a permutation). If we let $\mathrm{\Lambda }$ correspond to a permutation of the set C such that ${\lambda }_k\in C$ and ${\lambda }_1\prec {\lambda }_2\prec \dots \prec {\lambda }_N$then we write the probability density of a single ordering as

where ${\alpha }_{{\lambda }_k}$ is the proclivity associated with event ${\lambda }_k$and ${\mathrm{\Lambda }}^{\left(k\right)}=\left\{{\lambda }_k,{\lambda }_{k+1},\dots ,{\lambda }_N\right\}$ is the set of possible events after $k-1$ events have occurred.

Plackett-Luce mixtures

We hypothesized that there may be more than one set of copy number orderings present in our populationand so analysing all events in one ordering scheme may not be appropriate. Furthermorethe inhibition of AR-associated breakpoints implies that some CNAs may be found more frequently with a select set of otherswhich is in violation of Luce’s choice axiom. We therefore implemented a mixture modeling approach,^20,74 which reinstates Luce’s choice axiom as the selection of each CNA can be considered as independent conditional on the mixture component. Such a finite mixture model assumes that the population consists of a numberGof subpopulations. In this setting the probability of observing the ordering ${\mathrm{\Lambda }}_s$ for the $s^{\text{th}}$ sample is

where ${\omega }_g$ are the weight parameters (not to be confused with the weight matrix in the RBM) that quantify the probability that sample s belongs to subgroup g. The appropriate parameter values can be determined using maximum likelihood estimation via an EM algorithm.²⁰ The number of mixture components can be chosen using the Bayesian Information Criterion (BIC) estimationwhich is given by

where ${\mathrm{\Theta }}_{ML}$ is the parameter set that maximizes the log likelihood $\mathcal{l}(·)$N is the number of parametersand M is the number of samples.

Implementation

The general formulation of the Plackett-Luce model takes a matrix containing the sequence of events for each patient as its input. Howeverwe do not know the order in which these events occurredonly the presence and cancer cell fraction (CCF) of each CNA for each patient. As suchwe first estimate the phylogenetic trees for each patientand then determine the order of events from this. As we only have one tissue sample for each patientthere is often uncertainty in the tree topology and the possible sequence of eventsand so we use a Monte-Carlo sampling scheme in which we sample the trees and sequence of eventsand use these to estimate the distribution of possible orderings through the Plackett-Luce model. Samples with 0 or 1 CNA were not used in this analysis.

Another issue arises due to censoringwhich occurs when the sample is taken before all aberrations that would occur have occurredresulting in missing data. These are called partial-orderings in the Plackett-Luce frameworkand the general approach to addressing this is to reformulate the model so that all missing events are implicitly ranked lower than the observed data.^20,37 This may not be appropriate for our analysis as we may have multiple subgroupsand we anticipate that distinct aberrations may have similar or equivalent effects in each subtype and thus will rarely co-occur despite being indicative of the same type. For instancethe absence of a very early aberration may be due to the occurrence of another less frequent aberrationso including it at the bottom of the order would bias the rankings toward more frequent aberrations. As suchour algorithm works in two phases.

• (1)

  Determine the number of mixture components and assign patients to each component,

• (2)

  Estimate the ordering profiles of each component.

These are distinct as we treat the creation of the phylogenetic trees in a slightly different way in each of these processes to account for censoring. When estimating the number of componentswe calculate trees only using the observed CNAs. Howeverwhen estimating the full ordering profileswe introduce another sampling step into our Monte-Carlo scheme where we explicitly sample a number of additional CNAs with probability proportional to the subclonality of the aberration in tumors of each mixture component. Sampling in this way reduces the bias toward more frequent aberrations.

Assign samples to mixture components

In the first phasewe

• (1)

  Sample phylogenetic trees for each patient,

• (2)

  Sample sequence of events for each patient that are consistent with trees,

• (3)

  Calculate Bayesian Information Criterion (BIC) for 1–10 mixture components,

• (4)

  Repeat steps 1–3 1000 times,

• (5)

  Determine number of mixture components which consistently had lowest BIC score,

• (6)

  Assign patients to mixture components.

The phylogenetic trees are created by initially sorting the CNAs of each patient in descending order of CCF obtained from the output of the Battenberg algorithmiterating through them and sampling the possible parents with uniform probability. The CCF of a parent cannot be greater than the sum of the CCF of their childrenso viable parents are defined as ones where their CCF is greater than that of their current children plus the CCF of the CNA under consideration. The position in the sequence when the CNA occurred is sampled as any position after the parentwith uniform probability. The ordering estimates and assignment to the mixture components using the R package PLMIX as this incorporates mixture models and partial rankings (so the absence of a CNA from a sequence would not penalize its position in the ordering). A vector of assignments was retained for each sample runand the final assignment was determined by the most frequent assignment over the course of 1000 runs.

Estimate ordering profiles of each component

In the second phasewe

• (1)

  Sample phylogenetic trees for each patient,

• (2)

  Sample sequence of events for each patient that are consistent with trees,

• (3)

  Augment sequence with additional CNAs to alleviate censorship bias,

• (4)

  Calculate ordering profiles for each mixture component,

• (5)

  Repeat steps 1–4 1000 times,

• (6)

  Amalgamate results to determine final ordering profiles of each mixture component.

The phylogenetic trees and sequence of events were initially determined as before. Howeverinstead of utilizing partial rankings in the PL modelwe explicitly augmented the data with additional CNAs to account for those unobserved due to censorship. The probability of CNA being added to the sequence of events is equal to the proportion of subclonal occurrences relative to the total number of occurrences in the subpopulation defined by the mixture component. This can be written as

where $N_{sub}(·)$ and $N_{total}(·)$ denote the number of subclonal and total occurrences respectively of CNA $c_i$ in mixture component g. As events that are predominantly subclonal have a higher chance of being unobserved due to censorshipthis sampling scheme will mitigate this to a degree. Converselyevents that are predominantly clonal (i.e.early) may be unobserved due to factors other than censoringand these have a reduced chance of being imputed. Calculating these values using the patient samples for each mixture components rather than the entire population means that only CNA subclonality relevant to each subpopulation are considered. Imputation is performed by drawing a uniform random numberrfor each patient and including the CNA in the set of additional CNAs for each patient if $P\left({\tilde{c}}_{ig}\right)<r$. The set of additional CNAs for each patient are shuffled uniformly and added to the sequence. We then calculate the ordering for each mixture component individually using the Plackett-Luce model without partial ranking. This process is repeated 1000 times and the proclivity for each CNA is calculated and used to create an empirical distribution for proclivity for each CNAwhich are used to create the box-plots in MP Figure 3.

Synthetic data

We generated a synthetic dataset to evaluate our method. We simulated two subpopulationsA and B that each had exclusive sets of 5 ‘early’ and 5 ‘late’ CNAs as well as common sets of 5 early and 5 late CNAs. These can be written as

We used these sets to simulate the order of the events in each tumorand then created a set of CCF values consistent with this order. The main principle in the simulation is that early events will generally occur before late eventsbut there is no intrinsic order in the sets of early and late events themselves.

To obtain the set of events for each simulationYwe first sampled how many events occurred by the time of sampling from a Poisson distribution$e\sim Poiss\left(10\right)$ followed by the subpopulation S to draw fromwith $p\left(S=A\right)=p\left(S=B\right)=0.5$. To reflect the early/late ordering and exclusive/common nature of the CNAswe first sampled 5 events (or e events if $e<5$) from the pooled set of common and early events for that subpopulation$P_E=\left\{E_S\cup E_C\right\}$without replacementthat is $Y_E\subseteq P_E,\left|Y_E\right|=5$. If $e>5$we then pooled the events that had not been sampled from $P_E$ alreadydenoted here as $P_E^{\prime }$with the set of late events $P_L=\left\{L_S\cup L_C\right\}$ and sampled the remaining $e-5$ events from this set$Y_L\subseteq \left\{P_E^{\prime }\cup P_L\right\}$. We then randomly sampled the order in which the events in $Y_E$ occurredfollowed by those in $Y_L$. We then sampled how many of these events were clonal$e_c$uniformly at random as assigned these a CCF of 1. To obtain the CCFs values of the subclonal population we ranked the subclonal events$r_s\in \left\{1,2,\dots ,e_s\right\}$and assigned a CCF value as a linear function of their rank $CC{F}_k=1-r_s/\left(e_s+1\right)$.

We created synthetic CCF values for 200 tumors and used these as inputs into our algorithm described above. The BIC scores are shown in Figure S15Awhere two mixture components has the lowest BIC score and so the algorithm has identified the correct number of subpopulations. These subpopulations were used in to establish ordering profileswhich are shown in Figures S15B and S15C. We find that the algorithm has correctly identified all of the unique early and late events for each subpopulationas well as the common early and late events.

BIC scores for real data

Bayesian Information Criterion (BIC) scores were determined for each mixture component for each of the 1000 runs are shown in Figure S16. The BIC score was lowest for two mixture components for every sampled orderingand so this was taken as the value to use in subsequent analysis.

Statistical model of evotype convergence

We created a statistical model describing how the probability of convergence to the Canonical- or Alternative-evotypes changes as genetic alterations accumulate (MP Figure 4D). We assume that the accumulation of such aberrations in each individual tumor followed a stochastic process in which the order and relative timing of the aberrations occurred with some degree of randomness/stochasticity. Similar to the Ordering analysiswe utilized a statistical algorithm in which we simulated a number of possible aberrations consistent with the possible phylogenetic treesand then estimated the probability that tumors with these aberrations converged to the Canonical-evotype (the probability of convergence to the Alternative-evotype is 1 minus the probability of convergence to the Canonical-evotype). As many of the genetic alterations will occur clonally (i.e.with $CCF=1$) there is considerable uncertainty in the order in which they would have occurred. We incorporate this uncertainty into our approach using a method akin to propagation of errors through Monte Carlo simulation,⁷⁵ in which we sample (with uniform probability) from the space of possible trees for each tumor and calculate the probability of tumors with those genetic alterations converging to the Canonical-evotype. Repeating this random sampling many times provides a distribution of outputs that incorporates the uncertainty arising from our inputs in a principled fashionenabling downstream analysis.

Our algorithm is outlined in Figure S17. The algorithm iterates through an increasing number of aberrations (Loop i)performing several Monte-Carlo repeats of ordering samples (Loop j).

As individual evolutionary trajectories involve the stochastic accumulation of multiple genomic aberrationsis it impossible to specify each evolutionary route. Howeverwe can determine common modes of evolution by tracking the genetic alterations prevalent in tumors at the point of convergence to either evotype in our model. Through this we can identify paths in the probability density surface plot that correspond to the accumulation of these genetic alterations (black dashed linesMP Figure 4D)and the aberrations that distinguish them from each other.

Modeling the stochastic accumulation of genomic aberrations

We model the accumulation of aberrations in a tumor as a Poisson process.⁷⁶ For each iteration of Loop i we update the mean number of aberrations $x_i$this is then used as the input parameter to a Poisson random number generator to draw the number of aberrations to be samplednin each iteration of Loop j. We then identified those tumors with sufficient aberrations and selected one with uniform probabilityand used the data for these to sample a phylogenetic tree using the relative CCFs of the aberrations. We then used the phylogenetic trees to sample an order of occurrence for the aberrationsand retained the first n. The aberrations used were the SPOP mutations and the CNAs identified in the feature extraction; inter-intra chromosomal breakpointsETS status and chromothripsis are not included as these do not have associated CCFs and therefore cannot be used to determine the order of events.

We again use the data to calculate probability that tumors with the set of aberrations will be assigned to the Canonical-evotype. For a set of sampled aberrations$A_j=\left\{a_1,a_2,\dots ,a_n\right\}$we identified the patients for which $A_j\subseteq P_k$where $P_k$ denotes the full set of aberrations present in patient k. We can then identify which of these were assigned to the Canonical-evotype. We can now calculate the probabilities

where $N(·)$ denotes the number of tumors that obey the condition in brackets. We can now calculate the conditional probability

We performed 100,000 samples and thus obtained 100,000 values for each $p\left(Canonical|A_j\right)$. We input these values into a nonparametric density estimation scheme using Gaussian kernels with bandwidth 0.025. As we are estimating the probability density function of a set of probabilitieswhich are bound at $\left[0,1\right]$we ensured support only over this interval using the reflection method.⁷⁷ We performed this sampling step for $x_i\in \left\{0,0.01,0.02,\dots ,10\right\};i\in 1,2,\dots ,1000$.

Identifying genetic alterations in the convergent evolutionary trajectories

We used our model simulations to investigate the common evolutionary trajectories involved in convergence to each evotype (black dashed lines MP Figure 4D) as well as the aberrations that characterize them. In the modeling processwe recorded the order of genetic alterations for each of the trajectories used to calculate the pdf. We extracted each trajectory that had converged to the Canonical- or Alternative-evotypes (i.e.had a $p\left(Canonical|A_j\right)=0$ or 1) and assigned these into sets by the number of genetic alterations in the trajectories i.e.$\left\{A_1\right\},\left\{A_2\right\},\dots ,\left\{A_{10}\right\}$. We than ran a filtering step for each set where we removed any trajectories that had occurred in sets corresponding to fewer genetic alterationsmeaning we were left with trajectories that only converged to either evotype with the final genetic alteration for each set. We can then identify the position and frequency of occurrence of each genetic alteration in each set. The results of this are plotted in the bottom pane in Figures S7 and S8 for Canonical- and Alternative-evotype tumors respectively. Using this information we can calculate the pdf values for frequent combinations of genetic alterations in orderand use these to create the representative paths through the probability density (black dashed lines; MP Figure 4D).