The Germline and Somatic Origins of Prostate Cancer Heterogeneity

Abstract

This study uncovered 223 recurrently mutated driver regions using the largest cohort of prostate tumors to date. It reveals associations between germline SNPs, somatic drivers, and tumor aggression, offering significant insights into how prostate tumor evolution is shaped by germline factors and the timing of somatic mutations.

Figure 1.

Mutation rates of prostate tumors. A, Schematic roadmap of the key analyses conducted in this study, offering insights into the genomic origins of the ISUP GG. B, Distributions of somatic mutation frequency for SNVs, indels, mtSNVs, and GRs across 666 tumor–reference WGS pairs. INV, inversions; CTX, interchromosomal translocations. C, Distributions of somatic clonal and subclonal CNA frequency (number of base-pairs) for losses and gains. D, Clonal and subclonal CNA landscape of localized prostate cancer. Heatmaps represent CNA profiles for the cohort split by CNA occurring clonally or subclonally. Columns represent genes and rows represent patients, grouped by subtype, then sorted by ISUP GG, and clustered within the GG. The top two panels show clonal–subclonal differences in CNA frequency (FTrunk – FBranch) of the dominant CNA type and statistical significance −log₁₀(Q value) calculated using Pearson’s χ² test with FDR adjustment. To avoid confounding by subclonal whole-genome duplication, patients with subclonal PGA >80% (9/664) were excluded from this statistical analysis. (A, Created with BioRender.com.)

Figure 2.

Somatic driver mutations in localized prostate cancer. Driver mutation discovery in 666 localized prostate tumors. The top barplot shows the distribution of the number of drivers in patients; the covariates on the left show the region type and statistical significance from ActiveDriverWGS and GISTIC. The top heatmap shows drivers found in this study (rows) for each patient (columns). The bottom heatmap shows drivers found in an exome meta-analysis (14). Heatmaps are colored by mutation type. Right barplot shows the number of drivers per patient. Bottom covariate bars show clinical features of patients. Gene labels on the left are for rows identified by tick-marks on the axis. Schematics above each set of legends indicate the panels to which they apply in yellow. For example, the top set of legends apply to the regional covariate heatmaps on the left side of the plot.

Figure 3.

Functional characterization of driver mutations. A, Network diagrams represent multimodal pathway enrichment analysis of driver genes. Mutation types (i.e., the type of driver analysis) are indicated by shading of circles. Circle size represents the number of patients. Heatmaps show mutations in the cohort that affect genes contributing to two exemplar pathways dysregulated by multiple mutation types. The apoptosis pathway (GO:0042771 – intrinsic apoptotic signaling pathway in response to DNA damage by p53 class mediator), was identified to be significant only in the context of integrating statistical significance from all GR, and SNV, and indel driver analyses in coding and regulatory elements. The embryonic development pathway (GO:0043009 chordate embryonic development) was identified to be independently significant using SNVs and indels in regulatory elements, in coding elements, and in the Armenia and colleagues (14) dataset. Bottom covariates on the heatmaps show CNAs in pathway genes identified in driver CNA peaks (analyzed using GISTIC2). B, Associations between driver events and SBS signatures. Dot size and dot colors indicate median difference of signature activity. Background shading shows Q values from the Wilcoxon rank-sum test with FDR adjustment. Drivers and SBSs were ordered using hierarchical clustering. C, A summary of associations between driver events and mRNA abundance of driver genes. Dot size and colors indicate median difference of mRNA abundance. Background shading shows Q values from the Wilcoxon rank-sum test with FDR adjustment. D, Consensus clustering of 3,318 dysregulated mRNAs associated with driver mutations. Colors in the heatmap indicate median difference of mRNA abundance between patients with and without a specific driver mutation. Driver mutation type is on the left using colors from A. The right barplot shows the number of transcripts significantly associated with each driver mutation. E, Pathway enrichment analysis on the four mRNA subtypes from D. Clusters of biologically similar pathways are labeled and outlined for each subtype. The size of the pathway is indicative of the number of enriched genes. For (B and C) CNA drivers in patients with subclonal PGA >80% were excluded, and only driver events significantly associated with either SBS signatures or mRNA abundances are shown. (A, Created with BioRender.com.)

Figure 4.

Mutational subtypes of localized prostate cancer. A, Mutation densities (rows) differ by ETS fusion and NKX3-1 CNA status (columns). Dot size and color gives effect-size as a Z-score, scaled to ETS-negative, NKX3-1–neutral patients. The barplot on the right shows the FDR-adjusted P values from nonparametric Kruskal–Wallis tests. B, Comparison of log₁₀-transformed SNV mutation rate for patients divided by ETS fusion and NKX3-1 CNA status. P value is from a nonparametric Kruskal–Wallis test. C, Using a generalized linear model, eight driver mutations were identified whose frequency differed by ETS fusion and/or NKX3-1 CNA status after FDR adjustment for multiple-testing. Dot size and color indicate the difference in proportion, scaled to patients with ETS-negative, NKX3-1–neutral tumors. Background grayscale represents Q values from a proportion test. D, Co-occurrence and associations of driver region pairs across 666 localized prostate tumors. For each pair of driver regions, a hypergeometric test was used to assess whether more mutations were detected than expected by chance alone (co-occurrence) or fewer (mutual exclusivity) after FDR adjustment for multiple-testing (Q < 0.05). The bottom-left heatmap shows all driver pairs; the dotmap on the top right provides effect sizes (dots) and Q values for a subset. Yellow stars on the heatmap mark drivers which are the same gene in clonal and subclonal CNAs. Dot size reflects the difference in driver events, quantified as the observed number minus the expected number. Color indicates deviation direction. The red circle signifies more driver events than expected by chance, whereas the upside-down blue triangle indicates fewer events than expected. E, Clustering of driver regions identifies seven patient subtypes: IMS1–IMS7. Columns are patients. The bottom set of rows shows clinical characteristics, the second set shows mutation densities, and the third shows driver mutations whose frequency differs between subtypes (proportion test; Q < 0.05). The top barplot gives the number of mutated driver regions for each patient. F, Summary subtype profiles showing the proportion of patients in the subtype with certain aberrations. In the positive direction, the proportion of clonal CNA gains, subclonal CNA gains, select GRs, and select SNVs. The lollipops show the proportion of patients for GRs and SNVs. In the negative direction, the proportion of patients in the subtype with clonal and subclonal CNA losses is shown. For each subfigure, CNAs in patients with subclonal PGA >80% were excluded. INV, inversions.

Figure 5.

Mutational hallmarks of prostate cancer grade. A, A linear model was fit to relate each mutational density measure to ISUP GG using tumor and normal sequencing coverage as covariates. Dot size and color represents the effect size for each ISUP GG as a Z-score relative to ISUP GG 1. The barplot to the right shows the Q value from a nonparametric Kruskal–Wallis test. B, Distribution of the number of driver mutations per tumor in each ISUP GG; the median per GG is shown by a black dot. P value is from a one-way ANOVA. C, Genes whose mutation frequency is univariately associated with ISUP GG, ordered by the percentage of samples with mutations in ISUP GG 5 tumors. FDR-adjusted P values from the Pearson χ² test are shown. D, Two-sided Spearman correlation between clinical covariates and measures of genomic instability with dot size showing the magnitude of correlation and background color representing the statistical significance. E, Consensus clustering identified four groups of genes with similar patterns of change across age categories. For each gene cluster, the median mutation frequency for each age category is shown, along with the number of genes in each cluster. F, Venn diagram of the driver genes that were statistically associated with clinical features. G, Cox proportional hazard models were fit for the driver regions that were associated with clinical features. Significant regions after FDR adjustment are shown, as well as the driver type and clinical feature the region was associated with. H,MYC clonal and subclonal gains were associated with biochemical relapse. For each subfigure, CNAs in patients with subclonal PGA >80% were excluded. INV, inversions; WT, wild-type.

Figure 6.

dQTLs bias somatic mutational landscape. A, Schematic of dQTL detection. The PRS used was by Schumacher and colleagues (16). Linear local dQTLs were assessed within ±500 kbp around a driver. Spatial local dQTLs were evaluated using regions defined by RNA Pol-II ChIA-PET profiling in LNCaP, DU145, VCaP, and RWPE-1 cell lines and RAD21 ChIA-PET in LNCaP and DU145 cells. Enhancer regions were defined using H3K27ac HiChIP profiling in LNCaP cells. All discovered dQTLs were tested for replication in six replication cohorts. B, Summary of 26 dQTLs involving 25 unique variants. Dot size and color indicate the magnitude and direction of ORs between the SNP and somatic driver. Background shading indicates P values. Covariate on left indicates type of somatic mutation; the top covariate indicates the analysis strategy for the discovery cohort . C, Comparison of ORs in the discovery vs. replication cohort for tag dQTLs. Horizontal and vertical dotted lines represent OR = 1, and the diagonal line represents y = x. Halo around points indicates replication of direction, diamond around points indicates Q < 0.1 in the replication cohort, and dot color indicates the somatic driver. D and E, Contingency tables for rs11203152 association with clonal loss of TMPRSS2 in (D) discovery and (E) replication cohorts. F and G, Contingency tables of rs848048 associated with SNVs in FOXA1 3′ UTR in (F) discovery and (G) replication cohorts.

Figure 7.

Characterization of dQTLs. A, Summary of all 35 dQTLs involving 25 unique SNPs. Dot size and color indicate the magnitude and direction of association (as OR), and background shading indicates dQTL discovered strategy. B, Forest plot of OR and 95% confidence interval for dQTL associations across 1,991 prostate tumors. Background shading indicates Q < 0.1. The middle covariate indicates the driver mutation, and the right heatmap indicates cohorts included in the analysis. C, Summary of molecular and clinical characterization of dQTLs. Gray indicates dQTL was association with methylation (meQTL), RNA abundance (eQTL), protein abundance (pQTL), transcription factor–binding, histone modification, ISUP GG, BCR, or risk of prostate cancer diagnosis (PCa Risk). Left covariate indicates somatic drivers. D, rs11203152 is located within regulatory dense region. Tracks show chromatin looping anchored by RNAPII, RAD21, AR, or ERG in RWPE-1, LNCaP, VCaP, or DU145 cell lines. E, The number of chromatin loops was higher than expected by chance in LNCaP and VCaP cell lines. Barplots shows the number of anchors within one Mbp of rs11203152. Bottom covariate indicates cell line and target, whereas background shading indicates significant enrichment (Q < 0.05). The red X indicates the expected number of chromatin loop anchors based on 100,000 randomly sampled, equally sized regions. F, dQTLs may explain differences in somatic mutation frequencies across ancestries. Barplot shows the risk of acquiring a FOXA1 SNV or T2E in African (green) or Asian (purple) ancestry relative to European ancestry. The estimated percent of this risk explained by rs848048 (FOXA1) or rs11203152 (T2E) is indicated above the bar. The top covariate indicates ancestry: African in green and Asian in purple. Somatic mutation direction relative to European ancestry is indicated as higher (pink) vs. lower (teal). G, Schematic overview of primary prostate cancer evolution into ISUP GGs. The vertical blue arrow illustrates the temporal relationship between the germline context and driver mutations and the roles they play in tumor evolution. The germline SNPs (bottom) were found to be associated with driver acquisition illustrated by connecting lines. The somatic driver mutation frequency across IMSs is visualized using barplots. The horizontal arrow indicates increasing genomic instability across subtypes and ISUP GGs (9). The Sankey plot connects the IMSs and ISUP GGs, indicating the nondeterministic association between driver acquisition and clinical presentation, while noting the potential role of the tumor microenvironment (90) and epigenetics (45) (G, Created with BioRender.com.)