Widespread inference of weighted microRNA-mediated gene regulation in cancer transcriptome analysis

Abstract

MicroRNAs (miRNAs) comprise a gene-regulatory network through sequence complementarity with target mRNAs. Previous studies have shown that mammalian miRNAs decrease many target mRNA levels and reduce protein production predominantly by target mRNA destabilization. However, it has not yet been fully assessed whether this scheme is widely applicable to more realistic conditions with multiple miRNA fluctuations. By combining two analytical frameworks for detecting the enrichment of gene sets, Gene Set Enrichment Analysis (GSEA) and Functional Assignment of miRNAs via Enrichment (FAME), we developed GSEA-FAME analysis (GFA), which enables the prediction of miRNA activities from mRNA expression data using rank-based enrichment analysis and weighted evaluation of miRNA-mRNA interactions. This cooperative approach delineated a better widespread correlation between miRNA expression levels and predicted miRNA activities in cancer transcriptomes, thereby providing proof-of-concept of the mRNA-destabilization scenario. In an integrative analysis of The Cancer Genome Atlas (TCGA) multidimensional data including profiles of both mRNA and miRNA, we also showed that GFA-based inference of miRNA activity could be used for the selection of prognostic miRNAs in the development of cancer survival prediction models. This approach proposes a next-generation strategy for the interpretation of miRNA function and identification of target miRNAs as biomarkers and therapeutic targets.

Figure 1.

Outline of GFA. In step 1, GSEA is performed using C3 miRNA target gene sets (C3MIR), which contain potential target genes sharing a 3′UTR miRNA-binding motif for each miRNA, to assess whether each miRNA target gene set is enriched in group A or in group B. In step 2, leading-edge subsets, part of the members of each miRNA target gene set, which accounts for the enrichment of corresponding gene sets in group A or B in GSEA analysis, are collected for each gene set and assembled for groups A and B, respectively, to make the collection of overall miRNA target genes enriched in group A and B. As an option, similar procedures to steps 1 and 2 are performed using C3 transcription factor target gene sets (C3TFT) that contain genes sharing a transcription factor-binding site and the C3TFT target gene collection is subtracted from the C3MIR target gene collection subjected to FAME (step 2’). In step 3, FAME is applied to the C3MIR target gene collection in step 2, resulting in a ranked list of each miRNA activity and corresponding target genes accounting for this activity.

Figure 2.

Comparison of GSEA and GFA for miR-1 and miR-124 transfection data. GSEA (A) and GFA (B) were performed for the microarray data of HeLa cells transfected with miR-1, miR-124, mutant miR-124 (124mut5-6 and 124mut9-10) and chimeric miRNAs (chimiR-1/124 and chimiR-124/1) for 12 or 24 h. In GFA, we ran FAME for the two collection of leading-edge subsets, which were down-regulated (left panel) and up-regulated (right panel) in the cells transfected with miRNA relative to controls, in an enrichment direction of over-representation (enrichment, left panel) and under-representation (depletion, right panel), respectively. The distribution of GSEA normalized enrichment scores (NES) and GFA ranking [−log10(P-value)] are shown. GSEA NES and GFA ranking for miR-1 and miR-124 are indicated by red diamonds and blue circles, respectively.

Figure 3.

Widespread correlation of miRNA expression levels and miRNA activities assessed by GFA: DLBCL study. (A) Classification of DLBCL cases into ABC and GC subtypes. (B) Correlation between miRNA expression ranks and miRNA activity ranks analysed by GSEA, FAME [after the extraction of differentially expressed genes by the t-test (P < 0.05)] and GFA (using C3MIR and C3TFT gene sets) for the top 45 miRNAs with differential expression between ABC and GC subtypes. Differential miRNA expression was evaluated by −log10(P-value). (C) Spearman’s rank correlation coefficients for the correlation between miRNA expression ranks and miRNA activity ranks in (B). (D) RRHO analysis for the correlation between miRNA expression ranks and miRNA activity ranks in (B). Maximums of the Benjamini–Yekutieli-corrected RRHO map (left) and the representative RRHO heatmap (right) are shown. (E) Overlap between target genes of GFA-supported differentially expressed miRNAs and gene subsets of C4 cancer module (C4-CM) gene sets analysed by GSEA. In addition to GFA, GSEA was performed among ABC and GC subtypes using C4-CM gene sets. Red boxes represent overlap between target genes of GFA-supported differentially expressed miRNAs and leading-edge subsets of C4-CM gene sets enriched in either DLBCL subtypes.

Figure 4.

Availability of GFA for the selection of prognostic miRNAs in the development of cancer survival prediction models: TCGA glioblastoma study. (A) Analytical outline of survival prediction using TCGA glioblastoma data sets. (B) Examples of Kaplan–Meier plots representing survival probabilities according to low or high levels of risk scores developed by (i) the ‘expression level-based’ strategy and (ii) the ‘expression level/GFA-based’ strategy in a training set and a test set. (C) Likelihood ratio of survival prediction models developed by (i) the ‘expression level-based’ strategy and (ii) the ‘expression level/GFA-based’ strategy in training sets and test sets. Positive numbers indicate that the model fits the data better. Results with five randomizations are shown. (D) Time-dependent prediction errors of survival prediction models developed by (i) the ‘expression level-based’ strategy and (ii) the ‘expression level/GFA-based’ strategy in training sets and test sets. Average results with five randomizations are shown.