m6AVar: a database of functional variants involved in m6A modification.

Identifying disease-causing variants among a large number of single nucleotide variants (SNVs) is still a major challenge. Recently, N6-methyladenosine (m6A) 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 m6A modification, in order to gain a better understanding of them. Here, we report m6AVar (http://m6avar.renlab.org), a comprehensive database of m6A-associated variants that potentially influence m6A modification, which will help to interpret variants by m6A function. The m6A-associated variants were derived from three different m6A sources including miCLIP/PA-m6A-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 m6A-associated variants. We also integrated the RBP-binding regions, miRNA-targets and splicing sites associated with variants to help users investigate the effect of m6A-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 m6A-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. m6A) by changing the RNA sequences of the target sites or key flanking nucleotides.

N6-Methyladenosine (m6A) 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 m6A (known as Methylated RNA Immunoprecipitation Sequencing (MeRIP-Seq), Photo-Crosslinking-Assisted m6A Sequencing Strategy (PA-m6A-seq) and m6A individual-nucleotide-resolution cross-linking and immunoprecipitation sequencing (miCLIP)) has provided thousands of m6A sites and deep insights into the m6A machinery (10–13), revealing the essential regulatory roles of m6A in RNA splicing, miRNA function and RNA stability (7,8,14,15).

An increasing number of studies have revealed that dysregulation of m6A 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 m6A 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 m6A might be important in neurodevelopmental processes (10). To further investigate the potential pathogenesis of m6A modification, it is necessary to evaluate the effect of variants on m6A 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 m6A machinery, databases on m6A modification have emerged in recent years. In 2015, Liu et al. collected 74 samples from 22 different m6A-seq experiments and constructed MeT-DB, the first comprehensive m6A database of the mammalian transcriptome (20). Later, Sun et al. developed the RNA modification database called ‘RMBase’ that includes 226 000 m6A sites and 10 005 m5C sites (21). Although the above databases have greatly aided research of m6A functions, there is still no specific resource that would help study of influence of variants on m6A modification.

In this study, we present m6AVar (http://m6avar.renlab.org), a comprehensive database that allows the annotation, visualization and exploration of m6A-associated variants in humans and mice (Figure 1). A great number of the m6A-associated variants were derived from millions of germlines and somatic variants as well as three different m6A sources that included miCLIP experiments, PA-m6A-seq experiments, MeRIP-Seq experiments and transcriptome-wide predictions. We further annotated the m6A-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 m6A machinery and diseases.

Figure 1.

Overall design and construction of m6AVar.

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 m6A sites were derived from seven miCLIP experiments, two PA-m6A-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 m6A-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 m6A 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 m6A-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 m6A sites

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

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

The m6A 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 m6A 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 m6A sites using m6AFinder based on Random Forest algorithm from these consensus peaks (Supplementary method, Supplementary Figure S1).

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

Derivation of the m6A-associated variants

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

For m6A sites having a high confidence level, we retained the variants that located nearby the m6A sites and then looked for the variants that disrupt DRACH motif around the m6A 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.

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

For m6A sites with a low confidence level, we separately predicted the m6A 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 m6A sites in mutant sequence compared to reference sequence were defined as m6A-loss variants. In the opposite case, they were defined as m6A-gain variants.

Post-transcriptional regulation association analysis

First, m6A-associated variants were intersected with RNA-binding proteins (RBPs) regions for the same sample. In terms of miRNA targets, we matched all of the m6A-associated variants with miRNA targets to obtain the m6A-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 m6A-associated variants with these regions to obtain the splicing sites affected by the m6A-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 r2 > 0.8 in at least one of the four populations from CHB, CEU, JPT and TSI (19). Then we selected all the m6A-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 m6A-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 m6A-associated variants for human and mouse (Table 1). The m6A-associated variants with high confidence level were derived from miCLIP or PA-m6A-seq experiments. For human, there are 13 703 and 1494 high confidence level m6A-associated germline and somatic variants from dbSNP and TCGA, respectively. For mouse, there are 935 high confidence level m6A-associated germline variants from dbSNP. The m6A-associated variants with medium confidence level were derived from MeRIP-Seq experiments. For human, there are 54 222 and 7695 medium confidence level m6A-associated germline and somatic variants from dbSNP and TCGA, respectively. For mouse, there are 9,404 medium confidence level m6A-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 m6A sites. As a result, we obtained 296 933 and 29 982 low confidence level m6A-associated variants in human and mouse, respectively.

Table 1.

m6A-associated variants in m6AVar

	Human dbSNP147			Mouse dbSNP146			TCGA			Total
	Lossa variants	Gainb variants	All	Loss variants	Gain variants	All	Loss variants	Gain variants	All	Loss variants	Gain variants	All
miCLIP/PA-m6A-Seq (High confidence)	13 703	—	13 703	935	—	935	1 494	—	1 494	16 132	—	16 132
MeRIP-Seq (Medium confidence)	54 222	—	54 222	9404	—	9404	7695	—	7695	71 321	—	71 321
Prediction (Low confidence)	144 534	100 542	243 880	17 739	12 382	29 982	32 069	21 298	53 053	194 342	134 222	326 915
Total	212 360	100 542	311 706	28 065	12 382	40 308	41 243	21 298	62 227	281 668	134 222	414 241

aLoss variants were those variants resulting in loss of m6A sites in mutant sequence compared to reference sequence.

bGain variants were conversely formed.

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 m6A-associated variants from dbSNP and TCGA are related to 68 and 66 RBPs. 6371 and 338 m6A-associated variants from dbSNP and TCGA are related to 268 and 219 miRNAs. 158 469 and 39 605 m6A-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 m6A-associated variants recorded in ClinVar and GWAS, respectively. For mouse, there are 3370 m6A-associated variants related to 29 RBPs, 196 m6A-associated variants related to 173 miRNAs, and 13 382 m6A-associated variants related to the splicing sites of 7781 genes.

Table 2.

Statistics of associated data in m6AVar

	RBP-binding regions		miRNA Targets		Splicing sites		Disease-related variants
	Variants	RBPs	Variants	miRNAs	Variants	Genes	ClinVar	GWAS
Human dbSNP147	183 960 (59.02%)	68	6371 (2.04%)	268	158 469 (50.84%)	17 921	1919	372
Mouse dbSNP146	3370 (8.36%)	29	196 (0.49%)	173	13 382 (33.20%)	7781	0	0
TCGA	31 899 (51.26%)	66	338 (0.54%)	219	39 605 (63.65%)	12 273	178	168

Web interface and usage

m6AVar provides user-friendly web interfaces that enable users to browse, search and download all of the m6A-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 m6A modification in breast cancer using m6AVar. Through the ‘RsID’ search mode, users can check whether a variant of interest functionally affects the m6A status. In addition, m6A-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 m6A-associated variants in BRCA1 are presented as a table in the search results page (Figure 2B). Among them, 11, 26 and 65 m6A-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 m6A-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 m6A-associated variant. Users can select the tracks of interest to be shown, such as gene information, SNP site, m6A site, RBP binding regions, miRNA targets and the MeRIP-Seq peak level from the different samples (Figure 2F).

Figure 2.

A schematic workflow of the search interface in m6AVar. (A) m6AVar provides the four search modes of Rs ID, Gene, Region and Disease. (B) Snapshot of search results for ‘’ using the ‘Gene’ search mode. Basic information on all of the m6A-associated variants located in the ’ output presented as a table. (C and D) Distribution of m6A-associated variants in the different sources and databases. (E) Detailed information on m6A-associated variants related to post-transcriptional regulation and disease. (F) Visualization of specific m6A-associated variants with JBrowse.

Browse

The ‘Browse’ page displays: (i) Summary of m6A associated variants from three m6A 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 m6A-associated variants’ distribution in gene regions and gene types as well as other databases (Supplementary Figure S4). (iii) Browse m6A-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 m6A-associated variants that localize in the vicinity of m6A sites and potentially influence m6A modification in human and mouse. Currently, m6AVar holds ∼352 000 m6A-associated germline variants and ∼62 000 m6A-associated somatic variants, most of them were enriched in protein-coding genes (dbSNP147, 95.77%; dbSNP146, 92.12% and TCGA, 98.89%). The m6A-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 m6A-associated variants to identify the pathogenic variations contributing to dysregulation of m6A 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 m6A modification. (ii) m6AVar integrates somatic variants of 34 cancers from TCGA, which will help to reveal the potential mechanisms of m6A in cancer. (iii) m6AVar provides detail annotations and genomic coordinates for each variant and related m6A 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, m6A modification and other post-transcriptional regulation. (v) More than 2000 disease-related variants have been identified by linking the m6A-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 m6A-associated variants and search interested data by various criteria.

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

Click here for additional data file.