m6AVar: a database of functional variants involved in m⁶A modification

Abstract

Identifying disease-causing variants among a large number of single nucleotide variants (SNVs) is still a major challenge. Recently, N⁶-methyladenosine (m⁶A) has become a research hotspot because of its critical roles in many fundamental biological processes and a variety of diseases. Therefore, it is important to evaluate the effect of variants on m⁶A modification, in order to gain a better understanding of them. Here, we report m6AVar (http://m6avar.renlab.org), a comprehensive database of m⁶A-associated variants that potentially influence m⁶A modification, which will help to interpret variants by m⁶A function. The m⁶A-associated variants were derived from three different m⁶A sources including miCLIP/PA-m⁶A-seq experiments (high confidence), MeRIP-Seq experiments (medium confidence) and transcriptome-wide predictions (low confidence). Currently, m6AVar contains 16 132 high, 71 321 medium and 326 915 low confidence level m⁶A-associated variants. We also integrated the RBP-binding regions, miRNA-targets and splicing sites associated with variants to help users investigate the effect of m⁶A-associated variants on post-transcriptional regulation. Because it integrates the data from genome-wide association studies (GWAS) and ClinVar, m6AVar is also a useful resource for investigating the relationship between the m⁶A-associated variants and disease. Overall, m6AVar will serve as a useful resource for annotating variants and identifying disease-causing variants.

INTRODUCTION

Rapid improvement in high-throughput sequencing technology has resulted in the identification of millions of single nucleotide variants (SNVs) across multiple genomes. A major challenge in delineating these variants is to distinguish the functional variants from the rest. In recent years, numerous studies have been undertaken to explore disease-associated nonsynonymous SNVs that alter amino acid at the protein level (1). Nevertheless, there is growing evidence showing many synonymous SNVs, which do not alter the amino acid sequences of proteins and are considered ‘silent’ mutations, also affect the function of genes and cause various diseases, suggesting a role in transcriptional or post-transcriptional regulation (2). Many studies have shown that variants have the capacity to alter the secondary structure of RNA, influence RNA–protein interactions (3), and change the splicing sites of exonic splicing enhancers and silencers (4) as well as genetic information by means of RNA editing (5). We speculated that variants might also influence RNA modification (e.g. m⁶A) by changing the RNA sequences of the target sites or key flanking nucleotides.

N⁶-Methyladenosine (m⁶A) is a pervasive RNA modification in eukaryotes, that is involved in various biological processes such as embryonic development (6), cell apoptosis (7), spermatogenesis (8) and circadian rhythms (9). Recent development of the high-throughput sequencing techniques for m⁶A (known as Methylated RNA Immunoprecipitation Sequencing (MeRIP-Seq), Photo-Crosslinking-Assisted m⁶A Sequencing Strategy (PA-m⁶A-seq) and m⁶A individual-nucleotide-resolution cross-linking and immunoprecipitation sequencing (miCLIP)) has provided thousands of m⁶A sites and deep insights into the m⁶A machinery (10–13), revealing the essential regulatory roles of m⁶A in RNA splicing, miRNA function and RNA stability (7,8,14,15).

An increasing number of studies have revealed that dysregulation of m⁶A modification may impact various diseases. It has been found that the knockout of METTL3 in human cancer cells decreased the invasion of tumor cells (16). The activation of ALKBH5 in hypoxic breast cancer cells would promote cancer stem cell enrichment (17). In addition, previous studies have suggested that the m⁶A eraser FTO is related to metabolism dysfunction (18) and acts as an oncogenic role in Acute Myeloid Leukemia (19). Furthermore, a previous study in mice indicated that m⁶A might be important in neurodevelopmental processes (10). To further investigate the potential pathogenesis of m⁶A modification, it is necessary to evaluate the effect of variants on m⁶A modification. This will be helpful for both an understanding of the variants’ pathogenic molecular mechanisms and the identification of additional disease-causing variants.

As result of the intensifying researches and accumulating data on the m⁶A machinery, databases on m⁶A modification have emerged in recent years. In 2015, Liu et al. collected 74 samples from 22 different m⁶A-seq experiments and constructed MeT-DB, the first comprehensive m⁶A database of the mammalian transcriptome (20). Later, Sun et al. developed the RNA modification database called ‘RMBase’ that includes 226 000 m⁶A sites and 10 005 m⁵C sites (21). Although the above databases have greatly aided research of m⁶A functions, there is still no specific resource that would help study of influence of variants on m⁶A modification.

In this study, we present m6AVar (http://m6avar.renlab.org), a comprehensive database that allows the annotation, visualization and exploration of m⁶A-associated variants in humans and mice (Figure 1). A great number of the m⁶A-associated variants were derived from millions of germlines and somatic variants as well as three different m⁶A sources that included miCLIP experiments, PA-m⁶A-seq experiments, MeRIP-Seq experiments and transcriptome-wide predictions. We further annotated the m⁶A-associated variants by checking whether they localized in regions with RBP binding sites, as well as miRNA targets and splicing sites. Moreover, disease-associated data from GWAS and ClinVar database were also integrated into m6AVar, which allows users to explore the underlying relationship between the m⁶A machinery and diseases.

MATERIALS AND METHODS

Data resource

Germline and somatic variants were obtained from dbSNP and TCGA, respectively (Supplementary Table S1). We preserved those variants within the exonic regions for subsequent analysis. All of the m⁶A sites were derived from seven miCLIP experiments, two PA-m⁶A-seq experiments, 244 MeRIP-Seq experiments (Supplementary Table S2) and a transcriptome-wide prediction based on Random Forest algorithm. To identify the potential roles of m⁶A-associated variants in post-transcriptome regulation, the RBP binding sites from starBase2 (22) and CLIPdb (23) (Supplementary Table S3), the miRNA–RNA interactions from starBase2 and the canonical splice sites (GT-AG) from Ensembl annotations were collected. In addition, we obtained a large number of disease-associated SNPs from different data sets (GWAS catalog (24), Johnson and O’Donnel (25), dbGAP (26), GAD (27) and ClinVar (28)) to perform disease-association analysis. The detailed description and statistics for these data resources can be found in Supplementary Table S4.

Data preprocessing

As the raw data collected from the diverse databases utilized different data formats, it is essential to unify them under standard procedures. To do this, the genomic coordinates of all of the data resources were converted to GRCh37 for the human and GRCm38 for the mouse using the LiftOver (29). The location of each m⁶A site was then annotated by the transcript structure, including the CDS, 3′ UTR, 5′ UTR, start codon and stop codon etc. All the genomic information on the non-coding genes are from DASHR (30), miRBase (v21) (31), GtRNAdb (32) and piRNABank (33). Furthermore, all of the SNPs were annotated by ANNOVAR (updated to 1 February 2016) in two steps (34). First, we studied the conservation of evolutionary sequence in the m⁶A-associated variants using phastCons 100-way and 60-way gene conservation scores for the human and mouse respectively (35). Second, we measured the deleterious level of each variant by integrating the results from five predictors of variant function (SIFT (36), PolyPhen2 HVAR (37), PolyPhen2 HDIV (37), LRT (38) and FATHMM (39)). Each variant was scored from 0 to 5 scale by counting deleterious levels of the variants obtained by the above five methods according to their thresholds curated by the dbNSFP database (40).

Derivation of the m⁶A sites

The m⁶A sites in m6AVar were derived using three different strategies with confidence levels ranging from high to low as illustrated below:

1. The m⁶A sites having a high confidence level were extracted from the published single-nucleotide resolution m⁶A sites in the miCLIP experiments. Besides, we also obtained m⁶A sites that conformed to the DRACH (where D = A, G or U; R = G or A; H = A, C or U) motif from PA-m⁶A-seq experiments (10,11,41).

2. The m⁶A sites having a medium confidence level were predicted from the previously published MeRIP-seq data. We first downloaded all the MeRIP-Seq samples from the GEO database as raw data. Quality control was performed with FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc) and the sequencing adaptors were removed using Trimmomatic (v0.33) (42). A minimum of 25 nucleotides was required for unambiguous alignment. All qualified reads were mapped to reference genomes (GRCh37 for human and GRCm38 for mouse) by Tophat (v2.1.1) using default parameters (43). We applied three peak callers (MACS2 (44), MeTPeak (45) and Meyer's method (10)) to identify the m⁶A peaks separately. MSPC (46) was then applied to construct consensus peaks from the three methods (Supplementary method, Supplementary Table S5). We then predicted single-nucleotide resolution m⁶A sites using m⁶AFinder based on Random Forest algorithm from these consensus peaks (Supplementary method, Supplementary Figure S1).

3. In addition, to cover all potential m⁶A sites, we also obtain m⁶A sites with low confidence level by transcriptome-wide prediction using ‘m6AFinder’ (Supplementary method).

Derivation of the m⁶A-associated variants

We defined a variant as an m⁶A-associated variant by evaluating whether it has the potential to alter the DRACH motif or other sequence features essential for m⁶A modification. According to various levels of confidence, we extracted the corresponding m⁶A-associated variants as follows.

1. For m⁶A sites having a high confidence level, we retained the variants that located nearby the m⁶A sites and then looked for the variants that disrupt DRACH motif around the m⁶A sites, such as changing from D(A/G/U) to C, R(G/A) to C/T, A to C/G/U, C to G/A/U, H (A/C/U) to G.

2. For m⁶A sites with a medium confidence level, the m⁶A-associated variants were derived from the intersection between the variants and the m⁶A sites generated from MeRIP-Seq experiments. The Random Forest prediction model was subsequently applied to find the variants in the m⁶A site region that change the DRACH motif or other sequence features.

3. For m⁶A sites with a low confidence level, we separately predicted the m⁶A status for the sequence around the variants in both the reference sequence and mutant sequence by the Random Forest prediction model based on DRACH motif and other sequence features. The variants result in loss of m⁶A sites in mutant sequence compared to reference sequence were defined as m⁶A-loss variants. In the opposite case, they were defined as m⁶A-gain variants.

Post-transcriptional regulation association analysis

First, m⁶A-associated variants were intersected with RNA-binding proteins (RBPs) regions for the same sample. In terms of miRNA targets, we matched all of the m⁶A-associated variants with miRNA targets to obtain the m⁶A-associated variants which potentially impacted miRNA-target interactions. Additionally, we extracted 100 bp upstream from the 5′ splicing sites and 100 bp downstream from the 3′ splicing sites. Subsequently we matched all of the m⁶A-associated variants with these regions to obtain the splicing sites affected by the m⁶A-associated variants.

Disease association analysis

LD analysis was performed for each GWAS disease-associated SNP. We used Haploview to obtain its LD mutations with a parameter of r² > 0.8 in at least one of the four populations from CHB, CEU, JPT and TSI (19). Then we selected all the m⁶A-associated variants by mapping them with GWAS disease-associated SNPs and their LD mutations. Moreover, we also collected ClinVar data in order to annotate the m⁶A-associated variants with specific functions.

Database and web interface implementation

All the metadata in m6AVar were stored and managed in MySQL tables. The web interfaces were implemented in Hyper Text Markup Language (HTML), Cascading Style Sheets (CSS) and Hypertext Preprocessor (PHP). In order to provide visualization of all the analysis results, multiple statistical diagrams were shown by EChars and genome browser was implemented using Jbrowser (47).

RESULTS

Database content

m6AVar contains three different confidence levels of m⁶A-associated variants for human and mouse (Table 1). The m⁶A-associated variants with high confidence level were derived from miCLIP or PA-m⁶A-seq experiments. For human, there are 13 703 and 1494 high confidence level m⁶A-associated germline and somatic variants from dbSNP and TCGA, respectively. For mouse, there are 935 high confidence level m⁶A-associated germline variants from dbSNP. The m⁶A-associated variants with medium confidence level were derived from MeRIP-Seq experiments. For human, there are 54 222 and 7695 medium confidence level m⁶A-associated germline and somatic variants from dbSNP and TCGA, respectively. For mouse, there are 9,404 medium confidence level m⁶A-associated germline variants from dbSNP. In addition, a genome-wide prediction based on Random Forest algorithm was performed for the sequences around all the collected variants from dbSNP and TCGA to find the variants that cause potential gain or loss of m⁶A sites. As a result, we obtained 296 933 and 29 982 low confidence level m⁶A-associated variants in human and mouse, respectively.

m⁶A-associated variants in m6AVar

Moreover, m6AVar contains many associated data, such as RBPs, miRNA and splicing sites, as well as disease information (Table 2). For human, 183 960 and 31 899 m⁶A-associated variants from dbSNP and TCGA are related to 68 and 66 RBPs. 6371 and 338 m⁶A-associated variants from dbSNP and TCGA are related to 268 and 219 miRNAs. 158 469 and 39 605 m⁶A-associated variants from dbSNP and TCGA are related to the splicing sites of 17 921 and 12 273 genes. Moreover, there are 2097 and 540 disease-related m⁶A-associated variants recorded in ClinVar and GWAS, respectively. For mouse, there are 3370 m⁶A-associated variants related to 29 RBPs, 196 m⁶A-associated variants related to 173 miRNAs, and 13 382 m⁶A-associated variants related to the splicing sites of 7781 genes.

Statistics of associated data in m6AVar

Web interface and usage

m6AVar provides user-friendly web interfaces that enable users to browse, search and download all of the m⁶A-associated variants in the database.

Search

m6AVar provides four modes to query the database, i.e. by RsID, Gene, Chromosome region and Disease. Here, we illustrate an example to show how to utilize m6AVar by search function (Figure 2). We sought to undertake an investigation of m⁶A modification in breast cancer using m6AVar. Through the ‘RsID’ search mode, users can check whether a variant of interest functionally affects the m⁶A status. In addition, m⁶A-associated variants in known breast cancer-related genes may be obtained by using the ‘Gene’ mode (Figure 2A). Taking the human tumor suppressor gene BRCA1 as an example, 102 m⁶A-associated variants in BRCA1 are presented as a table in the search results page (Figure 2B). Among them, 11, 26 and 65 m⁶A-associated variants were derived from miCLIP (high confidence), MeRIP-Seq (medium confidence) and prediction (low confidence), respectively (Figure 2C). A statistical plot shows the number of germline and somatic m⁶A-associated variants (Figure 2D). Users may obtain the related RBPs, miRNA targets, splicing sites and diseases from the detailed information on each variant (Figure 2E). Furthermore, m6AVar also allows users to find more disease-related variants directly through ‘Disease’ mode. In order to facilitate follow-up experimental studies, it allows users to customize results with the advanced search and to sort the table by clicking on the column names. Furthermore, we applied the JBrowse Genome Browser to visualize every m⁶A-associated variant. Users can select the tracks of interest to be shown, such as gene information, SNP site, m⁶A site, RBP binding regions, miRNA targets and the MeRIP-Seq peak level from the different samples (Figure 2F).

Browse

The ‘Browse’ page displays: (i) Summary of m⁶A associated variants from three m⁶A sources (with a high, medium and low confidence level) (Supplementary Figure S2). (ii) Statistical graphs showing the overall functional gain and loss variants’ frequency distribution in a circular layout (Supplementary Figure S3), and m⁶A-associated variants’ distribution in gene regions and gene types as well as other databases (Supplementary Figure S4). (iii) Browse m⁶A-associated variants by gene types. To retrieve data more efficiently, various filters, such as gene types, associations and confidence levels are provided (Supplementary Figure S5).

Download

all data in the database can be downloaded from the ‘Download’ page, and a detailed introduction of m6AVar database as well as tutorial are available on the ‘Help’ page.

DISCUSSION

m6AVar is a comprehensive database of the m⁶A-associated variants that localize in the vicinity of m⁶A sites and potentially influence m⁶A modification in human and mouse. Currently, m6AVar holds ∼352 000 m⁶A-associated germline variants and ∼62 000 m⁶A-associated somatic variants, most of them were enriched in protein-coding genes (dbSNP147, 95.77%; dbSNP146, 92.12% and TCGA, 98.89%). The m⁶A-associated variants that can potentially affect RBP-binding regions, miRNA-targets and splicing sites were discovered by systematic association analyses. Furthermore, diseased-related variants from GWAS and ClinVar have been intersected with the m⁶A-associated variants to identify the pathogenic variations contributing to dysregulation of m⁶A modification.

m6AVar has the following advantages in comparison with MeT-DB and RMBase. (i) m6AVar is a specific database dedicated to the investigation of the functional association between variants and m⁶A modification. (ii) m6AVar integrates somatic variants of 34 cancers from TCGA, which will help to reveal the potential mechanisms of m⁶A in cancer. (iii) m6AVar provides detail annotations and genomic coordinates for each variant and related m⁶A site. This will help biologists determine its relevant biological features. (iv) m6AVar integrates the results from association analyses with RBP-binding regions, miRNA-targets and splicing sites, revealing the potential relationship among variants, m⁶A modification and other post-transcriptional regulation. (v) More than 2000 disease-related variants have been identified by linking the m⁶A-associated variants with GWAS and ClinVar data, which may assist the community in identifying the functional disease-causing variants. (vi) m6AVar is a user-friendly database with multiple statistical diagrams and genome browser through which users can browse all of the m⁶A-associated variants and search interested data by various criteria.

In conclusion, m6AVar provides useful information on m⁶A-associated variants to help experimental biologists interpret the disease-related variants by m⁶A function and explore the molecular mechanism of m⁶A modification. m6AVar will be continually updated whenever new high-throughput m⁶A sites data and variants data are made available in public databases.

SUPPLEMENTARY DATA

Supplementary Data are available at NAR Online.

FUNDING

National Natural Science Foundation of China [31771462, 81772614, 31471252, U1611261]; National Key Research and Development Program [2017YFA0106700]; Guangdong Natural Science Foundation [2014TQ01R387, 2014A030313181]; Science and Technology Program of Guangzhou, China [201604020003]. Funding for open access charge: National Natural Science Foundation of China [31471252].

Conflict of interest statement. None declared.

REFERENCES

Author notes