We downloaded the scRNA-seq data from Gene Expression Omnibus (GEO) (accession number: GSE107858). Following the filtering process described in the paper (3), we performed further analysis on 84 MCF-7 cells. RNAs with fragments per kilo base of transcript per million reads mapped (FPKM) >1 in at least 25% (84 × 0.25 = 21) of the cells were used for further analysis. The filtering parameter is referred to the paper (9). The dataset GSE107863 for T47D was an independent validation cohort to support the findings obtained using MCF-7 cells.

The ceRNA network was downloaded (10) (http://starbase.sysu.edu.cn/) using the Web API. The ceRNAs for all mRNAs, lncRNAs, and pseudogenes were downloaded using default parameters. The ceRNA network contained 308,266 ceRNA pairs composed of 18,942 RNAs. The complete information is presented in Table S1 .

We obtained the known cancer-related genes from the Cancer Gene Census (CGC) database (http://cancer.sanger.ac.uk/cosmic/census), which contains 576 genes (11). Other 876 cancer-related genes were also downloaded from the Genetic Association Database (GAD) database (http://geneticassociationdb.nih.gov/).

Functional gene sets “HALLMARK_ESTROGEN_RESPONSE_EARLY” and “HALLMARK_ESTROGEN_RESPONSE_LATE” were downloaded and extracted from the hallmark gene sets of MSigDB (https://www.gsea-msigdb.org/gsea/msigdb/, v7.2). REACTOME pathways and biological processes information of GO was also downloaded from MSigDB.

We downloaded the transcript annotation from Ensembl and obtained 215,307 annotations. The transcript ID, transcript type, and HUGO Gene Nomenclature Committee (HGCN) symbols were downloaded from Ensembl. Further, the annotations whose transcript type belonged to “lincRNA” or “antisense” were extracted as the lncRNAs. We obtained a total of 1,794 lncRNAs. In addition, we also downloaded the lncRNA annotation file lncipedia_5_0_hg19.gtf (full database) from LNCipedia (12). Considering that some lncRNAs had alternative names, we extracted the Ensembl ID, gene_alias and gene_id for each lncRNA. The lncRNA information from either Ensembl or LNCipedia was used to annotate the lncNRAs.

ER is the most important hormone receptor in breast cancer. We also screened the differentially expressed genes (DEGs) between ER-positive and ER-negative patients from public cohorts and denoted as ER_DEGs markers. The raw read counts for breast cancer was downloaded from The Cancer Genome Atlas (TCGA). The R package DESeq2 (13) was used for differential analysis. 22,946 DEGs were identified with adjusted p <0.05. The Z-score scaled expression profile of Molecular Taxonomy of Breast Cancer International Consortium (METABRIC) was also downloaded. T-test was used as the statistical method to calculate the p value of gene expression difference between ER+ vs ER− samples. The p value was then adjusted by fdr method using the R package fdrtool (14). As a result, we obtained 2,951 DEGs with fdr adjusted p <0.05.

The SSN method was developed by Liu et al. to construct a personalized network for individual patients based on their expression profiles (8). Briefly, a reference network can be constructed using the Pearson correlation coefficient (PCC) between molecules based on the expression data of the reference samples. After a new sample is added to the reference samples, the perturbed network can be similarly constructed. Then, the differential network is constructed by the edges with significantly changed correlation between the reference and perturbed networks.

The changed correlation follows a new type of distribution, which is called the “volcano distribution”. The tail areas of this distribution are similar to those of the normal distribution based on the Kolmogorov–Smirnov test with random sampling. Therefore, the statistical hypothesis Z-test was used to evaluate the significance level of the changed correlation because of the central limit theorem (15).

Liu et al. selected 8–17 normal samples as reference samples to construct an SSN. They also ensured that the SSN is robust and stable for the different reference sample sizes. Inspired by the SSN method, 20 MCF-7 cells captured at 0 h were selected as the reference cells in this study ( Figure 1A ). The reference network was constructed using the PCC. The RNA–RNA correlation was deemed significant with a p-value <0.05. For cells at 3, 6, or 12 h, we added one cell to the reference cells and recalculated the RNA–RNA correlation ( Figure 1B ). We retained a correlation network, named the perturbed network, containing significant RNA–RNA relationships with a p-value <0.05.

Next, we compared the significant RNA–RNA interactions in the perturbed network and reference network to keep only the edges with significantly changed correlations ( Figure 1C ). Figure 1D shows how to test the significance of the changed correlation (ΔCorrelation _N ) between a pair of RNAs. According to the SSN theory proposed by Liu et al. (8), the differential correlation followed a normal distribution, and the significance could be evaluated based on the Z-test.

where N is the number of reference cells.

From the expression perspective, the ceRNAs were positively correlated. Therefore, we considered only the positive and significant differential RNA–RNA interactions as candidate ceRNAs ( Figure 1E ). The ceRNA network from StarBase was further used to filter the ceRNA network to ensure its biological importance ( Figure 1F ).

The hypergeometric test was used to evaluate whether the genes in the CCN were significantly enriched in functional gene sets.

where M is the total number of genes in the background network, N is the number of genes in a functional gene set, m is the number of genes in the CCN, and n is the number of CCN genes shared by the functional gene set.

The R package igraph was used to calculate the topological features of RNAs in the CCN (R 4.0.2). The degree of RNA is the number of direct neighbors in the ceRNA network. RNAs with a high degree can be termed as hub RNAs, which play a pivotal role in maintaining the ceRNA–ceRNA relationships within CCN. The betweenness of RNA i can be calculated with the formula

where RNA s and t are RNAs in the CCN different from RNA i, d_st represents the number of the shortest paths from s to t, and δ_st(i) is the number of the shortest paths from s to t that i lies on. For RNA s and t, the ratio is the proportion of the shortest path that RNA i lies on. The sum of the ratios of all RNA pairs is the betweenness centrality of RNA i. The closeness coefficient is the average closeness of RNA i to other RNAs in the network. It is calculated as

where d_si represents the distance between RNA i and other RNAs.

A recently published paper (7) has provided a web tool, “ceRNA survival”, to perform survival analysis for ceRNA composed of lncRNA–miRNA–mRNA, based on The Cancer Genome Atlas (TCGA) datasets. The web tool was used to perform multivariate Cox regression analysis based on miRNA, mRNA, and lncRNA expression, without co-factors.

GSEA (16) was performed using Preranked utility implemented in the standalone version of the GSEA software (v 4.1.0). The RNA sequencing dataset of DLEU2 knockdown was downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/), with accession number GSE162677. We ranked the genes according to the fold change in expression (FPKM.siDLEU2/FPKM.siNC). The fold change was log2 transformed before GSEA.

The CCN was constructed following the workflow shown in Figure 1 (more details in Materials and Methods). Briefly, the edges with significant differential correlation between the reference network ( Figure 1A ) and perturbed network ( Figure 1B ) were used to construct a cell-specific RCN ( Figure 1C ). miRNA targets are negatively regulated by miRNAs. RNAs competitively bind to the same miRNAs as ceRNAs. Chen et al. generated a ceRNA network for each subtype of breast cancer, based on the principle of positive co-expression and shared miRNAs (17). Therefore, the positive RCN appeared to be a candidate ceRNA, based on ceRNA theory ( Figure 1E ). The ceRNA network from starBase was further used to reduce false-positive ceRNA relations ( Figure 1F ). As an example, we have demonstrated the ceRNA network for the cell MCF-7 12 h (RHM266) in Figure 1G .

ESR1 plays an important role in breast cancer. Therefore, we first examined its interactors in the cell-specific RCN, which is schematically shown in Figure 1E . We selected cells at 3, 6, and 12 h, in which ESR1 had a degree larger than 25. The known estrogen-regulated genes (ERGs), such as KLF4 ( Figures 2A–F ) and TSKU ( Figures 2A, C ), were significantly positively correlated with ESR1. Akaogi et al. reported high expression of KLF4 in ER-α-positive patients. KLF4 was found to bind to the DNA-binding region of ER-α and inhibit the binding of ER-α to estrogen response elements in the promoter regions (18). Known cancer-related genes from CGC, such as MYD88 ( Figures 2A, D ), DDX10 ( Figures 2B, F ), MLLT6 ( Figure 2D ), BCL10 ( Figure 2F ), and KAT6A ( Figure 2F ) also interacted with ESR1. The expression of MYD88 could be modulated in a single nucleotide polymorphism (SNP)- and estrogen-dependent fashion (19). Breast cancer-related genes from GAD were also identified as ESR1 interactors, including TP53BP1 ( Figure 2B ), SRA1 ( Figure 2D ), and BBS4 ( Figure 2E ). Low expression of TP53BP1 is associated with increased local recurrence in breast cancer patients treated with conserving surgery and radiotherapy (20). In addition to these protein-coding genes, ESR1 was also co-expressed with non-coding RNAs such as MIR302B ( Figure 2C ) and MIR4426 ( Figure 2E ). LncRNAs MIR181A1HG ( Figure 2A ), ATP1A1OS ( Figure 2B ), and LINC00094 ( Figure 2D ) were also shown to cross-talk with ESR1. ER_DEGs frequently interacted with ESR1 within individual cells ( Figures 2A–F ). These results indicate that the cell-specific RCN reflect the genes’ regulations of breast cancer.

We also constructed a cell-specific RCN based on the expression profiles of the T47D dataset. ESR1 was also co-expressed with known marker genes, including known ERGs, cancer-related genes from CGC, ER-DEGs, and breast cancer-related genes from GAD ( Figure S1 ).

The CCN was then constructed by integrating the RCN with the starBase ceRNA database, which was developed based on CLIP-Seq data. The average numbers of edges and lncRNAs in the CCN are shown in Table 1 . The average number of edges decreased after E2 stimulation, from 882 edges (3 h) to 662 edges (12 h). Meanwhile, the average number of lncRNAs decreased from 16 (3 h) to 14 (12 h). In contrast, there were more than 1,600 edges on average in the CCN from the T47D dataset. However, on average, five lncRNAs were involved in the CCN.

Cell type classification assumes high importance in single cell heterogeneity. Therefore, we defined cell subtypes by integrating CCN and the estrogen response hallmark. We retrieved estrogen early response and late response hallmarks from MSigDB. For each CCN, we extracted all the RNAs and performed a hypergeometric test to evaluate the extent of RNA enrichment occurring in these hallmark stages.

We selected one CCN at each time point as an example. The CCN was significantly enriched in the early response hallmark for MCF-7 3 h (RHM223, Figure 3A , p = 0.0002), MCF-7 6 h (RHM300, Figure 3B , p = 6.61E-8), and MCF-7 12 h (RHM265, Figure 3C , p = 3.69E-5). Among all 64 cells (3, 6, and 12 h), 41 were enriched in the estrogen early response hallmark ( Figure 3D , p < 0.05). Three out of 64 cells were enriched in the estrogen late response hallmark ( Figure S2 ). Similarly, the cells T47D 3 h (T47D_3 h_2B6), T47D 6 h (T47D_6 h_2H8), and T47D 12 h (T47D_12 h_D8) were enriched in the early response hallmark ( Figures S3A–C ). In the T47D dataset, seven cells were classified as early response cells ( Figure S3D ).

We further classified all 64 cells into two subtypes: early response ones and others. Traditionally, the nodes (RNAs in the network) were screened for biomarker identification. As a complex disease, breast cancer is induced by a set of dysregulated and synergetic genes rather than a single gene. Therefore, network biomarkers are more advantageous for characterizing disease states. Here, we explored the edge markers of both the cell subtypes. We selected the top 20 ceRNA–ceRNA relationships that only appeared in one cell subtype. From the heatmap, we could clearly distinguish the early response subtype from the other subtype ( Figure 3E ). A similar result is shown in Figure S3E for the T47 dataset.

Topological characterization of the CCN is crucial for identifying the pivotal genes that substantially contribute to gene regulation upon E2 stimulus. For all the CCNs, we analyzed the topological features, including degree, betweenness, and closeness coefficients. We focused on lncRNAs in the CCN with a high degree in the early response subtype. The top five lncRNAs are listed in Table 2 . The average degree for all lncRNAs in the early response subtype is in Table S3 .

The role of NEAT1 in breast cancer has been widely investigated. It is also a hub lncRNA in the CCN of early response cells with an average degree as high as 5.14. Cells with a degree ≥20 for NEAT1were selected. SMAD4 and NF1, the known cancer-related genes in the CGC database, are the common ceRNAs of NEAT1 in all six CCNs ( Figures 4A–F ). Another CGC gene, PTEN, appears in five of these six ceRNA networks. WWTR1, a CGC gene, is the ceRNA of NEAT1 in cell MCF-7 3 h (RHM254, Figure 4A ) and MCF-7 6 h (RHM271, Figure 4B ). The GAD genes PRKCA, PRLR, and POLK function as ceRNAs of NEAT1 in cell MCF-7 3 h (RHM254, Figure 4A ), MCF-7 6 h (RHM271, Figure 4B ), and MCF-7 12 h (RHM265, Figure 4C ), respectively. For cell MCF-7 3 h (RHM255, Figure 4D ), we identified the ERGs XRCC1 and RAPGEFL1 as ceRNAs of NEAT1. For cell MCF-7 6 h (RHM271), ERGs RAPGEFL1 and SLC7A2 are ceRNAs of NEAT1. TET2, a CGC gene, is the ceRNA of NEAT1 in cell MCF-7 6 h (RHM250, Figure 4E ).

Among NEAT1 ceRNAs in the six CCNs, we noted that the ceRNA–ceRNA relationship of NEAT1 and ZFX in MCF-7 12 h (RHM266, Figure 4F ) has been recently validated (21). According to a previous study (21), NEAT1 and ZFX competitively bind to miR-138-5p. Next, we performed multivariate Cox regression analysis for the NEAT1–ZFX–miR-138-5p regulation axis using the “ceRNA survival” tool of LnCeCell (7). We divided all breast cancer patients from TCGA into two groups, based on the median expression value of NEAT1–ZFX–miR-138-5p. Patients with high NEAT1–ZFX–miR-138-5p expression had worse overall survival than those with low expression ( Figure 4G ). Moreover, we curated from starBase that miR-493-5p and miR-513a-5p are significantly shared by NEAT1 and ZFX. Because the expression of miR-513a-5p is not available in the TCGA dataset of breast cancer, we tested the prognostic potential of NEAT1–ZFX–miR-493-5p. As shown in Figure 4H , the NEAT1–ZFX–miR-493-5p axis was also an unfavorable prognostic marker of breast cancer.

As NEAT1 is one of the hubs in CCN, we further used the web tool Kaplan–Meier Plotter (https://kmplot.com/analysis) to perform relapse-free and metastasis-free survival analysis for NEAT1. Three probes from the microarrays were mapped to NEAT1. The mean expression of the probes was used for the survival analysis of NEAT1. High and low NEAT1 expression levels were divided according to the median expression level. As shown in Figure 4I , NEAT1 was a prognostic marker for breast cancer, based on relapse-free survival analysis. However, NEAT1 was not a predictor of metastasis-free survival in breast cancer ( Figure 4J ).

We also constructed CCNs based on the T47D dataset. NEAT1 is not a hub lncRNA in the CCNs of T47D cells. The top five lncRNAs in the early response cells are shown in Table S2 for the T47D dataset. Topping the list is PITPNA-AS1. But there is no available siLncRNA dataset for PITPNA-AS1. Therefore, we focus on MALAT1, which has the second largest degree. MALAT1 has been widely investigated for its role in breast cancer. Cells with a degree ≥10 for MALAT1 were selected ( Figures S4A–C ). MALAT1 interacts with known ERGs, cancer-related genes from CGC, ER-DEGs, and breast cancer-related genes from GAD, which is consistent with the results from the MCF-7 dataset. The ceRNA–ceRNA relationships of MALAT1–ZFP36L2, MALAT1–PGRMC2, and MALAT1–PDS5B were shared among the three cells ( Figures S4A–C ). The ceRNA survival analysis revealed that they were unfavorable prognostic markers for breast cancer ( Figures S4D–F ). Survival analysis of the hub lncRNA MALAT1 demonstrated that MALAT1 was a prognostic marker of relapse-free survival ( Figure S4G ) but not metastasis-free survival ( Figure S4H ).

Function prediction and interpretation of lncRNAs are important factors to dissect their biological mechanisms. Therefore, we tried to infer the function of the hub lncRNA DLEU2, which has not been characterized well in breast cancer.

The silencing or overexpression is a commonly used measure of lncRNA function inference. We searched the GEO database for RNA-sequencing datasets generated by siLncRNA or overexpression of DLEU2. As a result, we found that the dataset for siDLEU2 (GSE162677) matched our criteria.

Differential expression analysis is a commonly used method to explore the function of lncRNAs, especially for in silico experiments. The significantly affected biological functions associated with DLEU2 expression could be theoretically identified based on functional enrichment of the genes affected by DLEU2. However, the dataset GSE162677 is generated from a cervical cell line. Thus, it is not suitable for the functional interpretation of DLEU2 in estrogen regulation in breast cancer.

In MCF-7 cells at 12 h (RHM227), DLEU2 had the highest number of ceRNAs. GSEA was used to enrich RNAs in the CCN. The RNAs in RHM227 were significantly enriched in the up-regulated genes after siDLEU2 ( Figure 5A ), which indicates that the genes in the CCN may have similar expression changes after siDLEU2 in breast cancer. Therefore, we further predicted DLEU2 function based on the genes in the CCN and used the hypergeometric test to explore the function of DLEU2 by functional enrichment of RNAs in the CCN of RHM227. The functional terms of the GO and REACTOME pathways were downloaded from MSigDB (v7.2). The top 10 most significant biological processes from GO and pathways from REACTOME are shown in Figures 5B, C , respectively. The most significant GO and REACTOME pathway was GPCR signaling. It should be noted that the biological process “ION_TRANSPORT” was also significantly enriched by RNAs in the CCN of RHM227.

Next, we focused on the ceRNAs of DLEU2 in the CCN of RHM227. SOS1 is involved in GPCR signaling from both GO and REACTOME, while GPR180 participates in GPCR signaling from GO. TSPAN13, a ceRNA of DLEU2, is also a known marker of late estrogen response. It synchronizes with other genes to facilitate ion transport. We retrieved the miRNAs shared by DLEU2, SOS1, GPR180, and TSPAN13. Multivariate Cox regression analysis of ceRNAs demonstrated their prognostic potential in breast cancer ( Figures 5D–F ). Furthermore, relapse-free and metastasis-free survival analysis using Kaplan–Meier Plotter (https://kmplot.com/analysis) revealed that DLEU2 can predict relapse-free survival ( Figure 5G ) but not metastasis-free survival ( Figure 5H ).

To validate the possibility of inferring lncRNA functions via the CCN, we also predicted the function of MALAT1, based on the CCN of T47D cells. The dataset GSE110239 is an mRNA profile generalized by RNA-sequencing for the mammary tumor mouse model PyMT. Mouse genes were mapped to human gene symbols using the R package biomaRt. The genes in the cell T47D_6 h_2C12 was significantly enriched in the up-regulated genes after MALAT1 KO ( Figure S5A ). Further functional analysis of RNAs in the CCN showed that they were enriched in the pathway of “fatty acid metabolism” (p = 0.015). The “fatty acid metabolism” has also been reported to be enriched by DEGs between BPA exposure and control in mouse liver. Several DEGs were key drivers, such as Apoa2, Akr1c12, and Malat1 (22). This provides additional support for the function inference of lncRNAs via the CCN.