Landscape of MicroRNA Regulatory Network Architecture and Functional Rerouting in Cancer

Genetic changes are a primary source of oncogenesis. Somatic mutations accumulate in cells during one's lifetime. Among them, some may contribute to a cell's selective advantage by altering activity in pathways associated with tumor formation and therapeutic resistance, and are therefore often referred to as “driver mutations.” Other variants may not have either phenotypic consequences or biological effects that are selectively advantageous to cancer cells, and are thus defined as “passenger mutations.” Many studies have paid attention to missense mutations in protein coding regions that can drive cancer development. A number of driver mutations have been identified in previous research efforts, such as BRAF (V600E), IDH1 (R132H), PIK3CA (E545K), EGFR (L858R), and KRAS (G12D; ref. 1). However, in addition to missense mutations that affect the protein coding components of human genome, miRNAs and their target binding sites occupy a significant proportion of the genome and can harbor somatic mutations, which play driver roles through miRNA-related pathways.

miRNAs are endogenous regulatory noncoding RNAs that are ∼22nt in length and act by targeting messenger RNAs (mRNA) for cleavage or translational repression. The diversity and abundance of miRNA targets contribute to the complexity of gene regulatory networks. Increasing lines of evidence have demonstrated that miRNAs play critical functions in various developmental, physiological, and pathological processes including cancer. Deregulation in the expression of miRNAs and their targets have been observed in various types of human cancer. In medulloblastoma, miR-217 was found to promote an oncogenic phenotype by blocking target genes such as histone deacetylase SIRT1 and SMAD family member SMAD7, inhibiting apoptosis (2). Furthermore, IGF2BP1 was found to interfere with miRNA-mediated mRNA suppression, enabling increased expression of mRNAs that contribute to tumorigenesis in ovarian cancer (3).

Mutations in miRNAs or their target genes may exert important effects on their deregulated expression. A number of studies have suggested that somatic mutations could impact miRNA–gene interactions (mGI) and are related to cancer development (4). Nucleotide variants in mRNAs within or near miRNA target regions may alter miRNA–mRNA binding by modifying mRNA structure nearby. For instance, a mutation within the E2F1 mRNA's miRNA target site was found to deregulate miRNA binding, contributing to oncogenesis in colorectal cancer (5). Beyond cancer, a gain-of-function mutation in miRNA-140 drastically altered its mRNA target profile, contributing to pathogenesis in human skeletal dysplasia (6). Taken together, mutations may potentially deregulate miRNA–mRNA interactions by several different mechanisms.

Several initial studies explored the mutational effect on miRNAs and their targets. SomamiR (7) mapped mutations to miRNAs and their targets, and developed a tool to calculate enrichment of mutated targets on KEGG pathways. PolymiRTS (8) used a TargetScan score to evaluate the effects changed by mutations on miRNA targets. However, mutations in specific cancer types and in individual patients were not considered in these studies. Many miRNA family members can form opposing associations with mRNAs, with some miRNAs exhibiting a positive correlation in expression levels whereas other members in the same family exhibit a negative correlation (9). Moreover, mRNAs found to bind many miRNA partners can change partners strongly based on the cellular environment (10). Such findings point to a need to assess miRNA–mRNA interactions independently in cancer-specific and patient-specific contexts.

miRNA targeting occurs not only in 3′ untranslated regions (3′ UTR), but also in coding regions and 5′ UTRs (11–16). Silent and nonsilent mutations in these regions have been shown to be relevant for cancer and other diseases (5, 11, 17). However, the mechanism of how mutations influence miRNA target genes was not clear in such studies. Several possible mechanisms have been elucidated, including the creation of novel splice variants, dosage control of oncogenes, and modifications in mRNA secondary structure (18–20). A study by Stegeman and colleagues used miRNA mimics and reporter gene assays and showed that mutations could change the expression of target genes in prostate cancer (21). Previously, miRNA-driven deregulation of cell-cycle genes and oncogenes has been documented to drive instances of human cytomegalovirus and pancreatic cancer (22, 23). In addition, compound mutations affecting both miRNA and mRNAs involved in an interaction are associated with a deregulated disease phenotype (24, 25). But it remains elusive the extent to which mutations could impair mGIs, and a global view of driver mutation-mediated gene regulatory network perturbations on a pan-cancer scale is needed.

Recently, with the development of high throughput sequencing projects, such as The Cancer Genome Atlas (TCGA) and International Cancer Genome Consortium (ICGC), numerous somatic mutations have been identified in various cancer types. Large-scale computational analysis has proven to be useful in characterizing mutations of unclear functional significance in cancer, suggesting that such a study with a focus on miRNA–mRNA interactions could lend novel insights (26–28). These mutations might create new miRNA binding sites or lose several binding sites, which further perturb miRNA–gene regulation. The quantity of binding sites is crucial to mGIs, as multiple binding sites for a single miRNA show stronger regulatory significance regardless of the sites’ location in the mRNA (29). To address these goals, we interrogated to what extent and how somatic mutations perturb miRNA–gene regulatory networks in cancer. In this study, we derived elaborate miRNA–target interaction networks in each cancer type by integrating empirically validated binding sites by CLIP-seq, and expression correlation between miRNAs and target genes (30). We demonstrate that somatic mutations are likely to occur selectively in miRNAs or target genes to perturb cancer hallmark-related functions. A mutually exclusive pattern is found for mutations in miRNAs and their targets. Driver mutations that significantly perturb miRNA regulatory networks in cancer are further identified. We show that driver mutations could exert their functions by impairing RNA binding affinity, resulting in alteration of target gene expression profiles. Intriguingly, driver gene targets that are significantly downregulated in cancer often function as tumor suppressors during cancer progression. Our study provides a valuable resource for systematic investigation of the functional impact of somatic mutations on miRNA regulation in cancer.

Argonaute (AGO) proteins are RNA-binding proteins (RBP) and essential components of the RNA-induced silencing complex (RISC), which is the molecular machinery for miRNA-induced silencing. High-throughput sequencing of immunoprecipitated RNAs after cross-link (CLIP-seq) to AGO proteins provides powerful ways to trace the footprints of miRNA binding sites (31). We downloaded peaks of 36 AGO CLIP-seq datasets from starBase v2.0 (32) and considered these binding sites as evidence of physical interactions between miRNAs and target genes (Fig. 1A). All the peaks of AGO CLIP-seq data were merged together based on their genomic locations to obtain a global pool of miRNA–gene physical interactions. Peaks that were overlapped with at least 1bp were merged together to provide evidence of permissive regions for miRNA binding. A total of 434,701 merged peaks were finally obtained after converting to hg38 genome assembly and ready for feeding to software prediction of precise sites.

Although the AGO CLIP-seq peaks could capture the footprints of miRNA binding sites, which miRNAs could bind to specific target sites still remains unclear. Agarwal and colleagues showed that many noncanonical sites detected by cross-linking method do not mediate repression despite binding the miRNAs whereas the vast majority of functional sites are canonical and can be identified by TargetScan (v7.0; ref. 30). To identify mGIs and functional binding sites on the mRNA, 16,347,639 interactions between all the miRNAs in miRBase release 21 and target genes in whole human genome were downloaded from TargetScan (v7.0; ref. 30). Binding sites of each miRNA–gene pair were further intersected with AGO-merged peaks and we finally obtained 279,924 mGIs among 2,586 mature miRNAs and 4,198 target genes as global reference in 3′ UTRs.

TargetScan (v7.0) predicts miRNA targets by searching for conserved sites six to eight nucleotides long that match to the seed region of 3′ UTR only (30), whereas RNA22 (v2) and RNAhybrid (v2.1.2) are seedless prediction tools available for detection of target sites in CDS, 5′ UTR, and 3′ UTR.

An assumption made by using TargetScan (7.0; ref. 30) to identify mGIs is that targets are only detected in 3′ UTRs, so miRNA targets in coding sequences (CDS) and 5′ UTRs are not detected by TargetScan (v7.0; ref. 30). To overcome this limitation, we downloaded ∼250M (million) predicted miRNA target sites from RNA22 (v2; ref. 33) for the 2,588 mature human miRNAs from miRBase (release 21; ref. 34), which targeted 140M sites in CDSs, 84M sites in 3′ UTRs, and 23M sites in 5′ UTRs from ENSEMBL78 (Fig. 1B; Table 1). To obtain more reliable sites, we overlapped these sites with AGO CLIP-seq sites collected from starBase (v2.0; ref. 32), which resulted in 2.7M sites in CDSs, 1.9M sites in 3′ UTRs, and 0.3M sites in 5′ UTRs (4.9M total sites; Fig. 1B; Table 1). These 4.9M target sequences and all 2,588 miRNA sequences were fed into RNAhybrid (v2.1.2; ref. 35) using default parameters with no helix constraint to ensure seedless predictions. Only the 0.15M sites participating in significant predictions (P < 0.05) by RNAhybrid (v2.1.2; ref. 35) and RNA22 (v2; ref. 33), and overlapping with AGO sites were kept.

Because RNA22 (v2) and RNAhybrid (v2.1.2) differ in that RNA22 (v2) is a pattern-based prediction tool that identifies miRNA targets through associating miRNAs with “target islands” (33), whereas RNAhybrid (v2.1.2) identifies miRNA targets through identifying energetically favorable hybridizations between miRNA and targets (35), we compared minimum free energies (MFE) of the 4.9M reliable sites overlapped by two algorithms in all types of sites (3′ UTR, CDS, and 5′ UTR) and calculated the Pearson correlation coefficients. We observed highly significant correlations (P < 2.2 × e⁻¹⁶) between RNA22 and RNAhybrid in all three types of sites, where the correlation is 0.53 for 3′ UTR sites, 0.58 for 5′ UTR sites and 0.52 for CDS sites (Fig. 1C). To further evaluate the consistency of predicted sites in CDS and 5′ UTR by two algorithms, we overlapped the significant sites detected by RNA22 (v2) and RNAhybrid (v2.1.2).

To obtain the functional mGIs in specific cancer type, we considered the expression correlation among samples for each mGI of global reference in the cancer type (Supplementary Table S1). The mGIs obtained from CDS and 5′ UTRs were then combined with mGIs from 3′ UTRs as the global reference for further analysis (Supplementary Table S2). Gene expression quantification by RNA-seq data and isoform expression quantification of mature miRNAs by miRNA-seq data of 33 cancer types were downloaded from the TCGA. These data were collected with written informed consent from the patients. This study was conducted in accordance with recognized ethical, and approved by an institutional review board. For each cancer type, datasets with more than 50 common samples (Supplementary Table S1) between RNA-seq and miRNA-seq data were extracted and Spearman's rank correlation of expression for each pair of miRNA and target gene was calculated.

where ρ_i is the rank correlation coefficient between the expression of miRNA_i and gene_i in the i_th mGI. The ranks of miRNA_i and gene_i expression in k_th sample are represented as rm_ik and rg_ik. N is the total number of samples in a cancer type. To control the FDR, P values of correlations for mGIs in each cancer type were further adjusted by Benjamini and Hochberg method. Finally, mGIs with ρ < 0 and FDR < 0.05 were considered as functional and reliable interactions in each cancer type (Supplementary Tables S3–S5).

To validate the mGIs identified from cancer samples, we built a benchmark dataset by collecting various sources of experimentally validated mGIs. Both validated positive and negative mGIs were collected from mirRecords (36), miRTarBase (37), TargetMiner (38), and TarBase (39). After mapping gene IDs to approved HGNC gene symbols and miRBase IDs (34) and integrating redundant mGIs from different sources, we finally retrieved 5,577 validated mGIs including 5,139 positive mGIs, which are validated as functional and 438 negative mGIs, which are validated as non-functional. Functional impacts of mutations were not considered during validation, as they would not be consistently common enough across cancer samples to skew the Spearman coefficient and affect mGI identification. Considering the incomplete nature for both benchmark and our predicted datasets, 166 common mGIs involved in our predictions were extracted for evaluation. Contingence table of prediction result in 3′ UTR was generated, and Fisher exact tests were calculated on the basis of the number of true positive (TP), true negative (TN), false positive (FP), and false negative (FN). Moreover, evaluations of performance were also calculated as below:

The cancer mutation dataset of 33 cancer types was achieved from “Multi-Center Mutation Calling in Multiple Cancers” (MC3) and was produced using six different algorithms on data from TCGA. A total of 3,518,261 mutations in 9,819 cancer samples were included in our analysis.

We mapped cancer mutations to the specified region by bedtools (40). To calculate the mutation probability in a genomic region, we first defined mutation rate of mut_i as: R(mut_i) = N(mut_i) / N, where N(mut_i) is the number of samples with mutation_i in a cancer and N is the number of all samples in a cancer. Then mutation probability of a region r was calculated as:

In brief, P(r) was calculated by summing up the mutation rate of all mutations in the region r and normalized by the length of region. In this way, we could compare the level of mutation probability for regions with different lengths and different number of mutations. To evaluate the level of mutation probability in mature miRNAs, the mature miRNA regions, seed regions (2–8 bp of each miRNA), upstream and downstream flanking 50bp regions of mature miRNAs were calculated for comparison.

For each mutated mGI, mutations could act on either the miRNA or target gene, or both of them. In this case, we classified mutated mGIs as miRNA-mutation mGIs, target-mutation mGIs, and dual-mutation mGIs. To find out whether the mutations on mGIs work in a synergetic way or in an exclusive way, we evaluated the mutation exclusivity of mGIs based on the occurrence of dual-mutation mGIs. In this way, mutation exclusivity on miRNA–gene pair was tested by hypergeometric distribution based on the number of all mutated mGIs, miRNA-mutation mGIs, target-mutation mGIs, and dual-mutation mGIs in each cancer type.

where N is the number of mutated mGIs, K is the number of mGIs with miRNA mutation, which equals m+d, M is the number of mGIs with target mutation, which equals t+d, and X is the number of dual mutation mGIs. In the case, m is the number of miRNA exclusively mutated mGIs, t is the number of target exclusively mutated mGIs, and d is the number of dual-mutation mGIs.

For each mGI_j identified in a type of cancer_i, first all the mutations within mature miRNA_j and binding sites of target gene_j were considered as candidate driver mutations. For each mutation_k on mGI_j, we searched samples with mutation_k in cancer_i and integrated them with gene expression of miRNA_j and gene_j in mutated samples. To identify the driver mutations that could affect miRNA–mRNA binding, we assumed the driver mutation on either miRNA or target site could alter the inhibitory role of miRNA_j and cause abnormal expression of target gene_j, which is isolated from the distribution of nonmutated control samples. Here two sets of control samples were considered, thus nonmutated cancer samples in 25 cancer types and nonmutated normal samples in five cancer types. For each mutated mGI_j in cancer_i, we used the linear regression model L_j to fit the expression distribution of control samples and calculated the prediction interval with a probability of 0.95.

where g_j and m_j are the expression of gene_j and miRNA_j among control samples. β₀ and β₁ are estimated parameters trained by expression of control samples, where only the models with significant P value (P < 0.05) are considered as successful. In this way, the control samples should fall into the predictive confident interval whereas samples with driver mutations on mGI_j should fall out of it (Supplementary Table S6).

We calculated the MFE for wild-type mGIs and driver-mutated mGIs (Supplementary Table S6) using RNAhybrid (41). The energy change of mGI altered by driver mutation is defined as:

The lower the ΔE is the more stable, as it is for the binding of miRNA–mRNA. Thus, it can be more easily for the miRNA to carry out its inhibitory function and downregulate the expression of target gene more effectively.

When calculating the proportions of mGI with changed or unchanged MFE, we defined mGIs with |ΔE | in top 30% as energy changed and all the other 70% as unchanged. Similarly for the gene expression, mGIs with |ΔExpression | in top 30% was defined as expression changed and the others as unchanged.

The differential impact of mutations in miRNAs compared with mRNA targets was further assessed using general linear mixed modeling. We defined our model L_k as:

where ΔE_k, M_k, and e_mk represent the change in free energy, the site of the mutation (miRNA vs. mRNA), and error for mGI_k retained after RNAhybrid analysis, respectively. β₀ and β₁ are estimated parameters trained on the mGIs, and b_C, 0c represents a control term for energy changes in the different cancer types. β₁ was used to assess which type of mutation was more destabilizing, as β₁ > 0 would suggest that mRNA mutations had a greater impact whereas β₁ < 0 suggested that mRNA mutations would have a reduced impact. Significance was assessed using the t value associated with β₁.

Certain cancer type exhibit differential mutation patterns between different sexes (42). To investigate the potential for such a bias in mGIs, we counted patients with mGI driver mutations for each gene in each cancer. To assess sex bias, we selected cancers that featured at least 1 male and 1 female patient. Within each cancer we selected genes that were part of an mGI driver mutation in at least 2 patients. With these criteria, 90 unique combinations of cancer types and genes were studied across 15 TCGA cancer cohorts.

Network visualization was conducted by the software Cytoscape. The layout of driver network was grouped by the classification of driver gene nodes, which are “mutated miRNAs,” “nonmutated miRNAs,” “mutated genes,” and “nonmutated genes.” Node degree is the number of interactions for each node in the global network.

To investigate the functional importance of the driver genes and mutations, functional enrichment analysis of the driver genes was carried out to investigate whether they were enriched in cancer hallmarks. The gene list of each cancer hallmark was obtained from one of the previous studies. We used hypergeometric test for exploring the statistical significance and the P value was computed as follow:

where N is the number of genes in the whole genome, of which, K genes were involved in the functional category under investigation (such as cancer hallmarks), and the number of candidate driver genes is M, of which, t genes were involved in the functional category.

Differential expression analysis was carried out for all genes in cancer types having both tumor and normal samples. Count of raw reads in each gene was downloaded in “HTSeq-Counts” format from TCGA. R package “DESeq2” was used to carry out differential analysis. Genes with reads count smaller than 10 in total samples were filtered out. Finally, genes with |fold change|>2 and adjusted P value <0.05 were considered as significantly differentially expressed.

For a given set of N genes, we calculated the proportion of significantly downregulated, upregulated, and nondifferentially expressed genes, which are denoted as D, U, and M. As the definition, D+U+M = 1. To test whether the observed values of D, U, and M are significant, we randomly picked N genes from all genes in all cancer types and calculated the proportions of them as d, u, and m in each random case. By repeating 10,000 times, we got the random distribution of variable d, u, and m. By comparing the D, U, and M with the random distributions of d, u, and m, we got the significance for each classification of differentially expressed genes.

Mutation information (MC3), gene expression RNA-seq data, miRNA-seq data as well as patient clinical features (survival and gender) in this study were obtained from The Cancer Genome Atlas Project (TCGA) project available at the Genomic Data Commons portal at https://portal.gdc.cancer.gov/. AGO CLIP-seq datasets are available from starBase v2.0 (32). miRNA details are available in miRBase (release 21; ref. 34), and miRNA target gene sequences are available in TargetScan (v7.0; ref. 30), RNA22 (v2; ref. 33), and RNAhybrid (v2.1.2; ref. 35). Validated positive and negative mGIs were compiled from mirRecords (36), miRTarBase (37), TargetMiner (38), and TarBase (39) to build a benchmark dataset. The data generated in this study are available in a user-friendly website resource, “mGI-map” (https://mgimap.dellmed.utexas.edu).

To analyze miRNA–gene regulation across cancer types, we proposed a computational method for identifying active miRNA–gene regulation in specific cancer types (Fig. 1A). In brief, high-throughput sequencing of immunoprecipitated RNAs after cross-link (CLIP-seq) to AGO proteins was used to identify endogenous genome-wide interaction maps for miRNAs (31) and bioinformatics tools were developed to infer specific target sites each miRNA could bind to (30, 33, 35). In a previous study, Agarwal and colleagues showed that many noncanonical sites detected by cross-linking methods did not mediate gene repression despite their ability to bind miRNAs whereas the vast majority of functional sites were canonical and could be identified by TargetScan (v7.0; ref. 30). In this study, AGO CLIP-seq data from starBase v2.0 (32) were integrated with TargetScan 7.0 (30) to get mGIs in 3′ UTR, and were integrated with RNA22 (v2; ref. 33) and RNAhybrid (v2.1.2; ref. 35) to get mGIs in CDS and 5′ UTR regions (Supplementary Table S7). We found predictions for either CDS or 5′ UTRs by RNA22 (v2) and RNAhybrid (v2.1.2) were highly similar (Supplementary Figs. S1A and S1B). To evaluate the consistency of predicted sites in 3′ UTR by different algorithms, we overlapped the significant sites in 3′ UTR detected by RNA22 (v2), RNAhybrid (v2.1.2), and TargetScan (v7.0). We found TargetScan (v7.0) sites covered almost all of the sites detected by RNA22 (v2) and RNAhybrid (v2.1.2; Supplementary Fig. S1C).

We derived 314,121 mGIs among 2,588 mature miRNAs and 5,207 target genes as a global reference network, including 279,924 mGIs in 3′ UTR, 27,016 mGIs in CDS, and 7,181 mGIs in 5′ UTR regions (Supplementary Table S2). After the inference of the global mGI network, we integrated the miRNA and gene expression profiles in specific cancer types to identify active miRNA–gene regulation. Only negatively correlated miRNA–gene pairs (FDR < 0.05) were kept for further analysis. In total, 114,469 unique mGIs were obtained from 30 of 33 cancer types, including 107,068 in 3′UTR, 5,938 in CDS, and 1,463 in 5′ UTR regions (Fig. 1D). The Spearman correlation coefficients ranged from −0.07 to −0.83 with a mean of −0.23, as the correlation associated with the FDR cutoff changed on the basis of the cancer cohort. Next, we validated 3′ UTR results using a benchmark dataset of experimental results collected from various sources. The performance of our predictions was: recall rate = 93.5%, precision = 94.1%, F1 score = 93.8%, accuracy = 88.6%. In addition, our predictions exhibited a significant overlap with the benchmark dataset (Supplementary Fig. S1D: OR, 6.36; P = 0.015). These results indicate that integration of multi-omics datasets enables us to identify cancer relevant miRNA-gene regulation at the systems level.

Comparing miRNA–gene regulatory networks across cancer types, we identified marked rewiring in the miRNA regulatory programs among different cancer types, with a unique “on/off” switch depending upon the cancer type. We found that ∼42.47% of 3′ UTR mGIs, 53.22% of CDS mGIs, and 56.87% of 5′ UTR mGIs occurred only in one cancer type, but only 1.55% of 3′ UTR mGIs, 0.72% of CDS mGIs, and 0.48% of 5′ UTR mGIs were observed in greater than 10 cancer types (Fig. 1D; Supplementary Figs. S2A and S2B). Notably, most of the active mGIs reside in 3′ UTR sites if compared with the number in CDS and 5′ UTR sites. We also observed that the distribution of cancer general (common) mGIs (N ≥ 10 cancer types) in 3′ UTR sites was more homogeneous than that of total mGIs (Supplementary Figs. S3A and S3B). The low conservation of miRNA–gene regulation could be explained in part by the cancer-specific expression of miRNAs or target genes. This indicates although the proportion of general mGIs is small, they are relatively stable across cancer types. In contrast, the number of cancer specific mGIs increased as the total number of mGIs grew (Supplementary Fig. S3C). Furthermore, hierarchical clustering analysis was conducted based on the occurrence of cancer frequent mGIs in 3′ UTR. Several cancer types with similar origins were clustered together (Supplementary Fig. S3D). OV (ovary) and CESC (cervix) were clustered into gynecologic system, showing low occurrences of cancer frequent mGIs. KIRP and KIRC (kidney) were clustered into urologic system, showing medium occurrences of frequent mGIs. COAD, READ (colorectal), and LUAD, LUSC (lung) were separately clustered into gastrointestinal and thoracic systems, respectively, showing relatively high occurrences of frequent mGIs. All these results indicate the intricate functional roles of miRNA regulation in cancer. Moreover, we performed functional enrichment analysis using the target genes from 3′ UTR that were observed in cancer frequent mGIs. We found that these genes were involved in general cancer-related functions, including FoxO signaling pathway, TGFβ signaling pathway and cell adhesion-related functions (Supplementary Fig. S3E). Taken together, our integrative analysis reveals cancer-specific miRNA–gene regulation, providing a valuable resource for mechanistic investigation of the function of miRNAs in cancer.

miRNAs play important roles in the development and progression of cancer. Several miRNA-associated mutations have been identified before, but systematic analysis of mutational landscape on miRNA regulation is still lacking. In this study, cancer somatic mutations obtained from the TCGA project were mapped to miRNAs and their target genes, or other regions in the human genome. We found that the mutation probability in miRNAs was greater than that in mRNAs (Fig. 2A, P < 2.2e−16, Wilcoxon rank-sum test). Moreover, we found that the genomic regions of miRNAs also had a higher mutation probability than upstream or downstream flanking regions (Fig. 2A, P < 2.2e−16). Specifically, the mutation probability of the seed regions within miRNAs was greater than other parts within miRNAs or the flanking regions (Supplementary Fig. S4A, P < 2.2e−16), suggesting the seed regions of miRNAs are required for miRNA regulation. In addition, we identified thousands of mutations located in target coding genes of miRNAs. Compared with the flanking regions, the miRNA binding sites in 3′ UTR, CDS, and 5′ UTR of target genes all exhibited a higher mutation probability (Fig. 2B–D; Supplementary Fig. S4B). These results indicate that mutations are likely to occur in the genomic regions that are critical for miRNA binding.

Next, we explored the signaling pathways that were likely perturbed by these mutations. As the hallmarks of cancer comprise biological capabilities acquired during the multistep development of human tumors, we evaluated whether these mutations could possibly perturb cancer hallmark related functions. As shown in Supplementary Fig. S4C, we found that the mutated miRNA target genes in 3′ UTRs were significantly enriched in various types of cancer hallmarks, such as regulation of cell proliferation and DNA repair. To further explore if these mutations could play a functional role in perturbing miRNA–gene regulation, we estimated the functional impact of somatic mutations identified in 3′ UTRs. Several methods have been proposed to assess the effects of mutations on protein function, and these methods are complementary to each other. We herein used Combined Annotation–Dependent Depletion (CADD; ref. 43), which is a framework that integrates multiple annotations into one metric, to explore the functional impact of mutations in target genes of miRNAs. For each cancer type, we randomly selected the same number of mutations as background controls. We then compared the CADD scores of mutations in target genes with those of randomly selected mutations. In the majority (23/31) of cancer types, mutations in target genes had significantly higher scores (Wilcoxon rank-sum test P < 0.05) than random controls (Supplementary Fig. S4D), suggesting that the identified miRNA target mutations could have deleterious effects in cancer. Evolutionary conservation of mutated residues has also been demonstrated to reflect the functional importance of the mutational event (44). We therefore explored the conservation feature of the mutations in target genes. We found that positions harboring target gene mutations were more likely to be conserved than positions harboring randomly selected mutations in most (29/31) cancer types (Supplementary Fig. S4D). All these results indicate that there are widespread mutations enriched within miRNA sites and their target genes, which could functionally perturb the binding of miRNAs, thus playing critical roles in cancer.

It has been observed that sets of genes that are co-involved in the same cancer pathways tend not to harbor mutations together in the same patient, which is called mutual exclusivity (45, 46). However, it remains elusive if the mutual exclusivity hypothesis could be extended in the context of miRNA regulation. We therefore assessed whether the same patients with cancer could simultaneously carry miRNA mutations and target gene mutations. To do so, we calculated the number of mGIs in 3′ UTR and CDS regions for three mutation perturbation models. The first one was that mutations located only in miRNAs to perturb the interactions; the second one was that mutations located only in target genes to perturb the interactions; and the third one was that mutations were located in both miRNAs and their target genes (Fig. 3A). We found that a vast majority of mGIs harbored miRNA mutations for 3′ UTR mGIs (Fig. 3B, left) or target mutations for CDS mGIs (Fig. 3B, right), and only a small fraction of mGIs had both miRNA and target mutations (Fig. 3B). 5′ UTR mGIs were too few to make statistical calculations. In addition, we found that patients in different cancer types exhibited different regulatory patterns (Fig. 3C and D; Supplementary Figs. S5A and S5B). For example, in 3′ UTR mGIs the majority of patients with ACC, LAML, PAAD, and READ carried miRNA mutations. In contrast, TGCT, THYM, and PCPG patients were likely to have target binding site mutations (Fig. 3C). In CDS mGIs, most of the cancer types carried target mutations except for ESCA, PAAD, READ, and THYM (Fig. 3D). In summary, we found that the mutual exclusivity of miRNA versus target mutations was more significant across cancer types than expected by chance (Fig. 3C and D, P < 0.05), which occurred in 23 of 25 cancer types for 3′ UTR mGIs and 14 of 23 cancer types for CDS mGIs.

Specifically, we identified mutually exclusive 3′ UTR mGIs that had been reported in various types of cancer. For example, the let-7 family is a conserved family of miRNAs and the deregulation of this family was observed in many cancer types, such as breast cancer, lung cancer, and melanoma. Here, several mutations were observed in let-7f-1–3p, let-7a-5p, and let-7b-3p. In addition, several mutations located within target genes were also observed in patients, such as LIFR, SPIRE1, and E2F3. However, no patient was found with both miRNA mutations and target gene mutations. These observations suggest that the mutations selectively target miRNAs or target genes to perturb miRNA–gene regulation in cancer, which is consistent with the pathway redundancy model in cancer.

Distinguishing driver mutations from passenger mutations in individual genomes is a big challenge in cancer genomic studies. In this study, we integrated the mGI regulatory networks with target gene transcriptomic networks to identify candidate driver mutations in each cancer type. Our hypothesis was that candidate driver mutations were likely to perturb mGIs, which could be reflected by a significant change in the expression levels of target genes (Fig. 4A). We thus used a linear regression model to fit the expression distribution among non-mutated control samples and identified driver mutations that elicited significant deviations from the normal correlation between miRNAs and their target genes (see details in Materials and Methods). In total, we identified 89 driver mutations for 3′ UTR mGIs, 51 driver mutations for CDS mGIs, and 8 driver mutations for 5′ UTR mGIs across 20 cancer types (Fig. 4B; Supplementary Figs. S6A and S6B; Supplementary Table S6). These driver mutations perturbed 154 mGIs in 3′ UTRs, 57 mGIs in CDS regions, and 8 mGIs in 5′ UTRs (Fig. 4B). Among the 89 driver mutations in 3′ UTRs, 75 of them occurred in target genes and 14 of them occurred in miRNAs. In addition, 9 of 75 driver mutations in target genes had indeed been characterized as cancer mutations in the Catalogue of Somatic Mutations in Cancer (COSMIC) database and in this study we inferred their possible functions through miRNA regulatory networks.

Mutations at mGI interfaces could impair the binding of miRNAs to target genes by changing the MFE of binding sites, thus further tuning the regulation of target gene expression. We calculated the changes of MFE using RNAhybrid (41) and changes of target gene expression between wild-type and mutated 3′ UTR and CDS mGIs. 5′ UTR mGIs were again too few to make statistical calculations. By comparative analysis of different types of mGIs, we found that the energy changes were significantly correlated with the changes in target gene expression levels by Fisher exact test for 3′ UTR mGIs (P = 0.002; OR, 3.3) but not for CDS mGIs (P = 0.793; OR, 1; Fig. 4C; Supplementary Table S6). We next speculated that if driver mutations disrupt miRNA–target binding, one would expect the repression of target gene expression by miRNAs to be reversed (i.e., increased). Consistently, we found that target genes in 3′ UTR mGIs with increased expression (e.g., top 5%) showed a significantly higher alteration of MFE than those with decreased expression (e.g., bottom 5%) by Wilcoxon rank-sum test (P = 0.038, Fig. 4D, left). However, we did not observe significantly different alteration of MFE for target genes in CDS mGIs (P = 0.575, Fig. 4D, right). This difference could be due to the repression efficiency of target sites in 3′ UTR and CDS regions. Although both 3′ UTR and CDS sites are functionally conserved across mammalian species, 3′ UTR sites are more conserved and appear to have more repressive effects than CDS sites (47). This may be a result of the complementary functions of 3′ UTR and CDS sites: 3′ UTR sites more efficiently trigger mRNA degradation whereas CDS sites more efficiently inhibit translation (48, 49). Fang and colleagues (50) have also investigated an additional synergistic function of CDS sites such that the presence of CDS sites appear to enhance the regulation of mRNAs at their 3′ UTR sites. Anyway, it is interesting that driver mutations in 3′ UTR showing greater influences through disruption of binding sites than those in CDS regions.

Furthermore, we investigated whether miRNA mutations and target mutations impact binding affinity in the same manner, we found that target mutations tended to have a greater impact on destabilizing free energy than miRNA mutations in 3′ UTR mGIs (P < 0.01, Supplementary Fig. S7). Together, these results suggest that driver mutations are likely to perturb mGIs in cancer through their functional influence on miRNA binding affinity to target genes.

Finally, recent studies yielded gender-specific regulatory networks in many human tissues (42). Our analysis of sexual dimorphism in 3′ UTR mGI mutation patterns yielded further insights. Across 15 TCGA cancer cohorts examined, we found that only EGFR exhibited a significant bias (P < 0.05) based on biological sex in LGG (Supplementary Fig. S8A; Supplementary Table S8). In this cohort, only females (9/285; 3.2%) possessed an EGFR mutation in mGI interaction sites. Although EGFR was only found to have a statistically significant bias in LGG, bias in EGFR mutation rates in SKCM and HNSC also ranked quite highly. Interestingly, the bias in EGFR-associated mGI mutations changes with the cancer type, as males were found to be at higher risk than females in SKCM (OR, 2.23) and lower risk than females in HNSC (OR, 0.27; Supplementary Fig. S8A). In addition, 8 of 9 mGI driver mutations in EGFR in LGG were missense mutations. Of these, five were clustered in the tyrosine kinase domain, and all 9 mutations were contained within or directly upstream or downstream of this domain (Supplementary Fig. S8B). Altogether, these findings indicate that sex bias associated with mGI mutations warrants further exploration.

To analyze the structure and properties of perturbed mGI networks by driver mutations, we extracted all mGIs from 3′ UTRs, CDS, and 5′ UTRs involving both miRNAs and target genes that harbored driver mutations. As shown in Fig. 5A, the bipartite mGI networks had four types of node components: mutated miRNAs, nonmutated target genes, nonmutated miRNAs, and mutated target genes. mGIs could be perturbed by mutations in the miRNAs (mostly interactions between mutated miRNAs and nonmutated targets), or perturbed by mutations in the target genes (mostly interactions between nonmutated miRNAs and mutated targets). We then examined the functional pathways impacted by the driver mutation-induced mGI network perturbations. Through analysis of cancer hallmark signatures, we found that perturbed target genes were enriched in “insensitivity to antigrowth signals” and “limitless replicative potential” (Fig. 5B). These two functions could help explain the increased growth and escaping from antigrowth signals in the cancer cells carrying these driver mutations. Furthermore, we investigated the expression profiles of all the target genes in the perturbed mGI networks. As shown in Fig. 5C and D, 60 mutated target genes in perturbed networks were differentially expressed in at least one cancer type, and 48 nonmutated target genes (but targeted by mutated miRNAs) were differentially expressed in at least one cancer type. These results indicate that the target genes perturbed by driver mutations are associated with cancer development.

Intriguingly, for the miRNAs carrying driver mutations (shown in top left of Fig. 5A), many of them were from the let-7 family including hsa-let-7a-2–3p, hsa-let-7a-3p, hsa-let-7a-5p, hsa-let-7b-3p, hsa-let-7c-3p, hsa-let-7c-5p, hsa-let-7f-1–3p, hsa-let-7f-2–3p, hsa-let-7f-5p, and hsa-let-7g-3p. As a well-known miRNA family with functions of tumor suppression, aberrant regulation of the let-7 family had been reported in many cancer types (51). These included increased cellular proliferation in lung cancer (52), increased proliferation and migration in liver cancer (53), increased invasion and metastasis in gastric adenocarcinoma, lymph node metastasis in breast cancer (54), and so on. Our study further recaptures the functional importance of the let-7 family in cancer by pointing out that the let-7 family could rewire molecular signaling through mutations within miRNAs themselves to dysregulate their targets.

To take a close investigation on the mutated 3′ UTR, we made a volcano plot for them in their matched cancer types and found 26.3% of the mutated 3′ UTR targets were downregulated (Fig. 6A). To explore the significance of these proportions, we randomly picked the same number of genes as the mutated target gene set in all cancer types 10,000 times and generated the random distribution as a background control. We found that the observed proportion of downregulated 3′ UTR target genes in reality was significantly higher than expected by chance (Fig. 6B, P < 0.001). Similarly, we found that 12.9% of the mutated CDS targets were downregulated (Fig. 6C), but this observed proportion of downregulated CDS target genes did not distinguish from random controls (Fig. 6D, P = 0.0925). In addition, the proportions of upregulated and nondifferentially expressed genes in 3′ UTRs were not significant (Supplementary Figs. S9A–S9C), which supports the conclusion that mutated 3′ UTR target genes are significantly downregulated. This result also matches with the different influences of driver mutations in 3′ UTR and CDS regions on target gene expression through disruption of binding sites. It indicates that 3′ UTR driver mutations are more likely to be involved in miRNA regulated pathways than CDS driver mutations.

By curation from literature, we found that most of the downregulated mutated 3′ UTR target genes were tumor suppressor genes (Supplementary Fig. S9D). This finding is corroborated by prior analyses that found a similar pattern (9). The top-mutated target genes were all cancer related and mostly annotated as inhibiting cancer progression. We further extracted the perturbed networks of these tumor suppressors (Supplementary Fig. S9E). Upstream miRNAs (red nodes) and their interactions with these tumor suppressors were perturbed by driver mutations. Although the other target genes of these upstream miRNAs were unperturbed, no significant functional pathways were enriched for unperturbed targets.

In the perturbed tumor suppressor networks, all the target genes involved and most of the miRNAs had previous evidence to be associated with cancer. For example, LIFR could confer a dormancy phenotype in breast cancer cells and also negatively regulate the metastasis of pancreatic cancer cells (55, 56). RBM24 was demonstrated to suppress cancer progression in nasopharyngeal carcinoma and also reported as a novel player in the p53 pathway (57, 58). KLF15 could inhibit the cell proliferation in gastric cancer (59). AKAP12 was reported as a tumor suppressor in glioblastoma and prostate cancer, and regulated by miRNA pathways (60, 61). CBX7 was reported as a tumor suppressor in colon, thyroid carcinomas, glioblastoma, and pancreatic cancer (60, 62, 63). It could play its role by modulating the KRAS gene and could be regulated by miRNAs (60, 62). For the upstream miRNAs of these tumor suppressors, hsa-miR-6881–3p was known to play its role through the p53-mediated ceRNAs network in hepatocellular carcinoma (64). hsa-miR-15b-5p was used to distinguish human ovarian cancer tissues from normal tissues (65). hsa-miR-193b-3p was found differentially expressed in hepatocellular carcinoma (66). hsa-miR-15a-5p was shown to decrease melanoma cell viability and could cause cell-cycle arrest at the G₀–G₁ phase (67). hsa-miR-769–3p was differentially expressed in patients with non–small cell lung cancer (NSCLC) harboring EGFR mutations (68). In conclusion, the mutated driver target genes perturbed through miRNA pathways play their roles probably as tumor suppressors and are downregulated in many cancer types thus accelerating the growth of cancer cells.

To help researchers apply the principles described in this work to any mGIs or mutations of interest, we developed a comprehensive and interactive web resource, mGI-map (https://mgimap.dellmed.utexas.edu). The features provided in the resource should serve as a guide for biologists interested in identifying miRNA–gene regulation in specific cancer types and understanding the consequences of mutations that perturb specific mGIs in cancer. All the data could also be downloaded for further analysis of the functional effect of mutations and miRNA–gene regulation (Supplementary Fig. S10).

In this study, we derived de novo cancer-specific mGIs with high quality by integration of in silico predictive models, AGO CLIP-seq data, and transcriptome networks. We built out these mGI networks of various cancer types as a framework for pan-cancer analysis of the functional consequences of somatic mutations on mGIs. Further analysis showed that the mutations on mGIs exhibited a mutually exclusive pattern. Mutual exclusivity was previously observed for mutations in gene members of the same pathways. Yeang and colleagues found that when a member of a pathway was altered, the selection pressure on the other members was diminished or even nullified (69). As a result, significantly less overlap in mutations of the other genes was expected, creating a mutually exclusive pattern between their alterations. Supporting this expectation, it was previously shown that some functionally related genes were mutated in a mutually exclusive manner in thyroid tumors (70) and in leukemia (71). In our study, we observed this pattern on mGIs in most of the cancer types.

To identify driver mutations that can perturb the mGIs, we compared the samples with mutations on mGIs with the fitted expression distribution in nonmutated samples and evaluated the deviations of each candidate mutations. After identification of driver mutations, we analyzed the alteration of MFE and target gene expression between driver mutant and wild-type mGIs. We found that the target gene expression was associated with the alteration of MFE, indicating that the driver mutations could perturb mGIs by influencing mGI binding sites.

The data and design of our study produce several limitations. Because AGO-CLIP data are derived from cell lines rather than live tumor samples, some mGIs may have been lost given the high context dependency of the interactions (72). In addition, our analysis grouped together the many miRNA isomers (isomiR) present in cancer as their source gene. This is significant due to prior evidence pointing to the differential involvement of isomiRs in cancer (73). Although this study offered increased statistical confidence in the identified mGIs, some resolution in the specificity of different isomiRs involved in the mGIs was lost. Telonis and colleagues found an abundant of isomiRs were also expressed in tissues and they were nearly 10 times higher than the reference sequences reported in miRBase from the analysis on TCGA dataset (73). In addition, around half of the most abundant isomiRs in various cancer types are not the reference miRNAs in the miRBase. The isomiR world makes the mGI regulation much more complicated and elaborated. Although we analyzed all the reference miRNAs in miRBase, we may still have missed ∼9/10 miRNA targets due to the fact of isomiRs. Furthermore, the elaborated forms of isomiRs raise more challenges in their target sites identification. The traditional CLIP protocols include an RNase digestion step and different RNase strengths would trim the endpoints of miRNA–target pair differently (74). Therefore the digestion step could affect the miRNA member of the pair at both 5′ and 3′ ends, which prevents the accurate identification of isomiR target sites. We do believe there was still an opening gap in the miRNA regulation world when considering the isomiRs, and we realize that has been exceeding our capacity in this research. These areas will offer key directions for further study from our presented mGIs on mGI-map.

The genetic changes that contribute to cancer tend to affect three main types of genes—proto-oncogenes, tumor suppressor genes, and DNA repair genes, which are called “drivers” of cancer. In our analysis, we have shown that miRNA-related driver mutations tend to play their roles on tumor suppressor genes by modulating the mGIs. Tumor suppressor genes are involved in controlling cell growth and division. Cells with certain alterations in tumor suppressor genes may divide in an uncontrolled manner. Concordantly, our analyses show that the identified target genes in the perturbed miRNA regulatory networks are enriched in the cancer hallmarks of “insensitivity to antigrowth signals” and “self-sufficiency in growth signals,” which are consistent with the function of tumor suppressors. In this study, although several known tumor suppressors have been revealed and perturbed by miRNA related driver mutations, such as LIFR, RBM24, PTGER3, several identified novel mutated driver genes such as LDLRAD4, AKAP5, KLHL6, and so on, could become potential biomarkers and therapeutic targets for tumor suppression.

J. Das reports grants from NIAID, NHGRI, and DoD during the conduct of the study and personal fees from Seromyx outside the submitted work. G.B. Mills reports grants, personal fees, nonfinancial support, and other support from Amphista, AstraZeneca, Chrysallis Biotechnology, GSK, ImmunoMET, Ionis, Lilly, PDX Pharmaceuticals, Signalchem Lifesciences, Symphogen, Tarveda, Turbine, Zentalis Pharmaceuticals, Catena, HRD assay to Myriad Genetics DSP to Nanostring, Adelson Medical Research Foundation, Breast Cancer Research Foundation, Komen Research Foundation, Ovarian Cancer Research Foundation, Prospect Creek Foundation, Nanostring Center of Excellence, Ionis (Provision of tool compounds), and Genentech during the conduct of the study. No disclosures were reported by the other authors.

X. Hua: Resources, data curation, formal analysis, validation, investigation, visualization, methodology, writing–original draft. Y. Li: Data curation, formal analysis, validation, visualization, methodology, writing–original draft. S.R. Pentaparthi: Resources, data curation, formal analysis, validation, investigation, writing–original draft. D.J. McGrail: Data curation, validation, writing–review and editing. R. Zou: Data curation, validation, investigation, visualization, writing–review and editing. L. Guo: Data curation, validation, investigation, writing–review and editing. A. Shrawat: Data curation, validation. K.M. Cirillo: Data curation, validation. Q. Li: Data curation, validation. A.V. Bhat: Resources, data curation. M. Xu: Resources, writing–review and editing. D. Qi: Resources, data curation, writing–review and editing. A. Singh: Resources, data curation, validation, investigation, visualization, writing–review and editing. F. McGrath: Resources. S. Andrews: Resources. K.L. Aung: Writing–review and editing. J. Das: Conceptualization, resources, data curation, formal analysis, supervision, funding acquisition, validation, investigation, methodology, writing–original draft, project administration, writing–review and editing. Y. Zhou: Data curation, writing–review and editing. A. Lodi: Writing–review and editing. G.B. Mills: Writing–review and editing. S.G. Eckhardt: Writing–review and editing. M.L. Mendillo: Writing–review and editing. S. Tiziani: Writing–review and editing. E. Wu: Conceptualization, resources, data curation, supervision, funding acquisition, validation, investigation, methodology, writing–original draft, project administration, writing–review and editing. J.H. Huang: Resources, data curation, writing–review and editing. N. Sahni: Conceptualization, resources, supervision, funding acquisition, investigation, writing–original draft, project administration, writing–review and editing. S.S. Yi: Conceptualization, resources, data curation, formal analysis, supervision, funding acquisition, validation, investigation, methodology, writing–original draft, project administration.

The authors are grateful to contributions from TCGA Research Network Analysis Working Group and also acknowledge the High-Performance Computing Cluster Facility at MD Anderson Cancer Center, the Biomedical Research Computing Facility at UT Austin, and Texas Advanced Computing Center (TACC) for assistance. This work was supported by the NIH grant R35GM133658 (to S. Yi), Komen Foundation grant CCR19609287 (to S. Yi), Early Career Investigator Grant by Ovarian Cancer Research Alliance Grant# 649968 (to N. Sahni), Department of Defense W81XWH-18-PRCRP-CDA CA181455 (to N. Sahni), NIGMS grant GM137836 (to N. Sahni), Cancer Prevention and Research Institute of Texas grant RR160093 (to S.G. Eckhardt), and NCI grants P50CA217685 and U01CA217842 (to G.B. Mills). S. Yi was also supported by 2021–2022 and 2022–2023 program in Oncological Data and Computational Sciences sponsored by The Joint Center for Computational Oncology between the Oden Institute, MD Anderson, and Texas Advanced Computing Center (TACC). G.B. Mills also received a kind gift from the Miriam and Sheldon Adelson Medical Research Foundation and Breast Cancer Research Foundation. N. Sahni is a CPRIT Scholar in Cancer Research with funding from the Cancer Prevention and Research Institute of Texas (CPRIT) New Investigator Grant RR160021. D.J. McGrail was supported by NCI grant K99CA240689. The authors would also like to acknowledge the support from BioRender in figure generation.

The publication costs of this article were defrayed in part by the payment of publication fees. Therefore, and solely to indicate this fact, this article is hereby marked “advertisement” in accordance with 18 USC section 1734.