Methods for cancer prognosis

Abstract

We describe a method and apparatus for stratifying subjects into one of two prognostic groups. The first prognostic group may be referred to as an alternative evotype group, and the second prognostic group may be referred to as a canonical evotype group. The canonical evotype group includes tumors that evolve along the same or different trajectories into a form of cancer that may be considered to have a standard form. The alternative evotype group includes tumors that evolve along the same or different trajectories into a non-standard or alternative form of cancer. One of the methods includes analyzing a biological sample obtained from a subject with cancer or metastatic disease using DNA sequencing and / or RNA sequencing, identifying genetic abnormalities in the biological sample, and classifying the subject into a first prognostic group based on the presence of one or more genetic abnormalities selected from Table 1 and a second prognostic group based on the presence of one or more genetic abnormalities selected from Table 2.

Description

[Technical field]

[0001] The present invention relates to a method for stratifying cancer patients into prognostic groups, in particular prostate cancer patient groups. [Background technology]

[0002] Tumor evolution is a dynamic, progressive process (1) involving the accumulation of genetic alterations that result in a pathological phenotype (2). In several cancer types, the genomic and expression alterations that characterize these phenotypes have been used to develop clinically-actionable classification frameworks (3-6). In localized prostate cancer, stratification based on the presence of specific molecular alterations (7), combinations of alterations (8), or gene expression profiles (9) have been proposed.

[0003] However, detailed investigations of the prostate cancer genome (10, 11, 12) have shown substantial heterogeneity in the occurrence of genomic variants that preclude the clinical utility of these simple classification schemes (13). More recent studies have distinguished between events that may occur early or late in the evolution of prostate cancer that may be informative for early-onset disease (14) and aggressive disease (15). However, the role of evolution in the development of disease types remains largely unclear.

[0004] Therefore, there is a need for methods that can classify cancer types in a clinically useful manner. [Prior art documents] [Non-patent literature]

[0005] [Non-Patent Document 1] Green and Sambrook et al., Molecular Cloning: A Laboratory Manual, 4th ed., Cold Spring Harbor Laboratory Press, Cold Spring Harbor, NY (2012). Summary of the Invention [Problem to be solved by the invention]

[0006] Cancer development is an evolutionary process, but the factors that drive the emergence of different disease types are still poorly understood. We applied three statistical and machine learning methods to genomic measurements from 159 prostate cancer patients, each of which identifies a different aspect of tumor evolution. Taken together, these results reveal that tumors follow an evolutionary trajectory that converges into two forms of the disease, defined as canonical and alternative evotypes. Canonical evotype tumors evolve into a standard form of the disease with a normal prognosis. Alternative evotype tumors evolve into a different form of the disease with a poor prognosis.

[0007] Statistical modeling revealed multiple routes to each evotype, dependent on the stochastic acquisition of complementary genetic alterations. Thus, evotype classification reflects the influence of several interacting factors and provides a powerful new paradigm for cancer stratification. [Means for solving the problem]

[0008] Therefore, according to the present invention there is provided a method as set out in the accompanying claims.Other features of the present disclosure will be apparent from the dependent claims and from the following description.

[0009] Stratification methods According to a first aspect of the present invention, there is provided a method for stratifying a subject into one of two prognostic groups, the method comprising: analyzing a biological sample obtained from a subject with cancer or metastatic disease using DNA and / or RNA sequencing; determining the location of a double-stranded DNA break relative to an androgen receptor binding site (ARBS); and classifying the cancer patient into a first prognostic group when the determined location is less frequently proximal to the androgen receptor binding site (ARBS) than expected; and classifying the cancer patient into a second prognostic group when the determined location is more frequently proximal to the androgen receptor binding site (ARBS) than expected. The method may further comprise classifying the cancer patient into the second prognostic group when there is no statistically significant difference between the closeness of the determined location to the ARBS and the expected closeness of the location to the ARBS.

[0010] The first prognostic group may be called the alternative evotype group, and the second prognostic group may be called the canonical evotype group. The canonical evotype group includes tumors that evolve along the same or different trajectories to a form of cancer that may be considered to have a standard morphology. The alternative evotype group includes tumors that evolve along the same or different trajectories to a non-standard or alternative form of cancer. Tumors that evolve into alternative forms have a worse prognosis than tumors that evolve into standard forms. An unfavorable prognosis is characterized by a lower likelihood of progression-free survival, determined through the time to biochemical recurrence, defined as a prostate-specific antigen (PSA) level in the blood above 0.2 ng / mL in the period after definitive treatment. Alternatively, a poor prognosis may be characterized by the observation of metastasis or death. Definitive treatment may include radical prostatectomy or radiation therapy. Tumors that evolve into standard forms have a standard prognosis, characterized by a higher likelihood of progression-free survival. For example, a subject may have a progression free survival of up to 120 months, or 100 months, or 80 months, or 60 months, or 40 months, or 24 months, or 12 months.

[0011] By using ARBS as explained above, this method of classification can be called ARBS classification. ARBS is the binding site for the androgen receptor protein. "Androgen receptor" (AR) is a DNA-binding transcription factor that regulates gene expression. AR is widely expressed in many cells and tissues, and thus AR has diverse biological actions, including important roles in the development and maintenance of the reproductive, musculoskeletal, cardiovascular, immune, nervous, and hematopoietic systems. AR signaling may also be involved in the development of tumors in the prostate, bladder, liver, kidney, and lung. AR has also been identified as playing an important role in prostate cancer, especially castration-resistant prostate cancer. AR is a member of the steroid hormone receptors and a group of steroid-induced transcription factors that share a consensus DNA-binding motif.

[0012] Methods for identifying transcription factor binding sites, such as androgen receptor binding sites (ARBS), are known to those skilled in the art, for example, chromatin immunoprecipitation assay combined with sequencing (ChIP-seq) is a common method for identifying genome-wide DNA binding sites for transcription factors. According to ChIP protocol, DNA binding proteins are immunoprecipitated using specific antibodies. The bound DNA is then co-precipitated, purified and salt sequenced. Sequencing can be performed using next-generation sequencing (NGS). In one example, ARBS can be identified using processed ChIP-seq data targeting AR for 13 primary prostate cancer tumors from Gene Expression Omnibus (Accession No. GSE70079, MM Pomerantz et al., Nature Genetics 47, 1346 (2015)). This ChIP-seq data can be merged to be used as the arrangement of ARBS.

[0013] The determined alignment closeness to the ARBS can be compared to the closeness of a baseline (or expected) distribution of breakpoint alignments to the ARBS to obtain an ARBS score that is used to classify patients into one of the first and second prognostic groups. The baseline distribution can be a random distribution of double-stranded DNA breakpoint alignments across the genome. The baseline distribution can be identified by a permutation approach. For example, the baseline distribution can be determined from a number of samples, such as those identified above in Nature Genetics. The observed breakpoints in the sample data across the genome (e.g., GRCh37) can be randomly shuffled 1000 times using the R package RegioneR (67) and masked for, for example, assembly gaps (AGAPS mask) and intra-contig ambiguity (AMB mask) to obtain a baseline distribution.

[0014] A double-stranded DNA break may be considered to be relatively proximal to the ARBS when the break is less than a threshold number of base pairs (e.g., 20,000 bps) from the ARBS, and may be considered to be distal when the break is greater than or equal to a threshold number of base pairs (e.g., 20,000 bps) from the ARBS. The method may include determining a percentage of the determined (e.g., observed) configurations that are less than a number of base pairs from the ARBS. The method may also include obtaining a percentage of configurations in a baseline distribution that are less than a number of base pairs from the ARBS. The method may include normalizing the determined percentage by the obtained percentage to obtain an ARBS score, and determining whether the determined configuration is more or less frequently proximal to the androgen receptor binding site than expected. When the percentage of the determined configurations that are relatively proximal is greater than an upper threshold (e.g., 97.5%) of the percentage of configurations in the baseline distribution that are relatively proximal, the determined configuration may be considered to be more frequently proximal than the given configuration (i.e., the tumor may be classified as enriched). When the ARBS score is used, the tumor may be classified as enriched when the ARBS score is higher than an upper threshold. When the proportion of the determined location that is relatively proximal is less than a lower threshold (e.g., 2.5%) of the proportion of locations in the baseline distribution that are relatively proximal, the determined location may be considered to be less frequently proximal than the given location (i.e., the tumor may be classified as depleted). When the ARBS score is used, the tumor may be classified as depleted when the ARBS score is lower than a lower threshold. If neither of these conditions are met, it may be considered that there is no statistically significant difference (i.e., the tumor may be classified as indeterminate).

[0015] Without wishing to be bound by theory, it is hypothesized that in alternative evotypes, some genetic alterations cause alterations in androgen receptor binding, promote DNA breaks distal to the ARBS, resulting in different copy number changes and giving rise to mechanistically distinct forms of cancer (e.g., prostate cancer) with poor prognosis. In contrast, in canonical evotypes, genetic alterations occur with double-stranded DNA breaks proximal to the androgen receptor binding site, resulting in standard progression of prostate cancer. In other words, the ARBS score (i.e., double-stranded DNA breaks) can be used for prostate cancer prognosis.

[0016] The ARBS score can be considered to represent genomic abnormalities in the sample. The inventors have identified additional genomic abnormalities characteristic of the first and second prognostic groups, namely the alternative and canonical evotypes of prostate cancer. Thus, the term genomic abnormality as used herein can be defined as any alteration of genomic sequence, such as deletion, insertion, inversion, duplication, loss of heterozygosity, gain of heterozygosity, DNA breakage, gene fusion, any other chromosomal mutation, or a measure of such alteration, such as PGA (percentage of genomic alteration), number of breakpoints, ARBS score. Genetic abnormalities that convey complementary information that can be used to distinguish prognostic groups can be identified using neural networks (e.g., restricted Boltzmann machines or autoencoders) and grouped together as features. Thus, each feature can represent multiple individual genetic abnormalities, and thus a complete set of abnormalities can be reconstructed with respect to these features and further analysis carried out directly on them.

[0017] The method according to any one of the preceding claims includes identifying further genomic abnormalities present in the sample and classifying the cancer patients into a first prognostic group based on the presence of one or more genomic abnormalities selected from Table 1 and into a second prognostic group based on the presence of one or more genomic abnormalities selected from Table 2. It will be understood that the absence of some or all of the genomic abnormalities in Table 1 may also indicate the second prognostic group. Similarly, the absence of some or all of the genomic abnormalities in Table 2 may also indicate the first prognostic group. This classification may thus combine the ARBS score with the presence of one or more genomic abnormalities to determine the most likely prognostic group. The probability that the classification is a correct classification may also be output together with the classification.

[0018] [Table 1]

[0019] [Table 2]

[0020] The number of breakpoints and / or ARBS scores are omitted in Tables 1 and 2 above, which focus on additional genomic abnormalities. It will be understood that combinations including two, three, four or more genomic abnormalities may be used in the classification step. Genomic abnormalities may be selected based on their importance to classification and / or the ease with which they can be identified in the sample. For example, genomic abnormalities in the target region may be included in such combinations in preference to genomic abnormalities that require whole genome testing.

[0021] The method detects loss of heterozygosity in one or more of the following regions: 2q14.3-2q23.3, 5q15-5q23, 5q11.1-5q14.1 (IL6ST, PDE4D), 6q12-6q22.32 (MAP3K7, ZNF292), 10q23.1-10q25, 16q12.1-16q24.3, 17p, 18q, 3q21.2-3q29, the entire chromosome 7, 8p23.3-8p22, 8q, 9q1 The method may further include identifying additional genomic abnormalities in the biological sample selected from the group including: gain of heterozygosity in one or more of the regions of 2.9-9q21.11 and entire chromosome 19, ratio of intrachromosomal to interchromosomal linkage variants, kataegis, ETS, percentage of genomic alterations (subclonal component), and percentage of genomic alterations (clonal component). Subjects were randomly assigned to receive genomic DNA samples containing loss of heterozygosity in one or more of the following regions: 2q14.3-2q23.3, 5q11.1-5q14.1 (IL6ST, PDE4D), 5q15-5q23, 6q12-6q22.32 (MAP3K7, ZNF292), 18q, loss of heterozygosity in one or more of the following regions: 3q21.2-3q29, entire chromosome 7, 8p23.3-8p22, 8q, 9q12.9-9q21.11. and / or a combination of genomic abnormalities selected from 2q14.3-2q23.3, 6q12-6q22.32 (MAP3K7, ZNF292), loss of heterozygosity in the region of 18q, and gain of heterozygosity in one or more of the entire chromosome 7 and regions of 8q, more specifically, in the first prognostic group.The subject may be classified into a second prognostic group based on the presence of one or more genomic abnormalities selected from the group including: loss of heterozygosity in one or more of the regions of 10q23.1-10q25, 16q12.1-16q24.3, 17p, gain of heterozygosity in one or more of the regions across chromosome 19, ratio of intrachromosomal to interchromosomal linkage variants, ETS, percentage of genomic alterations (subclonal component), and percentage of genomic alterations (clonal component), more specifically based on the presence of a combination of genomic abnormalities selected from the ratio of intrachromosomal to interchromosomal linkage variants, loss of heterozygosity in one or more of the regions of 10q23.1-10q25, 17p, ETS, and percentage of genomic alterations (subclonal component).

[0022] For ease of reference, these features are listed in the table below, ranked in order of the importance or significance of the features to the classification: Thus, the classification may be based on the most important combination of features, for example the top 5, or even the top 3.

[0023] [Table 3]

[0024] [Table 4]

[0025] In the above table, the abnormality may occur in any region of the chromosome referenced. In some cases, a specific gene is listed in brackets, which refers to a gene / genes that are present within the chromosomal region where the abnormality may occur, but the abnormality is not limited to occurring within this gene / genes.

[0026] Probabilities are returned with the identified genomic abnormalities, including ARBS scores, which indicate the probability that the tumor belongs to the first or second prognostic group (i.e., belongs to one of the evotypes) based on each identified genomic abnormality. It will be understood that the features listed for the first prognostic group are strongly negative for the second prognostic group, and vice versa. The classification can be based on consideration of all features, including features that, when present, strongly indicate a particular classification and features that, when present, strongly indicate not that particular classification. A probability threshold can be used to assign the classification, for example, when the probability is higher than p=0.5. It will be understood that the presence of individual genetic alterations will result in a smaller change in the probability of converging to a particular evotype than when a combination is present. Thus, probabilities are returned with the overall classification, and probabilities can be calculated based on the combination of genomic abnormalities that have been identified.

[0027] Previously collected samples with known clinical outcomes are analyzed and the samples may be classified into two or more clusters based on a combination of features indicative of a prognostic group. New samples may then be classified into prognostic groups based on the presence of features corresponding to the features of the clustered samples indicative of a prognostic group. This method of stratification may be considered to use a clustering classification (e.g., hierarchical clustering), where the classification may be selected from metacluster A or metacluster B. Classification as metacluster A may indicate a first prognostic group, i.e., alternative evotype. Classification as metacluster B may indicate a second prognostic group, i.e., canonical evotype. Metacluster B may be subdivided into two subclasses, i.e., metacluster B1 and metacluster B2.

[0028] Individual features characteristic of metacluster A (i.e., alternative evotype) may include combinations of genetic abnormalities including intrachromosomal structural variants (SVs), SPOP mutations, chromosoliposis, and loss of heterozygosity in regions 5q15-5q23.1 (spanning CHD1) and 6q14.1-6q22.32 (MAP3K7, ZNF292). A tumor may be classified as metacluster A when at least some of the abnormalities in this group (also referred to as clusters) are present. Features characteristic of metacluster B1 (i.e., canonical evotype) may include combinations of genetic abnormalities including ETS gene fusions and loss of heterozygosity in regions 17p (TP53) and 19p13.3-19p13.2 and 22q11.21-22q11.22. A tumor may be classified as metacluster B1 when at least some of these abnormalities are present. Features characteristic of metacluster B2 (i.e., the canonical evotype) may include a combination of genetic abnormalities including ETS gene fusions, interchromosomal linked structural variants (cSVs), and loss of heterozygosity in the regions of 5q11.1-5q14.1 (IL6ST, PDE4D), 10q23.1-10q25.1 (PTEN), and 17p (TP53). A tumor may be classified as metacluster B2 when at least some of these abnormalities are present.

[0029] For ease of reference, these features are listed in the table below, ranked in order of the importance of the features to the classification: Thus, classification may be based on a combination of the most important features, for example the top 5, or even the top 3.

[0030] [Table 5]

[0031] [Table 6]

[0032] This clustering classification can be used together with ARBS classification or separately from ARBS classification to stratify patients.Therefore, according to another aspect of the present invention, a method for stratifying cancer patients into one of two prognostic groups can be provided, which comprises: analyzing biological samples obtained from subjects suffering from cancer using DNA sequencing; identifying genomic abnormalities in the biological samples; and classifying cancer patients into a first prognostic group based on the presence of one or more genetic abnormalities selected from a set of genetic abnormalities comprising intrachromosomal structural variants, SPOP mutations, chromosoliposis, and loss of heterozygosity in the region 5q15-5q23.1 (spanning CHD1) and 6q14.1-6q22.32 (MAP3K7, ZNF292) by using clustering classification. Cancer patients may be classified into a second prognostic group based on the presence of one or more genetic abnormalities selected from a first set of genetic abnormalities including ETS gene fusions and loss of heterozygosity (LOH) in regions 17p (TP53) and 19p13.3-19p13.2 and 22q11.21-22q11.22 or a second set of genetic abnormalities including combinations of ETS fusions and interchromosomal linked structural variants (cSVs), as well as LOH affecting 17p (TP53), 10q23.1-10q25.1 (PTEN) and 5q11.1-5q14.1 (IL6ST, PDE4D).

[0033] The method may further comprise determining the order in which genomic abnormalities occur. This comprises performing bulk cell sequencing and determining the percentage of cells containing each genetic abnormality. It is determined that the abnormalities present in a higher percentage of cells occur before the abnormalities present in a lower percentage of cells. The ordering may be a consensus ordering, which may be determined using a statistical ranking method, such as the Plackett-Luce model. The cancer patients are classified based on the identified order, and such classification may be referred to as an ordered classification. Thus, according to another aspect of the present invention, a method is provided for stratifying cancer patients into one of two prognostic groups, the method comprising: providing a biological sample from a subject suffering from prostate cancer; analyzing the biological sample using bulk cell DNA sequencing; identifying genomic abnormalities present in the biological sample; determining the percentage of cells in which the genomic abnormalities exist; identifying the order in which the genomic abnormalities occur by determining that the genomic abnormalities present in a larger percentage of cells occur before the genomic abnormalities present in a smaller percentage of cells; and classifying the cancer patients into one of a first and a second prognostic group based on the identified order.

[0034] The ordering for the first prognostic group can be called ordering II, and the ordering for the second prognostic group can be called ordering I. Genomic abnormalities that can indicate ordering classification include some or all of the following: loss of heterozygosity, SPOP mutation, and ETS fusion in one or more of the following regions: 5q15-5q23.1 (spanning CHD1), 6q14.1-6q22.32 (MAP3K7, ZNF292), 8p (NKX3.1), 10q23.1-10q25.1 (PTEN), 13q12.3-13q21.1 (RB1, BRCA2) and 13q21.1-13q33.1 (EDNRB), 16q12.1-16q24.1 (CDH1) and 17p (TP53). Cancer patients are classified based on the identified order, and such classification can be called ordering classification. In other words, the accumulation of genetic abnormalities may provide further insight into the evolutionary trajectory of a cancer, for example, whether the cancer converges to an alternative evotype or a canonical evotype (i.e., ordering II or ordering I, respectively). Instead of (or in addition to) identifying the order in which genomic abnormalities occurred based on the percentage of cells in which genomic abnormalities are present, the method may include analyzing a second biological sample obtained from the cancer patient at a later time point using DNA sequencing, identifying genomic abnormalities present in the second sample, and comparing the genomic abnormalities identified in the second sample with the genomic abnormalities identified in the first sample to identify the order in which genomic abnormalities occurred. Additional subsequent samples may be obtained and analyzed at multiple time intervals, providing more detailed information about the order in which genomic abnormalities occur. As such, the method may be used in methods of monitoring disease progression and selecting a therapy for cancer treatment, or for therapy monitoring over time.

[0035] This ordering classification may be used together with the ARBS classification and / or the clustering classification, or separately from the ARBS classification and / or the clustering classification, to stratify patients.

[0036] Cancer patients can be classified as ordering II when loss of heterozygosity in one or more of the regions 6q14.1-6q22.32 (MAP3K7, ZNF292), 13q12.3-13q21.1 (RB1, BRCA2), and 13q21.1-13q33.1 (EDNRB) occurs early in the sequence of genomic abnormalities. Loss of heterozygosity in 5q15-5q23.1 (spanning CHD1) and SPOP mutations can also occur early in the sequence of genomic abnormalities for cancer to be classified as ordering II. High frequency of copy number gain can also indicate classification as ordering II. Late gain on chromosome 19 can also indicate classification as ordering II.

[0037] Cancer patients can be classified as ordering I when loss of heterozygosity or ETS fusion in 8p region (NKX3.1) occurs early in the sequence of genomic abnormalities. Loss of heterozygosity in one or more of the regions 10q23.1-10q25.1 (PTEN), 13q12.3-13q21.1 (RB1, BRCA2), 16q12.1-16q24.1 (CDH1), and 17p (TP53) that occurs later in the sequence of genomic abnormalities than loss of heterozygosity or ETS fusion in region 8p (NKX3.1) can also indicate ordering I. Late gains (i.e., only present in later samples or occurring in a small percentage of cells in a single sample) on chromosome 19 can also indicate classification as ordering I. Occasionally, very early loss of heterozygosity in region 1q42.12-42.3 can also indicate ordering I.

[0038] For ease of reference, the most relevant features for ordering classification are listed in the table below, ranked in order of the importance of the features to the classification. Thus, classification may be based on a combination of the most important features, for example the top 5, or even the top 3.

[0039] [Table 7]

[0040] [Table 8]

[0041] The three classifications may be applied separately or in combination to stratify the patient into one of two prognostic groups. When there is a combination of two or more classifications selected from the group of ARBS classification, clustering classification, and ordering classification, the selected classification may be considered as an intermediate classification. An overall classification may be determined based on the intermediate classification. For example, when all three intermediate classifications are used, an overall classification as a first prognostic group (alternative evotype) may be provided when at least two of the intermediate classifications classify the patient into the first prognostic group. In other words, the tumor has at least two intermediate classifications selected from the classification as metacluster MC-A, the ARBS classification of exhaustion, and the ordering classification of ordering II. The tumor may be assigned to a second prognostic group (canonical evotype) based on a similar majority approach, i.e., classifying the patient into the second prognostic group, when at least two of the classifications show the canonical evotype. For example, the clustering classification as metacluster MC-B1 or B2, the ARBS classification as enriched or indeterminate, and the ordering classification of ordination I indicate the overall classification as a canonical evotype.

[0042] Thus, according to another aspect of the present invention, there may be provided a method for stratifying cancer patients into one or two prognostic groups, the method comprising: analyzing a biological sample obtained from a subject having cancer or metastatic disease using bulk cell sequencing and DNA and / or RNA sequencing; determining in the biological sample the location of double stranded DNA breakpoints relative to the androgen receptor binding site (ARBS); identifying genomic abnormalities in the biological sample; and determining whether or not a genomic abnormality exists in the region 5q15-5q23. and determining the percentage of cells in the biological sample that have one or more genetic abnormalities selected from loss of heterozygosity, SPOP mutations, and ETS fusions in one or more of: 1 (spanning CHD1), 6q14.1-6q22.32 (MAP3K7, ZNF292), 8p (NKX3.1), 10q23.1-10q25.1 (PTEN), 13q12.3-13q21.1 (RB1, BRCA2) and 13q21.1-13q33.1 (EDNRB), 16q12.1-16q24.1 (CDH1) and 17p (TP53). The method includes obtaining a first intermediate classification using the ARBS classification, wherein the cancer patient is classified into a first prognostic group when the determined location is proximal to the androgen receptor binding site less frequently than expected, and into a second prognostic group when the determined location is proximal to the androgen receptor binding site more frequently than expected.The method includes obtaining a second intermediate classification using a clustering classification, wherein the cancer patients are classified into a first prognostic group based on the presence of one or more genetic abnormalities selected from a set of genetic abnormalities including intrachromosomal structural variants, SPOP mutations, chromosoliposis, and loss of heterozygosity in regions 5q15-5q23.1 (spanning CHD1) and 6q14.1-6q22.32 (MAP3K7, ZNF292), ETS gene fusions, and loss of heterozygosity in regions 17p (TP53) and and classifying into a second prognostic group based on the presence of one or more genetic abnormalities selected from a first set of genetic abnormalities including loss of heterozygosity at 19p13.3-19p13.2 and 22q11.21-22q11.22 or a second set of genetic abnormalities including ETS gene fusions, interchromosomal linked structural variants, and loss of heterozygosity in regions 5q11.1-5q14.1 (IL6ST, PDE4D), 10q23.1-10q25.1 (PTEN), and 17p (TP53). The method includes obtaining a third intermediate classification using an ordering classification, wherein the cancer patient is classified into one of the first and second prognostic groups based on the identified order in which the genomic abnormalities occurred, and identifying the order in which the genomic abnormalities occurred includes determining that the genomic abnormalities present in a large proportion of cells occurred before the genomic abnormalities present in a small proportion of cells. The method includes determining an overall classification as a first prognostic group when at least two of the first, second, and third intermediate classifications classify the patient into a first prognostic group, and as a second prognostic group when at least two of the first, second, and third intermediate classifications classify the patient into a second prognostic group.

[0043] The pathway to alternative evotypes may result from a sequence of steps in which several genetic alterations cause altered AR binding, which promotes DNA breaks at a different set of locations, which results in different copy number changes and gives rise to mechanistically distinct forms of disease. In contrast, canonical evotype tumors show genetic alterations that progress down the "default route" and have breakpoints near normal AR binding sites. These accumulate to the extent that the alternative route is closed (perhaps the cell becomes nonviable or progression on the pathway to the alternative route is now too far away). Thus, alterations to AR binding may be considered important in determining classification as canonical versus alternative evotypes. It will be appreciated that there are other effects of this AR cistrome modification that can also be used to examine evotypes (e.g., occurrence of point mutations in open chromatin regions associated with alternative AR binding).

[0044] The probability that the classification is the correct classification may also be output with the classification. For example, when an SPOP mutation occurs first, this gives a high probability (~0.91) of progression to the alternative evotype. As explained above, other routes leading to the alternative evotype involve the accumulation, in any order, of multiple separate LOH events involving genes such as MAP3K7, CHD1, or EDNRB. LOH of IL6ST or gain of the region 8p23.3-8p22 strongly influences the convergence to the alternative evotype after multiple abnormalities have already accumulated. Conversely, classification into a second prognostic group, for example the canonical evotype, may have a higher probability when several important abnormalities are identified, for example early TP53 loss or ETS gene fusions almost certainly ensure fixation to the canonical evotype. Loss of the region covering PTEN or CDH1 may also indicate the canonical evotype. For the canonical evotype, there were numerous aberrations that occurred late but ensured convergence, and thus in many cases the final step, notably LOH of 19p13.3-19p13.2, and gains of chromosome 19 and the region 22q11.1-22q11.23.

[0045] The above various classifications are based on the presence or absence of genomic abnormalities. Instead of considering each classification separately, subjects can be classified into one of the prognostic groups based on the identified genomic abnormalities. The probability that the classification is correct can also be output with the classification. Genomic abnormalities are ranked in order of importance according to the significance determined by the proportion of tumors that have this characteristic.

[0046] [Table 9]

[0047] It should be noted that the ARBS score does not have a positive association with the first prognostic group, however, a low ARBS score is highly indicative of the first prognostic group and may be used in combination with the presence of listed genomic abnormalities.

[0048] [Table 10]

[0049] Classification may be based on the combination of the most important features, for example the top 5 or even the top 3. Alternatively, the most important features may be selected by ranking the features based on their importance in subclassification. For example, features that are important to all three subclassifications may be preferably included when considering combinations. Features that are important to two subclassifications may be optionally included. Also, the time and effort involved in carrying out the test for the presence of a feature may affect whether the feature is included in the genomic abnormality to be identified before classification. Thus, features that reflect changes in specific regions of chromosomes may be selected for classification.

[0050] Genomic abnormalities can be identified using known methods. For example, structural variants can be detected using Brass (31). Somatic mutations in the SPOP gene can be determined using CaVEMan. Copy number alterations (i.e., LOH, HD, and gain) can be determined using the Battenberg algorithm on whole genome sequencing, or using ADTEx, CoNVEX SeqCNV on exome or targeted sequencing. ETS fusions can be detected using BRASS. Chromosomycosis can be identified by identifying copy number breakpoints and segmenting the inter-breakpoint distances along the genome using piecewise constant fitting (pcf from the R package copynumber v1.22.0). Regions with a density of greater than one breakpoint per 3 Mb can be flagged as being high-density regions. Chromosophthisis regions may be defined as regions with a high density of up to three allele-specific copy number states covering a proportion of the region greater than min(1,0.006N+1.1) with a number of copy number breakpoints N>15, a nonrandom segment size distribution (Kolmogorov-Smirnov test for exponential distribution, P<0.05), and the proportion of each type of structural variant is random with equal probability PTD=PDel=PH2Hi=PT2Ti=0.25 (multinomial test P>0.01), where TD=tandem duplication, Del=deletion, H2Hi=head-to-head inversion, and T2Ti=tail-to-tail inversion.

[0051] For example, kataegis can be identified using SeqKat https: / / github.com / cran / SeqKat. DNA breakpoints associated with linkage events can be identified using Chainfinder (http: / / archive.broadinstitute.org / cancer / cga / chainfinder) version 1.01. Clonal / subclonal ratio can be used to quantify the number of SNVs that are in all cancer cells (clonal) or only a subset (subclonal) in a sample, i.e., SNVs with cancer cell fraction (CCF)=1 and CCF<1, respectively. Genomic alteration percentage can be calculated as the total percent of the genome affected by CNAs (copy number alterations) (37). We also recorded the percentage of affected clonal and subclonal CNAs (i.e., CNAs with CCF=1 and CCF<1, respectively). "Number of breakpoints" is the number of DNA breakpoints that can be determined by BRASS. Inter- / intra-chromosomal breakpoints are determined using Chainfinder, which can identify breakpoints that occurred as part of the same event (e.g., chromosomes co-locating in a transcription factory, splitting, being misplaced, etc.). If these events involve different, i.e., non-homologous, chromosomes (e.g., inter-chromosomal translocations), they are classified as inter-chromosomal breakpoints, and if they only involve the same chromosome (e.g., deletions), they are classified as intra-chromosomal breakpoints. The ratio between the number of breakpoints in the two categories can then be determined. "Gene fusion" refers to a hybrid gene formed from two previously separate genes. Multiple genetic events, such as gene translocations, deletions, etc., can cause gene fusions.

[0052] The term "loss of heterozygosity" may refer to a chromosomal event in which a gene or chromosomal region is lost, a common form of allelic imbalance in which one of the two alleles is lost, causing a heterozygous cell to become homozygous. Loss of heterozygosity can cause somatic loss of the wild-type allele, and this form of chromosomal instability is sufficient to provide a selective growth advantage and has been recognized as a major cause of tumorigenesis. It will be understood that the term "gain" is used to indicate that a region of a chromosome or an entire chromosomal arm has been duplicated once or multiple times.

[0053] Gene loss can occur through many mechanisms, including large deletions, which are often the result of nonsense mutations or frameshifts. The former is the result of a standard mutation from one nucleotide to another, causing a premature stop codon. Frameshift mutations are the result of insertions or deletions, usually small, not a multiple of three, within the coding region that change the way the codon is translated into amino acids. These processes can occur, resulting in the loss of protein-coding genes. However, similar effects can be achieved for the loss of non-coding genes or regulatory regions.

[0054] The term "subject" or "patient" refers to an animal that is the object of treatment, observation, or experiment. By way of example only, a subject includes, but is not limited to, a mammal, which includes, but is not limited to, a human or non-human primate, non-human mammal such as a murine, bovine, equine, canine, ovine, feline, etc.

[0055] DNA and / or RNA sequencing is used to analyze the sample. Next generation sequencing, high throughput sequencing methods can be used. Whole genome sequencing, whole exome sequencing, targeted gene sequencing, RNA-seq, transcriptome sequencing, methylation sequencing, bisulfite sequencing, or combinations thereof can be performed. Bulk cell sequencing and / or single cell sequencing methods can be used. The methods used for DNA and RNA extraction are known in the art and can be used to obtain DNA or RNA from various samples to perform sequencing.

[0056] The cancer may be prostate cancer.

[0057] The biological sample may be a biopsy or a whole blood specimen from a tumor.The biological sample may be obtained during a radical prostatectomy or biopsy or transurethral electroresection performed on the subject.The sample may be fresh frozen or formalin fixed.

[0058] Samples may be obtained from subjects suffering from acinar adenocarcinoma, ductal adenocarcinoma, transitional cell carcinoma (urothelial carcinoma), squamous cell prostate carcinoma, small cell prostate carcinoma, large cell prostate carcinoma, mucinous adenocarcinoma, signet cell prostate carcinoma, basal cell prostate carcinoma, leiomyosarcoma, or rhabdomyosarcoma.

[0059] Complementary method As mentioned above, the first prognostic group may be called alternative evotypes and the second prognostic group may be called canonical evotypes. The alternative evotypes are associated with a poorer prognosis than the canonical evotypes. Thus, the method of stratifying cancer patients into one of the two prognostic groups described above may also be used to predict the clinical outcome of the patient.

[0060] The above method may include the additional step of treating the patient with a therapy. Alternatively, as another aspect of the present invention, a method for treating cancer in a subject may be provided, comprising stratifying the subject into one of two prognostic groups according to the method described above, and further comprising administering a cancer therapy to the subject.

[0061] There are many therapies recommended for the treatment of prostate cancer. Treatments are recommended depending on the stage of progression of the disease. Radiation therapy, hormone therapy, and chemotherapy are the three options often used in the treatment of prostate cancer. A single treatment or a combination of treatments may be used.

[0062] Chemotherapy is often used to treat prostate cancer that has spread to other organs of the body (metastatic prostate cancer). Chemotherapy destroys cancer cells by interfering with their growth. Chemotherapy may be used to control prostate cancer and reduce symptoms, thus minimizing the impact on daily life. Patients may receive hormone therapy before undergoing chemotherapy to increase the chances of successful treatment.

[0063] Radiation therapy can be used to treat localized and locally advanced prostate cancer. Radiation therapy can also be used to slow the progression and relieve symptoms of metastatic prostate cancer. Hormonal therapy may also be recommended after radiation therapy to reduce the chance of recurrence.

[0064] Hormone therapy is often used in combination with radiation therapy. Hormone therapy should not usually be used alone to treat localized prostate cancer in men who are fit and willing to undergo surgery or radiation therapy. Hormone therapy can also be used to slow the progression and relieve symptoms of advanced prostate cancer. Hormones control the growth of cells in the prostate. In particular, prostate cancer requires the hormone testosterone to grow. The goal of hormone therapy is to block the action of testosterone, either by stopping the production of testosterone or by stopping the patient's body from using testosterone.

[0065] Other therapies that may be used to treat prostate cancer include surgery, such as radical prostatectomy, high intensity focused ultrasound, cryotherapy, brachytherapy, patient monitoring, transurethral resection of the prostate, gene therapy, viral therapy, RNA therapy, bone marrow transplantation, nanotherapy, targeted anti-cancer therapy, or oncolytic drugs. Examples of therapeutic agents include steroids, antibodies targeting prostate specific membrane antigen (PSMA), checkpoint inhibitors, anti-tumor drugs, immunogens, attenuated cancer cells, tumor antigens, tumor-derived antigens, or antigen-presenting cells such as dendritic cells pulsed with nucleic acids, immune stimulating cytokines (e.g., IL-2, IFNa2, GM-CSF), targeted small molecules and biological molecules (such as agents that bind to tumor-specific antigens, including components of signal transduction pathways, such as modulators of tyrosine kinases and inhibitors of receptor tyrosine kinases, EGFR antagonists), anti-inflammatory agents, cytotoxic agents, radiotoxic agents, or immunosuppressants and cells transfected with genes encoding immune stimulating cytokines or steroids (e.g., GM-CSF).

[0066] The method allows clinicians to select appropriate treatment for prostate cancer based on evotype (i.e., based on prognostic group). For example, alternative evotypes are associated with poor prognosis, and therefore more aggressive treatments can be selected. Aggressive treatments can be selected from one or more of external beam radiation, brachytherapy, radical prostatectomy, hormone therapy, and / or chemotherapy. In particular, external beam radiation and brachytherapy can be used in combination. Canonical evotypes are associated with standard prognosis, and therefore standard treatments involving patient monitoring, radiation therapy, hormone therapy, chemotherapy, and / or combinations thereof can be adopted.

[0067] The progression of prostate cancer can be determined using the methods of the invention. The therapeutic efficacy of prostate cancer treatment can be monitored using the methods of the invention.

[0068] kit Kits for use in the methods of the invention are provided herein.

[0069] The kit may include reagents for performing DNA and / or RNA sequencing, and tools for detecting DNA double-strand breaks.The reagents for DNA and / or RNA sequencing may be for whole genome sequencing, whole exome sequencing, targeted gene sequencing, RNA-seq, transcriptome sequencing, methylation sequencing, bisulfite sequencing, or combinations thereof.Double-strand breaks may be identified by probes such as fluorescent nucleotides or antibodies.Double-strand DNA breaks may also be identified in sequencing data using computational tools such as BRASS.

[0070] The kit may also include instructions for use.

[0071] The kit may also include probes for the detection of one or more of the genomic abnormalities listed in Table 1, Table 2, or any one of Tables A to F. These probes may be DNA probes for annealing to specific target sequences in the sample DNA, for example FISH (fluorescence in situ hybridization) may be used.

[0072] The kit may include one or more probes for detecting one or more genomic abnormalities selected from loss of heterozygosity in regions 17p (TP53), 19p13.3-19p13.2, and 21q22.2-21q22.3 (ERG), and instructions for use. Additionally or alternatively, the kit may comprise one or more probes for detection of one or more of genomic abnormalities selected from loss of heterozygosity in any one of the regions 1q42.12-1q42.13, 2q14.3-2q23.3, 5q11.1-5q14.1 (IL6ST, PDE4D), 5q15-5q23.1 (spanning CHD1), 6q14.1-6q22.32 (MAP3K7, ZNF292), 13q12.3-13q21.1 (RB1, BRCA2), 13q21.1-13q33.1 (EDNRB), and / or gain in any one of the regions 3q21.2-3q29, chromosome 7, 8p23.3-8p22, and 8q (MYC).

[0073] The kit detects the following regions: 1q42.12-1q42.13, 2q14.3-2q23.3, 5q11.1-5q23.1 (IL6ST, PDE4D), 5q15-5q23.1 (CHD1), 6q12-6q22.32 (MAP3K7, ZNF292), 13q12.3-13q21.1 (BRCA2, RB1), 13q13.3-13q33.1 (EDNRB), 17p (TP53), 19p13.3-19p13.2 (LKB), 21q22.2-21q22.3 ( The present invention may further comprise probes for detection of one or more of the following genetic abnormalities selected from: loss of heterozygosity in chromosome 3q21.2-3q29, chromosome 7, 8p23.3-8p22, gain of heterozygosity in 8q (MYC), intrachromosomal structural variants, SPOP mutations, kataegis, chromosomal thromboplastin, ETS gene fusions, PGA clonals, high number of breakpoints, high ratio of interchromosomal / intracromosomal breakpoints. In one embodiment, the kit may further comprise probes for detection of one or more of genetic abnormalities selected from intrachromosomal structural variants, SPOP mutations, chromosoliposis, ETS gene fusions, loss of heterozygosity in the regions 5q15-5q23.1 (spanning CHD1), 6q14.1-6q22.32 (MAP3K7, ZNF292), 17p (TP53), 19p13.3-19p13.2, 22q11.21-22q11.22, 5q11.1-5q14.1 (IL6ST, PDE4D), and / or 10q23.1-10q25.1 (PTEN), and instructions for use. In one embodiment, the kit comprises a probe for detection of one or more genetic abnormalities selected from 5q15-5q23.1 (spanning CHD1), 6q14.1-6q22.32 (MAP3K7, ZNF292), 8p (NKX3.1), 10q23.1-10q25.1 (PTEN), 13q12.3-13q21.1 (RB1, BRCA2), 13q21.1-13q33.1 (EDNRB), 16q12.1-16q24.1 (CDH1) and 17p (TP53), SPOP mutations, ETS fusions, and instructions for use. [Brief description of the drawings]

[0074] [Figure 1a] 1 is a schematic flow chart of computer-implemented steps for discovering canonical or alternative evotypes. [Figure 1b-1] FIG. 2 illustrates different types of data used as input data. [Figure 1b-2] FIG. 2 illustrates different types of data used as input data. [Figure 1b-3] FIG. 2 illustrates different types of data used as input data. [Figure 2a] FIG. 1b is a schematic diagram of a neural network engine implementing the latent feature model used in the method of FIG. 1a. [Figure 2b] FIG. 1 is a schematic diagram of a latent feature model. [Diagram 3] 13 is a graph showing the frequency with which a particular number of features are estimated across 200 network runs using a subset of the data. [Figure 4] FIG. 2b is a schematic diagram of the steps taken when combining multiple weighting matrices to arrive at the fixed weighting matrix of FIG. 2a; [Diagram 5] 5a, b are heatmaps showing the relationship between the input data and patients, and between the reduced set of feature data and patients. [Figure 6] FIG. 1 is a dendrogram showing the probability of observing the listed feature in each cluster along with a discrimination score that quantifies the relevance of each feature in predicting recurrence. [Figure 7a] FIG. 13 plots the normalised proportion of DNA breakpoints for each sample ordered by normalised proportion. [Figure 7b-1] 7b is a heatmap of genomic features for each sample using the ordering from FIG. 7a. [Figure 7b-2] 7b is a heatmap of genomic features for each sample using the ordering from FIG. 7a. [Figure 7c]FIG. 1 is a dendrogram showing the proportion of CNAs in each ARBS group. [Figure 8] A plot of the BIC score against the number of mixture components showing the mean score as well as the individual scores. [Figure 9a] FIG. 13 is a plot of the Proportion of samples versus the Plackett-Luce coefficient for Ordering-I and Ordering-II. [Figure 9b] FIG. 1 shows a plot of copy number alteration versus Plackett-Luce coefficient for Ordering-I and Ordering-II. [Figure 10a] FIG. 1 plots progression free survival versus time for patients with tumors classified as either evotype. [Figure 10b] This figure plots the proportion of each evotype in each tumor stage. [Figure 10c] The figure plots the percentage of each ISUP Gleason grade group and each evotype. [Figure 10d] The figure shows the percentage of each evotype plotted according to PSA (ng / ml). [Figure 10e] 1 is a bar graph showing the prevalence of each gene abnormality in each evotype. [Figure 11a] 1 is a flowchart of a statistical algorithm for obtaining the probability of convergence to a canonical or alternative evotype based on the accumulation of genetic modifications. [Figure 11b] 1 is a surface plot showing the probability density that a tumor is assigned to a canonical evotype versus the number of aberrations. [Figure 12a] FIG. 2D surface plot showing the probability density that all canonical evotype tumors are assigned to the canonical evotype with increasing numbers of aberrations. [Figure 12b] Graph showing the percentage of lineages that converged to the canonical evotype at each number of genetic modifications. [Figure 12c-1] 1 is a bar graph showing the relative proportion of genetic alterations and the positions where they occurred for lineages that converged to the canonical evotype. [Figure 12c-2] 1 is a bar graph showing the relative proportion of genetic alterations and the positions where they occurred for lineages that converged to the canonical evotype. [Figure 13a] FIG. 2D surface plot showing the probability density that all alternative evotype tumors are assigned to alternative evotypes with increasing number of aberrations. [Figure 13b] 1 is a graph showing the percentage of lineages that converged to alternative evotypes at each number of genetic modifications. [Figure 13c-1] 1 is a bar graph showing the relative proportion of genetic alterations and the positions where they occurred for lineages that converged to alternative evotypes. [Figure 13c-2] 1 is a bar graph showing the relative proportion of genetic alterations and the positions where they occurred for lineages that converged to alternative evotypes. [Figure 14a] FIG. 1 shows the abnormalities present in canonical evotype tumors when divided into ETS− and ETS+. [Figure 14b] Kaplan-Meier plot of ETS+ and ETS- tumors classified into canonical evotypes. [Figure 15a] 1 is a schematic flow chart of computer-implemented steps for classifying tumors as canonical or alternative evotypes. [Figure 15b]13 is a plot of the relevance of features for classifying tumors into a particular cluster, metacluster A. [Figure 15c] Plot of the relevance of features for classifying tumors into a particular cluster, metacluster B. [Figure 15d] Plot of feature relevance for classifying tumors into alternative evotype (depleted tumors) or canonical evotype (enriched tumors). [Figure 15e] Plot of feature relevance for classifying tumors into alternative evotype (depleted tumors) or canonical evotype (enriched tumors). [Figure 15f] 1 is a plot of the relevance of features for classifying tumors into a particular ordering, ordering II. [Figure 15g] 1 is a plot of the relevance of features for classifying tumors into a particular ordering, ordering I. [Figure 15h] 13 is a plot of feature relevance for directly classifying tumors as canonical or alternative evotype tumors. [Figure 16-1] 13 is a plot of feature relevance for directly classifying tumors as canonical or alternative evotype tumors using RNA sequencing. [Figure 16-2] 13 is a plot of feature relevance for directly classifying tumors as canonical or alternative evotype tumors using RNA sequencing. [Figure 17] FIG. 1 is a schematic diagram of an associated system for performing computer-implemented aspects of these methods. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0075] The present invention is now further described. In the following text, different aspects of the present invention are defined in more detail. Each aspect thus defined may be combined with any one or more other aspects, unless otherwise specified. In particular, any feature shown to be preferred or advantageous may be combined with any other feature or features shown to be preferred or advantageous. The implementation of the present invention employs conventional techniques of immunology, molecular biology, cell biology, chemistry, biochemistry and recombinant DNA technology, which are within the skill of the art, unless otherwise specified. Such techniques are fully described in the literature, see, for example, Green and Sambrook et al., "Molecular Cloning: A Laboratory Manual, 4th ed.", Cold Spring Harbor Laboratory Press, Cold Spring Harbor, NY (2012).

[0076] FIG. 1a is a schematic flow chart of steps in a method for discovering canonical and alternative evotypes, which may be a computer-implemented method. As shown, the first step is to receive a dataset of information collected from a patient's tumor sample (step S100). The dataset may be collected by using DNA or RNA sequencing, for example, in this case, whole genome sequencing (target depth: 50X) on the sample along with a matched blood control. The dataset includes a large number (e.g., more than 120, perhaps as many as 140) of summary measurements or genomic features, and each sample may be represented in terms of these features (step S102). These measurements may include some or all of the following: number of single nucleotide variants (SNVs), number of indels (insertions, deletions, or complexes), number of structural variants including genomic rearrangements (inversions, deletions, tandem duplications, and translocations) and mutational signatures, percentage genome alterations (PGAs), DNA breakpoints (and whether they involve linked intrachromosomal or interchromosomal structural variants), telomere length, number of gene fusions, whole genome duplications (WGDs), presence or absence of any one of kataegis, ETS+ status, and chromosomal thromboplasty, presence or absence of significant driver mutations, copy number alterations (divided into CNA, loss of heterozygosity (LOH), homozygous deletions (HD), and gains). Data may be collected using any known technique, including those described below with respect to the data used to develop the classification.

[0077] The next step S104 is where the input data is reformulated as a reduced set of features that encapsulate the underlying relationships between the original inputs. As explained in detail below, this may involve using an adapted unsupervised neural network to perform feature learning on the dataset to identify associations between the inputs to obtain a reduced set with, for example, 30 features (shown in the table below).

[0078] [Table 11A]

[0079] [Table 11B]

[0080] In this approach, the relationship between the inputs and features is more easily interpretable and named with the corresponding genomic abnormality. Thus, the term genomic abnormality as used herein may be defined as any alteration of the genomic sequence (i.e., genetic abnormality), such as a deletion, insertion, inversion, duplication, loss of heterozygosity, DNA breakage, gene fusion, any other chromosomal mutation, or a measure of such genomic alteration, such as PGA (percentage of genomic alteration), number of breakpoints, ARBS score. When a feature reflects attributes of multiple genomic inputs, the attributes are separated by semicolons in the feature name. We expressed each sample in terms of these genomic features, which formed the basis of our findings described below. As set out in step S106, the next step is to quantify the discriminatory ability of each feature in predicting disease recurrence to identify patterns of genomic abnormalities indicative of adverse clinical outcomes.

[0081] Using the information from step S106 together with the feature presentation as input to the two-stage clustering method resulted in the identification of two different metaclusters characterized by different sets of abnormalities. Thus, as set out in step S108 and described in more detail below, tumors can also be classified as belonging to metacluster A (MC-A), metacluster B1 (MC-B1) or metacluster B2 (MC-B2). Tumor samples showing a combination of intrachromosomal structural variants (SVs), SPOP mutations, chromosomal schizophrenia, and loss of heterozygosity (LOH) in regions 5q15-5q23.1 (spanning CHD1) and 6q14.1-6q22.32 (MAP3K7, ZNF292) can be classified as metacluster A (MC-A). Tumor samples that exhibit a combination of ETS fusions and LOH affecting 17p (TP53) and regions 19p13.3-13.2 and 22q11.21-22q11.22 are classified as metacluster B1 (MC-B1).Tumor samples that exhibit a combination of ETS fusions and interchromosomal linked structural variants (cSVs) and LOH affecting 17p (TP53), 10q23.1-10q25.1 (PTEN), and 5q11.1-5q14.1 (IL6ST, PDE4D) are classified as metacluster B2 (MC-B2).

[0082] The next step is to examine the influence of the androgen receptor (AR) on the involved DNA breakpoints. AR is known to precipitate DNA double-strand breaks (DSBs) in conjunction with topoisomerase II-beta, and AR-associated breakpoints are frequently found in early-onset prostate cancer. As indicated in step S110, tumors can be classified as enriched when breakpoints occur significantly more frequently than expected proximal to AR-binding sites (ARBS), as depleted when breakpoints occur significantly less frequently than expected proximal to AR-binding sites (ARBS), or as indeterminate if they do not show a statistically significant association. By examining the ARBS group in conjunction with previously identified features, as described in detail below, depleted tumors were associated with multiple CNAs, chromosomal thrombosis, and SPOP mutations. Enriched tumors were associated with CNAs affecting 16q12.1-16q24.3 (CDH1) and 17p (TP53), high interchromosomal / intrachromosomal cSV ratios, and ETS fusions. Further clustering work also confirms the association of these CNAs with ARBS distal breakpoint prevalence.

[0083] The next step was to adapt a Plackett-Luce mixture model to extract a consensus ordering of the CNAs identified in the genomic features. Bayesian model selection determined that two separate ordering profiles were optimal. For ordering classification, each tumor was classified as belonging to one of two orderings -- ordering-I and ordering-II -- (step S112). The two profiles showed notable differences. Tumors classified as ordering-I frequently underwent early 8p LOH (spanning NKX3.1) and ETS fusions, and lack of LOH of regions covering RB1, BRCA2, CDH1, TP53, or PTEN genes could also occur. Very early LOH of 1q42.12-42.3 was also possible for tumors in this ordering. Tumors classified as ordering-II often showed early LOH events covering MAP3K7 and 13q (EDNRB, RB1, BRCA2), as well as copy number gains. Early mutations in the SPOP gene and LOH covering CHD1 are also possible but much rarer. Both sequencings showed late gains of chromosome 19.

[0084] The agreement of these three classification methods revealed a remarkable relationship. We introduce the term evotype to describe tumors linked by a common evolutionary mode resulting in similar disease characteristics. The metacluster MC-A is, to a large extent, a subset of the depletion group, and both are almost entirely subsets of ordering-II. We can therefore infer that there exists a subset of tumors that exhibit all the corresponding characteristics, i.e., evolutionary trajectories (ordering-II), breakpoint mechanisms (ARBS classification of depletion), and characteristic patterns of abnormalities (metacluster: MC-A). The term evotype may be used to describe tumors linked by a common evolutionary mode resulting in similar disease characteristics. Tumors that are assigned to at least two of MC-A, depletion, or ordering-II may be classified as alternative evotypes. Similarly, a tumor that is classified into at least two of the clustering classifications as metacluster MC-B1 or B2, the ARBS classification of enriched or indeterminate, and the ordering classification of Ordering-I indicates an overall classification as a canonical evotype. This match can be used to assign the tumor to one of the two evotypes (step S114).

[0085] Patient samples and data used in developing the classification The clustering, ARBS, and ordination classifications used above are based on applying three statistical and machine learning methods to genomic measurements collected from 159 samples. Data were collected from cancer samples from 205 patients treated at the Royal Marsden NHS Foundation Trust in London, Addenbrooke's Hospital in Cambridge, Oxford University Hospitals NHS Trust, and Changhai Hospital in Shanghai, China, as previously described (31, 32). Ethical approval was obtained from the respective local ethical committees and The Trent Multicentre Research Ethics Committee. All patients consented to ICGC standards. 159 of these samples passed rigorous quality control for copy number profiles and structural variants and were used in this study.

[0086] DNA was extracted from frozen tumor tissues and whole blood samples (matched controls) and quantified using a ds-DNA assay (UK-Quant-iT™ PicoGreen® dsDNA Assay Kit for DNA) with a Fluorescence Microplate Reader (Biotek SynergyHT, Biotek) according to the manufacturer's instructions. Acceptable DNA had a concentration of at least 50 ng / μl in TE (10 mM Tris / 1 mM EDTA) and an optical density 260 / 280 (OD ) of 1.8-2.0. 260 / OD 280) ratios were shown. Whole-genome sequencing (WGS) was performed at Illumina (Illumina Sequencing Facility, San Diego, CA, USA) or BGI (Beijing Genome Institute, Hong Kong) to a target depth of 50x for cancer samples and 30x for matched controls as previously described (31, 32) (31). Burrows-Wheeler Aligner (33) (BWA) was used to align the sequencing data to the GRCh37 reference human genome.

[0087] Sequencing data generated for this study have been deposited in the European Genome-phenome Archive under accession code EGAS00001000262. Alignments and variant calling were performed using the Wellcome Trust Sanger Institute's Cancer Genome Project (CGP) analysis pipeline, available at https: / / github.com / cancerit / dockstore-cgpwgs. The Battenberg algorithm (34) was used to call clonal and subclonal copy number alterations (CNAs) in all samples; see https: / / github.com / Wedge-Oxford / battenberg. The resulting copy number profiles were subjected to quality control.

[0088] A total of 123 summary measures were generated, including some or all of the following: number of single nucleotide variants (SNVs), number of indels (insertions, deletions, or complex), number of genomic rearrangements (inversions, deletions, tandem duplications, and translocations) and structural variants including mutational signatures, percentage of genome alterations (PGA), DNA breakpoints (and whether they involve linked intra- or inter-chromosomal structural variants), telomere length, number of gene fusions, whole genome duplications (WGDs), presence or absence of any one of kataegis, ETS+ status, and chromosomal thromboplasty, presence or absence of significant driver mutations, copy number altered CNAs (divided into loss of heterozygosity (LOH) and homozygous deletions (HD)), and CNA gains.

[0089] Figure 1b illustrates the input data used as training and validation data. Where applicable, the top number on the Y-axis corresponds to the highest value of the data (e.g., 7887 SNVs) and the dashed line indicates the median. Bar graphs show some of the measured data, namely, number of SNVs, number of indels, number of structural variants including genomic rearrangements and mutational signatures, PGA, DNA breakpoints, telomere length, number of gene fusions. There are grid plots showing the presence or absence of WGD, kataegis, ETS+ status, and chromosophiasis. Finally, there are heatmaps showing the presence or absence of significant driver mutations, as well as CNAs (split into LOH and HD) and CNA gains.

[0090] The summary measurements detailed above form the dataset for further analysis. However, it contains many different data types (binary, categorical, ordinal, continuous), is highly dimensional in terms of number of patients, and will undoubtedly contain highly correlated, co-occurring, or equivalent, events that may confound the analysis. To address this, feature extraction preprocessing is performed prior to the analysis, as described above. Since our downstream analysis is to investigate genomic patterns that exhibit evolutionary behavior, it is important that the results of these analyses can be easily interpreted. This requires a method in which links between input variables corresponding to features are identifiable.

[0091] We outline how each of our summary measures was generated; unless otherwise noted, default parameters were used.

[0092] Number of SNVs, indels, and structural variants SNVs, insertions, and deletions were detected using the Cancer Genome Project Wellcome Trust Sanger Institute pipeline as previously described (31). Briefly, SNVs were detected using CaVEMan with a cutoff “somatic” probability of 0.95. Insertions and deletions were called using a modified version of Pindel (35). Variant allele frequencies of all indels were corrected by locally realigning unmapped reads to the mutant sequence. Structural variants were detected using BRASS (31). The total number of SNVs per sample was calculated, as were the total and type of indels (insertions, deletions, and complex), structural variants (large insertions or deletions, tandem duplications, translocations).

[0093] Clonal and subclonal SNVs Clonal / subclonal quantifies the number of SNVs that are present in all cancer cells (clonal) or only in a subset (subclonal) of the samples, i.e., SNVs with cancer cell fraction (CCF) = 1 and CCF < 1, respectively. These were calculated as previously described by calculating the proportion of reads with the SNV compared to the total number of reads covering that position, followed by adjustments for tumor purity and copy number obtained through the Battenberg algorithm ( 36 ).

[0094] Genome modification percentage This was calculated as the total percentage of the genome affected by CNAs (37). We also recorded the percentage affected by clonal and subclonal CNAs (i.e., CNAs with CCF = 1 and CCF < 1, respectively).

[0095] ploidy We adopted the same approach as previously described ( 36 ), and whole-genome duplication samples were those with an average ploidy, as identified by the Battenberg algorithm, greater than 3. These samples were designated as tetraploid; otherwise, samples were diploid.

[0096] Kataegisu Kataegis were identified using SeqKat https: / / github.com / cran / SeqKat.

[0097] ETS Status A positive ETS status was assigned if a breakpoint between ERG, ETV1, ETV3, ETV4, ETV5, ETV6, ELK4, or FLI1 and a partner DNA sequence was detected and the fusion was in-frame.

[0098] Gene fusion We reported a number of in-frame gene fusions, as well as those affecting only the ETS genes or only TMPRSS2 / ERG.

[0099] Breakpoints Breakpoints were identified using Chainfinder (http: / / archive.broadinstitute.org / cancer / cga / chainfinder) version 1.01. The total number of breakpoints, as well as the total number of linked breakpoints (i.e., when breakpoints are interdependent (38)), the number of linkages, the number and percentage of breakpoints involved in linked events, the number of breakpoints in the longest linkage, and the mean, median, and maximum number of chromosomes involved in linkage. Information on the type of breakpoint, including the number of deletion bridges, intra- and inter-chromosomal events, and the inter- to intra-chromosomal ratio, was also recorded.

[0100] Mutation driver genes A set of driver genes was identified from our previous publication ( 36 ). Using CaVEMan output, we determined nonsynonymous mutations within the exonic regions of these genes as positive events in our dataset.

[0101] Copy number alteration We followed a previous approach (36) to identify consistently aberrant regions. A permutation test was developed in which the CNAs detected from each sample were randomly placed across the genome, and then the total number of hits of the region by each type of CNA in this random assignment was compared to the number of hits of the region in the real data. This process was repeated 100,000 times, and recurrent (or enriched) regions were defined as those with a false discovery rate (FDR) less than 0.05. This was performed separately for gain, loss of heterozygosity (LOH), and homozygous deletion (HD). We first identified small regions, which were merged into larger regions, defined as the union of adjacent regions that all had an FDR less than 0.05. For each sample, the respective datum was set to 1 if a breakpoint corresponding to a gain, LOH, or HD occurred in each region, and 0 otherwise.

[0102] Telomere length Telomere length was estimated as described in a previous publication. (39) A mean correction was applied to the batches to compensate for the effects of chemistry variations in the project.

[0103] Chromosipsis Identified copy number breakpoints were segmented according to inter-breakpoint distance along the genome using piecewise constant fitting (pcf from R package copynumber v1.22.0). Regions with a density higher than 1 breakpoint per 3 Mb were flagged as being dense regions. Chromosophthisis regions may then be defined as regions with a high density of up to three allele-specific copy number states covering a proportion of the region greater than min(1,0.006N+1.1), with the number of copy number breakpoints N>15, a non-random segment size distribution (Kolmogorov-Smirnov test for exponential distribution, P<0.05), and the proportion of each type of structural variant being random with equal probability PTD=PDel=PH2Hi=PT2Ti=0.25 (multinomial test P>0.01), where TD=tandem duplication, Del=deletion, H2Hi=head-to-head inversion, and T2Ti=tail-to-tail inversion.

[0104] Clustering Classification FIG. 2a illustrates a model trained as described below to generate a reduced set of features and score the features with respect to the occurrence of recurrence. In the example of FIG. 2a, the model is a modified restricted Boltzmann machine (44) (RBM) neural network. Latent feature (or latent variable) analysis provides a method to reformulate input data into a reduced set of features that encapsulate the underlying relationships between the original inputs. This framework can be described using a graphical model, such as FIG. 2b, where two latent features contribute to five observed variables. Note that the lack of connections between the latent features indicates that they are conditionally independent. Downstream analysis can then be performed directly on the latent features.

[0105] Many latent feature models have been proposed, each with an associated method of inference on the features (a process called feature learning). These included methods such as nonnegative matrix factorization (41), Bayesian nonparametric methods (42), and neural networks (43). However, none of these known models could meet all of our requirements, so we created a custom RBM neural network.

[0106] RBMs are extensible to multiple data types (45, 46) and can provide interpretable hidden units with appropriate modifications (47). A basic RBM unit consists of only two layers, known as the visible and hidden layers, and one weight matrix used to update both the visible and hidden layers. As shown in Figure 2a, both these layers and the weight matrix are present in the custom-made RBM. The information about the transformation from the visible units (input representation) to the hidden units (feature representation) is encapsulated in the weight matrix. Therefore, we also refer to it as the input-to-feature map. In Figure 2a, the weight matrix is ​​described as fixed, but as explained in detail below, the fixing of the weight matrix is ​​done partway through the training process. Thus, Figure 2a represents the trained state of the RBM.

[0107] The custom RBM in Figure 2a is adapted to compute a discriminative score for each feature, as described in detail below. Thus, the RBM includes an extra classification layer, fully connected to the hidden layer. There is another set of weights that indicate the strength of the connections between the hidden layer and the classification layer.

[0108] The hidden layer typically has fewer units than the visible layer. The basic RBM is formulated as a probabilistic network, which means that each unit represents a random variable rather than a fixed value. All units can only take values ​​of 1 or 0 (active or inactive, respectively), and the input to each unit represents the probability that the unit is active. Thus, the visible layer represents the distribution of observed data values, and the hidden layer represents the distribution of hidden units. Note that the RBM needs to be trained as described below, and the training is done in stages, with each layer being sampled in turn and used to update the weights and parameters of the other layers (the input data is used in the first hidden unit sample). The biases that adjust the baseline activation probability of each unit are not shown for ease of understanding.

[0109] Just as background, note that RBMs are functionally similar to another type of neural network architecture called an autoencoder.(43) The hidden layers of an RBM perform a similar function to the coding layers in an autoencoder, albeit in a probabilistic representation. It has also been shown that RBMs are equivalent to a graphical model of factor analysis.(48) Thus, each hidden unit can be interpreted as a latent feature.

[0110] The standard RBM formulation (44) considers all visible units v = {v i} and hidden unit h={h i} Bernoulli random variables, where v i , h j ∈{0,1}, and each bias is a={a i}, b={b j}, and a i , bj ∈(-∞,∞), weight matrix W, w ij ∈(-∞,∞). Training of an RBM is based on minimizing the free energy of the visible units, since low free energy corresponds to states where the data is well explained through the parameterization of the model. The energy-based probability distribution is of the form

[0111]

number

[0112] Take where E(v,h) is the energy function and Z is a normalization factor, which is the probability of observing the joint v,h pair. The energy function of RBM is E(v,h)=-a T vb T hv T Wh (2) is given as:

[0113] In this formulation Z=Σ v Σ h e -E(v,h) , (3) and This is difficult to compute due to the number of possible combinations of v and h.

[0114] Training is performed on the energy in visible units, so we consider a dataset D={d k , k=1, 2, ..., K} k To calculate the likelihood of observing a visible unit corresponding to , we need to marginalize on h in Eq. 1. The likelihood is

[0115]

number

[0116] It is calculated using where θ∈{{ai},{b j},{w ij}} is the complete parameter set. To simplify notation, we use L(θ|v=d k ) to L(d k ) for each parameter we want to update. k )) / ∂θ. The partial derivative of the logarithm in Equation 4 has the form

[0117]

number

[0118] Take.

[0119] We then run

[0120]

number

[0121] Calculate the expectation using This can be used to update the model parameters via gradient descent.

[0122]

number

[0123] corresponds to the expected energy state that would be invoked from observing a data sample, and the term

[0124]

number

[0125] are the expected energy states of the model configuration, and both depend on the current model parameters. So, they are, respectively,

[0126]

number

[0127] and

[0128]

number

[0129] When we calculate the partial derivatives with respect to the parameters, we get

[0130]

number

[0131] is obtained, This is the update equation

[0132]

number

[0133] is used to set up the learning rate ν and η.

[0134]

number

[0135] The value can be easily estimated by taking the arithmetic mean.

[0136] term

[0137]

number

[0138] The terms are generally hard to compute since they involve summing over all possible configurations of v and h. An alternative method is to perform Gibbs sampling using conditional probabilities since these are fairly easy to compute due to the conditional independence between units in the same layer. We can estimate the conditional probability of the hidden layer values ​​from the visible layer and vice versa, as follows: (h|v)=Π j P(h j |v) (14) P(v|h)=Π i P(v i |h) (15)

[0139] P(h j |v) and P(v i The form of |h) depends on the activation function, which takes as input the product of a unit in a layer with its corresponding weight, and outputs the probability that the unit is active. In this work, we use the logistic sigmoid (or simply "sigmoid") function, which is

[0140]

number

[0141] is given by where x depends on the layer we sample from, and so the individual hidden and visible probabilities are P(h j |v) = σ(b j +Σ i v i w ij ) (17) P(v i |h) = σ(a i +Σ j h j w ij ) (18) can be written as:

[0142] The sample is appropriately j|v) or P(v i The hidden units are extracted by setting the corresponding units to 1 with a probability given by the value of |h). These can then be used to compute estimates for P(v) and P(h) by marginalization over the condition variables. In practice, full Gibbs sampling at every update iteration would be prohibitively slow, and therefore we used an approximation called contrastive divergence (44), in which a 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(1)), where the data (or a mini-batch of data) is presented as a matrix and used to sample the hidden unit values, which 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) (49). Thus, information travels in both directions with respect to the weights in these early stages.

[0143] During training, the results of these updates are stored in three matrices (H, W, and V) that correspond to the weights as well as the network representation of the tumor data in the visible and hidden layers. These matrices correspond to the network reconstruction of the data (visible layer, V), the latent feature representation of the data (hidden layer, H), and the input feature mapping (weights, W). When the network is trained, these can be extracted and utilized in analysis.

[0144] A number of simple modifications were made to the standard RBM to ensure that the representation is interpretable, generalizable, stable, and reproducible. These modifications include data integration, use of non-negative weights, pruning of hidden units, sparsity, avoidance of overfitting, and convergence to a global solution. All modifications are incorporated, as described below, but it will be understood that alternative versions may incorporate some but not all of the modifications. Data Integration

[0145] Our data consist of multiple different modalities, and unlike traditional multiomics approaches that have a large number of data points from a small number of sources, we have a small number of data points from a large number of sources. Therefore, data integration had to be carefully considered. RBMs can be modified to incorporate inputs of multiple modalities, sometimes through modifications of the energy function (50, 51). However, we decided to avoid this complication and standardize all our inputs by ranking all integer and continuous variables before rescaling them to [0,1]. For example, the specific transformations incorporated are: Binary -- set as {0, 1}. Category--One Hot Encoding*. Integer--rank and scale to [0,1]. Continuous -- ranked and scaled to [0,1].

[0146] For the integer and continuous cases, we used ranking because it decouples the values ​​from the distribution of inputs and, after scaling to [0,1], the new value can be interpreted as the probability that the corresponding visible unit is active. Therefore, all inputs are treated equally in the processing of the RBM. These transformations do not affect the hidden units, which are the Bernoulli random variables h i ∈{0,1}. In one-hot encoding, categorical variables are replaced by vectors of length equal to the number of categories. The values ​​of the vector are all zero except for a 1 in the nth position indicating membership in the mth category.

[0147] Non-negative weights Neural networks are considered a black-box approach because the transformations they perform are highly complex. To improve the interpretability of the network's behavior, we impose non-negativity constraints on the weight updates, specifically by penalizing negative values. We use an approach in which for each negative weight, a quadratic gate function is subtracted from the likelihood (47). Mathematically, this is

[0148]

number

[0149] It is written as follows: Here, α represents the strength of the penalty.

[0150]

number

[0151] It becomes.

[0152] This is the update rule

[0153]

number

[0154] Leads to.

[0155] W {-} is a matrix that contains the negative entries of W, with other parts being zero. This formulation is equivalent to an L2-norm penalty for negative weights, and therefore imposes a greater penalty on weights that are more strongly negative. When used in a training scheme, this forces the network weights to be non-negative, simplifying the interpretation of the input feature maps. This can be thought of as a non-linear extension of non-negative matrix decomposition (41), which can similarly be used to represent the underlying structure of data in terms of parts that are features in machine learning terms.

[0156] Since the weights can no longer be traded off against each other with countervailing weights of opposite sign, this means that the lowest free energy states correspond to those with the least redundancy, and thus during training the hidden units compete to convey information about a single input (52). This means that the input is represented by only a small number of latent variables, and thus when the initial number of hidden units is of a similar order of magnitude to the number of data inputs, this results in some of the biases or weights converging to negligibly small values, and the corresponding hidden layer activations converging to arbitrary fixed values. The latter are then called dead units. This is of fundamental importance to our method, as it can be used as an estimate of the intrinsic dimensionality of the data.

[0157] Hidden unit pruning During training, we prune dead units to improve the speed of the algorithm. However, in probabilistic networks such as RBMs, it is not easy to determine dead units because the values ​​in the network at each state change probabilistically. To get around this, we apply L 1 / 2 We apply a norm penalty, which penalizes nonzero activation values ​​(53). This forces the values ​​for all patient samples to be zero rather than some arbitrary value, which can be easily identified and removed by a thresholding approach. This penalty function is computed over all training data samples, so, for consistency with Equation 4, we compute the likelihood for each sample as

[0158]

number

[0159] This can be formulated as follows: Here, f(y k )=P(h|v k), and β is a parameter that describes the strength of this penalty. We compute the gradient of an additional likelihood term with respect to the bias of each hidden unit, which is given by the formula

[0160]

number

[0161] is given by:

[0162] We then define the vector of gradients for all hidden unit biases as

[0163]

number

[0164] Therefore, the corresponding update rule is

[0165]

number

[0166] It can be written as:

[0167] In our learning algorithm, we prune dead units every 50 iterations after the first 1000 iterations.

[0168] Sparsity Sparsity is a desirable property for a latent space representation, since it means that information is conveyed in a concise form. The penalty measure defined in Equation 22 introduces sparsity by penalizing highly active hidden units, forcing the network into a sparse configuration (53). No further sparsity measures were used in training, since the weight matrices, which define the input to the feature mapping, are heavily filtered in later stages.

[0169] Avoiding overfitting A concern with neural network formulations is the tendency to overfit the data, which in this application results in a feature set that does not represent the true underlying structure and is therefore not generalizable. To mitigate this, we have taken several measures, for example: 1. DropConnect, 2. Max-norm regularization, 3. Bootstrap aggregation, 4. Early termination.

[0170] In DropConnect(54), at each training iteration, a given percentage of weights in the network are randomly set to zero with uniform probability. This helps to prevent overfitting by temporarily decoupling correlations between features, thus increasing the likelihood of learning features that are independent of the state of other features.

[0171] When using max-norm regularization (55), we set an absolute value on the norm of each weight vector that forms the input to a single hidden unit. If the vector becomes too large, we rescale it to obey the constraint. It is possible that non-negative weights will continue to grow throughout training, because the binary nature of some inputs means that they are already within the maximum output of the sigmoid activation function when present, and so it does not matter whether the values ​​are accurate or not. Max-norm regularization prevents this from happening, facilitating comparisons between weight matrices from different runs.

[0172] For bootstrap aggregation(56) (bagging), multiple networks with the same initial architecture are trained on subsets of the data and the outputs are merged. In our feature learning representation, we extract the weight matrices from each of the networks and merge them according to the cosine distance between the features, as shown in Figure 4 and described in more detail below.

[0173] Finally, when implementing early stopping (57), we need to compare the performance of the network on the training set with its performance on an unseen validation set. If the network performs similarly on the training and validation sets, it is a good indicator that it returns generalizable outputs. We start with the subsets extracted for ensemble learning and use the data left out when the subsets were sampled as a validation set, which is propagated through the network. Since RBMs are formulated as energy-based models, early stopping is predicted by comparing the free energy in the training set with the free energy in the validation set (58). If the free energy resulting from the training set becomes consistently lower than that of the validation set, overfitting has occurred and training is stopped.

[0174] Convergence to a global solution When training multiple networks and merging their results, it is important that each network converges to a global solution, otherwise the results will be inconsistent. Furthermore, because RBMs are trained by stochastic gradient descent, it is possible that the algorithm may get stuck in a local optimum. To minimize this possibility, we used a cyclic learning rate scheme (59) in which the learning rate for each of the variables oscillates between zero and a maximum value throughout training. Since the maximum value undergoes decay, the maximum learning rate decreases to zero throughout training. This approach has been shown to aid convergence to a global solution, and has the advantage that the learning rate parameter does not need to be tuned (59).

[0175] We trained 2000 networks using 75% of the data as a training set (selected uniformly at random). The remainder of the data was used as a validation set for early stopping. If early stopping occurred, the entire network was discarded (as it may not have had time to converge to an accurate feature representation) and another network was trained in its place, and this was repeated until training was successfully completed. Figure 3 illustrates that the dimensionality of the extracted features trained above is consistent with a mean of 26.30 and a standard deviation of 1.51.

[0176] As explained above, multiple networks were trained with this data, with each individual network run providing a similar, though not identical, weight matrix. Therefore, the weight matrices from each network run were merged and filtered to form the final input feature map. Because the number of features, the inputs they represent, their magnitudes, and order do not necessarily occur the same in each network, we constructed an algorithm based on cosine similarity.

[0177] Figure 4 illustrates an outline of the steps of the algorithm. Each heatmap in Figure 4 shows the relative magnitude of the network weights corresponding to the map between each input and each feature. The individual weight matrices on the left of the figure are concatenated to form the large matrix in the middle of the figure. For each input, co-occurring inputs and their relative magnitudes are calculated to form a preliminary feature set. Cosine distances are then calculated pairwise between new features and used to merge features that were within a threshold of similarity.

[0178] An example of pseudocode for merging weight matrices is shown in the following algorithm. Algorithm 1: Input: A set of weight matrices from each network run Concatenate the weight matrices to form matrix W, Set small weights to zero, Set the similarity threshold τ = 0.5. Initialize the matrix M with the same number of rows and columns as the number of inputs, Initialize an empty feature matrix F, for i=1 to number of inputs do Set the i-th row of M equal to the average of all rows of W with i-th weight >0 end Compute the pairwise cosine similarity matrix I from M, while number of lines in S>0 do Read the first row of S as the current similarity vector s, s j Identify all j for which > τ; Add the mean of M[j,:] as a row to F, Remove the jth row from S and M, end Rescale all rows in the feature matrix F by max-norm,

[0179] Low magnitude weights were weights less than 50% of the maximum weight value for each hidden unit. The combined weight matrix has 30 features as opposed to 22-31 in each individual run. This is primarily due to low frequency inputs not being consistently represented after the data is subsetted for cross-validation.

[0180] Referring again to Fig. 2a, the merged weight matrix is ​​used as the fixed weight matrix. A feature representation of the data is then obtained using the network of Fig. 2a, but utilizing all patient samples and the fixed weight matrix. This was done by initializing the weights to the merged weight matrix and setting the weight learning rate to zero. Bias learning was enabled since it may be different from the bias in the previous network due to the removal of low magnitude weights.

[0181] After the remaining network parameters converge during training, performing further iterations is equivalent to sampling the hidden units / representations for each patient. Therefore, we averaged the hidden unit values ​​taken every 10th iteration during the final 1000 iterations to obtain the final representation.

[0182] FIG. 5a illustrates a heatmap showing the relationship between patients and input data. FIG. 5b illustrates how the data is converted to a feature representation as described above. FIG. 5b shows that the 123 inputs have been reduced to 30 features, as set forth in the table showing the reduced set of features and described with respect to FIG. 1a.

[0183] Two-stage clustering The dimensionality of the feature representation in Figure 5b is still rather large relative to conventional clustering techniques. Therefore, we adopted a two-stage approach, first clustering by the features most informative for clinical outcome, computing the centroids of these first-stage clusters for all features, and then clustering these in a second stage of clustering to produce the results shown in Figure 6. Details of the identification of informative features using discriminative scores and the clustering method used are described below.

[0184] Discrimination score Several methods have been proposed to quantify the relative importance of units in neural networks (60). However, most of these are generally formulated to find the inputs that are important in discriminating the output (61, 62). For our application, we wish to quantify the discriminatory power of each of the features (hidden layer) with respect to clinical outcomes. Because we utilize nonnegative weights to determine the relevance of inputs to hidden units in feature extraction, a similar approach can be adopted to determine the importance of hidden units to outcomes.

[0185] As briefly explained above with reference to Fig. 2a, the architecture of the base RBM is similar to ClassRBM (63), modified so that a discriminative score for each feature can be obtained. Thus, as shown in Fig. 2a, our RBM includes an extra classification layer fully connected to the hidden layer, whose units contain class values. We want to uncover the underlying relationships in the data (encapsulated by the features) in an unbiased way, and then determine to what extent these features are related to clinical outcomes. Therefore, we enforced that the classification weights are unidirectional, and the information used in training is passed only from the hidden layer to the classification weights. This means that the latent structure encapsulated by the hidden units remains unbiased by knowledge of clinical outcomes, and the algorithm for feature learning can still be considered unsupervised. (In contrast, in ClassRBM there is another set of weights that indicate the strength of the connections between the hidden layer and the classification layer, and these are trained in the same bidirectional manner as the input weights.)

[0186] Furthermore, we enforced non-negativity constraints on these class weights, as well as on the input weights. Thus, when trained, the relative magnitudes of these class weights quantify how important each corresponding feature is for discriminating the corresponding clinical outcome, in a manner similar to standard non-negative matrix factorization. We define the discrimination score s as the absolute value of the weight corresponding to recurrence minus the weight corresponding to no recurrence, which can be expressed as s=|c r -c r' | (26) It can be expressed as where c r is the class weight associated with the recurrence, and c r' is the class weight associated with no recurrence.

[0187] These s-values ​​can be considered heuristic and quantify the importance of the corresponding feature to the clinical output, similar to how component loadings quantify the explained variance of the corresponding principal components in principal component analysis (PCA). Since there is no prescribed rule for determining the number of features, we followed an approach similar to that traditionally used in PCA and selected the number of features using cumulative distribution. We selected a cutoff value of 0.9 of the cumulative discrimination score sum, which resulted in 14 features out of 30 features being selected in the initial clustering stage. These 14 features are listed below and highlighted (red) in Figure 6.

[0188] [Table 12]

[0189] Clustering Tumor clustering was performed on the latent feature representations in a two-step process to facilitate the identification of clusters associated with clinical outcomes. Because the feature representation for each patient can be thought of as a vector containing the probability that the corresponding feature is active, it is appropriate to use a distance measure that quantifies the distance between the probabilities. Therefore, we calculated the average Jensen-Shannon (JS) divergence (64) between tumors pairwise.

[0190] Hidden layer h A and h B For a pair of patients A and B represented by the latent feature representation of

[0191]

number

[0192] can be written as Where:

[0193]

number

[0194] h A and h B The additive terms in the brackets in Equation 27 represent the Kullback-Leibler divergence between each element of the latent feature representation for any patient and the corresponding element of the midpoint vector m.

[0195] As we do not use a Euclidean distance metric, clustering via k-means is not appropriate, and therefore in the first stage we use k-medoid clustering, which is similar to k-means but selects a representative data point (medoid) as the centroid for each cluster instead of the mean. We determined that 11 clusters were optimal using the silhouette method (65). In the second stage of clustering, the medoids themselves were clustered using hierarchical clustering (again using JS divergence), which was used to generate and order the clusters according to the dendrogram shown in Figure 6.

[0196] The discrimination scores quantifying the relevance of each feature in predicting recurrence are shown as a green heatmap in Figure 6. The 14 features (listed in red in the table of features for the initial clustering stage) are used as input for k-medoid clustering with 11 clusters (determined by the silhouette method).

[0197] Medoids of each cluster were used as input for hierarchical clustering using all features, which revealed two main metaclusters, MC-A and MCB, with distinct profiles. Metacluster MC-B was further divided into MC-B1 and MC-B2 as shown in the dendrogram. The main heatmap shows the medoid feature values ​​for patients in each cluster, ordered by the hierarchical clustering (right scale). The colors of the metaclusters are indicated in the text above the dendrogram.

[0198] Thus, metacluster A (MC-A) can be identified by samples with intrachromosomal structural variants, SPOP mutations, chromosomal smears, and loss of heterozygosity (LOH) in regions 5q15-5q23.1 (spanning CHD1) and 6q14.1-6q22.32 (MAP3K7, ZNF292). Metacluster B1 (MC-B1) can be identified by samples with ETS fusions and loss of heterozygosity (LOH) affecting 17p (TP53) and regions 19p13.3-19p13.2, 22q11.21-22q11.22. Metacluster B2 (MC-B2) can be identified by samples with frequent ETS fusions, interchromosomal linked structural variants and loss of heterozygosity (LOH) affecting 17p (TP53) and regions 5q11.1-5q14.1 (IL6ST, PDE4D) and 10q23.1-10q25.1 (PTEN).

[0199] ARBS classification (classification based on the proximity of the DNA breakpoint to the androgen receptor binding site) To examine the proximity of DNA breakpoints to androgen receptor binding sites (ARBS), we designed a permutation approach to quantify deviations from a random distribution of breakpoints across the genome. We downloaded processed ChIP-seq data targeting AR for 13 primary prostate cancer tumors from Gene Expression Omnibus (accession no. GSE70079) (66) and merged them to use as an alignment of ARBS.

[0200] To detect significant departures from a uniform random distribution, we use the observed and permuted data (B obs and B. perm The percentage of breakpoints within 20,000 base pairs (bp) of the ARBS was calculated for each gene. obs >P 97.5% (B perm ), the tumor is classified as enriched; otherwise, Bobs <P 2.5% (B perm ), the tumor was classified as depleted. Otherwise, the difference was not significant and the tumor was classified as undefined. The levels of enrichment or depletion of breakpoints near the ARBS used in Figure 7a were calculated using the formula

[0201]

number

[0202] was estimated according to

[0203] The method was validated using the same data used to train the modified RBM described above. Figure 7a shows the results of calculating the proportion of DNA breakpoints within 20 kilobases (kb) of the AR binding site for each patient in 159 samples. For each of our 159 samples, we randomly shuffled the observed breakpoints across the genome (GRCh37) masked for assembly gaps (AGAPS mask) and intracontig ambiguity (AMB mask) 1000 times using the R package RegioneR (67). In Figure 7a, the number of breakpoints is normalized by the number of proximal breakpoints expected by chance. Each of the tumor samples is ordered according to this normalized proportion. The class (enriched, depleted, or indeterminate) was determined based on whether the tumor exhibited more proximal breakpoints than expected (enriched), fewer proximal breakpoints than expected (depleted), or no statistically significant difference (indeterminate).

[0204] Figure 7b shows a heatmap of genomic features for each patient using the ordination from Figure 7a. The genomic features include genetic alterations associated with features previously identified from the modified RBM. As shown, the depleted tumors had the highest genomic alteration percentage (PGA), and the highest frequency of multiple CNAs, chromosomal thrombosis, kataegis, and SPOP mutations (Relationship column, Figure 7b). Enriched and indeterminate tumors did not show significant differences for any CNAs, but both showed a higher frequency of CNAs covering PTEN and TP53 than the depleted group (Relationship column, Figure 7b). For ETS fusions and interchromosomal / intrachromosomal cSV ratios, the enriched group showed higher enrichment than the intermediate group, which in turn showed higher enrichment than the depleted group. Both enriched and depleted tumors showed a higher number of breakpoints than the indeterminate tumors. Associations with ARBS pairs were established with a one-sided Mann-Whitney U test with P<0.05.

[0205] In Figure 7b, statistically significant relationships for the three classes are shown in the "Relationship" column, with E, D, and I indicating enrichment, depletion, or indeterminate, respectively. Brackets {.,.} indicate no relationship between the enclosed classes, but they both indicate significant differences for the remaining classes. Relationships are ordered such that the leftmost class shows a significantly greater proportion of genetic modifications. For Bernoulli variables, significance was determined by chi-squared tests followed by Fisher's exact tests for each pairwise relationship, and for continuous variables, Kruskal-Wallis tests with Tukey's HSD method were used (adjusted P<0.005 for all tests).

[0206] Figure 7c shows the ARBS groups in two additional datasets compared to 159 samples from the ICGC UK dataset. The aim of this analysis was to validate the ARBS findings within additional datasets. The first set is a set of low- to intermediate-risk tumors from the Canadian Prostate Cancer Genome Network (CPC-GENE) (12), and the second set is a set of high-risk tumors from the Melbourne Prostate Cancer Research Group, Australia (unpublished). The top left bar graph shows the proportion of each ABRS group in each country's data. The main figure shows the results of clustering these groups by the proportion of CNAs. We found that the depleted groups clustered together across all datasets (P < 0.0337, approximately unbiased multiscale bootstrap method).

[0207] ARBS clusters were identified by a bespoke permutation test with multiple testing correction. For example, an agglomerative hierarchical clustering of ARBS groups across the Australian, Canadian, and UK datasets was generated using the R package pvclust (68) v2.0.0 using the ward.D2 clustering method with squared Euclidean distance (100000 iterations). This package also allowed the estimation of approximately unbiased multiscale bootstrap (AU) P values ​​for the depleted groups. These clustering results were confirmed by a partitional clustering approach using the R packages cluster v2.1.0 and factoextra v1.0.5.

[0208] Classification by Ordering A consensus ordering of events has been previously determined by estimating a phylogenetic tree from the cancer cell fractions (CCFs) containing each aberration and applying the Bradley-Terry model to determine the most consistent ordering of events (36). This approach has many sources of uncertainty. In particular, we are often unable to estimate the true phylogenetic tree for each patient, and furthermore, it is not possible to determine the relative timing of events on parallel branches. However, we can estimate a set of possible trees using the relative cancer cell fractions (CCFs) of the involved genomic aberrations, and from these we can estimate a set of possible orderings. Thus, we created an algorithm that samples a single possible tree from the data and uses this to sample feasible orderings of events for each patient. This is repeated multiple times, thereby encapsulating the uncertainty in these estimates in the output distribution. This type of algorithm is called Monte Carlo simulation, and emphasizes the use of randomness in the procedure.

[0209] In this application, we have adopted an extension of the Bradley-Terry model known as the Plackett-Luce model (69, 70) as the basis for our ordering analysis. This model is used to construct a probability distribution over the relative ranking of a finite set of items, whose parameters can be estimated from a large number of individual rankings. This can be used to quantify the expected rank of each item with respect to other items across a population. In our application, the items correspond to events, i.e., the emergence and fixation of novel copy number alterations (CNAs) as specified in the extracted features. The ranking of these events is therefore related to the order in which they are expected to occur. We also utilize a Plackett-Luce mixture model, which allows for subpopulations of data with different orderings.

[0210] Plackett-Luce model Given a set of CNA expression events for each patient with associated subclonality, we would like to infer the order in which these events typically occur. To do this, we used the Plackett-Luce model, which is formulated as a ranking method and returns a value that quantifies the ranking preference. We use a different interpretation, namely, ordering, which is defined as the inverse of the ranking preference (71). Like the Bradley-Terry model, the Plackett-Luce model does not return any temporal information other than the expected order of events.

[0211] We define a set of N copy number events of interest as C = {c1, c2, ?, c N} (29) having We can apply Luce's axiom of choice (69), which states that the probability of selecting one event from a set of events over another does not depend on the presence or absence of other events in the set. Thus, we can define the probability of observing event i as

[0212]

number

[0213] can be written as Here, {α i} is a coefficient that quantifies the relative probability of observing the i-th event. To reflect the ordering aspect of our application, we call this value proclivity. Plackett (70) used this formalism to construct a generative model in which all N events are randomly sampled without replacement (i.e., permutations) from C. Let us denote Λ as a permutation of a set C and λ as k ∈C and

[0214]

number

[0215] If we make the correspondence such that

[0216]

number

[0217] Write it like this: Here, α λk is the event λ k is the proclibity associated with Λ (k) ={λ k ,λ k+1 ,...,λ N} is the set of possible events after k-1 events have occurred.

[0218] Plackett-Luce mixture We hypothesized that there may be multiple sets of copy number orderings present in a population, and therefore analyzing all events with one ordering scheme may not be appropriate. Furthermore, inhibition of AR-associated breakpoints means that some CNAs are found more frequently together with a selected set of other CNAs, which contradicts Luce's selection axioms. Therefore, we implemented a mixture modeling approach (71, 72) that revives Luce's selection axioms, since the selection of each CNA can be considered conditionally independent of the mixture components. In such a finite mixture model, we assume that the population consists of G subpopulations. In this setting, the ordering Λ for the sth sample is s The probability of observing

[0219]

number

[0220] and Here, ω gare weighting parameters (not to be confused with the weighting matrix described above) that quantify the probability that sample s belongs to subgroup g. Appropriate parameter values ​​can be determined using maximum likelihood estimation via the EM algorithm (72).

[0221] The number of mixture components may be selected using a Bayesian Information Criterion (BIC) estimate, which is BIC=Nlog(M)-2l(Θ ML ) (33) is given by Here, Θ ML is the parameter set that maximizes the log-likelihood l(·), N is the number of parameters, and M is the number of samples.

[0222] The general formulation of the Plackett-Luce model takes as input a matrix containing the sequence of events for each patient. However, we do not know the order in which these events occurred, only the presence of each CNA and the cancer cell fraction (CCF) for each patient. Therefore, we first estimate a phylogenetic tree for each patient and then determine the order of events from this. Since we have only one tissue sample for each patient, there is often uncertainty in the tree topology and possible sequences of events, and therefore we use a Monte Carlo sampling method to sample the trees and sequences of events and 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.

[0223] Another problem arises from censoring, which occurs when samples are taken before all possible anomalies have occurred, resulting in missing data. These are called partial orderings in the Plackett-Luce framework, and a common approach to deal with this is to reformulate the model so that all missing events are implicitly ranked lower than the observed data (72, 73). This may not be appropriate for our analysis since there may be multiple subgroups, and we expect that different anomalies have similar or comparable effects in each subtype and therefore rarely co-occur despite representing the same type. For example, the absence of a very early anomaly could be due to the occurrence of another, less frequent anomaly, so including it at the bottom of the ordering would bias the ranking towards the more frequent anomalies. Therefore, our algorithm works in two stages: 1. Determine the number of mixture components and assign patients to each component. 2. Estimate the ordering profile for each component. These differ because we handle the creation of phylogenetic trees in each of these processes slightly differently to account for censoring. When estimating the number of components, we compute the tree using only the observed CNAs. However, when estimating the full ordering profile, we introduce another sampling step into the Monte Carlo method, where we explicitly sample a number of additional CNAs with a probability proportional to the subclonality of the anomalies in each mixture component tumor. Sampling in this way reduces the bias toward more frequent anomalies.

[0224] Assigning samples to mixture components In the first phase, we will do the following: 1. Sample the dendrogram for each patient. 2. Sample the sequence of events for each patient that matches the tree. 3. Calculate the Bayesian Information Criterion (BIC) for 1 to 10 mixture components. 4. Repeat steps 1 to 3 1000 times. 5. Determine the number of mixture components that consistently had the lowest BIC scores. 6. Assign patients to mixture components.

[0225] The phylogenetic tree is constructed by first sorting the CNAs of each patient in descending order of CCF obtained from the output of the Battenberg algorithm, and then iterating over them, sampling possible parents with uniform probability. Since the CCF of a parent cannot be less than the sum of the CCFs of their children, viable parents are defined as parents whose CCF is greater than or equal to the CCF of their current child plus the CCF of the CNA under consideration. The position in the sequence where the CNA occurred is sampled with uniform probability as any position after the parent. For the estimation of the ordering and assignment to mixture components, the R package PLMIX was used as it incorporates mixture models and partial ranking (so the absence of a CNA from a sequence does not penalize its position in the ordering). A vector of assignments was kept for each sample run, and the final assignment was determined by the most frequent assignment over the course of 1000 runs.

[0226] Bayesian Information Criterion (BIC) scores were determined for each mixture component for each of the 1000 runs, which are shown in Figure 8. The y-axis shows the BIC score calculated for each ordering when there are 1 to 10 mixture components as indicated on the x-axis. Each individual score is shown as a cross (blue) and the average of the scores for each component is shown as a line (red). Since the BIC score was lowest for the two mixture components for all orderings sampled, this was considered the value to use in subsequent analyses.

[0227] Estimate the ordering profile of each component In the second phase, we will: 1. Sample the dendrogram for each patient. 2. Sample the sequence of events for each patient that matches the tree. 3. Augment sequences with additional CNAs to reduce censoring bias. 4. Compute the ordering profile for each mixture component. 5. Repeat steps 1 to 5 1000 times. 6. The results are merged to determine the final ordering profile for each mixture component.

[0228] The phylogenetic tree and sequence of events were first determined as before. However, instead of utilizing partial rankings in the PL model, we explicitly augmented the data with additional CNAs to account for those not observed due to censoring. The probability that a CNA is added to the sequence of events is equal to the proportion of subclonal occurrences with respect to the total number of occurrences in the subpopulations defined by the mixture components. This means that

[0229]

number

[0230] can be written as Here, N sub (·) and N total (·) is the CNA c in mixture component g i The subclonal and total occurrences of 10 ...

[0231] Calculating these values ​​using patient samples for each mixture component, rather than the entire population, means that only the CNA subclonalities associated with each subpopulation are considered. Imputation is done by drawing a uniform random number r for each patient,

[0232]

number

[0233] This is done by including the CNA in a set of additional CNAs for each patient if . The set of additional CNAs for each patient is uniformly shuffled and added to the sequence. Imputation helps to mitigate truncation. We then compute orderings for each mixture component separately using a Plackett-Luce model without partial ranking. This process is repeated 1000 times and the Plackett-Luce coefficients for each CNA are calculated and used to create an empirical distribution of the Plackett-Luce coefficients for each CNA, which is used to create the boxplots in Figure 9a.

[0234] Figure 9a shows the proportions of 159 samples for the Plackett-Luce coefficient for ordering I and ordering II. Phylogenetic trees from individual tumors were used to estimate two ordering profiles using a Plackett-Luce (PL) mixture model, as explained above. Tumors were assigned to ordering I (top) or ordering II (bottom). Horizontal box and whisker plots (5th / 25th / 75th / 95th percentiles) show the negative Plackett-Luce coefficient α for the i th genetic alteration (x-axis). i represents the bootstrap estimate of i The lower the value of , the more likely it is that a genetic alteration will occur early. The y-axis indicates the proportion of samples in the mixture component in which a genetic alteration is observed. Genetic alterations with a proportion higher than 0.25 have chromosomal regions annotated with notable driver genes within the region shown in brackets. The color of the box and whiskers indicates the chromosome in which the abnormality occurred.

[0235] Figure 9a shows that the two orderings show striking differences. Tumors corresponding to ordering I frequently underwent early LOH of 8p (spanning NKX3.1) and ETS fusions. Less frequent LOH of regions covering RB1, BRCA2, CDH1, TP53, or PTEN genes could also occur. The profile showed occasional very early LOH of 1q42.12-42.3. Tumors corresponding to ordering-II consistently showed early LOH events covering MAP3K7 and 13q (EDNRB, RB1, BRCA2), as well as copy number gains. However, the earliest events, mutations in the SPOP gene and LOH covering CHD1, were less frequent. Both orderings showed late gains of chromosome 19.

[0236] Figure 9b shows the variation in the order of copy number alterations among individuals from 159 samples. When comparing the occurrence of aberrations among individuals in each ordering, we found that the relative order of alterations is highly variable, indicating that it occurs stochastically. The value at the left end of each bar is the lowest Plackett-Luce (PL) coefficient among all CNAs that must have occurred after the genetic alteration named at the left (i.e., were found to have occurred subclonally (CCF<1) when the named CNA was observed in all sampled cells (CCF=1)). The value at the right end of each bar is the highest PL coefficient among all CNAs that must have occurred before the genetic alteration named at the left (i.e., were observed in all sampled cells when the named CNA occurred subclonally (CCF=1)). The black dots represent the PL coefficients of the CNAs named at the left. CNAs are ordered from top to bottom by their PL coefficient values.

[0237] Comparison of three classification methods The table below demonstrates the agreement of the three classification methods described above by showing which of the 159 samples were assigned to each category.

[0238] [Table 13]

[0239] [Table 14]

[0240] The above table reveals a remarkable relationship: MC-A is a large subset of the depleted group (22 / 27), and both are near-complete subsets of ordering II (26 / 27 and 30 / 32, respectively). We can therefore infer that there exists a subset of tumors that exhibit all the corresponding characteristics, i.e., evolutionary trajectories (ordering-II), breakpoint mechanisms (ARBS: depleted), and characteristic patterns of aberrations (meta-clusters:MC-A). Thus, to classify by evotype, we adopted a majority-vote approach and defined tumors assigned to at least two of MC-A, depleted, or ordering II as belonging to alternative evotypes, distinguishing these from canonical evotype tumors that can evolve via trajectories involving canonical AR processes.

[0241] Figure 10a plots progression-free survival versus time for patients with tumors classified as either evotype. The plot is a Kaplan-Maier plot, with P-value (0.0218) and hazard ratio (HR) calculated using the log-rank method. HR is estimated with a 5th to 95th percentile range of -2.26 (0.964 to 5.3). As shown, patients with alternative evotype tumors had a poor prognosis. The endpoint is time to biochemical recurrence.

[0242] This poorer prognosis is perhaps surprising given that other clinical characteristics such as tumor stage, ISUP Gleason grade group, and PSA (ng / ml), plotted in each of Figures 10b-d, indicate that no statistical differences are observed between the two classifications. The p-values ​​of the chi-square test are P=0.5968 for the results in Figure 10b, P=0.0586 for the results in Figure 10c, and P=0.191 for the results in Figure 10d. All clinical signs were collected at the time of prostatectomy.

[0243] Figure 10e is a bar graph showing the prevalence of each genetic abnormality in each evotype. Evotype classification is determined using majority consensus. Abnormalities with significant differences between evotypes (P<0.05 using Fisher's exact test) are listed below (colored in red for alternative evotypes and blue for canonical evotypes). Thus, each evotype is characterized by a different propensity for several abnormalities in combination, but note that no single abnormality was necessary or sufficient for assignment to any evotype.

[0244] [Table 15]

[0245] [Table 16]

[0246] Statistical model of evotype convergence Figure 11a is a flow chart of a statistical algorithm for obtaining the probability of convergence to a canonical or alternative evotype based on the accumulation of genetic alterations. An example output from this algorithm is shown in Figure 11b. We assume that the accumulation of such anomalies in individual tumors followed a stochastic process in which the order and relative timing of the anomalies occurs with some degree of randomness / stochasticity. Similar to the ordering analysis (described above), we utilized a statistical algorithm that simulates a large number of possible anomalies consistent with possible phylogenetic trees and then estimates the probability that tumors carrying these anomalies will converge to the canonical evotype (the probability of convergence to an alternative evotype is 1 minus the probability of convergence to the canonical evotype). The algorithm is repeated with an increasing number of anomalies (loop i) and several Monte Carlo iterations of the ordering samples (loop j).

[0247] The accumulation of aberrations in tumors is modeled as a Poisson process (74). Figure 11a shows that the first step in each iteration of the first loop, i, determines the average number of aberrations across all patients in the i-th iteration, x i (step S1000). This is then used as an input parameter to a Poisson random number generator to draw the number of anomalies n to be sampled in each iteration of loop j. In other words, the number of anomalies n considered in this iteration is randomly sampled (step S1002).

[0248] We then identified tumors with enough anomalies and selected one with uniform probability (step S1004). The data for the selected tumor is then used to sample a phylogenetic tree using the relative CCF of the anomalies (step S1006). The phylogenetic tree is sampled from the anomaly data. We then used the phylogenetic tree to sample the order of occurrence for the anomalies (step S1008), retaining the first n (hence the output from step S1002 is used in this step, as illustrated by the connecting arrows). In other words, using both n and the phylogenetic tree obtained in the previous step, a set of anomalies is sampled that is consistent with the possible ordering of events allowed by the phylogenetic tree. The set of anomalies generated in this step is A j ={a1,a2,_,a n The abnormalities used are SPOP mutations and CNAs identified in feature extraction; intra / inter chromosomal breakpoints, ETS conditions, and chromosomal thromboplasties are not included as they do not have associated CCFs and cannot be used to determine the order of events.

[0249] The sampled set of anomalies is then used to calculate the proportion of tumors that have these anomalies classified as canonical evotypes (step S1010). j The probability that a tumor with p is assigned to a canonical evotype j To perform this calculation, anomaly data is used. A set of sampled anomalies A j ={a1,a2,...,a n}, we consider A j ⊆P k We identified patients who were k denotes the complete set of abnormalities present in patient k. We can then identify which of these were assigned to the canonical evotype. We then define the probability

[0250]

number

[0251] can be calculated, where N(·) denotes the number of tumors subject to the conditions in brackets. Next, we calculate the conditional probability

[0252]

number

[0253] can be calculated. The final step in each inner loop (Monte Carlo loop) is to determine whether further iterations should be performed (step S1012). If further iterations should be performed, the method loops back to step S1002 of randomly selecting the number of anomalies, and steps S1004, S1006, S1008, and S1010 are repeated.

[0254] If no further iterations of the inner loop should be performed, i.e. no further samples are to be considered, the results obtained so far are collated (step S1014). i The set of probabilities for all selected samples of the average number s i =[p1,p2,...,p j ], where p1 is the probability calculated in the first iteration, p2, p3, ... until the jth iteration is completed. The collated set of probabilities is used to obtain a non-parametric density estimate (step S1016), pdf(canonical|x i ) is an anomaly x i is the probability density function of tumors being assigned to a canonical evotype for the average number of iThe values ​​of x are passed to a nonparametric density estimation scheme using a Gaussian kernel with bandwidth 0.025. Because we are estimating the probability density function for a set of probabilities bounded in [0,1], we use reflection techniques to ensure support only over this interval (75).

[0255] In this example, we ran 100,000 samples, so each p(canonical | A j ) for 100,000 values. The next step is to determine whether further iterations of the outer loop i should be performed (step S1018). The outer loop i performs a step of iterating each number of average anomalies, e.g., x i ∈{0, 0.01, 0.02, ..., 10}, i∈1, 2, ..., 1000. If all iterations are not yet completed, the method loops back to step S1000, where the average number of anomalies is updated in the first step S1000, and then the inner loop j is repeated. If no further iterations of the outer loop should be performed, the obtained results are collated (step S1020).

[0256] In summary, loop i iterates with increasing numbers of average anomalies, and loop j runs multiple samples, randomly selecting patients and sampling the sequence of events for possible phylogenetic trees and consistency with the current number of average anomalies. The samples are collated and used to estimate a probability density function for each average number of anomalies.

[0257] Figure 11b shows an example output from the algorithm in Figure 11a, generated using data from 159 samples as previously described. Figure 11b is a surface plot showing the probability density of tumors being assigned to a canonical evotype with respect to the number of aberrations. Since individual evolutionary trajectories involve the stochastic accumulation of multiple genomic aberrations, it is not possible to specify each evolutionary route. However, linking regions with increasing density as the number of aberrations increases can indicate a common mode of evolutionary progress. In other words, we can determine the common mode of evolution by tracking the genetic alterations that are prevalent in tumors at the time of convergence to any evotype in this model. Through this, we can identify paths in the probability density surface plot that correspond to the accumulation of these genetic alterations.

[0258] Initially, the probability density is centered at ~0.78, the fraction of canonical evotype tumors in our sample set. As the number of aberrations increases, the density diverges and accumulates to 1 (corresponding to unambiguous assignment to the canonical evotype) and 0 (alternative evotype). Individual tumors follow trajectories through this probabilistic landscape that depend on the type and order of aberrations, favoring regions of high probability density that are not necessarily contiguous. Examples of such routes (or pathways) are illustrated by the black dashed lines in Fig. 11b. These include Canonical:Rapid, Canonical:Moderate, Canonical:Punctuated, Alternative:Rapid, and Alternative:Incremental. There are also two balanced routes involving LOH in NKX3.1 or IL6ST, or LOH in RB1 and BRCA2. Labels include their likely evotypes, behavioral descriptions, and notable driver genes affected by the aberrations prevalent in regions along the pathway.

[0259] Canonical:Rapid is indicated by early TP53 loss or ERG gene fusion fixation leading to the canonical evotype. Alternatively, loss of regions covering PTEN or CDH1 forces progression to the canonical evotype, and this evolutionary trajectory is referred to as Canonical:Moderate. For the canonical evotype, there are numerous abnormalities that were often the last steps in convergence, notably LOH of 19p13.3-19p13.2, and gain of chromosome 19 and region 22q11.1-22q11.23, and this trajectory is referred to as Canonical:Punctuated.

[0260] When an SPOP mutation first occurs, it confers a high probability (~0.91) of progression to the alternative evotype, which is termed Alternative:Rapid. The other route to the alternative evotype involves the accumulation, in any order, of multiple separate LOH events involving genes such as MAP3K7, CHD1, or EDNRB; this trajectory is termed Alternative:Incremental. LOH of IL6ST or gain of the region 8p23.3-8p22 strongly influences convergence after multiple abnormalities have already accumulated, which is termed Alternative:Abrupt.

[0261] In other words, the model simulations from Fig. 11a can be used to investigate the common evolutionary trajectories involved in the convergence to each evotype (black dashed lines in Fig. 11b). The anomalies that characterize the common evolutionary process can be further investigated, as detailed in Figs. 12a to 13c. In the modeling process, we recorded the order of genetic modifications for each of the trajectories used to calculate the pdf. We calculated the number of trajectories that converged to the canonical or alternative evotype (i.e., p(canonical|A j ) = 0 or 1) and set these by the number of genetic alterations in the locus, i.e., {A1},{A2},...,{A 10}. We then performed a filtering step on each set to remove trajectories that occurred in the set corresponding to fewer genetic alterations, i.e., only trajectories that converged to any evotype with a final genetic alteration for each set were left. We can then identify the location and frequency of occurrence of each genetic alteration in each set. Using this information, we can calculate pdf values ​​for frequent combinations of genetic alterations in turn and use these to create representative paths through the probability density (dashed black lines, shown in Fig. 11b).

[0262] Figure 12a is a 2D surface plot showing the probability density of all canonical evotype tumors being assigned to the canonical evotype with increasing number of aberrations. Figure 12b is a graph showing the percentage of lineages that converged to the canonical evotype at each number of genetic alterations. Figure 12c is a bar graph showing the relative percentage of genetic alterations and the position where they occurred for lineages that converged to the canonical evotype. The bar graph shows the relative percentage for each number of genetic alterations (e.g., 2, 3, ..., 10). For example, for two genetic alterations, the relative percentage of each of the first and second genetic alterations is shown, with the largest percentages shown for TP53 and ERG.

[0263] Figure 13a is a 2D surface plot showing the probability density that all alternative evotype tumors are assigned to the alternative evotype with increasing number of aberrations. Figure 13b is a graph showing the percentage of lineages that converged to the alternative evotype at each number of genetic modifications. Figure 13c is a bar graph showing the relative percentage of genetic modifications and the position where they occurred for lineages that converged to the alternative evotype. The bar graph shows the relative percentage for each number of genetic modifications (e.g., 2, 3, ..., 10). For example, for two genetic modifications, the relative percentage of each of the first and second genetic modifications is shown, with the largest percentage occurring in SPOP.

[0264] Taken together, the findings described above reveal prostate cancer subtypes that arise as a result of distinct trajectories of a stochastic evolutionary process in which different alterations may shift the landscape toward either outcome. The definition of evotypes provides additional context to the relationship between individual abnormalities reported in previous studies. Previously identified co-occurring abnormalities can be associated with specific evotypes. For the canonical evotype, this includes LOH events affecting PTEN and CDH1 (20), or PTEN and TP53 (21). Conversely, CHD1 loss, previously observed in conjunction with SPOP mutations (22, 23), but with LOH affecting MAP3K7 (24) and 2q22 (25), all of which are associated with alternative evotypes.

[0265] The most widely used basis for genomic prostate cancer subtyping is ETS status, where tumors are classified as ETS+ and ETS- depending on the presence or absence of ETS fusions, respectively (7, 8, 10, 11). Figures 14a and 14b illustrate the methods described above and some comparative data for the ETS data. For alternative evotype tumors, 94% were ETS-. Furthermore, alterations such as SPOP mutations and CHD1 LOH characteristic of this evotype have previously been associated with the ETS- subtype (10, 26). In contrast, for tumors of the canonical evotype, there is a relatively even balance of ETS- and ETS+ tumors.

[0266] Figure 14a shows the aberrations present in canonical evotype tumors when divided into ETS- (n=42 or 44%) and ETS+ (n=83 or 66%). Continuous variables were converted to binary by setting values ​​equal to or greater than the median to 1 and values ​​less than the median to 0. Samples were ordered by hierarchical clustering using Hamming distance average. No aberrations were significantly associated with either ETS group (Q>0.05, Fisher's exact test). Figure 14b shows the Kaplan-Meier plot of ETS+ and ETS- tumors classified into canonical evotypes. P-values ​​(0.909) and hazard ratios (1.06 (0.413 to 2.7)) were calculated using the log-rank method, and HRs were estimated in the range from the 5th to the 96th percentile. The end point in Figure 14b is the time to biochemical recurrence. As shown in these figures, there were no significant differences in risk or prevalence of genomic features between ETS+ and ETS− tumors of the canonical evotype consistent with their definition as distinct disease types.

[0267] Application of Evotype to the classification of new tumors Having established the existence of evolutionary pathotypes, we can directly classify metaclusters and even evotypes from the feature set using classification methods such as, but not limited to, neural networks, random forests, and boosted decision trees. Figure 15a illustrates a possible method for classifying tumors. In the first step, a dataset is received (step S200). Figure 15a shows that three different classifications can be performed based on the received dataset. These classifications are a clustering classification based on clustering, an ARBS classification based on the proximity of DNA breakpoints to the androgen receptor binding site (ARBS), and an ordering classification based on the ordering of events. These classifications can be performed in parallel or sequentially. In a preferred arrangement, all three classifications are performed, but the overall classification of the tumor can be based on one, two, or three of the classifications (which may be referred to as intermediate classifications). The received data must be relevant to the classification method used. For example, in a cluster-based classification, gene sequencing can be used to extract relevant information.

[0268] When using the metacluster classification shown in the first branch, the trained neural network is used to process raw inputs from new samples to generate feature representations that can be used to assign them to one of the metaclusters. Alternatively, we can develop a simple ML classifier that uses the raw inputs (or a subset of them) to directly classify by metaclusters. For example, the SHapley Additive explanation (SHAP) value can be used to quantify the relative importance of a feature when performing classification using the gradient boosting decision tree method XGBoost. A value of zero indicates that the feature is not necessary to perform the classification. Figure 15b illustrates the SHAP value for each feature when classifying a sample as belonging to metacluster A, and Figure 15c illustrates the SHAP value for each feature when classifying a sample as belonging to metacluster B. When generating these figures, the ARBS score has been omitted as a feature for consistency with the other figures below.

[0269] For each of Figures 15b and 15c, each feature is ranked by its feature value, and not surprisingly, there are similar rankings in each figure. Each feature with a high positive ranking for metacluster A has a high negative ranking for metacluster B, and vice versa, with positive values ​​indicating that the feature is strongly suggestive of belonging to a particular metacluster:

[0270] [Table 17]

[0271] [Table 18]

[0272] In a first step, samples can optionally be described in terms of genomic features (S204), such as those identified above. Tumors can then be classified into specific clusters (S206). We are able to achieve 95.60% accuracy in distinguishing metacluster 1 (MC-A) from metaclusters 2 (MC-B1) and 3 (MC-B2).

[0273] Referring again to FIG. 15a, the next classification illustrated is ARBS classification, where the first step is to obtain the location of DNA breakpoints with respect to the androgen receptor binding site (ARBS) from the input data (step S214). To classify the tumor, the closeness of the obtained location to the ARBS was compared to the closeness of the location in the expected distribution to the ARBS, thereby determining the ARBS score. The expected distribution of breakpoints may be referred to as the baseline distribution. The class (enriched, depleted, or indeterminate) was determined based on whether the tumor exhibited more proximal breakpoints than expected (enriched), fewer proximal breakpoints than expected (depleted), or no statistically significant difference (indeterminate) (step S214). The baseline distribution may be defined as the distribution that would be expected if the DNA breakpoints were uniformly distributed over the genome. Alternatively, another baseline distribution may be used.

[0274] As an example, the baseline (or expected) distribution may include permutation data that may be generated by simulating 1000 data sets in which the location of DNA breakpoints is permuted to new locations in the genome with a uniform distribution. The distance to the nearest ARBS was calculated for each simulated breakpoint in each data set. Similarly, the distance to the nearest ARBS was calculated for each observed or obtained breakpoint. A double-stranded DNA break may be considered relatively proximal to an ARBS when the break is less than a threshold number of base pairs (e.g., 20,000 bps) from the ARBS. The ARBS score may be calculated by normalizing the number of relatively proximal DNA breaks by the number of proximal breakpoints expected by chance. Thus, the percentage of breakpoint locations that are relatively proximal may be calculated for both the observed data (B_obs) and the permutation data (B_perm). If the observed proportion of relatively proximal breakpoints (B_obs) is higher than the upper threshold (e.g., the 97.5% percentile) of the proportion of relatively proximal breakpoints in the permutation data (B_perm), i.e., B_obs>P _ 97.5% (B_perm). In other words, a tumor is classified as enriched if the ARBS score is higher than an upper threshold (e.g., 97.5%). A tumor is classified as enriched if the observed proportion of relatively proximal breakpoints (B_obs) is lower than a lower threshold (e.g., the 2.5% percentile) of the proportion of relatively proximal breakpoints in the permutation data (B_perm). obs <P 2.5% (B perm )), the tumor is classified as depleted. In other words, if the ARBS score is higher than a lower threshold (e.g., 2.5%), the tumor is classified as depleted. Otherwise, the difference is not significant and the tumor is classified as indeterminate.

[0275] When using the ARBS classification shown in the second branch, the feature representation is used along with the ARBS score itself, so that each sample can be assigned to one of two classifications: enriched or depleted. The SHapley Additive explanation (SHAP) value can be used to quantify the relative importance of features when performing classification using the gradient boosting decision tree method XGBoost. Figure 15d illustrates the SHAP value for each feature when classifying a sample as belonging to a depleted tumor, and Figure 15e illustrates the SHAP value for each feature when classifying a sample as belonging to an enriched tumor. When generating these figures, the ARBS score has been omitted as a feature for consistency with the other figures.

[0276] For each of Figures 15d and 15e, each feature is ranked by its feature value, and not surprisingly, there are similar rankings in each figure. Each feature with a high positive ranking for depleted tumors has a high negative ranking for enriched tumors and vice versa, with positive values ​​indicating that the feature is strongly suggestive of belonging to a particular type of tumor:

[0277] [Table 19]

[0278] [Table 20]

[0279] The next classification illustrated is an ordered classification. This can be done by inferring the order of genetic modifications (step S224). The order of genetic modifications can be inferred by performing bulk cell sequencing and determining the percentage of cells that contain each genetic modification. It is determined that the abnormality present in a higher percentage of cells occurred before the abnormality present in a lower percentage of cells. For example, if we estimate that CHD1 LOH appears in 90% of cancer cells and PTEN LOH appears in 40% of cancer cells, there must be cells that contain CHD1 LOH alone and CHD1 LOH and PTEN LOH. Therefore, CHD1 LOH occurs first. The abnormalities can be ranked in order of percentage.

[0280] The tumor may then be classified based on the determined order of anomalies (step S226). Similar to Figures 15b to 15e, the SHapley Additive explanation (SHAP) value may be used to quantify the relative importance of features when performing classification using the gradient boosting decision tree method XGBoost. Figure 15f illustrates the SHAP value for each feature when classifying samples as belonging to ordering I, and Figure 15g illustrates the SHAP value for each feature when classifying samples as belonging to ordering II. When generating these figures, the ARBS score has been omitted as a feature for consistency with the other figures.

[0281] For each of Figures 15f and 15g, each feature is ranked by its feature value, and not surprisingly, there are similar rankings in each figure. Each feature with a high positive ranking for ordering I has a high negative ranking for ordering II, and vice versa, with positive values ​​indicating that the feature strongly suggests belonging to a particular ordering:

[0282] [Table 21]

[0283] [Table 22]

[0284] As suggested in the figure and table above, genomic abnormalities that may indicate ordering classification include loss of heterozygosity, SPOP mutations, and some or all of ETS fusions in one or more of the following regions: 5q15-5q23.1 (spanning CHD1), 6q14.1-6q22.32 (MAP3K7, ZNF292), 8p (NKX3.1), 10q23.1-10q25.1 (PTEN), 13q12.3-13q21.1 (RB1, BRCA2) and 13q21.1-13q33.1 (EDNRB), 16q12.1-16q24.1 (CDH1) and 17p (TP53). These are summarized in Table B below.

[0285] [Table 23]

[0286] The next step (step S230) may be to combine one or more of the clustering, ARBS, and ordering classifications to provide an overall classification for the tumor. A tumor classified as an alternative evotype indicates a poor prognosis. Each of the clustering classification as metacluster MC-A, the ARBS classification of depletion, and the ordering classification of ordering II indicates an overall classification as an alternative evotype. If all three intermediate classifications are used, an overall classification as an alternative evotype is provided when the tumor has at least two intermediate classifications selected from the classification as metacluster MC-A, the ARBS classification of depletion, and the ordering classification of ordering II. Similarly, each of the clustering classification as metacluster MC-B1 or B2, the ARBS classification of enrichment or indeterminate, and the ordering classification of ordering I indicates an overall classification as a canonical evotype. A tumor may be assigned to a canonical evotype based on a similar majority approach when at least two of the intermediate classifications indicate a canonical evotype.

[0287] Instead of proceeding with each of the classifications separately, it is also possible to classify tumors directly as either canonical or alternative evotypes based on the presence of a combination of genomic abnormalities. This can be used in combination with one or more of the classifications. Alternatively, as shown by the dotted line, the method proceeds directly from receiving the dataset in step S200 to identifying genetic abnormalities (step S232). Similar to metacluster classification, a trained neural network is used to process raw inputs from new samples to generate feature representations that can be used to assign to one of the evotypes. Alternatively, we can develop a simple ML classifier that uses the raw inputs (or a subset thereof) to directly classify by evotype as shown in step S234. For example, the SHapley Additive explanation (SHAP) value can be used to quantify the relative importance of features when performing classification using the gradient boosting decision tree method XGBoost. A value of zero indicates that the feature is not necessary to perform the classification. Figure 15h illustrates the SHAP value for directly classifying evotypes. Comparing these features with those shown in Figure 10e, there is a significant overlap with the highest ranking SHAP values. The overlapping features are listed in the table below according to their ranks shown in Figure 15h. For the ARBS score, as shown in Figure 15h, the SHAP value is significantly higher for this score, and therefore this is likely to be the most useful feature. As explained above, this can be used to indicate whether a tumor is a canonical or alternative evotype by considering a threshold value.

[0288] [Table 24]

[0289] Tumors can then be classified into specific evotypes using these features. It is likely to focus on methods and / or kits that target combinations of specific regions as described above. More general genomic tests, such as tests to determine whether chromosilysis or PGA is present, may be omitted from the kit and / or subject classification method to provide a quicker and more convenient test / method. We can achieve an accuracy of 94.97% when directly classifying canonical and alternative evotypes. The classification is then output in step S236, optionally with an associated probability that the assignment to the classification is correct.

[0290] As can be seen, there is an overlap between the features considered for each sub-classification (ARBS, clustering, and ordering) as well as direct classification. These are compared in the table below and ranked using the rankings in Figure 15h. Anomalies are marked with a Y in the table if they have a positive value in the corresponding figure and their SHAP score is within 99% of the cumulative total.

[0291] [Table 25]

[0292] As shown in Table 1 above, genomic abnormalities were present in all three subclassifications, and these were LOH at 6q12-6q22.32 (MAP3K7, ZNF292), LOH at 5q15-5q23.1 (CHD1), LOH at 2q14.3-2q23.3, chromosomal thrombocytopenia, and LOH at 1q42.12-1q42.13. Thus, the presence of a combination of some or all of these features can also be used to classify a subject into a first prognostic group, in particular where the combination includes at least two of the highest ranked features targeting specific regions in the genome, for example at least the top four, i.e., LOH at 6q12-6q22.32 (MAP3K7, ZNF292), LOH at 5q15-5q23.1 (CHD1), LOH at 2q14.3-2q23.3, and LOH at 1q42.12-1q42.13, more particularly at least the top two, i.e., LOH at 6q12-6q22.32 (MAP3K7, ZNF292) and LOH at 5q15-5q23.1 (CHD1). Similarly, genomic abnormalities exist in at least two sub-classifications, which are: kataegis, LOH at 5q11.1-5q14.1 (IL6ST, PDE4D), gain of the entire chromosome 7, gain of 8p23.3-8p22, LOH:18q, LOH at 12p12.32-12p12.3, LOH at 13q12.3-13q21.1, SPOP, gain of 8q (MYC).Therefore, the presence of a combination of some or all of these features can also be used to classify subjects into a first prognostic group, in particular, the combination includes at least the highest-ranked features that target specific regions, which are LOH at 5q11.1-5q14.1 (IL6ST, PDE4D) and gain of the entire chromosome 7.There can also be three sub-classifications and a combination of two sub-classifications.For example, the presence of a combination including at least the highest ranked features targeting particular regions, e.g., LOH at 5q11.1-5q14.1 (IL6ST, PDE4D), LOH at 6q12-6q22.32 (MAP3K7, ZNF292), and LOH at 5q15-5q23.1 (CHD1), can be used to classify into a first prognostic group.

[0293] [Table 26]

[0294] As shown in Table 2 (Table 26) above, genomic abnormalities exist in all three subclassifications, which are ETS gene fusions and LOH at 17p. Therefore, the presence of some or all of these features in combination can also be used to classify subjects into a second prognostic group, and in particular, the combination includes at least a feature that targets a specific region, such as LOH at 17p. Similarly, genomic abnormalities exist in at least subclassifications, which are interchromosomal / intrachromosomal breakpoint ratio and LOH at 21q22.2-21q22.3 (ERG). Thus, the presence of at least these features in combination can also be used to classify subjects into a second prognostic group. There can also be a combination of three subclassifications and two subclassifications. For example, the presence of a combination that includes at least a feature that targets a specific region, such as LOH at 17p and LOH at 21q22.2-21q22.3 (ERG), can also be used to classify tumors into a second prognostic group.

[0295] Instead of using the combinations described above, combinations based on the ranking of the proportion of tumors with the features shown in Fig. 10e can also be used.For example, a combination that includes at least two of the features that target specific regions, such as LOH at 6q12-6q22.32, LOH at 13q21.1-13q33.1, and LOH at 13q12.3-13q21.1, can also be used to classify subjects as belonging to a first prognosis group.For example, a combination that includes at least two of the features that target specific regions, such as LOH at 17p and LOH at 21q22.2-21q22.3 (ERG), can also be used to classify subjects as belonging to a second prognosis group.It will be understood that these selections are included only as examples, and the top three, four, five, or more features can also be included.

[0296] In the various tables and descriptions above, there are acronyms for genes, which are listed below along with the full gene name.

[0297] [Table 27A]

[0298] [Table 27B]

[0299] Figure 16 illustrates that evotypes can also be classified using other techniques. For example, we have RNA-seq from tumor and adjacent normal tissues for 136 of the 159 samples used to derive evotypes. Performing discriminative gene expression analysis using the EdgeR package reveals that there are 588 genes that are significantly differentially expressed (adjusted P-value < 0.05) between the canonical and alternative evotypes. This set can potentially be used as a basis for classification by evotype. By performing classification with XGBoost, we can see that we obtain a classification accuracy of 84.56%. Calculating the SHAP value for this classifier, we find 77 variables with non-zero SHAP values, as shown in Figure 16.

[0300] Since this is still a fairly large number of parameters to explore in the XGBoost algorithm, we attempt to further reduce the number of inputs and optimize the classification by finding the set of transcripts with the highest SHAP value that maximizes the classification accuracy. Through this method, we find that we can obtain a maximum classification accuracy of 91.91% when classifying using the top 18 transcripts. These are listed in the table below. Table of features from RNMA expression that can be used for classification

[0301] [Table 28]

[0302] Therefore, we conclude that information from RNA expression can directly classify evotypes. Furthermore, using the full set of 77 transcripts, we obtain an accuracy of 94.12% in classifying tumor and benign samples using XGBoost.

[0303] FIG. 17 is a schematic diagram of a relevant system for carrying out computer-implemented aspects of the methods described above (both discovery and classification). The system includes a computing device 10, which may be a portable handheld device that the clinician can carry from patient to patient, and an app for performing the prediction may be loaded onto the device. The computing device 10 comprises standard components such as a processing unit or processor 20, a user interface unit 22 for allowing a user to input information, and a memory 24. The user interface may be a display 24 for displaying information or, alternatively, displaying information to the user, e.g., treatment suggestions as described above. There may also be a communication module 28 to communicate with other devices and / or access the cloud, e.g., for processing data as described below.

[0304] The computing device 10 also has a discrimination score module 30 for calculating a discrimination score, a clustering module 32 for determining a clustering classification, a DNA breakpoint analysis module 34 for analyzing the placement of breakpoints within the base sequence for use in determining an ARBS classification, and an ordering module 36 for determining an ordering classification as described above. Each of the modules may be stored in the memory 24 or in a separate storage device (not shown) of the device. The modules may also be located remotely from the computing device 10, for example in the cloud.

[0305] The system shown in this schematic diagram may be fabricated in part or in whole using dedicated hardware. Terms such as "module" or "unit" as used herein may include, but are not limited to, hardware devices that perform specific tasks or provide related functionality, such as circuits in the form of discrete or integrated components, field programmable gate arrays (FPGAs) or application specific integrated circuits (ASICs). In some embodiments, the described elements may be configured to reside on tangible, persistent, addressable storage media and configured to execute on one or more processors. These functional elements may include, in some embodiments, components such as, for example, software components, object-oriented software components, class components and task components, processes, functions, attributes, procedures, subroutines, segments of program code, drivers, firmware, microcode, circuits, data, databases, data structures, tables, arrays, and variables. Although the exemplary embodiments are described with reference to the components described herein, such functional elements may be combined into fewer elements or separated into additional elements.

[0306] summary As described above, a comprehensive analysis of genomic measurements from 159 prostate cancer patients has been performed using three statistical and machine learning methods. In this analysis, two distinct forms of prostate cancer evolutionary types, referred to herein as "evotypes", can be characterized by various characteristics. First, evotypes can be characterized by the arrangement of double-stranded DNA breaks with respect to androgen receptor binding sites (ARBS classification as described above). Second, evotypes can be characterized by some genetic abnormalities and some combinations of genetic abnormalities (e.g., using clustering classification or ordering classification as described above). Evotypes can be characterized by the arrangement of DNA double-stranded breaks and combinations of genetic abnormalities.

[0307] Stratification by evotype may have epidemiological implications. For example, non-Caucasian groups may have a higher incidence of many alternative evotype aberrations (27-29) and therefore a higher predisposition to this phenotype. Conversely, cancers occurring in younger patients have been reported to have enrichment for ARBS proximal breakpoints (17) and develop via an evolutionary progression similar to the canonical evotype (14, 17). It may also be possible to tailor therapeutic strategies to each evotype. In particular, cancers with alternative evotype aberrations have been shown to be sensitive to ionizing radiation (22) and to have a good response to treatment with PARP inhibitors (30) and androgen deprivation (23). Thus, our model of prostate cancer evolutionary phenotypes provides a conceptual framework that unifies the results of many previous studies and has far-reaching implications for our understanding of disease progression, prognosis, and treatment.

[0308] Unless otherwise defined herein, scientific and technical terms used in connection with this disclosure shall have the meanings commonly understood by those skilled in the art. The foregoing disclosure provides a general description of the subject matter encompassed within the scope of the present invention, including how to make and use the invention, as well as the best mode thereof, but the following examples are provided to further enable those skilled in the art to practice the invention and provide a complete written description thereof. However, those skilled in the art will understand that the specific content of these examples should not be read as limiting the invention, the scope of which should be understood from the claims appended to this disclosure and their equivalents. Various further aspects and embodiments of the present invention will be apparent to those skilled in the art in light of the present disclosure.

[0309] All of the features disclosed in this specification (including the accompanying claims, abstract, and drawings), and / or all of the steps of the methods or processes so disclosed, may be combined in any combination, except combinations in which at least some of such features and / or steps are mutually exclusive.

[0310] Each feature disclosed in this specification (including the accompanying claims, abstract, and drawings), unless expressly stated otherwise, may be replaced by alternative features serving the same, equivalent, or similar purpose. Thus, unless expressly stated otherwise, each feature disclosed is only an example of a generic series of equivalent or similar functions.

[0311] All documents and references that list gene / protein accession numbers mentioned herein are incorporated herein by reference in their entirety.As used herein, "and / or" is deemed to be a specific disclosure of each of the two specified features or components, regardless of the presence or absence of the other.For example, "A and / or B" is interpreted as a specific disclosure of each of (i) A, (ii) B, and (iii) A and B, as if each were individually described herein.Unless otherwise indicated by the context, the description and definition of features described above are not limited to a particular aspect or embodiment of the present invention, but are equally applicable to all aspects and embodiments described. [Explanation of symbols]

[0312] 10. Computing Devices 20 Processing Unit or Processor 22 User Interface Unit 24 Memory 24 Display 28 Communication Module 30 Discrimination Score Module 32 Clustering Module 34 DNA Breakpoint Analysis Module 36 Sequencing Module

Claims

1. 1. A method of stratifying a subject into one of two prognostic groups, said method comprising: analyzing a biological sample obtained from said subject with cancer or metastatic disease using DNA sequencing and / or RNA sequencing; determining the location of double-stranded DNA breakpoints relative to the androgen receptor binding site (ARBS) in the biological sample; obtaining an ARBS score for the sample by comparing the closeness of the determined configuration to the ARBS with the closeness of a baseline distribution of double-stranded DNA breakpoints to the ARBS; The ARBS score determines whether the subject a first prognostic group when the ARBS score indicates that the determined location is proximal to an androgen receptor binding site less frequently than expected; and classifying the determined location into a second prognostic group when the ARBS score indicates that the determined location is proximal to an androgen receptor binding site more frequently than expected.

2. 2. The method of claim 1, further comprising classifying the subject into the second prognostic group when the ARBS score indicates that there is no statistically significant difference between the proximity of the determined location to the androgen receptor binding site and the proximity of the predicted placement of the breakpoint to the androgen receptor binding site.

3. 10. The method of claim 1, further comprising defining a baseline distribution of breakpoint placements by randomly shuffling observed breakpoints in the sample data.

4. calculating the ARBS score determining the percentage of the determined alignments that are less than a threshold number of base pairs from an androgen receptor binding site; obtaining a percentage of the breakpoint locations within the baseline distribution that are less than a threshold number of base pairs from an androgen receptor binding site; 2. The method of claim 1, further comprising the step of normalizing the determined proportion by the obtained proportion to obtain the ARBS score, and determining whether the determined location is more or less frequently proximal to an androgen receptor binding site than expected.

5. 5. The method of claim 4, further comprising classifying the subject into the second prognostic group when the ARBS score is greater than an upper threshold and classifying the subject into the first prognostic group when the ARBS score is less than a lower threshold.

6. identifying additional genomic abnormalities present in the sample; Furthermore, the object classified into a first prognostic group based on the presence of one or more genomic abnormalities selected from Table 1; and and classifying the patient into a second prognostic group based on the presence of one or more genomic abnormalities selected from Table 2.

7. classifying the subject into the first prognostic group based on the ARBS score and the presence of a combination of genomic abnormalities including loss of heterozygosity in at least the regions 6q12-6q22.32 (MAP3K7, ZNF292) and 5q15-5q23.1 (CHD1); and classifying the subject into the second prognostic group based on the ARBS score and the presence of a combination of genomic abnormalities comprising loss of heterozygosity in regions 16q12.1-16q24.3 and 17p.

8. Regions 1q42.12-1q42.13, 2q14.3-2q23.3, 5q11.1-5q14.1 (IL6ST, using a clustering classification to classify the subject into the first prognostic group based on the presence of one or more genomic abnormalities selected from a set of genomic abnormalities including: loss of heterozygosity in 12p12.32-12p12.3 and 18q, gain in the entire chromosome 7 and region 8q, kataegis, SPOP mutations, and percentage of genomic alterations (clonal component), more specifically based on the presence of a combination of genomic abnormalities including loss of heterozygosity in regions 2q14.3-2q23.3, 5q11.1-5q14.1 (IL6ST, PDE4D), 5q15-5q23.1 (spanning CHD1), kataegis, and percentage of genomic alterations (clonal component); and ETS gene fusions, the ratio of interchromosomal / intrachromosomal breakpoints, the percentage of genomic alterations (subclonal component), losses of heterozygosity in regions 17p(TP53), 16q12.1-16q24.3 and 22q11.21-22q11.22, and gains in regions 9q12.9-9q21.11 and the entire chromosome 19, more specifically based on the presence of a combination of genomic abnormalities including ETS gene fusions, the ratio of interchromosomal / intrachromosomal breakpoints, the percentage of genomic alterations (subclonal component), losses of heterozygosity in regions 17p(TP53) and 16q12.1-16q24.

3.

9. 1. A method of stratifying a subject into one of two prognostic groups, said method comprising: analyzing a biological sample obtained from a subject with cancer or metastatic disease using DNA sequencing and / or RNA sequencing; identifying genomic abnormalities in the biological sample; By clustering classification, the object The first prognostic group is based on the presence of one or more genomic abnormalities selected from a set of genomic abnormalities including loss of heterozygosity in regions 1q42.12-1q42.13, 2q14.3-2q23.3, 5q11.1-5q14.1 (IL6ST, PDE4D), 5q15-5q23.1 (spanning CHD1), 6q14.1-6q22.32 (MAP3K7, ZNF292), 12p12.32-12p12.3 and 18q, gain of heterozygosity in the entire chromosome 7 and region 8q, kataegis, SPOP mutations, and genomic alteration percentage (clonal component). and classifying the patient into a second prognostic group based on the presence of one or more genomic abnormalities selected from a set of genomic abnormalities including ETS gene fusions, interchromosomal / intrachromosomal breakpoint ratios, genomic alteration percentage (subclonal component), loss of heterozygosity in regions 17p(TP53), 16q12.1-16q24.3 and 22q11.21-22q11.22, and gains in regions 9q12.9-9q21.11 and the entirety of chromosome 19.

10. analyzing the biological sample obtained from the subject with cancer or metastatic disease using bulk cell sequencing; determining the proportion of cells in the biological sample that have one or more genomic abnormalities; identifying the order in which the genomic abnormalities occurred by determining that the genomic abnormalities present in a large proportion of cells occurred before the genomic abnormalities present in a small proportion of cells; and classifying the subject into one of the first and second prognostic groups using an ordered classification based on the identified order; The genomic abnormalities were located in the regions 5q15-5q23.1 (spanning CHD1), 6q14.1-6q22.32 (MAP3K7, ZNF292), 8p (NKX3.1), 10q23.1-10q25.1 (PTEN), 13q12.3-13q21.1 (RB1, BRCA2), 13q21.1-13q33.1 (EDNRB), 16q12.1-16q24.1 (CDH1), and 17p (TP 53), and at least one or more of loss of heterozygosity in one or more of 21q22.2-21q22.3, gain in one or more of regions 8p23.3-8p22 and 9q12.9-9q21.11, genomic alteration percentage (subclonal component), ratio of intrachromosomal / interchromosomal linked structural variants, and ETS fusions.

11. A method for stratifying a subject into one of two prognostic groups, said method comprising: analyzing a biological sample obtained from a subject suffering from cancer or metastatic disease using bulk cell sequencing; determining the proportion of cells in the biological sample that have one or more genomic abnormalities; identifying the order in which the genomic abnormalities occurred by determining that the genomic abnormalities present in a large proportion of cells occurred before the genomic abnormalities present in a small proportion of cells; and classifying the subject into one of the first and second prognostic groups by an ordered classification based on the identified order; The genomic abnormalities were located in the regions 5q15-5q23.1 (spanning CHD1), 6q14.1-6q22.32 (MAP3K7, ZNF292), 8p (NKX3.1), 10q23.1-10q25.1 (PTEN), 13q12.3-13q21.1 (RB1, BRCA2), 13q21.1-13q33.1 (EDNRB), 16q12.1-16q24.1 (CDH1), and 17q15.1-17q26.1 (SEQ ID NO: 1). p(TP53), and loss of heterozygosity in one or more of 21q22.2-21q22.3, gain in one or more of regions 8p23.3-8p22 and 9q12.9-9q21.11, genomic alteration percentage (subclonal component), ratio of intrachromosomal / interchromosomal linked structural variants, and at least one or more of ETS fusions.

12. the genomic abnormalities include at least one or more of: loss of heterozygosity in one or more of the regions 5q15-5q23.1 (spanning CHD1), 6q14.1-6q22.32 (MAP3K7, ZNF292), 13q12.3-13q21.1 (RB1, BRCA2), 13q21.1-13q33.1 (EDNRB); gain in one or more of the regions 8p23.3-8p22 and 9q12.9-9q21.11; percentage of genomic alterations (subclonal component); and ratio of intra- / inter-chromosomal linked structural variants.

12. The method of claim 10 or claim 11, wherein the method sometimes comprises classifying the subject into the first prognostic group, more specifically classifying the subject into the first prognostic group based on the presence of a combination of genomic abnormalities including loss of heterozygosity and a genomic alteration percentage (subclonal component) in the regions 5q15-5q23.1 (spanning CHD1), 6q14.1-6q22.32 (MAP3K7, ZNF292), 13q12.3-13q21.1 (RB1, BRCA2), 13q21.1-13q33.1 (EDNRB).

13. 12. The method of claim 10 or claim 11, comprising classifying the subject into the second prognostic group when the genomic abnormalities include at least one or more of loss of heterozygosity in one or more of regions 8p (NKX3.1), 10q23.1-10q25.1 (PTEN), 16q12.1-16q24.1 (CDH1), 17p (TP53) and 21q22.2-21q22.3, and ETS fusions, more particularly classifying the subject into the second prognostic group based on the presence of a combination of genomic abnormalities including loss of heterozygosity in regions 8p (NKX3.1), 10q23.1-10q25.1 (PTEN), 16q12.1-16q24.1 (CDH1), 17p (TP53) and 21q22.2-21q22.

3.

14. The method described in claim 10, wherein when the subject is classified into the first prognostic group by at least two of the ARBS, clustering and ordering classification, the subject has an overall classification as the first prognostic group, and when the subject is classified into the second prognostic group by at least two of the ARBS, clustering and ordering classification, the subject has an overall classification as the second prognostic group.

15. identifying additional genomic abnormalities present in the sample; Loss of heterozygosity in one or more of the regions 2q14.3-2q23.3, 5q11.1-5q14.1 (IL6ST, PDE4D), 5q15-5q23, 6q12-6q22.32 (MAP3K7, ZNF292), 18q, gain of heterozygosity in one or more of the regions 3q21.2-3q29, the entire chromosome 7, 8p23.3-8p22, 8q, 9q12.9-9q21.11, and kataegis. classifying the subject into the first prognostic group based on the presence of one or more genomic abnormalities, more specifically based on the presence of a combination of genomic abnormalities including loss of heterozygosity in one or more of the regions 2q14.3-2q23.3, 6q12-6q22.32 (MAP3K7, ZNF292), 18q, and gain of heterozygosity in one or more of the following regions: the entire chromosome 7, and 8q; classifying the subject into the second prognostic group based on the presence of one or more genomic abnormalities selected from the group comprising: loss of heterozygosity in one or more of the regions 10q23.1-10q25, 16q12.1-16q24.3, 17p; gain of heterozygosity in one or more of the following regions: entire chromosome 19; ratio of intrachromosomal / interchromosomal linked structural variants; ETS; percentage of genomic alterations (subclonal component); and percentage of genomic alterations (clonal component); more specifically, based on the presence of a combination of genomic abnormalities comprising ratio of intrachromosomal / interchromosomal linked structural variants; loss of heterozygosity in one or more of the regions 10q23.1-10q25, 17p; ETS; and percentage of genomic alterations (subclonal component); 12. The method of claim 1, 9, 10, or 11, comprising:

16. 1. A method of stratifying a subject into one of two prognostic groups, said method comprising: analyzing a biological sample obtained from a subject with cancer or metastatic disease using DNA sequencing and / or RNA sequencing; identifying genomic abnormalities in the biological sample; The object is classified into a first prognostic group based on the presence of one or more genomic abnormalities selected from Table 1; and and classifying the patient into a second prognostic group based on the presence of one or more genomic abnormalities selected from Table 2.

17. 17. The method of claim 16, comprising: classifying the subject into the first prognostic group based on the ARBS score and the presence of a combination of genomic abnormalities comprising loss of heterozygosity in at least the regions 6q12-6q22.32 (MAP3K7, ZNF292) and 5q15-5q23.1 (CHD1); and classifying the subject into the second prognostic group based on the ARBS score and the presence of a combination of genomic abnormalities comprising loss of heterozygosity in the regions 16q12.1-16q24.3 and 17p.

18. 17. A kit for use in the method of claim 1, 9, 10, 11 or 16, comprising reagents for whole genome sequencing and probes for the detection of DNA double strand breaks used to calculate the ARBS score.

19. 17. The method of claim 1, 9, 10, 11 or 16, wherein subjects stratified into the first prognostic group are identified for treatment selected from one or more of external beam radiation, brachytherapy, radical prostatectomy, hormone therapy, and / or chemotherapy, and subjects stratified into the second prognostic group are selected for patient surveillance.

20. 17. The method of claim 1, 9, 10, 11 or 16, wherein the subject has prostate cancer.

Images