Full Text Article

Introduction

Breast cancer (BC) is a common and aggressive malignant tumor that originates from breast tissue, especially in women. Every year, nearly 1.3 million people are diagnosed with BC, and 450,000 people died from BC in the world [1]. The incidence of BC is getting younger, while the morbidity and mortality rates remain high. Breast cancer is highly heterogeneous, with numerous subtypes, and clinical tests typically assess the expression of four proteins, including Ki67 (a protein associated with breast cancer proliferation), Estrogen Receptor (ER), Human epidermal growth factor receptor 2 (HER2), and Progesterone Receptor (PR). Depending on the expression, if cancer cells are negative for the receptor, this indicates that the receptor is not present in the tumor tissue, e.g. cells with PR receptor are PR-negative or PR-; conversely, if cancer cells are positive for the receptor, this indicates that this receptor is present in the tumor tissue, e.g. cells with PR receptor are PR-positive or PR+. According to these criteria, breast cancer was classified into four subtypes in clinical studies[2-4], as shown in Table 1. Similarly, breast cancer is highly metastatic, with 20-30% of breast cancer patients developing metastases after diagnosis and treatment, with sites including bone, lung, liver and brain[5,6]. Whereas metastasis leads to a decrease in overall patient survival, the 5-year overall survival rate is 22.8% for bone metastases and 16.8% for lung metastases, and brain metastases are the most risky factor among patients with metastases[7-9]. Currently, BC diagnosis and treatment methods include mammography, computed tomography (CT), magnetic resonance imaging (MRI) and ultrasound imaging, and other imaging methods[10-13].Trials have shown that mammography can reduce breast cancer mortality by 20%, and ultrasound can also detect cancer patients, but the positive predictive value of ultrasound is only 3-8%[14, 15]. MRI, although proven to have a high cancer detection rate in breast cancer imaging modalities, is neither applicable nor cost-effective in screening large numbers of women at high risk for breast cancer [16]. Therefore, defects in sensitivity and low specificity at the time of diagnosis can cause patients to miss the best time for treatment. In addition, the high metastasis and invasion of BC also put forward more stringent requirements for clinical treatment. Considering the complexity and the unclear pathogenesis of the clinical and pathological symptoms of BC, it is still necessary to conduct in-depth research so as to fully grasp the regulatory mechanism of BC.

As a type of RNA that does not encode protein, NcRNA includes long non-coding RNAs (lncRNAs), microRNAs (miRNAs), and circular RNAs (circRNAs) and so on [17].. More and more evidences have shown that non-coding RNA (ncRNA) has played an important role in many human tumors including breast cancer by involving cell development, proliferation, differentiation and apoptosis [18]. Firstly, proposed the competing endogenous RNA hypothesis, named ceRNA, suggesting that the presence of MRE can relieve the inhibition of miRNA to the target mRNA by combining with miRNA. An increasing number of studies suggested that lncRNA can act as ceRNA or sponge to indirectly regulate mRNA expression by blocking miRNA [19]. In addition, ceRNA regulatory network analysis can be used for mining potential biomarkers, providing valuable information for exploring the mechanism of tumor occurrence and development. For example, LINC0092 and C2orf71 has been used as new prognostic indicators for estrogen receptor positive (ER+) and ER negative (ER-) breast cancer [20], RP11-434D9.1, LINC00052, BC016831 and IGKV are related to the occurrence of TNBC(triple negative breast cancer), which can be used as a potential biomarker for clinical diagnosis or treatment targets of TNBC [21]. Besides, lncRNA can be used as a new biomarker to determine the subtype of breast cancer [22]. These studies demonstrated that ceRNA network disorder is related to the occurrence of breast cancer, and the construction of ceRNA regulatory network is an effective way to study the occurrence of cancer. However, the current research on ceRNA regulatory networks related to breast cancer is scarce.

In the current study, we aimed to conduct research on BC targeted genes. Through the GEO database, we identified differentially expressed RNA by bioinformatics methods and constructed lncRNA-miRNA-mRNA regulatory network. And through molecular docking technology to find the relevant target information that is helpful for early prevention and diagnosis of breast cancer provides theoretical support for new drug development, and also opens the way for precision treatment.

Materials and methods

The Cancer Genome Atlas Program

GEO(https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi) is a public open gene expression profile database established and maintained by the National Center of Biotechnology Information (NCBI) in 2000, which can store high-throughput gene expression data sets [23]. Potential reasonable RNA-Seq databases were selected based on the following criterion:

(1) The research organism is homo sapiens;
(2) Specimens with RNA;
(3) Human BC and control Expression data from breast;
(4) Expression profiling by array;
(5) Supported by GEO2R analysis.

Use the keywords: (("Breast Cancer") AND ("RNA" OR "Expression profiling by array") AND ("Homo sapiens")) to find potential datasets. Then, the databases of GSE5764, GSE7904, GSE8977 and GSE10810 were selected, and the samples were divided into BC group and control group for the next study, respectively, according to the predefined datasets inclusion criteria. Gene expression data from breast cancer are downloaded from TCGA (The Cancer Genome Atlas Program, https://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga) [24] with clinical data and screened and organized for analysis.

Identification of Differentially Expressed RNA

In order to have a better view of RNA with significant expression changes in different samples, we identified differentially expressed RNAs(DERNAs) through GEO2R, a useful online analysis tool that comes with the GEO database [25]. Namely, enter a GSE number into the GEO accession box. After all samples (GSMs) were assigned to groups we wanted to compare (BC vs. Normal), we performed the test with default parameters and clicked “Save all results” to obtain all the DERNAs. The threshold for DERNAs is P.Value<0.05 and |log FC|≥2. The Venn diagram in FunRich_v1.3.1 [26] was used to show the common DERNAs identified by GEO2R in the four datasets.

Bioinformatics Analysis of Differentially Expressed RNA

In order to obtain a comprehensive understanding of mechanism and progression in BC patients, we performed Gene Ontology including biological process, molecular function, cellular component and Kyoto Encyclopedia of Genes and Genomes (KEEG) analysis on DERNAs using FunRich_v1.3.1 and DAVID to predict their possible functions and signaling pathways [27].

We adopted online tool STRING[28] (https://string-db.org/) to develop protein-protein interaction (PPI) network of differentially expressed mRNAs (DEmRNAs). Briefly,PPI network was showed by inputing all gene symbol into the search box and selecting correct organism (homo sapiens). Cytoscape_v3.6.1 is used to visualize the network [29].

Construction of ceRNA Regulatory Network

The combination of MiRNA with lncRNA was predicted by databases starbase [30], miRNA targeting genes intersection common mRNA were obtained through online databases mirDIP [31]. Organized and sort out the corresponding data of lncRNA-miRNA-mRNA. The established ceRNA network were visualized with Cytoscape_v3.6.1.

Statistical Methods

The genes in the ceRNA network were individually subjected to monofactor analysis using the R programming language, and the mRNAs were screened again using least absolute shrinkage and selection operator (Lasso) analysis, followed by cox regression of the screened variables to obtain a risk model. The Kaplan-Meier method was used for survival analysis, and the predictive performance of the risk model was assessed by the time-dependent receiver operating characteristic (ROC) of the "SurvivalROC" R package. In addition, the Mann-Whitney test was used to assess differences between sample groups, and the Kruskal-Wallis test and Dunn's multiple test were used for multiple comparisons to calculate the median, with P<0.05 being considered a statistically significant difference.

Molecular Docking

After downloading the protein structure corresponding to the 12 common DEmRNAs from RCSB database [32] (https://www.rcsb.org/pdb), we consulted the relevant literatures to get related treatment drugs of BC (Tamoxifen, Lapatinib, Eribulin and Afimoxifene) as ligands for molecular docking (table 2). Both ligands and protein receptors were prepared using AutoDock_v4.0 [33]. Before the molecular docking, we pre-treated protein receptors (PDB ID: 3FBY, 1XED, 1O7Y) and drug molecules (hydrogenation, charging, etc.) and calculated the grid energy. the GirdBox grid region covered the studied target proteins was used as the center (the active site is unknown). In the docking process, we employed the genetic algorithm for regional search, the initial population was set to 150, the maximum number of energy assessments was 2.500000, and other parameters were selected by default. After the docking, binding energy (BE) was used as the criterion to select the optimal conformation. The greater energy indicated the stronger interaction between molecules. We selected this conformation as the optimal result of docking. Finally, the results was imported into PyMOL_v2.3.3 to visualize the conformation of the complex and analyze the interaction between protein and drug [34].

Results

Identification of Differentially Expressed RNA

Through searching and inclusion criterion, we found four databases of GSE5764, GSE7904, GSE8977 and GSE10810G from GEO (Table 3). All RNA-Seq data are from GPL570 platform. After processed by GEO2R, we obtained DERNAs with the criterion (P. Value<0.05 and |log FC|≥2) (Figure 2). Meanwhile, 12 DEmRNAs (SFRP1, PLPP4, PIGR, MT1M, MMP11, INHBA, EMP1, CXCL10, COMP, COL11A1, COL10A1 and C2orf40) identified by GEO2R were used for subsequent drug target analysis (Figure. 2).

Bioinformatics Analysis of Differentially Expressed Gene

The functional enrichment analysis of DEmRNAs was performed using DAVID and FunRich. DAVID enrichment factor with P.value<0.05 were thought to be significant. The top ten terms were visualized by R language. The results showed that DEmRNAs are related to cytokine activity and protein binding and other biological functions, and can affect cell growth and communication (Figure 3 and Figure4).

We also performed enrichment analysis on 12 common DEmRNAs with DAVID (Figure 5) and found that these DEmRNAs are related to the two pathways of protein digestion and absorption and ECM-receptor interaction, and are related to endoderm cell differentiation, negative regulation of cell growth, and biological functions such as the development of multicellular organisms, and can also predict drug response.

The DEmRNAs identified from the GSE122063 database was imported into the SRTING to construct PPI network. The isolated protein does not participate in the construction of the network, and hides it to obtain the protein interaction map of the DERNAs. From this PPI network,it was found that the numbers of node and edge are 683 and 6236, respectively. average node degree is 18.3, average local clustering coefficient is 0.44(Figure 6b). After importing 12 common DEmRNAs into the SRTING database, only 4 mRNAs were found to involve in the formation of protein network, namely COL10A1, COMP, MMP11 and COL11A1. From this PPI network, it was found that the number of nodes is 12, the number of edges is 4, the average node degree is 0.667, and the average local clustering coefficient is 0.278(Figure 6a).

Construction of ceRNA regulatory network

Through the online databases mirDIP and ENCORI database, we got 15 lncRNA, 2 miRNAs, and 72 mRNA. Based on the relationship between miRNA and mRNA, lncRNA, we obtained ceRNA regulatory network and visualize it through Cytoscape (Figure 7).

Building a risk model based on ceRNA regulatory network

Building a risk model based on mRNA

There were 1076 samples matching in TCGA database, 72 mRNA expression matrices in ceRNA network were extracted, and 18 genes were screened using lasso analysis: LRP1, CLDN11, AXIN2, NKAPD1, PPP1R3A, ZNF800, BMPR1B, EDN3, HMGB3, CALN1, GEMIN6, EGR3, ZIC1, SPIB, EDNRB, RCAN1, EGR1, BCAP29, for which we also set a 10-fold cross-validation to construct the optimal risk model as Risk Score=0.088*LRP1-0.042*CLDN11+0.053*AXIN2-0.084* NKAPD1-0.163*PPP1R3A-0.036*ZNF800+0.002*BMPR1B-0.021*EDN3+0.008*HMGB3-0.096*CALN1-0.096*GEMIN6-0.008*EGR3+0.028*ZIC1-0.049*SPIB- 0.0189*EDNRB-0.022*RCAN1-0.039*EGR1-0.025*BCAP29. In addition, a multifactorial survival analysis based on model genes identified three genes, LRP1, ZIC1, and AXIN2, as risk genes with a risk ratio greater than 1, which are closely associated with breast cancer prognosis (Figure 8A). We used the R package max stat to calculate the optimal cutoff value of Risk Score, setting the minimum grouping sample size greater than 25% and the maximum sample size grouping less than 75%, and finally obtained the optimal cutoff value of -0.539204150358237, based on which the patients were divided into two groups, high and low, and the overall survival rate of patients in the low risk group was better than that in the high risk group (Figure 8B), and in addition the ROC curves at 1, 3 and 5 years showed the more accurate predictive ability of the model with the area under the curve of 0.60, 0.63 and 0.64, respectively (Figure 8C). We then analyzed the relationship between different risk scores and the follow-up time, events and changes in the expression of each gene in patients, and it was observed that the survival rate of patients decreased significantly with increasing risk scores, and as expected, all three genes, LRP1, ZIC1 and AXIN2, were risk factors, and the expression showed an up-regulation trend with increasing risk scores (Figure 8D). We also calculated the immune infiltrating cell score for each sample based on the expression profile of breast cancer tumor samples using the MP Counter method selected by the R package IOBR and found six types of immune cells in patients in the high-risk and low-risk groups: T cells, CD8+ T cells, cytotoxic lymphocytes, NK cells, B lymphocytes, myeloid dendritic cells, and two stromal cells, endothelial cells and fibroblasts, were significantly different in abundance (Figure 8E).

Independence of the Constructed Risk Model

Clinical characteristics of breast cancer were downloaded from TCGA and selected age, sex, histological type, TNM stage, stage, new event, cancer status, breast carcinoma estrogen receptor status, breast carcinoma progesterone receptor status, lymph node examined count, margin status, number of lymph nodes positive by he, and race along with the risk model for univariate survival analysis. In univariate analysis, age, cancer status, lymph node examined count, M, new event, race, and risk score were statistically significant, and then using them in multivariate survival analysis, risk Score, cancer status, new event were still statistically significant (Table 4), and the HR of risk score reached 6.4, which was an important influencing factor for patient prognosis.

When the risk score was analyzed together with cancer status and new event, it was found that there was no significant difference in risk score among patients with different status of cancer status and new event and the risk score could also classify patients with different status into high-risk and low-risk groups, and there was a significant prognostic difference between the two groups (Figure 9). These results suggest that the risk model we constructed has good independence in predicting breast cancer prognosis.

Construction and Calibration of an Integrated Monogram

Two TBX2-AS1 and MUC20-OT1 lncRNAs that were significant for patient prognosis, with HR values of 1.21 and 1.19, respectively. Finally, to more accurately predict the prognosis of breast cancer patients, a risk model was constructed with clinical characteristics including cancer status, new event, TBX2-AS1, and MUC20-OT1. A nomogram was constructed, as shown in Figure 10. It was assigned specific scores based on the contribution of risk scores and pathological features to the prognosis of breast cancer. The overall C-index of the model was: 0.625,95% CI (0.603-0.646), P value=1.415e-29, and the observed overall survival was better matched with the actual survival at 3 and 5 years (Figure 10B), in addition the model was able to better classify patients into two groups of high and low risk, and the area under the curve reached 0.60, 0.67, and 0.69, suggesting that the risk model constructed based on ceRNA network can reliably and accurately predict the prognosis of breast cancer patients.

Molecular Docking

According to the common mRNA obtained in the GEO database, three corresponding target proteins (COMP, PIGR and CXCL10) downloaded from the RCSB database were used as receptors for molecular docking (Table 5). In molecular docking, the ligands are Tamoxifen, Lapatinib, Eribulin and Afimoxifene. The docking results showed the poor docking effect of CXCL10 and Lapatinib and the strong docking effects of other target proteins with drug ligands, this result indciated that the three proteins can be used as potential sites for breast cancer therapy (Table 6). We imported the best conformation into PyMOL, and analyze the three-dimensional effect of the docking results (Figure 11).

Discussion

BC is a common malignant tumor with extremely high mortality. Due to the lack of specific diagnostic and prognostic biomarkers, the survival rate of BC patients is very low. In order to study the clinical results and explore the regulatory mechanism of the initiation and progression of BC, the RNA expression data in the GEO database was used to construct a ceRNA regulatory network by conducting gene enrichment and protein interaction analysis using bioinformatics methods. Meanwhile, through molecular docking, the interaction between target protein and related therapeutic drugs was studied to find potential drug targets for the treatment of BC.

In this study, we conducted an in-depth analysis of 3 target proteins (COMP, PIGR and CXCL10). Cartilage oligomeric matrix protein (COMP) is a soluble pentameric protein that is expressed in cartilage. It can lead to the severity of the disease by increasing the invasiveness and viability of tumor cells through metabolic conversion [35]. Studies have found that the imbalance of COMP expression is significantly related to the aggressiveness of breast cancer, which may affect the survival rate and recurrence rate of patients [36,37]. The polymeric immunoglobulin receptor (pIgR) is a key component of the mucosal immune system that mediates the transcytosis of immunoglobulin epithelial cells. In addition, up-regulation of pIgR expression was reported to relate to several human cancers [38,39], but there is currently little evidence in clinical trials. CXCL10 is a CXC chemokine [40,41], which mainly exerts its biological effects through CXCR3 [42], and is related to many human diseases. Bioinformatics and luciferase reporter gene analysis indicate that miR-34a can be used as a target of CXCL10, and over-expressed miR-34a can significantly inhibit breast cancer cell growth, invasion, migration and induce apoptosis [43,44]. Furthermore, LncRNA UFC1 can regulate the biological activity of breast cancer through the miR-34a/CXCL10 axis, which provide new diagnostic biomarkers and potential therapeutic targets for breast cancer [44].

There are four kinds of drugs as molecular receptors: Tamoxifen, Lapatinib, Eribulin and Afimoxifene. Among them, Tamoxifen, Lapatinib and Eribulin have been approved by the FDA as clinically effective drugs for the treatment of BC, and Afimoxifene, the active metabolite of tamoxifen, is under clinical phase II study [45]. Tamoxifen, an estrogen receptor antagonist, is the most commonly used drug for the treatment of estrogen receptor-positive breast cancer, or as an adjuvant therapy in the early stage of the disease [46-48]. Lapatinib is a dual tyrosine kinase inhibitor (TKI), which can block HER1 and HER2 tyrosine kinase activity by binding to the ATP-binding site of the receptor's intracellular domain [49]. Erebrin is highly effective for metastatic breast cancer (MBC), and combined treatment studies are currently undergoing combined therapy research [50].

Using molecular docking, we analyzed the interaction between the target protein and the drug and found no interaction between CXCL10 and Lapatinib. The docking effects of the other target proteins and drug were good, but the docking effects were inconsistent. According to the sequence of binding energy for COMP, Eribulin> Afimoxifene> Tamoxifen> Lapatinib, for PIGR, Eribulin> Lapatinib> Tamoxifen> Afimoxifene, for CXCL10, Afimoxifene>Eribulin> Tamoxifen. According to this ranking, it is found that the Leu residue of COMP and the GLN residue of PIGR can interact with Eribulin, and CXCL10 residues GLN interacts with Afimoxifene. Based on this finding, we suspect that these related amino acid residues of COMP, PIGR and CXCL10 are related with the efficiency of certain drugs for the treatment of breast cancer.

In addition, in the current study we constructed a risk model based on the constructed ceRNA network, and the model achieved AUC values of 0.60, 0.63, and 0.64 at 1, 3, and 5 years, and the model well divided the patients into two groups of high and low risk. 6 immune cells in both groups: T cells, CD8+ T cells, cytotoxic lymphocytes, NK cells, B lymphocytes, myeloid dendritic cells, and two stromal cells, endothelial cells and fibroblasts, were significantly different in abundance, suggesting that dysregulation of proliferation of these eight cell types may be associated with breast cancer patients.

Along with risk modeling we also identified three risk genes with risk ratios greater than 1, LRP1, ZIC1, and AXIN2. LRP1 (low-density lipoprotein receptor-related protein 1) is a large multifunctional endothelial cell surface receptor that recognizes numerous ligands and is associated with the regulation of multiple LRP1 has been shown to regulate tumor growth and progression in a variety of tumors and has been identified as a hub in a network of biomarkers used to predict clinical outcomes in multiple cancers [51], a conclusion supported by our study.ZIC1 (zinc finger of the cerebellum 1) is a multi-tumor It inhibits cancer cell growth by reducing Akt and Erk phosphorylation [52], and the high expression of ZIC1 has become a potential biomarker for good prognosis [53], and it has been demonstrated that ZIC1 expression is decreased in breast cancer samples as well as in breast cancer cell lines, while elevated ZIC1 expression inhibited the growth of breast cancer cells and xenograft tumors, while significantly downregulating survivin expression through the Akt/mTOR/P70S6K pathway [54], which presents opposite results to our results, and breast carcinogenesis may be associated with ZIC1 dysregulation, and the specific mechanism of ZIC1 in tumorigenesis still needs further study.AXIN2 (axis formation inhibitor 2) is an important regulator of the Wnt/β-catenin signaling pathway and is involved in cellular functions such as cell proliferation, cell chemotaxis, migration, and apoptosis [55]. Although AXIN2 is known as an oncogene, recent studies have reported that AXIN2 plays an oncogene role in cancers such as colorectal, hepatocellular and gastric cancers [56], and in our study, AXIN2 presented the same situation in breast cancer.

Similarly, when the joint model and clinical features were analyzed, it was found that the risk model was independent of the traditional clinical features and the clinical features with prognostic ability, lncRNA and risk model were used to construct the nomogram, and it was found that the model had the highest percentage score, but the nomogram had better predictive ability, which could provide a theoretical basis for the later multi-target and multi-mechanism combination of breast cancer.

Previous research shows that the patients with estrogen receptor (ER)-positive tumors, tumors with high CXCL10 level showed improved effect of tamoxifen treatment in terms of local recurrence-free survival and CXCL10 can be used as a prognosis biomarker for tamoxifen [57].

Until now, the pathogenesis of BC is still unclear, and there are few effective drugs approved by the FDA for clinical treatment, which poses a great challenge for future research. On the whole, the current analysis of drug targets is mostly based on existing databases, and the clinical experimental data is lacking. To further clarify the pathogenesis of BC and effective targets for treatment, more advanced biological technologies need to be used in basic and clinical research on BC and other tumor-related drug targets. This article mainly studies the related targets of BC from the theoretical level. Therefore, experiment is necessary to further verify the accuracy of prediction results.

Data Availability Statement: The datasets provided in this study can be found in an online repository.

Conflicts of Interest: The authors declare that the study was conducted without any commercial or financial relationships that could be interpreted as potential conflicts of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Funding: This study was supported by Chongqing Natural Science Foundation (cstc2021jcyj-msxmX0834) and the Research Program of Chongqing University of Posts and Telecommunications.

Acknowledgments

K.X.S, Y.F.X. and Y.N.Z conceived and designed the experiments. H.X.R., Y.Z., Z.L.L. and Z.Y.X. assisted the experiments. K.X.S, Y.F.X. and Y.N.Z discussed the results and contributed to the manuscript. Y.F.X. and Y.N.Z wrote the manuscript. All authors have read and agreed to the published version of the manuscript.