CVm6A: A Visualization and Exploration Database for m⁶As in Cell Lines

Abstract

N6-methyladenosine (m⁶A) has been identified in various biological processes and plays important regulatory functions in diverse cells. However, there is still no visualization database for exploring global m⁶A patterns across cell lines. Here we collected all available MeRIP-Seq and m⁶A-CLIP-Seq datasets from public databases and identified 340,950 and 179,201 m⁶A peaks dependent on 23 human and eight mouse cell lines respectively. Those m⁶A peaks were further classified into mRNA and lncRNA groups. To better understand the potential function of m⁶A, we then mapped m⁶A peaks in different subcellular components and gene regions. Among those human m⁶A modification, 190,050 and 150,900 peaks were identified in cancer and non-cancer cells, respectively. Finally, all results were integrated and imported into a visualized cell-dependent m⁶A database CVm6A. We believe the specificity of CVm6A could significantly contribute to the research for the function and regulation of cell-dependent m⁶A modification in disease and development.

1. Introduction

As one important post-transcriptional modification, N6-methyladenosine (m⁶A) was largely discovered by high throughput sequencing in recent years [1,2,3]. m⁶A was identified with consensus sequence surrounding m⁶A site RRACH (R=G or A, H=A, C or U) and conserved in human, mouse, chimpanzee and even in plants [1,4,5]. m⁶A was also found to exist in bacterial and archaeal species [6]. The abundance of m⁶A is reported as being correlated with evolutionarily conserved region of genome [2]. m⁶A modification was a reversible status mediated by methyltransferases METTL3/ METTL14/ WTAP complex [7], demethylases FTO/ALKBH5 [2,8] and recognized by m⁶A binding proteins YTH (YT521-B homology) domain family/HNRNPA2B1 [9,10], which were called writer, eraser and reader, respectively.

m⁶A can regulate the multiple biological functions in spatial and temporal [11]. m⁶A methyltransferase complex controls the neuronal functions and fine-tuning sex determination in Drosophila [12]. m⁶A also acts as a regulator at molecular switches in murine naive pluripotency for proper lineage priming and differentiation [13]. The existence of m⁶A in lncRNA XIST mediated the gene silencing on X chromosome. Knockdown of m⁶A methytransferase METTL3 can impair XIST-mediated gene silencing [14]. m⁶A RNA can recruit DNA polymerase k (Pol k) to facilitate repairing of ultraviolet-induced DNA damage [15]. Furthermore, m⁶A could alter RNA structure to affect RNA-protein interactions in cells [16]. The m⁶A-driven gene network was already constructed and the dynamic interactions between m⁶A related methyltransferases and demthylases were established [17]. The deficiency of m⁶A modification led to various diseases, such as obesity, cancer, type 2 diabetes mellitus, infertility and developmental arrest, etc. [18].

In previous researches, m⁶A was discovered mainly located near stop codons, large internal exons and 3’UTR (3’-Untranslated region), as well as in CDS (Coding sequence), transcriptional start sites and intron [1,2,19]. Dynamic m⁶A modification could affect translation status and lifetime of mRNA in Hela [20]. Several lncRNAs also hosted m⁶A modification [1,2] and long intergenic noncoding RNAs (lincRNAs) established significantly higher level than mRNAs in B-cell lymphoblastoid cell line GM12878 [21]. In CD4T, m⁶A modification happened on HIV-1 RNA could regulate viral infection [22].

Though m⁶A patterns were analyzed in different cells independently, the global patterns across those cells were not well summarized. Several databases collected and detected m⁶A from public datasets, such as RMBase [23] and MeT-DB [24]. However, RMBase and MeT-DB were not focused on cell-dependent m⁶A. For examples, MeT-DB only included m⁶A datasets from a portion of wild type cell lines, and RMBase included m⁶A sites from various samples without indicating the cell sources. To better understand the function of m⁶A in cellular biological processes, a more specific database is needed for exploring and comparing the distribution and patterns of m⁶A in different cell lines. Here, using latest public datasets, we collected MeRIP-Seq and m⁶A-CLIP-Seq datasets from 23 human cell lines and eight mouse cell lines from previous work, and inspected the global patterns of m6As across those cell lines, including the distribution and abundance of m⁶A modification in lncRNA or mRNA, different subcellular location and gene regions. The m⁶A patterns from cancer or non-cancer cell lines were also classified. Moreover, validated m⁶A sites from previous experiments were also collected and summarized. All results were imported into a cell-dependent m⁶A database CVm6A (http://gb.whu.edu.cn:8080/CVm6A) providing a visualization interface for searching and comparing the m⁶A patterns in different cell lines, which could contribute to the function and regulation research of m⁶A in disease and development.

2. Data Collection and Database Content

2.1. Cell Line Samples in CVm6A

Previous studies showed that MeRIP-Seq (Methylated RNA Immunoprecipitation sequencing) [20], miCLIP-Seq (m⁶A individual-nucleotide-resolution cross-linking and immunoprecipitation sequencing) [25] and PA-m⁶A-Seq (Photo-crosslinking-assisted m⁶A-seq) [26] could be used for detecting m⁶A modification in transcriptomic level. Therefore, we collected all available MeRIP-Seq, miCLIP-Seq and PA-m⁶A-Seq datasets with total RNA or PolyA enriched library construction from NCBI GEO database (http://www.ncbi.nlm.nih.gov/GEO). In total, 47 samples from 23 human cell lines and 22 samples from 8 mouse cell lines were collected (Table S1).

2.2. Identification of Cell m⁶A Peaks

For MeRIP-Seq datasets, both reads from IP (Immunoprecipitation) and Input samples were mapped to human (hg38 version) and mouse (mm10 version) genome separately via Hisat2 [27]. Mapped reads with MAPQ <30 were filtered by samtools [28], and removed PCR duplicates using Picard (http://broadinstitute.github.io/picard, v2.16.0). Then m⁶A peaks were called and enrichment score of each peak was calculated by MeTPeak [29]. m⁶A sites from miCLIP-Seq and PA-m⁶A-Seq were collected from previous works [19,25,26]. Gene annotation of GENCODE (GRCh38 release 28 and GRCm38 release M20) including 35,048 human genes and 31,237 mouse genes were used to annotate m⁶A sites or peaks. Detailed pipeline was included in Supplemental Method. In all cell lines, total 340,950 m⁶A peaks from 16,950 human genes and 179,201 m⁶A peaks from 14,360 mouse genes were identified. In human cell lines, we retrieved 6345 (H1299) ~ 23,052 (A549) m⁶A peaks, and 2562 (H1299) ~ 6838 (GSC-11) genes with m⁶A modification (Figure 1A). In mouse cell lines, 6833 (3T3-L1) ~ 20,892 (iPSC) m⁶A peaks and 2882 (SC) ~ 7125 (NSC) genes with m⁶A modification were identified (Figure 1B).

2.3. Prediction of m⁶A lncRNA and mRNA

All m⁶A peaks were mapped to mRNA and lncRNA using GENCODE gene annotation (GRCh38 release 28 and GRCm38 release M20) [30] via Bedtools [31]. To view the similarities and differences of m⁶A modification in lncRNA and mRNA, those all m⁶A genes were separated into lncRNA or mRNA groups. In human cell lines, there were 225 (HEK293) ~ 2627 (LCL) peaks from lncRNA, while 6044 (H1299) ~ 22,630 (A549) peaks were from mRNA (Figure 1C). We also checked the enrichment score of those m⁶A peaks from lncRNA and mRNA. The average enrichment scores of lncRNA were from 2.37 (U2OS) ~ 9.55 (ESC), while in mRNAs, the scores were from 2.36 (U2OS) ~ 11.15 (iPSC) (Figure S1A). In mouse cell lines, there were 101 (3T3-L1) ~ 1,243 (iPSC) peaks from lncRNA, while 6732 (3T3-L1) ~ 19,649 (iPSC) peaks were from mRNA (Figure 1D). The average enrichment scores of lncRNA and mRNA from mouse cell lines were also established (Figure S1B).

2.4. Prediction of Subcellular Location

To view the location of m⁶A in subcellular component, we acquired the subcellular location information of genes from public database Hum-mPloc 3.0 [32] and Euk-mPLoc [33] and classified m⁶A peaks into subcellular components according to their annotated genes. Total 309,137 peaks were located into different components in all cell lines. Across human cell lines, 2439 (HEK293) ~ 8866 (A549) of peaks were located in Cytoplasm, 2772 (H1299) ~ 9195 (A549) of peaks located in Nucleus and 447 (H1299) ~ 2878 (A549) of m⁶A peaks located in Plasma membrane. Other peaks were distributed in other subcellular location: Centrosome, Cytoskeleton, Endoplasmic reticulum, Endosome, Golgi apparatus, Lysosome, Mitochondrion, Peroxisome, and Extracellular (Figure 1E). The average enrichment scores of those peaks from different components were also calculated (Figure S1C). Similar distributions were observed in mouse cell lines (Figure 1F and Figure S1D).

2.5. Prediction of Gene Regions

Previous works revealed m⁶A modification were not uniform on different gene regions [1,2,20]. In CVm6A, we separated all annotated genes into 6 regions: TSS, 5′UTR, CDS, Stop codon, 3′UTR and Intron and located m⁶A peaks in these regions. The middle site of each peak from MeRIP-Seq or the precise m⁶A site from m⁶A-CLIP-Seq was used for the localization. For TSS and stop codon, a 200 bp window (±100 bp surrounding the coordinates) were allowed to locate the m⁶A site according to previous work [1,2]. Across human cell lines, 2668 (H1299) ~ 11,209 (A549) of m⁶A peaks were distributed in CDS, while 1517 (H1299) ~ 10,706 (A549) were located in 3′UTR and 710 (CD8T) ~ 4608 (A549) were distributed in Stop codon (Figure 1G). We also checked the average enrichment scores distributed in those regions. In all cell lines, stop codon had average scores from 2.34 (U2OS) to 12.5 (iPSC), while the average scores in CDS were 2.41 (U2OS) ~ 11.7 (iPSC) and the average scores in 3′UTR were 2.32 (U2OS) ~ 11.6 (iPSC) (Figure S1E). The m⁶A distribution of gene regions in mouse cell lines was also established (Figure 1H and Figure S1F).

2.6. Classification of Cancer and Non-Cancer m⁶A

To view the association of m⁶A and diseases, all human m⁶A peaks were classified into cancer and non-cancer groups. Overall, 190,050 m⁶A peaks from 14,628 genes were identified in 12 cancer cell lines and 150,900 peaks from 14,346 genes were identified in 11 non-cancer cell lines (Figure S1G).

2.7. Validated m⁶A Sites

Previous research had validated several m⁶A modification in cell lines. To enhance the usability of CVm6A, we collected those m⁶A sites from m⁶A-RIP or m⁶A-CLIP experiment from previous literatures. Totally, CVm6A contains validated m⁶A modification in 96 genes, which were identified in 25 cell lines from human, mouse, zebrafish and Drosophila (Table S2).

3. Database Organization and Web Interface

All the analyzed results, including peak region, gene type, gene region, subcellular location, conditions and library types associated with m⁶A peaks were integrated into a set of interactive MySQL tables. Laravel–an open-source web framework based in PHP (https://laravel.com) and JavaScript library were used to construct the CVm6A database. The web interface of CVm6A is summarized in Figure 2.

3.1. Browse Page

On this page, users can browse all m⁶A peaks from 23 human and eight mouse cell lines. All information, including the peak region, strand, enrichment score, gene symbol, gene type (mRNA/lncRNA), gene region (CDS, 3’UTR, et al.), subcellular location (plasma_membrane, nucleus, et al.), cell line, condition (cancer/non-cancer) and library type for each peak are displayed in the table (Figure 2A). Peak region and gene symbol are linked to the visualization page of m⁶A peaks located in this gene (Figure 2B).

3.2. Visualization Page

On this page, users can view all m⁶A peaks distributed in selected annotated gene. All peaks in current gene are displayed with dark yellow color. In the top left corner, user can also select special cell line in the search box. All peaks in the selected cell line are displayed with brown color both in the figure and table. The summit of each peak is determined according to the enrichment score. The gene structure with exon (blue box) and intron (gray line) is displayed below peak figures. If the selected gene has more than one transcript, all transcripts are displayed. The relative location of each peak and the annotated gene are placed according to genomic coordinates. Mouse hover over on or click each peak can display the peak region, enrichment score, gene region, subcellular location, cell line, condition and library type of this peak. The table below the figures includes the coordinates of annotated transcripts and detailed information of all peaks in the selected gene (Figure 2B).

3.3. Search Page

On this page, users can search m⁶A peaks by gene symbol, cell line or genomic coordinates. While users select gene symbol or cell line, all peaks in this gene or cell line will be displayed in tables below and can be exported into files. Batch search by gene symbol is also provided. The search function by genomic coordinates supports fuzzy search. All peaks surrounding the input coordinates will be displayed (Figure 2C).

4. Summary and Future Directions

CVm6A collects available MeRIP-Seq and m⁶A-CLIP-Seq datasets in human and mouse cell lines, and provides a visualized m⁶A database to benefit functional studies of m⁶A in cell lines. Those samples include the most frequently used cell lines in previous researches. For example, users working on stem cells could explore the m⁶A modification in ESC, MSC, NPC and iPSC and can compare the distribution of m⁶A peaks in nucleus and other subcellular components, as well as 3′UTR and other gene regions. Users working on the immune cell lines could inspect the distribution in CD4T, CD8T and LCL. Moreover, more than ten cancer cell lines are included in CVm6A, which allow researchers to study the potential function of m⁶A in cancers. CVm6A also predicts the enrichment score of each peak, which allow users to check the abundance of m⁶A. Due to the limited m⁶A-Seq datasets for total RNA, only two total RNA datasets are included in the current version, which cannot thoroughly establish the distribution of m⁶A on lncRNA and other non-polyA RNAs. We will update CVm6A when more sequencing data from other library types becomes available.