1. Introduction

N6-Methyladenosine (m6A) is an abundant modification in mRNA and is found within some viruses, and most eukaryotes including mammals, insects, plants and yeast. It is also found in tRNA, rRNA, and small nuclear RNA (snRNA) as well as several long non-coding RNA.

m6A and m6Am are the two 6mA (6-methyladenine)-containing nucleotides in mRNA. 6mA is found in two distinct epitranscriptomic modifications: m6A (N6-methyladenosine) and m6Am (N6,2'-Odimethyladenosine). The methyl groups that are added enzymatically to alter the function of the nucleotide are shown in red. m6A is found in 5'UTRs, coding sequences, and 3'UTRs. m6Am is found in one place in mRNAs, the first encoded nucleotide, which is adjacent to the m7G (N7-methylguanosine) cap.

The methylation of adenosine is directed by a large m6A methyltransferase complex containing METTL3 as the SAM-binding sub-unit. In vitro, this methyltransferase complex preferentially methylates RNA oligonucleotides containing RRACH motif and a similar preference was identified in vivo in mapped m6A sites in Rous sarcoma virus genomic RNA and in bovine prolactin mRNA. Wilms' tumor 1-associating protein (WTAP) may function to regulate recruitment of the m6A methyltransferase complex to mRNA targets. More recent studies have characterized other key components of the m6A methyltransferase complex in mammals, including METTL14, KIAA1429(Vir), ZFP217, RMB15, RBM15b, HAKAI(CBLL1), ZC3H13 and other proteins. Following a 2010 speculation of m6A in mRNA being dynamic and reversible, the discovery of the first m6A demethylase, fat mass and obesity-associated protein (FTO) in 2011 confirmed this hypothesis and revitalized the interests in the study of m6A by He, Chuan in The University of Chicago. A second m6A demethylase alkB homolog 5 (ALKBH5) was later discovered as well.

The biological functions of m6A are mediated through a group of RNA binding proteins that specifically recognize the methylated adenosine on RNA. These binding proteins are named m6A readers. The YT521-B homology (YTH) domain family of proteins (YTHDF1, YTHDF2, YTHDF3, YTHDC1 and YTHDC2) have been characterized as direct m6A readers and have a conserved m6A-binding pocket. Insulin-like growth factor-2 mRNA-binding proteins 1, 2, and 3 (IGF2BP1–3) are reported as a novel class of m6A readers. IGF2BPs use K homology (KH) domains to selectively recognize m6A-containing RNAs and promote their translation and stability. These m6A readers, together with m6A methyltransferases (writers) and demethylases (erasers), establish a complex mechanism of m6A regulation in which writers and erasers determine the distributions of m6A on RNA, whereas readers mediate m6A-dependent functions. m6A has also been shown to mediate a structural switch termed m6A switch.

A review named Dynamic transcriptomic m6A decoration: writers, erasers, readers and functions in RNA metabolism published on Cell Research described the current understanding of the m6A modification, particularly the functions of its writers, erasers, readers in RNA metabolism, with an emphasis on its role in regulating the isoform dosage of mRNAs.

2. Materials and Methods

2.1 m6A Immunoprecipitation and library construction

Total RNA was extracted using Trizol reagent (Invitrogen, CA, USA) following the manufacturer's procedure. The total RNA quality and quantity were analysis of Bioanalyzer 2100 and RNA 6000 Nano LabChip Kit (Agilent, CA, USA) with RIN number >7.0. Approximately more than 50 ug of total RNA was subjected to isolate Poly (A) mRNA with poly-T oligo attached magnetic beads (Invitrogen). Following purification, the poly(A) mRNA fractions is fragmented into ~100-nt-long oligonucleotides using divalent cations under elevated temperature. Then the cleaved RNA fragments were subjected to incubated for 2h at 4℃ with m6A-specific antibody (No. 202003, Synaptic Systems, Germany) in IP buffer (50 mM Tris-HCl, 750 mM NaCl and 0.5% Igepal CA-630) supplemented with BSA (0.5 μg μl−1). The mixture was then incubated with protein-A beads and eluted with elution buffer (1 × IP buffer and 6.7mM m6A). Eluted RNA was precipitated by 75% ethanol. Eluted m6A-containing fragments (IP) and untreated input control fragments are converted to final cDNA library in accordance with a strand-specific library preparation by dUTP method. The average insert size for the paired-end libraries was ~100±50 bp. And then we performed the paired-end 2×150bp sequencing on an Illumina Novaseq™ 6000 platform at the LC-BIO Bio-tech ltd (Hangzhou, China) following the vendor's recommended protocol.

2.2 m6A-seq data analysis

Firstly, Cutadapt[1] and perl scripts in house were used to remove the reads that contained adaptor contamination, low quality bases and undetermined bases. Then sequence quality was verified using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). We used bowtie[2] to map reads to the genome of Homo sapiens (Version: GRCh38.p10) with default parameters. Mapped reads of IP and input libraries were provided for R package exomePeak[3], which identifies m6A peaks with bed or bam format that can be adapted for visualization on the UCSC genome browser or IGV software (http://www.igv.org/). MEME[4] and HOMER[5] were used for de novo and known motif finding followed by localization of the motif with respect to peak summit by perl scripts in house. Called peaks were annotated by intersection with gene architecture using ChIPseeker[6]. Then StringTie[7] was used to perform expression level for all mRNAs from input libraries by calculating FPKM (FPKM=[total_exon_fragments/mapped_reads(millions)×exon_length(kB)]). The differentially expressed mRNAs were selected with log2 (fold change) >1 or log2 (fold change) <-1 and p value < 0.05 by R package edgeR[8].

Reference:

[1] Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. Embnet Journal 17 (2011).

[2] Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nature Methods 9, 357–359 (2012).

[3] Meng, J., et al. A protocol for RNA methylation differential analysis with MeRIP-Seq data and exomePeak R/Bioconductor package. Methods. 69(3), 274-281. (2014).

[4] Timothy L. Bailey, et al. MEME Suite: tools for motif discovery and searching. Nucleic Acids Research. 37.Web Server issue:202-8. (2009).

[5] Heinz, S., et al. Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and b cell identities. Molecular Cell. 38(4), 576. (2010).

[6] Yu, G., et al. ChIPseeker: an R/Bioconductor package for ChIP peak annotation, comparison and visualization. Bioinformatics. 31.14:2382. (2015).

[7] Pertea, M. et al. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nature Biotechnology. 33, 290–295 (2015).

[8] Robinson, M. D. et al. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 26(1), 139. (2010).

3. Project information

3.1 Sample information

Species: Human

Latin name: Homo sapiens

Specimens: METTL14 overexpress HEK293T cell line

3.2 Database

3.3 Bioinformatics software

4. Results

4.1 Quality control and statistics of Reads

Statistics of reads from 2 libraries was shown below. Results of FastQC, PCA analysis and sample correlation were deposited in summary/3_peak_calling and summary/9_quality_control/.

document location: summary/1_raw_data/data_stat.xlsx

4.2 Genome mapping

4.2.1 Statistics of mapped reads

document location:summary/2_mapped_stat/mapped_stat_out.xlsx

4.2.2 Statistics of mapped region

document location:

summary/2_mapped_stat/mapped_region_stat.png

summary/2_mapped_stat/mapped_region_stat.xlsx

4.3 Peak calling on whole genome

The m6A-seq approach is analogous to ChIP-seq—both are based on global identification of regions of signal enrichment. Indeed, the distribution of sequence tags in m6A-seq data is quite similar to that of ChIP-seq: in both cases, enriched regions are discrete and form sharp peaks along the genome or transcriptome. As a marriage of ChIP-Seq and RNA-Seq, m6A-seq has the potential to study the transcriptome-wide distribution of various post-transcriptional RNA modifications. Compared with other relatively well studied data types such as ChIP-Seq and RNA-Seq, the study of m6A-seq data is still at very early stage, and existing protocols are not optimized for dealing with the intrinsic characteristic of m6A-seq data. We utilized an easy-to-use protocol of using exomePeak R/Bioconductor package along with other software programs such as bowtie, ChIPseeker, MEME and HOMER for analysis of m6A-seq data, which covers raw reads alignment, RNA methylation site detection, motif discovery, differential RNA methylation analysis, and functional analysis.

4.3.1 Peak calling

Definition of table:

document location: summary/3_peak_calling/*/*_peak.xlsx

4.3.2 Accumulation of m6A-IP reads

document location:

summary/3_Peak_Calling/heatmap.png

summary/3_Peak_Calling/distProfile.png

4.4 Differential Methylation

4.4.1 Diff peak analysis

Definition of table:

document location:summary/4_diff_peak/*VS*/*VS*_diffPeak.xlsx

peak distribution within different gene contexts:

M14_OEVSM14_NC_diffPeak_pie.png

document location:summary/4_diff_peak/*VS*/*VS*_diffPeak_pie.png

4.5 Motif identification

Besides differential methylation sites, R package-exomePeak also reports all the detected RNA methylation sites in BED format, based on which the motifs of RNA methylation sites can be detected. The sequence motifs, conjectured to have a biological significance, can be identified in both ChIP-Seq and m6A-seq. The main computational difference is whether to search the reverse complement strand. While a DNA motif may appear on either strand of the two, RNA motif should appear only on the strand where the transcript is located, and thus the strand information should be kept at all times. When the strand information of the RNA fragment is lost in MeRIP-Seq unstranded library construction, it may still be derived from the information of transcripts to which the fragments (reads) are mapped. We use MEME (http://meme-suite.org/) and HOMER (http://homer.ucsd.edu/homer/motif/) to discover RRACH motif of m6A and other noval motifs in methylated sequences.

RRACH motif identification:

de novo motif prediction

document location: summary/5_motif/

4.6 Overview of gene expression

StringTie was utilized to assemble all trascripts from mapped reads (bam format) by mapping software Bowtie. Also both gene level and transcript level were estimated by StringTie using FPKM (FPKM=[total_exon_fragments/mapped_reads(millions)×exon_length(kB)]).

4.6.1 Gene expression profiling

Definition of table:

document location:summary6_transcript_expression/4_genes_fpkm_expression.xlsx

4.6.2 Transcript expression profiling

Definition of table:

Definition of class code:

document location: summary6_transcript_expression/5_isoforms_fpkm_expression.xlsx

4.6.3 Density and boxplot of gene expression

Density:

DEMO_gene_expression_density.pngDEMO_transcript_expression_density.png

Boxplot:

DEMO_gene_expression_boxplot.pngDEMO_transcript_expression_boxplot.png

document location

summary6_transcript_expression/*_gene_expression_density.png

summary6_transcript_expression/*_transcript_expression_density.png

summary6_transcript_expression/*_gene_expression_boxplot.png

summary6_transcript_expression/*_transcript_expression_boxplot.png

4.7 Differential expression analysis of mRNA

Exploratory analysis, visualization and statistical modeling are the next steps after assembling and quantifying genes and transcripts from StringTie. The assembly of input reads reconstructs the exon–intron structure of genes and their isoforms. R package edgeR takes both transcript and gene expression levels from StringTie and applies rigorous statistical methods to determine which transcripts and genes are differentially expressed between two or more experiments. edgeR also includes plotting tools as part of the R/Bioconductor package that help visualize the results.

4.7.1 Profiling of differentially expressed genes

Definition of table:

document location: summary/7_differential_expression/*VS*/*_Gene_differential_expression.xlsx

4.7.2 Profiling of differentially expressed transcripts

Definition of table:

document location: summary/7_differential_expression/*VS*/*_Transcript_differential_expression.xlsx

4.7.3 Volcano of differentially expressed genes

M14_OEVSM14_NC_Genes_volcano.pngM14_OEVSM14_NC_Transcripts_volcano.png

document location:

summary/7_differential_expression/*VS*/*_Genes_volcano.png

summary/7_differential_expression/*VS*/*_Transcripts_volcano.png

4.7.4 Heatmap of differentially expressed genes

Heatmap is a graphical representation of data where the individual z values (repeated) or Log10 (non-repeated) contained in a matrix are represented as color scheme. Calculation formula of z value is shown below. Sample-i means expression level (fpkm value).

M14_OEVSM14_NC_genes_heatmap.pngM14_OEVSM14_NC_transcripts_heatmap.png

document location:

summary/7_differential_expression/*VS*/*_genes_heatmap.png

summary/7_differential_expression/*VS*/*_transcripts_heatmap.png

4.7.5 GO enrichment analysis of differentially expressed genes

Gene ontology (GO) (http://www.geneontology.org) is a major bioinformatics initiative to unify the representation of gene and gene product attributes across all species. More specifically, the project aims to: 1) maintain and develop its controlled vocabulary of gene and gene product attributes; 2) annotate genes and gene products, and assimilate and disseminate annotation data; and 3) provide tools for easy access to all aspects of the data provided by the project, and to enable functional interpretation of experimental data using the GO, for example via enrichment analysis.

GO is part of a larger classification effort, the Open Biomedical Ontologies (OBO).

Although gene nomenclature itself aims to maintain and develop controlled vocabulary of gene and gene products, the Gene Ontology extends the effort by using markup language to make the data (not only of the genes and their products but also of all their attributes) machine readable, and to do so in a way that is unified across all species (whereas gene nomenclature conventions vary by biologic taxon).

Significant GO terms were calculated by Hypergeometric equation as shown below. TB gene number=number of total genes; TS gene number=number of differentially expressed genes in total genes; B gene number=total number of gene in all GO terms; S gene number=number of differentially expressed genes in this GO term. Those GO terms with p value<0.05 were defined as significant GO terms.

Significant GO terms:

document location: summary/7_differential_expression/*VS*/GO_Enrichment/*VS*_GO_enrichment_Gene.xlsx

Barplot of enriched GO terms:

M14_OEVSM14_NC_GO_enrichment.png

document location: summary/7_differential_expression/*VS*/1_GO_Enrichment/*VS*_GO_enrichment.png

Scatterplot of enriched GO terms:

Number of differentially expressed genes enriched in GO terms, p value and rich factor are shown in scatterplot. Rich factor=(number of differentially expressed genes in GO term)/(total number of genes in GO term). The larger rich factor is, the higher enrichment is.

M14_OEVSM14_NC_GO_enrichment_scatterplot.png

document location: summary/7_differential_expression/*VS*/1_GO_Enrichment/*VS*_GO_enrichment_scatterplot.png

4.7.6 KEGG enrichment analysis of differentially expressed genes

KEGG (Kyoto Encyclopedia of Genes and Genomes) (http://www.kegg.jp/) is a collection of databases dealing with genomes, biological pathways, diseases, drugs, and chemical substances. KEGG is utilized for bioinformatics research and education, including data analysis in genomics, metagenomics, metabolomics and other omics studies, modeling and simulation in systems biology, and translational research in drug development.

Significant KEGG pathways were calculated by Hypergeometric equation as shown below. TB gene number=number of total genes; TS gene number=number of differentially expressed genes in total genes; B gene number=total number of genes in KEGG pathways; S gene number=number of differentially expressed genes in this KEGG pathway. Those KEGG pathways with p value<0.05 were defined as significant KEGG pathways.

Significant KEGG pathways:

document location: summary/7_differential_expression/*VS*/2_KEGG_Enrichment/*VS*_KEGG_enrichment_Gene.xlsx

Scatterplot of enriched KEGG pathways:

M14_OEVSM14_NC_KEGG_enrichment_scatterplot.png

document location: summary/7_differential_expression/*VS*/2_KEGG_Enrichment/*VS*_KEGG_enrichment_scatterplot.png

KEGG pathway maps:

Graphical map objects

The KEGG pathway map is a moleculalr interaction/reaction network diagram represented in terms of the KEGG Orthology (KO) groups, so that experimental evidence in specific organisms can be generalized to other organisms through genomic information. Each map is manually drawn with in-house software called KegSketch, which generates the KGML+ file. This file is an SVG file containing graphics objects that are associated with KEGG objects (see KEGG object identifiers). Basic graphics objects in the reference KEGG pathway maps are:

boxes - ortholog (KO) groups identified by K numbers and, in metabolic maps, reactions identified by R numbers as well

circles - other molecules, usually chemical compounds identified by C numbers, but including glycans identified by G numbers

lines - reactions identified by R numbers in metabolic maps; ortholog (KO) groups identified by K numbers in global metabolism maps

and in organism specific pathway maps that are computationally generated:

boxes - genes or gene products identified by the combination of the KEGG organism code and gene identifiers

These map objects can be searched in the search box at the top of the KEGG PATHWAY page, in the search box in each pathway map, and by the KEGG Mapper tools.

Convention of map number prefix

Each pathway map is identified by the combination of 2-4 letter code and 5 digit number (see KEGG Identifiers). The prefix has the following meaning:

map - Reference pathway

ko - Reference pathway (KO)

ec - Reference pathway (EC)

rn - Reference pathway (Reaction)

org - Organism-specific pathway map

Only the first reference pathway map is manually drawn; all other maps are computationally generated. For metabolic pathways, each box (or line) in the reference map is linked to the K number (KO identifeir), the EC number, and the R number (reaction identifier). The KO, EC, and reaction maps are linked to only one of them. For all metabolic and non-metabolic maps, K numbers are converted to gene identifiers in each organism to generate organism-specific pathways.

As shown above, "map" pathways are not colored, "ko/ec/rn" pathways are colored blue, and organism-specific pathways are colored green, where coloring indicates that map objects exist and are linked to corresponding entries. For global metabolism maps, "map" pathways are fully colored, so that "ko/ec/rn" pathways and organism-specific pathways are generated by reducing the coloring indicating the absence of corresponding entries.

document location: summary/7_differential_expression/*VS*/3_KEGG_Pictures

4.8 Association analysis

4.8.1 Correlation between RNA-seq and m6A-seq

Definition of table:

document location: summary/8_diffPeak_differential_expression/*VS*/*_diffPeak_Gene_differential_expression.xlsx

4.8.2 GO enrichment analysis of diff genes/peaks

This section is similar with 4.7.5 GO enrichment analysis of differentially expressed genes

document location: summary/8_diffPeak_differential_expression/*VS*/1_GO_Enrichment

4.8.3 KEGG enrichment analysis of diff genes/peaks

This section is similar with 4.7.6 KEGG enrichment analysis of differentially expressed genes

document location: summary/8_diffPeak_differential_expression/*VS*/2_KEGG_Enrichment

5. Contact us

2575 West Bellfort Street

Suite 270

Houston, TX

77054 USA

Local (713) 664-7087

Toll Free: 1-888-528-8818

Fax: (713) 664-8181

http://www.lcsciences.com/

---

---