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" 指定文件路径。