1. 目的
最近在某个分析的时候,需要BSgenome object,然而我的物种不是常见的物种,是自己组装出来的多倍体物种,因此没法用已有的BSgenome参考基因组,需要自己构建BSgenome参考基因组,通过BSgenome相关介绍说明的时候,发现有一个How to forge a BSgenome data package说明,这说明了我们能够自己通过BSgenome包,构建自己的参考基因组,下面就是构建BSgenome参考基因组的过程
2. 输入fasta的准备
构建BSgenome参考基因组的方法有两种,一种是不需要masked sequences,另外一种是需要masked sequences文件,这里我选择第一种进行构建。其实两者最主要的差别在于是否提供masked sequences文件。
构建BSgenome参考基因组具体需要什么文件呢?说明文档中也有详细的说明,其中最主要的文件就是一个fasta文件,这里的fasta文件可以是所有染色体合并一起的一个文件,也可以是单条染色体的fasta文件,文件格式可以是压缩文件格式,后续的处理几乎都是在fasta文件的基础上进行处理。
前面提到了fasta文件需要进行处理,究竟如何处理呢?由于输入的fasta文件是TwoBit文件,因此需要下载faToTwoBit工具,对fasta文件进行转换,具体命令如下:
wget -c "http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/faToTwoBit"
chmod 755 faToTwoBit
faToTwoBit gene.fa gene.2bit
这里faToTwoBit运行需要一定的时间,主要是与参考基因组的大小相关,生成的是二进制文件,生成twoBit后,就可以进行下一步的准备了。
3. seed文件的准备
构建BSgenome参考基因组另外一个需要准备的文件就是seed文件,其实这里的seed文件就是在开发R包的DESCRIPTION文件, 这里相当于构建了一个R包,R包需要添加一下DESCRIPTION信息,因此seed的文件,可以参考DESCRIPTION的格式说明,具体格式要求如下:
上图可以看出,seed文件,主要包括Package/Title/Description/Version/Author/Maintainer/License,其他的选填项可以不用填写。seed文件格式有很多种,可根据实际需要进行撰写;另外,也可以根据已经安装的BSgenome包来查看seed格式,具体的查看命令如下:
library(BSgenome)
seed_files <- system.file("extdata", "GentlemanLab", package="BSgenome")
tail(list.files(seed_files, pattern="-seed$"))
#[1] "BSgenome.Sscrofa.UCSC.susScr11-seed" "BSgenome.Sscrofa.UCSC.susScr3-seed" #"BSgenome.Sscrofa.UCSC.susScr3.masked-seed" "BSgenome.Tguttata.UCSC.taeGut1-seed"
#这里可以看到有很多seed文件,可以随便cp一个文件,然后在此基础上就行即可。
seed文件其他的一些可选项也有很多,比如这里的organism、common_name、genome、provider、release_date、source_url、organism_biocview等等,这些是为了更详细的记录一些参考基因组信息,一遍后续使用或者构建的时候了解相关信息,其他的一些说明,这里不再进行详细的介绍,可以直接通过How to forge a BSgenome data package进行相关信息的查询。
4. forgeBSgenomeDataPkg
seed文件生成以后,就可以直接构建BSgenome包了,这里构建BSgenome包的具体命令如下:
seed_file="gene.物种.seed"
#seqs_srcdir;destdir 序列文件所在以及输出的位置
forgeBSgenomeDataPkg(seed_file, seqs_srcdir=getwd(), destdir=getwd(), verbose=TRUE,unlink=TRUE)
#forgeBSgenomeDataPkg(seed_file, verbose=TRUE)
#unlink参数表示是否overwrite已有的目录,seqs_srcdir是twoBit的目录,destdir为生成包的目录,这里需要一定的时候。
5. 包的构建以及安装
前面forgeBSgenomeDataPkg运行成功以后,会在目的目录中生成R包的几个基础文件,如下图所示:
这是一个基本的R包的结果目录,然后在此基础上,在linux命令进行包的build、check、install,具体的命令如下:
R CMD build <pkgdir>
R CMD check <tarball>
R CMD INSTALL <tarball>
一般来说,build没有问题的话,就不用再check了,我这边构建的包后进行check,结果报错,报错的是LaTeX有问题,然后尝试不管他,再直接安装,结果能够正常安装,并且正常导入包。因此这里的报错可以不用考虑,这应该是要生成pdf的说明文件报错,对后续其他的结果没有影响。
6. 后记
这里构建BSgenome参考基因组是选择的最简单的方式,没有提供masked sequences文件,如果需要masked sequences文件,那么可以直接参考方说明文档。可以通过此方法,构建任何自己的参考基因组,进行其他相关的分析。
参考文档
2021年1月11日