Introduction

Breast cancer is the second most frequent and the fifth cause of cancer-associated mortality1. This type of cancer is associated with dysregulation of several genes (including both coding and non-coding ones) and signaling pathways2. Breast cancer is a molecularly heterogeneous disorder which is classified to five subtypes including luminal A, luminal B, basal-like, HER2-enriched and normal-like. This classification is based on the presence/ abundance of estrogen receptor (ER), progesterone receptor (PR), HER2 and Ki673,4. However, several recent studies have indicated significance of other genes and signaling pathways in determination of overall survival (OS) of patients2,5. Among the recently appreciated genes in this regard are long non-coding RNAs (lncRNAs)6. These transcripts are involved in the regulation of fundamental cell survival pathways and have functional interactions with proteins and other non-coding RNAs that participate in the pathogenesis of breast cancer7. Identification of such networks is an important step towards design of targeted therapies in breast cancer.

In the current study, we have retrieved data of two microarray datasets (GSE65194 and GSE45827) from the NCBI Gene Expression Omnibus database (GEO). R package was used for identification of differentially expressed genes (DEGs), assessment of gene ontology (GO) and pathway enrichment evaluation. The DEGs were integrated to construct a protein–protein interaction (PPI) network. Next, hub genes were recognized using the Cytoscape software and lncRNA–mRNA co-expression analysis was performed to evaluate the potential roles of lncRNAs.

Methods

In this study, we used a system biology approach for data mining of two microarray datasets of normal and malignant breast tissue (GSE65194 and GSE45827). We aim to identify differentially expressed genes (DEGs) and lncRNAs and construct an mRNA–lncRNA network based on co-expression analysis. Figure 1 shows summary of the steps accomplished in bioinformatics strategy.

Figure 1
figure 1

Study design flowchart.

Full size image

Gene expression profile data collection

Two gene expression profiles associated with breast cancer (GSE65194 and GSE45827) were obtained from the NCBI Gene Expression Omnibus database (GEO, https://www.ncbi.nlm.nih.gov/geo/). A chip-based platform GPL570 (HG-U133_Plus_2) Affymetrix Human Genome U133 Plus 2.0 Array was applied for both datasets. The GSE65194 included 130 breast cancer samples (41 Triple negative, 30 HER2 positive, 29 Luminal A, 30 Luminal B) and 11 normal breast tissue samples8. Similarly, the GSE45827 contained 130 tumor tissue specimen (41 Triple negative, 30 HER2 positive, 29 Luminal A, 30 Luminal B) as well as 11 normal tissues samples9.

Data preprocessing and DEGs identification

All raw data files were subjected to quantile normalization and background correction using Robust Multichip Average (RMA)10. RMA is an effective tool in the affy Bioconductor package for both mRNA and lncRNA profiling data11. Quality Control assessment was done with AgiMicroRna Bioconductor Package12. We conducted a dimensional reduction analysis by performing Principal component analysis (PCA)13 with the purpose of finding similarities between each group of samples using ggplot2 package of R software14. The linear models for microarray data (LIMMA) R package15 in Bioconductor (https://www.bioconductor.org/)16 were used to perform differential expression gene analysis (DEGA) between breast cancer and normal breast samples. The Student's t-test was applied and DEGs with false discovery rate (FDR) < 0.01 and a |log2FC (fold change)|> 2 were screened.

Functional enrichment analysis

To identify the role of DEGs in breast cancer, KEGG Pathway and GO function enrichment analysis in 3 functional ontologies namely biological process (BP), cellular component (CC) and molecular function (MF) were performed using the DAVID system (https://david.ncifcrf.gov/). The adjusted P < 0.05 was considered as statistically significant17.

PPI network construction, cluster analysis and key gene identification

To predict interactive relationships among common DEGs encoding proteins, we constructed a PPI network using online STRING database (https://string-db.org/)18. The minimum interaction score > 0.4 was required to construct the PPI network. Cytoscape software version 3.7.1 (https://www.cytoscape.org/) was applied to visualize the PPI networks and analyze the hub genes19. We used Molecular COmplex DEtection (MCODE) algorithm (version 1.5.1) to find PPI subnetwork and the highly interconnected clusters within the PPI network. MCODE is a Cytoscape plug-in in which we set maximum depth = 100, node score = 0.2, and K-core = 2 as threshold parameters20. CytoHubba (version 1.6)21 and CytoNCA (version 2.1.6)22 are two other plug-in in which provide multiple algorithms to detect hub genes in the network. In addition, identified key genes were selected for additional expression analysis on 1104 cancer and 113 normal samples from the TCGA project in The Encyclopedia of RNA Interactomes (ENCORI) database (https://starbase.sysu.edu.cn/panCancer.php). Pearson correlation coefficient was assessed between hub genes. The correlation coefficients were also checked on TCGA dataset by using Gene Expression Profiling Interactive Analysis (GEPIA) database (https://gepia.cancer-pku.cn/).

Prediction of lncRNAs function

LncRNA–mRNA co-expression analysis was performed to evaluate the potential roles of lncRNAs. The full list of lncRNA genes with approved HUGO Gene Nomenclature Committee (HGNC) symbols was downloaded from (https://www.genenames.org/)23. The list of lncRNA gene names was compared to our dataset gene symbols and overlapped genes were chosen. Then, differentially expressed lncRNAs were selected according to (|logFC|) > 0.5 and the adjusted P value < 0.01 cutoff criteria. The reason for application of easier selection criteria was the lower expression level of lncRNAs compared with mRNAs. Then, the Pearson correlation coefficient was calculated between the differentially expressed lncRNA and 2 key protein-coding genes that were obtained from the previous steps based on functional annotation and co-expression analysis (MAD2L1 and CCNA2) in our dataset. LncRNAs with correlation coefficients higher than 0.6 or lower than − 0.6 were chosen as the lncRNAs that co-expressed with MAD2L1 and CCNA2. In order to uncovering the importance of these candidate genes in different molecular subtypes of breast cancer, the expression of these genes was also examined in four breast cancer subtypes, including luminal A, luminal B, basal-like and HER2-enriched.

Survival analysis

Survival analysis was carried out on these candidate hub genes to check out their effects on breast cancer survival. Recurrence free survival (RFS) analysis and overall survival (OS) analysis were performed based on expression data from 6234 breast cancer patients by Kaplan Meier plotter (kmplot.com/) that can evaluate the effect of gene expression on survival in 21 cancer types24. We split patients by Mean. In other words, the groups were divided with low expression level and high expression level based on Mean in the survival analysis. The hazard ratio was calculated for both RFS and OS and the P value was determined applying log-rank tests.

Results

DEGs screening

Before performing differentially expressed gene analysis (DEGA), background correction and normalization were done and we removed batch effect. We used AgiMicroRna Bioconductor Package for Quality Control assessment. Degradation plots which indicate the quality of RNA hybridization along the probe sets was drawn and the RNA quality was good. Furthermore, box plots for gene expression data were created to assess the distribution of data after normalization. In the box plots the different arrays had the similar median expression level. This result indicated correction was performed properly. Additionally, a PCA plot was drawn to illustrate the spatial distribution of the samples before and after batch effect correction (Supplementary Figure S1a). Principal components analysis (PCA) provides information about the structure of the analyzed dataset. It can be used to find similarities between samples. We found two samples from the normal group which is spatially far from other normal samples. As a consequence, we removed these two samples. Furthermore, a heatmap was drawn to illustrate the correlation between samples using Pheatmap package of R 3.6.1 software25 (Supplementary Figure S1b). After correction, removing the batch effects and performing data normalization, 887 DEGs including 730 upregulated and 157 downregulated DEGs were screened between breast cancer and normal samples from GSE65194 and GSE45827 according to |logFC|> 2 and FDR < 0.01 as cut-off criteria. The list of upregulated and downregulated DEGs are indicated in Supplementary Tables S1 and S2, respectively. Figure 2a is a Venn diagram which illustrates the overlap between 2 datasets. Moreover, to visualize the overall gene expression levels of the DEGs, a Volcano plot was created with log2 FC score and log10 P values in R software (Fig. 2b).

Figure 2
figure 2

(a) A Venn diagram of 887 overlapping DEGs in GSE65194 and GSE45827. (b) Volcano plot of significant DEGs with |logFC|> 2. (c) Volcano plot of significant differentially expressed lncRNA with |logFC|> 0.5.

Full size image

KEGG and GO enrichment analysis

To further examine the role of common DEGs in breast cancer, we performed GO and KEGG pathway enrichment analysis26,27. We found 10 dysregulated pathways based on the adjusted P < 0.05. Up-regulated DEGs were enriched in six pathways including ‘Cell cycle’, ‘Oocyte meiosis’ and ‘Focal adhesion’. Down-regulated DEGs were enriched in five pathways including ‘Peroxisome-proliferator-activated receptors (PPAR) signaling pathway’, ‘Metabolism of xenobiotics by cytochrome P450’, ‘Adipocytokine signaling pathway’ and ‘Cytokine-cytokine receptor interaction’ pathways (Fig. 3a). The results for each GO functional analysis are presented in Fig. 3b–d. The genes enriched in KEGG pathway and GO enrichment analysis have shown in Supplementary Tables S3–S6.

Figure 3
figure 3figure 3

KEGG and GO enrichment analysis. (a) KEGG pathways (based on KEGG PATHWAY database26,27). (b) GO for DEGs, Biological process. (c) GO for DEGs, Cellular component. (d) GO for DEGs, Molecular function.

Full size image

PPI network construction and module analysis

The interactive information among DEGs and the PPI network was obtained using the STRING online database. Among the total of common DEGs, 887 DEGs (730 up-regulated and 157 down-regulated) were filtered into the PPI network with 887 nodes and 10,398 edges, at a combined score > 0.4. Finally, Genes with a combined score > 0.9 were selected as key DEGs to be imported into Cytoscape. The Cytoscape software was applied to evaluate the interactive relationships between the candidate proteins. Afterward, two clusters consist of 65 nodes and 23 nodes were screened with a cut-off k-score = 12 depend on the MCODE scoring system (Supplementary Figure S2). The CytoNCA and the CytoHubba are two Cytoscape plug-in for centrality analysis and give us some insight into the most influential nodes or edges in a network. We ran CytoHubba application and extracted data from four calculations methods (EPC, MCC, MNC, and Stress). The top 100 nodes ranked by these four methods were selected (Supplementary Table S7). Moreover, four algorithms from CytoNCA application (Degree, Eigenvector, Betweenness, and Closeness) were employed and the top 100 nodes based upon these four approaches were obtained (Supplementary Table S7). Besides, a Venn diagram was created to identify the significant hub genes that are similar between all groups. The result of Venn diagram is mentioned in Supplementary Table S8. Eventually, through overlapping analysis, we identified a list of 26 key genes most of them belonged to MCODE cluster 1 (Supplementary Table S8). Since highly interconnected proteins in a network accumulate in a cluster, we chose only 20 genes from our list that belonged to cluster 1 (Table 1). All the selection steps are illustrated in Fig. 4a.

Table 1 Key differentially expressed genes acquired by centrality analysis.
Full size table
Figure 4
figure 4

(a) A Venn diagram of 26 overlapping genes between different calculation methods of Cytohubba and CytoNCA. (b) Heatmap correlation plot for 20 candidate genes (depicted using Pheatmap package of R software25).

Full size image

Key genes functional annotation and co-expression analysis

GO enrichment and KEGG pathway analysis on these 20 genes indicated that four pathways were enriched, including cell cycle, progesterone-mediated oocyte maturation, oocyte meiosis, and p53 signaling pathway. CCNA2, CDK1, MAD2L1, and CCNB1 were significantly enriched in some biological aspects such as cell cycle, mitosis, nuclear division, M phase, cell cycle and progesterone-mediated oocyte maturation pathways. In particular, by checking the expression data of 1104 cancer and 113 normal samples from the TCGA project in ENCORI database, we found that these four genes showed strong expression in the breast cancer specimens as compared to their expression in normal breast tissue ( Including : MAD2L1, Fold change: 4.28, Adjusted P value: 1.4e−70; CCNA2, Fold change: 6.88, Adjusted P value: 3.2e−91; CCNB1, Fold change: 5.63, Adjusted P value: 1.8e−111; CDK1, Fold change: 8.54, Adjusted P value: 5.3e−121). Additionally, we calculated the Pearson correlation for these 20 candidate genes and found a strong and significant correlation between them (Supplementary Table S9, Fig. 4b). Interestingly, CCNA2 and MAD2L1 which are two important genes in the cell cycle pathway and some crucial biological processes related to cell division, were highly correlated genes with a correlation coefficient higher than 0.9 in our analysis. Furthermore, these two genes correlation in TCGA dataset in the GEPIA database was consistent with our analysis.

Identification of differentially expressed lncRNAs and co-expression analysis

After downloading the list of lncRNA genes from HGNC database, lncRNAs genes symbols were extracted from the GSE65194 and GSE45827. A total of 334 lncRNA probes were identified in these two datasets by using this approach. Finally, 159 lncRNAs probe ID with |logFC|> 0.5 and adjusted P value < 0.01 among 20 normal samples and 258 breast tissue samples were picked out (Fig. 2c). Among these lncRNAs, 77 lncRNAs were up-regulated (Supplementary Table S10) and 80 lncRNAs were down-regulated (Supplementary Table S11) in breast cancer. We calculated Pearson correlation coefficient between differentially expressed lncRNAs and MAD2L1 and CCNA2 based on their expression value. LncRNA with Pearson correlation coefficient ≥ 0.6 or ≤  − 0.6 were selected as key lncRNA which co-expressed with MAD2L1 and CCNA2. Totally, 12 lncRNAs meet this criterion (Table 2). Additionally, Table 3 indicates the expression of these genes in four breast cancer subtypes. Our selected genes appear to be more important in more aggressive sub-types (basal-like and HER-enriched). However, deregulation of CARMN, PRINS and MEG3 may be crucial in all subtypes of breast cancer.

Table 2 Key lncRNAs which co-expressed with MAD2L1 and CCNA2.
Full size table
Table 3 Relative expression of our candidate genes in different molecular subtypes of breast cancer and healthy breast tissue in GSE65194 and GSE45827.
Full size table

Survival analysis of candidate hub genes

Associations between expression of candidate hub genes and RFS and OS of the breast cancer patients were evaluated using KM method to estimate the prognostic importance of the hub genes in our study. The results indicated that low expression of MAD2L1, CCNA2 and NCK1-DT lead to higher RFS rate than high expression. Inversely, high expressions of MEG3, RAD51-AS1, PRINS, LINC01089, LINC02256, FUT8-AS1, LINC01279, CARMN, EPB41L4A-AS1, EIF3J-DT and TNFRSF14-AS1 result in a significantly longer RFS time among patients with breast cancer. The results showed that MAD2L1, CCNA2, RAD51-AS1 and LINC01089 have the most prediction potential based on RFS among all candidate hub genes. Besides, hazard ratio was also calculated for OS. High expression of MAD2L1, CCNA2 and FUT8-AS1 lead to lower OS rate than low expression. On the other hand, low expressions of LINC01279, RAD51-AS1 and CARMN were correlated with significantly worse OS in breast cancer patients. Other candidate hub genes expressions were not significantly relevant to OS (Table 4).

Table 4 Recurrence free survival (RFS) and overall survival (OS) of candidate hub genes.
Full size table

Discussion

In the present study, we used a bioinformatics strategy to identify key genes and signaling pathways in breast cancer pathogenesis with a focus on the role of lncRNAs and their interactions with protein-coding genes. Such interactions can be assessed using experimental approaches which are costly and laborious. Bioinformatics methods for such purpose fall into two groups: strategies that use sequence, structural data and physicochemical features, and methods that are based on network construction. The latter can provide the inherent characteristics of topological configuration of associated biological networks which is often disregarded by the former strategies6. In the present work, we used GPL570 which is a good platform to evaluate the expression level of lncRNAs in tumorigenesis11,28. We identified 730 upregulated and 157 downregulated DEGs between breast cancer and normal samples. Up-regulated DEGs were enriched in ‘Cell cycle’, ‘Oocyte meiosis’ and ‘Focal adhesion’. A previous bioinformatics study using topological characteristics of genes in breast cancer has identified these pathways as hub subnetworks29. The role of these pathways has been acknowledged in the pathogenesis of another hormone related cancer namely prostate cancer30. We also detected down-regulated DEGs were enriched in ‘PPAR signaling pathway’, ‘Metabolism of xenobiotics by cytochrome P450′, ‘Adipocytokine signaling pathway’ and ‘Cytokine-cytokine receptor interaction’ pathways. PPARs are nuclear hormone receptors which participate in modulation of different aspects of tumorigenesis such as cell proliferation, survival and apoptosis31. Xenobiotic metabolizing enzymes are also involved in the tumorigenesis and response of cancer patients to therapeutic options. Integration of expression data of these genes with eQTL data and allele frequency data from the 1000 Genomes project has shown considerable inter-population differences in the related pathways which might influence cancer prognosis and response to treatment32. Adipocytokines can also influence cell proliferation and survival, and malignant phenotypes of breast cancer cells through regulation several cellular and molecular pathways thus aggravating survival of patients33. Cytokine signaling has important functions in formation, proliferation, and migration of breast cancer, thus modulating invasiveness, angiogenesis and metastatic potential of these cells34.

Our in silico analyses revealed that CCNA2, CDK1, MAD2L1 and CCNB1 were significantly enriched in several biological pathways. These four genes showed strong expression in breast cancer samples as compared to their expression in normal breast tissue. Notably, these four genes have been among the top dysregulated genes in small cell lung cancer as revealed by GO, KEGG analysis and construction of PPI network35. Such similarity between these two different types of cancers implies fundamental role of these genes in the carcinogenesis process and potentiates them as therapeutic targets. MAD2L1 form a complex with the APC/C and CDC20 and subsequently stimulate the M-A checkpoint to halt the transition of cell at this stage in the presence of anomalous segregation of chromatin. Yet, over-expression of E2F1 in atypical cells affects the formation of the mentioned complex leading to cell cycle transition even in the presence of abnormal chromosomes36. CDK1/cyclin B is a maturation-promoting factor37 and the checkpoint for G2/M transition38,39, so it is expected to be involved in the process of cell cycle regulation and tumorigenesis. We also identified 12 lncRNAs with significant correlation with MAD2L1 and CCNB1 genes. As expected from KEGG analysis, KM analysis indicated that low expression of MAD2L1, CCNA2 and NCK1-DT lead to higher RFS rate than high expression. Inversely, high expressions of MEG3, RAD51-AS1, PRINS, LINC01089, LINC02256, FUT8-AS1, LINC01279, CARMN, EPB41L4A-AS1, EIF3J-DT and TNFRSF14-AS1 result in a significantly longer RFS time among patients with breast cancer. Additionally, hazard ratio was also calculated for OS. High expression of MAD2L1, CCNA2 and FUT8-AS1 and low expressions of LINC01279, RAD51-AS1 and CARMN were correlated with significantly worse OS in breast cancer patients, while Other candidate hub genes expression were not significantly relevant to OS.

According to previous studies MEG3 is down-regulated in breast cancer tissues40,41,42. Recently, Zhang et al. showed MEG3 ability in promoting breast cancer growth and induction of apoptosis by activating ER stress, NF-κB and p53 pathways in breast cancer cell line43. RAD51-AS1, also known as TODRA is transcribed from upstream of RAD51 in a divergent manner. Gazy et al. identified a conserved E2F1 binding site in the promoter region of RAD51-AS1 and considered this lncRNA as a target gene of E2F1 in breast cancer. RAD51-AS1 negatively regulates RAD51 expression and higher expression of RAD51-AS1 has been associated with a less aggressive tumor phenotype44. PRINS (Psoriasis susceptibility-related RNA Gene Induced by Stress) is a stress induced lncRNA which regulates apoptosis45,46. Min Yu et al. considered PRINS as a HIF-1α dependent lncRNA due to its significant over-expression in hypoxic conditions in renal tubular cells47. Moreover, increased levels of PRINS have been observed in colorectal adenocarcinoma cells. This lncRNA interacts with trefoil factor 3 (TFF3), AKT/PI3K signaling pathway and miR-491-5p48. PRINS levels were down-regulated in MCF-7 and MDA-MB-231 cell lines following exposure to the apoptotic and anti-proliferative agent CCT13769049. LINC01089 (also known as LncRNA Inhibiting Metastasis; LIMT) is an EGF regulated lncRNA which is down-regulated in breast cancer tissues and cell lines, especially in aggressive subtypes of breast cancer50. Yuan et al. have reported a significant correlation between low expression of LINC01089 and lymph node metastasis and poor prognosis of breast cancer. LINC01089 is modulates breast tumorigenesis by inhibiting β-catenin transcription and consequently blocking Wnt/β-catenin signaling51. LINC02256 (ENSG00000261064) is a validated novel long intergenic non-protein coding RNA with 2 transcripts which is located on 15q13.3. Based on GTEx (Release v6) results, it has ubiquitous expression in breast and other tissues. Potential contribution of this lncRNA in breast cancer should be evaluated in future studies. FUT8‐AS1 was up-regulated in endometrioid endometrial cancer patients in association with poor survival52. According to another TCGA data mining study on glioblastoma, FUT8-AS1 over-expression has been associated with poor patients outcomes53. LINC01279 was significantly upregulated in patients with endometriosis. Based on Liu et al. study, there is a strong association between this lncRNA and cell cycle-dependent kinase-14 and CXC motif chemokine ligand-12. Hence, LINC01279 might contribute in the pathogenesis of endometriosis54. In another bioinformatics analysis of differential gene expression in breast cancer LINC01279 was significantly down-regulated55. However, further studies should be done to elucidate its function in breast cancer. CARMN (also known as MiR143HG) is recognized as a tumor suppressor in bladder cancer. Xie et al. observed down-regulation of CARMN in bladder cancer tissues compared with normal tissues. Moreover, there was an association between CARMN over-expression and a high survival rate in bladder cancer patients. CARMN/miR‐1275/AXIN2 axis takes part in bladder tumorigenesis by interacting with the Wnt/β‐catenin pathway56. Furthermore, there is an association between down-regulation of CARMN and poor survival in endometrial carcinoma57. CARMN was significantly down-regulated in hepatocellular carcinoma (HCC) tissues and cells. Over-expression of CARMN was associated with good prognosis. Generally, this gene contributes to development and progression of HCC by blocking the MAPK and Wnt signaling pathways58. EPB41L4A-AS1 (also known as TIGA1) is a p53-regulated gene. EPB41L4A-AS1 was down-regulated in many human cancers in correlation with poor prognosis. EPB41L4A-AS1 acts as a repressor of the Warburg effect and is involved in cancer metabolic reprogramming59. In a recent study in early stage breast cancer, significant down-regulation of EPB41L4A-AS1 was observed in tumor tissues60. EIF3J-DT is a novel lncRNA with no publication reporting its biological function in breast cancer up to now. Based on one study in HCC, EIF3J-DT might have potential prognostic value61. In addition, EIF3J-DT regulates multi-drug resistance by interacting with autophagy in gastric cancer62. Based on He et al. study, TNFRSF14-AS1 might have a prognostic value in breast cancer but this result needs to further confirmation63. Based on the available literature, the identified lncRNAs in the current study has putative roles in the pathogenesis of breast cancer and other types of cancer.

Taken together, in the present study, we intended to introduce a precise method to discover and prioritize the most probable candidate genes involved in breast cancer. Gene expression analysis in different molecular subtypes indicated the importance of our chosen genes in more aggressive subtypes. On the other hand, CARMN, PRINS and MEG3 probably have an important role in pathogenesis of all subtypes of breast cancer. We also added several evidences from literature regarding the role of candidate genes in the pathogenesis of cancer. Although this study provides some impressive evidence for future differential expression studies in breast cancer, the limitation of this study is lack of experimental evaluation of the candidate genes. our in silico method identified a number of hub genes and related lncRNAs which are possibly involved in the pathogenesis of breast cancer and patients' prognosis, so can be used as therapeutic targets or biomarkers for this malignancy.