exomePeak2 使用筆記之自建 BSgenome Object
最近需要使用 exomePeak2 這個(gè)包來(lái)對(duì)玉米的 RNA-seq 和 RIP-seq 數(shù)據(jù)進(jìn)行 peak calling，但其中的 bsgenome 參數(shù)需要導(dǎo)入對(duì)應(yīng)物種的 BSgenome object，所以需要自主構(gòu)建參考基因組的 BSgenome data package。
BSgenome 數(shù)據(jù)包
BSgenome 數(shù)據(jù)包是基于 Biostrings 構(gòu)建的全基因組數(shù)據(jù)序列數(shù)據(jù)包。BSgenome 數(shù)據(jù)包中包含物種的數(shù)據(jù)并且具有相似的數(shù)據(jù)結(jié)構(gòu)，可以用統(tǒng)一的方式進(jìn)行處理，由 BSgenome 軟件包提供。
首先安裝 BSgenome 包。
if (!requireNamespace("BiocManager", quietly = TRUE))?
? install.packages("BiocManager")?
BiocManager::install("BSgenome")??
然后導(dǎo)入包，并查看官方已有的 BSgenome 數(shù)據(jù)包
library(BSgenome)
ag <- available.genomes()
但當(dāng)前版本中沒(méi)有玉米的數(shù)據(jù)包，所以需要自己構(gòu)建。
1. 構(gòu)建玉米的 BSgenome 數(shù)據(jù)包
1) 首先將 fa 文件轉(zhuǎn)換為 2bit 文件。 此時(shí)需要用到 faToTwoBit。
faToTwoBit zmB73.fa zmB73.2bit
2) 提取玉米的染色體名稱(chēng)
less -S zmB73.fa | grep ">" |awk '{print $1}' | sed 's/^>//g' > zmB73.chromName.txt
3) 書(shū)寫(xiě)玉米的 seed 文件 ZmaysZm4_seed。此處我使用的參考基因組版本是 Zea Mays RefGen_v4 （AGPv4）。
根據(jù) BSgenome 的參考文件?https://www.bioconductor.org/packages/devel/bioc/vignettes/BSgenome/inst/doc/BSgenomeForge.pdf?第二章： 2.2?Prepare the BSgenome data package seed file 提示編寫(xiě)。
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 中構(gòu)建 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)
運(yùn)行完畢后會(huì)在指定目錄或工作目錄中產(chǎn)生輸出文件夾：BSgenome.Zmays.Ensemble.zmv4
2. 在 linux 中構(gòu)建 R 包
1） 使用 build 命令構(gòu)建 R 包
R CMD build BSgenome.Zmays.Ensemble.zmv4 #運(yùn)行后會(huì)產(chǎn)生相應(yīng)的 tar.gz 壓縮文件
此時(shí)發(fā)生報(bào)錯(cuò)：
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.
根據(jù)提示查閱?Writing R Extensions?文檔的?The DESCRIPTION file 章節(jié)，并對(duì)比自己的 description 文件。
發(fā)現(xiàn)是 Version 這個(gè)參數(shù)寫(xiě)錯(cuò)了，規(guī)定的范圍值應(yīng)該是 0.5-1，但這里我自己的文檔寫(xiě)了 41。 這里應(yīng)該直接從 seed 文件修改。
先把 seed 文件的 Version 改為 1.0， 然后重復(fù) 1. 4的步驟，重新構(gòu)建 R 包。再次執(zhí)行 R CMD build 命令，此時(shí)就成功了。
2） 用 check 命令查看是否產(chǎn)生 tar.gz 壓縮文件
R CMD check?BSgenome.Zmays.Ensemble.zmv4_1.0.tar.gz
3） 導(dǎo)入 R 包
R CMD INSTALL BSgenome.Zmays.Ensemble.zmv4_1.0.tar.gz
exomePeak2 進(jìn)行 peak calling
1. 導(dǎo)入所需要的 R 包
沒(méi)有的包需要先用 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
構(gòu)建玉米的 TxDb object 數(shù)據(jù)。
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")
運(yùn)行 exomePeak2 函數(shù)
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)
過(guò)程大概需要幾個(gè)小時(shí)，視數(shù)據(jù)大小而定。結(jié)束后會(huì)自動(dòng)生成輸出文件夾 exomePeak2_output，也可通過(guò)參數(shù)?save_dir = "exomePeak2_output" 指定文件路徑。
