M6A2Target: a comprehensive database for targets of m⁶A writers, erasers and readers

Abstract

N⁶-methyladenosine (m⁶A) is the most abundant posttranscriptional modification in mammalian mRNA molecules and has a crucial function in the regulation of many fundamental biological processes. The m⁶A modification is a dynamic and reversible process regulated by a series of writers, erasers and readers (WERs). Different WERs might have different functions, and even the same WER might function differently in different conditions, which are mostly due to different downstream genes being targeted by the WERs. Therefore, identification of the targets of WERs is particularly important for elucidating this dynamic modification. However, there is still no public repository to host the known targets of WERs. Therefore, we developed the m⁶A WER target gene database (m6A2Target) to provide a comprehensive resource of the targets of m⁶A WERs. M6A2Target provides a user-friendly interface to present WER targets in two different modules: ‘Validated Targets’, referred to as WER targets identified from low-throughput studies, and ‘Potential Targets’, including WER targets analyzed from high-throughput studies. Compared to other existing m⁶A-associated databases, m6A2Target is the first specific resource for m⁶A WER target genes. M6A2Target is freely accessible at http://m6a2target.canceromics.org.

Introduction

RNA modifications were first described approximately 60 years ago [1]. To date, more than 100 kinds of RNA modifications have been discovered and determined to be involved in various kinds of cellular biological processes [2]. Among these hundreds of RNA modifications, one of the most abundant and widely distributed internal modifications in mRNA is N⁶-methyladenosine (m⁶A) [3]. As the most abundant epigenetic modification, m⁶A has been shown to be associated with the regulation of mRNA export [4], stability [5], splicing [6] and translation [7], while malfunction of m⁶A will result in changes in multiple cell functions, leading to disease occurrence. For example, FTO, a well-known eraser of m⁶A, has long been suggested to be associated with obesity in several population studies [8, 9] and plays a vital role in leukemogenesis [10]. Downregulation of METTL14, a writer for m⁶A modification, leads to poor prognosis in hepatocellular carcinoma and is significantly associated with tumor metastasis [11]. In addition, m⁶A modification is also considered to be related to type 2 diabetes mellitus [12], infertility [13] and growth retardation [14].

M⁶A modification is reversible and modulated by ‘writers’ (such as METTL3, METTL14 and WTAP) [15, 16], ‘erasers’ (such as ALKBH5 and FTO) [17, 18] and ‘readers’ (such as YTHDF1, YTHDF2 and YTHDF3) [19–22]. Although these ‘writers’, ‘erasers’ and ‘readers’ (WERs) are functional through m⁶A modifications, different WERs may execute entirely different functions [17–25], and even the same WER may have distinct functions under different conditions [26–28]. These differences are mostly due to different sets of downstream genes targeted by different WERs. Different WERs may target different sets of genes under the same conditions, and the same WER may also target different sets of genes under different conditions. Therefore, identification of WER targets is particularly important for further characterization of this dynamic modification.

Growing studies have focused on the identification of the targets of m⁶A WERs. Various kinds of next-generation sequencing technologies [29–33], mass spectrometry [34] and some binding experiments, such as coimmunoprecipitation [35], have been applied to infer the relationship between m⁶A WERs and their targets. To provide a comprehensive and convenient platform for researchers to obtain related information about WERs and their targets, we established the m⁶A WER target gene database (m6A2Target), which is freely accessible at http://m6a2target.canceromics.org. The pipeline for the database construction is described in Figure 1. M6A2Target integrated basic information of WERs and their targets, including gene annotation information from authority sites, such as NCBI, Ensemble and UniProt. In addition, m6A2Target also provided information on the associations between WERs and their targets, such as protein–DNA, protein–RNA or protein–protein interactions, and changes in gene expression level, m⁶A modification level or translation efficacy and alternative splicing events followed by perturbation of WERs. Growing studies have revealed vital roles of m⁶A in some fundamental cellular functions, such as messenger RNA stability [36], mRNA splicing [28] and RNA translation efficiency [19], and in various kinds of diseases, including cancers. Therefore, we believe that such a database collecting the targets of m⁶A WER in a cell-specific manner will notably facilitate further downstream functional and mechanistic studies in the m⁶A research field.

Materials and methods

Data sources

To establish m6A2Target, we first collected all reported m⁶A WERs from the literature in humans and mice. For human species, we obtained 10 writers, 2 erasers and 10 readers, and for mice, we obtained 6 writers, 2 erasers and 6 readers. We collected all the studies and datasets related to m⁶A WERs from NCBI PubMed, Gene Expression Omnibus (GEO) [37] and the Sequence Read Archive (SRA) [38] updated to April 2019. The quality of the raw data was checked using fastqc (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/), and the raw data that did not pass the criteria of fastqc were not used. In total, 77 datasets for 80 kinds of cell lines or mouse models were obtained. Genome Reference Consortium Human Build 38 (hg38) and Genome Reference Consortium Mouse Build 38 (mm10) were used for genome reference for humans and mice, respectively. To conduct downstream biological applications using our database, we collected gene expression data, gene mutation data and clinical information of 33 cancer types from The Cancer Genome Atlas (TCGA) [39].

Identification of validated WER targets from low-throughput studies

To obtain validated WER targets, we collected all m⁶A-related articles from the publicly available literature (Supplementary Table 1). We then manually curated WER targets from low-throughput experiments, such as immunoprecipitation assays, luciferase reporter assays and perturbation experiments. A gene was considered a validated WER target if it was shown to interact with any WER according to such experiments as immunoprecipitation and luciferase reporter, or its expression level, methylation level, translation efficiency, stability or other biological manifestations were altered when any WER was perturbed (such as knocked down, knocked out, mutated and overexpressed).

Identification of potential targets with binding evidence

To identify potential targets with binding evidence, we collected high-throughput sequencing data about protein–DNA interactions and protein–RNA interactions and mass spectrometry data about protein–protein interactions for all WERs (Supplementary Table 2).

To identify potential targets with protein–protein interactions with WERs, we summarized the mass spectrometry results from the supplementary data of published papers.

For identification of potential targets having protein–DNA interactions with WERs, we collected raw ChIP-seq data for all WERs from GEO [37] and SRA [38]. These raw ChIP-seq data were cleaned using fastp [40], mapped to the reference genome using bowtie2 [41], then called peaks using MACS2 [42]. The peaks were annotated to genes using HOMER (http://homer.ucsd.edu/homer/ngs/annotation.html) [43]. All genes with ChIP-seq peaks were considered potential WER targets with protein–DNA binding evidence.

For identification of potential targets having protein–RNA interactions with WERs, we collected raw sequencing data, such as RIP-seq data and various types of CLIP-seq data, for all WERs from GEO [37] and SRA [38]. For iCLIP-seq data, HITS-CLIP-seq data and eCLIP-seq data, the CLIP Tool Kit [44] was used to identify the potential binding sites, and HOMER (http://homer.ucsd.edu/homer/ngs/annotation.html) [43] was used for annotation. For PAR-CLIP-seq, after cleaning with fastp [40], bowtie [45] was first used for read alignment, PARalyzer v1.5 [46] was subsequently used to identify potential binding sites and HOMER (http://homer.ucsd.edu/homer/ngs/annotation.html) [43] was used for annotation. For RIP-seq data, with fastp [40] for read cleaning, RSEM [47] was used for RNA level quantification followed by FPKM (Fragments Per Kilobase Million) normalization, and genes with fold change between IP and INPUT greater than 2 (IP/INPUT >2) were considered potential binding targets.

Identification of potential targets inferred from WER perturbation

The potential targets inferred from WER perturbation have four levels: RNA, m⁶A, translation and alternative splicing. To identify potential targets at the RNA level, we collected RNA-seq data from studies with wild-type and perturbed WERs of interest. The perturbation condition could be knocked down, knocked out, mutated or overexpressed. RSEM [47] was used to generate gene counts, and DESeq2 [48] was used to obtain differentially expressed genes (P-value <0.05, fold change >1.5). The differentially expressed genes were considered potential WER targets whose RNA level could be altered by perturbation of the WER of interest.

To identify potential targets at the m⁶A level, we collected m⁶A-specific methylated RNA immunoprecipitation sequencing (MeRIP-seq) data from studies with wild-type and perturbed WERs of interest. The following pipeline was used for MeRIP-seq data analysis. First, using fastp [40] for read cleaning, STAR [49] was used for read alignment. Second, MACS2 [42] and MeTPeak [50] were used for m⁶A peak calling. Third, the relative m⁶A level was obtained by calculating the ratio between the IP RPKM value and the input RPKM value for each peak. The differential m⁶A peaks (P-value <0.05, fold change >1.5) between the two conditions were considered potential WER targets whose m⁶A level could be altered by perturbation of the WER of interest. The significance of the differential m⁶A methylation was measured using the model from Audic et al. [51].

To identify potential targets at the translation level, we collected ribosome profiling data from studies with wild-type and perturbed WERs of interest. The raw ribosome profiling data were first cleaned using fastp [40] and then filtered for mitochondrial DNA and ribosome RNA using Bowtie2 [41]. The remaining reads were mapped to the genome using STAR [49]. RSEM [47] was used to obtain the read count for each gene followed by FPKM normalization. The translation efficiency difference between the two conditions was calculated as log2 (FPKM of the treatment group/FPKM of the control group). Genes with different translation efficiencies (fold change >1.5) were considered potential WER targets whose translation efficiency could be altered by the perturbation of the WER of interest.

To identify potential targets with differential alternative splicing events, we used the cleaned RNA-seq data mentioned above with STAR [49] for read alignment and then used rMATS [52] for alternative splicing event detection. The differential alternative splicing events (P-value <0.05) were considered to be potential WER targets whose alternative splicing level could be altered by perturbation of the WER of interest.

Web interface implementation

All the analysis results and sample information were managed using MySQL tables (https://www.mysql.com/cn/). The web interfaces were implemented using Vue.js, Element UI and JavaScript. Furthermore, we inlayed Genomic Browser and allowed the genomic coordinates of the regions of interest as an input to search for the corresponding target location information of the user-defined genes.

Application of m6A2Target in cancer study

We collected and downloaded all mutation data, gene expression data and clinical data of 33 cancer types from TCGA (https://portal.gdc.cancer.gov/) [39]. The mutation rates of each WER in all cancer types and the oncoplot were obtained by Maftools [53]. Limma [54] was used to obtain the P-value and log2 fold change of differential KIAA1429 expression between tumor and normal samples across all cancer types. The hazard ratio of overall survival probability between the high KIAA1429 expression group and the low KIAA1429 expression group was analyzed by Cox proportional hazards regression analysis. KOBAS 3.0 [55–57] was used for gene set enrichment analysis, and the top 10 enrichment pathways are shown. Genes belonging to the PI3K-AKT signaling pathway were obtained from Kyoto Encyclopedia of Genes and Genomes (KEGG) [58]. IGV [59] was used for the visualization of the genomic distribution of two randomly selected genes. All analyses were performed using R 3.5.2 (https://www.r-project.org/).

Result

Overview of m6A2Target

For validated targets, m6A2target contains 668 targets for 22 WERs in humans and 366 targets for 14 WERs in mice. For potential targets with binding evidence, m6A2Target contains 112,592 WER-target pairs in different cell lines from humans and 12,180 WER-target pairs in cell lines from mice. For potential targets inferred from WER perturbation, in human cell lines, there were 98,046 WER-target pairs with RNA level changes, 54,146 WER-targets with m⁶A level changes, 17,878 WER-targets with translation efficiency changes and 61,370 WER-target pairs with differential alternative splicing events; in mice, there were 35,010 WER-target pairs with RNA level changes, 27,438 WER targets with m⁶A level changes and 36,214 WER-target pairs with differential alternative splicing events.

Validated WER targets from low-throughput studies

The number of m⁶A-related studies has been increasing annually since 2012, indicating that the m⁶A has become an increasingly popular research topic in the past several years (Figure 2A). According to the validated targets recorded in m6A2Target, the m⁶A WERs shared many targets and formed a complex network, indicating that these WERs may cooperate with each other to exert biological functions through m⁶A (Figures 2B and S1A). The validated targets could be classified into four major categories according to the alteration in cellular functions followed by perturbation of WER, including gene expression level, m⁶A level, translation efficiency and alternative splicing, and other functions, such as RNA export, RNA localization and RNA editing, were also reported (Figures 2C and S1B). Most of the validated targets were protein-coding genes, and microRNAs were the second most common gene type (Figures 2D and S1C).

Potential targets with binding evidence revealed by high-throughput studies

The potential targets with binding evidence were primarily obtained from the analysis of CLIP-seq, RIP-seq, ChIP-seq and mass spectrum data of all available WERs (Figure 3A). We classified these targets into three categories: the targets inferred from RIP-seq and CLIP-seq were regarded as ‘Protein–RNA’ associations, targets inferred from ChIP-seq were regarded as ‘Protein–DNA’ associations and targets summarized from WER mass spectrometry results were regarded as ‘Protein–Protein’ associations (Figure 3B). Among these potential targets with binding evidence per WER, the targets with m⁶A peaks were considered m⁶A-dependent targets, while others were considered to be m⁶A-independent. Approximately half of the targets did not contain m⁶A peaks, indicating that some m⁶A WER targets might be functional independent of m⁶A (Figures 3C and S2A). Most of these potential targets were protein-coding RNAs, and the second were long noncoding RNAs (lncRNAs) (Figures 3D and S2B).

We next performed several analyses to confirm the reliability of our database. For instance, YTHDF1 is an RNA-binding protein that is known to be functional through m⁶A. We found that the distribution of the ‘Protein–RNA’ targets of YTHDF1 along the genes was similar to the general m⁶A distribution, a peak in the stop-codon region, indicating their close association with m⁶A methylation (Figure 3E) FTO, a known m⁶A eraser. We calculated the distance between the ‘Protein–RNA’ targets of FTO and FTO-associated m⁶A peaks and found that the smaller the distance between FTO targets and m⁶A peaks is, the higher the density is, further indicating that FTO binding targets may be strongly influenced by m⁶A (Figure 3F). In addition, we chose the mass spectrometry results of an m⁶A writer, WTAP, for further analysis. These ‘Protein–Protein’ targets of WTAP can be enriched into dozens of signaling pathways, such as ribosome, spliceosome, mRNA surveillance pathway and RNA transport, illustrating the vital roles of m⁶A WERs in various fundamental cellular functions (Figure 3G).

Potential targets inferred from m⁶A WER perturbation followed by high-throughput sequencing

The other kinds of potential targets were inferred from m⁶A WER perturbation followed by high-throughput sequencing, such as RNA-Seq, m⁶A-Seq and ribosome profiling (Figures 4A and S3A). For example, after METTL3 was perturbed, 26,945 potential targets were altered at the RNA level, 14,881 potential targets were altered at the m⁶A level and 12,507 potential targets were altered at the translation level; the targets of the three levels exhibited considerable overlap, indicating the complex roles of METTL3 (Figure 4B). All the WER targets constituted a complicated network, suggesting that these WERs function in a cooperative manner (Figures 4C and S3B). These targets can also be grouped into various kinds of gene types, demonstrating their vital functions with WERs through different kinds of mechanisms (Figures 4D and S3C).

Web interface and usage

M6A2Target is publicly available at http://m6a2target.canceromics.org. Users can browse, search and download all related data of m⁶A WERs and their targets through our web interface. M6A2Target is mainly composed of five modules, including ‘Validated Targets’, ‘Potential Targets’ and ‘Genomic Browser’, ‘Search’ and ‘Download.’

In the ‘Validated Targets’ module (Figure 5A), users can browse validated targets in specific species, specific cell lines or specific types of WER. Moreover, users can browse the validated targets for a specific WER by clicking the name of the WER of interest. By clicking the ‘Detail’ button for a specific target in the interactive table, users can obtain detailed information of the selected WER-target pair, including detailed information of the WER and target gene, as well as detailed information for the interaction between the WER and target (Figure S4A).

The ‘Potential Targets’ module includes two submodules, named ‘Binding’ and ‘Perturbation.’ Similar to the ‘Validated Targets’ module, users can browse all potential targets with binding evidence in the ‘Binding’ module (Figure 5B) and browse all the potential targets with perturbation evidence in the ‘Perturbation’ module (Figure S4B). Specifically, the page below in this module will display some statistical information about targets of this user selected WER, including target numbers per RNA type (Figure S4C), target numbers per cell type (Figure S4D) and the overview of these targets (Figure S4E). In addition, we also provided a pathway enrichment function in each module (Figure S4F), and users can conveniently obtain enriched pathways of user-selected WER target genes.

In the ‘Genomic browser’ module (Figure 5C), we provided an IGV genomic browser to visualize the targets. The ‘Search’ page provided search functions to help users to explore m6A2Target. Users could search any gene of interest to check whether it is a validated or potential target of any WER (Figure 5D). We also provided a ‘Download’ page (Figure 5E) to allow users to download all validated targets and potential targets in human or mouse species. Genome Reference Consortium Human Build 38 (hg38) was used as the human species reference genome, and for mouse species, Genome Reference Consortium Mouse Build 38 (mm10) was used. ‘Validated Target’ contains all targets in the validation submodule, ‘Potential Target (binding)’ contains all targets in the binding submodule and ‘Potential Target (perturbation)’ corresponds to the perturbation submodule.

Application of m6A2Target in cancer studies

Through analysis of all WER mutation rates in 33 cancer types in TCGA, we found that a well-known writer, KIAA1429, also known as VIRMA, had the highest mutation rate among all m⁶A WERs (Figure 6A). The expression level of KIAA1429 was significantly different between tumor samples and normal samples in several kinds of cancer, such as LIHC and HNSC (Figure 6B). Furthermore, the expression level of KIAA1429 was tightly associated with overall survival in multiple cancers, where high expression of KIAA1429 indicated better survival in BRCA, KICH, KIRP, PAAD, BRAC and UCEC but worse survival in KIRC and THYM (Figure 6C). Collectively, we considered that KIAA1429 played important roles in cancer progression. We then applied our database m6A2Target to explore the targets of KIAA1429 to obtain clues for downstream mechanistic studies of KIAA1429 in cancer progression. To obtain a reliable list of targets of KIAA1429, the targets with binding evidence and the targets with perturbation evidence (including m⁶A level targets and gene expression level targets) overlapped, which enabled 130 targets to be obtained (Figure 6D). The 130 targets were enriched in many cancer-related pathways, such as focal adhesion and the PI3K-AKT signaling pathway (Figure 6E). For example, the m⁶A levels of ITGAV and YWHAE in the PI3K-AKT signaling pathway were significantly hypomethylated in KIAA1429 knockdown cells compared to wild-type cells (Figure 6F). The above results suggest that m6A2Target may help researchers to investigate the m⁶A-related functions and mechanisms of genes of interest.

Discussion and conclusions

As the most abundant internal RNA modification, m⁶A is currently one of the most heavily researched subjects, and many outstanding researchers have verified its vital roles in the regulation of many fundamental biological processes and the occurrence and development of some diseases. m6A2Target provides targets of all identified WERs in human and mouse species and some related detailed comprehensive information. Currently, m6A2Target contains 1034 WER-target pairs validated by experiments, 124,772 WER-target pairs interacting with WERs and 330,102 WER-target pairs with differential expression levels, differential methylation levels, differential translation efficiencies or differential alternative splicing levels.

Compared to other m⁶A-related databases, m6A2Target is the first comprehensive and specific resource focused on the targets of m⁶A writers, erasers and readers. This resource covers human and mouse species and contains various abundances of information for targets of all identified WERs. The m6A2Target website is convenient and user-friendly, and users can browse all m⁶A-associated genes and search their targets of interest by various criteria in m6A2Target. Users can also obtain detailed information about WERs, targets and WER-target associations by simply clicking the ‘detail’ button. If users want to obtain more detailed information about this WER or its target gene, m6A2Target also provides links that help our users go to the corresponding official website, including Ensembl (http://asia.ensembl.org/index.html), HGNC (https://www.genenames.org/), UniProt (https://www.uniprot.org/) and NCBI (https://www.ncbi.nlm.nih.gov/).

m6A2Target may help researchers to characterize m⁶A modifications and their clinical significance more thoroughly. For instance, we observed that target genes of METTL14 are enriched in pathways related to the PI3K-AKT-mTOR signaling pathway and endometrial cancer pathway, and target genes of METTL3 are also enriched in endometrial cancer. He et al. proved the crucial roles of m⁶A methylation in the proliferation and tumorigenicity of endometrial cancer by regulating AKT activity [60]. These researchers found that approximately 70% of endometrial tumors presented reduced m⁶A methylation levels, which may be due to either METTL14 mutation or reduced METTL3 expression; this finding is highly consistent with the pathway enrichment results obtained with our database.

In addition, m6A2Target also has a limitation. At present, we only collected m⁶A WER-related data for the two most frequently used species in m⁶A studies, that is, human and mouse species, which may cause our results to be slightly limited. We will attempt to solve this problem in the updated versions of our database in the future.

Overall, the results of this study indicate that m6A2Target may enable researchers in the m⁶A community to obtain more detailed information about m⁶A writers, erasers, readers and their targets, thereby establishing a foundation for further research investigating m⁶A. We will update the data regularly as m⁶A-related studies are published.

Conflict of interest statement

The authors have declared that no competing interests exist.

Shuang Deng, Hongwan Zhang, Kaiyu Zhu, Xingyang Li, Ying Ye, Rui Li, Xuefei Liu, Zhixiang Zuo, Jian Zheng Sun Yat-sen University Cancer Center, State Key Laboratory of Oncology in South China and Collaborative Innovation Center for Cancer Medicine, Guangzhou, China.

Dongxin Lin Sun Yat-sen University Cancer Center, State Key Laboratory of Oncology in South China and Collaborative Innovation Center for Cancer Medicine, Guangzhou, China; Department of Etiology and Carcinogenesis, National Cancer Center/Cancer Hospital, Chinese Academy of Medical Sciences and Peking Union Medical College, Beijing, China.

References

Author notes