An in-silico method leads to recognition of hub genes and crucial pathways in survival of patients with breast cancer

Breast cancer is a highly heterogeneous disorder characterized by dysregulation of expression of numerous genes and cascades. In the current study, we aim to use a system biology strategy to identify key genes and signaling pathways in breast cancer. We have retrieved data of two microarray datasets (GSE65194 and GSE45827) from the NCBI Gene Expression Omnibus database. R package was used for identification of differentially expressed genes (DEGs), assessment of gene ontology and pathway enrichment evaluation. The DEGs were integrated to construct a protein–protein interaction 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. Finally, the clinical importance of the obtained genes was assessed using Kaplan–Meier survival analysis. In the present study, 887 DEGs including 730 upregulated and 157 downregulated DEGs were detected between breast cancer and normal samples. By combining the results of functional analysis, MCODE, CytoNCA and CytoHubba 2 hub genes including MAD2L1 and CCNB1 were selected. We also identified 12 lncRNAs with significant correlation with MAD2L1 and CCNB1 genes. According to The Kaplan–Meier plotter database MAD2L1, CCNA2, RAD51-AS1 and LINC01089 have the most prediction potential among all candidate hub genes. Our study offers a framework for recognition of mRNA–lncRNA network in breast cancer and detection of important pathways that could be used as therapeutic targets in this kind of cancer.


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.
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 ( 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 data 11 . Quality Control assessment was done with AgiMicroRna Bioconductor Package 12 . 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 software 14 . The linear models for microarray data (LIMMA) R package 15 in Bioconductor (https ://www.bioco nduct or.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 .ncifc rf.gov/). The adjusted P < 0.05 was considered as statistically significant 17 .
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 ://strin g-db.org/) 18  www.nature.com/scientificreports/ 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 parameters 20 . Cyto-Hubba (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 ://starb ase.sysu.edu.cn/panCa ncer.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 .cance r-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.genen ames.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 HER2enriched.

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 types 24 . 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 software 25 (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).

KEGG and GO enrichment analysis.
To further examine the role of common DEGs in breast cancer, we performed GO and KEGG pathway enrichment analysis 26,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-proliferatoractivated 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 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 www.nature.com/scientificreports/ 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.

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   . 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. www.nature.com/scientificreports/ 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.
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).

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 strategies 6 .
In the present work, we used GPL570 which is a good platform to evaluate the expression level of lncRNAs in tumorigenesis 11,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 subnetworks 29 . The role of these pathways has been acknowledged in the pathogenesis of another hormone related cancer namely prostate cancer 30 . We also detected down-regulated DEGs were enriched in 'PPAR signaling pathway' , 'Metabolism of xenobiotics by cytochrome P450′, ' Adipocytokine signaling pathway' and  www.nature.com/scientificreports/ '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 apoptosis 31 . 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 treatment 32 . 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 patients 33 . Cytokine signaling has important functions in formation, proliferation, and migration of breast cancer, thus modulating invasiveness, angiogenesis and metastatic potential of these cells 34 .
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 network 35 . 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 chromosomes 36 . CDK1/cyclin B is a maturation-promoting factor 37 and the checkpoint for G2/M transition 38,39 , so it is expected to be involved in  www.nature.com/scientificreports/ 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 tissues [40][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 line 43 . 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 phenotype 44 . PRINS (Psoriasis susceptibility-related RNA Gene Induced by Stress) is a stress induced lncRNA which regulates apoptosis 45,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 cells 47 . 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-5p 48 . PRINS levels were down-regulated in MCF-7 and MDA-MB-231 cell lines following exposure to the apoptotic and anti-proliferative agent CCT137690 49 . 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 cancer 50 . 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 signaling 51 . LINC02256 (ENSG00000261064) is a validated novel long intergenic nonprotein 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 survival 52 . According to another TCGA data mining study on glioblastoma, FUT8-AS1 over-expression has been associated with poor patients outcomes 53 . 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 endometriosis 54 . In another bioinformatics analysis of differential gene expression in breast cancer LINC01279 was significantly down-regulated 55 . 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 pathway 56 . Furthermore, there is an association between down-regulation of CARMN and poor survival in endometrial carcinoma 57 . 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 pathways 58 . EPB41L4A-AS1 www.nature.com/scientificreports/ (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 reprogramming 59 . In a recent study in early stage breast cancer, significant down-regulation of EPB41L4A-AS1 was observed in tumor tissues 60 . 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 value 61 . In addition, EIF3J-DT regulates multi-drug resistance by interacting with autophagy in gastric cancer 62 . Based on He et al. study, TNFRSF14-AS1 might have a prognostic value in breast cancer but this result needs to further confirmation 63 . 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.