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

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 Study design flowchart. Study design flowchart. 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 samples 8 . 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 samples 9 . 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.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)  2 were screened. 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 significant 17 . 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 genes 19 . 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 parameters 20 . 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/ ). 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 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.

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 S1 a). 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 S1 b). 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  2 a 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.  2 b). 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. ( 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. 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-proliferator-activated receptors (PPAR) signaling pathway’, ‘Metabolism of xenobiotics by cytochrome P450’, ‘Adipocytokine signaling pathway’ and ‘Cytokine-cytokine receptor interaction’ pathways (Fig.  3 a). The results for each GO functional analysis are presented in Fig.  3 b–d. The genes enriched in KEGG pathway and GO enrichment analysis have shown in Supplementary Tables S3 – S6 . Figure 3 KEGG and GO enrichment analysis. ( a ) KEGG pathways (based on KEGG PATHWAY database 26 , 27 ). ( b ) GO for DEGs, Biological process. ( c ) GO for DEGs, Cellular component. ( d ) GO for DEGs, Molecular function. KEGG and GO enrichment analysis. ( a ) KEGG pathways (based on KEGG PATHWAY database 26 , 27 ). ( b ) GO for DEGs, Biological process. ( c ) GO for DEGs, Cellular component. ( d ) GO for DEGs, Molecular function. 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.  4 a. Table 1 Key differentially expressed genes acquired by centrality analysis. Gene logFC adj. P .Val MCODE Centrality analysis by CytoNCA Centrality analysis by CytoHubba MCODE Score Betweenness Closeness Degree Eigenvector EPC MCC MNC Stress CDK1 3.729327 3.03E−28 46.020339 6576.003525 0.551637 131 0.1437205 55.3 9.22E+13 158 356,658 CCNB1 3.115417 1.59E−27 46.020339 3303.18478 0.534146 116 0.1417335 55.9 9.22E+13 133 213,424 CCNA2 2.219921 2.57E−14 46.020339 2534.615585 0.52019 105 0.1342352 53.4 9.22E+13 122 163,390 CDC20 2.756389 1.63E−14 46.020339 2758.724299 0.496036 103 0.1345755 51.8 9.22E+13 116 116,042 MAD2L1 3.346366 4.29E−26 46.020339 1131.260005 0.48079 94 0.1333797 52.9 9.22E+13 110 72,556 KIF11 3.207367 4.39E−28 46.020339 1464.21417 0.487751 92 0.1298102 50.4 9.22E+13 104 95,940 CENPA 2.75168 8.64E−15 47.094949 1704.200531 0.481319 92 0.1282006 50.3 9.22E+13 100 74,860 PCNA 2.453675 1.03E−37 43.857039 2862.973827 0.489385 92 0.1150848 47.9 9.22E+13 106 108,266 EZH2 2.942149 1.48E−27 42.305272 3357.588136 0.511085 91 0.1046336 47.1 9.22E+13 112 243,314 KIF23 2.687969 7.39E−20 46.020339 1657.240201 0.493799 89 0.1254939 49 9.22E+13 97 99,424 TOP2A 4.595822 5.72E−29 46.020339 1122.678983 0.493243 88 0.1293264 51.4 9.22E+13 104 94,808 UBE2C 3.062068 1.20E−22 46.020339 1455.353604 0.466951 88 0.1219277 48.8 9.22E+13 101 67,434 BIRC5 2.792394 5.43E−15 46.020339 1678.399139 0.493799 88 0.1287475 48.9 9.22E+13 99 108,288 KIF2C 2.389307 1.51E−15 46.94026 1182.010234 0.487751 88 0.1259678 48 9.22E+13 97 85,072 RRM2 4.481084 1.30E−31 46.020339 1745.829126 0.497727 86 0.1228496 50.9 9.22E+13 101 110,592 RACGAP1 2.400852 1.87E−23 46.020339 1361.211732 0.493799 86 0.1225891 49.8 9.22E+13 94 93,686 KIF4A 2.206078 2.43E−13 46.94026 1008.406749 0.478689 80 0.1192849 46 9.22E+13 88 65,850 KPNA2 3.428273 1.09E−42 46.880102 961.5679185 0.47454 78 0.111362 46.1 9.22E+13 86 91,340 TYMS 2.878076 1.89E−21 47.774118 1340.60442 0.473514 76 0.1156782 48.1 9.22E+13 92 99,082 RRM1 2.268661 2.13E−25 40.445411 1145.682149 0.45768 69 0.0981711 43.2 9.22E+13 80 69,510 logFC, log2 fold change; adj. P .Val, adjusted P value; EPC, Edge Percolated Component; MCC, Maximal Clique Centrality; MNC, Maximum Neighborhood Component. 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 software 25 ). Key differentially expressed genes acquired by centrality analysis. logFC, log2 fold change; adj. P .Val, adjusted P value; EPC, Edge Percolated Component; MCC, Maximal Clique Centrality; MNC, Maximum Neighborhood Component. ( 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 software 25 ). 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.  4 b). 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. 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.  2 c). 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. Symbol MA2L1 CCNA2 Pearson correlation P value Pearson correlation P value MEG3 − 0.68 < 2.2e−16 − 0.73 < 2.2e−16 LINC01279 − 0.66 < 2.2e−16 − 0.68 < 2.2e−16 EIF3J-DT − 0.64 < 2.2e−16 − 0.68 < 2.2e−16 LINC01089 − 0.64 < 2.2e−16 − 0.67 < 2.2e−16 RAD51-AS1 − 0.64 < 2.2e−16 − 0.63 < 2.2e−16 TNFRSF14-AS1 − 0.61 < 2.2e−16 − 0.63 < 2.2e−16 LINC02256 − 0.60 < 2.2e−16 − 0.62 < 2.2e−16 PRINS − 0.60 < 2.2e−16 − 0.62 < 2.2e−16 EPB41L4A-AS1 − 0.60 < 2.2e−16 − 0.62 < 2.2e−16 CARMN − 0.60 < 2.2e−16 − 0.60 < 2.2e−16 FUT8-AS1 − 0.60 < 2.2e−16 − 0.60 < 2.2e−16 NCK1-DT 0.65 < 2.2e−16 0.66 < 2.2e−16 CCNA2 0.93 < 2.2e−16 1 < 2.2e−16 MAD2L 1 1 < 2.2e−16 0.936209 < 2.2e−16 Table 3 Relative expression of our candidate genes in different molecular subtypes of breast cancer and healthy breast tissue in GSE65194 and GSE45827 . Symbol Basal-like HER2-enriched Luminal A Luminal B logFC adj. P .Val logFC adj. P .Val logFC adj. P .Val logFC adj. P .Val MAD2L1 4.4 3.02E−61 3.57 6.58E−44 1.709 3.53E−13 3.297 4.45E−37 CCNA2 3.91 9.31E−50 2.78 1.50E−28 0.644 0.00823 1.986 1.09E−15 RAD51-AS1 − 1.6 1.47E−33 − 1.6 1.57E−30 − 0.9 6.18E−12 − 1.27 5.54E−21 LINC01089 − 1 2.44E−18 − 0.8 6.34E−11 − 0.32 0.00964 − 0.71 5.79E−09 EIF3J-DT − 0.9 1.02E−13 − 0.8 1.66E−10 − 0.04 0.77445 − 0.46 0.00031 LINC02256 − 0.9 1.13E−21 − 0.9 3.56E−22 − 0.44 3.20E−06 − 0.57 1.22E−09 TNFRSF14-AS1 − 0.8 8.62E−26 − 0.8 6.44E−25 − 0.34 6.84E−06 − 0.5 7.35E−11 CARMN − 3.2 6.13E−22 − 3 4.42E−18 − 1.91 4.53E−08 − 2.41 5.23E−12 EPB41L4A-AS1 − 1.6 1.78E−20 − 1.3 2.99E−13 − 0.7 0.0002 − 1.24 3.19E−11 LINC01279 − 2.1 2.58E−08 − 0.9 0.022916 0.484 0.27445 − 0.66 0.10494 MEG3 − 2.8 7.98E−40 − 2.1 2.33E−25 − 1.55 5.77E−14 − 1.91 8.24E−20 FUT8-AS1 − 0.8 2.61E−31 − 0.8 7.15E−26 − 0.49 5.63E−11 − 0.55 1.80E−13 PRINS − 1.2 4.18E−44 − 1.2 2.64E−38 − 0.92 4.31E−24 − 1.09 1.10E−31 NCK1-DT 1.58 1.22E−22 0.94 4.64E−09 0.625 0.00022 0.897 8.21E−08 logFC, log2 fold change; adj. P .Val, adjusted P value. Key lncRNAs which co-expressed with MAD2L1 and CCNA2. Relative expression of our candidate genes in different molecular subtypes of breast cancer and healthy breast tissue in GSE65194 and GSE45827 . logFC, log2 fold change; adj. P .Val, adjusted P value. 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. Gene name Multivariate analysis for RFS Univariate analysis for RFS Multivariate analysis for OS Univariate analysis for OS Probe ID Symbole HR CI logrank P HR CI logrank P HR CI logrank P HR CI logrank P 210794_s_at MEG3 0.8 0.71–0.89 0.0001 0.73 0.65–0.81 1.30e−08 0.93 0.74–1.17 0.5396 0.82 0.66–1.02 0.0713 227061_at LINC01279 0.74 0.63–0.87 0.0003 0.69 0.59–0.81 4.2e−06 0.73 0.53–1.01 0.0553 0.72 0.53–0.99 0.0391 235124_at EIF3J-DT 0.57 0.48–0.68 0 0.5 0.42–0.58  < 1e−16 0.76 0.54–1.07 0.1189 0.62 0.45–0.85 0.0029 226369_at LINC01089 0.5 0.42–0.58 0 0.46 0.39–0.54  < 1e−16 0.84 0.6–1.17 0.2993 0.72 0.53–0.99 0.0433 1560081_at RAD51-AS1 0.45 0.38–0.53 0 0.41 0.35–0.49  < 1e−16 0.7 0.5–0.99 0.0427 0.61 0.45–0.84 0.0023 232190_x_at TNFRSF14-AS1 0.63 0.54–0.74 0 0.58 0.5–0.68 1.4e−11 0.86 0.62–1.18 0.3494 0.76 0.55–1.04 0.0868 234423_x_at LINC02256 0.63 0.53–0.74 0 0.58 0.49–0.68 5.2e−12 1.26 0.91–1.75 0.1708 1.06 0.78–1.46 0.6976 216051_x_at PRINS 0.7 0.63–0.78 0 0.69 0.62–0.77 2.4e−11 0.91 0.73–1.13 0.3928 0.86 0.69–1.06 0.1658 225698_at EPB41L4A-AS1 0.77 0.65–0.91 0.0023 0.65 0.55–0.76 3.3e−08 0.91 0.65–1.28 0.5966 0.75 0.55–1.03 0.0781 1558828_s_at CARMN 0.62 0.53–0.73 0 0.59 0.51–0.69 5.3e−11 0.61 0.44–0.84 0.0028 0.56 0.41–0.77 0.0003 242889_x_at FUT8-AS1 0.8 0.68–0.93 0.0046 0.75 0.64–0.88 0.00028 1.74 1.25–2.43 0.0011 1.54 1.11–2.12 0.0083 228799_at NCK1-DT 1.06 0.91–1.25 0.4498 1.18 1.01–1.38 0.035 0.95 0.68–1.32 0.751 1.09 0.8–1.49 0.5982 203418_at CCNA2 1.24 1.09–1.42 0.0011 1.84 1.64–2.05  < 1e−16 1.36 1.05–1.77 0.0193 1.55 1.25–1.93 5.0e−05 203362_s_at MAD2L1 1.66 1.47–1.87 0 1.86 1.67–2.08  < 1e−16 1.8 1.42–2.28 0 2.02 1.62–2.51 1.8e−10 HR, hazard ratio; CI, confidence interval; RFS, recurrence free survival; OS, overall survival. A multivariate analysis was performed for MKI67, ESR1, and HER2 (ERBB2). Recurrence free survival (RFS) and overall survival (OS) of candidate hub genes. HR, hazard ratio; CI, confidence interval; RFS, recurrence free survival; OS, overall survival. A multivariate analysis was performed for MKI67, ESR1, and HER2 (ERBB2).

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 ‘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 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 – 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 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 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 (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.

Breast cancer is the second most frequent and the fifth cause of cancer-associated mortality 1 . This type of cancer is associated with dysregulation of several genes (including both coding and non-coding ones) and signaling pathways 2 . 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 Ki67 3 , 4 . However, several recent studies have indicated significance of other genes and signaling pathways in determination of overall survival (OS) of patients 2 , 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 cancer 7 . 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.

Supplementary Information 1. Supplementary Information 2. Supplementary Information 3. Supplementary Information 4. Supplementary Information 5. Supplementary Information 6. Supplementary Information 7. Supplementary Information 8. Supplementary Information 9. Supplementary Information 10. Supplementary Information 11. Supplementary Information 12. Supplementary Information 1. Supplementary Information 2. Supplementary Information 3. Supplementary Information 4. Supplementary Information 5. Supplementary Information 6. Supplementary Information 7. Supplementary Information 8. Supplementary Information 9. Supplementary Information 10. Supplementary Information 11. Supplementary Information 12.