Funm6AViewer
Identification and visualization of functional differential m6A methylation genes (FDmMGenes) and single base DmM sites. We also developed the Funm6AViewer webserver which is freely available from http://funm6aviewer.rnamd.com/ .
1. Installation
Funm6AViewer depends on GenomicFeatures, Guitar, trackViewer, DESeq2, STRINGdb, TxDb.Hsapiens.UCSC.hg19.knownGene, org.Hs.eg.db R packages and please make sure they are installed before installing Funm6AViewer. An R version >= 3.6 is suggested.
Install the required packages
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install(c("GenomicFeatures", "GenomicAlignments", "Rsamtools", "Guitar", "trackViewer", "DESeq2", "apeglm", "STRINGdb",
"TxDb.Hsapiens.UCSC.hg19.knownGene", "org.Hs.eg.db"), version = "3.10")
Note that if your are using an R version lower than 3.6, please install the corresponding Bioconductor version and if you installed a different version of Bioconductor packages rather than version 3.10, you should check the STRINGdb supported version and assign it to version parameter of funm6aviewer function. For example, by default, the version of funm6aviewer is set as version = "10" which is corresponding to STRINGdb version = 10 and Bioconductor version = 3.10; if you installed Bioconductor version = 3.11, the corresponding STRINGdb version should be "11", then the version should be set as version = "11".
Install Funm6AViewer
if (!requireNamespace("devtools", quietly = TRUE))
install.packages("devtools")
devtools::install_github("NWPU-903PR/Funm6AViewer")
Test the installation
To test the installation, please run the following toy example:
library(Funm6AViewer)
dminfo <- system.file("extdata", "DMinfo_toy.xls", package="Funm6AViewer")
dminfo <- read.table(dminfo, header = TRUE, stringsAsFactors = FALSE)
bamreadsgr <- system.file("extdata", "bamgrlist_toy.RData", package="Funm6AViewer")
load(bamreadsgr)
re <- coverageplot(dminfo = dminfo, grlist = grlist, intrested_gene = "MYC")
2. Data required
Funm6AViewer adopted FunDMDeep-m6A to idenify functional DmM genes which required 4 PPI networks. Then to run Funm6AViewer, users need to firstly download this data from http://180.208.58.66/Funm6AViewer/Download/Funm6AViewer_data.zip
3. One step usage of Funm6AViewer
funm6aviewer takes single base DmM sites information, gene DE information as input and output:
-
Functional DmM genes (FDmMGenes);
-
DmM sites distribution on RNA;
-
Counts of DmM sites on different RNA regions;
-
DmM sites and reads coverage on interested gene;
-
Function enrichment of FDmMGenes;
-
Context specific function of interested genes;
-
DmMGene's MSB score along with DE score;
-
Network of FDmMGene's MSB neighbors.
Following is an example to achieve these using funm6aviewer:
library(Funm6AViewer)
Get input data:
dminfo <- system.file("extdata", "DMinfo_toy.xls", package="Funm6AViewer")
deinfo <- system.file("extdata", "DEinfo_toy.xls", package="Funm6AViewer")
dminfo <- read.table(dminfo, header = TRUE, stringsAsFactors = FALSE)
deinfo <- read.delim(deinfo, header = TRUE, stringsAsFactors = FALSE)
dminfo contains the position annotation and log2 foldchange of DmM sites. It can be extracted from the result of DMDeep-m6A package using summarydmdeepm6A (see 9.2 for more details). Alternatively, users can use any other method to make it as the following formate:
head(dminfo)
## chr chromStart chromEnd name score strand log2fd
## 1 chr1 155160832 155160833 4582 0.9137714 - 2.7950754
## 2 chr1 171505224 171505225 23215 0.9431386 + 0.7589479
## 3 chr1 241767682 241767683 23596 0.8125095 - 1.0680801
## 4 chr1 243418399 243418400 9859 0.9652586 - 1.7805035
## 5 chr1 8073372 8073373 54206 0.8477583 - 2.5678589
## 6 chr1 8073689 8073690 54206 0.8089832 - 1.1375522
The 'name' column can be entrez gene ID or gene symbol.
deinfo contains the differential expresion p-value and fdr for genes. It can be made using makegrreadsfrombam and getdeinfo (see 9.2 for more details), or users can use any other method to make it as the following formate:
head(deinfo)
## name pval padj log2FoldChange
## 1 1 7.578860e-01 8.990406e-01 -0.07382034
## 2 100 6.958592e-01 8.695820e-01 0.07649361
## 3 1000 4.155368e-06 6.420489e-05 -0.73757946
## 4 10000 2.043250e-02 9.424864e-02 0.40075777
## 5 100009676 4.524888e-01 7.148682e-01 -0.16186352
## 6 10001 6.708161e-01 8.566831e-01 0.08459812
The 'name' column can be entrez gene ID or gene symbol.
bamreadsgr can be generated using makegrreadsfrombam from the MeRIP-Seq data in bam formate (see 9.2 for more details).
bamreadsgr <- system.file("extdata", "bamgrlist_toy.RData", package="Funm6AViewer")
load(bamreadsgr)
siggene <- c("CCNT1", "MYC", "BCL2")
permutime <- 1000
The datapath is the filepath where the required PPI data saved and the enrich_input_directory is the filepath passed to string_db, the GO and KEGG function annotation data will be downloaded to this path. All these required data can be downloaded from http://180.208.58.66/Funm6AViewer/Download/Funm6AViewer_data.zip
datapath <- "F:/Funm6A_package/data"
enrich_input_directory <- "F:/Funm6A_package/data"
savepath <- getwd()
re <- funm6aviewer(dminfo, deinfo, grlist, intrested_gene = siggene, permutime = permutime, version = "10",
datapath = datapath, enrich_input_directory = enrich_input_directory, savepath = savepath)
The version parameter is passed to STRINGdb, it depends on which version is supported by STRINGdb. The results will be saved to savepath.
If you are using other genomes, you need to install txdb and orgdb annotation for the corresponding genome. Taking mouse mm9 genome as an example, you should firstly install the genome annotation:
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install(c("TxDb.Mmusculus.UCSC.mm9.knownGene",
"org.Mm.eg.db"))
And assign the annotation to txdb, orgdb and orgsymbol when running funm6aviewer:
library(TxDb.Mmusculus.UCSC.mm9.knownGene)
library(org.Mm.eg.db)
re <- funm6aviewer(dminfo = dminfo,
deinfo = deinfo,
grlist = grlist,
intrested_gene = siggene,
txdb = TxDb.Mmusculus.UCSC.mm9.knownGene,
orgdb = org.Mm.eg.db,
orgsymbol = org.Mm.egSYMBOL,
permutime = permutime,
datapath = datapath,
enrich_input_directory = enrich_input_directory,
savepath = savepath)
4. DmM sites plot for interested gene
Users can use dmsiteplot to visaulize the DmM sites on their interested genes and their isoforms.
Get input:
dminfo <- system.file("extdata", "DMinfo_toy.xls", package="Funm6AViewer")
dminfo <- read.table(dminfo, header = TRUE, stringsAsFactors = FALSE)
siggene <- c("MYC")
Make plot:
re <- dmsiteplot(dminfo = dminfo, intrested_gene = siggene)
5. Reads coverage plot for interested gene
Users can use coverageplot to visaulize the reads coverage of DmM sites on their interested genes.
Get input:
dminfo <- system.file("extdata", "DMinfo_toy.xls", package="Funm6AViewer")
dminfo <- read.table(dminfo, header = TRUE, stringsAsFactors = FALSE)
bamreadsgr <- system.file("extdata", "bamgrlist_toy.RData", package="Funm6AViewer")
load(bamreadsgr)
siggene <- c("MYC", "CCNT1")
Make plot:
re <- coverageplot(dminfo = dminfo, grlist = grlist, intrested_gene = siggene)
Users can zoom the reads coverage near the differential m6A sites by setting the zoom_region as following:
re <- coverageplot(dminfo = dminfo, grlist = grlist, intrested_gene = "CCNT1", zoom_region = 100)
Then 100 bp up- and down-stream of the sites on gene CCNT1 will be zoomed.
6. FunDMDeep-m6A
If users only hase a list of DmM genes and their DE information, then they can use fdmdeepm6A to identify functional DmM genes (FDmMGenes).
DMgene is a group of DmM genes, can be gene symbol or entrez gene ID.
dminfo <- system.file("extdata", "DMinfo_toy.xls", package="Funm6AViewer")
deinfo <- system.file("extdata", "DEinfo_toy.xls", package="Funm6AViewer")
dminfo <- read.table(dminfo, header = TRUE, stringsAsFactors = FALSE)
deinfo <- read.delim(deinfo, header = TRUE, stringsAsFactors = FALSE)
DMgene <- unique(dminfo$name)
descore <- getdescore(deinfo)
The datapath is the filepath where the required PPI data saved and it can be downloaded from http://180.208.58.66/Funm6AViewer/Download/Funm6AViewer_data.zip
datapath <- "F:/Funm6A_package/data"
permutime <- 1000
Identify FDmMGenes:
re <- fdmdeepm6A(DMgene = DMgene, descore = descore, datapath = datapath, permutime = permutime)
Plot interested genes' MSB score:
siggene <- c("CCNT1", "MYC", "BCL2")
siggenescoreplot(fdmgene = re, siggene = siggene)
7. Context specific function annotation of interested FDmMGenes
Users can visualize the context specific function of interested FDmMGenes identified by FunDMDeep-m6A using siggenepathplot.
fdmgene is a group of identified FDmMGnes. siggene is interested FDmMGne.
dminfo <- system.file("extdata", "DMinfo_toy.xls", package="Funm6AViewer")
dminfo <- read.table(dminfo, header = TRUE, stringsAsFactors = FALSE)
fdmgene <- unique(dminfo$name)
siggene <- c("MYC")
The input_directory is the filepath passed to string_db, the GO and KEGG function annotation data will be downloaded to this path. Users can also donwloaded the annotation data previously from http://180.208.58.66/Funm6AViewer/Download/Funm6AViewer_data.zip and set the input_directory as where you save the data.
input_directory <- "F:/Funm6A_package/data"
re <- siggenepathplot(fdmgene = fdmgene, intrested_gene = siggene,
version = "11", input_directory = input_directory)
8. MSB net plot for interested FDmMGenes
Users can visualize the MSB neighbours of interested FDmMGenes identified by FunDMDeep-m6A using msbnetplot.
dmgene is a group of DmM gnes. siggene is interested FDmMGnes.
dminfo <- system.file("extdata", "DMinfo_toy.xls", package="Funm6AViewer")
deinfo <- system.file("extdata", "DEinfo_toy.xls", package="Funm6AViewer")
dminfo <- read.table(dminfo, header = TRUE, stringsAsFactors = FALSE)
deinfo <- read.delim(deinfo, header = TRUE, stringsAsFactors = FALSE)
dmgene <- unique(dminfo$name)
descore <- getdescore(deinfo)
siggene <- c("CCNT1", "MYC", "BCL2")
The datapath is the filepath where the required PPI data saved and it can be downloaded from http://180.208.58.66/Funm6AViewer/Download/Funm6AViewer_data.zip
datapath <- "F:/Funm6A_package/data"
Plot for one FDmMGne:
re <- msbnetplot(genesymbol = siggene[1], dmgene = dmgene, descore = descore, datapath = datapath)
plot for several FDmMGnes:
re <- msbnetplot(genesymbol = siggene, dmgene = dmgene, descore = descore, datapath = datapath,
savename = "InterestedGene")
9. Functional m6A analysis pipline from bam files
We introduced this pipeline using a pretend example MeRIP-Seq data which contains 2 replicates in bam format for each condition (untreated and treated), named as:
"untreated_input_rep1.bam","untreated_ip_rep1.bam",
"untreated_input_rep2.bam", "untreated_ip_rep2.bam",
"treated_input_rep1.bam", "treated_ip_replicate1.bam",
"treated_input_rep2.bam", "treated_ip_rep2.bam".
These bam format files can be obtained by aligning the raw sequenced MeRIP-Seq data to genome, i.e., human hg19 using splice junction mapping tools for RNA-Seq reads, like TopHat, HISAT or STAR.
9.1 Calling differential m6A methylation sites using DMDeepm6A
Users can firstly call single base differential m6A methylation (DmM) sites using DMDeepm6A R package (https://github.com/NWPU-903PR/DMDeepm6A1.0) as following:
library(DMDeepm6A)
ip_bams <- c("treated_ip_rep1.bam"," treated_ip_rep2.bam",
"untreated_ip_rep1.bam ", "untreated_ip_rep2.bam")
input_bams <- c("treated_input_rep1.bam","treated_input_rep2.bam",
"untreated_input_rep1.bam", "untreated_input_rep2.bam")
sample_condition <- c("treated", "treated", "untreated", "untreated")
Assuming the aligned genome is hg19, users can call DmM sites as following:
output_filepath <- getwd()
re <- dmdeepm6A(ip_bams = ip_bams,
input_bams = input_bams,
sample_conditions = sample_condition,
output_filepath = output_filepath,
experiment_name = "DMDeepm6A_out")
See ?dmdeepm6A for more details if you are using other genomes. The experiment_name is the name of the file folder where saved the result of dmdeepm6A.
9.2 Making input for Funm6AViewer
Funm6AViewer takes single base DmM sites information, gene DE information as input and a list of GRanges converted from MeRIP-Seq bam files using makegrreadsfrombam if users would like to see the reads coverage of interested DmM sites on gene. Single base DmM sites information named dminfo contains the position annotation and log2 foldchange of DmM sites. It can be extracted from the result of DMDeepm6A using summarydmdeepm6A as following:
dminfo <- summarydmdeepm6A(dmpath = " DMDeepm6A_out", sigthresh = 0.05)
dminfo contains the following information:
head(dminfo)
## chr chromStart chromEnd name score strand log2fd
## 1 chr1 155160832 155160833 4582 0.9137714 - 2.7950754
## 2 chr1 171505224 171505225 23215 0.9431386 + 0.7589479
## 3 chr1 241767682 241767683 23596 0.8125095 - 1.0680801
## 4 chr1 243418399 243418400 9859 0.9652586 - 1.7805035
## 5 chr1 8073372 8073373 54206 0.8477583 - 2.5678589
## 6 chr1 8073689 8073690 54206 0.8089832 - 1.1375522
The ‘name’ column can be entrez gene ID or gene symbol.
Gene DE information named deinfo contains the differential expresion p-value and fdr for genes. It can be made using makegrreadsfrombam and getdeinfo from MeRIP-seq input samples as following:
ip_bams <- c("treated_ip_rep1.bam"," treated_ip_rep2.bam",
"untreated_ip_rep1.bam ", "untreated_ip_rep2.bam")
input_bams <- c("treated_input_rep1.bam","treated_input_rep2.bam",
"untreated_input_rep1.bam", "untreated_input_rep2.bam")
sample_condition <- c("treated", "treated", "untreated", "untreated")
savepath <- getwd()
grlist <- makegrreadsfrombam(IP_bams = ip_bams,
Input_bams = input_bams,
condition = sample_condition,
minimal_alignment_MAPQ = 30,
txdb = TxDb.Hsapiens.UCSC.hg19.knownGene,
savepath = savepath)
deinfo <- getdeinfo(grlist = grlist,
txdb = TxDb.Hsapiens.UCSC.hg19.knownGene,
savepath = savepath)
The converted grlist will be saved to savepath named as " bamgrlist.RData ". Users can directly load it to use next time. If you are using other genomes please install the corresponding txdb annotation similar to TxDb.Hsapiens.UCSC.hg19.knownGene of the genome and assign it to txdb.
deinfo contains the following information:
head(deinfo)
## name pval padj log2FoldChange
## 1 1 7.578860e-01 8.990406e-01 -0.07382034
## 2 100 6.958592e-01 8.695820e-01 0.07649361
## 3 1000 4.155368e-06 6.420489e-05 -0.73757946
## 4 10000 2.043250e-02 9.424864e-02 0.40075777
## 5 100009676 4.524888e-01 7.148682e-01 -0.16186352
## 6 10001 6.708161e-01 8.566831e-01 0.08459812
The ‘name’ column can be entrez gene ID or gene symbol.
9.3 One step usage of Funm6AViewer
Following is an example to achieve Functional m6A analysis using funm6aviewer with previously generated dminfo, deinfo and grlist:
siggene <- c("CCNT1", "MYC", "BCL2")
permutime <- 100*length(unique(dminfo$name))
The datapath is the file path where the required PPI data saved and the enrich_input_directory is the file path passed to string_db, the GO and KEGG function annotation data will be downloaded and saved to this path. All these required data can be downloaded from http://180.208.58.66/Funm6AViewer/Download/Funm6AViewer_data.zip
datapath <- "~./Funm6AViewer_data"
enrich_input_directory <- "~./Funm6AViewer_data"
savepath <- getwd()
re <- funm6aviewer(dminfo = dminfo,
deinfo = deinfo,
grlist = grlist,
intrested_gene = siggene,
permutime = permutime,
datapath = datapath,
enrich_input_directory = enrich_input_directory,
savepath = savepath)
The results will be saved to savepath.