1. m6A RNA methylation in the bone marrow of childhood B-ALL
ALKBH5, FTO, METTL14, YTHDF2 and YTHDF3 as the key regulators of m6A RNA methylation are expressed in the bone marrow of childhood B-ALL by qRT-PCR test (Figure 1).
The 5 genes in the datasets from Oncomine database were screen and the P values were collected and analyzed (Figure 2 and Additional file: Data S1).
2. Transcriptome-wide detection of m6A modification in childhood B-ALL
We obtained 4594662, 4388434, 8565686 and 6755462 high-quality reads from group A, B, C and D respectively after filtering out low-quality data (table.1).
3. The distribution of m6A modification exhibits a distinct topology in the childhood B-ALL
To understand the preferential location of m6A in transcripts, we next investigated the metagene profiles of m6A peaks in the entire transcriptome. We observed that m6A peaks were markedly correlated with two distinct coordinates: immediately following the start of the 5′untranslated regions (5′UTRs), and near the end of the coding sequence (CDS) and the start of the 3′untranslated regions (3′UTRs) in the four group (Fig. 3A); meanwhile, the end CDS of peaks were more pronounced than the start of the 5′UTRs peaks. This dominant m6A enrichment was not found in the remaining mRNA. To assess the enrichment methodically, we assigned each m6A peak to one of six no overlapping transcript segments: 5′UTRs, start codon, CDS, stop codon, 3′UTR, and unmapped (Fig. 3B). Most (~80%) of the m6A peaks were within genic regions, and almost 40% of genic peaks were localized near the stop codon. The topological patterns distributing within genes were highly similar in the four groups, suggesting that recognition of motif for m6A methylation was conserved among human tissues.
4. The General features of m6A modification in childhood B-ALL
To determine whether our identified m6A peaks shared the conserved RRACH motif (where R represents purine, A is m6 A and H is a non-guanine bas), a motif search was performed to confirmed motifs enriched in regions Surrounding m6A peaks. Clustering of significantly enriched sequences perfectly recapitulated the previously established m6A RRACH consensus sequence in four groups (Fig.4). The identification of a strong consensus reinforces the authenticity of the discovered m6A peaks, and supports the existence of predominant methylation machinery. There is about 70% methylated transcripts exhibited one or two m6A sites. The percentage of more than 2 m6A modifications in each mRNA is approximately 30%. And over 20% contained four or more sites (Fig. 5), which was much higher than previously reported in human.
5. The commonly and selectively methylated peak between groups.
5.1 Commonly and selectively methylated genes in Ph positive ALL and Ph negative ALL
Comparing the methylated genes intra the Ph positive ALL or between the Ph positive and the Ph negative ALL, we discovered that more than 12,000 tissue specific methylated genes (TSMGs) were methylated intra the Ph positive ALL and between the Ph positive ALL and the Ph negative ALL (Fig.6). Further investigation revealed the commonly methylated is more than the selectively methylated (SMG) in either Ph positive or negative ALL. Obviously, selectively methylated genes (SMG) were only modified in certain tissues. By statistical test analysis, we found the numbers of SMGs in Ph negative ALL were more than the Ph positive ALL were differentially expressed between the two groups.
5.2 m6A-containing genes are associated with transcriptional factors and involved in important biological pathways
To further predict general functional insights about m6A in the bone marrow from childhood B-ALL at the time of diagnosis. We systematically screened these commonly m6A modified genes which appear on four groups and identified the enriched gene ontology (GO) terms with the help of the GO consortium database. The commonly methylated genes intra the Ph positive ALL or between the Ph positive and the Ph negative ALL, we found the major molecular functions (MF) were mainly involved in a variety of binding activity, transferase, transcription regulatory, and transferase activity; the cellular components (CC) were responsible for intracellular, nucleus and membrane-bounded organelle and the biological processes (BP) were responsible for metabolic process (Fig.7A-CC, Fig.7B-MF, Fig.7C-BP). Collectively, these data demonstrated that m6A -containing RNAs were involved in a variety of biological pathways relevant to tissue development or cellular signaling.
5.3 The KEGG pathway analysis
Based on KEGG pathway analysis, certain transcripts presenting higher methylation was related to pathways of infection, shigellosis, and ubiquitin mediated proteolysis and chronic myeloid leukemia (Fig.8).
6. The differential m6A methylation between Ph positive and Ph negative
6.1 the peak number between groups
6.2 The heat map
Differentially expressed peak between Ph+ (A/B/C) and Ph- (D) childhood ALL at diagnosed. The up regulations of specific m6A peaks are highlighted in green, the down regulations are highlighted in red (Fig.9).
7. Correlation of m6A modification with gene features and transcription levels
Next, we performed MeRIP-Seq and RNA-Seq analysis of Ph+ and Ph- B-ALL bone marrow (Table 4, Fig.10 and Additional file: Data S2). Intriguingly, 9471 transcripts were overlapped in the transcriptome RNA sequencing and MeRIP-seq data (Fig.10)
We found that the number of transcripts displayed m6A methylation changes was almost the same as the number of genes that were differentially expressed. Among genes with m6A modification, more are upregulated than downregulated in Ph positive and Ph negative childhood B-ALL. The data obtained from the B-ALL bone marrow samples suggest that m6A RNA methylation tends to positively correlate with gene expression in a large fraction of transcripts in Ph positive and Ph negative childhood B-ALL (Table.5 and Fig.11).
Furthermore, To confirm the DEC expression in our microarray data, three upregulated (PCDHGB2, ZNF99 and PCDHGB4) and two downregulated (SHISA2 and TRAC) genes were randomly chosen for validation of our analysis above by qRT-PCR analysis (Figure 12).
8. GO and KEGG analysis of changes in gene expression
We performed GO and KEGG pathway analyses to investigate the biological significance of these genes with differentially methylated m6A peaks and differential expression between the Ph positive (CK1) and Ph negative (CK2) B-ALL. Go analysis we found the CC were responsible for plasma membrane; the MF were mainly involved in a variety of binding activity and the biological processes were responsible for metabolic process. The KEGG pathway analysis revealed that the changes in gene expression mainly occur in microRNA, ECM-receptor interaction and secretion process (Figure 13).