Differentially expressed miRNAs screening
We screened out 12 common miRNAs by using |log2FC| > 1 and p < 0.05 threshold between 5 miRNA expression arrays, GSE113740, GSE113486, GSE83270, GSE73002, and GSE57897. All of them were upregulated (Figure 1).
Weighted co-expression network construction and key modules identification
GSE113740 (Figure 2A), GSE113486 (Figure 3A), and GSE57897 (Figure 4A) were selected for WGCNA analysis. pickSoftThreshold function specified β = 18 for GSE113740, β = 7 for GSE113486 and β = 7 for GSE57897 choose for reaching scale-free topology fit index (scale free R2 = 0.9). Hierarchical clustering of GSE113740 represented five differently modules, including blue (120 miRNAs), brown (260 miRNAs), grey (503 miRNAs), turquoise (112 miRNAs), and yellow (181 miRNAs) (Figure 2B and Figure 2C). Hierarchical clustering of GSE113486 represented six modules including yellow, blue, brown, green, turquoise and grey (Figure 3B and Figure 3C). Hierarchical clustering of GSE57897 represented six modules including turquoise, blue, yellow, and brown, green and grey (Figure 4B and Figure 4C). The significant association of modules and clinical status (disease status) were identified. Then eigengene and clinical traits correlation was demonstrated on a color-coded table. Also, miRNAs and clinical traits correlation were displayed on scatterplot. We choosed turquoise module of GSE113740 (correlation= 0.76 and P.value <1e-200) (Figure 2D). For GSE113486, turquoise module was selected (correlation = 0.57 and P.value =2.2e-16) (Figure 3D) and for GSE57897, yellow module was selected (correlation = 0.21 and P.value = 0.017) (Figure 4D). At last, we merge the common miRNAs of selected modules with DEMs and 11 core miRNAs were identified, Hsa-miR-151a-5p, Hsa-miR-185-5p, Hsa-miR30c-2-3p, Hsa-miR-22-3p, Hsa-miR-181c-5p, Hsa-miR-501-3p, Hsa-miR-34a-5p, Hsa-miR-532-5p, Hsa-miR-450b-5p, Hsa-miR-1307-3p, and Hsa-miR-30a-5p. All of them were upregulated and considered for other analysis.
Survival analysis of core miRNAs
We combined the results of co-expression networks and DEMs and then extracted the survival plot of 11 core miRNAs including Hsa-miR-151a-5p, Hsa-miR-185-5p, Hsa-miR30c-2-3p, Hsa-miR-22-3p, Hsa-miR-181c-5p, Hsa-miR-501-3p, Hsa-miR-34a-5p, Hsa-miR-532-5p, Hsa-miR-450b-5p, Hsa-miR-1307-3p, and Hsa-miR-30a-5p via OncoLnc database. As it is clear from plots 6 out of 11 miRNAs show meaningful P.Value, therefore are good candidates for next steps of analysis. 6 core miRNAs are hsa-MiR-450b-5p, hsa-MiR-1307-3p, hsa-MiR-34a-5p, hsa-MiR-151a-5p, hsa-MiR-501-3p, and hsa-MiR-532-5p (Figure 5).
Evaluation of core miRNAs expression using TCGA
To confirm the efficacy of 6 core miRNAs, we detected the expression of them by mirtv database (Figure 6). As it is clear from the figure the expression of six miRNAs including Hsa-miR-151a-5p, Hsa-miR-34a-5p, Hsa-miR-1307-3p, Hsa-miR-450b-5p, Hsa-miR-501-3p, and Hsa-miR-532-5p are concordant with the previous results.
In order to verify the differential expression of the core miRNAs in breast cancer patients experimentally, by applying 100 BC serum samples and 100 healthy controls, Real-Time PCR was performed that showed significant expression difference in breast cancer patients compare to the healthy controls (Figure 7).
Efficacy of diagnosis of hub miRNAs
Using ROC curve, we evaluated the diagnostic efficacy of six hub miRNAs based on Real-time PCR result. The AUC showed that miRNAs indicated excellent diagnostic efficiency for breast cancer in comparison with healthy samples (Figure 8). As it is clear from the figure P.value of all six hub miRNAs is meaningful (P.value < 0.05) and the AUC of six hub miRMAs including Hsa-miR-151a-5p (0.93), Hsa-miR-34a-5p (0.67), Hsa-miR-1307-3p (1.00), Hsa-miR-450b-5p (0.76), Hsa-miR-501-3p (0.84), and Hsa-miR-532-5p (0.88) show excellent diagnosis efficacy.
Validation of core miRNAs based on METABRIC database
To further validation of hub miRNAs, we evaluated the survival plot of them using METABRIC database that present an integrated analysis of copy number and gene expression in a discovery and validation set of 1262 primary breast tumors, with long-term clinical follow-up (Figure 9). As, it is clear from plots five of them including Hsa-miR-151a-5p, Hsa-miR-34a-5p, Hsa-miR-1307-3p, Hsa-miR-501-3p, and Hsa-miR-532-5p have meaningful P.Value, while Hsa-miR-450b-5p has P.Value = 0.086 which is meaningless.
MiRNA-mRNA network of hub miRNAs
We apply six hub miRNAs and constructed miRNA-mRNA network by CyTargetLinker, based on miRTarBase, Targetscan, and TransmiR databases (Figure 10). All hub miRNAs are upregulated which are presented in red while the target genes of hub miRNAs are downregulated and are presented in green.
Target genes expression evaluation by other datasets
The expression status of target genes was evaluated using 12 GEO series including GSE10797, GSE20711, GSE29431, GSE31448, GSE3165, GSE38959, GSE42568, GSE45827, GSE57297, GSE65216, GSE65194, and GSE61724 that in total cover the expression data of 1600 BC patients. The result was shown as heatmap (Figure 11). Target genes of hub miRNAs are downregulated and their expression is concordant with the expectations.
Enrichment of target genes
Target genes were enriched, and then analyzed through the cytoscape plug-in ClueGo. According to the KEGG database, these genes are predicted to play role in several cellular processes like p53 signaling pathway, insulin resistance, melanoma, central carbon metabolism in cancer, prostate cancer, longevity regulating pathway, amphetamine addiction, and cocaine addiction, and adipocytokine signaling pathway (Figure 12).