MiRNA target enrichment analysis of co-expression network modules reveals important miRNAs and their roles in breast cancer progression

2.1 Data source and pre-processing

RNA-seq, miRNA-seq, and clinical data of TCGA-BRCA project were downloaded by R package TCGAbiolinks [18], and imported to R programming environment. To eliminate possible effects of treatments on molecular mechanisms, those patients who have received any pharmaceutical or radiation therapy were removed from dataset. From those treatment-naive patients, tumor samples which have both RNA-seq and miRNA-seq data were chosen for further analysis. In addition, samples without stage information were removed. To filter out lowly expressed features, mRNAs and miRNAs which have less than 50 and 10 total counts number respectively among all samples were removed. Then sequencing data were normalized and log2 transformed. To find targets of miRNAs, Tarbase [19] and mirtarbase [20] which are both curated databases of mRNA-miRNA interactions were used. Both databases were merged into a single matrix and mRNAs that exist in the merged matrix were chosen and used to construct co-expression network. Before network construction, the quality of mRNA data was assessed by hierarchical clustering analysis.

2.2 Co-expression network construction

WGCNA [13], which is a powerful R package was used to construct signed gene co-expression network. To do so, first Pearson’s correlation coefficients between genes were calculated and used to construct adjacency matrix as follow:

where a is adjacency matrix, cor(x _i , x _j ) is Pearsion’s correlation coefficient between gene i and j, and β is soft threshold which eliminate weak correlation. To find appropriate soft threshold, expression matrix was transformed to similarity matrix and different values of β were tested by pickSoftThreshold function. Finally the smallest value which ensured constructed network has scale-free topology by R² > 0.9 was chosen. Next, topology overlap measure (TOM) was calculated according to adjacency matrix and different gene modules have been distinguished by using dynamicTreeCut method for dendrogram of dissimilarity (1-TOM). Minimum required size of detected modules were set to 30 (minModuleSize = 30) and modules with less than 0.25 height in dendrogram were merged together (MergeCutHeight = 0.25). Other parameters of blockwiseModules function of WGCNA package, were set as follow: genes were split to blocks with maximum number of 5,000 (maxBlockSize = 5,000); Pearson was used as correlation type (corType = “pearson”); deepSplit = 2 was used for moderate sensitivity over module splitting; minimum connectivity of 0.5 between module’s eigengenes and 10 members was set as threshold to maintain detected modules (minCoreKME = 0.5, minCoreKMESize = 10); also module’s member with connectivity less than 0.3 to corresponding module’s eigengene were omitted from the module (minKMEtoStay = 0.3).

For each module, a representative value (eigengene) was calculated. To find important modules in cancer progression, correlations between module eigengenes and stage-related traits were calculated. Those modules that have |correlation| bigger than 0.2 and p-value less than 0.05 with at least two traits related to progression (clinical stage and TNM) were considered as important modules.

2.3 MiRNA target enrichment and correlation analysis

To test whether targets of specific miRNA are enriched in specific network modules, hypergeometric test was performed as follows:

P_HT: hypergeometric test p-value

x: the number of genes in a module

n: the number of targets of a specific miRNA

z: the number of module’s members and miRNA’s targets

N: total number of genes in the co-expression network

The total number of genes were 17,463 (N = 17,463) and the test was performed for 770 miRNAs against significant modules. Because multiple tests were performed, p-value was corrected by Benjamini–Hochberg method and False Discovery Rate (FDR) < 0.05 was chosen as the first threshold. As the second criterion, we used correlation between miRNAs and modules’ eigengenes. Here we observed Pearson’s correlation coefficients ranging from −0.56 to 0.89 and because miRNAs mostly have negative regulatory effect on gene expression, correlation less than −0.3 was set as the second threshold. Thus miRNA was flagged as significant if it had both significant number of targets in the specific module according to hypergeometric test result (FDR < 0.5) and negative correlation (coefficient < −0.3) with that module.

2.4 MiRNAs’ targeted genes GSVA score

Gene set variation analysis (GSVA) [21] is a popular and open-source R package and estimates variation of gene set in each sample of the population. For each significant miRNA, their target genes in targeted module were considered as a gene set and their GSVA score was calculated and compared among different stages of breast cancer.

2.5 Pathway and gene ontology enrichment analysis

DAVID database [22] and Cytoscape plugin, ClueGO [23] were used to performed functional enrichment analysis. When using DAVID, we chose three types of GO terms (biological process [BP], cellular compartment [CC], and molecular function [MF]) and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways as resources. While biological process (BP) gene ontology was used to construct network by ClueGO.

2.6 PPI network construction and cluster analysis

In order to find molecular interaction and possible functional module in targeted gene list of each miRNA, protein-protein interaction network was constructed by CluePedia [24], a powerful plugin for Cytoscape [25]. STRING [26] was used as an interaction resource and the confidence cut-off was set to 0.7 (high confidence). Molecular complex detection (MCODE) algorithm [27] was used to find clusters in PPI networks. Node density cut-off = 0 and node score cut-off = 0.2 were set as MCODE parameters.

3.1 Data pre-processing

After selecting treatment-naive patients and removing samples with missing mRNA, miRNA, or stage data, 106 out of 1098 TCGA-BRCA samples were included. Minimum, average, and maximum age of included patients was 30, 64, and 90 years respectively. Chosen data comprised 21 stage I, 58 stage II, 23 stage III, and 4 stage IV samples. Initial mRNA count matrix had 56602 Ensembl gene ID and after removing lowly expressed genes, it has 17463 ID in common with the merged tarbase and mirtarbase interaction data which had information of 20480 Ensembl gene ID. There was no missing data in count matrix and no outlier was seen in the data. Hierarchical clustering analysis and stage and TNM heatmap can be seen in Figure 1A.

Figure 1:

Dendrogram and heatmap represent distance between mRNA data samples and their corresponding stage and TNM data respectively (A) plots show scale free topology fitting R² and mean connectivity of different powers. Power = 16 fulfills R² > 0.9 and was chosen for network construction. (B) Dendrogram shows modules’ eigengenes and their distances. (C) Heatmap demonstrates important modules and correlations between their eigengenes and stage and TNM. P-values are also mentioned in parentheses (D).

3.2 Co-expression network construction and module detection

To construct co-expression network, correlation coefficients were raised to β = 16. As shown in Figure 1.B, it fitted scale-free topology by R² = 0.915 and has mean connectivity = 8.11. Forty distinct modules were found in constructed network. Grey module which consists of unassigned genes was removed in further analysis. The smallest module (mediumpurple3) had 30 genes and the largest one (turquoise) had 1,539 genes. Figure 1c shows hierarchical clustering of modules’ eigengenes. Correlations between eigengenes and stage and TNM were calculated and among all modules, four modules had |correlation| > 0.2 with at least two of stage-related traits. Figure 1D shows stage_modules correlation heatmap and their corresponding p-values and top GO term and pathways for important modules can be seen in Table 1. Saddlebrown and tan modules have negative correlation with cancer progression and top results of their biological process gene ontology enrichment analysis were cilium movements and retinol metabolic process respectively. There were also two other important modules which show positive correlation with cancer progression. Turquoise module was one of them and mitochondrial processes such as mitochondrial translation and oxidative phosphorylation were among its top GO terms and pathways. Another one was grey60 module and it consisted of small group of genes that exist in cytoplasm and nucleoplasm.

3.3 miRNA targets enrichment analysis

After performing hypergeometric test and calculating false discovery rate (FDR), 29 miRNA passed FDR < 0.05 threshold. Moreover, correlations between normalized miRNA expression and modules’ eigengenes were calculated and correlation < −0.3 was used as the second level of threshold and three miRNA passed both threshold. Hsa-let-7b-5p and hsa-mir-23b-3p had 471 and 372 targets in the turquoise module respectively and hsa-mir-30a-5p had 28 targets in grey60 module. The turquoise module had 1,539 and the grey60 module had 92 genes and both modules were positively correlated with stage-related traits which shows suppression effect of mentioned miRNAs on the breast cancer progression. Important miRNAs’ target enrichment results and their correlation with important modules are depicted in Figure 2A.

Figure 2:

Heatmap shows important miRNAs and their target enrichment analysis result for important modules. Corresponding FDR values are shown and correlation between miRNA expression and modules’ eigengenes are mentioned in parentheses. (A) For each miRNA, its targeted genes were considered as a single gene list and GSVA score was calculated for them. Boxplots show GSVA scores of miRNAs’ targeted gene lists against stages. Coefficients and corresponding p-value regarding correlation between GSVA scores and stages have been mentioned above each plot (B–D).

3.4 Important miRNA targets GSVA score

For each important miRNA, their targets were considered as a gene set and gene set variation analysis (GSVA) score was calculated and compared among different stages. As shown in Figure 2B–D, targeted gene set of each important miRNA (hsa-mir-23a-3p, hsa-let-7b-5p, and hsa-mir30a-5p) had a higher GSVA score in higher stages.

3.5 Pathway and gene ontology enrichment analysis of miRNAs’ target gene sets

Because hsa-mir-23b-3p and hsa-let-7b-5p targets were significantly enriched in the same module (turquoise), their functional roles were investigated together and their targets were treated as a single gene list. As ClueGO result (Figure 3A) demonstrate, mir-23b and let-7b target related cellular activities as expected. Major cellular activities such as mRNA metabolism and processing, mitochondrial gene expression, and oxidative phosphorylation enhance in higher breast cancer stages and it can be due to reducing inhibitory effects of mir-23b and let7b. These miRNAs seem to have complementary roles in breast cancer. For instance, mRNA metabolism is mostly controlled by mir-23b and mRNA processing is affected by let-7b; or mitochondrial gene expression is regulated by let-7b and mir-23b regulates oxidative phosphorylation.

Figure 3:

ClueGO result shows functional enrichment analysis for mir-23b and let-7b targets. Nodes show biological process GO terms with p-value < 0.05 and their colors indicate percentage of each miRNA’s targets. Blue nodes show terms consist of more than 60 % let-7b targeted genes and red nodes show terms with more than 60 % genes targeted by mir-23b. (A) MCODE first two modules are shown; blue and red nodes indicate let-7b and mir-23b targeted genes respectively and gray nodes are genes which are targeted by both. (B, C) DAVID pathway enrichment analysis result for second MCODE modules are sorted by log10(FDR) and plotted as bar-chart (D).

On the other hand, because hsa-mir-30a-5p targeted separated module (grey60), its functional roles were assessed individually. As represented in DAVID analysis (Table 2), targeted genes of mir-30a were mostly in cytosol and autophagy was a frequent term that had been seen in its functional analysis.

3.6 PPI network was constructed, and important modules were found

To investigate targets of mir-23b and let-7b holistically, PPI network was constructed by Cytoscape. Because targeted genes were originally co-expressed, it is possible to see different functional protein modules which work together in the same activity or even different genes which are subunits of one complex protein. To find these modules, MCODE was performed and two important modules with the highest scores were chosen. The first module which its MCODE score was 21.818 (Figure 3B), consists of 21 mithochondrial ribosomal proteins, mitochondrial translation initiation factor (MTIF1) and mitochondrial inner membrane protein OXA1L that interact with mitochondrial ribosomes and require for insertion of proteins to mitochondrial inner membrane. Second module with score = 18.182 is represented in Figure 3C. As DAVID analysis results indicated (Figure 3D), its genes are mostly involved in spliceosome and there are several genes related to transcription in this module such as several RNA polymerase subunits.

To examine mir-23b and let-7b targets in more detail, we also constructed separated PPI networks for each of them and run MCODE algorithm separately to find their important network’s clusters. For let-7b, the first module comprised mitochondrial ribosomal proteins and the second module consisted of spliceosomal proteins as same as modules in union PPI network of mir-23b and let-7b targets. Modules of let-7b targets’ PPI networks can be seen in Figure 4A and B and DAVID analysis result of second module is shown in Figure 4C. As shown in Figure 4D and E, first module of mir-23b targets mostly consists of proteosomal genes and its DAVID analysis result shows its genes mostly participate in ubiquitin activity in cell cycle and ubiquitin-dependent catabolism processes. Second module of mir-23b targets and their enrichment analysis are depicted in Figure 4F and G, and most of its nodes are spliceosomal protein and it indicates mRNA splicing as important target of both mir-23b and let-7b.

Figure 4:

MCODE first and second modules of let-7b targets’ PPI networks are shown. (A, B) DAVID analysis result of second module is depicted as bar-chart (C). Proteins and their interaction, and functional analysis result of mir-23b targets’ first MCODE module (D, E) and the second module (F, G).

4.1 Hsa-mir-30a-5p

One of the important miRNAs we found was hsa-mir-30a-5p and pathway enrichment analysis of its target showed autophagy as an important term. Several studies have shown mir-30a is a tumor suppressor miRNA in breast cancer and revealed its suppression effect on progression and metastasis 28], [29], [30], [31. Also its higher level of expression has been shown to be associated with better overall survival of patients [32]. Moreover there are also several studies that show mir-30a has an inhibitory effect on autophagy in different types of cancers including breast cancer [33, 34]. Most of these studies focused on beclin-1 as autophagy-related target of mir-30a, but here our results showed that the inhibitory effect of mir-30a goes behind beclin-1 and there are several autophagy-related genes which are targeted by mir-30a in breast cancer. Autophagy is a crucial process in different types of cancer and it plays an important role in breast cancer tumorigenicity [35]. There are also several studies indicated that it can have both anti- and pro-metastatic effects on breast cancer cells. One of the possible explanations is autophagy can help unattached cells to survive before dormancy in new environments [36] and it can explain why the lower expression of mir-30a is associated with growth in unattached breast cancer cells [37].

4.2 Hsa-let-7b-5p

Hsa-let-7b-5p was another miRNA that we found to be important in breast cancer progression. Let-7 microRNA family have been shown to be dysregulated in several malignancy in both cancerous tissues and blood sample of patients [38]. This family have been indicated to regulate JAK-STAT3 pathway via STAT3 targeting, and MYC pathway by regulating expression of c-MYC. Also there are evidences which shows the regulatory role of let-7 family on breast cancer metastasis and stemness mostly by targeting oncogenes such as RAS, HMGA2 [39]. Regarding to let-7b, a member of let-7 family, there is also an in vivo evidence of their downregulation in breast cancer [40]. Dysregulation of let-7b have been also shown in vitro, and it’s role in the apoptosis process has been indicated by breast cancer cell line transfection [41]. Moreover, a study indicated the effect of let-7b on cell growth and metastasis through reprogramming aerobic glycolysis via modulation of Hexokinase 2 (HK2) [42]. This had been followed by another evidence which indicated the role of metabolic reprogramming via let-7 family to be associated with tumor progression and autophagy [43].

Here we have found the targets of let-7b had been enriched in one of the important module in the breast cancer co-expression network (turquoise module) which had positive correlation with progression traits. Pathway enrichment analysis showed mitochondrial gene expression is in this module and targeted by let-7b. Moreover PPI network analysis of let-7b targets in the turquoise module, indicated this process maybe targeted mostly by targeting mitochondrial ribosomal gene expression. Localization of precursor and mature form of let-7 family microRNA have been indicated previously by in situ hybridization [44]. our results suggested metabolic reprogramming by let-7b could be due to upstream modulation of mitochondrial genes’ expression indeed by controlling regulation of mitochondrial ribosomal proteins. Also mitochondrial ribosomal genes belonged to breast cancer progression related module (turquoise). This suggest the importance of mitochondrial ribosome in breast cancer cell growth and metastasis. There are studies that showed the importance of mitochondrial ribosomes in breast cancer and their subunits are differentially expressed in this disease [45, 46]. Here we found this dysregulation maybe caused by let-7b.

4.3 Hsa-mir-23b-3p

Also we found hsa-mir-23b-3p as another important miRNA and together with let-7b, they have significant number of targets in the same module. Mir-23b have been shown to by dysregulated in several malignancies. Mir-23b has been shown to have dual role in different cancer types and evidences have indicated both anti-tumor effects and tumor promoting effects for this microRNA [47]. Also in breast cancer, there are several studies show mir-23b dysregulation but there is also duality of its roles in breast cancer formation and progression. Some studies emphasized on oncogenic effects of mir-23b and their positive effect on breast cancer progression. For instance, knockout of mir-23b and mir-27b have been shown to reduce cell growth in breast cancer cell line (MCF7) but had no effect on migratory behavior of breast cancer cell line [48]. Another study indicated elevation in mir-23b expression in malignant breast cancer tissue comparing to normal tissue; but decreasing in higher stages comparing to lower stages [49]. Other studies addressed their tumor-suppressing roles especially on progression and metastasis. Using mir-23b transfection in breast cancer cell line, one study have shown regulatory role of mir-23b on several cytoskeletal genes’ expression which cause cytoskeleton remodeling and increase motility of cancer cells [50]. Also another evidence indicated mir-23b stabilized epithelial phenotype and reverse Epithelial Mesenchymal Transition (EMT) which is an important behavior in metastasis [51].

Here we found mir-23b may have tumor-suppressing role in breast cancer and inhibits progression and metastasis. We have found mir-23b targets have been enriched in turquoise module and they were mostly involved in oxidative phosphorylation process and NF-KB signaling pathway. Previous study showed suppression of mir-23b expression by c-myc can change metabolism program and enhances glutamine metabolism. Considering converting glucose to lactate in anaerobic metabolism in cancer cells (Warburg effect), glutamine metabolism is pivotal to maintaining mitochondrial metabolism in cancer cells [52]. This phenomenon (glutamine metabolism) have been shown to be an important behavior in breast cancer cells like other malignancy [53]. Here we have shown oxidative phosphorylation in mitochondria is correlated to breast cancer progression and mir-23b maybe the key regulatory element in this process. Also there are a line of evidences that have shown NF-KB pathway can change the expression of mir-23b by increasing its expression [54, 55]. Interestingly our result indicated the regulatory effect of mir-23b on NF-KB pathway. This have been also demonstrated in other studies [56, 57]. This would be suggesting feedback loop between the expression of mir-23b and NF-KB pathway. Further study is necessary to investigate relationship between mir-23b and NF-KB to unravel their regulatory effects on each other in breast cancer.

Also our results show that let-7b and mir-23b affect mRNA metabolism and processing such as mRNA splicing in breast cancer. Spliceosomes play important roles in breast cancer and their components have been shown to be related to breast cancer patients’ survival [58]. To our knowledge, no study have demonstrated complementary role of let-7b and mir-23b in breast cancer. Also we showed several new processes which are targeted by these two miRNAs in breast cancer.

In summary, we identified mir-30a, mir-23b, and let-7b as important miRNAs related to breast cancer progression by enrichment analysis of WGCNA modules. They regulate expression of not just few specific genes but two progression-related modules in breast cancer which may contain even several hundreds of genes. We found mir-23b and let-7b targeted module contains metabolic-related genes and mir-30a targets several autophagy-related genes. Our finding revealed new mechanisms of these miRNAs in breast cancer and their possible roles in breast cancer progression.