We used data from the PPCG consortium of primary prostate cancer samples from a total of 1001 prostate cancer donors. Informed written ethical consent was obtained at clinical follow-up, and was consistent with local research ethics and International Cancer Genome Consortium (ICGC) guidelines (https://icgc.org/). Ethical approval was obtained from local research ethical committees.
The full list of approving boards/committees:
Human Research Ethics Committee at Melbourne Health (HREC 2006.073; HREC 2011.009; HREC 2012.275; HREC 2012.220; HREC 2016.087) (Australia).
University Health Network Research Ethics Board and CHU de Québec–Université Laval Research Ethics Board (UHN 06‑0822‑CE; UHN 11‑0024‑CE; CHUQc‑UL 2012‑913:H12‑03‑192) (Canada).
UBC Clinical Research Ethics Board (UBC REB H21‑03722) (Canada).
The National Committee on Health Research Ethics (reference no. 1302791) and notification to the Danish Data Agency (no. 1‑16‑02‑330‑13) (Denmark).
CPP Ile de France IV Institutional Review Board (IRB 00003835) (France).
Ethik‑Kommission der Ärztekammer Hamburg (PV3552; PV4445; PV4679) (Germany).
NHS East of England–Cambridge Research Ethics Committee (REC 3/0180) (UK).
NHS East Midlands Research Ethics Committee (01/4/061) (UK).
NHS London Research Ethics Committee (CCR 2075) (UK).
University of Pretoria Human Research Ethics Committee (FWA00002567; IRB00002235; IORG0001762; approval #43/2010) (South Africa).
St Vincent’s Hospital Human Research Ethics Committee (#SVH/12/231; #SVH/15/227) (Australia).
Human Research Protection Office, US Army Medical Research and Development Command (approvals E02371 and E03280 for genomic interrogation of SAPCS samples) (USA).
Dataset
The study leverages data from the Pan-Prostate Cancer Group (PPCG) cohort, a harmonised, clinically annotated dataset comprising 2,021 prostate cancer donors across seven countries. Multi-modal molecular profiling includes whole genome sequencing (WGS), 450 K DNA methylation arrays, RNA-seq, and where available, long-read Oxford Nanopore Technologies (ONT) sequencing. Detailed clinical annotations including Gleason Grade Group, PSA levels, treatment, and survival outcomes are described in Jakobsdottir et al. 29.
Of the 2021 donors, 1001 had WGS (1209 samples), 1490 had 450 K methylation arrays (1928 samples), and 1292 had RNA-seq (1594 samples) available after quality control.
The details of the PPCG cohort are described in the accompanying marker resource paper29.
WGS data and variant calling
Whole-genome sequencing (WGS) data were generated across multiple international centres as part of the PPCG consortium and processed using standardised, containerised pipelines to ensure reproducibility and harmonisation across sites. Two centralised workflows for somatic and germline variant calling were applied uniformly, as described in Jakobsdottir et al. 29.
Allele-specific somatic copy number alterations (CNAs) were inferred using the Battenberg algorithm (v2.2.9, https://github.com/Wedge-lab/battenberg), which phases heterozygous SNPs using the 1000 Genomes Project reference panel. B-allele frequency (BAF) and logR ratios were then jointly modelled to identify clonal and subclonal CNAs. The method incorporates statistical testing to assess clonality and estimates cancer cell fraction (CCF), tumour purity, and ploidy for each sample.
CNA profiles were generated for all tumour–normal WGS pairs. CNA call sets were manually reviewed by two analysts and samples were classified as satisfactory (n = 1128), unsatisfactory (n = 24), over-fragmented (n = 13), or samples without detectable copy-number alterations (n = 301). Only high-confidence CNA profiles (satisfactory samples) were retained for downstream analyses in Epi2Hit.
Hemizygous loss on autosomes were defined as segments with total copy number = 1 and CCF ≥ 0.9, ensuring high-confidence clonal deletion events. These calls rely on Battenberg’s allele-specific modelling and purity/ploidy-aware segmentation, providing robust and consistent annotation of deletion events for integration with methylation and gene expression data.
Further processing and harmonisation of the PPCG WGS data are described in detail in the accompanying marker resource paper29
DNA methylation
450 K array-based methylation data were processed using a standardised and publicly available RnBeads-based workflow100,101 (https://github.com/panprostate/PPCG_DNA_methylation), which integrates preprocessing, normalisation, and QC steps designed to harmonise large multi-cohort data, as described in Jakobsdottir et al. 29. Briefly, Raw IDAT files were imported and preprocessed using RnBeads v2.15.1 with background subtraction, dye-bias correction, and bead count-based filtering. Low-quality samples (median bead count <3 or high probe-wise detection p-values > 0.05) were flagged and excluded. Probes were filtered out if they overlapped known SNPs (MAF > 0.05), were known to be cross-reactive102, or had detection p-values ≥ 0.05 in ≥1% of samples. Methylation data were normalised using a single-sample approach implemented via the SeSAMe method103, re-implemented in the RnBeads framework. Background correction was performed using the “noobsb” method, which fits a sample-specific background distribution based on out-of-band signal intensities from Type I probes and regresses this from foreground intensities. Dye-bias correction was conducted using internal channel scaling. This normalisation strategy enabled robust cross-cohort comparison without requiring batch metadata. Technical consistency across countries and cohorts was evaluated using density plots of 65 SNP-tagging CpG probes. Additionally, beta-value distributions were visually inspected, and principal component analysis (PCA) of normalised values confirmed that primary axes of variation reflected biological rather than technical differences. To further mitigate population structure effects, genotype-derived SNP principal components were included in downstream analyses.
Reference methylation data for major immune and stromal cell types were compiled from public datasets and in-house FACS-sorted prostate epithelial populations. To distinguish among 11 major cell types, the 200 most differentially methylated CpG sites per cell type (100 hypermethylated and 100 hypomethylated), along with CpGs within 50 bp, were selected104, resulting in a set of 2851 CpGs. Cell-type proportions were then estimated using the Houseman algorithm105,106 as implemented in RnBeads, and normalised to sum to 100% per sample for downstream adjustment in methylation-expression analyses. Estimated cell-type proportions were used as covariates in downstream methylation-expression modelling within the Epi2Hit framework.
Further processing and harmonisation of the PPCG Methylation data are described in detail in the accompanying marker resource paper29.
Principal component analysis is shown in Supplementary Fig. 1a
RNA sequencing data
RNA-seq data were uniformly processed across samples using a harmonised pipeline based on STAR for alignment to the GRCh37 reference genome, followed by transcript-level quantification with Salmon (v1.4.0) using GENCODE v19 annotations (https://github.com/panprostate/RNA).
RNA-seq data from the PPCG cohort, generated across several centres using diverse chemistries and platforms, underwent a multi-step normalisation strategy to correct for both known and unknown sources of technical variation. Lowly expressed genes were filtered, retaining 22760 transcripts for downstream analysis. Tumour purity was estimated using ESTIMATE107,108 and singscore methods.
Normalisation was performed using the RUV-III method109 with pseudo-replicates of pseudo-samples (PRPS), implemented in both supervised and unsupervised modes. In the supervised approach, known biological covariates (ERG status, tumour purity) guided the selection of negative control genes. In the unsupervised approach, mutual nearest neighbours (MNN) and KNN-based algorithms were used to generate pseudo-replicates and identify control genes based on low biological variability and high association with technical noise. RUV-III-PRPS outperformed other tested methods (CPM, TMM, upper-quartile, VST, ComBat) in preserving biological signals while removing technical artefacts. All normalisation code is publicly available at https://github.com/RMolania/RUVprps.
Further processing and harmonisation of the PPCG RNA data are described in detail in Jakobsdottir et al. 29.
Principal component analysis is shown in Supplementary Fig. 1a
Multi-omic data integration strategy
Due to the sparsity data representation, a stepwise integration strategy was employed to maximise the use of available data while minimising loss of statistical power. First, regulatory CpGs were identified using the full set of available DNA methylation profiles (n = 1268) together with ATAC-seq data from an independent cohort (GSE188797) and further annotated with external Hi-C-based chromatin loops (GSE164347). Next, expression-correlated CpGs (eCpGs) were determined within the subset of samples for which both methylation and RNA-seq data were available (n = 1292). Finally, candidate biallelic inactivation events were evaluated by integrating methylation data with WGS-derived copy number profiles (n = 1001). This modular approach enabled robust cross-modal analysis without restricting the study to only the samples profiled in all three data types.
Hi-C data
Hi-C data from prostate cancer were retrieved from the Gene Expression Omnibus (GEO) from data sets GSE16434737.
For integration with Oxford Nanopore (ONT) long-read sequencing data, which was processed using the hg38 genome build, we performed liftover of the relevant Hi-C loop regions from hg19 to hg38 using pyliftover v0.4.1 to enable visualisation and plotting in the hg38 coordinate space.
ChIP-seq data
Processed ChIP-seq data for ZFHX3 was obtained from the Gene Expression Omnibus (GEO) (GSE49402, GSM1208743; Yan et al. 2013)72 and from human neural stem cells (Pérez Baca et al., 2024)71,72. To provide prostate-specific evidence, we additionally used published ChIP-PCR data in C4-2B prostate cancer cells (Hu et al., 2019)42,43, where primers targeting three regions of the MYC promoter (Regions A–C) were designed. These primer sequences were validated in silico using the UCSC isPcr tool against the hg38 assembly, which confirmed single, specific amplicons within the MYC promoter (chr8:127,735,377–127,735,930). The identified regions overlapped with ZFHX3 binding peaks detected in neural stem cells and colon cells, confirming promoter-proximal binding across distinct cellular contexts.
liftover from hg18 to hg19 and from hg19 to hg38 was performed using pyliftover 0.4.1.
Epigenetic marks (H3K4me3, H3K27ac) and CTCF binding sites at the loci were derived for prostate cancer cells (PC-3, C4-2B) from the ENCODE project110
ATAC sequencing data
To identify CpG sites located within open chromatin regions, processed ATAC-seq data generated from human tumour samples was retrieved from the Gene Expression Omnibus (GEO) from data set GSE18879736
Long-read sequencing
Tumour and matched normal genomic DNA from six prostate cancer donors were prepared for Oxford Nanopore long-read whole-genome sequencing using the Ligation Sequencing Kit V14 kit (SQK-LSK114). The libraries were sequenced on a PromethION P24 instrument with R10.4.1 flow cells (FLO-PRO114M). Average library sequencing depth was 30x and 18x, range 19.5x-48.6x and 6.0x-48.0x, for tumour and normals, respectively. Basecalling was performed with Dorado v1.0.0 with model hac@v5.2.0, while 5mC and 5hmC modifications were called concurrently with model sup@v5.2.0_5mC_5hmC@v1. The data was then aligned to hg38 and phased with the epi2me-labs wf-human-variation workflow v2.7.1, performing genome-wide allele specific CpG methylation calling (5mC and 5hmC) with modkit v0.3.3.
For ZFHX3, phased haplotypes allowed us to distinguish the deleted allele (allele 1, copy number = 0) from the retained allele (allele 2). The expected pattern for biallelic inactivation is (i) absence of methylation signal on the deleted allele and (ii) strong hypermethylation on the retained allele.
Identification of regulatory CpG sites (eCpGs)
DNA methylation data
Methylation data was obtained from the 450k array methylation dataset provided by the PPCG. The data included β-values corresponding to the methylation levels at various CpG sites across multiple samples. Samples that failed quality control were excluded from the analysis. The remaining 1268 primary tumour samples were used for the analysis.
Annotation of CpG sites with gene information
Each site was annotated with the nearest gene(s) using the hg19 reference genome. This annotation was performed using the pybedtools package and the ‘closest’ function. The reference gene annotations were retrieved from the GENCODE v19 database. The CpG sites were matched to the closest genes, considering both strand orientation and distance. The output included all potential gene associations for each CpG site, facilitating further exploration of gene-specific methylation patterns.
Filtering based on methylation correlations
The dataset was subjected to a series of filtering steps to prioritise CpG sites with specific methylation patterns:
-
Correlation with immune and stromal components: CpG sites were filtered based on their correlation with immune and stromal cells abundance. Pearson correlation coefficients were calculated for each site using methylation beta-values against a methylation-derived estimation of stromal and immune cell type proportions (described in detail in Jakobsdottir et al. 29. Sites showing low correlation (correlation coefficient <0.4) with these components were retained, ensuring that the selected CpG sites primarily reflect methylation changes in tumour cells rather than in the surrounding stromal or immune cells.
-
Correlation with global loss of methylation: Additionally, CpG sites were filtered based on their correlation with global loss of methylation estimated by the EpiCMIT-hypo score111 (as described in Jakobsdottir et al.29.), reflecting epigenetically-determined Cumulative MIToses which result in stochastic loss of methylation and the formation of’partially methylated domains (PMDs).
ChromHMM annotation for regulatory elements
Filtered CpG sites were annotated with chromatin states using prostate-specific ChromHMM data from Pomerantz et al. 38. ChromHMM annotations were intersected with CpG sites using pybedtools, and states were mapped to corresponding categories (e.g., Active Prostate Lineage-Specific Promoter). This step identified CpG sites within enhancer and promoter regions.
Merging CpG sites by methylation correlation
CpG sites annotated within the same regulatory region were merged if their methylation levels exhibited a Pearson correlation coefficient above 0.6. Merged sites were represented by the average methylation value and spanned the combined genomic coordinates of the individual CpG sites.
Hi-C data integration for chromatin loops
Hi-C data from prostate cancer (GEO; GSE164347) was used to annotate regulatory CpG sites with chromatin loop information. CpG sites overlapping loop anchor regions were identified, and corresponding loop IDs were assigned to these sites. Sites associated with loops were indicated by a ‘Ls_’ prefix in their identifiers.
Linear regression analysis to identify expression-associated CpG sites. We performed linear regression to identify CpG sites associated with gene expression (eCpGs) by regressing methylation levels against RNA-seq-based gene expression from the PPCG group.
-
Regression modelling: ordinary Least Squares (OLS) regression was performed using the statsmodels package, treating methylation as the independent variable and gene expression as the dependent variable. Regression coefficients, p-values, confidence intervals, and R-squared values were extracted.
-
Multiple testing correction: P-values were adjusted for multiple comparisons using the Benjamini-Hochberg False Discovery Rate (FDR) method from the statsmodels package (Ruiz-Arenas et al. 112; Advani et al., Nat Commun 2024113; Ando et al., Nat Commun 2019114), Supplementary Data 2.
-
eCpG identification: CpG sites with significant negative associations (adjusted p-value < 0.05) were classified as eCpGs. To ensure biological relevance, we further applied a regression coefficient threshold of < −0.4, which corresponds to an approximate 25% reduction in gene expression across the methylation range. Since gene expression was modelled as log2(TPM + 1), a coefficient of −0.4 implies a moderate-to-strong transcriptional repression effect per unit increase in methylation. This level of effect size is consistent with enhancer and promoter methylation–expression associations previously reported in cancer cohorts99 (Supplementary Data 2).
Sensitivity analyses for tumour purity and global genomic instability
Purity-adjusted methylation–expression regression
To assess the robustness of methylation-expression associations to tumour purity, we re-fit the regression models used to define eCpGs. In the baseline analysis, for each CpG–gene pair we modelled gene expression as a function of methylation using ordinary least squares (OLS):
$${expression}\,=\,\alpha \,+\,\beta \cdot {{\rm{\cdot }}}{methylation}\,+\,\varepsilon$$
(1)
$${expression}\,=\,\alpha \,+\,\beta \cdot {{\rm{\cdot }}}{methylation}\,+\,\gamma \cdot {{\rm{\cdot }}}{purity}\,+\,\varepsilon$$
(2)
For all CpG-gene pairs considered in the main analysis we extracted the methylation coefficients (β) and R2 values from both models and compared their distributions. The direction of the methylation effect was preserved and the association remained significant after Benjamini–Hochberg false discovery rate (BH-FDR) correction (q < 0.05). Genome-wide concordance of methylation effect sizes between baseline and purity-adjusted models was summarised using Pearson and Spearman correlations (Pearson r = 0.996, Spearman ρ = 0.996) (Supplementary Fig. 1e).
To further account for global genomic instability, we additionally fit a two-covariate model including both tumour purity and fraction genome altered (FGA):
$${expression}\,=\,\alpha \,+\,\beta \cdot {{\rm{\cdot }}}{methylation}\,+\,\gamma \cdot {{\rm{\cdot }}}{purity}\,+\,\delta \cdot {{\rm{\cdot }}}{FGA}\,+\,\varepsilon$$
(3)
Methylation coefficients from this model were compared to the baseline model in the same way (Pearson r = 0.969, Spearman ρ = 0.964) (Supplementary Fig. 1f).
For loci highlighted in Fig. 3a (ZFHX3, PDE4D, PSD3, KLF5, NAT1, APC and RBPMS), we report locus-level regression results with and without purity adjustment, including methylation coefficients, p-values, FDR-adjusted q-values and R2 (Supplementary Table 2), to illustrate the stability of effect direction and magnitude at key Epi2Hit loci.
Categorisation of methylation Levels
Methylation levels for each feature were categorised into three groups: ‘Not/Low_methyl’, ‘Medium_methyl’, and ‘Highly_methyl’. The classification was based on predefined thresholds (down_threshold=0.2 and up_threshold=0.7)115,116,117,118 and the distribution of methylation values. Features with all values below the down_threshold were categorised as ‘Not/Low_methyl’, while those with values above the up_threshold were labelled ‘Highly_methyl’. Features with methylation values spanning these thresholds were further analysed using kernel density estimation (KDE) and peak detection to identify potential multimodal distributions.
For features with intermediate methylation levels, a Gaussian kernel density estimation (KDE) was applied to estimate the probability density function (PDF) of the methylation values.
The KDE was performed using the scipy.stats.gaussian_kde function.
Peaks in the KDE were identified using the scipy.signal.find_peaks function. Depending on the number of detected peaks, the feature’s methylation values were categorised accordingly:
To ensure reproducibility and robustness, we conducted a sensitivity analysis across three key KDE parameters: distance, height, and prominence, which define the minimum distance between peaks, the KDE height threshold, and the peak prominence threshold respectively.
A total of 305 CpG sites were selected from the prostate cancer cohort, each visually confirmed to exhibit bimodal or trimodal methylation distributions. We then systematically evaluated 125,000 parameter combinations generated from a three-dimensional grid (50 linearly spaced values from 0.01–0.5 for each parameter; see Supplementary Data 3).
Performance was measured by the number of CpGs for which the algorithm correctly detected 2–3 peaks. The optimal parameter range was defined as:
-
distance_frac: 0.01–0.11
-
height_quantile: 0.01–0.33
-
prominence_frac: consistently optimal at 0.01
The midpoint values (distance_frac=0.06, height_quantile=0.17, prominence_frac=0.01) were selected as defaults for downstream analyses and recovered 94% (287/305) of the benchmark regions.
Pairwise KDE parameter performance is visualised in Supplementary Fig. 2d, where warmer (red) regions indicate higher peak detection accuracy across the parameter grid.
This parameter configuration was further validated on an independent luminal breast cancer cohort (TCGA-BRCA119), where an additional 450 manually selected CpG sites were analysed. Using the same settings, the KDE algorithm correctly identified 1, 2 or 3 peaks in 443/450 (98.4%) regions, confirming the sensitivity of the approach (Supplementary Data 4-5).
Depending on the number of detected peaks, the feature’s methylation values were categorised accordingly:
-
Single peak: a Gaussian Mixture Model (GMM) with three components was fitted to the methylation values using sklearn.mixture.GaussianMixture(https://scikit-learn.org/stable/modules/generated/sklearn.mixture.GaussianMixture.html). The resulting clusters were used to assign methylation categories based on the cluster means, with thresholds calculated from the median and minimum of these means.
We chose GMM for unimodal distributions based on its successful use in prior studies116, where it stratified promoter methylation into binary states. In contrast, our method applies GMM at the CpG level and partitions the distribution into three states: ‘Not/Low_methyl’, ‘Medium_methyl’, and ‘Highly_methyl’, capturing partial methylation patterns frequently observed in cancer. GMM is only used when KDE fails to detect more than one peak, ensuring method selection is driven by the distributional shape.
A visual overview of the KDE vs. GMM logic, is shown in Supplementary Fig. 2a.
-
Two peaks: features displaying a bimodal distribution were manually reviewed, with the lower and upper peak values serving as thresholds for categorising the methylation levels into ‘Not/Low_methyl’, ‘Medium_methyl’, and ‘Highly_methyl’.
-
Three peaks: features with a trimodal distribution were also manually reviewed, using the median and minimum of the peak values as thresholds for classification.
If no peaks were detected or distributions were ambiguous, conservative fallback quantile thresholds (quant_low=0.3, quant_high=0.7) were used.
To identify genes that are biallelic inactivated we used our previously published methodology14,39,120, expanding it to DNA methylation data.
Loss/Methylation: A heterozygous loss combined with high methylation (som_loss/methyl).
Validation of Epi2Hit in TCGA-PRAD
TCGA-PRAD data acquisition and preprocessing
For external validation, we applied Epi2Hit to an independent prostate cancer cohort from TCGA-PRAD. DNA methylation (450k array β-values), RNA-seq expression, and copy-number alteration (CNA) profiles for primary tumours (n = 493) were downloaded from cBioPortal(https://www.cbioportal.org/), and structural variant (SV) calls were obtained from TCGA-PRAD whole-genome sequencing samples in the NCI Genomic Data Commons (GDC(https://portal.gdc.cancer.gov/); accessed 22 July 2025). Only samples with matched methylation, expression and CNA data were retained. Gene expression values were transformed to log2(TPM + 1).
Re-implementation of the Epi2Hit workflow in TCGA
We re-implemented the complete Epi2Hit pipeline on TCGA-PRAD using the same parameter settings as for PPCG. Briefly, CpG sites were annotated to genes, filtered (removal of CpGs correlated with immune/stromal estimates or global PMD-associated hypomethylation, restriction to prostate ChromHMM promoters/enhancers, intersection with prostate tumour ATAC-seq peaks, and prioritisation of CpGs at Hi-C loop anchors), and tested for association with gene expression to define eCpGs. CpG-level p-values were adjusted using the BH-FDR, and within-gene multiplicity was controlled using Simes-based aggregation121 of eCpG p-values to obtain per-gene q-values. For each gene, methylation distributions at eCpGs were then modelled with KDE/GMM as described above to obtain discrete methylation states, which were integrated with CNAs and SVs to call Epi2Hit events (heterozygous loss plus hypermethylation).
Comparison of biallelic inactivation patterns across cohorts
For the most frequently disrupted tumour suppressor genes (TSGs), we quantified the number of samples in each cohort harbouring genomic biallelic hits versus Epi2Hit events These counts were summarised as bar plots for TCGA-PRAD and PPCG (Supplementary Fig. 3c, Supplementary Table 4-5).
MethylMix validation
Construction of gene-level methylation profiles
To compare Epi2Hit with a MethylMix framework for methylation-driven genes, we used the pre-eCpG set defined in the PPCG discovery cohort (prior to KDE/GMM-based peak detection).
Gene-level methylation and expression matrices were analysed with the MethylMix algorithm (R package MethylMix 2.0(https://bioconductor.org/packages/release/bioc/html/MethylMix.html)) using recommended settings.
Overlap with Epi2Hit and visualisation of methylation states
To quantify concordance, we intersected the list of methylation-driven genes with the set of genes harbouring Epi2Hit events in PPCG. Overall, 47% of methylation-driven genes were also identified by Epi2Hit, including key tumour suppressors such as APC, KLF5, NAT1 and ZFHX3. For these four genes, we visualised the MethylMix beta mixture models as density curves overlaid on β-value histograms (Supplementary Fig. 3d), illustrating the presence of 2–3 distinct methylation states per gene.
Standardised mean difference analysis (Cohen’s d)
We first used publicly available ChIP-seq data (GSE49402) to identify genes with ZFHX3 binding sites located within ±2 kb of their annotated transcription start sites, representing putative promoter-proximal targets of ZFHX3. From this list, we selected genes previously reported in the literature to have tumour-suppressive or oncogenic functions in prostate cancer. This final set of genes was used for expression analysis.
To quantify the differences in gene expression between ZFHX3 wild-type (WT) and epigenetic biallelic inactivation (Epi2Hit) prostate cancer samples, we computed the standardised mean difference (Cohen’s d) for the preselected genes. For each gene, we calculated the mean expression in the WT and Epi2Hit groups and divided their difference by the pooled standard deviation across the two groups. This approach allows comparison of expression differences across genes on a standardised scale, independent of absolute expression values.
The pooled standard deviation (sp) was computed as:
$${s}_{p}=\,\sqrt{\frac{({n}_{1}-\,1){s}_{1}^{2}\,+\,({n}_{2}\,-\,1){\,s}_{2}^{2}}{{n}_{1}\,+\,{n}_{2}\,-\,2}}$$
(4)
where n1 and n2 are the sample sizes of the WT and Epi2Hit groups, respectively, and s1 and s2 are their respective standard deviations.
Cohen’s d is then:
$$d=\frac{{\underline{x}}_{1}-{\underline{x}}_{2}}{{s}_{p}}$$
(5)
Positive values of d indicate higher expression in the WT group, and negative values indicate higher expression in the Epi2Hit group.
Essential genes analysis
To investigate the association between essential genes and tumour suppressor genes (TSGs) that are rarely or never inactivated by bi-allelic loss, we utilised data from the DepMap(https://depmap.org/portal/) database. Genes with scores below this threshold (<−1.2) were considered essential.
Calculation of proximity to essential genes
The distance between each TSG and the closest essential gene was calculated as the absolute genomic distance in megabases (Mb), considering both upstream (5′) and downstream (3′) directions. For each TSG, the minimum of the two distances (5′ or 3′) to the nearest essential gene was used. The presence of an essential gene near a TSG was annotated as True if the distance between the TSG and the nearest essential gene was less than the first quantile (Q1) of the mean distance (dTSG
For comparative analyses, we included only TSGs affected by Epi2Hit or homozygous deletion in more than two tumour samples. Differences in the distance to essential genes between these groups were assessed using the two-sided Mann–Whitney U test.
To explore functional consequences, we also compared average gene expression (log2(TPM + 1)) between essential genes and TSGs using RNA-seq data from the same cohort. Statistical significance was evaluated using the Mann–Whitney U test (Fig. 5e).
Structural covariates and multivariable models
To assess whether the proximity signal could be explained by structural genomic features and to address potential confounding, we assembled a unified gene-level dataset for all TSGs that included:
-
Gene length (bp, converted to kb),
-
Local gene density in a ± 250 kb window (gene_density_500kb),
-
Distance to the centromere (dist_to_centromere, in kb),
-
Proximity to essential genes at several DepMap CERES thresholds.
For each DepMap essentiality cutoff X∈{CERES < − 0.2, < − 0.5, < − 0.7, < − 1.0, < − 1.2, < − 1.5}, we computed:
-
dist_CERES_X: distance (kb) to the nearest essential gene,
-
count_CERES_X: number of essential genes within 500 kb,
-
has_dist_CERES_X_within_500kb: indicator of ≥1 essential neighbour within 500 kb.
All distance-like variables were rescaled to kilobases for interpretability. These variables were used as predictors in multivariable logistic regression models with outcome Epi2Hit vs homozygous loss for each CERES threshold. Each model included gene length, local gene density, distance to centromere, number of nearby essential genes (count_CERES_X), distance to the nearest essential gene (dist_CERES_X) and the 500 kb proximity indicator. Odds ratios (ORs) and 95% confidence intervals are summarised in Supplementary Fig. 5a, b and Supplementary Table 8. At the main common-essential threshold (CERES <–1), both the number of essential neighbours and, to a lesser extent, their proximity remained significantly associated with Epi2Hit status after adjustment for gene length, gene density and chromosomal context.
Prostate-specific essentiality
To evaluate whether these patterns hold in a prostate-specific context, we repeated the analysis using only prostate cancer cell lines (22RV1, LNCaP, DU145) from DepMap. Essentiality scores were recalculated across the same CERES thresholds, and prostate-specific proximity measures (count_CERES_X, dist_CERES_X, has_dist_CERES_X_within_500kb) were recomputed for each TSG. Multivariable logistic models with the same covariates (gene length, gene density, distance to centromere and essential-gene proximity) were then refitted. The number and proximity of prostate-essential neighbours remained strong predictors of Epi2Hit versus homozygous loss across thresholds (Supplementary Fig. 5c, d, Supplementary Table 9), supporting a genuine constraint by essential neighbours’ effect rather than structural artefacts or dependence on a particular essentiality definition.
Copy number loss frequency (ideogram-level)
For each ideogram interval (chromosome, start, end; hg19 autosomes), we calculated the proportion of tumour samples (n = 991) with overlapping homozygous loss (norm_total = –2) or hemizygous loss (norm_total = –1) in the segmentation data. Overlap was defined as segment_start region_start. Homozygous loss (hom %) and hemizygous-only loss (het %; excluding samples that also carried a homozygous loss in that region) were computed separately. Their sum (loss %) represented the total fraction of samples with any loss in the interval, without double-counting. Only intervals with loss % > 10% of samples were retained for ideogram plotting.
The chromosome ideogram was generated using the Python package quyuan (v1.1.3; https://github.com/tcztzy/quyuan).
Copy number loss frequency (chromosome-level)
For each autosomal chromosome, we determined the proportion of tumour samples (n = 991) with at least one segment meeting the loss criterion anywhere on that chromosome. Homozygous loss (hom %, norm_total = –2) and hemizygous loss (het %, norm_total = –1) were calculated separately, and the union of these sample sets (loss %) represented the total percentage of tumours with any loss on that chromosome, without double-counting those carrying both types of loss.
Exclusivity score of copy number alterations
To quantify the tendency of driver genes to occur independently rather than co-occur with other driver genes across cancer samples, we applied the exclusivity score metric, as described by Gerhauser et al. (Cancer Cell, 2018)14. This score estimates the degree to which a given gene alteration is exclusive or co-occurs with other gene alterations.
For each driver gene, the exclusivity score is calculated as the mean inverse of the number of co-altered driver genes in samples where the gene is present. Specifically, for each sample with an alteration of the target gene, the reciprocal sum of driver genes with the alteration is computed and the exclusivity score is calculated as the arithmetic mean of the reciprocal values across all relevant samples.
Formally, the exclusivity score is defined as:
$${Exclusivity}\,{Score}\,=\,\frac{1}{n}{\sum }_{i=1}^{n}\frac{1}{{c}_{i}}$$
(6)
Where:
\(n\) is the number of samples containing the target gene,
\({c}_{i}\) is the number of driver genes altered in a sample
The exclusivity score ranges from 0 to 1. Higher values indicate greater exclusivity (occurring as the only driver alteration), while lower values suggest frequent co-occurrence with other driver genes.
Survival analysis
Lifelines package (lifelines.KaplanMeierFitter and lifelines.CoxPHFitter) (https://github.com/lifelines/lifelines) was used for survival analysis across subgroups.
Statistics & reproducibility
All statistical tests were two‑tailed, unless otherwise specified. For parametric data, we used tests assuming normal distribution. For non‑parametric data, we used distribution‑free tests. Significance thresholds and specific tests are detailed in the figure legends and Methods.
The Epi2Hit pipeline was implemented in Python and comprises modular steps for methylation processing, expression-methylation correlation, probe filtering and merging, gene annotation, chromatin state mapping, and integration with Hi-C, ATAC-seq, and long-read sequencing data. Default parameters and thresholds (CNV, KDE, regression covariates) are detailed in Supplementary Table 11.
The pipeline supports parallel processing through the joblib library (e.g., n_jobs=16) to accelerate steps such as correlation-based probe merging. All analyses were run on a shared Linux server with 64 CPU cores and 512 GB of RAM. Jobs were typically submitted with a memory request of up to 40 GB RAM. Runtime for the full pipeline applied to the PPCG cohort (~1000 samples and ~450,000 probes) was approximately 3-4 h, depending on data modalities and filtering settings. Due to its modular and parallelizable structure, the pipeline scales efficiently with increasing cohort size and is readily applicable to other tumour types and multi-omic datasets.
Epi2Hit is distributed as a Docker image available at https://hub.docker.com/r/al0ka/epi2hit. The container includes all necessary dependencies and project scripts and is available for both x86_64 (linux/amd64) and ARM64 (linux/arm64) architectures. Users can obtain the image by executing “docker pull al0ka/epi2hit:latest” and then run the workflow by following the included documentation. The image was tested using the included example data on an Ubuntu 24.04.3 LTS machine with 4GB of RAM and 4 CPU cores. The analysis consisted of three steps: the initial step (eCpG identification – run-epi2hit_py) took approximately 7 min, the second step (generating metadata.yaml) was nearly instantaneous, completing in about 1 s, and the final step (running biallelic_py) took approximately 1 min.
All data analyses were performed using Python 3.8 (pandas v1.1.5, numpy v1.19.5, lifelines 0.27.7, scipy v1.5.4, and scikit-learn v0.23.2, statsmodels v0.13.2, joblib v1.2.0, tqdm v4.64.1, plotnine v0.10.1, pybedtools v0.10.0, pyBigWig v0.3.22, statannot v0.2.3) and and R 4.2.2 (MethylMix v2.41.0).
The plots and visualisations were generated using the seaborn v0.11.0, matplotlib v3.6.3, and additional figures were created using BioRender (BioRender.com) and quyuan (v1.1.3).
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
