构建植物（棉花）BSgenome 参考基因组

构建非常见物种BSgenome参考基因组
1.前期工作，在linux下工作
1）将fa格式的文件转换成twoBit。BSgeome不识别fa的格式. 需要下载ucsc的faToTwoBit工具进行格式转化。
下载地址：http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/faToTwoBit

wget -c http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/faToTwoBit
faToTwoBit dm6.fa dm6.2bit
faToTwoBit genome.Ghir.HAU.fa genome.Ghir.HAU.2bit
1
2
3

2）获得每一条染色体的名字

samtools faidx genome.Ghir.HAU.fa |awk '{print $1}' > genome.Ghir.HAU.chromName.txt
# 或 
grep ">" genome.Ghir.HAU.fa|awk '{print $1}' | sed 's/^>//g' > genome.Ghir.HAU.chromName.txt
1
2
3

3.安装 BSgenome 包
if (!require(“BiocManager”, quietly = TRUE)) install.packages(“BiocManager”)
BiocManager::install(“BSgenome”)

4.进入工作目录，准备 seed 文件。
seed 文件可以理解为一种配置文件或者说明文件，包含许多必要的和非必要的词条，详细内容参考Bioconductor - BSgenome。在 BSgenome 包安装目录有示例文件，可以直接复制到工作目录进行修改。
cd /biodata/02.software/R-4.1.2/library/BSgenome/extdata/GentlemanLab

5.下面是seed 文件的一些bug及注意事项。
首先是 Package，必须由四部分组成，每个部分必须由逗号（.）分开。part 1：BSgenome；part 2：物种缩写，如 ovis aries 写成Oaries；part 3：参考基因组来源，一般是NCBI，Ensembl 或者 UCSC；part 4：参考基因组版本。必须要注意的是！！！在package词条不能出现多余的逗号（.）和其他的分隔符（-或_）。比如 part 4中，最容易理解的方式是写成 rambouillet.v1，rambouillet_v1 或 rambouillet-v1，然而写成这样会报错，报错内容是“malformed package name”！！！

第二个 bug 是在 seed 文件中必须提供 genome 词条，在这里我用的是 genome: bosTau9，是的，这是一个牛的参考基因组名称。第一，BSgenome 只支持 UCSC 和 NCBI 数据库下载的参考基因组，而这个绵羊的参考基因组是在 Ensembl 下载的；第二，在 UCSC 和 NCBI 找不到可用的、同版本的参考基因组，因此，没有办法，只能随便找一个代替。

第三个bug是seqnames 构建包时在R环境中可行，但是安装的时候会出现找不到该对象，解决的办法是如下写法（注：必须为字符向量）：

seqnames: as.character(read.table("/biodata/05.standard_analysis/cotton_AJFX2220520037_scATAC_20220530/new_analysis/Ghir_HAU/genome.Ghir.HAU.chromName.txt")$V1)
1

下面是我的seed文件供参考

Package: BSgenome.Ghau.Ensembl.GhirHAUv1
Title: Full genome sequences for Ghir HAU(Ensembl version GhirHAU_v1)
Description: Full genome sequences for Ghir HAU(cottom) as provided by Ensembl (GhirHAU_v1) and stored in Biostrings objects.
Version: 1.0.0
organism: Ghir HAU
common_name: cottom
genome: TAIR10
provider: Ensembl
release_date: July 2022
release_name: IGENEBOOK
source_url: http://hgdownload.cse.ucsc.edu/goldenPath/rn4/bigZips/
organism_biocview: Ghir_HAU
BSgenomeObjname: HAU
#seqnames: paste("chr", c(1:20, "X", "M", "Un", paste(c(1:20, "X", "Un"), "_random", sep="")),sep="")
#seqnames: set the sequence names in R first: seqnames_all
#seqnames: seqnames1
seqnames: as.character(read.table("/biodata/05.standard_analysis/cotton_AJFX2220520037_scATAC_20220530/new_analysis/Ghir_HAU/genome.Ghir.HAU.chromName.txt")$V1)
# circ_seqs: "dmel_mitochondrion_genome"
# SrcDataFiles: chromFa.tar.gz from http://hgdownload.cse.ucsc.edu/goldenPath/rn4/bigZips/
SrcDataFiles: genome.Ghir.HAU.2bit, transferred by faToTwoBit
# seqs_srcdir: /fh/fast/morgan_m/BioC/BSgenomeForge/srcdata/BSgenome.Rnorvegicus.UCSC.rn4/seqs
#seqfile_name: the full of sequence name
seqfile_name: genome.Ghir.HAU.2bit
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

6.构建R包，在R中的工作

library(BSgenome)
#读取前面构建的染色体名字的文件，注意 seqnames应该和seed文件中的关键字一致

seqnames=read.table("dm6.chromName.txt")
seqnames=as.character(seqnames$V1)
#染色体名有<211000022280427>,这种纯数字的，这种chr是不能识别的
seqnames1=seqnames[grep("^211",seqnames,invert=T)]

#给出seed文件
dm6_seed="BSgenome.Dmelanogaster.Ensemble.dm6-seed"

#seqs_srcdir;destdir 序列文件所在以及输出的位置
forgeBSgenomeDataPkg(dm6_seed, seqs_srcdir=getwd(), destdir=getwd(), verbose=TRUE)
forgeBSgenomeDataPkg(dm6_seed, verbose=TRUE)
1
2
3
4
5
6
7
8
9
10
11
12
13
14

7.Linux 系统下运行如下命令

tree BSgenome.Ghau.Ensembl.GhirHAUv1
1

BSgenome.Ghau.Ensembl.GhirHAUv1
├── DESCRIPTION
├── inst
│ └── extdata
│ └── single_sequences.2bit
├── man
│ └── package.Rd
├── NAMESPACE
└── R
└── zzz.R

BSgenome.Oaries.Ensembl.rambouilletv1 同级目录运行

R CMD build ./BSgenome.Ghau.Ensembl.GhirHAUv1/  #生成BSgenome.Ghau.Ensembl.GhirHAUv1_1.0.0.tar.gz 文件
#如果发现Bug查找问题及时修复，只要 check for installation 不报错，基本就没有太大问题
R CMD check BSgenome.Ghau.Ensembl.GhirHAUv1_1.0.0.tar.gz
R CMD INSTALL BSgenome.Ghau.Ensembl.GhirHAUv1_1.0.0.tar.gz
1
2
3
4

8.参考文章
https://blog.csdn.net/flashan_shensanceng/article/details/125155320