Integrative network analysis identifies cell-specific trans regulators of m6A

Abstract

N6-methyladenosine (m6A) is a reversible and dynamic RNA modification in eukaryotes. However, how cells establish cell-specific m6A methylomes is still poorly understood. Here, we developed a computational framework to systematically identify cell-specific trans regulators of m6A through integrating gene expressions, binding targets and binding motifs of large number of RNA binding proteins (RBPs) with a co-methylation network constructed using large-scale m6A methylomes across diverse cell states. We applied the framework and successfully identified 32 high-confidence m6A regulators that modulated the variable m6A sites away from stop codons in a cell-specific manner. To validate them, we knocked down three regulators respectively and found two of them (TRA2A and CAPRIN1) selectively promoted the methylations of the m6A sites co-localized with their binding targets on RNAs through physical interactions with the m6A writers. Knockdown of TRA2A increased the stabilities of the RNAs with TRA2A bound near the m6A sites and decreased the viability of cells. The successful identification of m6A regulators demonstrates a powerful and widely applicable strategy to elucidate the cell-specific m6A regulators. Additionally, our discovery of pervasive trans-acting regulating of m6A provides novel insights into the mechanisms by which spatial and temporal dynamics of m6A methylomes are established.

Figure 1.

Schematic flow chart demonstrating the computational framework to identify trans regulators of m⁶A.

Figure 2.

Analyses of m⁶A methylomes of multiple cell lines. (A) Normalized distributions of m⁶A peaks across 5′ UTR, CDS and 3′ UTR for representative cell lines. (B) The unsupervised hierarchical clustering and heatmap of the m⁶A ratios for the m⁶A peaks with the largest CVs across all cell lines. The technical information is indicated above the heatmap. (C) Box plot representing the CVs of m⁶A ratios for the peaks located at different regions of mRNAs. (D) Normalized distributions of variable and stable m⁶A peaks across 5′ UTR, CDS and 3′ UTR. (E) Densities of logarithm transformed lengths of the internal exons with variable m⁶A peaks and stable m⁶A peaks. The P-value of Wilcoxon test is indicated. (F) Tracks showing the read coverage of the IPs, inputs and the merged m⁶A peaks of the representative cell lines as well as the HEK293 m⁶A sites from miCLIP-seq data on ILF2. The tracks are shown for optimal viewing. The variable and stable m6A peaks are highlighted, respectively. (G) Box plots representing the logarithm transformed TPMs of variable and stable m⁶A peaks. (H), Box plot representing the maximum m⁶A ratios across all cell lines of variable and stable m⁶A peaks. (I) Bar plot showing the percentages of variable peaks and stable peaks that overlap with m⁶A sites obtained from miCLIP-seq. ‘n.s.’ denotes non-significant.

Figure 3.

Classification and analyses of co-methylated m⁶A modules. (A) Classification of co-methylated m⁶A modules through dynamical cutting of the clustering dendrogram of all variable m⁶A peaks. (B) Heatmap representing the m⁶A indexes of all the 41 co-methylation modules across all the cell lines. (C) Normalized distributions of m⁶A peaks in different combined modules across 5′ UTR, CDS and 3′ UTR. (D) Densities of logarithm transformed lengths of the internal exons with m⁶A peaks in different combined co-methylation modules.

Figure 4.

Discovery of a co-methylation module specifically methylated in cancer cell lines. (A) Heatmaps representing the Z-scores of m⁶A ratios (upper panel) and gene expressions (low panel) of the peaks and corresponding genes across all cell lines. The types of cell lines are indicated at the top of the upper panel. (B) Box plot representing the gene expression indexes of the genes corresponding to cancer-specific module for cancer and normal samples of 14 cancer types in TCGA. (C) Cox correlations between the gene expression indexes of the genes corresponding to cancer-specific module and the survival of cancer patients of 14 cancer types in TCGA. OS: overall survival; DFS: disease-free survival; HR: hazard ratio. (D) Tracks representing the gene expression indexes of the genes corresponding to the cancer-specific module and genetic alteration spectrum of the key markers as well as clinical phenotypes of the breast cancer patients from TCGA. The patient samples are sorted according to the gene expression indexes of the cancer-specific module.

Figure 5.

Systematic identification of m⁶A regulators. (A) Q-Q plot comparing the distributions of expected P-values and observed P-values of the correlations between the gene expressions of RBPs and m⁶A indexes of the co-methylation modules. (B) Heatmaps representing the m⁶A ratios of the m⁶A peaks within the module M15 (upper panel) and the gene expressions of the RBPs that significantly correlated with the m⁶A indexes of M15 (lower panel). The cell lines are sorted according to the m⁶A indexes of M15. (C andD) Barplot representing the percentages of the pairs of RBPs and modules that enriched for CLIP-seq binding sites (C) or motifs (D) of the RBPs out of the pairs that showing significant and non-significant (top 1000 least significant) correlations between gene expressions of the RBPs and the m⁶A indexes of the modules. (E) Venn diagram demonstrating the identification of 32 high-confidence m⁶A regulators. (F) Scatter plot showing the correlation between TRA2A gene expression and m⁶A indexes of module M15 across all cell lines. (G) Representative motifs enriched in module M15 and TRA2A CLIP-seq targets.

Figure 6.

Experimental validation of selected m⁶A regulators. (A) Plot of cumulative fraction of log2 fold change of m⁶A ratios upon TRA2A knockdown for the m⁶A peaks overlap or non-overlap with TRA2A CLIP-seq targets. P-value of two-tailed Wilcoxon test is indicated. (B) Plot of cumulative fraction of log2 fold change of m⁶A ratios upon TRA2A knockdown for the m⁶A peaks within or not within the co-methylation modules correlated with TRA2A. P value of two-tailed Wilcoxon test is indicated. (C) Tracks displaying the read coverage of IPs normalized by inputs as well as the miCLIP-seq m⁶A sites and TRA2A CLIP-seq peaks in HepG2 cells on the long non-coding gene MALAT1. The m⁶A peak with down-regulated m⁶A ratio in shTRA2A is highlighted. The dashed lines indicate the summits of the peaks. (D) Western blots showing the interaction between TRA2A and METTL3 with and without RNase treatment respectively. * indicates a non-specific band (see Supplementary Figure S5B). (E) Plot of cumulative fraction of log2 fold change of m⁶A ratios upon CAPRIN1 knockdown comparing the m⁶A peaks in the correlated modules and overlap with CAPRIN1 CLIP-seq targets versus the peaks not within the correlated module or overlap with CAPRIN1 CLIP-seq targets. P-value of two-tailed Wilcoxon test is indicated. (F) Western blots showing the interactions of CAPRIN1 with METTL3 and METTL14 with and without RNase treatment respectively. * indicates a non-specific band (see Supplementary Figure S5B). (G) Plot of cumulative fraction of log2 fold change of m⁶A ratios upon MOV10 knockdown comparing the m⁶A peaks in the correlated modules and overlap with MOV10 CLIP-seq targets versus the peaks not within the correlated module or overlap with MOV10 CLIP-seq targets. P-value of two-tailed Wilcoxon test is indicated.