小麦BSgenome构建

安装R包BSgenome

#R
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("BSgenome")

输入文件1：2bit格式的基因组序列文件

#shell
#从UCSC中下载faToTwoBit
wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/faToTwoBit
chmod 755 faToTwoBit

#下载小麦基因组序列文件
wget https://ftp.ensemblgenomes.ebi.ac.uk/pub/plants/release-56/fasta/triticum_aestivum/dna/Triticum_aestivum.IWGSC.dna.toplevel.fa.gz

#将fasta格式转为2bit格式
faToTwoBit Triticum_aestivum.IWGSC.dna.toplevel.fa.gz Triticum_aestivum.IWGSC.dna.toplevel.2bit

输入文件2：seed文件

BSgenome包内置了某些基因组的seed文件，存在~/anaconda3/envs/python3.8/lib/R/library/BSgenome/extdata/GentlemanLab/路径下，我们将其中拟南芥的seed文件复制一份，修改成小麦的

小麦seed文件展示

#shell
#cat Triticum_aestivum.IWGSC-seed
Package: BSgenome.Triticum.aestivum.IWGSC
Title: Full genome sequences for Triticum aestivum
Description: Full genome sequences for Triticum aestivum as provided by IWGSC and stored in Biostrings objects.
Version: 1.0.0
organism: Triticum aestivum
common_name: bread wheat
genome: TAIR10.1#这里要求的genome名是UCSC或者NCBI中的，但小麦的基因组是从EnsemblPlants下载的，因此用拟南芥的代替，不会影响后续使用
provider: EnsemblPlants
release_date: 2018/03/15#直接用拟南芥的，未做改动
source_url: https://ftp.ensemblgenomes.ebi.ac.uk/pub/plants/release-56/fasta/triticum_aestivum/dna/
organism_biocview: Triticum.aestivum
BSgenomeObjname: Triticum.aestivum.IWGSC
SrcDataFiles: Triticum_aestivum.IWGSC.dna.toplevel.fa.gz from https://ftp.ensemblgenomes.ebi.ac.uk/pub/plants/release-56/fasta/triticum_aestivum/dna/
PkgExamples: genome[["1A"]]#从EnsemblPlants下载的小麦基因组的染色体号为1A、1B、1D...，因此这里将原本拟南芥的genome[["1"]]改为genome[["1A"]]
seqs_srcdir: ~/genome/ChineseSpring/BSgenome#指定2bit文件路径
seqfile_name: Triticum_aestivum.IWGSC.dna.toplevel.2bit#指定2bit文件名

在R中构建小麦BSgenome

#R
library(BSgenome)
forgeBSgenomeDataPkg("./Triticum_aestivum.IWGSC-seed",seqs_srcdir=getwd(),destdir=getwd(),verbose=T)
#seqs_srcdir指定2bit文件路径，destdir指定输出文件路径，空值则代表当前路径

结束后会在当前路径中生成BSgenome.Triticum.aestivum.IWGSC目录，目录结构如下

#shell
BSgenome.Triticum.aestivum.IWGSC
├── DESCRIPTION
├── inst
│   └── extdata
│       └── single_sequences.2bit
├── man
│   └── package.Rd
├── NAMESPACE
└── R
    └── zzz.R

5 directories, 5 files

将BSgenome.Triticum.aestivum.IWGSC下的文件打包成R包

#shell
R CMD build BSgenome.Triticum.aestivum.IWGSC
R CMD check BSgenome.Triticum.aestivum.IWGSC_1.0.0.tar.gz
R CMD INSTALL BSgenome.Triticum.aestivum.IWGSC_1.0.0.tar.gz

打包完后就可以直接在R中调用了

#R
library("BSgenome.Triticum.aestivum.IWGSC")
#显示染色体号
seqnames(BSgenome.Triticum.aestivum.IWGSC)
#获取染色体长度信息
seqlengths(BSgenome.Triticum.aestivum.IWGSC)
#获取指定染色体的序列，因为R不识别数字开头的变量，所以$1A要写成$"1A"，而$Un则可以直接写成$Un
A1=BSgenome.Triticum.aestivum.IWGSC$"1A"
Un=BSgenome.Triticum.aestivum.IWGSC$Un