exomePeak2 使用笔记之自建 BSgenome Object
最近需要使用 exomePeak2 这个包来对玉米的 RNA-seq 和 RIP-seq 数据进行 peak calling，但其中的 bsgenome 参数需要导入对应物种的 BSgenome object，所以需要自主构建参考基因组的 BSgenome data package。
BSgenome 数据包
BSgenome 数据包是基于 Biostrings 构建的全基因组数据序列数据包。BSgenome 数据包中包含物种的数据并且具有相似的数据结构，可以用统一的方式进行处理，由 BSgenome 软件包提供。
首先安装 BSgenome 包。
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("BSgenome")
然后导入包，并查看官方已有的 BSgenome 数据包
library(BSgenome)
ag <- available.genomes()
但当前版本中没有玉米的数据包，所以需要自己构建。
1. 构建玉米的 BSgenome 数据包
1) 首先将 fa 文件转换为 2bit 文件。 此时需要用到 faToTwoBit。
faToTwoBit zmB73.fa zmB73.2bit
2) 提取玉米的染色体名称
less -S zmB73.fa | grep ">" |awk '{print $1}' | sed 's/^>//g' > zmB73.chromName.txt
3) 书写玉米的 seed 文件 ZmaysZm4_seed。此处我使用的参考基因组版本是 Zea Mays RefGen_v4 （AGPv4）。
根据 BSgenome 的参考文件 https://www.bioconductor.org/packages/devel/bioc/vignettes/BSgenome/inst/doc/BSgenomeForge.pdf 第二章： 2.2 Prepare the BSgenome data package seed file 提示编写。
Package: BSgenome.Zmays.Ensemble.zmv4
Title: Genome sequences for Zea mays (Ensemble AGPv4)
Description: A BSgenome package containing the full genome sequences for Zea mays (Maize) as provided by Ensemble (B73 AGPv4, Sept. 2020) and stored in Biostrings objects.
Version: 1.0
organism: Zea mays
common_name: Maize
provider: Ensemble
provider_version: zmv4
release_date: Sept. 2020
release_name: Maize Genome Sequencing B73 RefGen_v4.0
#source_url: ftp://ftp.ensemblgenomes.org/pub/plants/release-49/fasta/zea_mays/dna/
circ_seqs: c("Mt","Pt")
organism_biocview: Zea_mays
BSgenomeObjname: Zmays
seqs_srcdir: /data/zhanmin/preparation/index/
seqfiles_suffix: .fa
seqnames: seqnames1 # or you can use the R code directly: as.character(read.table(" zmB73.chromName.txt ")$V1)
seqfile_name: zmB73.2bit
4) 在 R 中构建 R 包
library(rtracklayer)
library(BSgenome)
## import the sequence names of maize
## If you already use the R code in seed file to read the sequence names, this part code is not necessary.
seqnames=read.table("zmB73.chromName.txt")
seqnames1=as.character(seqnames$V1)
# import seed file
zmB73_seed="BSgenome.Zmays.Ensemble.zm4-seed"
# make the BSgenome data package for maize. seqs_srcdir: the dirctory of sequence file; destdir: the directory of output
forgeBSgenomeDataPkg(zmB73_seed, seqs_srcdir=getwd(), destdir=getwd(), verbose=TRUE)
运行完毕后会在指定目录或工作目录中产生输出文件夹：BSgenome.Zmays.Ensemble.zmv4
2. 在 linux 中构建 R 包
1） 使用 build 命令构建 R 包
R CMD build BSgenome.Zmays.Ensemble.zmv4 #运行后会产生相应的 tar.gz 压缩文件
此时发生报错：
ERROR
cannot change to directory ‘BSgenome.Zmays.Ensemble.zmv4’
(python3.7) user@amax:/data/user/preparation/index/BSgenome.Zmays.Ensemble.zmv4$ cd ..
(python3.7) user@amax:/data/user/preparation/index$ R CMD build BSgenome.Zmays.Ensemble.zmv4
* checking for file ‘BSgenome.Zmays.Ensemble.zmv4/DESCRIPTION’ ... OK
* preparing ‘BSgenome.Zmays.Ensemble.zmv4’:
* checking DESCRIPTION meta-information ... ERROR
Malformed package version.
See section 'The DESCRIPTION file' in the 'Writing R Extensions'
manual.
根据提示查阅 Writing R Extensions 文档的 The DESCRIPTION file 章节，并对比自己的 description 文件。
发现是 Version 这个参数写错了，规定的范围值应该是 0.5-1，但这里我自己的文档写了 41。 这里应该直接从 seed 文件修改。
先把 seed 文件的 Version 改为 1.0， 然后重复 1. 4的步骤，重新构建 R 包。再次执行 R CMD build 命令，此时就成功了。
2） 用 check 命令查看是否产生 tar.gz 压缩文件
R CMD check BSgenome.Zmays.Ensemble.zmv4_1.0.tar.gz
3） 导入 R 包
R CMD INSTALL BSgenome.Zmays.Ensemble.zmv4_1.0.tar.gz
exomePeak2 进行 peak calling
1. 导入所需要的 R 包
没有的包需要先用 BioManager::install() 命令安装。
library(AnnotationHub) # help to build the annotation object
library(biomaRt)
library(GenomicFeatures)
library(exomePeak2)
library(rtracklayer)
library(BSgenome)
library(BSgenome.Zmays.Ensemble.zmv4) # import the Maize BSgenome data pacakage
构建玉米的 TxDb object 数据。
maize_txdb<- makeTxDbFromBiomart(biomart = "plants_mart",dataset = "zmays_eg_gene",host = "http://plants.ensembl.org")
saveDb(maize_txdb, file="maize_v4.sqlite")
maize_txdb <- loadDb("maize_v4.sqlite")
运行 exomePeak2 函数
exomePeak2(bam_ip = c("IP_1.bam","IP_2.bam"),
bam_input = c("Input_1.bam","Input_2.bam"),
txdb = maize_txdb,
bsgenome = BSgenome.Zmays.Ensemble.zmv4,
paired_end = TRUE)
过程大概需要几个小时，视数据大小而定。结束后会自动生成输出文件夹 exomePeak2_output，也可通过参数 save_dir = "exomePeak2_output" 指定文件路径。
