Exploration of prognosis-related microRNA and transcription factor co-regulatory networks across cancer types.

The study of cancer prognosis serves as an important part of cancer research. Large-scale cancer studies have identified numerous genes and microRNAs (miRNAs) associated with prognosis. These informative genes and miRNAs represent potential biomarkers to predict survival and to elucidate the molecular mechanism of tumour progression. MiRNAs and transcription factors (TFs) can work cooperatively as essential mediators of gene expression, and their dysregulation affects cancer prognosis. A panoramic view of cancer prognosis at the system level, considering the co-regulation roles of miRNA and TF, remains elusive. Here, we establish 12 prognosis-related miRNA-TF co-regulatory networks. The characteristics of prognostic target genes and their regulators in the network are depicted. Although the target genes and co-regulatory patterns exhibit cancer-specific properties, some miRNAs and TFs are highly conserved across cancers. We illustrate and interpret the roles of these conserved regulators by building a model associated with cancer hallmarks, functional enrichment analysis, network community detection, and exhaustive literature research. The elaborated system-level prognostic miRNA-TF co-regulation landscape, including the highlighted roles of conserved regulators, provides a novel and powerful insights into further biological and medical discoveries.

Introduction

Prognostic information is important for clinicians treating patients with cancer; it may inform decisions about reasonable medical interventions and strategies for precision medicine [1,2]. Owing to recent advances in next-generation sequencing technology and its emerging application in various clinical settings, a number of signatures associated with survival outcomes have been extensively investigated. These markers, either genetic or epigenetic, carry various indicative features and clues for further biological and clinical discoveries [3,4].

The regulation of gene expression controls developmental, physiological, and pathophysiological processes in eukaryotic organisms. Associated dysfunction is tightly related to tumorigenesis and progression [5-7]. In the fine-tuned modulation at multiple levels, transcription factors (TFs) and microRNAs (miRNAs) have been recognized to play important roles at the transcriptional and post-transcriptional levels, respectively. The transcriptional program determines cancer phenotype and prognosis by shaping the gene signature in cancer cells [8]. Detectable dysregulated miRNAs in tumour biopsies have readily emerged as promising diagnostic, prognostic and therapeutic indicators [9-11]. In particular, increasing evidence suggests the existence of cooperation and crosstalk between miRNAs and TFs, mainly to buffer gene expression and/or adjust signalling [12]. Specifically, miRNAs and TFs can coordinatively regulate shared target genes in feed-forward loops (FFLs) [13]. Indeed, as recent studies have shown, perturbations of the interwoven regulatory system involving miRNAs and TFs may trigger global alterations in gene expression and affect cancer prognosis (Fig. S1). For example, in colorectal cancer, Mullany et al. found the expression of TFs and their related miRNAs together influence survival [14] and Wang et al. pointed out abnormal expression of two miRNAs (hsa-mir-25 and hsa-mir-31), one TF (BRCA1), and two other genes (ADAMTSL3 and AXIN1) affected patient survival [15]. Fulciniti et al. exhibited the existence of a novel miRNA-TF FFL with a critical role in growth and survival in multiple myeloma [16]. Kong et al. identified an interwoven network of miRNAs and TFs that regulates CD147, a known risk factor for breast cancer associated with poor prognosis in breast cancer patients [17].

The biological network is an integrated and system-level lens through which researchers may uncover the mechanism underlying disease [18]. At the network level, miRNA-TF FFLs are major network motifs (i.e., interconnection patterns that occur more often by chance in biological networks), forming the basic building blocks of the miRNA-TF co-regulatory network [19-21]. Despite substantial efforts to identify the prognostic signatures and potential roles of FFLs in prognosis, an integrative and system-level analysis remains lacking. Hence, we seek to investigate prognostic signatures and the regulatory mechanism behind them in the context of the miRNA-TF co-regulatory network.

In this study, we discerned prognostic FFLs and constructed 12 prognosis-related miRNA-TF co-regulatory networks by integrating cancer genomics data, prognostic signature findings, and interactome data. The characteristics and features of the co-regulatory network are summarized. Through a detailed analysis, we found some miRNAs and TFs are common and conserved during co-regulation of the prognostic network. After considering the conserved regulators, we designed a hierarchical model associated with cancer hallmarks to elucidate the regulatory mechanism affecting cancer prognosis. Detailed enrichment analysis revealed a common theme among prognostic signatures in the co-regulatory network. Several potential prognostic modules were identified inside the miRNA-TF co-regulatory network. A comprehensive survey of the conserved regulators, including literature consultation, was performed to validate and highlight their pan-cancer prognostic functionality. The investigation of prognostic miRNA-TF co-regulatory networks provides novel insight into cancer outcomes, elucidates the commonality among regulatory mechanisms, and offers implications for clinical biomarkers and the study of therapeutics.

Results

Construction of prognostic miRNA-TF co-regulatory networks for human cancers

We developed a five-step pipeline based on the framework in our previous studies [22] in order to build prognosis-related miRNA-TF co-regulatory networks across cancers types.

First, we obtained prognostic genes and miRNAs from the Human Pathology Atlas (HPA) and OncomiR, respectively (Fig. 1A). We then used the information provided by 10 established regulatory databases (Fig. 1B). In total, 72,801 TF-gene, 178,689 TF-miRNA, 305,858 miRNA-gene interactions were collected. We identified prognostic regulatory interactions whose target nodes or regulator nodes are known to be relevant to prognosis and formed a combinatorial network by merging all interaction types (Fig. 1C). Using the network motif detection algorithm, we identified three types of FFLs (TF-FFLs, miRNA-FFLs, and composite-FFLs) in the combinatorial network (Fig. 1D). We then constructed the co-regulatory network, which comprises three types of FFLs (Dataset S1) and incorporated expression data from The Cancer Genome Atlas (TCGA) to discern more precise FFL patterns in each network (Fig. 1E) (Dataset S2). The final information for each prognosis-related network is shown in Table 1 and Figure S2.

Figure 1.

An overview of the computational approach to build prognostic miRNA-TF co-regulatory networks in human cancers. (a) We collected prognosis-related genes and miRNAs for 12 cancers by referring to established databases. (b) Regulatory relationships were obtained from 10 public interactome resources. (c) We screened out prognosis-related regulatory interactions whose target nodes or regulator nodes are known to be relevant to prognosis, forming an entirely synthetic network by merging all interaction types. (d) We then identified three types of FFLs from the combinatorial network using a network motif detection algorithm. (e) We constructed the co-regulatory network which comprises three types of FFLs and incorporated expression data from TCGA to filter out more precise FFL patterns in each network.

Table 1.

Summary of FFLs in 12 prognosis-related co-regulatory networks.

Cancer	TF-FFLs	miRNA-FFLs	Composite-FFLs	TF-Gene	TF-miRNA	miRNA-Gene	nTF	nmiRNA	nGene
BRAC	77	27	15	93	63	106	32	27	59
CMM	29	14	5	44	28	49	22	16	31
CXSCC	101	27	20	129	59	140	39	31	93
ENAC	542	121	96	571	285	672	93	116	394
HCC	1498	959	188	2048	589	1996	183	194	1107
HNSCC	269	51	36	271	170	294	60	63	179
OVAC	23	3	19	41	24	54	9	18	39
PAAC	487	107	61	469	346	583	101	117	307
PRAC	87	22	43	102	57	150	20	29	85
STAC	206	25	26	176	144	207	36	47	107
THYCA	149	35	25	129	131	183	52	63	77
UC	727	127	119	677	322	867	69	124	518

The landscape of prognostic miRNA-TF co-regulatory networks

To assess the topological structure of all 12 co-regulatory networks, we examined the degree distribution of each network (Table S1). The results showed that each prognosis-related co-regulatory network followed the power-law distribution, indicating that these co-regulatory networks had scale-free characteristics, a common feature of most types of biological networks [23].

We investigated the components of co-regulatory networks by examining similarities in FFLs, genes, TFs, and miRNAs across cancer types. In order to measure the pairwise overlap of FFLs, genes, TFs and miRNAs across cancer types, we used two metrics: Fisher’s p-value and the Jaccard index (Fig. 2 and Fig. S3). Little similarity of prognostic FFLs with genes was observed across 12 cancers, suggesting a relatively limited number of common FFLs and genes. In contrast, significant overlap was observed among regulators (i.e., miRNAs and TFs) in the networks. The low conservation of target genes in the prognostic co-regulatory network was concordant with previous studies finding that prognostic genes themselves lack cross-cancer conservation [24,25], which also led to the lack of cross-cancer conservation of FFLs that comprise prognostic genes. However, when we focused on miRNAs and TFs that regulate target genes, some miRNAs and TFs played roles in multiple prognostic co-regulatory networks. For example, ETS1 had a regulatory role in 12 co-regulatory networks; ESR1, MYC, and GATA2 appearing in 11 networks, respectively. These results indicate that conserved regulators impact non-conserved prognostic targets. Based on this FFL pattern, some conserved miRNAs and TFs may influence the cancerous and clinical outcomes for multiple cancers.

Figure 2.

Heat map showing the Fisher’s -log (p-value) for the pairwise overlap of FFLs, genes, miRNAs and TFs between the prognostic co-regulatory networks.

Regulators common to multiple prognostic co-regulatory networks

To systematically investigate regulators that acted across networks, miRNAs and TFs were grouped into two categories based on the extent to which a regulator was common across prognostic-related co-regulatory networks. We analysed the distribution of target genes, miRNAs, and TFs in different cancers (Fig. 3A–C). Only 0.19% of target genes were present in >6 cancers; 3.37% miRNAs and 6.09% TFs occurred in ≥7 cancers. We defined ‘common’ regulators (including common miRNAs and common TFs) as miRNAs or TFs that occurred in ≥7 cancer networks (18 TFs and 11 miRNAs). All others were considered to be ‘specific’ regulators. The 18 common TFs and 11 common miRNAs were shown in Fig.3D–E, and the lists of prognostic target genes and ‘specific regulators’ were recorded in Supplementary datasets S3-S5.

Figure 3.

The conserved regulators in prognostic miRNA-TF co-regulatory networks. (a-c) Occurrence of prognostic (a) genes, (b) miRNAs and (c) TFs in 12 co-regulatory networks. (d-e) The identification of common (d) miRNAs and (e) TFs occurred in ≥7 co-regulatory networks.

Using this classification, we divided the regulatory elements in each prognostic co-regulator network into two categories. Nodes with high degree (regarded as hub nodes) are known to play important roles in networks. We compared the degree of common vs. specific regulators in each network (Fig. 4A). In nine cancers, common TFs were hub nodes with significant high degree, rather than non-common TFs. Common miRNAs had a significant high degree in five cancers. These results highlight the pivotal function of common regulators conserved across prognostic co-regulatory networks. Compared to common TFs, such tendency of common miRNAs was weaker, which may be explained by the relatively weaker regulatory function of miRNAs.

Figure 4.

Common regulators may govern and maintain prognostic co-regulatory networks’ architecture across cancers. (a) Common regulators tend to have a higher degree in each network. The left plot shows the node degree comparison of common TFs vs. specific TFs, and the right plot makes a comparison of common miRNAs. (b) Common regulators controlled broad FFLs (left panel) and target genes (right panel) in each network.

For each co-regulatory network, we further investigated the number of FFLs and target genes containing common regulators (Fig. 4B). It is noteworthy that common regulators controlled broad FFLs (p-value <0.01, Wilcoxon test, paired) and target genes (p-value <0.01, Wilcoxon test, paired) in each network. This suggests that common regulators not only affect multiple cancers but also govern broad targets through FFLs in each co-regulatory network. Both of these findings indicate that common regulators may govern and maintain prognostic co-regulatory networks’ architecture across cancers.

A hierarchical model to illustrate the roles of common regulators in cancer prognosis

Despite the daunting complexity and remarkable diversity of neoplastic diseases, several cancer hallmarks contribute to the development of human tumours [26,27]. Each cancer hallmark represents the biological capability for oncogenic progress that underlies [28] the tumour phenotype. Cancer hallmark genes also have a tight and subtle relationship with cancer prognosis [29]. After considering the important role of common regulators in the co-regulatory network, we designed a hierarchical model related to prognosis in order to elucidate the contribution of common regulators to the development of cancer and the effect of prognostic genes on cancer biology (Fig. 5A).

Figure 5.

A hierarchical model associated with cancer hallmarks. (a) A hierarchical model considering cancer hallmarks to comprehend the functions of common regulators in cancer prognosis. (b) Heatmaps containing the number of target genes (under the control of common regulators), common TFs and common miRNAs that link to cancer hallmarks across cancer types.

This model lays out a hierarchy for 12 major cancers, common regulators, the cooperative regulators of common regulators, target prognostic genes, annotated GO terms related to biological progress, and nine cancer hallmarks. Common regulators (top layer), govern prognostic target genes directly or with the help of their cooperative miRNAs/TFs, in an FFL pattern. These prognostic genes are enriched in specific GO (p-value <0.05, Dataset S6) and GO terms are correlated with various cancer hallmarks, demonstrating the vivid interaction between cancer hallmarks and clinical outcomes.

We next enumerated the common regulators and target genes linked to each GO term related to cancer hallmarks, across diverse cancer types (Fig. 5B). We found that common regulators and target prognostic genes tended to appear in cancer hallmarks for ‘sustaining proliferative signaling’ and for ‘tissue invasion and metastasis’. Sustaining proliferative signalling plays a fundamental role in cancer. Cancer cells display autonomous, chaotic growth because of dysregulated growth signals. Tissue invasion and metastasis are known to be associated with the progression of carcinoma [30,31]. Our model suggests that the production and release of growth-promoting signals and the invasion-metastasis cascade may play key driver roles, reflecting the organizing principle of common regulators with respect to cancer prognosis.

We sought to further elucidate the biological role of common regulators and target genes in the model. Pathway enrichment analysis was performed for common TFs, common miRNAs, and prognostic genes under the control of common regulators (p-value <0.05, Fig S4 and S5). In addition to pathways related to specific cancers, the cell cycle is the most shared enriched pathway for common regulators and their targets. Disturbance of the cell cycle has already been proven significant in the prognosis of several cancers [32,33]. Furthermore, pathways that target genes enriched, such as focal adhesion, TGF-ß and p53 may be significant for prognosis [34-37]. While pathways such as ERBB signalling and regulation of the actin cytoskeleton are previously less characterized pathways in prognosis. The results of pathway enrichment analysis may be used to identify additional regulators and genes related to cancer prognosis.

miRNA-TF cooperative modules as prognostic biomarkers in multiple cancer types

Networks present modular structure, and decomposition of the network is beneficial for the elucidation of complex systems [38]. Compared with individual genes, module biomarkers are more powerful predictors of prognosis [39]. Given the characteristics of pan-cancer and core roles inside each network, we are interested in finding modules that comprise common regulators associated with the survival of cancer patients.

We used the GLay community detection algorithms to decompose each prognostic co-regulatory network. In total, 115 modules were identified (Dataset S7). The results of survival analysis for each module identified 32 modules that could predict overall survival (log-rank test, p < 0.05); 23 out of 32 of these modules comprised common regulators (Dataset S8).

We focused on two significant modules of cervical cancer (CXSCC), a common gynaecological cancer. The two modules comprised the common regulators MYC and GATA2 (Fig S6), and the elements contained in the two modules are closely related to cancer progression. Amplification and overexpression of MYC are related to CXSCC progression and GATA2 mutations cause a multifaceted disorder [40,41].In the first module, GATA2 regulates both hsa-mir-30e and the target genes IL1A and ITGA5, with the hsa-mir-30e repressing the target genes. The target IL1A can promote tumour growth, invasion and migration [42] and the ITGA5 expression is induced in transformed epithelial cells during epithelial to mesenchymal transition (EMT) process which fuels metastasis by endowing cells with enhanced migratory and invasive potential [43-45]. In the second module, MYC and hsa-mir-342 control a joint target TFRC. Hsa-mir-342 has the potential to suppress cell proliferation, migration and invasion of human cervical cells [46]. Clinical data have shown that high TFRC expression in cervical cancers is related to advanced clinical pathologic characteristics, and the TFRC is also an independent predictor for survival in cervical cancer [47].

The above results showed that larger clusters of FFLs may play a role in prognosis stratification. Furthermore, common regulators may affect prognosis in larger modules.

A comprehensive survey of common regulators

After establishing the significant roles of common regulators, we investigated the association between common regulators and pan-cancer signatures (Fig.6A, B). We collected four consensus lists of pan-cancer gene data [48-51]. Notably, the overlap result showed that 14 out of 18 common TFs have been identified as pan-cancer genes. The most notable TF is EZH2, which occurred in all four data sets. Although current studies do not provide enough data on pan-cancer miRNAs, the overlap between common miRNA and two pan-cancer-related miRNA datasets [52,53], namely ‘Pan-cancer miRNA superfamily’ and ‘SDEmiRNA’, showed that common miRNA hsa-mir-93 is oncogenic in multiple cancers. These results suggest the dual function of several TFs and miRNAs in oncogenesis and prognosis.

Figure 6.

Overlap between common regulators and pan-cancer signatures. (a) A Venn diagram showing the overlap between common miRNAs and two pan-cancer miRNAs datasets, namely pan-cancer miRNA superfamily and SDEmiRNA. (b) A Venn diagram showing the overlap between common TFs and four pan-cancer gene datasets, namely Cosmic, Cancer5000, Netsig5000 and IntOGen.

We then carried out a detailed literature survey of common regulators. We searched the PubMed database with keywords including ‘prognosis’, ‘prognostic’, ‘survival’, and ‘clinical outcome’ for each common regulator that we found. We manually extracted analyses related to cancer. As a result, we consulted about 130 published studies describing the associations between common regulators and cancer prognosis (Table 2, 3). Common TFs are reportedly associated with the prognosis of 3–7 cancers. Remarkably, the most common TFs reported are ETS1 and EZH2; the former is a common TF that plays roles in 12 co-regulatory networks, while the latter is the significant one who has a dual function in both oncogenesis and prognosis as noted before. The most heavily studied prognostic common miRNA is hsa-mir-34a, which is related to 18 cancers. Common miRNAs such as hsa-mir-9–2, hsa-mir-23b, and hsa-mir-361 have not previously been investigated, and further study will be necessary to verify their pan-cancer prognostic potential. These findings validate and support the pan-cancer prognostic functionality of conserved regulators.

Table 2.

Published studies describing the associations between common TFs and cancer prognosis.

Common TFs	Reported cancer	PMID
AR	Breast cancer	26526356
	Osteosarcoma	28262798
	Colorectal cancer	25376484
	Prostate cancer	30105831
BRCA1	Breast cancer	24258259
	Ovarian cancer	25398451
	Non-small-cell lung cancer	27179511
CEBPA	Prostate cancer	30430607
	Cervical squamous cell carcinoma	24913332
	Acute myeloid leukemia	15746035
	Hepatocellular carcinoma	25363290
CREB1	Gastric cancer	29870889
	Breast cancer	17786359
	Prostate cancer	26743006
	Colorectal cancer	27046651
	Ovarian cancer	22596241
E2F1	Cervical cancer	28559983
	Lung carcinoma	29754146
	Ovarian cancer	28667302
	Gastric cancer	28569791
E2F4	Breast cancer	25440089
	Bladder cancer	26032289
	Lung cancer	29754146
ESR1	Thyroid carcinoma	28124274
	Breast cancer	29482551
	Ovarian cancer	24368280
ETS1	Breast cancer	26392377
	Cervical cancer
	Colorectal cancer
	Gastric cancer
	Lung cancer
	Oral cancer
	Ovarian cancer
EZH2	Cervical cancer	27697098
	Head-and-neck squamous cell carcinoma	26604082
	Renal clear cell carcinoma	30405850
	Non-small cell lung carcinoma	24097870
	Hepatocellular carcinoma	27920552
	Colorectal cancer	29061982
	Oral squamous cell carcinomas	18619895
GATA1	Clear cell renal cell carcinoma	25230694
	Acute erythroid leukemia	27086927
	Breast cancer	22020876
GATA2	Clear cell renal cell carcinoma	25230694
	Acute erythroid leukemia	19097174
	Colorectal cancer	26287967
	Hepatocellular carcinoma	24498120
HIF1A	Hepatocellular carcinoma	26115041
	Non-small cell lung cancer	24631267
	Oral cancer	19449077
	Pancreatic cancer	18362831
MYC	Breast cancer	24316975
	Gastric cancer	25618371
	Acute myeloid leukaemia	26856970
	Colorectal cancer	24503701
	Lung adenocarcinoma	21148746
RELA	Pancreatic cancer	17622249
	Chronic lymphocytic leukemia	19124804
	Non-small cell lung cancer	18215193
SP1	Hepatocellular carcinoma	28028181
	Glioma	21469139
	Colorectal cancer	22821729
	Gastric cancer	15217947
STAT3	Gastric cancer	27938379
	Diffuse large B-cell lymphoma	21806788
	Cervical cancer	19638983
	Ovarian cancer	17063503
TFAP2A	Bladder cancer	21489314
	Nasopharyngeal carcinoma	24335623
	Breast cancer	21375726
	Gastric adenocarcinoma	21966377
TP53	Glioblastoma	24248532
	Breast cancer	26910472
	Colorectal cancer	22038927
	Thymic carcinoma	25299233
	Head and neck squamous cell carcinoma	25108461
	Pancreatic cancer	25428385
	Hepatocellular carcinoma	21616106

Table 3.

Published studies describing the associations between common miRNAs and cancer prognosis.

Common miRNAs	Reported cancer	PMID
hsa-mir-9–2	Hepatocellular carcinoma	23364900
		26046780
hsa-mir-23b	Ovarian cancer	24997860
	Colorectal cancer	26269151
hsa-mir-30d	Ovarian cancer	30095616
	Prostate cancer	28241827
	Hepatocellular carcinoma	26046780
hsa-mir-34a	Acute myeloid leukemia	29945348
	Cholangiocarcinoma	30050323
	Chronic lymphocytic leukaemia	30111844
	Hepatocellular carcinoma	29303511
	Colorectal adenocarcinoma	28624481
	Cervical cancer	28615991
	Ewing sarcoma	25015333
	Bladder cancer	25556547
	Laryngeal squamous cell carcinoma	27450916
	Sinonasal squamous cell carcinoma	22624980
	Glioma	23529798
	Gastric lymphomas	24232982
	Non-small-cell lung cancer	19736307
	Colon cancer	23243217
	Ovarian cancer	21516127
	Renal cell carcinoma	29104726
	Prostate cancer	25053345
	Breast cancer	22439831
hsa-mir-93	Non-small cell lung cancer	29309884
	Gastric cancer	28842285
	Breast cancer	28518139
	Colon cancer	23354160
	Lung cancer	24037530
	Cervical cancer	30098344
hsa-mir-150	Hepatocellular carcinoma	28811864
	Pancreatic cancer	25906450
	Prostate cancer	25778313
	Esophageal squamous cell carcinoma	23013135
	Colorectal cancer	22052060
hsa-mir-155	Clear cell renal cell carcinoma	30278113
	Cervical cancer	27470551
	Lung cancer	16530703
	Prostate adenocarcinoma	25938433
	Acute myeloid leukemia	25428263
	Head and neck squamous cell carcinoma	28347920
	Oral squamous cell carcinoma	27035278
	Bladder cancer	27035278
	hepatocellular carcinoma	27035278
	Glioblastoma	23302469
	Colorectal cancer	29361687
hsa-mir-221	Multiple myeloma	28168095
	hepatocellular carcinoma	22009537
	Breast cancer	30110679
	Ovarian cancer	28350128
	Bladder cancer	29181884
	Colon cancer	25932237
	Glioma	25636684
	Thyroid cancer	28061868
	Gastric carcinoma	27712596
	Prostate carcinoma	19676045
hsa-mir-335	Glioma	22644918
	Gallbladder carcinoma	24250228
	Gastric cancer	29075357
	Breast cancer	24132943
hsa-mir-361	Non-small cell lung cancer.	28051257
	Breast cancer	27959953
hsa-let-7g	Head and neck squamous cell carcinoma	30171046
	Non-small cell lung cancer	23820752
	Oral cavity squamous cell carcinoma	25050621

Discussion

In the present study, integrated data and network-based methods were used to identify miRNA-TF cooperative events for cancer prognosis. Twelve prognosis-related co-regulatory networks were identified by our multi-step pipeline. Since the incorporation of multi-omics data, prognostic signatures, mechanistic regulatory information, and careful refinement in the pipeline of network construction, the prognostic miRNA-TF co-regulatory network is powerful and reliable. MiRNAs and TFs may jointly regulate gene expression in the form of FFLs, which impact many aspects of cellular processes and disease progression. The miRNA-TF co-regulatory network brings a system-level heuristic view of gene expression regulation to cancer prognosis. A panoramic view of the functional networks may help to characterize prognostic targets and conserved regulators. We observed that common regulators maintain the structure of the co-regulatory network. This motivates us to move the focus from heterogeneous prognostic genes to their regulators – especially conservative regulators. MiRNAs and TFs have been treated as diagnostic, prognostic, and therapeutic objects. Evidence shows that therapies targeting TFs constitute an important part of the most commercially successful drugs approved by the US Food & Drug Administration. A broader effect was found when therapies targeting miRNAs and TFs were compared with those targeting a single gene [54].

Although the prognosis of cancer is as complex as cancer itself, the results presented above pertaining to common regulators demonstrate the cascaded regulatory principle among regulators, prognostic targets, and cancer hallmarks. We note that ‘sustaining proliferative signaling’ and ‘tissue invasion and metastasis’ are vital to cancer prognosis as cancer hallmarks. Pathways enrichment analysis and network community detection were used to elucidate the biological and topological roles of conserved regulators. Finally, we conducted a comprehensive survey of common regulators, emphasizing evaluation of the pan-cancer prognostic function of conserved regulators.

Conclusion

In this study, we investigated 12 miRNA-TF co-regulatory networks in the context of cancer prognosis, in order to elucidate prognostic signatures and the regulatory mechanism behind them. This network-based study highlights conservative regulators (beyond the prognostic genes that vary across cancers), highlighting the clinical importance of regulatory mechanisms in prognosis. We hope our work opens new avenues for the study of cancer prognosis and accelerates the development of precision medicine.

Materials and methods

Collection of prognostic miRNAs and genes

In this study, we focused our analysis on 12 tumour types. We used open resources to obtain the genes and miRNAs for 12 major clinical cancer outcomes. Prognostic genes were identified using the Human Pathology Atlas [29]. MiRNAs associated with clinical outcomes were identified using OncomiR [55] (significance threshold:0.05) (Table S2).

Regulatory relationships among miRNAs, TFs, and target genes

Regulatory relationships among miRNAs, TFs, and target genes were determined from public databases, as follows: (i) TF-gene: ITFP, TRED, TRRUST, HTRIdb [56-59]; (ii) miRNA-gene: miRTarBase, miR2Disease, miRecords [60-62]; (iii) TF-miRNA: mirTrans, PuTmiR, TransmiR [63-65] (Table S3). In this paper, the term ‘TF’ refers specifically to TF genes; the term ‘gene’ includes both TF and non-TF genes; ‘target gene’ refers to only non-TF genes. We unified TF/miRNA/target gene symbols in the regulatory relationship by referring the Hugo Gene Nomenclature Committee (HGNC) [66], the approval TF list [67,68], and miRBase [69].

Omics data across multiple cancers

We capitalized on expression data from TCGA [70] to filter more precise co-regulatory interactions. The TCGA clinical data were used for survival analysis of network clusters (Table S4). The R package RTCGAToolbox [71] was used to assess TCGA data (run date Jan. 2016) provided in the Firehose data repository and to perform survival analysis. The RNA-Seq expression values were transformed by the log2(x + 1) transformation, where ‘x’ is the original expression value. These values were used for subsequent analyses.

Network motif detection

Based on the collected prognostic signatures and regulatory relationships, we filtered out prognostic regulatory interactions whose target nodes or regulator nodes are known to be relevant to prognosis. We then pooled the prognostic regulatory relationships including TF-gene, TF-miRNA and miRNA-gene, generating a combinational network. Using a high-efficiency FANMOD algorithm for network motif detection [72], we detected three types of three-node FFLs: TF-FFLs, miRNA-FFLs, and composite-FFLs and formed the raw co-regulatory network.

Prognostic co-regulatory network construction and refinement

The preliminary co-regulatory network comprises three types of FFLs. A single TF-FFL has a master TF that regulates a partner miRNA and their joint target. A miRNA-FFL contains a master miRNA regulator, which represses its partner TF and their joint target gene. In a composite-FFL, the miRNA and TF regulate each other, thereby controlling their joint target.

We used TCGA expression data to discern more precise co-regulatory interactions in the raw co-regulatory network. We calculated pairwise Spearman correlation values among TFs, miRNAs, and genes for each FFL in the raw co-regulatory network. For TF-genes and TF-miRNA pairs, we retained p < 0.05 as the level of statistical significance. For miRNA-gene, we retained p < 0.05 and correlation coefficient <0, because most miRNAs are assumed to inhibit the expression of their targets. We removed less-significant FFLs from the raw co-regulatory networks in order to yield the final prognosis-related co-regulatory networks.

Classification of network regulators

For systematic analysis of the regulators in co-regulatory networks, we split the intra-network miRNAs and TFs into two groups: (i) common regulators: miRNAs or TFs that occurred in ≥7 co-regulatory networks; (ii) specific regulators: regulators with frequency <7 across 12 networks.

Network visualization, topological measurements, and identification of network modules

The miRNA-TF co-regulatory networks and the cancer hallmark-associated model were visualized with Cytoscape 3.7.0 [73]. Topological measurements of the networks were obtained using the NetworkAnalyzer plugin for Cytoscape. We used the GLay community clustering algorithm to generate clusters from large complex networks [74].

Functional annotation and enrichment analysis

The R/Biocondutor software ClusterProfiler [75] package was used for enrichment analysis. Data from MSigDB(v6.2) [76], miEAA (release date Apr.2016) [77] were utilized in enrichment analysis. Specifically, we selected MSigDB GO gene sets for GO enrichment analysis, and chose MSigDB KEGG gene sets for pathway enrichment analysis. We utilized miEAA pathway annotation to perform miRNA pathway enrichment analysis. To build a linkage between target genes and cancer hallmarks, we referred to a previous study [78] to determine a list of GO terms related to the hallmarks (Table S5).