Reconstructing protein interactions at enhancer-promoter regions in prostate cancer

doi:10.21203/rs.3.rs-2252897/v1

Download PDF

Article

Reconstructing protein interactions at enhancer-promoter regions in prostate cancer

https://doi.org/10.21203/rs.3.rs-2252897/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 06 Dec, 2023

Read the published version in Nature Communications →

Version 1

posted

You are reading this latest preprint version

DNA-binding proteins (DBPs) and in particular transcription factors interact with enhancers and their target genes through enhancer-promoter (E-P) interactions. Technological advancements such as chromosome conformation capture allow to identify E-P interactions, but the protein networks involved have not yet been characterized. Most importantly, the role of nuclear protein networks in human diseases has been so far poorly investigated. Prostate cancer (PrCa) heritability is associated with variations in enhancers that affect specific gene expression. Here, we introduce a novel approach, called Promoter-ENhancer-GUided Interaction Networks (PENGUIN), to identify protein-protein interactions (PPI) in E-P interactions and apply it to our PrCa dataset. PENGUIN integrates chromatin interactions between a promoter and its enhancers defined by high-coverage H3K27ac-HiChIP data, with a tissue-specific PPI network inferred from DNA-binding motifs and refined with gene expression. Among a total of 4,314 E-P networks, PENGUIN performed unsupervised clustering. We functionally validated this clustering procedure by searching for enrichments of specific biological features. We confirmed PENGUIN structural classification of E-P networks by showing a clear differential enrichment of the architectural protein CTCF. Next, and directly related to our PrCa case study, we observed that one of our 8 main clusters, containing 273 promoters, is particularly enriched for PrCA associated single nucleotide polymorphisms (SNPs) and oncogenes. Our approach proposes a mechanistic explanation for 208 PrCa SNPs falling either inside the binding sites of DNA-binding proteins (DBPs) or within genes encoding for intermediate proteins bridging E-P contacts. PENGUIN not only confirmed the relevance of key regulators in PrCa, but also identified new candidates for intervention, opening up new directions to identify molecular targets for disease treatment.

Biological sciences/Computational biology and bioinformatics/Genome informatics

Biological sciences/Genetics/Cancer genetics

Biological sciences/Molecular biology/Transcription

Prostate cancer (PrCa) is the 2nd most common cancer in men ¹. Its distinct hormone-dependent nature is characterized by high expression and frequent genetic amplification of AR, regulator of homeostasis and proteases transcription, such as KLK3 encoding PSA (Prostate-Specific Antigen), and principal therapeutically targeted oncogene in PrCa ². Increased genetic instability resulting in chromosomal rearrangements and high frequency of mutations are deemed indicative of PrCa aggressiveness ³ for which there is need of ad hoc treatments ⁴. Recurrent mutations in FOXA1, involved in prostate organogenesis and regulator of AR transcription, have been observed in several populations ^5,6. Hundreds of PrCa-associated single nucleotide polymorphisms (SNPs) have been identified by genome-wide association studies (GWAS), including genomic regions within tumor suppressor genes and oncogenes, such as MYC ⁷. However, the functional relationship between most of these SNPs and PrCa pathophysiology is unknown. This missing part of the picture, together with the growing evidence of abnormal transcriptional programs driven by genetic instability, led us to investigate the role of chromatin architecture in PrCa. In particular, we studied the nuclear proteins potentially involved in transcriptional regulation through the interaction of promoters and non-coding regulatory elements, enhancers.

Enhancer-promoter (E-P) interactions have an important role in gene regulation. DNA-binding proteins (DBPs), such as transcription factors (TFs), regulate gene expression by binding promoters to enhancers sometimes through intermediate proteins. This chromatin interaction between E-P can be mediated by co-activators (e.g., mediators), chromatin structural proteins (e.g., cohesin), noncoding RNA-binding proteins, and others. Disruption of these interactions centered in a single promoter (either protein-protein or protein-DNA interactions), which we collectively call E-P protein-protein Interaction Networks (EPIN), is increasingly linked to diseases such as cancer ^8–10. Enhancers are often the target of sequence and structural variation in cancer, specifically deregulation of TFs and chromatin modifiers ¹¹, and represent promising pharmacological targets in PrCa ¹². To date, techniques such as HiC or, more specifically, its recent derivative HiChIP ¹³, HiC combined with ChIP-seq, are able to identify specific chromatin interactions between a promoter and its enhancers. For instance, H3K27ac-HiChIP is designed to specifically detect and amplify E-P interactions and it has been recently employed by our lab to uncover susceptibility genes in PrCa ¹⁴.

Although both protein-protein and protein-DNA interactions play central roles in E-P interactions, previous analyses focused on one or the other aspect. Indeed, a number of studies have focused on DBPs networks ^15,16. A more recent analysis exploited PPIs with chromosome conformation capture assays to facilitate prioritization of functional interactions ¹⁷. At present, the characterization of context specific intermediate protein-protein interactions (PPIs) involved in disease pathways and their association with DBPs remains largely unanswered ¹⁸.

To characterize protein interactions that take place at the E-P contacts in PrCa, we developed the Promoter-ENhancer-GUided Interaction Networks (PENGUIN) approach. For each promoter annotated in the genome and covered by at least one HiChIP interaction, PENGUIN builds an EPIN by integrating several sources of information: (1) high-resolution chromatin interaction maps enriched for a marker of active promoter-enhancer activity (H3K27ac-HiChiP); (2) gene expression (RNA-sequencing); (3) tissue-specific physical nuclear PPIs; (4) high-quality curated binding motifs of protein-DNA interactions. The PrCa specificity of this dataset is given by the H3K27ac-HiChiP of androgen-sensitive human prostate adenocarcinoma cells (LNCaP) and RNA-seq form the same cell line, and prostate-specific PPIs and DNA binding motifs extracted from publicly-available datasets (Methods).

PENGUIN is a software that can be applied to other diseases. The method can identify clusters of PPI networks dependent on the potential binding proteins in promoters and enhancers that are enriched for particular annotations like GWAS and CTCF. It offers a comprehensive framework to integrate disease data and provides insights into protein networks found in E-P interactions for further molecular validation. For instance, PENGUIN can be used to identify trans-acting factors (e.g., interaction cascades of TFs and chromatin regulators) that could be targeted by drugs, or cis-acting factors (e.g., DBPs with binding motifs in regulatory elements) whose DNA binding affinity could be modified through knock-outs via CRISPR for therapeutic intervention. Moreover, unlike traditional TF enrichment analysis which detects general enrichments of particular proteins, PENGUIN can help identify the specific protein cascade potentially disrupted at enhancer loci for the disease under study. This methodology identifies new directions in the molecular characterization of chromatin interactions as well as in the definition of potential targets for molecular screening towards disease treatment.

In this paper, we propose two main applications of PENGUIN. First, we used it to facilitate the identification of key factors that may play a role in transcriptional regulation of PrCa. By clustering together promoters with similar EPINs, we identified 273 promoters whose genes are enriched in PrCa GWAS, known PrCa oncogenes, and ChIP-validated binding sites of transcriptional repressor CTCF. The proteins that populate such EPINs constitute putative PrCa-related bioentities, some of which have not been previously described to be associated with PrCa SNPs or oncogenes. Second, the EPINs detected by PENGUIN enable the characterization of distinct molecular cascades potentially affected by PrCa SNPs at E-P contacts. These represent potentially new molecular targets in PrCa that cannot be identified through conventional analytical procedures, such as E-P contacts and GWAS overlap. Finally, we provide a dedicated web server to explore the results at https://penguin.life.bsc.es/.

PENGUIN identifies clusters of protein interactions based on chromatin contacts

We leveraged 24,547 E-P contacts (30,416 after refinement and prioritization, Methods) identified using H3K27ac-HiChIP data in LNCaP, 810 binding motifs from 639 DNA-binding proteins, and 31,944 prostate-specific, experimentally validated, physical and nuclear PPIs (filtering out proteins from unexpressed genes, Methods) to construct 4,314 EPINs using the PENGUIN pipeline outlined in Fig. 1 (Methods). Each EPIN is centered around one promoter that we found to be contacted by a median of 4 enhancers, with a maximum of 93 enhancers for the promoter of the gene CRNDE (Table S1). Altogether, the 4,314 EPINs contain a total of 8,215 interactions (edges) among a total of 885 proteins (nodes) that are expressed in LNCaP (Methods). A mean of 36% proteins found in these EPINs are encoded by genes differentially expressed in LNCaP compared to LHSAR (an immortalized prostate epithelial line overexpressing androgen receptor; Methods and Table S1).

Each EPIN is composed of three types of nodes (Table S2): proteins with DNA binding motifs in the promoter (promoter-bound nodes), proteins with DNA binding motifs in the enhancers (enhancer-bound nodes), and proteins interacting with promoter-bound or enhancer-bound nodes but without DNA binding motifs onto the promoter or the enhancers (intermediate nodes). Overall, 751 out of the 885 proteins represent intermediate nodes, with 127 of them acting both as intermediate and as DNA-bound nodes in different EPINs. 261 unique DNA-binding proteins have predicted binding sites in at least one of the anchors of enhancers and promoters. A mean of 32.8 (s.d. 11.5) distinct DBPs were identified per promoter anchor with SP1, EGR1, SP2 being the most represented; and a mean of 24.8 (s.d. 7.69) were predicted per enhancer anchor with SP1, IRF1 and TFAP2A being the most represented.

PENGUIN reconstructs an EPIN by grouping enhancers interacting with the same promoter based on the well-known partial redundancy of enhancers ¹⁹. A mean of 1.43 (normalized) promoters (0.88 s.d.) are shared among enhancers, with a maximum of 15 promoters for the same enhancer. To identify communalities and differences among the 4,314 EPINs in PrCa, we performed an unsupervised, hierarchical clustering based on edge composition (Ward’s linkage method, Methods). Using this approach, we identified 8 clusters of promoters with specific networks (Table S1, Table S3, Fig. 2 and Figure S1A).

Protein networks of PrCa risk-associated regulatory elements

We characterized the 8 clusters using PrCa specific annotations. We used the 95% credible set of SNPs (henceforth referred to as PrCa SNPs) across 137 PrCa-associated regions fine-mapped from the largest publicly available GWAS summary statistics [N = 79,148 cases and 61,106 controls ²⁰]. By comparing each cluster with all other clusters, we found a significant enrichment of PrCa SNPs in one specific cluster (cluster 8 or GWAS + cluster; Fisher’s exact test, Methods, Figures S3 B-C). Interestingly this enrichment is exclusively due to SNPs in enhancers (Table 1). Our results show that E-P interactions containing PrCa SNPs are clustered together (red branches in Fig. 2A) indicating that they have similar characteristics in the way their PPI networks are wired.

We found that most edges (67.5%, or 5,550 out of 8,215 edges) are shared among all clusters. We identified the protein interactions that are enriched in each cluster and estimated the significance of overrepresentation of each edge in a cluster compared to all others (Methods). GWAS + cluster (cluster 8 in Fig. 2) exhibits the lowest number of promoters and distinctive network characteristics (Table S3, Figure S1B). Per promoter, it displays the largest number of edges (p-value < 1e-16) and intermediate nodes (p-value < 1e-16), in line with its greater number of enhancers per promoter (p-value < 1e-16), see Figure S2. Moreover, the EPINs of the GWAS + cluster have the lowest values of node-level centrality measures, namely betweenness and degree (Figure S1C). The degree of a node measures the amount of connections it has, while the betweenness centrality measures the amounts of shortest paths that pass through it. Low values of betweenness and degree indicate a lower amount of connections among different components and regions of the network. Betweenness and degree are significantly different across clusters (Kruskall-Wallis test p-value < 1e-16), but not with respect to the ensemble of all EPINs, which indicates that, despite the high number of shared pairwise interactions (67.5% of edges), the wiring of the cluster-specific EPINs are distinctive.

Since CTCF is a major actor in the formation and maintenance of transcriptionally productive E-P interactions ^22,23, we tested the clusters identified by PENGUIN for enrichment in CTCF binding. For this analysis we used CTCF ChIP-seq peaks, from the same cell line (LNCaP), from the ENCODE project instead of predictions based on DNA-binding motifs (Methods). We found that the interactions with CTCF peaks, that we call CTCF+, cluster together (red branches in Fig. 2B, Figure S3A) suggesting that the presence of CTCF in chromatin interactions results in the formation of characteristic PPI networks between the promoter and its enhancers.

CTCF + clusters overlap the GWAS + cluster (Fig. 2, Figures S3A to S3C), also suggesting that CTCF-mediated interactions could be more functionally relevant to PrCa. In particular, GWAS + cluster (representing 6% of the total number of promoters considered) is the only one presenting the unique and significant enrichment in CTCF binding, PrCa SNPs, and oncogenes (Table 2, Table S3, Figures S3D and S3E). This cluster is enriched in the Hippo signaling pathway (KEGG:04390) (Bonferroni-corrected p-value = 0.012) and the Signal transduction pathway (Reactome REAC:R-HSA-162582) with genes such as FOXA1, MYC, FOS (Bonferroni-corrected p-value = 0.047) (Methods, Table S5).

We note that the GWAS + cluster, or any other cluster did not significantly stand out in terms of overall expression level (Figure S4A) or, particularly in terms of fraction of differentially expressed genes (Figure S4B).

We next assessed whether our results are specific to PrCa by comparing the association with SNPs from other diseases. We compared across 17 GWAS from other diseases in the NHGRI GWAS Catalog ²⁴ (Methods). We found that the highest enrichment of SNPs is for GWAS in PrCa and Multiple Sclerosis (Table S6). Interestingly, a possible direct co-morbidity between multiple sclerosis and PrCa has been previously reported ^25,26.

In conclusion, PENGUIN enabled the identification of a cluster of E-P contacts whose EPINs are uniquely enriched in PrCa SNPs, ChIP-seq CTCF binding sites, and oncogenes (a.k.a. GWAS + cluster or cluster 8, Fig. 2 and Table 2). We note that using a different set of GWAS for PrCa, we could still identify cluster 8 as enriched (Methods, Table S6).

Comparison with baseline analytical procedures

Among the 273 promoters belonging to the identified GWAS + cluster (cluster 8 in Fig. 2A), 11 belong to known oncogenes, FOXA1, ZFHX3, CDKN1B, KDM6A, BRCA2, CDH1, CCND1, NKX3-1, BAG4, MYC, GATA2 (Methods). We compared enrichment of PrCa functional annotations in the reconstructed networks with and without inclusion of intermediate proteins. Including intermediates allows increasing the number of retrieved PrCa-related oncogenes in cluster GWAS + from 6 to 11 and increasing significance of enrichment indicating improved specificity (Table S4). We then compared our results with the simple overlap of the genomic regions of E-P contacts and known oncogene promoters. Table S1 reports on the overlaps of E-P contacts with CTCF binding sites (in both enhancers and promoters, see Methods), and PrCa SNPs (in enhancers). Only 30 promoters (12 overlapping the GWAS + cluster) would be identified that overlap both PrCa SNPs and CTCF binding sites. Of these, just 3 are promoters of known oncogenes (and only one, ZMYM3, is not in the GWAS + cluster). We therefore conclude that PENGUIN, and the integration of intermediate PPI network, increases the number of promoters of candidate PrCa-related genes with the additional and unique information on their specific interactome.

Involvement of E-P protein interactomes in tumor-related functional processes

We analyzed the functional enrichment of the set of 885 proteins composing the universe of nodes used to create PPI networks in EPINs. 43 out of these 885 proteins are encoded by one of the 122 known PrCa oncogenes (32 intermediates, 7 DBPs among which MGA, ETV4, ETV1, GATA2, ETV3, ERF, NKX3-1, and 4 of both types among which TP53, MYC, FOXA1, AR; Methods and Table S2). In total, 11 out of 885 have been targeted by PrCa-specific drugs (source: DrugBank; protein targets: ESR2, ESRRA, AR, PARP1, NFKB2, NFKB1, NCOA2, NCOA1, AKT1, TOP2A, TOP2B; drugs: Estramustine, Genistein, Flutamide, Nilutamide, Bicalutamide, Enzalutamide, Olaparib, Custirsen, Amonafide); and 190 out of 885 are targets of non-prostate drugs indicating the possibility of re-purposing.

Considering 500 out of 885 proteins with annotations for KEGG pathways retrieved using g:Profiler ²⁷, 45 were enriched in the prostate cancer pathway (KEGG:05215) (adjusted p-value = 2.54e-27) (Methods and Table S7). We also studied specific protein enrichments in the nodes of the EPINs of each identified cluster (Table S8). Although intermediates are ubiquitous and generally shared among all clusters, we could identify 22 significantly specific proteins enriched in the GWAS + cluster (Methods). Functional enrichment analysis of these 22 proteins revealed significant relationships with tumorigenic processes (Table S9). KEGG Prostate cancer pathway (KEGG:05215) appears highly enriched (adjusted p-value = 4.38e-5) together with other pathways related to tumors such as Colorectal cancer (KEGG:05210, adjusted p-value = 3.45e-9) Pancreatic cancer (KEGG:05212, adjusted p-value = 1.28e-6) and Breast cancer (KEGG:05224, adjusted p-value = 2.05e-6). KEGG pathway KEGG:04919 (Thyroid hormone signaling pathway) appears as the third most enriched pathway (adjusted p-value = 4.02⁻e-8). Thyroid hormones have been previously described as modulators of prostate cancer risk ^28–31. Pathway KEGG:05200 (called Pathways in cancer) appears as the fourth most enriched KEGG concept (adjusted p-value = 1.30e-7). Other classical tumorigenic pathways, such as Wnt signaling pathway (KEGG:04310, adjusted p-value = 3.48e-4) and TGF-beta signaling pathway (KEGG:04350, adjusted p-value = 3.30e-5) appear also enriched. In this regard, recent studies analyzed the involvement of Wnt signaling in the proliferation of prostate cancer cells ^32,33, as well as the involvement and TGF-beta signaling ^34,35.

SNPs path analysis in the E-P protein interactomes

Next, we sought to perform a SNPs analysis of the paths in an EPIN promoter (Methods). In this analysis, a path in a network is a sequence of edges joining a sequence of nodes and going from enhancer to promoter (Fig. 1A). We distinguish between two possible scenarios: (1) PrCa SNPs falling in the DNA binding motifs found in enhancers, indicating a possible dysregulation of TFs binding and activity (Fig. 1B); (2) PrCa SNPs in the genomic regions of the genes that encode for the intermediate nodes of the EPINs, indicating a possible alteration of the PPIs (Fig. 1C). The first analysis aims to identify the location of enhancers that could be targeted by genetic perturbation techniques such as CRISPRi. The second analysis aims to identify the proteins that are potentially affected by mutations so as to enhance our understanding of prostate cancer biology. Overall, we characterized 188 PrCa SNPs falling within any path that connects enancers to promoters.

Network paths with PrCa SNPs in enhancer binding motifs

We sought to detect SNPs located in the DNA binding motifs found in the enhancers of the EPINs. Based on previous evidence ^36,37, our hypothesis is that SNPs in enhancers could disrupt the binding of proteins such as TFs having an impact on their interactome (Fig. 3B).

In Table S10 we list the 36 PrCa SNPs falling within 60 DBP motifs in enhancer regions linking 34 different promoters whose EPINs include 5,184 edges. Among these, we identified 17 PrCa SNPs falling within 16 promoter EPINs (1,894 edges) belonging to the GWAS + cluster that had at least one PrCa SNP in their enhancers. Several of these EPINs promoters were also found differentially expressed (such as DLL1, STOM and SEC11C in the tumor/normal dataset; ID2, RPS27, SEC11C, CASZ1, CRTC2, C5 and STOM in the LNCaP/LHSAR dataset). At the level of intermediate proteins, we also found a group encoded by genes reported to be differentially expressed. We observed that the mean proportion of intermediates that are differentially expressed is on average 40% (Figure S4). We tested whether promoters belonging to the GWAS + cluster were significantly enriched for intermediate protein encoding for differentially expressed genes (Methods). Among the 16 EPINs belonging to the GWAS + cluster that had at least one PrCa SNP in their enhancers, 11 contain expression data to study potential direct effects of the SNPs. In this subset we found 4 EPINs differentially expressed in promoters (3 also differentially expressed in intermediates: CASZ1, ID2, SEC11C), and 4 EPINs only differentially expressed in intermediates: MIIP, MRPL14, MYC, TMEM63B (Table S1). The differential expression of intermediates makes it easier to identify interesting and potentially novel cases. For instance, MYC is not differentially expressed but it has differentially expressed intermediates.

Network paths with PrCa SNPs in the genes coding for intermediate proteins

In this analysis, we identify EPINs with PrCa SNPs falling within genes that encode for intermediate nodes (Table S11), indicating a potential alteration of PPIs involved in E-P contacts. We found that the GWAS + cluster has the highest proportion of PrCa SNPs in the intermediate nodes compared to all other clusters (mean = 53.2, SE = 18.0, p-value < = 0.01, Table S12). The EPINs of STK40 and GATA2 promoters GWAS + cluster display the highest fraction of intermediate proteins with PrCa SNPs in their corresponding genes encoding them (Table S1).

We used the SNP paths to link 172 PrCa SNPs falling within the gene bodies of 26 genes of which 7 are known oncogenes (MAP2K1, CHD3, AR, SETDB1, ATM, CDKN1B, USP28). We identify edges that are most enriched in our GWAS + cluster which could be pointing to essential links between the gene encoding for the intermediate and containing a PrCa predisposing SNP at a particular EPIN. For example, we identified the link between MDM4 containing SNP rs35946963 (PrCa p-value 1e-24) and TP53 ³⁸ and between KDM2A containing SNP rs12790261 (PrCa p-value 1e-7) and BCL6 ³⁹ and ARNT continuing SNP rs139885151 (PrCa p-value 3e-13) and HIF1A ⁴⁰.

We integrated information from pQTL associations between the 172 PrCa SNPs and protein levels (Methods). Two intermediate protein levels (CREB3L4, MAP2K1) were associated with PrCa SNPs falling within the gene encoding for them (p-value of association with protein levels were 7.75e-86 for CREB3L4 and 2.40e-5 for MAP2K1). We identified 3 out of 26 promoter EPINs (TRIM26, MEIS1, POU2F2) with evidence (p-value < 1e-5) of association between the PrCa SNP with the PENGUIN-linked promoter EPIN, pointing to the cancer promoting mechanistic action of these variants: gene with SNPs in POU2F2 linked to the EPIN promoter of gene PHGDH (SNP with lowest p-value rs113631324 = 3.80e-8); gene with SNPs in TRIM26 and EPIN promoter of gene RRM2 (SNP with lowest p-value rs2517606 = 2.69e-7); gene with SNPs in MEIS1 and EPIN promoter of gene STOM (SNP with lowest p-value rs116172829 = 8.19e-6).

A case study: the SNP paths connected to the MYC promoter

Considering HiChIP data, the MYC promoter (chr8:128747814–128748813) is in contact with 73 enhancer regions among which one holds the SNP rs10090154 (p-value of association with PrCa = 1.4e-188). This SNP is located in the binding motif of the transcription factor FOXA1. The integration of PrCa SNPs information highlights paths in the EPIN of MYC that are particularly compelling in the context of the disease (red line in Fig. 4A). The promoter region of MYC binds 8 proteins TFAP2C, KLF5, RBPJ, SP1, ZBTB14, ATF6, ZBTB7A, PRDM1 and contains 17 protein interactors (dots in Fig. 4A) that might be affected by the possible disruption of its binding motif, namely, FOXA1, HMGA1, RCC1, TFAP4, NFIC, PBX1, HOXB9, NFIX, NACC1, RARA, PIAS1, RPA2, H2AFY, RECQL, SATB2, CREB1, AR. The gene encoding for FOXA1 is differentially expressed, along with others of its interactors (Table S10; Methods). Interestingly, 24 correlated PrCa SNPs fall within the genomic region of AR (marked by an asterisk next to the gene name), all with p-value of association with PrCa below 1e-11 (Table S11). AR is targeted by several drugs used in the treatment of prostatic neoplasms, such as apalutamide, bicalutamide, diethylstilbestrol, enzalutamide, flutamide, and nilutamide (triangle in the Fig. 4A, source: DrugBank).

Notably, mutations in FOXA1 enhancers were previously shown to alter TF bindings in primary prostate tumors ³⁷. And, also in line with our observations, FOXA1 enhancer region has been previously reported to be coupled to MYC ⁴¹ and has been shown to have a strong binding of AR ⁴².

Here we used PPI networks to uncover a set of genes implicated in PrCa that could not be identified with other methods. We observed that intermediate nodes of the PPI networks carry intrinsic properties to be exploited for the classification of E-P chromatin loops. We have shown that, without any prior information, our approach is able to group genes according to their implication in PrCa. Most importantly, we found that genes identified in the cluster enriched in SNP associated to PrCa (GWAS + cluster) are candidate oncogenes or partners of oncogenes sharing onco-enhancers.

Our study opens a new path towards the understanding and identification of biological markers in a disease. The method proposed here, PENGUIN, assumes that the PPI network structure between a promoter and its enhancers can be used as a signature to be associated to a specific function or disease. This assumption builds on top of previous works showing the relation between loop 3D topology and chromatin state or expression ⁴³. With PENGUIN we aim to provide a molecular explanation for the observed distinct structural features of loops, and directly relate them to disease. Indeed, we observed that interactions belonging in the GWAS + cluster, but not carrying a GWAS SNP, can still be related to PrCa. We note that our method does not consider the presence of specific SNPs or oncogenes that are only used in the post-processing of data to validate the defined clusters.

Previous methods combined GWAS hits with PPI networks. For example, Ratnakumar and colleagues identified proteins enriched for PPI with GWAS hits ⁴⁴. Recently, Dey and colleagues have demonstrated advantage in using strategies capturing both distal and proximal gene regulation to prioritize genes for disease ⁴⁵. Also, some authors have exploited information from 3D chromatin interactions and GWAS SNPs to relate intergenic SNPs to gene regulation and to cancer ^17,46,47. Differently from these approaches, our method is completely blind to the presence of SNPs and combines information from PPI network and E-P information from enhancer-promoter conformation (H3K27ac-HiChIP) into one framework.

PENGUIN finds clusters of EPINs significantly related to cancer, while classical gene expression analysis did not reveal any significant trend in PrCa. This observation is not in contradiction with our definition of EPIN clusters or even of our core definition of EPIN. In fact, given the level of evidence brought by our analysis, we believe that PENGUIN allows the detection of associations with cancer with a greater sensitivity than a differential expression analysis. Furthermore, our analysis can identify specific links that are disrupted by PrCa-predisposing variants, such as CTCF sites linking a promoter to its enhancers, or intermediate structural proteins involved in E-P network. To further investigate SNPs pathways involved in prostate cancer we provide a web interface at https://penguin.life.bsc.es/.

We note that in this work we employed intermediate nodes through PPI to inform clustering of E-P into networks (intermediate = 1 in Table S4). Yet, PENGUIN’s clustering procedure allows the use of a variety of inputs. First, it can run without grouping the data by promoter into EPINs, considering only single, pairwise, E-P interactions. Alternatively, exploiting our definition of EPINs, PENGUIN allows building networks either: providing the E-P interactions, without information from PPI (intermediate = 0 in Table S4); or only giving the intermediate PPI-network (not considering DBPs for the clustering) without the DNA-binding proteins bound to enhancer or to the promoter (data not shown); or increasing the number of intermediates (e.g. intermediates = 2, data not shown). In any of the above-mentioned configurations PENGUIN is able to yield a clustering that significantly segregates CTCF-enriched and GWAS-enriched EPINs. Here, we used PPI information (intermediates = 1) and observed that the presence of intermediate PPI networks increases segregation significance compared to not using PPI (Table S4).

Our analysis has also caveats. We considered only E-P interactions detected through HiChIP, protein-DNA interactions retrieved from FIMO, and tissue-specific protein-protein interactions available from the integrated interactions database (IID). The completeness of this information depends on the limits of the databases and methods used. Additionally, we considered only focused on protein networks (i.e., only proteins having edges). For sake of visualization purposes, we reduced the number of investigated proteins and study one intermediate only (i.e., 1 edge away). Yet, the clustering is even more pronounced when including three intermediate proteins (data not shown). We should also considere that we used E-P in a stable environment (LNCaP cells), which represents just a snapshot in time. This is acceptable considering that the literature to date supports minimal and quantitative small changes in E-P interactions.

In conclusion, the networks obtained with PENGUIN enable a fine characterization of the molecular associations at play in cancer chromatin and are suitable to train sophisticated machine learning models such as graph neural networks. We were able to prioritize PrCa intermediates that interact in E-P networks in cancer cells (LNCaP) and are amenable to therapeutic intervention. High-throughput functional studies would validate the impact of genetic perturbation of thousands of enhancers at a time. For example, CRISPR-Cas9 technology could allow the targeted editing of specific genomic regions. The PENGUIN approach applied to study PrCa in LNCaP cells can be exploited to investigate any other human disease provided similar data and to study other scenarios (e.g., cell type/GWAS combination) of interest. For example, using an E-P set from another prostate cancer cell line could be exploited to identify target genes regulated by enhancers from different models. These can be prioritized using a genome-wide set of risk SNPs from a disease of interest.

Conformation capture and E-P interactions

We used Hi-C followed by chromatin immunoprecipitation (HiChIP) targeting H3K27Ac in LNCaP cells (androgen-sensitive prostatic carcinoma cell line) across 5 biological replicates including 1 billion reads as previously described ¹⁴. HiChIP is an efficient protein-mediated chromatin-conformation assay ¹³. H3K27Ac is a marker of active enhancers and promoters. Briefly, we used HiC-Pro ⁴⁸ to map HiCHiP reads and extract unique interactions; FitHiChIP ⁴⁹ was used to identify significant interactions with a predefined set of peaks from H3K27ac ChIP-seq in LNCaP to refine accurate anchor ranges. We used q-value < 0.01 and a 5 kb resolution and considered only interactions between 5 kb and 3 Mb. In this analysis, we restricted to a stringent global background estimation to reduce as much as possible the number of false-positive interactions. The corresponding FitHiChIP specifications used were “IntType = 3” (the peak-to-all) for the foreground, meaning at least one anchor to be in the H3K27 peak, and “UseP2PBackgrnd = 1” (the peak-to-peak (stringent)) for the global background estimation of expected counts and contact probabilities for each genomic distance for learning the background and spline fitting. We identified 49,565 significant interactions (FitHiChIP, FDR < 0.01).

We categorized interactions by overlapping anchors with transcription start sites (TSS) and enhancers identified by H3K27ac ChIP-seq as previously described ¹⁴. Briefly, we first extended anchors by 5 kb on either side; we defined promoter regions around the TSS (+/- 500 bases) using RefSeq hg19 (see Data Availability); we defined enhancer regions using 49,638 regions from H3K27ac LNCaP in regular media (union of narrow and broad peaks). Out of the 49,565 significant interactions, we considered only the 24,547 E-P interactions. Specifically, we labeled the promoters and enhancer regions that overlap either right or left anchors and considered E-P if only one anchor overlaps a promoter and the other an enhancer region. The enhancer anchors at this stage of the analysis are of length 15 kb (5 kb resolution of the HiChIP data analysis and additional 5 kb padding added to anchors on either side).

We further prioritized E-P interactions to 1 kb regions and discarded from enhancers the 1 kb bins with fewer HiChIP interactions with the promoter (see E-P HiChIP prioritization section). The 15 kb original E-P interactions dataset contained a mean of 1.6 (1.3 s.d.) promoter anchors per enhancer anchor (after prioritization of enhancer anchor to 1 kb region, mean of 1.4 (0.9 s.d.) promoters per enhancer). There were 11,127 (17,683 prioritized 1 kb regions) enhancer anchors in total; 7,341 (12,385 prioritized 1 kb regions) enhancer anchors are contacted by one promoter anchor with a maximum of 21 promoter anchors (15 using prioritized enhancer regions) sharing the same enhancer.

E-P HiChIP prioritization

In order to reduce experimental artifacts in the context of our EPINs, we developed a specific prioritization method. This prioritization start by normalizing the data assuming, as most used capture-C normalizations (ICE ⁵⁰, Vanilla, or KR ⁵¹) that all biases (e.g. GC content, number of restriction sites, mappability, or in the case of HiChIP, immunoprecipitation bias) can be corrected together. For this normalization step, we assume that there is a specific bias per any 1 kb genomic loci (ꞵ_x for loci x; see Figures S5A and S5B). This bias causes the difference between a theoretical expected number of interactions (E_XY between loci X and Y) and the observed number of interactions (O_XY between loci X and Y). In this representation we can define a system of 9 equations involving three 1 kb loci in the promoter and three 1 kb loci on the enhancer side. This system of equations is then solved using Sequential Quadratic Programming (SQP) ⁵². The procedure is repeated in an overlapping window manner along the 15 kb of the enhancer, alway against the target 1 kb of the promoter and its two 1 kb neighboring loci.

Before the normalization step, we observed a different interaction pattern for interactions below 10 kb (Figure S5C) due, in part, to the contiguity of restriction-enzyme fragments or chromatin persistence length. As these interactions may also be a source of bias in the construction of a PPI network, we removed them from our study.

We applied the normalization to the remaining interactions and observed a better correlation between genomic distance and interaction count (Figures S5D and S5E).

The normalized profile of interactions was finally used to prioritize most interacting 1 kb loci on the 15 kb enhancer (Figure S5F). The selected 1 kb regions are referred to as prioritized enhancer regions.

DNA binding motifs

DNA binding motifs were retrieved from JASPAR, an open-access database of curated, non-redundant binding profiles of DBPs (a.k.a. motifs) stored as position frequency matrices (PFMs). To detect the binding motifs, we used FIMO from the MEME-suite software, with p-value < = 1e-4 and q-value < = 5e-2 cutoffs. JASPAR contains 810 DNA binding motifs of 640 proteins that overlap the E-P contacts identified with HiChIP.

Gene expression data

We assayed RNA sequencing (RNA-seq) in the cell line LNCaP for two replicates using the VIPER pipeline as previously described ¹⁴, and fragments per kilobase of transcript per million mapped reads (FPKM) values were calculated for 20,114 RefSEQ genes. Genes with expression levels above the threshold of 0.003 in both replicates were considered in the entire analysis (Figure S4C).

Protein-protein interaction network

We obtained protein-protein interactions (PPIs) from the Integrated Interactions Database (IID) ⁵³. To better contextualize the interactome information, we combined the annotations of the PPIs from IID database with the LNCaP gene expression data. As for the IID annotations, we applied the following selection criteria. First we selected interactions annotated as “experimental” in the “evidence type” field and identified by at least two independent biological assays reported in the “methods” field. Then, we filtered only for interactions in the prostate or in prostate cancer cells and between nuclear proteins. Finally, we retain proteins whose gene expression levels were FPKM > 0.003 in both replicates (this cut-off removes ~ 30% of the genes). In total, 14,221 proteins from a pool of 20,111 human protein coding genes meet the gene expression criteria. The combination of the above filtering criteria (gene expression, using only nuclear, prostate cancer or prostate and experimentally by 2 methods) resulted in an unweighted network of 31,944 prostate-specific nuclear PPIs among 4,295 proteins.⁵³

PENGUIN pipeline

We set up graph-based approach, called Promoter-ENhancer-GUided Interaction Networks (PENGUIN), to reconstruct individual networks of protein interactions that might occur between one promoter (P) and its contacting enhancers (E), that we call E-P protein-protein Interaction Networks (EPINs). To reconstruct the EPINs, PENGUIN integrates information about chromatin contacts, protein-DNA binding, and protein-protein interactions (PPIs). For the case under study in this work (prostate cancer, PrCa), chromatin contacts information comes from H3K27Ac HiChIP of LNCaP cells (4,314 promoters and 5,789 enhancer regions; see Methods, “Conformation capture and E-P interactions''), protein-DNA binding information ^51,52 comes from the JASPAR database (810 DNA binding motifs of 640 proteins; see Methods, “DNA binding motifs”), and PPIs information comes from the IID database (31,944 prostate-specific nuclear PPIs among 4,295 proteins; see Methods, “Protein-protein interaction network”).

The reconstruction of EPINs follows these steps: for each E-P contact, (1) DNA binding motifs are detected in the corresponding sequences of promoter and enhancer regions; (2) a subnetwork of PPIs is selected containing all promoter-bound proteins, all enhancer-bound proteins, and all their intermediate interactors, with a maximum of 1 intermediate node between enhancer and promoter bound DNA binding proteins; (3) intermediate interactors are discarded if they only connect promoter-bound proteins or enhancer-bound proteins. Using the provided PrCa information, PENGUIN reconstructed 4,314 EPINs consisting of a total of 9,141 PPIs among 885 proteins of which 751 are intermediate proteins linking promoter-bound and enhancer-bound proteins.

Node centrality measures

We employed two measures of node centrality, namely betweenness and degree. Betweenness is a measure of centrality in a graph based on shortest paths. For every pair of nodes in a connected graph, there exists at least one shortest path between the vertices such that the number of edges that the path passes through is minimized. The degree of a node in a network is the number of connections it has to other nodes and the degree distribution is the probability distribution of these degrees over the whole network.

Clustering EPIN

We defined clusters by considering their edge content. We collected the full universe of edges using all existent edges between all promoter EPINs (the union graph) and we performed clustering using a binary representation that encodes this particular edge list. Clustering was performed over the overlap index distance matrix, by calculating Euclidean distance and using Ward’s linkage method. Each leaf in the obtained cluster is a promoter EPIN.

Identification of transcription factors binding directly enhancer to promoter

The TFs binding in the Enhancer region (TFs.E) and the TFs binding in the Promoter region (TFs.P) were identified. The PPIs networks were then searched to see if there is a path linking the (TFs.E) and the (TFs.P).

Identifying enriched clusters using functional annotations

We performed Fisher tests on the obtained clustering for every single branch of the dendrogram. We examined if there is an enrichment of any feature (CTCF, GWAS) for the leaves under the branch of interest compared to those in the rest of the tree. For the CTCF binding feature, we require a CTCF binding (see CTCF ChiP-Seq peaks) to Promoter and at least one of the Enhancer regions of the promoter EPIN. For the GWAS feature, we require the presence/overlap of a GWAS SNP (see Genome-wide association data) in at least one of the Enhancers of the promoter EPIN.

Fisher tests were used to calculate the odds ratio (OR) and enrichment p-values for presence of PrCa annotations within the identified clusters.

Druggability information

We extracted information for target druggability from DrugBank ⁵⁴. The use of each drug was extracted from the Therapeutic Target Database ⁵⁵. We annotated as PrCa druggable each protein node that is target for drugs that are assigned as Approved or under Clinical Trials (Phase 1, 2, 3) or Investigable for Prostate Cancer.

CTCF ChiP-Seq peaks

CTCF ChiP-seq peaks were retrieved from ENCODE project (https://www.encodeproject.org/) for the same Genome assembly, hg19 (Experiment ID: ENCSR315NAC, file ID: ENCFF155SPQ). These narrow peaks were mapped on the enhancer regions using the python package pyranges (see “E-P contacts” section). We found overlaps of the CTCT binding sites with enhancer and promoter anchors allowing a 10 kb gap between them.

PrCa SNPs

To explore enrichment of GWAS across the identified clusters, and to identify the SNP paths, we used the 95% credible set from fine-mapping of the largest PrCa genome-wide association studies (GWAS) collected from Schumacher et al. 2019 (N = 79,148 cases and 61,106 controls), which includes 20,370,946 SNPs, as previously described ⁵⁶. Briefly, we fine-mapped 137 GWAS regions using PAINTOR ⁵⁷, a Bayesian statistical method, with no functional annotations and specifying a maximum of 1 causal SNP. We then constructed a 95% credible set by taking the cumulative sum of the posterior probability until a cumulative 95% posterior probability was reached. This set was composed of 5412 distinct SNPs (rsid). We will refer to these 95% credible set SNPs as PrCa SNPs for brevity. Note that this set also includes SNPs that do not reach genome-wide-filters of p-value significance. We mapped the SNP location to prioritized enhancer regions anchor locations with a window of 10kb. 518 out of 5412 overlap our prioritized enhancer regions; 18 of them overlap our promoter regions. In total 218 prioritized enhancers and 14 promoters overlap a PrCa SNP.

SNP paths (PrCa SNPs in enhancer binding motifs)

A path in a network is a sequence of edges joining a sequence of nodes. We detected PrCa SNPs located in the DNA binding motifs in the enhancers and identified the corresponding SNP paths (linked edges and nodes) for each EPIN promoter. For SNP paths analyses and the web-browser, we used all PrCa SNPs in the 95% credible set. There were 36 PrCa SNPs falling in enhancer binding motifs across clusters 3, 4, 5, 6, 7, 8. To report the most interesting cases in the Tables and Results, we used the subset of those passing genome-wide significance of p-value for PrCa association < 5e-8. There were 15 PrCa SNPs falling in enhancer binding motifs across clusters 3, 5, 6, 7, 8.

SNP paths (PrCa SNPs in intermediate proteins)

We detected PrCa SNPs falling within genes that encode for intermediate nodes, and identified the corresponding SNP paths (linked edges and nodes) for each EPIN promoter. For SNP paths analyses and the web-browser, we used all PrCa SNPs in the 95% credible set. To report the most interesting cases in the Tables and Results, we used the subset of those passing genome-wide significance of p-value for PrCa association < 5e-8.

PrCa GWAS enrichment using GWAS Catalog and comparison with other diseases

This analysis had two aims: 1) explore whether we could replicate our finding and identify the GWAS enriched cluster using a different source for the GWAS; 2) to compare the GWAS signal for different diseases. We estimated enrichment of SNPs overlapping the enhancers in each of the identified clusters by exploring the NHGRI GWAS Catalog associations ²⁴. First, we retrieved GWAS data and filtered the traits according to their “umlsSemanticTypeName” as defined in DisGeNet database ⁵⁸ to one of the following: "Mental or Behavioral Dysfunction", "Neoplastic Process", "Disease or Syndrome", "Congenital Abnormality; Disease or Syndrome", "Disease or Syndrome; Congenital Abnormality", "Disease or Syndrome; Anatomical Abnormality". We considered only traits with at least 10 genome-wide-significant SNPs (unadjusted p-value < 5e-8). We mapped the SNP location to prioritized enhancer anchor locations with a window of 10kb. 104 diseases had SNPs overlaps and 17 of them have more than 10 SNP overlapping (Table S5). For each cluster, we tested enrichment of disease-associated SNPs using Fisher tests and considered significant p-value < 0.01 and OR > 1.

Oncogenes Gene list

We used a previously identified list of 122 Genes ("PrCa_GeneList_Used.csv") known to be somatically mutated in PrCa oncogenesis (37 out of 4,314 promoters considered). As previously described ¹⁴, the 122 oncogenes are a set of prostate cancer–genes curated from three large-scale PrCa studies that show evidence of somatically acquired mutations, at both localized and advanced prostate cancer, known and recurrently altered in localized prostate cancer and metastatic prostate cancer.

Enriched edges within each cluster

Fisher tests were used to compute odds ratios and p-values of the edges and nodes in the eight different clusters. Specifically, each edge or node was tested for presence/absence in a cluster compared to all others. Therefore, one edge or node can be enriched in one or more than one cluster, it cannot be enriched in all clusters.

Enriched intermediate nodes within each cluster

We computed protein importance for each cluster in terms of two network centrality measures: betweenness and degree. For each protein we obtain both betweenness and degree specificity ratios in order to equitably quantify internal protein centrality differences between the clusters. For each of the found clusters we independently estimated the specificity of the observed protein centrality measures (“Betweenness” and “Degree”). For a given protein (P_i) in a particular cluster (C_j), we define the specificity as the ratio between the mean centrality value of P_i inside the fraction of networks belonging to C_j ; divided by the mean centrality value of P_i for the fraction of networks outside of the cluster C_j.

Specificity ratio (P_i, C_j) = (mean (Pi centrality in C_j networks) + 1) / (mean (P_i centrality in non C_j networks) + 1)

We assessed protein specificity ratio significance for each cluster upon random network cluster generation. Aiming to assess the significance of the different specificity ratios for the proteins within each cluster, we developed a significance analysis test based on random cluster subsamplings. In order to compute the significance of a given protein specificity ratio (P_i) within a particular cluster of analysis (C_j), we performed 1000 random network samplings to produce random network clusters containing the same amount of networks as the real cluster being analyzed (i.e. if the real cluster contains 100 networks, the random clusters generated will contain 100 random networks out of the 4,314 clustered networks). Within each of those 1000 random clusters, we compute the corresponding protein specificity ratios, with the p-value representing the probability of finding the protein specificity ratio to be higher or equal to the real value computed for the particular cluster of interest (C_j).

We also performed Fisher tests to assess enrichment for the presence of the node in the cluster (Fisher test p-value < 0.01). EP300 was excluded from the enrichment test as the presence of that node was not significantly enriched (Fisher test p-value < 0.01). 22 proteins (SMAD2, KAT5, NCOR2, MAPK8, SMAD4, CREBBP, CTNNB1, PGR, HDAC3, HDAC2, GSK3B, UBA52, UBE2I, JUND, PIAS1, XRCC5, CDK6, XRCC6, MAPK1, FOS, HIF1A and MAPK3) were found to be significantly specific for both betweenness and degree ratios (p-value < 0.01 for both centrality measures and over-represented in this cluster using Fisher tests).

Functional gene set enrichment analysis

Functional enrichment analysis was performed using the g:GOST module from g:Profiler, a web tool to perform simultaneous gene set enrichment analysis across multiple biomedical databases ²⁷. g:GOST performs cumulative hypergeometric tests of an input geneset against preprocessed database-specific gene sets. We run the web service considering only annotated genes for the statistical domain scope. Reported adjusted p-values correspond to Benjamini-Hochberg correction for multiple testing, with adjusted p-values ≤ 0.05 considered to be significant. Gene set enrichment analysis results are provided for KEGG pathways, Reactome, Gene Ontology, Wikipathways, TRANSFAC, miRTarBase, Human Protein Atlas, CORUM and Human Phenotype Ontology.

For the enrichment analysis of significantly specific proteins of the GWAS + cluster, we provided as input the 22 previously described proteins. For the enrichment analysis of the GWAS + cluster, we provided as input all genes associated with the EPIN promoters in cluster GWAS+.

Differential Gene Expression

We integrated data from EPIN promoters with differential gene expression (DE) from two sources.

DE analysis on prostate cancer tumor versus normal was downloaded from GEPIA: http://gepia2.cancer-pku.cn/#degenes, which use the TCGA and GTEx projects databases to compare gene expression between tumor and normal tissues under Limma, both under and over expressed. We used the default thresholds of log2FC of 1 and qvalue cut-off of 0.01. These data covered 84 out of 885 genes encoding for intermediates in PENGUIN and 413 out of 4,314 promoter EPINs.

DE analysis of RNA-Seq on LHSAR (an immortalized prostate epithelial line overexpressing androgen receptor) versus LNCaP was performed as previously described. Briefly, RNA-seq data were processed using the VIPER pipeline ⁵⁹. Reads were aligned to the hg19 human genome built with STAR. FPKM values were calculated with Cufflinks for 20,114 RefSEQ genes included in the VIPER repository. Differential expression analysis was performed with the DESeq2 R package ⁶⁰. 15,650 genes with DE data covered 884 of the 885 genes encoding for intermediates in PENGUIN and 3,286 genes out of 4,314 promoter EPINs.

We annotated whether the EPIN promoters themselves and the genes encoding the intermediate proteins in our data were DE using either of the two databases. We considered as DE those genes passing |log2 fold change| > 1 and adjusted p-value < = 0.01. For the LNCAP/LHSAR dataset, we could compute a Fisher test of enrichment of differentially expressed genes encoding for intermediate proteins within each EPIN promoter versus within the SNP paths (we could not compute this for the GEPIA since we did not have the full dataset of covered genes). The genes that were not passing these filters were considered non-DE and the genes not covered by the two datasets were excluded from the enrichment analysis described next. For each EPIN we calculated the fraction of DE intermediates within the SNP paths and we estimated the enrichment of those compared to the fraction of DE intermediates in the full EPIN network.

To find the enrichment of DE genes in SNP paths (PrCa SNPs in intermediate proteins) compared to those in the entire EPIN, we computed as enrichment the ratio of Fraction1 / Fraction2, where:

Fraction1 = (number of DE intermediates within SNP paths) / (number of covered intermediates within SNP paths), and

Fraction2 = (number of DE intermediates the EPIN) / (number of covered intermediates in the EPIN).

We report the EPIN genes passing enrichment (“enrichment_DE_deseq_SNP.bs.TF.path”) > 1.

pQTL look-up

We downloaded summary statistics with genome-wide association between SNPs and 4907 proteins reported in the deCODE study ⁶¹ and annotated with pQTL association the SNPs we identified falling in either in enhancer binding sites or in node genomic locations. The deCODE pQTL summary statistics data contained information on 4,907 proteins and 186 (201 PrCa SNPs out of the 213 PrCa SNPs we looked up were in the data and 186 also matched by alleles). 808 out of the 4,314 genes promoters ("Gene_network") and 278 out of the 885 gene intermediates (in total 997 out of 4,918 genes promoters and coding for intermediates in our networks) have information on associations with their respective coded proteins covered by the pQTL deCODE data.

Acknowledgments

The authors are grateful to José María Fernández González (Barcelona Supercomputing Center) for the crucial guidance with the PENGUIN web server. They also thank the Biola Javierre’s lab at the Josep Carreras Leukaemia Research Institute for the support, the ‘RNA initiative’ at IIT and all the members of Tartaglia’s lab at CRG, Sapienza and IIT.

FUNDING

The research leading to these results has been supported by the European Research Council [RIBOMYLOME_309545 and ASTRA_855923] and the H2020 projects [IASIS_727658 and INFORE_825080]. CG has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 754490 – MINDED project.

Data Availability

RefSeq hg19 from UCSC Genome Browser is available at the following URL:

http://genome.ucsc.edu/cgi-bin/hgTables?hgsid=694977049_xUU5i1QkIJ50dj5miBt9wkAYuxN3&clade=mammal&org=&db=hg19&hgta_
group=genes&hgta_track=knownGene&hgta_table=knownGene&hgta_regionType=genome&position=&hgta_outputType=
selectedFields&hgta_outFileName=knownGene.gtf

Ensembl hg19 data for overlaps of SNPs with intermediates: biomart / ensembl from bioaRt package

All EPINs and related statistics can be downloaded through the PENGUIN web service at https://penguin.life.bsc.es/

Rebello, R. J. et al. Prostate cancer. Nat Rev Dis Primers 7, 9 (2021).
Tan, M. H. E., Li, J., Xu, H. E., Melcher, K. & Yong, E.-L. Androgen receptor: structure, role in prostate cancer and drug discovery. Acta Pharmacol. Sin. 36, 3–23 (2015).
Cancer Genome Atlas Research Network. The Molecular Taxonomy of Primary Prostate Cancer. Cell 163, 1011–1025 (2015).
de Bono, J. et al. Olaparib for Metastatic Castration-Resistant Prostate Cancer. N. Engl. J. Med. 382, 2091–2102 (2020).
Adams, E. J. et al. FOXA1 mutations alter pioneering activity, differentiation and prostate cancer phenotypes. Nature 571, 408–412 (2019).
Parolia, A. et al. Distinct structural classes of activating FOXA1 alterations in advanced prostate cancer. Nature 571, 413–418 (2019).
Ahmadiyeh, N. et al. 8q24 prostate, breast, and colon cancer risk loci show tissue-specific long-range interaction with MYC. Proc. Natl. Acad. Sci. U. S. A. 107, 9742–9746 (2010).
Dekker, J. & Misteli, T. Long-Range Chromatin Interactions. Cold Spring Harb. Perspect. Biol. 7, a019356 (2015).
Norton, H. K. & Phillips-Cremins, J. E. Crossed wires: 3D genome misfolding in human disease. J. Cell Biol. 216, 3441–3452 (2017).
Krumm, A. & Duan, Z. Understanding the 3D genome: Emerging impacts on human disease. Semin. Cell Dev. Biol. 90, 62–77 (2019).
Sur, I. & Taipale, J. The role of enhancers in cancer. Nat. Rev. Cancer 16, 483–493 (2016).
Chen, X., Ma, Q., Shang, Z. & Niu, Y. Super-enhancer in prostate cancer: transcriptional disorders and therapeutic targets. NPJ Precis Oncol 4, 31 (2020).
Mumbach, M. R. et al. HiChIP: efficient and sensitive analysis of protein-directed genome architecture. Nat. Methods 13, 919–922 (2016).
Giambartolomei, C. et al. H3K27ac HiChIP in prostate cell lines identifies risk genes for prostate cancer susceptibility. Am. J. Hum. Genet. 108, 2284–2300 (2021).
Zhang, K., Li, N., Ainsworth, R. I. & Wang, W. Systematic identification of protein combinations mediating chromatin looping. Nat. Commun. 7, 12249 (2016).
Wang, R. et al. Hierarchical cooperation of transcription factors from integration analysis of DNA sequences, ChIP-Seq and ChIA-PET data. BMC Genomics 20, 296 (2019).
Liu, N. et al. Seeing the forest through the trees: prioritising potentially functional interactions from Hi-C. Epigenetics Chromatin 14, 41 (2021).
Deng, W. & Blobel, G. A. Manipulating nuclear architecture. Curr. Opin. Genet. Dev. 25, 1–7 (2014).
Kvon, E. Z., Waymack, R., Gad, M. & Wunderlich, Z. Enhancer redundancy in development and disease. Nature Reviews Genetics vol. 22 324–336 Preprint at https://doi.org/10.1038/s41576-020-00311-x (2021).
Schumacher, F. R. et al. Author Correction: Association analyses of more than 140,000 men identify 63 new prostate cancer susceptibility loci. Nat. Genet. 51, 363 (2019).
Huerta-Cepas, J., Serra, F. & Bork, P. ETE 3: Reconstruction, Analysis, and Visualization of Phylogenomic Data. Mol. Biol. Evol. 33, 1635–1638 (2016).
Zuin, J. et al. Cohesin and CTCF differentially affect chromatin architecture and gene expression in human cells. Proc. Natl. Acad. Sci. U. S. A. 111, 996–1001 (2014).
Pugacheva, E. M. et al. CTCF mediates chromatin looping via N-terminal domain-dependent cohesin retention. Proc. Natl. Acad. Sci. U. S. A. 117, 2020–2031 (2020).
Buniello, A. et al. The NHGRI-EBI GWAS Catalog of published genome-wide association studies, targeted arrays and summary statistics 2019. Nucleic Acids Res. 47, D1005–D1012 (2019).
Bosco-Lévy, P. et al. Incidence and risk of cancer among multiple sclerosis patients: A matched population-based cohort study. Eur. J. Neurol. 29, 1091–1099 (2022).
Schwartz, G. G. Multiple sclerosis and prostate cancer: what do their similar geographies suggest? Neuroepidemiology 11, 244–254 (1992).
Raudvere, U. et al. g:Profiler: a web server for functional enrichment analysis and conversions of gene lists (2019 update). Nucleic Acids Res. 47, W191–W198 (2019).
Mondul, A. M. et al. Circulating thyroxine, thyroid-stimulating hormone, and hypothyroid status and the risk of prostate cancer. PLoS One 7, e47730 (2012).
Hsieh, M.-L. & Juang, H.-H. Cell growth effects of triiodothyronine and expression of thyroid hormone receptor in prostate carcinoma cells. J. Androl. 26, 422–428 (2005).
Lehrer, S., Diamond, E. J., Stone, N. N. & Stock, R. G. Serum thyroid-stimulating hormone is elevated in men with Gleason 8 prostate cancer. BJU Int. 96, 328–329 (2005).
Hellevik, A. I. et al. Thyroid function and cancer risk: a prospective population study. Cancer Epidemiol. Biomarkers Prev. 18, 570–574 (2009).
Ma, F. et al. Autocrine canonical Wnt signaling primes noncanonical signaling through ROR1 in metastatic castration-resistant prostate cancer. Cancer Res. (2022) doi:10.1158/0008-5472.CAN-21-1807.
Wei, X. et al. Paracrine Wnt signaling is necessary for prostate epithelial proliferation. Prostate 82, 517–530 (2022).
Natani, S. et al. Activation of TGF-β - SMAD2 signaling by IL-6 drives neuroendocrine differentiation of prostate cancer through p38MAPK. Cell. Signal. 91, 110240 (2022).
Xi, X. et al. High expression of small nucleolar RNA host gene 3 predicts poor prognosis and promotes bone metastasis in prostate cancer by activating transforming growth factor-beta signaling. Bioengineered 13, 1895–1907 (2022).
Speedy, H. E. et al. Insight into genetic predisposition to chronic lymphocytic leukemia from integrative epigenomics. Nat. Commun. 10, 3615 (2019).
Zhou, S. et al. Noncoding mutations target cis-regulatory elements of the FOXA1 plexus in prostate cancer. Nat. Commun. 11, 441 (2020).
Mejía-Hernández, J. O. et al. Targeting MDM4 as a Novel Therapeutic Approach in Prostate Cancer Independent of p53 Status. Cancers 14, (2022).
Liu, L., Liu, J. & Lin, Q. Histone demethylase KDM2A: Biological functions and clinical values (Review). Exp. Ther. Med. 22, 723 (2021).
Mandl, M. & Depping, R. ARNT is a potential direct HIF-1 target gene in human Hep3B hepatocellular carcinoma cells. Cancer Cell Int. 17, 77 (2017).
Sur, I., Tuupanen, S., Whitington, T., Aaltonen, L. A. & Taipale, J. Lessons from functional analysis of genome-wide association studies. Cancer Res. 73, 4180–4184 (2013).
Jia, L. et al. Functional enhancers at the gene-poor 8q24 cancer-linked locus. PLoS Genet. 5, e1000597 (2009).
Galan, S., Serra, F. & Marti-Renom, M. A. Identification of chromatin loops from Hi-C interaction matrices by CTCF-CTCF topology classification. NAR Genom Bioinform 4, lqac021 (2022).
Ratnakumar, A., Weinhold, N., Mar, J. C. & Riaz, N. Protein-Protein interactions uncover candidate ‘core genes’ within omnigenic disease networks. PLoS Genet. 16, e1008903 (2020).
Dey, K. K. et al. SNP-to-gene linking strategies reveal contributions of enhancer-related and candidate master-regulator genes to autoimmune disease. Cell Genom 2, (2022).
Javierre, B. M. et al. Lineage-Specific Genome Architecture Links Enhancers and Non-coding Disease Variants to Target Gene Promoters. Cell 167, 1369–1384.e19 (2016).
López de Maturana, E. et al. A multilayered post-GWAS assessment on genetic susceptibility to pancreatic cancer. Genome Med. 13, 15 (2021).
Servant, N. et al. HiC-Pro: an optimized and flexible pipeline for Hi-C data processing. Genome Biol. 16, 259 (2015).
Bhattacharyya, S., Chandra, V., Vijayanand, P. & Ay, F. Identification of significant chromatin contacts from HiChIP data by FitHiChIP. Nat. Commun. 10, 4221 (2019).
Imakaev, M. et al. Iterative correction of Hi-C data reveals hallmarks of chromosome organization. Nat. Methods 9, 999–1003 (2012).
Rao, S. S. P. et al. A 3D map of the human genome at kilobase resolution reveals principles of chromatin looping. Cell 159, 1665–1680 (2014).
Virtanen, P. et al. Author Correction: SciPy 1.0: fundamental algorithms for scientific computing in Python. Nat. Methods 17, 352 (2020).
Kotlyar, M., Pastrello, C., Sheahan, N. & Jurisica, I. Integrated interactions database: tissue-specific view of the human and model organism interactomes. Nucleic Acids Res. 44, D536–41 (2016).
Wishart, D. S. et al. DrugBank 5.0: a major update to the DrugBank database for 2018. Nucleic Acids Res. 46, D1074–D1082 (2018).
Zhou, Y. et al. Therapeutic target database update 2022: facilitating drug discovery with enriched comparative data of targeted agents. Nucleic Acids Res. 50, D1398–D1407 (2022).
Schumacher, F. R. et al. Association analyses of more than 140,000 men identify 63 new prostate cancer susceptibility loci. Nat. Genet. 50, 928–936 (2018).
Kichaev, G. et al. Integrating functional data to prioritize causal variants in statistical fine-mapping studies. PLoS Genet. 10, e1004722 (2014).
Piñero, J. et al. The DisGeNET knowledge platform for disease genomics: 2019 update. Nucleic Acids Res. 48, D845–D855 (2020).
Cornwell, M. et al. VIPER: Visualization Pipeline for RNA-seq, a Snakemake workflow for efficient and complete RNA-seq analysis. BMC Bioinformatics 19, 135 (2018).
Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550 (2014).
Ferkingstad, E. et al. Large-scale integration of the plasma proteome with genetics and disease. Nat. Genet. 53, 1712–1721 (2021).

Table 1. Enrichment of PrCa SNPs in cluster 8 (GWAS+) when considering SNPs overlapping enhancers, promoters, either or both.

PrCa SNPs overlaps	Odds Ratio (OR)	p-value
Only enhancers	11.329	1.80e-12
Only promoters	1.139	0.600
Either enhancers or promoter	8.551	2.68e-11
Both enhancers and promoter	0	1

Table 2. Enrichment of PrCa SNPs, CTCF ChIP-seq binding sites (“CTCF” in the header), and other PrCa annotations (oncogene promoters and PrCa SNPs from GWAS Catalog) across the eight clusters identified by PENGUIN. Cluster 8 is enriched in CTCF binding, PrCa SNPs, and oncogenes. Symbols code: +, enriched; -, depleted; =, neutral. OR: Fisher’s exact test Odds Ratio.

Cluster	Number of genes	CTCF	OR CTCF	P-value CTCF	PrCa SNPs	OR PrCa SNPs	P-value PrCa SNPs	Number of oncogene promoters	OR oncogenes	P-value oncogenes
1	825	-	0.617	1.91e-9	-	0.28	2.46e-2	8	1.17	0.67
2	399	-	0.613	9.65e-6	-	0.00	2.00e-2	5	1.54	0.38
3	544	+	1.348	1.35e-3	=	0.80	8.27e-1	2	0.39	0.31
4	491	=	1.084	4.09e-1	=	0.51	3.60e-1	4	0.94	1.00
5	465	=	0.841	9.12e-2	=	0.75	8.14e-1	1	0.23	0.17
6	641	-	0.664	4.24e-6	=	0.38	1.03e-1	1	0.16	0.03
7	676	+	1.655	2.12e-9	=	1.42	3.18e-1	5	0.84	1.00
8	273	+	3.287	3.64e-20	+	11.33	1.80e-12	11	6.48	1.04e-5

There is NO Competing Interest.

3972520supp7061944rkyzhz.xlsx
Supplementary Tables
SupplementaryFigures.docx

Download PDF

Journal Publication

published 06 Dec, 2023

Read the published version in Nature Communications →

Version 1

posted

You are reading this latest preprint version

Reconstructing protein interactions at enhancer-promoter regions in prostate cancer

Status:

Journal Publication

Version 1

Abstract

Figures

Introduction

Results

PENGUIN identifies clusters of protein interactions based on chromatin contacts

Protein networks of PrCa risk-associated regulatory elements

Comparison with baseline analytical procedures

Involvement of E-P protein interactomes in tumor-related functional processes

SNPs path analysis in the E-P protein interactomes

Network paths with PrCa SNPs in enhancer binding motifs

Network paths with PrCa SNPs in the genes coding for intermediate proteins

A case study: the SNP paths connected to the MYC promoter

Discussion

Methods

Declarations

Acknowledgments

FUNDING

Data Availability

References

Tables

Additional Declarations

Supplementary Files

Status:

Journal Publication

Version 1