Constructing methylation-driven ceRNA networks unveil tumor heterogeneity and predict patient prognosis

Abstract

Cancer development involves a complex interplay between genetic and epigenetic factors, with emerging evidence highlighting the pivotal role of competitive endogenous RNA (ceRNA) networks in regulating gene expression. However, the influence of ceRNA networks by aberrant DNA methylation remains incompletely understood. In our study, we proposed DMceNet, a computational method to characterize the effects of DNA methylation on ceRNA regulatory mechanisms and apply it across eight prevalent cancers. By integrating methylation and transcriptomic data, we constructed methylation-driven ceRNA networks and identified a dominant role of lncRNAs within these networks in two key ways: (i) 17 cancer-shared differential methylation lncRNAs (DMlncs), including PVT1 and CASC2, form a Common Cancer Network (CCN) affecting key pathways such as the G2/M checkpoint, and (ii) 24 cancer-specific DMlncs construct unique ceRNA networks for each cancer type. For instance, in LUAD and STAD, hypomethylation drives DMlncs like PCAT6 and MINCR, disrupting the Wnt signaling pathway and apoptosis. We further investigated the characteristics of these methylation-driven ceRNA networks at the cellular level, revealing how methylation-driven dysregulation varies across distinct cell populations within the tumor microenvironment. Our findings also demonstrate the prognostic potential of cancer-specific ceRNA relationships, highlighting their relevance in predicting patient survival outcomes. This integrated transcriptomic and epigenomic analysis provides new insights into cancer biology and regulatory mechanisms.

Graphical Abstract

competitive endogenous RNA (ceRNA), DNA methylation, long noncoding RNA, cancer

Introduction

Tumor heterogeneity, a formidable challenge in cancer research, reflects the inherent diversity in genetic and phenotypic characteristics both within and among different tumors [1]. This intricate variability arises from complex interactions between genetic and epigenetic factors, particularly those involving long non-coding RNAs (lncRNAs) and DNA methylation [2]. These interactions ultimately result in the abnormal expression of oncogenes and tumor suppressor genes [3]. Notably, DNA methylation at cytosine-phosphate-guanine (CpG) dinucleotide sites has been widely reported to regulate gene expression, maintain genome stability, and influence various cellular processes, including growth, proliferation, differentiation, and apoptosis [4]. Concurrently, lncRNAs involved in gene regulation function as competing endogenous RNA (ceRNA), exerting significant influence on multiple facets of cancer biology [5].

Numerous studies have identified significant roles for lncRNAs via the ceRNA mechanism in various cancer processes. For instance, DLEU2L competes with BRCA2 for binding to miR-210-3, inhibiting the proliferation, invasion, and migration of pancreatic cancer cells [6, 7]. Additionally, HOXD-AS1 binds competitively with miR-130a-3p, preventing miRNA-mediated degradation of SOX4, thereby activating the expression of EZH2 and MMP2 and promoting hepatocellular carcinoma (HCC) metastasis [8]. Moreover, the expression of lncRNAs is directly influenced by DNA methylation, and some lncRNAs can regulate DNA methylation in cancer by recruiting chromatin modifiers [9]. For example, the expression of the hypomethylation-driven lncRNA HUMT is upregulated in the promoter region, promoting tumor proliferation and lymph node metastasis in triple-negative breast cancer (TNBC) by activating FOXK1 [10]. Furthermore, lncRNAs can influence DNA methylation status by controlling the accumulation of S-adenosylhomocysteine (SAH), thereby regulating the binding of DNA methyltransferases to target gene promoter regions [11]. These findings underscore the combined role of lncRNAs and DNA methylation in cancer. Therefore, a systematic analysis of the characteristics of DNA methylation-driven mechanisms of lncRNA-mediated ceRNA in cancer is urgently needed.

In this study, we proposed a computational method to comprehensively construct methylation-driven ceRNA networks, leveraging the integration of transcriptome and epigenome data to provide a thorough analysis of cancer mechanisms. We characterize the common and cancer-specific methylation-driven ceRNA networks at tissue and cellular levels, revealing the contribution of aberrant methylation-driven ceRNA mechanisms to cancer dysfunction and highlighting the heterogeneity among cancers. This research enhances our understanding of the complex networks that drive cancer and paves the way for targeted treatments and personalized medicine.

Results

Constructing methylation-driven ceRNA regulatory networks across cancers with DMceNet

To characterize the effects of DNA methylation on competitive endogenous RNA (ceRNA) regulatory mechanisms in various cancers, we proposed a computational method, DMceNet (Fig. 1a). DMceNet constructs methylation-driven ceRNA networks by: (i) assessing the methylation levels and dividing samples into hypo- and hypermethylation groups, (ii) identifying differential methylation lncRNAs and mRNAs (DMlncs and DMms) based on differential expression, and (iii) calculating correlations between DMlncs and DMms, and compiling ceRNA pairs that share miRNA regulation.

Figure 1

Identification and analysis of methylation-driven ceRNAs. (a) Framework for DMceNet. (b) Upset plots display the overlap of DMlnc (top) and DMm (bottom), with the bar plot illustrating the count of specific and recurrent DMlnc and DMm. (c) Bar plots depict identified DMlnc and DMm across cancers; the point on the diagonal indicate cancer-specific cases and others denote recurrence across multiple cancers. (d) Heatmap shows methylation level of DMlnc across eight cancers. (e) ceRNA networks are cancer-specific. Light rhomboid nodes represent hypomethylated DMlnc, dark rhomboid nodes represent hypermethylated DMm, with light and dark circular nodes indicating hypo- and hypermethylated DMlnc, respectively. Node size correlates with degree.

Employing the DMceNet across eight prevalent cancer types, we identified 200-400 DMlncs and ~1000 DMms per cancer, with significant recurrence across multiple cancers. A total of 2198 DMms and 486 DMlncs were found in multiple cancers, with 164 DMms and 26 DMlncs common to all eight, while 786 DMms and 298 DMlncs were cancer-specific. THCA and KIRP showed the highest numbers of cancer-specific DMms and DMlncs (134 and 64, respectively; Table 1 and Fig. 1b). STAD had significant overlap with other cancers (Jaccard coefficients of 0.313 for DMlncs and 0.398 for DMms), while HNSC and CESC had the highest molecular similarity, likely due to their squamous epithelial origin (Fig. 1c). Additionally, the average methylation levels of DMlncs across cancers revealed consistent hypomethylation patterns, particularly for DMlncs like CAS2 and PVT1, suggesting their key roles in cancer dysregulation (Fig. 1d).

Table 1

Number of differential methylation lncRNAs and mRNAs in eight cancers.

Cancer	Abbreviation	Differentially expressed molecules		Molecules in ceRNA network
		DMm	DMlnc	DMm	DMlnc
Cervical squamous cell carcinoma and endocervical adenocarcinoma	CESC	1260	216	293	18
Uterine Corpus Endometrial Carcinoma	UCEC	1230	297	221	30
Kidney renal papillary cell carcinoma	KIRP	1011	298	331	32
Head and Neck squamous cell carcinoma	HNSC	1483	334	706	35
Lung adenocarcinoma	LUAD	1366	286	223	27
Thyroid carcinoma	THCA	988	203	239	21
Colon adenocarcinoma	COAD	1017	189	274	29
Stomach adenocarcinoma	STAD	1551	297	373	33

Importantly, we constructed eight DNA-methylation-driven ceRNA networks and assessed the involvement of DMlncs and DMms in ceRNA mechanisms across various cancers, all exhibiting a scale-free structure with power-law exponents (α) consistently above 2, peaking at 2.87 (Fig. S1). Approximately 25% of DMms and 10% of DMlncs formed ceRNA networks, with HNSC showing the most extensive involvement (705 DMms and 35 DMlncs), indicating a complex regulatory landscape (Fig. 1e). In conclusion, DMceNet effectively unveiled methylation-driven ceRNA regulatory relationships across eight major cancers, identifying numerous DMlncs and DMms with both shared and cancer-specific patterns. The consistent hypomethylation of certain DMlncs underscores their potential roles in cancer dysregulation, offering valuable insights into cancer development.

Differential methylation lncRNA dominates the complex methylation-driven ceRNA networks

To characterize the methylation-driven ceRNA networks identified across various cancers, we first explore the nodes within the networks. Although differential methylation lncRNAs (DMlncs) represented less than 10% of the nodes, they exhibited significantly higher degrees than differential methylation mRNAs (DMms) across all eight cancers (Wilcoxon test; Fig. 2a and b). Nodes with higher degrees, particularly the top 10%, were identified as hubs. DMlnc hubs outnumbered DMm hubs across multiple cancers, highlighting their potential as central players in the complex ceRNA networks (Fig. 2c). Given their centrality, we investigated the functional roles and molecular mechanisms of DMlnc hubs. DMlnc hubs exhibited lower methylation levels compared to other DMlncs (Fig. 2d), suggesting higher expression levels. This elevated expression implies a crucial role for DMlnc hubs in cancer processes, making them important targets for further investigation.

Figure 2

Patterns of methylation-driven DMlnc in cancer. (a) Bar plots display the number of DMlnc/DMm in ceRNA networks across cancer types. (b) Boxplots compare the degree of hub DMlnc and hub DMm, ^**P < 0.01; ^***P < 0.0001 with t-tests. (c) Bar plots show the number of hub DMlnc/DMm in ceRNA network across cancer types. (d) Violin plots illustrate the methylation level of hub DMlnc vs. other DMlnc. ^*P < 0.05; ^**P < 0.01; ^***P < 0.0001. (e) Violin plots show the length of DMlnc within and outside ceRNA network across cancers. ^*P < 0.05; ^**P < 0.01; ^***P < 0.0001. (f) Circos diagram shows DMlnc-DMm relationships in HNSC, with edges representing inter- and intra-chromosomal interactions. (g) Bar plots show the number of DMlnc-DMm pairs on the same or different chromosomes. (h) Heatmaps depict ceRNA network similarity across cancer types, with the Jaccard index representing network similarity.

We then examined the influence of lncRNA length on the ceRNA regulatory mechanism by comparing the transcript length of DMlncs involved in ceRNA networks with those not involved. Except for CESC, DMlnc in the networks had significantly longer transcripts than those not involved (Wilcoxon test, P-value < 0.05), indicating that longer lncRNAs are more likely to function as molecular sponges in gene regulation (Fig. 2e). Additionally, we analyzed the chromosomal distribution of DMlnc-DMm competing pairs (Figs 2f, g, and S2). While a small subset of these pairs shared the same chromosomes (marked in red), most were located on different chromosomes (marked in gray). This finding aligns with a previous study by Zhang et al. [12], supporting the idea that methylation-influenced lncRNA-mRNA competing pairs often operate through distal regulation. Jaccard index analysis revealed low similarity between ceRNA networks across cancers, indicating distinct methylation-driven regulatory mechanisms in each cancer (Fig. 2h). In addition, DMlncs were more frequently shared across cancers than DMms, suggesting a more conserved role for DMlncs across cancers, consistent with their inherent characteristics.

In conclusion, our study reveals key characteristics of methylation-driven ceRNA networks, highlighting the dominant role of hypomethylated, long DMlncs in cancer regulation. These findings enhance our understanding of ceRNA dynamics and their influence on cancer progression.

Pan-cancer ceRNA network analyses reveal cancer-common perturbed pathways

Cancer development often shares common characteristics, such as uncontrolled cell division and growth [13, 14]. We investigate the influence of conserved methylation-driven ceRNA mechanisms constituted by differential methylation lncRNAs (DMlncs) shared across different cancers on common dysfunction mechanisms. We identified 17 cancer-shared DMlncs, including well-known cancer-related ncRNAs like PVT1 and CASC2, based on their presence as hubs in four or more cancers (Fig. 3a). The proportion of these common hubs varied across cancers, ranging from 70% in CESC to 28% in KIRP. These DMlncs and their associated mRNAs formed a comprehensive Common Cancer Network (CCN) (Fig. 3b and Tables 3 and S1).

Figure 3

Regulatory mechanisms of cancer common DMlnc. (a) Heat map showing the identified DMlnc in cancers. Each row represents a cancer type and each column represents a DMlnc, the DMlnc identified as a hub in the cancer ceRNA network was labeled. (b) Cancer common network (CCN), constructed with cancer-shared DMlnc. (c) Enrichment analysis for genes regulated by DMlnc in the CCN related to Hallmark geneset, or KEGG pathways. (d) Common network with known cancer-related genes in CCN, as CCN2. (e) Cancer common regulation mechanism involving PVT1 and CASC2. (f) Top and bottom: Upset plot showing the overlap of DMm and Hallmarker function across various cancers, respectively. Middle: Sankey plot showing the DMm and Hallmarker influencing by PVT1 across various cancers.

Table 3

Hyper- and hypo-methylation lncRNA mediated ceRNA pairs.

Cancer	Common		Specific
	hypermethylation	hypomethylation	hypermethylation	hypomethylation
CESC	6	192	0	3
COAD	21	332	8	0
HNSC	82	709	521	100
KIRP	0	411	79	16
LUAD	39	171	0	30
STAD	50	418	0	6
THCA	0	167	43	0
UCEC	8	113	27	14

To explore the common dysfunction mechanism influenced by these DMlncs, we performed hallmark [15] and KEGG pathway enrichment analysis on their regulated genes (Methods, Fig. 3c). Key pathways, including the G2/M checkpoint, cell cycle, p53 signaling pathway, and Hippo signaling pathway were dysregulated across cancers. These DMlncs also regulated genes like MYC oncogenes and TGFb1, crucial in cancer onset and progression. Among the cancer-shared DMlncs, we discoed six well-established cancer-related lncRNAs: CASC2, SNHG7, PVT1, TTN-AS1, RHPT1-AS1, and LEF1-AS1. A subnetwork (CCN2) with these six lncRNAs as core components revealed their significant impact on multiple cancers (Fig. 3d). Within this network, CASC2, and PVT1 were identified for their wide effect in multiple cancers, and their interaction genes were enriched in cancer-related pathways. CASC2, driven by hypomethylation in seven cancer types, competes with PIK3CA, FGFR2, and MET, forming ceRNA pairs (Fig. 3e bottom). These oncogenes are enriched in cancer-associated pathways like carbon metabolism and EGFR tyrosine kinase inhibitor resistance, which are associated with cancer progression and poor prognosis [16]. PVT1, also driven by hypomethylation in six cancers, competes with several miRNAs to regulate TRAF5, TRAF1, and MYC, all of which are associated with critical cancer-related pathways, including the NF-κB and TNF signaling pathways (Fig. 3e top). TRAF5, TRAF1, and MYC are crucial for prognosis, survival signaling, and cell proliferation, respectively [17–19]. In short, PVT1 upregulation due to hypomethylation enhances miRNA competition, disrupting key pathways and driving cancer progression. These findings highlight the potential of PVT1 and CASC2 as therapeutic targets.

To understand the diversity in regulated patterns of shared DMlncs, we analyzed their impact across different cancers (Figs 3f and S3). While DMlncs recur across cancers, the affected differential methylation mRNAs (DMms) vary significantly. For instance, PVT1 influences different DMms across cancers, with the highest disparity in HNSC (44 unique DMms), followed by STAD and COAD. Four pathways, including reactive oxygen species, were uniquely enriched in HNSC, while one and three pathways were uniquely affected in STAD and COAD, respectively. However, nine pathways, including G2M checkpoint and E2F target, were affected by PVT1 in most cancers. This suggests that while cancer-shared lncRNA, such as PVT1, engages in ceRNA mechanisms with different DMms in diverse cancers, it ultimately impacts similar functional mechanisms (Fig. 3f).

In summary, we identified 17 shared DMlncs, including PVT1 and CASC2, forming a Common Cancer Network that affects key pathways such as G2/M checkpoint, p53 signaling, and Hippo signaling. Abnormal methylation of these DMlncs drives mRNA dysregulation, leading to cancer development. Hypomethylation of CASC2 and PVT1 disrupts key signaling pathways, underscoring their therapeutic potential.

Cancer-specific ceRNA networks reveal tumor heterogeneity

To explore the regulated patterns of differential methylation lncRNAs (DMlncs) in specific cancers, we identified 24 cancer-specific DMlncs, including eight proven cancer-related lncRNAs, and constructed cancer-specific methylation-driven ceRNA networks uniquely associated with particular cancers (Tables 2 and 3).

Table 2

Cancer-specific differential methylation lncRNAs. Proven cancer-related lncRNAs are marked in bold black.

Cancer	DMlnc in specific cancer
CESC	EPHA1-AS1
UCEC	FLJ37453, LINC02878, LINC01001, AL162586.1
KIRP	PSMG3-AS1, ZFAND2A − DT, AC008040.1, MIR3936HG, AC092535.2
HNSC	LINC00339, FEZF1-AS1, DLG1-AS1, CTBP1-DT, LINC01703, TM4SF19-AS1, LNCOC1, AP000866.1
LUAD	PCAT6, TBX5-AS1
THCA	PINK1-AS, HEIH
COAD	LINC01806
STAD	MINCR

In THCA, two significant lncRNAs, HEIH and PINK1-AS1, were identified as potential cancer biomarkers. DNA hypermethylation drives the downregulation of HEIH and PINK1-AS1, which reduces the ability to compete for miRNA binding. For example, PINK1-AS1 forms ceRNA interactions with 14 oncogenes (purple nodes in the network) through 6 miRNAs, while HEIH interacts only with PDE7B by competitively binding hsa-miR-429. As lncRNAs are hypermethylated and expression decreases, the efficiency as ceRNAs diminishes. This reduction inhibits the expression level of related oncogenes. Then it leads to dysfunction of cell apoptosis and nucleosome assembly, finally causing the occurrence of cancer (Fig. 4a). In KIRP, PSMG3-AS1, driven by hypermethylation, similarly experiences downregulated expression. PSMG3-AS1 predominantly regulates oncogenes via interactions with hsa-miR-210-3p and hsa-miR-143-3p, and decreased expression due to hypermethylation reduces its effectiveness as a ceRNA. This dysregulation impacts critical pathways such as the G2M checkpoint and TGF-β signaling pathway (Fig. 4b). In HNSC, we identified seven oncogenes, such as CASP8 and LAMC1, forming ceRNA relationships with the hypermethylation-driven DMlnc CTBP1-DT via interactions with 22 miRNAs. Another ceRNA relationship involves the oncogene PKP4 and the DMlnc LINC00339 through hsa-miR-218-5p. These oncogenes significantly impact cancer-related functions. LAMC1 is linked to immune cell infiltration [20], CASP8 is an apoptosis-related cysteine peptidase relevant to HNSC risk [21], and PKP4 regulates cell adhesion [22]. DNA hypermethylation drives the reduced expression of CTBP1-DT and LINC00339, impairing the ability to compete for miRNA binding. This reduction leads to decreased expression of genes associated with cancer, resulting in dysfunctions such as aberrant cell proliferation and immune infiltration. This finding also suggests a potential association of CTBP1-DT with head and neck squamous cell carcinoma (Fig. 4c). In COAD, LINC01806 competed with a colon oncogene CASP8 for binding hsa-miR-145a-3p and hsa-miR-194-5p. CASP8 is integral to apoptosis and closely relevant to colon cancer development [23]. The hypermethylation of LINC01806 leads to its dysregulation expression, which in turn disrupts the regulation of CASP8. This disruption impairs crucial cellular processes like apoptosis, ultimately contributing to colon cancer development (Fig. 4d).

Figure 4

Regulations of cancer-specific DMlnc in ceRNA network. The marked circular nodes in the network represent oncogenes related to the given cancer, and the diamond node represent cancer-related DMlnc.

In contrast, in cancers like LUAD and STAD, hypomethylation upregulates lncRNAs such as PCAT6 and MINCR. Hypomethylation enhances the expression and ceRNA efficiency of lncRNAs, promoting oncogenic pathways. For instance, in LUAD, the hypomethylation-driven DMlnc PCAT6 competes with the oncogene TRAF2 through hsa-miR-330-5p to form a ceRNA relationship, which disrupts necrotic apoptosis [24]. Another ceRNA relationship involves the lung cancer gene APC, and the lncRNA TBX5-AS1, affecting the Wnt signaling pathway and apoptosis [25]. DNA hypomethylation in lung cancer drives the expression of PCAT6 and TBX5-AS1, resulting in abnormal expression of associated oncogene and disrupting the Wnt signaling pathway and apoptosis. (Fig. 4e). In STAD, hypomethylation upregulates the DMlnc MINCR, enhancing its ability to compete for miRNA binding and upregulate oncogenes TIMM17A and SET. This disrupts the normal expression of oncogenes, contributing to STAD development (Fig. 4f). Although no studies have confirmed an association with cancer, the specific DMlncs identified in UCEC and CESC may provide directions for follow-up studies (Fig. S4).

In summary, hypermethylation tends to downregulate lncRNAs, reducing the ceRNA efficiency and impairing the ability to suppress oncogene expression, thereby promoting cancer pathways. Conversely, hypomethylation upregulates lncRNAs, enhancing their ceRNA function and facilitating oncogenic processes. These findings suggest that the methylation status of DMlncs significantly impacts their role in cancer pathways, either by decreasing or increasing the efficiency as ceRNAs.

Cellular DNA methylation-driven ceRNA networks display unique characteristics

We extended our analysis to characterize methylation-driven ceRNA mechanisms across cell types. In tumor cells, ceRNA networks driven by DNA methylation are primarily involved in promoting oncogenic signaling, cell proliferation, and survival pathways. For example, in COAD, tumor cells exhibited strong interactions within ceRNA networks that regulate genes involved in cell cycle progression and metastasis. In contrast, cancer-associated fibroblasts (CAFs) are pivotal in shaping the TME by supporting tumor growth, invasion, and immune evasion. In CESC, UCEC, and KIRP, methylation-driven ceRNA networks were more active in CAFs. These interactions primarily regulate extracellular matrix remodeling, angiogenesis, and secretion of pro-tumorigenic factors, contributing to creating a pro-tumorigenic microenvironment (Fig. S5). In HNSC and THCA, methylation-driven ceRNA networks were evenly distributed across cell types, indicating complex tumor microenvironment interactions. Tumor cells, cancer-associated fibroblasts, and T cells were the primary interacting cells within these networks. This indicates that the methylation-driven ceRNA network regulates various aspects of cancer formation, development, and immunity (Fig. S5, and Table S2.).

HNSC, with the most intricate methylation-driven ceRNA networks, was chosen for in-depth analysis. We curated single-cell transcriptome data from the GEO database and annotated 13 cell types using canonical markers (Fig. 5a and b): Tumor cells (KRT7 and KRT17), Salivary cells (STATH), Cancer-associated fibroblasts (CAFs) (COL1A2 and MMP2), Myofibroblasts (ACTA2), Endothelial cells (PECAM1), T cells (PTPRC and CD3E), B cells (CD79A), Natural killer (NK) cells (NKG7 and XCL2), Plasma cells (IGHG1), Mast cells (TPSAB1), Mature dendritic cells (DCs) (LAMP3), plasmacytoid dendritic cells (pDCs) (LILRA4), Tumor-associated macrophages (TAMs) (CD163) [26]. We delineated the methylation-driven ceRNA landscape at single-cell resolution by examining DMlnc-DMm pair expression across cell types (Fig. 5c and d). In HNSC, apart from a subset of DMlnc-DMm pairs that exhibited high expression across most cells, the expression of DMlnc-DMm pairs in different cell types displayed pronounced specificity. For instance, FGD5-AS1 predominantly formed ceRNA interactions with HSPG2 in endothelial cells, whereas PTAM predominantly formed ceRNA interactions in tumor and immune cells, consistent with MAPK5’s features. We found that PVT1, a lncRNA associated with cancer and driven by hypomethylation, predominantly regulated the ribosome-related gene RPL12 in tumor cells and T cells, with stronger expression observed in T cells (Fig. 5e). This suggests a potential role for PVT1 in promoting immune evasion of tumor cells by upregulating RPL12 expression in HNSC. In summary, our mapping of the methylation-driven ceRNA network at the single-cell level, combined with the analysis of its cellular origins, revealed distinct regulatory characteristics and mechanisms of action across different cell types.

Figure 5

Characterization of ceRNA networks driven by methylation at single-cell resolution for HNSC. (a) t-SNE representation for HNSC. Color-coded for cell types. (b) Bubble heatmap showing expression levels of canonical cell markers for various cell types. (c) Heatmap showing expression levels of DMlnc-DMm pairs identified from bulk RNA-seq data. (d) Methylation-driven ceRNA networks at HNSC. Pie charts show the corresponding cell type. (e) Expression levels of specific DMlnc-DMm pairs. Color coded for cell types.

Methylation-driven ceRNAs predict clinical prognosis

Based on the methylation level of the cancer-specific ceRNA network, we performed survival analysis across multiple cancers to delve into their direct impact on cancer prognosis. By examining the methylation levels within specific ceRNA networks, we uncovered compelling evidence in COAD, KIRP, and THCA, where ceRNA interactions tied to distinct molecular mechanisms were available in effectively dividing cancer samples into high-risk and low-risk groups. This stratification had significant clinical implications, as high-risk groups were associated with notably poorer prognostic outcomes (log-rank test, P < 0.05) (Fig. 6a). Moreover, one remarkable aspect of our findings was the superiority of ceRNA-based risk prediction over single-gene or single-lncRNA approaches. Such as COAD, hypermethylation driven LINC01806 competitive binding CASP8, and disruption impairs crucial cellular processes like apoptosis. This methylation-driven ceRNA pair significance influences the survival risk (P < 0.05) compared to lncRNA (P = 0.44) or mRNA (P = 0.17) alone. Traditional risk regression analysis relying on individual mRNA or lncRNA was insufficient for accurately distinguishing the survival status of cancer samples. In contrast, using the complex ceRNA interactions as a comprehensive risk prediction strategy resulted in accurate differentiation of patients’ prognostic outcomes. This highlights the strength and promise of the risk score obtained from ceRNA networks as a dependable and robust prognostic factor in cancer (Fig. 6b).

Figure 6

Survival analysis of ceRNA relationships. (a) Kaplan–Meier plot for overall survival of cancer (COAD, KIRP, THCA) patients stratified by methylation of ceRNA relationships. (b) Results of multivariable cox regression model for the interactions and single gene in ceRNA.

These findings emphasize the intricate nature of cancer regulatory networks and underscore the clinical significance of their disruption. The ceRNA interactions, especially those intricately connected to specific molecular processes, emerge as promising candidates for the creation of prognostic tools within cancer research and clinical application. Their capacity to decode the complex landscape of cancer prognosis represents a substantial advancement in comprehending and potentially addressing this intricate disease.

Discussion

lncRNAs function as molecular sponges in competitive endogenous RNA (ceRNA) networks, sequestering miRNAs and regulating mRNA expression, playing a key role in disease etiology [27]. Additionally, DNA methylation plays a crucial role in shaping the expression of both coding and non-coding genomic elements, especially in complex diseases like cancer [2]. In this study, we integrated transcriptomic and epigenomic data to explore the impact of DNA methylation on ceRNA dysregulation mechanisms in cancer. Our approach sheds light on the complex mechanisms underlying cancer development and highlights the extensive heterogeneity across different types of cancers.

We have proposed a computational method for constructing methylation-driven ceRNA networks and systematically applied it across eight distinct cancers. Our investigations revealed that, despite minor variations in network complexity among cancer types, the regulatory mechanisms driven by differential methylation lncRNAs (DMlncs) were remarkably consistent. DNA hypomethylation primarily drives lncRNAs with long transcripts to form ceRNA mechanisms with distal mRNAs, thereby regulating cancer development. Given the pivotal role of DMlncs, our primary research focus gravitated towards understanding their mechanisms, particularly in conjunction with well-established cancer-related lncRNAs. Additionally, we explore the methylation-driven ceRNA network within cells, elucidating its unique characteristics in different cell types. Unlike traditional ceRNA studies that primarily focus on broad miRNA-lncRNA-mRNA interactions within the transcriptome [28, 29], our approach integrates both epigenome and transcriptome data to dig deeper into the underlying mechanisms by which methylation drives ceRNA dysfunction. Our method uncovers the mechanisms by which epigenetic modifications drive the expression of lncRNAs and mRNAs, thereby shaping ceRNA networks in cancer. This deeper, multi-omics perspective offers a more comprehensive understanding of regulatory dynamics by linking epigenetic modifications, like methylation, to ceRNA regulation. Furthermore, by integrating molecular regulatory information data from five databases, we minimize false positives and enhance the robustness of our results, laying a strong foundation for further experimental validation.

However, our analysis represents only the beginning of uncovering the many regulatory mechanisms within these complex networks. Our findings suggest that the regulatory mechanisms driven by DMlncs are remarkably consistent across cancer types and highlight further exploration to evaluate the framework’s applicability beyond these specific cancers. Future studies could explore the generalizability of this framework to other cancers and its potential relevance to non-cancerous diseases where epigenetic regulation is known to play a significant role, such as neurological disorders, cardiovascular diseases, and autoimmune conditions. Moreover, as research on cancer-related lncRNAs progresses and single-cell DNA methylation sequencing improves in quality and quantity, we expect our findings to provide a more comprehensive and detailed understanding of these regulatory landscapes. This ongoing exploration promises to deepen our understanding of cancer biology and ultimately enhance our ability to develop targeted therapeutic strategies.

In this study, we embarked on a comprehensive analysis of cancer ceRNA mechanisms by ingeniously integrating transcriptomic and epigenomic datasets. Our pioneering exploration into the role of DNA methylation and ceRNA has not only illuminated previously uncharted territories in cancer biology, but has also set the stage for deeper dissection of the intricate molecular machinery underlying tumorigenesis. Our study offers a new perspective on integrating additional epigenetic modifiers and multi-omics data. Utilizing advanced approaches, these interdisciplinary studies have the potential to revolutionize our understanding of cancer etiology and facilitate the development of highly personalized treatment strategies based on individual molecular profiles.

Materials and methods

Collection and preprocessing of datasets

Transcriptome expression data

Genome-wide mRNA and lncRNA expression profiles across 8 cancer types were obtained from The Cancer Genome Atlas Program (TCGA). This includes cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC), uterine corpus endometrial carcinoma (UCEC), kidney renal papillary cell carcinoma (KIRP), head and neck squamous cell carcinoma (HNSC), lung adenocarcinoma (LUAD), thyroid carcinoma (THCA), colon adenocarcinoma (COAD), stomach adenocarcinoma (STAD). Genes expressed in at least 70% of samples were selected. After imputing missing values using the k-nearest neighbor method, expression profiles were log2-transformed for data regularization. Using GENCODE (v36), we classified lncRNAs and protein-coding genes.

DNA methylation data

450 k methylation microarray data for the same eight cancer types were obtained from TCGA and processed using the ChAMP package. We filtered data with (i) detection p-values greater than 0.01, (ii) non-CpG sites, (iii) all SNP-related, (iv) mapped to multiple loci, and (v) probes on the X chromosomes and Y chromosome, etc. Low-quality samples were removed, and only samples containing both DNA methylation profiles and gene expression profiles were selected for subsequent analysis.

miRNA–mRNA/lncRNA interaction pairs

To ensure comprehensive and reliable miRNA targets, we integrated data from lncACTdb [30], starBase [31], miRSponge [32], miRTarBase [33] and LncRNASNP [34] databases. This integration resulted in 9282 miRNA-lncRNA interaction pairs and 375 865 miRNA-mRNA interaction pairs.

Cancer-related lncRNA and mRNA

Known cancer-related lncRNAs were obtained from LncRNADisease [35] and Lnc2Cancer [36] databases. Cancer-related mRNA data were obtained from the Network of Cancer Gene (http://ncg.kcl.ac.uk/index.php).

DMceNet: Construction of methylation-driven ceRNA network across cancers

We developed a computational method to construct DNA methylation-driven ceRNA networks by integrating DNA methylation and gene expression data. We hypothesize that significantly differentially expressed mRNA/lncRNA in different methylation groups are driven by methylation. By analyzing methylation levels in the promoter regions, we divide the samples into hypo- and hypermethylation groups, extract highly expressed genes, and define methylation-driven ceRNA pairs through shared miRNA and significant co-expression.

Identification of DNA methylation-driven lncRNA and mRNA. For each gene, promoter regions were defined as the ±2 kb regions around transcription start sites (TSS), and the average β value of all methylation sites mapped to its promoter region was calculated to determine methylation level.

$$ \beta =\frac{M}{M+U} $$

Where U represents unmethylated signal strength and M represents methylated signal strength.

Samples were ranked based on their methylation levels, with the top 40% classified as hypermethylated and the bottom 40% as hypomethylated. Significantly differentially expressed lncRNAs and mRNAs were extracted and labeled as DMlnc and DMm, respectively (t-test). We focused on genes highly expressed in the hypomethylated group, reflecting the negative regulatory relationship between methylation and expression (log2FC > 0 and FDR < 0.05).

Identification of DMlnc and DMm shared targeting. We collected lncRNA-miRNA and mRNA-miRNA regulatory information from five molecular regulation databases, retaining pairs present in more than two databases. After integration, we obtained 3 468 717 lncRNA-mRNA relationship pairs that share at least one miRNA target. DMlnc-DMm pairs with shared miRNA targets were selected for further analysis.

Calculation of the Pearson correlation between DMlncs and DMms. For each DMlnc (⁠|$X$|⁠) and DMm (⁠|$Y$|⁠), the Pearson correlation coefficient |$r$| was calculated, and Student’s t-test was used for statistical comparison.

$$ r=\frac{\sum_{i=1}^n\left({X}_i-\overline{X}\right)\left({Y}_i-\overline{Y}\right)}{\sqrt{\sum_{i=1}^n{\left({X}_i-\overline{X}\right)}^2}\sqrt{\sum_{i=1}^n{\left({Y}_i-\overline{Y}\right)}^2}} $$

A DMlnc-DMm pair was defined as a methylation-driven ceRNA pair if it met two conditions: (1) shared targeting by at least one miRNA and (2) significant co-expression (Pearson correlation coefficient > 0.3 and P-value < 0.05) within cancer. These pairs formed the basis of the methylation-driven ceRNA network, visualized using Cytoscape.

Network similarity

To compare ceRNA networks across cancer types, the Jaccard index for shared DMlncs and DMms was calculated, quantifying the similarity between networks. The Jaccard index is defined as:

$$ J\left(A,B\right)=\frac{N_A\cap{N}_B}{N_A\cup{N}_B} $$

Where A and B represent two networks. |${N}_A$| and |${N}_B$| represent the number of DMlnc or DMm in the networks A and B, respectively.

Functional enrichment analysis

We used hallmark gene sets from MSigDB (https://www.gsea-msigdb.org/gsea/msigdb), Gene Ontology (GO) terms [37, 38], and KEGG pathways [39] for functional enrichment analysis. Hallmark gene sets were used to identify functional characteristics for cancer-shared nodes in the ceRNA network, while KEGG pathways characterized ceRNA mechanisms. GO terms were primarily used for other functional analyses. Enrichment analysis was performed using the cumulative hypergeometric distribution.

$$ P=1-\sum_{k=0}^m\frac{\left(\genfrac{}{}{0pt}{}{n}{k}\right)\left(\genfrac{}{}{0pt}{}{N-n}{M-k}\right)}{\left(\genfrac{}{}{0pt}{}{N}{M}\right)} $$

Where N represents the total number of genes in the genome, M is the number of genes included in a specific pathway/term/hallmark, n represents the number of mRNAs interacting with a specific lncRNA, and m is the number of mRNAs interacting with the lncRNA and annotated to the respective pathway/term/hallmark.

Single-cell RNA-seq data processing across cancer types

The single-cell RNA-seq data for eight different cancer types were obtained from the Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) and ArrayExpress (https://www.ebi.ac.uk/biostudies/arrayexpress) (Table 4), The count matrix for each cancer type was downloaded and subjected to data preprocessing and standardization using the R package “Seurat” v4.2.0. Cells with fewer than 200 detected genes or with more than 15% mitochondrial content were removed for quality control. Following filtering, the functions “NormalizeData” and “ScaleData” were applied to normalize gene expression across cells, with a scale factor of 10 000. The top 2000 variable genes were identified using the “FindVariable” function, and principal component analysis was performed on them, the first 20 principal components were used with the functions “FindNeighbors” and “FindClusters” for clustering. The clusters were annotated using canonical cell markers obtained from recent research [26, 40–46] and the Cellmarker 2.0 database [47]. Visualization was achieved using t-distributed stochastic neighbor embedding (t-SNE). Subsequently, after normalizing the transcriptome expression profiles, we mapped each DMlnc-DMm pair across different cancers based on its expression levels.

Table 4

Single-cell RNA-seq data resources.

Cancer	Accession No	Resources	PMID
CESC	E-MTAB-11948	Arrayexpress	37427448
COAD	GSE161277	GEO	35221332
HNSC	GSE188737	GEO	36973261
KIRP	GSE152938	GEO	34722263
LUAD	GSE189357	GEO	36434043
STAD	GSE183904	GEO	34642171
THCA	GSE191288	GEO	35515000
UCEC	GSE173682	GEO	34739872

Survival analysis

Clinical data for eight cancer types were obtained from TCGA. Patients were divided into two groups based on risk scores, calculated as:

$$ R\mathrm{isk}\ \mathrm{Score}={\sum}_{i=1}^n{r}_i Me(i) $$

where n represents the number of DMlnc/DMm in cancer-specific molecular mechanism, |${r}_i$| denotes the Cox regression coefficient of the DMlnc/DMm, and |$Me(i)$| represents the methylation level of the DMlnc/DMm in the patient.

Survival analysis was conducted using the Cox regression model, and survival curves were drawn. Patients were categorized into high-risk and low-risk groups based on mean risk scores, and survival differences were assessed using the Kaplan–Meier estimation method and log-rank test.

Statistical analysis

Unless stated otherwise, all statistical analyses and presentations were conducted using R Project (version 4.2.1), and statistical significance was defined as a P-value less than 0.05. The significance of differences between the two groups was determined by the Wilcoxon rank sum test. Bar plots, dot plots, and box plots were produced using the R base package “ggplot2”. The circus plots were generated using OmicStudio (https://www.omicstudio.cn/tool) Networks were visualized by Cytoscape (version 3.9.1).

Author contributions

Conceptualization: Jing Bai; Data curation: Xinyu Li, Chuo Peng, Hongyu Liu; Formal analysis: Xinyu Li, Hongyu Liu; Funding acquisition: Jing Bai; Investigation: Xinyu Li, Chuo Peng, Hongyu Liu; Methodology: Xinyu Li, Hongyu Liu, Weixin Liang; Project administration: Jing Bai; Resources: Jing Bai; Software: Xinyu Li, Chuo Peng, Hongyu Liu; Supervision: Jing Bai; Validation: Mingjie Dong, Shujuan Li, Weixin Liang; Visualization: Xinyu Li, Mingjie Dong, Shujuan Li, Weixin Liang; Writing—original draft: Hongyu Liu, Xinyu Li; Writing—review & editing: Jing Bai, Xia Li.

Conflict of interest statement: The authors declared that they have no competing interests.

Funding

This work was supported by the STI2030-Major Projects [2021ZD0202400]; National Natural Science Foundation of China [U23A20166, 62102123, 32070673]; Hainan Provincial Natural Science Foundation [822RC698]; Heilongjiang Touyan Innovation Team Program.

Data availability

All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. All unique identifiers, or web links for publicly available datasets are described in the paper. All codes are available upon reasonable request, and the major codes for DMceNet are available on GitHub (https://github.com/superlaoyuzi/DMceNet).

Author notes

Xinyu Li, Chuo Peng and Hongyu Liu contributed equally to this work.