Somatic mutations are a major source of cancer development, and many driver mutations have been identified in protein coding regions. However, the function of mutations located in miRNA and their target binding sites throughout the human genome remains largely unknown. Here, we built detailed cancer-specific miRNA regulatory networks across 30 cancer types to systematically analyze the effect of mutations in miRNAs and their target sites in 3′ untranslated region (3′ UTR), coding sequence (CDS), and 5′ UTR regions. A total of 3,518,261 mutations from 9,819 samples were mapped to miRNA–gene interactions (mGI). Mutations in miRNAs showed a mutually exclusive pattern with mutations in their target genes in almost all cancer types. A linear regression method identified 148 candidate driver mutations that can significantly perturb miRNA regulatory networks. Driver mutations in 3′UTRs played their roles by altering RNA binding energy and the expression of target genes. Finally, mutated driver gene targets in 3′ UTRs were significantly downregulated in cancer and functioned as tumor suppressors during cancer progression, suggesting potential miRNA candidates with significant clinical implications. A user-friendly, open-access web portal (mGI-map) was developed to facilitate further use of this data resource. Together, these results will facilitate novel noncoding biomarker identification and therapeutic drug design targeting the miRNA regulatory networks.

Significance:

A detailed miRNA–gene interaction map reveals extensive miRNA-mediated gene regulatory networks with mutation-induced perturbations across multiple cancers, serving as a resource for noncoding biomarker discovery and drug development.

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.

Identification of global mGIs in 3′ UTRs

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.

Figure 1.

The landscape of miRNA–gene regulatory networks across cancer types. A, The workflow to identify genome-wide miRNA–gene regulation in cancer. Firstly, global mGIs were identified by integrating sequence alignment and AGO-CLIP-seq datasets. Then miRNA and mRNA expression profiles were integrated to identify cancer-specific interactions. Somatic mutations were mapped to mGIs to identify candidate driver mutations. B, Target sites distribution in 3′ UTR, CDS, and 5′ UTR regions from RNA22 raw data; RNA22 and AGO binding overlapped sites; RNA22 significant (P < 0.05), RNAhybrid significant (P < 0.05), and AGO binding overlapped sites. C, MFE correlation between RNA22 and RNAhybrid in CDS, 5′ UTR, and 3′ UTR regions at all AGO binding sites. D, The proportion of active mGIs from 3′ UTR, 5′ UTR, and CDS regions in different numbers of cancer types. (A, Created with BioRender.com.)

Figure 1.

The landscape of miRNA–gene regulatory networks across cancer types. A, The workflow to identify genome-wide miRNA–gene regulation in cancer. Firstly, global mGIs were identified by integrating sequence alignment and AGO-CLIP-seq datasets. Then miRNA and mRNA expression profiles were integrated to identify cancer-specific interactions. Somatic mutations were mapped to mGIs to identify candidate driver mutations. B, Target sites distribution in 3′ UTR, CDS, and 5′ UTR regions from RNA22 raw data; RNA22 and AGO binding overlapped sites; RNA22 significant (P < 0.05), RNAhybrid significant (P < 0.05), and AGO binding overlapped sites. C, MFE correlation between RNA22 and RNAhybrid in CDS, 5′ UTR, and 3′ UTR regions at all AGO binding sites. D, The proportion of active mGIs from 3′ UTR, 5′ UTR, and CDS regions in different numbers of cancer types. (A, Created with BioRender.com.)

Close modal

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.

Identification of global mGIs in coding sequence and 5′ 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.

Table 1.

miRNA targets by different algorithms.

Algorithms5′ UTR3′ UTRCDSTotal
RNA22 raw 22558002 83939429 145391270 251888701 
RNA22 raw and overlapped AGO sites 329761 1928045 2673408 4931214 
RNA22 and RNAhybrid significant and overlapped AGO sites 12089 64839 68381 145309 
Algorithms5′ UTR3′ UTRCDSTotal
RNA22 raw 22558002 83939429 145391270 251888701 
RNA22 raw and overlapped AGO sites 329761 1928045 2673408 4931214 
RNA22 and RNAhybrid significant and overlapped AGO sites 12089 64839 68381 145309 

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−16) 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).

Identification of cancer type–specific mGIs

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.

formula

where ρi is the rank correlation coefficient between the expression of miRNAi and genei in the ith mGI. The ranks of miRNAi and genei expression in kth sample are represented as rmik and rgik. 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).

Validation of identified mGIs in 3′ UTR

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:

formula
formula
formula
formula

Cancer mutations

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.

Mutation probability

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 muti as: R(muti) = N(muti) / N, where N(muti) is the number of samples with mutationi in a cancer and N is the number of all samples in a cancer. Then mutation probability of a region r was calculated as:

formula

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.

Significance of mutation exclusivity on mGI

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.

formula

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.

Identification of driver mutations on mGI

For each mGIj identified in a type of canceri, first all the mutations within mature miRNAj and binding sites of target genej were considered as candidate driver mutations. For each mutationk on mGIj, we searched samples with mutationk in canceri and integrated them with gene expression of miRNAj and genej 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 miRNAj and cause abnormal expression of target genej, 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 mGIj in canceri, we used the linear regression model Lj to fit the expression distribution of control samples and calculated the prediction interval with a probability of 0.95.

formula

where gj and mj are the expression of genej and miRNAj among control samples. β0 and β1 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 mGIj should fall out of it (Supplementary Table S6).

Minimum-free energy change by driver mutations on mGI

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:

formula

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 Lk as:

formula

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

Sex bias across driver mutations in mGIs

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.

mGI network analysis

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.

Cancer hallmark and enrichment

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:

formula

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

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.

Random permutation test for differentially expressed genes

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.

Data availability

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).

Integrative method for construction of pan-cancer miRNA regulatory networks

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.

miRNA mutations and target gene mutations 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. 2BD; Supplementary Fig. S4B). These results indicate that mutations are likely to occur in the genomic regions that are critical for miRNA binding.

Figure 2.

miRNA and target gene mutation profiles. A, The distribution of mutation probability for miRNA-related regions, including upstream, downstream of miRNA regions, miRNAs, seed regions, and mRNA regions. B, The distribution of mutation probability for 3′ UTR binding sites, including upstream, downstream, miRNA binding sites in targets, and all target mRNA regions. C, The distribution of mutation probability for CDS binding sites, including their upstream, downstream, target binding sites, and all target mRNA regions. D, The distribution of mutation probability for 5′ UTR binding sites, including their upstream, downstream, target binding sites, and all target mRNA regions.

Figure 2.

miRNA and target gene mutation profiles. A, The distribution of mutation probability for miRNA-related regions, including upstream, downstream of miRNA regions, miRNAs, seed regions, and mRNA regions. B, The distribution of mutation probability for 3′ UTR binding sites, including upstream, downstream, miRNA binding sites in targets, and all target mRNA regions. C, The distribution of mutation probability for CDS binding sites, including their upstream, downstream, target binding sites, and all target mRNA regions. D, The distribution of mutation probability for 5′ UTR binding sites, including their upstream, downstream, target binding sites, and all target mRNA regions.

Close modal

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.

Mutually exclusive mutations in miRNA–gene regulation

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.

Figure 3.

The mutual exclusivity of mutations in mGIs. A, Three models to illustrate the mutual exclusivity of mutations in mGIs. B, The proportion of 3′ UTR and CDS mGIs with different types of mutations. Purple, mGIs with miRNA mutations; green, mGIs with target mutations; blue, mGIs with both miRNA and target mutations. C and D, Pie plots show the proportion of samples with different types of mutated 3′ UTR mGIs (C) and CDS mGIs (D). I, the proportion of samples with miRNA mutations; II, the proportion of samples with target mutations; III, the proportion of samples with both miRNA and target mutations. Dot plots at the bottom indicate the significance level for mutual exclusivity test.

Figure 3.

The mutual exclusivity of mutations in mGIs. A, Three models to illustrate the mutual exclusivity of mutations in mGIs. B, The proportion of 3′ UTR and CDS mGIs with different types of mutations. Purple, mGIs with miRNA mutations; green, mGIs with target mutations; blue, mGIs with both miRNA and target mutations. C and D, Pie plots show the proportion of samples with different types of mutated 3′ UTR mGIs (C) and CDS mGIs (D). I, the proportion of samples with miRNA mutations; II, the proportion of samples with target mutations; III, the proportion of samples with both miRNA and target mutations. Dot plots at the bottom indicate the significance level for mutual exclusivity test.

Close modal

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.

Identification of driver mutations on mGIs

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.

Figure 4.

The candidate driver mutated mGIs across cancer types. A, The linear regression model is used to identify the candidate driver mutated mGIs in cancer. B, The statistics of driver mutations (left) and perturbed mGIs (right) in 3′ UTR, CDS, and 5′ UTR regions across cancer types. C, The proportion of mGIs in 3′ UTR (left) and CDS (right) regions with energy changed or not for expression changed/unchanged mGIs. Orange, energy changed; green, energy unchanged. D, The correlation between energy change and expression change for perturbed mGIs from 3′ UTR (left) and CDS (right) regions.

Figure 4.

The candidate driver mutated mGIs across cancer types. A, The linear regression model is used to identify the candidate driver mutated mGIs in cancer. B, The statistics of driver mutations (left) and perturbed mGIs (right) in 3′ UTR, CDS, and 5′ UTR regions across cancer types. C, The proportion of mGIs in 3′ UTR (left) and CDS (right) regions with energy changed or not for expression changed/unchanged mGIs. Orange, energy changed; green, energy unchanged. D, The correlation between energy change and expression change for perturbed mGIs from 3′ UTR (left) and CDS (right) regions.

Close modal

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.

mGI network perturbations by driver mutations

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.

Figure 5.

Network perturbation by driver mutations. A, Gene regulatory networks perturbed by driver mutations in 3′ UTR, CDS, and 5′ UTR mGIs. Yellow edges, perturbed mGIs; gray edges, unperturbed mGIs. Red nodes, mutated miRNAs; orange nodes, nonmutated miRNAs; blue nodes, mutated target genes; cyan nodes, nonmutated genes. B, Cancer hallmark enrichment of driver target genes in 3′ UTR, CDS, and 5′ UTR mGIs. Hallmarks with OR (x-axis) >1 and P value (y-axis) <0.05 are significantly enriched. C and D, Differential expression of mutated driver target genes (C) and nonmutated driver genes (D) from 3′ UTR, CDS, and 5′ UTR mGIs. Each row represents a driver gene that is differentially expressed in at least one cancer type. Each column represents a cancer type. The occurrence of each point indicates the gene is significantly differentially expressed in that cancer, which means the fold change of expression is larger than 2 or smaller than 1/2 and the P value calculated by DESeq2 is smaller than 0.05. Red indicates the gene is upregulated in cancer samples versus normal samples, whereas blue indicates downregulation.

Figure 5.

Network perturbation by driver mutations. A, Gene regulatory networks perturbed by driver mutations in 3′ UTR, CDS, and 5′ UTR mGIs. Yellow edges, perturbed mGIs; gray edges, unperturbed mGIs. Red nodes, mutated miRNAs; orange nodes, nonmutated miRNAs; blue nodes, mutated target genes; cyan nodes, nonmutated genes. B, Cancer hallmark enrichment of driver target genes in 3′ UTR, CDS, and 5′ UTR mGIs. Hallmarks with OR (x-axis) >1 and P value (y-axis) <0.05 are significantly enriched. C and D, Differential expression of mutated driver target genes (C) and nonmutated driver genes (D) from 3′ UTR, CDS, and 5′ UTR mGIs. Each row represents a driver gene that is differentially expressed in at least one cancer type. Each column represents a cancer type. The occurrence of each point indicates the gene is significantly differentially expressed in that cancer, which means the fold change of expression is larger than 2 or smaller than 1/2 and the P value calculated by DESeq2 is smaller than 0.05. Red indicates the gene is upregulated in cancer samples versus normal samples, whereas blue indicates downregulation.

Close modal

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.

Driver mutations in 3′ UTRs primarily target tumor suppressor genes

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.

Figure 6.

Tumor suppressors are often targeted in perturbed miRNA regulatory networks. A and C, Volcano plots of 48 mutated 3′ UTR target driver genes (A) and 31 mutated CDS target driver genes (C) in the cancer type where the driver mutation occurs. X-axis shows the log2 value of expression fold change and y-axis shows the value of –log10 (FDR). Each point indicates the gene in the cancer type where the driver mutation is found. The color of each point indicates the cancer type. B, Mutated driver genes in 3′ UTR mGIs are significantly downregulated. Bar plot shows the random distribution of the proportion of downregulated genes in mutated driver genes. The smoothed curve is the probability density of random distribution. The dashed vertical line indicates the observed proportion of downregulated genes. D, Mutated driver genes in CDS mGIs are not significantly downregulated (P > 0.05). Legend is the same as B.

Figure 6.

Tumor suppressors are often targeted in perturbed miRNA regulatory networks. A and C, Volcano plots of 48 mutated 3′ UTR target driver genes (A) and 31 mutated CDS target driver genes (C) in the cancer type where the driver mutation occurs. X-axis shows the log2 value of expression fold change and y-axis shows the value of –log10 (FDR). Each point indicates the gene in the cancer type where the driver mutation is found. The color of each point indicates the cancer type. B, Mutated driver genes in 3′ UTR mGIs are significantly downregulated. Bar plot shows the random distribution of the proportion of downregulated genes in mutated driver genes. The smoothed curve is the probability density of random distribution. The dashed vertical line indicates the observed proportion of downregulated genes. D, Mutated driver genes in CDS mGIs are not significantly downregulated (P > 0.05). Legend is the same as B.

Close modal

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 G0–G1 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.

mGI-map: a web-based resource for mutation-mediated mGI network perturbations in cancer

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.

Note: Supplementary data for this article are available at Cancer Research Online (http://cancerres.aacrjournals.org/).

1.
McGranahan
N
,
Favero
F
,
de Bruin
EC
,
Birkbak
NJ
,
Szallasi
Z
,
Swanton
C
.
Clonal status of actionable driver events and the timing of mutational processes in cancer evolution
.
Sci Transl Med
2015
;
7
:
283ra54
.
2.
Kumar
V
,
Chaudhary
AK
,
Coulter
DW
,
McGuire
T
,
Mahato
RI
.
Impact of miRNA-mRNA profiling and their correlation on medulloblastoma tumorigenesis
.
Molecular Therapy-Nucleic Acids
2018
;
12
:
490
503
.
3.
Muller
S
,
Bley
N
,
Glass
M
,
Busch
B
,
Rousseau
V
,
Misiak
D
, et al
.
IGF2BP1 enhances an aggressive tumor cell phenotype by impairing miRNA-directed downregulation of oncogenic factors
.
Nucleic Acids Res
2018
;
46
:
6285
303
.
4.
Liu
CC
,
Rennie
WA
,
Carmack
CS
,
Kanoria
S
,
Cheng
JJ
,
Lu
J
, et al
.
Effects of genetic variations on microRNA: target interactions
.
Nucleic Acids Res
2014
;
42
:
9543
52
.
5.
Lopes-Ramos
CM
,
Barros
BP
,
Koyama
FC
,
Carpinetti
PA
,
Pezuk
J
,
Doimo
NTS
, et al
.
E2F1 somatic mutation within miRNA target site impairs gene regulation in colorectal cancer
.
PLoS One
2017
;
12
:
e0181153
.
6.
Grigelioniene
G
,
Suzuki
HI
,
Taylan
F
,
Mirzamohammadi
F
,
Borochowitz
ZU
,
Ayturk
UM
, et al
.
Gain-of-function mutation of microRNA-140 in human skeletal dysplasia
.
Nat Med
2019
;
25
:
583
90
.
7.
Bhattacharya
A
,
Cui
Y
.
SomamiR 2.0: a database of cancer somatic mutations altering microRNA-ceRNA interactions
.
Nucleic Acids Res
2016
;
44
:
D1005
10
.
8.
Bhattacharya
A
,
Ziebarth
JD
,
Cui
Y
.
PolymiRTS Database 3.0: linking polymorphisms in microRNAs and their target sites with human diseases and biological pathways
.
Nucleic Acids Res
2014
;
42
:
D86
91
.
9.
Dhawan
A
,
Scott
JG
,
Harris
AL
,
Buffa
FM
.
Pan-cancer characterisation of microRNA across cancer hallmarks reveals microRNA-mediated downregulation of tumour suppressors
.
Nat Commun
2018
;
9
:
5228
.
10.
Plotnikova
O
,
Baranova
A
,
Skoblov
M
.
Comprehensive analysis of human microRNA-mRNA interactome
.
Front Genet
2019
;
10
:
933
.
11.
Lal
A
,
Kim
HH
,
Abdelmohsen
K
,
Kuwano
Y
,
Pullmann
R
Jr
,
Srikantan
S
, et al
.
p16(INK4a) translation suppressed by miR-24
.
PLoS One
2008
;
3
:
e1864
.
12.
Duursma
AM
,
Kedde
M
,
Schrier
M
,
le Sage
C
,
Agami
R
.
miR-148 targets human DNMT3b protein coding region
.
RNA
2008
;
14
:
872
7
.
13.
Tay
Y
,
Zhang
J
,
Thomson
AM
,
Lim
B
,
Rigoutsos
I
.
MicroRNAs to Nanog, Oct4 and Sox2 coding regions modulate embryonic stem cell differentiation
.
Nature
2008
;
455
:
1124
8
.
14.
Shen
WF
,
Hu
YL
,
Uttarwar
L
,
Passegue
E
,
Largman
C
.
MicroRNA-126 regulates HOXA9 by binding to the homeobox
.
Mol Cell Biol
2008
;
28
:
4609
19
.
15.
Tsai
NP
,
Lin
YL
,
Wei
LN
.
MicroRNA mir-346 targets the 5'-untranslated region of receptor-interacting protein 140 (RIP140) mRNA and up-regulates its protein expression
.
Biochem J
2009
;
424
:
411
8
.
16.
Lee
I
,
Ajay
SS
,
Yook
JI
,
Kim
HS
,
Hong
SH
,
Kim
NH
, et al
.
New class of microRNA targets containing simultaneous 5'-UTR and 3'-UTR interaction sites
.
Genome Res
2009
;
19
:
1175
83
.
17.
Maqbool
R
,
Ismail
R
,
Hussain
M
.
Mutations in MicroRNA genes and their binding sites are infrequently associated with human colorectal cancer in the kashmiri population
.
Microrna
2014
;
2
:
219
24
.
18.
Supek
F
,
Minana
B
,
Valcarcel
J
,
Gabaldon
T
,
Lehner
B
.
Synonymous mutations frequently act as driver mutations in human cancers
.
Cell
2014
;
156
:
1324
35
.
19.
Shabalina
SA
,
Spiridonov
NA
,
Kashina
A
.
Sounds of silence: synonymous nucleotides as a key to biological regulation and complexity
.
Nucleic Acids Res
2013
;
41
:
2073
94
.
20.
Sauna
ZE
,
Kimchi-Sarfaty
C
.
Understanding the contribution of synonymous mutations to human disease
.
Nat Rev Genet
2011
;
12
:
683
91
.
21.
Stegeman
S
,
Amankwah
E
,
Klein
K
,
O'Mara
TA
,
Kim
D
,
Lin
HY
, et al
.
A large-scale analysis of genetic variants within putative miRNA binding sites in prostate cancer
.
Cancer Discov
2015
;
5
:
368
79
.
22.
Grey
F
,
Tirabassi
R
,
Meyers
H
,
Wu
GM
,
McWeeney
S
,
Hook
L
, et al
.
A viral microRNA down-regulates multiple cell cycle genes through mRNA 5 ' UTRs
.
PLoS Pathog
2010
;
6
:
e1000967
.
23.
Zhou
HL
,
Rigoutsos
I
.
MiR-103a-3p targets the 5 ' UTR of GPRC5A in pancreatic cells
.
RNA
2014
;
20
:
1431
9
.
24.
Brest
P
,
Lapaquette
P
,
Souidi
M
,
Lebrigand
K
,
Cesaro
A
,
Vouret-Craviari
V
, et al
.
A synonymous variant in IRGM alters a binding site for miR-196 and causes deregulation of IRGM-dependent xenophagy in Crohn's disease
.
Nat Genet
2011
;
43
:
242
U24
.
25.
Shah
SP
,
Morin
RD
,
Khattra
J
,
Prentice
L
,
Pugh
T
,
Burleigh
A
, et al
.
Mutational evolution in a lobular breast tumour profiled at single nucleotide resolution
.
Nature
2009
;
461
:
809
U67
.
26.
Gartner
JJ
,
Parker
SCJ
,
Prickett
TD
,
Dutton-Regester
K
,
Stitzel
ML
,
Lin
JC
, et al
.
Whole-genome sequencing identifies a recurrent functional synonymous mutation in melanoma
.
Proc Nat Acad Sci U S A
2013
;
110
:
13481
6
.
27.
Seiler
M
,
Peng
SY
,
Agrawal
AA
,
Palacino
J
,
Teng
T
,
Zhu
P
, et al
.
Somatic mutational landscape of splicing factor genes and their functional consequences across 33 cancer types
.
Cell Rep
2018
;
23
:
282
96
.
28.
Biankin
AV
,
Waddell
N
,
Kassahn
KS
,
Gingras
MC
,
Muthuswamy
LB
,
Johns
AL
, et al
.
Pancreatic cancer genomes reveal aberrations in axon guidance pathway genes
.
Nature
2012
;
491
:
399
405
.
29.
Moretti
F
,
Thermann
R
,
Hentze
MW
.
Mechanism of translational regulation by miR-2 from sites in the 5 ' untranslated region or the open reading frame
.
RNA
2010
;
16
:
2493
502
.
30.
Agarwal
V
,
Bell
GW
,
Nam
JW
,
Bartel
DP
.
Predicting effective microRNA target sites in mammalian mRNAs
.
Elife
2015
;
4
:
e05005
.
31.
Konig
J
,
Zarnack
K
,
Luscombe
NM
,
Ule
J
.
Protein-RNA interactions: new genomic technologies and perspectives
.
Nat Rev Genet
2012
;
13
:
77
83
.
32.
Li
JH
,
Liu
S
,
Zhou
H
,
Qu
LH
,
Yang
JH
.
starBase v2.0: decoding miRNA-ceRNA, miRNA-ncRNA and protein-RNA interaction networks from large-scale CLIP-Seq data
.
Nucleic Acids Res
2014
;
42
:
D92
7
.
33.
Miranda
KC
,
Huynh
T
,
Tay
Y
,
Ang
YS
,
Tam
WL
,
Thomson
AM
, et al
.
A pattern-based method for the identification of MicroRNA binding sites and their corresponding heteroduplexes
.
Cell
2006
;
126
:
1203
17
.
34.
Kozomara
A
,
miRBase
G-JS
.:
annotating high confidence microRNAs using deep sequencing data
.
Nucleic Acids Res
2014
;
42
:
D68
73
.
35.
Rehmsmeier
M
,
Steffen
P
,
Hochsmann
M
,
Giegerich
R
.
Fast and effective prediction of microRNA/target duplexes
.
RNA
2004
;
10
:
1507
17
.
36.
Xiao
FF
,
Zuo
ZX
,
Cai
GS
,
Kang
SL
,
Gao
XL
,
Li
TB
.
miRecords: an integrated resource for microRNA-target interactions
.
Nucleic Acids Res
2009
;
37
:
D105
D10
.
37.
Chou
CH
,
Chang
NW
,
Shrestha
S
,
Hsu
SD
,
Lin
YL
,
Lee
WH
, et al
.
miRTarBase 2016: updates to the experimentally validated miRNA-target interactions database
.
Nucleic Acids Res
2016
;
44
:
D239
D47
.
38.
Bandyopadhyay
S
,
Mitra
R
.
TargetMiner: microRNA target prediction with systematic identification of tissue-specific negative examples
.
Bioinformatics
2009
;
25
:
2625
31
.
39.
Papadopoulos
GL
,
Reczko
M
,
Simossis
VA
,
Sethupathy
P
,
Hatzigeorgiou
AG
.
The database of experimentally supported targets: a functional update of TarBase
.
Nucleic Acids Res
2009
;
37
:
D155
D8
.
40.
Quinlan
AR
,
Hall
IM
.
BEDTools: a flexible suite of utilities for comparing genomic features
.
Bioinformatics
2010
;
26
:
841
2
.
41.
Kruger
J
,
Rehmsmeier
M
.
RNAhybrid: microRNA target prediction easy, fast and flexible
.
Nucleic Acids Res
2006
;
34
:
W451
W4
.
42.
Lopes-Ramos
CM
,
Chen
CY
,
Kuijjer
ML
,
Paulson
JN
,
Sonawane
AR
,
Fagny
M
, et al
.
Sex differences in gene expression and regulatory networks across 29 human tissues
.
Cell Rep
2020
;
31
:
107795
.
43.
Kircher
M
,
Witten
DM
,
Jain
P
,
O'Roak
BJ
,
Cooper
GM
,
Shendure
J
.
A general framework for estimating the relative pathogenicity of human genetic variants
.
Nat Genet
2014
;
46
:
310
5
.
44.
Watson
IR
,
Takahashi
K
,
Futreal
PA
,
Chin
L
.
Emerging patterns of somatic mutations in cancer
.
Nat Rev Genet
2013
;
14
:
703
18
.
45.
Babur
O
,
Gonen
M
,
Aksoy
BA
,
Schultz
N
,
Ciriello
G
,
Sander
C
, et al
.
Systematic identification of cancer driving signaling pathways based on mutual exclusivity of genomic alterations
.
Genome Biol
2015
;
16
:
45
.
46.
Thomas
RK
,
Baker
AC
,
Debiasi
RM
,
Winckler
W
,
Laframboise
T
,
Lin
WM
, et al
.
High-throughput oncogene mutation profiling in human cancer
.
Nat Genet
2007
;
39
:
347
51
.
47.
Liu
G
,
Zhang
R
,
Xu
J
,
Wu
CI
,
Lu
X
.
Functional conservation of both CDS- and 3'-UTR-located microRNA binding sites between species
.
Mol Biol Evol
2015
;
32
:
623
8
.
48.
Hausser
J
,
Syed
AP
,
Bilen
B
,
Zavolan
M
.
Analysis of CDS-located miRNA target sites suggests that they can effectively inhibit translation
.
Genome Res
2013
;
23
:
604
15
.
49.
Hafner
M
,
Landthaler
M
,
Burger
L
,
Khorshid
M
,
Hausser
J
,
Berninger
P
, et al
.
Transcriptome-wide identification of RNA-binding protein and microRNA target sites by PAR-CLIP
.
Cell
2010
;
141
:
129
41
.
50.
Fang
Z
,
Rajewsky
N
.
The impact of miRNA target sites in coding sequences and in 3'UTRs
.
PLoS One
2011
;
6
:
e18067
.
51.
Wang
T
,
Wang
G
,
Hao
D
,
Liu
X
,
Wang
D
,
Ning
N
, et al
.
Aberrant regulation of the LIN28A/LIN28B and let-7 loop in human malignant tumors and its effects on the hallmarks of cancer
.
Mol Cancer
2015
;
14
:
125
.
52.
Takamizawa
J
,
Konishi
H
,
Yanagisawa
K
,
Tomida
S
,
Osada
H
,
Endoh
H
, et al
.
Reduced expression of the let-7 microRNAs in human lung cancers in association with shortened postoperative survival
.
Cancer Res
2004
;
64
:
3753
6
.
53.
Ji
J
,
Zhao
L
,
Budhu
A
,
Forgues
M
,
Jia
HL
,
Qin
LX
, et al
.
Let-7g targets collagen type I alpha2 and inhibits cell migration in hepatocellular carcinoma
.
J Hepatol
2010
;
52
:
690
7
.
54.
Hu
X
,
Guo
J
,
Zheng
L
,
Li
C
,
Zheng
TM
,
Tanyi
JL
, et al
.
The heterochronic microRNA let-7 inhibits cell motility by regulating the genes in the actin cytoskeleton pathway in breast cancer
.
Mol Cancer Res
2013
;
11
:
240
50
.
55.
Johnson
RW
,
Finger
EC
,
Olcina
MM
,
Vilalta
M
,
Aguilera
T
,
Miao
Y
, et al
.
Induction of LIFR confers a dormancy phenotype in breast cancer cells disseminated to the bone marrow
.
Nat Cell Biol
2016
;
18
:
1078
89
.
56.
Ma
D
,
Jing
X
,
Shen
B
,
Liu
X
,
Cheng
X
,
Wang
B
, et al
.
Leukemia inhibitory factor receptor negatively regulates the metastasis of pancreatic cancer cells in vitro and in vivo
.
Oncol Rep
2016
;
36
:
827
36
.
57.
Hua
WF
,
Zhong
Q
,
Xia
TL
,
Chen
Q
,
Zhang
MY
,
Zhou
AJ
, et al
.
RBM24 suppresses cancer progression by upregulating miR-25 to target MALAT1 in nasopharyngeal carcinoma
.
Cell Death Dis
2016
;
7
:
e2352
.
58.
Jiang
Y
,
Zhang
M
,
Qian
Y
,
Xu
E
,
Zhang
J
,
Chen
X
.
Rbm24, an RNA-binding protein and a target of p53, regulates p21 expression via mRNA stability
.
J Biol Chem
2014
;
289
:
3164
75
.
59.
Sun
C
,
Ma
P
,
Wang
Y
,
Liu
W
,
Chen
Q
,
Pan
Y
, et al
.
KLF15 inhibits cell proliferation in gastric cancer cells via up-regulating CDKN1A/p21 and CDKN1C/p57 expression
.
Dig Dis Sci
2017
;
62
:
1518
26
.
60.
Ko
HK
,
Akakura
S
,
Peresie
J
,
Goodrich
DW
,
Foster
BA
,
Gelman
IH
.
A transgenic mouse model for early prostate metastasis to lymph nodes
.
Cancer Res
2014
;
74
:
945
53
.
61.
Peixoto
P
,
Blomme
A
,
Costanza
B
,
Ronca
R
,
Rezzola
S
,
Palacios
AP
, et al
.
HDAC7 inhibition resets STAT3 tumorigenic activity in human glioblastoma independently of EGFR and PTEN: new opportunities for selected targeted therapies
.
Oncogene
2016
;
35
:
4481
94
.
62.
Wu
W
,
Zhou
X
,
Yu
T
,
Bao
Z
,
Zhi
T
,
Jiang
K
, et al
.
The malignancy of miR-18a in human glioblastoma via directly targeting CBX7
.
Am J Cancer Res
2017
;
7
:
64
76
.
63.
Ni
S
,
Wang
H
,
Zhu
X
,
Wan
C
,
Xu
J
,
Lu
C
, et al
.
CBX7 suppresses cell proliferation, migration, and invasion through the inhibition of PTEN/Akt signaling in pancreatic cancer
.
Oncotarget
2017
;
8
:
8010
21
.
64.
Zhang
Y
,
Kang
R
,
Liu
W
,
Yang
Y
,
Ding
R
,
Huang
Q
, et al
.
Identification and analysis of P53-mediated competing endogenous RNA network in human hepatocellular carcinoma
.
Int J Biol Sci
2017
;
13
:
1213
21
.
65.
Wang
L
,
Zhu
MJ
,
Ren
AM
,
Wu
HF
,
Han
WM
,
Tan
RY
, et al
.
A ten-microRNA signature identified from a genome-wide microRNA expression profiling in human epithelial ovarian cancer
.
PLoS One
2014
;
9
:
e96472
.
66.
Jiang
Y
,
Gao
H
,
Liu
M
,
Mao
Q
.
Sorting and biological characteristics analysis for side population cells in human primary hepatocellular carcinoma
.
Am J Cancer Res
2016
;
6
:
1890
905
.
67.
Alderman
C
,
Yang
Y
.
The anti-melanoma activity and oncogenic targets of hsa-miR-15a-5p
.
RNA Dis
2016
;
3
:
e1450
.
68.
Ma
Y
,
Pan
X
,
Xu
P
,
Mi
Y
,
Wang
W
,
Wu
X
, et al
.
Plasma microRNA alterations between EGFR-activating mutational NSCLC patients with and without primary resistance to TKI
.
Oncotarget
2017
;
8
:
88529
36
.
69.
Yeang
CH
,
McCormick
F
,
Levine
A
.
Combinatorial patterns of somatic gene mutations in cancer
.
FASEB J
2008
;
22
:
2605
22
.
70.
Xing
M
,
Cohen
Y
,
Mambo
E
,
Tallini
G
,
Udelsman
R
,
Ladenson
PW
, et al
.
Early occurrence of RASSF1A hypermethylation and its mutual exclusion with BRAF mutation in thyroid tumorigenesis
.
Cancer Res
2004
;
64
:
1664
8
.
71.
Carbuccia
N
,
Trouplin
V
,
Gelsi-Boyer
V
,
Murati
A
,
Rocquain
J
,
Adelaide
J
, et al
.
Mutual exclusion of ASXL1 and NPM1 mutations in a series of acute myeloid leukemias
.
Leukemia
2010
;
24
:
469
73
.
72.
Ro
S
,
Park
C
,
Young
D
,
Sanders
KM
,
Yan
W
.
Tissue-dependent paired expression of miRNAs
.
Nucleic Acids Res
2007
;
35
:
5944
53
.
73.
Telonis
AG
,
Magee
R
,
Loher
P
,
Chervoneva
I
,
Londin
E
,
Rigoutsos
I
.
Knowledge about the presence or absence of miRNA isoforms (isomiRs) can successfully discriminate amongst 32 TCGA cancer types
.
Nucleic Acids Res
2017
;
45
:
2973
85
.
74.
Kishore
S
,
Jaskiewicz
L
,
Burger
L
,
Hausser
J
,
Khorshid
M
,
Zavolan
M
.
A quantitative analysis of CLIP methods for identifying binding sites of RNA-binding proteins
.
Nat Methods
2011
;
8
:
559
64
.
This open access article is distributed under the Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International (CC BY-NC-ND 4.0) license.