Effects of Differentially Methylated CpG Sites in Enhancer and Promoter Regions on the Chromatin Structures of Target LncRNAs in Breast Cancer

Abstract

Aberrant DNA methylation plays a crucial role in breast cancer progression by regulating gene expression. However, the regulatory pattern of DNA methylation in long noncoding RNAs (lncRNAs) for breast cancer remains unclear. In this study, we integrated gene expression, DNA methylation, and clinical data from breast cancer patients included in The Cancer Genome Atlas (TCGA) database. We examined DNA methylation distribution across various lncRNA categories, revealing distinct methylation characteristics. Through genome-wide correlation analysis, we identified the CpG sites located in lncRNAs and the distally associated CpG sites of lncRNAs. Functional genome enrichment analysis, conducted through the integration of ENCODE ChIP-seq data, revealed that differentially methylated CpG sites (DMCs) in lncRNAs were mostly located in promoter regions, while distally associated DMCs primarily acted on enhancer regions. By integrating Hi-C data, we found that DMCs in enhancer and promoter regions were closely associated with the changes in three-dimensional chromatin structures by affecting the formation of enhancer-promoter loops. Furthermore, through Cox regression analysis and three machine learning models, we identified 11 key methylation-driven lncRNAs (DIO3OS, ELOVL2-AS1, MIAT, LINC00536, C9orf163, AC105398.1, LINC02178, MILIP, HID1-AS1, KCNH1-IT1, and TMEM220-AS1) that were associated with the survival of breast cancer patients and constructed a prognostic risk scoring model, which demonstrated strong prognostic performance. These findings enhance our understanding of DNA methylation's role in lncRNA regulation in breast cancer and provide potential biomarkers for diagnosis.

Keywords: DNA methylation; Hi-C; breast cancer; enhancer; lncRNA; promoter.

Figure 1

Probe annotation and methylation features of lncRNAs. (A) Number of DNA methylation probes and lncRNAs for five lncRNA categories. (B) Expression levels of five lncRNA categories across tumor and normal samples. (C) Global methylation value distribution of five lncRNA categories across tumor and normal samples. (D) DNA methylation patterns of protein-coding genes, lncRNAs, and miRNAs across the gene body and ±5 kb flanking regions of the gene body. (E) DNA methylation patterns of five different categories of lncRNAs across the gene body and ±5 kb flanking the gene body. Analyses utilized TCGA-BRCA DNA methylation and gene expression data. Statistical significance was assessed using the t-test (* p < 0.05, ** p < 0.01, *** p < 0.001).

Figure 2

Genomic location of proximal and distally associated CpGs of lncRNAs according to ChromHMM and TF binding regions (A) Bar plot showing the enrichment of distally associated CpGs of lncRNAs across functional regulatory regions based on MCF-7 ChromHMM annotation. The bars represent the ratio of observed to expected frequencies for distal CpGs. (B) Enrichment of the proximal CpGs of lncRNAs in functional regulatory regions. (C) Functional annotation of the lncRNAs with DMCs in the relevant enhancer regions. Top 20 clusters with their representative enriched terms (colored by cluster ID), where nodes that share the same cluster ID are typically close to each other. (D) Network plot of the lncRNAs with DMCs in promoter regions. Data obtained from MCF-7 ChromHMM annotation and Metascape analysis.

Figure 3

Assessment of prognostic values in the TCGA-BRCA, GSE20711, and GSE20685 datasets. Comparison of the OS status of breast cancer patients with varying risk scores, and K–M curves in (A) TCGA-BRCA, (B) GSE20711, and (C) GSE20685 datasets. The AUC values for the time-dependent ROC curves depict the OS prediction values for the (D) TCGA-BRCA, (E) GSE20711, and (F) GSE20685 datasets. Statistical significance was determined using the log-rank test for K–M curves and the concordance index (C-index) for AUC calculations.

Figure 4

Determination of the independent prognostic value of risk model (A) A signature-based nomogram was applied to estimate 1-, 3-, and 5-year overall survival probabilities.(* p < 0.05, ** p < 0.01, *** p < 0.001). (B) Calibration plots of the nomogram for predicting the 1-, 3-, and 5-year OS. (C) The C-index was utilized to assess and compare the prognostic accuracy between clinical factors and the risk score. (D) Subgroup analysis of the Kaplan–Meier (K–M) survival curves was conducted using the log-rank test based on factors such as age, histopathological grade, and clinical stage. All results were based on the TCGA-BRCA dataset.

Figure 5

Effect of DMCs on enhancer–promoter looping and DP-lncRNA TMEM220-AS1 expression. Hi-C contact maps for the regions of chromosome 17 (10.224 Mb–11.224 Mb) at 10 kb resolution in the (A) MCF-7 cell lines and (B) HMEC cell lines; purple points indicate chromatin loops in the HMEC and MCF-7 cell lines, respectively; red rectangles represent TMTM22–AS1; green rectangles indicate significant CpG sites. (C) Genome browser snapshots of TMEM220-AS1 in the MCF-7 cell lines, the purple vertical bars highlight the important loop anchor regions which are associated with DP-lncRNAs. The promoter of TMEM220-AS1 co-localize with the anchor of important chromatin loop, DNase-seq peaks, and histone (H3K27ac, H3K4me1) ChIP-seq peaks, and (D) genome browser snapshots of TMEM220-AS1 in the HMEC cell lines. Hi-C data were obtained from the ENCODE database, and snapshots were generated using the Washu Epigenome Browser.

Figure 6

DNA methylation affects 3D genome architecture through distal enhancers. (A) Hi-C contact maps of C9orf163 in the MCF-7 cell lines and (B) HMEC cell lines. (C) Hi-C contact maps of HID1-AS1 in the MCF-7 cell lines and (D) HMEC cell lines. Hi-C data were obtained from the ENCODE database. (Purple point: chromatin loops; Red: lncRNAs; Green: significantly different CpG sites, Blue: enhancer regions.).

Figure 7

Relationships between immune cell infiltration, the TME, and the key methylation-driven lncRNA signature. (A) Heat map of immune responses among the high- and low-risk groups based on the CIBERSORT, ESTIMATE, and ssGSEA algorithms. (B) Comparison of TME scores in both risk groups via the ESTIMATE algorithm. (C) Box plot comparing 13 immune-linked functions in both risk groups. (D) The CIBERSORT algorithm was used to quantify the distribution of 22 tumor-infiltrating immune cells in all HNSCC patients. (E) Violin plot showing the fraction of 22 immune cells in both risk groups. Statistical significance was assessed using the t-test (* p < 0.05, ** p < 0.01, and *** p < 0.001).

Figure 8

Classification of lncRNAs based on their genomic proximity to neighboring transcripts [50].