LncRNA RP11-10E18.7 cooperates with lncRNA RP11-481C4.2 to affect the overall survival of breast cancer patients: a TCGA-based retrospective study

Introduction

Breast cancer (BC) is one of the leading causes of cancer-related mortality, with approximately 2.3 million new cases (11.7%) per year (1). Despite advances in surgery, chemotherapy, radiotherapy, and neoadjuvant therapy, the mortality rate of BC remains a global challenge (2). In addition, BC also shows a high recurrence rate (3). Due to the pathogenic complexity of BC, many biomarkers have been found, such as estrogen receptor (ER), progesterone receptor (PR), and human epidermal growth factor 2 (HER2) (4,5). However, the prognosis is still a complex problem, indicating the possible existence of new prognosis-influencing molecular mechanisms (6).

Long non-coding RNAs (lncRNAs) are a large class of transcripts from non-protein-coding regions and a group of non-coding RNAs (ncRNAs) of size >200 nucleotides in length (7). LncRNAs are involved in diverse biological processes, such as cell cycle, cell growth, proliferation, migration, invasion, and apoptosis (8). In recent years, lncRNAs have emerged as crucial players in cancer, with specific lncRNAs identified as potential biomarkers for prognosis and therapeutic targets in various cancer types, such as liver, lung, and ovarian cancers (9-12). Further, some lncRNAs may serve as prognosis biomarkers in BC (13-15).

Over the years, mounting studies have revealed that the interactive effects of lncRNA-microRNA (miRNA), lncRNA-messenger RNA (mRNA), and lncRNA-miRNA-mRNA are prognostic signatures in predicting the survival of BC patients (16-18). However, a gap in our understanding pertains to the statistical correlation between interactions among lncRNAs themselves and their impact on the overall survival (OS) of BC. To address this gap, we conducted an analysis using data from the Cancer Genome Atlas (TCGA) data to explore the potential clinical significance of interactive effects between two lncRNAs in the context of BC prognosis. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-23-1941/rc).

Methods

Study populations

The fragments per kilobase per million (FPKM) of lncRNA expression data and corresponding clinical information of BC were derived from the TCGA database (https://portal.gdc.cancer.gov), which includes different types of BC, such as infiltrating ductal carcinomas and infiltrating lobular carcinomas. LncRNA annotation information came from gencode.v22 (https://www.gencodegenes.org/). The study sample was collected between May 2021 and October 2022. The requirement for approval from our Ethics Committee was waived since the study data was acquired from TCGA. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Quality control (QC) procedures for lncRNA expression data

FPKM data and gencode.v22 annotations were matched to obtain the expression of 15,900 lncRNAs. In order to obtain reliable biomarkers in lncRNA expression data, we performed QC, after which we performed statistical association analysis. The exclusion criteria of lncRNAs were when all gene expression values equaled to 0 or the proportion of missing values was greater than 10%. Further, samples with any missing clinical variables and male patients were also removed. Finally, 890 female samples with 4,636 lncRNAs remained in the subsequent statistical association analysis. LncRNA was logarithmically transformed before analysis. The average age was 57.86 years for patients (Table 1). The validation sample was enrolled by the same method, including 157 female samples with 817 lncRNAs.

mRNA expression data

Expression and mRNA sequencing data for all 890 patients were also downloaded from TCGA. Gene expression was measured by RNA sequencing (RNA-seq). Data processing and QC were completed by the TCGA workgroup. RNA-Seq expectation maximization (RSEM) was adopted to normalize raw counts. We downloaded the Level-3 (gene-level) gene quantification data from TCGA and further checked the data quality. Gene expression data were extracted and transformed on log₂ scale before statistical association analysis.

Genome-wide lncRNA-lncRNA interaction analysis

The analysis workflow is displayed in Figure 1. We applied a multivariate Cox proportional hazards model adjusted for age, race, clinical stage, neoadjuvant therapy, ER, and PR to test lncRNA-lncRNA interaction items, by using the R package survival (https://cran.r-project.org/web/packages/survival/index.html). The effect size was measure by hazard ratio (HR) and its 95% confidence interval (CI) according to per 5% expression increment. The significance level accounting for multiple comparison was set based on the Bonferroni method, where the significance level was defined as 0.05 divided by the number of tests. Therefore, we controlled the overall type I error at the 0.05 level. The significance level of lncRNA-lncRNA interaction study was set to be 4.65×10⁻⁹=0.05/(4,636×4,635/2).

Figure 1 Flow chart of study design and statistical analyses. R version 3.6.1 was used for creation of the figure. BC, breast cancer; lncRNA, long non-coding RNA.

Gene-gene interaction analysis

For a better exploration of the function of the significant lncRNA-lncRNA items, the related mRNAs were identified by co-expression methods based on the Pearson correlation, respectively. The related mRNAs were screened according to |COR| >0.35, Bonferroni adjusted P<0.0001. For the significant co-expression genes, we also applied a multivariate Cox proportional hazards model adjusted for the aforementioned covariates to test gene-gene interactions.

Statistical analysis

Continuous variables were described as mean ± standard deviation and compared by Student’s t-test using R package t-test (19), and these categorized variables were summarized as frequency (n) and proportion (%) and compared by chi-square test. All statistical analyses were performed using R version 3.6.1 (The R Foundation for Statistical Computing, Vienna, Austria), unless otherwise specified. P values less than 0.5 were considered as statistically significance.

Results

For the lncRNA-lncRNA interaction analysis, only RP11-10E18.7×RP11-481C4.2 (HR_interaction =1.04, 95% CI: 1.03–1.06, P=3.35×10⁻⁹) was significantly associated with BC OS. As presented in Figure 2A, with the increased expression level of RP11-481C4.2, there was the elevated risk for RP11-10E18.7 on BC OS. Therefore, RP11-481C4.2 was a modifier of the association between RP11-10E18.7 and BC survival. To illustrate the modification effect, patients were categorized into low and high groups based on the intersection (0.69) of RP11-481C4.2 expression in Figure 2A. The effect of RP11-10E18.7 varied across patients with different RP11-481C4.2 expressions. For patients with a low level of RP11-481C4.2, high expression of RP11-10E18.7 had significantly better OS (HR_{highvs. low} =0.54, 95% CI: 0.34–0.84, P=7.13×10⁻³) (Figure 2B,2C). Conversely, high expression of RP11-10E18.7 had poor effect on patients with high level of RP11-481C4.2 (HR_{highvs. low} =2.36, 95% CI: 1.06–5.25, P=3.37×10⁻³) (Figure 2B,2D). The results showed that the effect direction of RP11-10E18.7 on BC OS was opposite at the low and high levels of RP11-481C4.2. To further evaluate our findings regarding RP11-10E18.7 and RP11-481C4.2, we repeated the interaction analysis in the validation dataset of 817 gene expression samples, and the results were similar (HR_interaction =1.02, 95% CI: 1.01–1.05, P=3.01×10⁻⁹).

Figure 2 RP11-10E18.7 and RP11-481C4.2 interaction on overall survival of BC patients. (A) HR of RP11-10E18.7 estimated based on the expression level of RP11-481C4.2. The shallow area represents 95% CI, with grey and blue areas indicating low and high expression, respectively. Histogram on the top shows the distribution of RP11-481C4.2 expression. (B) Forest plots of the effects of RP11-10E18.7 among BC patients with low or high expression of RP11-481C4.2. P_{heterogeneity} was used to evaluate the heterogeneity of HRs across groups. (C) Kaplan-Meier survival curves of low and (D) high expression of RP11-10E18.7 among BC patients with varying RP11-481C4.2 expression levels. R version 3.6.1 was used for creation of the figure. HR, hazard ratio; CI, confidence interval; BC, breast cancer.

Based on the BC-related mRNA expression data from TCGA, co-expression analysis was performed on RP11-10E18.7 and RP11-481C4.2 respectively. The results showed that 14 and 139 genes were closely related to the 2 lncRNAs, respectively, that were then used for gene-gene interaction analysis (Figure 3). Finally, 4 pairs of genes were identified with Bonferroni adjusted P≤0.05 (Table 2). However, as shown in Figure 4, there were high correlations among the expression of these genes, and the subsequent analysis focused on the FOXA1×U2SURP that had the lowest P value (HR_interaction =1.49, 95% CI: 1.28–1.73, P=2.16×10⁻⁷). The associations between RP11-10E18.7 and FOXA1 (r=−0.37, P=3.70×10⁻³⁰), RP11-481C4.2 and U2SURP (r=0.38, P=1.37×10⁻³²) were all significant (Figure S1).

Figure 3 Co-expression analysis of BC patients from TCGA cohort. (A) RP11-10E18.7 related mRNAs. (B) RP11-481C4.2 related mRNAs. Blue points represent P values of correlation between gene expression and lncRNA, ordered by genomic position. Grey lines represent significant connections with |COR| >0.35 and Bonferroni adjusted P≤0.05. R version 3.6.1 was used for creation of the figure. BC, breast cancer; TCGA, The Cancer Genome Atlas; mRNA, messenger RNA; lncRNA, long non-coding RNA.

Figure 4 Correlations among the expression of genes from significant gene-gene interactions. R version 3.6.1 was used for creation of the figure.

The similar interactive pattern was exhibited between FOXA1 and U2SURP (Figure 5). FOXA1 had an elevated risk on BC OS with the increased expression level of U2SURP (Figure 5A). Again, after classifying the patients by intersection (3.47), we observed that the effect of FOXA1 opposite in patients with different U2SURP expression level (HR_{highvs. low} =0.58, 95% CI: 0.34–0.99, P=0.046 in low expression of U2SURP; HR_{highvs. low} =1.56, 95% CI: 1.18–2.87, P=0.029 in high expression of U2SURP) (Figure 5B-5D).

Figure 5 FOXA1 and U2SURP interaction on survival of BC patients. (A) HR of FOXA1 estimated based on the expression level of U2SURP. The shallow area represents 95% CI, with gray and blue areas indicating low and high expression, respectively. Histogram on the top shows the distribution of U2SURP expression. (B) Forest plots of the effects of FOXA1 among BC patients with low or high expression of U2SURP. P_{heterogeneity} was used to evaluate heterogeneity of HRs across groups. (C) Kaplan-Meier survival curves of low and (D) high expression of FOXA1 among BC patients with varying U2SURP expression levels. R version 3.6.1 was used for creation of the figure. BC, breast cancer; HR, hazard ratio; CI, confidence interval.

Discussion

In this study, we used transcriptional data from TCGA to identify RP11-10E18.7×RP11-481C4.2 interaction that had an impact on BC OS. Then, through co-expression analysis, 4 pairs of gene-gene interactions, including FOXA1×U2SURP, were also found to be associated with BC OS.

Approximately 93% of DNA can be transcribed as RNA in the human genome, 98% of which is known as ncRNAs (20). Among them, lncRNAs are RNAs longer than 200 nucleotides in length that have been demonstrated to play important roles in epigenetic control and the regulation of transcription and translation (21). As a new type of gene regulator, lncRNA is related to the occurrence, development, and prognosis of human diseases, especially cancer (22,23). In recent years, lncRNAs have been demonstrated to be engaged in BC development, progression, invasion, and metastasis (24-27).

Several studies have demonstrated that the abnormal expression of some lncRNAs may be associated with the poor prognosis of BC by testing their main effects (28,29). To the best of our knowledge, this may be the first attempt to explore the relationship between lncRNA-lncRNA interactions and BC OS at the population level. Interactions have been shown to provide important clues to the biologic mechanisms of complex diseases (30). Besides this, interactions could increase the power to detect associations and then be leveraged for the identification of new biomarkers (31). Previous studies have identified some gene-gene interactions related to BC survival (32,33). Our results found that biomarkers with RP11-10E18.7×RP11-481C4.2 interaction and FOXA1×U2SURP interaction significantly affected the prognosis of BC.

LncRNA RP11-10E18.7 was found to be associated with miRNA hsa-miR-181a-5p that regulated gene SRPK2 in blood (34). Forkhead box A1 (FOXA1) is a forkhead box transcription factor expressed in mammary luminal epithelial cells (LECs) (35). FOXA1 mutations are a hallmark of ER⁺ BC and have been widely regarded as a determinant of breast tumor response to endocrine therapy as well as a marker for favorable patient prognosis (36). Previous study has indicated that elevated FOXA1 expression level is associated with better outcome in BC (37). However, our results showed that high level of FOXA1 expression had protective effect on BC survival only when U2SURP expression was low, meaning that U2SURP modified the relationship between FOXA1 and BC OS. Meanwhile, U2SURP was also found to be significantly associated with BC survival (38).

Our study has several advantages. First, to our knowledge, this is the first study of the relationship between lncRNA-lncRNA interaction and the OS of BC. The RP11-10E18.7×RP11-481C4.2 interaction provides potential evidence that complex disease is driven by intricate association patterns. Second, our study used co-expression analysis to find gene-gene interactions related to lncRNAs and impacting BC OS, which was a comprehensive evaluation of lncRNAs. Third, for the lncRNA-lncRNA interaction analysis, we applied the most conservative Bonferroni correction to control for false positives.

The study has some limitations. First, our study lacked independent validation. Additional available public databases and further researches are warranted. Second, biological evidence requires further functional experiments of lncRNAs, not just our statistical evidence. This association should be interpreted with caution.

Conclusions

Our study identified RP11-10E18.7×RP11-481C4.2 and FOXA1×U2SURP interactions as potential biomarkers for BC prognosis. Our findings elucidated potential molecular mechanisms of BC progression under complex association patterns and provided potential dynamic and reversible therapeutic targets for BC patients.

Acknowledgments

Funding: This study was funded by Outstanding Young Medical Talents Program of Pudong Health Bureau of Shanghai (No. PWRq2020-40).

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-23-1941/rc

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-23-1941/prf

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-23-1941/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Sung H, Ferlay J, Siegel RL, et al. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin 2021;71:209-49. [Crossref] [PubMed]
2. Valachis A, Nearchou AD, Lind P. Surgical management of breast cancer in BRCA-mutation carriers: a systematic review and meta-analysis. Breast Cancer Res Treat 2014;144:443-55. [Crossref] [PubMed]
3. Yang B, Chou J, Tao Y, et al. An assessment of prognostic immunity markers in breast cancer. NPJ Breast Cancer 2018;4:35. [Crossref] [PubMed]
4. Senkus E, Kyriakides S, Ohno S, et al. Primary breast cancer: ESMO Clinical Practice Guidelines for diagnosis, treatment and follow-up. Ann Oncol 2015;26:v8-30. [Crossref] [PubMed]
5. Zhu A, Yun Z, You M, et al. Surgical reduction in chest wall disease to prolong survival in breast cancer patients: a retrospective study. Gland Surg 2022;11:1015-25. [Crossref] [PubMed]
6. Arnedos M, Vicier C, Loi S, et al. Precision medicine for metastatic breast cancer--limitations and solutions. Nat Rev Clin Oncol 2015;12:693-704. [Crossref] [PubMed]
7. Matsui M, Corey DR. Non-coding RNAs as drug targets. Nat Rev Drug Discov 2017;16:167-79. [Crossref] [PubMed]
8. Bhan A, Soleimani M, Mandal SS. Long Noncoding RNA and Cancer: A New Paradigm. Cancer Res 2017;77:3965-81. [Crossref] [PubMed]
9. Lin C, Yang L. Long Noncoding RNA in Cancer: Wiring Signaling Circuitry. Trends Cell Biol 2018;28:287-301. [Crossref] [PubMed]
10. Parasramka MA, Maji S, Matsuda A, et al. Long non-coding RNAs as novel targets for therapy in hepatocellular carcinoma. Pharmacol Ther 2016;161:67-78. [Crossref] [PubMed]
11. Pan J, Fang S, Tian H, et al. lncRNA JPX/miR-33a-5p/Twist1 axis regulates tumorigenesis and metastasis of lung cancer by activating Wnt/β-catenin signaling. Mol Cancer 2020;19:9. [Crossref] [PubMed]
12. Xu D, Guo Y, Lei S, et al. Identification and Characterization of TF-lncRNA Regulatory Networks Involved in the Tumorigenesis and Development of Adamantinomatous Craniopharyngioma. Front Oncol 2022;11:739714. [Crossref] [PubMed]
13. Crudele F, Bianchi N, Reali E, et al. The network of non-coding RNAs and their molecular targets in breast cancer. Mol Cancer 2020;19:61. [Crossref] [PubMed]
14. Tomar D, Yadav AS, Kumar D, et al. Non-coding RNAs as potential therapeutic targets in breast cancer. Biochim Biophys Acta Gene Regul Mech 2020;1863:194378. [Crossref] [PubMed]
15. Chen J, Ling C. Construction of a predictive model for breast cancer metastasis based on lncRNAs. Transl Cancer Res 2023;12:387-397. [Crossref] [PubMed]
16. Müller V, Oliveira-Ferrer L, Steinbach B, et al. Interplay of lncRNA H19/miR-675 and lncRNA NEAT1/miR-204 in breast cancer. Mol Oncol 2019;13:1137-49. [Crossref] [PubMed]
17. Jiang YZ, Liu YR, Xu XE, et al. Transcriptome Analysis of Triple-Negative Breast Cancer Reveals an Integrated mRNA-lncRNA Signature with Predictive and Prognostic Value. Cancer Res 2016;76:2105-14. [Crossref] [PubMed]
18. Fan CN, Ma L, Liu N. Systematic analysis of lncRNA-miRNA-mRNA competing endogenous RNA network identifies four-lncRNA signature as a prognostic biomarker for breast cancer. J Transl Med 2018;16:264. [Crossref] [PubMed]
19. Ramsden JD. Angiogenesis in the thyroid gland. J Endocrinol 2000;166:475-80. [Crossref] [PubMed]
20. Mattick JS. Non-coding RNAs: the architects of eukaryotic complexity. EMBO Rep 2001;2:986-91. [Crossref] [PubMed]
21. Carrieri C, Cimatti L, Biagioli M, et al. Long non-coding antisense RNA controls Uchl1 translation through an embedded SINEB2 repeat. Nature 2012;491:454-7. [Crossref] [PubMed]
22. Liang C, Qi Z, Ge H, et al. Long non-coding RNA PCAT-1 in human cancers: A meta-analysis. Clin Chim Acta 2018;480:47-55. [Crossref] [PubMed]
23. Chen Y, Li Z, Chen X, et al. Long non-coding RNAs: From disease code to drug role. Acta Pharm Sin B 2021;11:340-54. [Crossref] [PubMed]
24. Huang QY, Liu GF, Qian XL, et al. Long Non-Coding RNA: Dual Effects on Breast Cancer Metastasis and Clinical Applications. Cancers (Basel) 2019;11:1802. [Crossref] [PubMed]
25. Zhang Y, Huang X, Liu J, et al. New insight into long non-coding RNAs associated with bone metastasis of breast cancer based on an integrated analysis. Cancer Cell Int 2021;21:372. [Crossref] [PubMed]
26. Sha R, Wu Z, Xu Y, et al. Predictive value of lncRNA LOC100505851 in breast cancer in the neoadjuvant setting. Gland Surg 2021;10:1899-909. [Crossref] [PubMed]
27. Zhang Y, He W, Zhang S. Seeking for Correlative Genes annd Signaling Pathways with Bone Metastasis from Breast Cancer by Integrated Analysis. Front Oncol 2019;9:138. [Crossref] [PubMed]
28. Yang J, Zhang X, Ye Y, et al. Postmastectomy radiation therapy can improve survival for breast cancer patients with 1–3 positive axillary lymph nodes: a retrospective cohort study using the SEER database. Transl Cancer Res 2021;10:1984-2001. [Crossref] [PubMed]
29. Li W, Jia G, Qu Y, et al. Long Non-Coding RNA (LncRNA) HOXA11-AS Promotes Breast Cancer Invasion and Metastasis by Regulating Epithelial-Mesenchymal Transition. Med Sci Monit 2017;23:3393-403. [Crossref] [PubMed]
30. Trerotola M, Relli V, Simeone P, et al. Epigenetic inheritance and the missing heritability. Hum Genomics 2015;9:17. [Crossref] [PubMed]
31. Cordell HJ. Detecting gene-gene interactions that underlie human diseases. Nat Rev Genet 2009;10:392-404. [Crossref] [PubMed]
32. Chen X, Theobard R, Zhang J, et al. Genetic interactions between INPP4B and RAD50 is prognostic of breast cancer survival. Biosci Rep 2020;40:BSR20192546. [Crossref] [PubMed]
33. Dai X, Fagerholm R, Khan S, et al. INPP4B and RAD50 have an interactive effect on survival after breast cancer. Breast Cancer Res Treat 2015;149:363-71. [Crossref] [PubMed]
34. Su L, Wang C, Zheng C, et al. A meta-analysis of public microarray data identifies biological regulatory networks in Parkinson's disease. BMC Med Genomics 2018;11:40. [Crossref] [PubMed]
35. Carroll JS, Liu XS, Brodsky AS, et al. Chromosome-wide mapping of estrogen receptor binding reveals long-range regulation requiring the forkhead protein FoxA1. Cell 2005;122:33-43. [Crossref] [PubMed]
36. Hurtado A, Holmes KA, Ross-Innes CS, et al. FOXA1 is a key determinant of estrogen receptor function and endocrine response. Nat Genet 2011;43:27-33. [Crossref] [PubMed]
37. Shou J, Lai Y, Xu J, et al. Prognostic value of FOXA1 in breast cancer: A systematic review and meta-analysis. Breast 2016;27:35-43. [Crossref] [PubMed]
38. An J, Luo Z, An W, et al. Identification of spliceosome components pivotal to breast cancer survival. RNA Biol 2021;18:833-42. [Crossref] [PubMed]