R/BioC序列处理之四：BSgenome简介

一、BSgenome和BSgenome数据包

Bioconductor提供了某些物种的全基因组序列数据包，这些数据包是基于Biostrings构建的，称为BSgenome数据包。不同物种的BSgenome数据包都有类似的数据结构，可以用统一的方式进行处理。但是BSgenome数据包仅包含有数据，它们的处理的方法由另外一个软件包提供，即BSgenome包。
先安装BSgenome包（如果没有安装）：

library(BiocInstaller)
biocLite("BSgenome")

载入BSgenome包，并查看当前版本提供的BSgenome数据包：

library(BSgenome)
(ag <- available.genomes())

##  [1] "BSgenome.Alyrata.JGI.v1"
##  [2] "BSgenome.Amellifera.BeeBase.assembly4"
##  [3] "BSgenome.Amellifera.UCSC.apiMel2"
##  [4] "BSgenome.Athaliana.TAIR.04232008"
##  [5] "BSgenome.Athaliana.TAIR.TAIR9"
##  [6] "BSgenome.Btaurus.UCSC.bosTau3"
##  [7] "BSgenome.Btaurus.UCSC.bosTau4"
##  [8] "BSgenome.Btaurus.UCSC.bosTau6"
##  [9] "BSgenome.Celegans.UCSC.ce10"
## [10] "BSgenome.Celegans.UCSC.ce2"
## [11] "BSgenome.Celegans.UCSC.ce6"
## [12] "BSgenome.Cfamiliaris.UCSC.canFam2"
## [13] "BSgenome.Cfamiliaris.UCSC.canFam3"
## [14] "BSgenome.Dmelanogaster.UCSC.dm2"
## [15] "BSgenome.Dmelanogaster.UCSC.dm3"
## [16] "BSgenome.Drerio.UCSC.danRer5"
## [17] "BSgenome.Drerio.UCSC.danRer6"
## [18] "BSgenome.Drerio.UCSC.danRer7"
## [19] "BSgenome.Ecoli.NCBI.20080805"
## [20] "BSgenome.Gaculeatus.UCSC.gasAcu1"
## [21] "BSgenome.Ggallus.UCSC.galGal3"
## [22] "BSgenome.Ggallus.UCSC.galGal4"
## [23] "BSgenome.Hsapiens.UCSC.hg17"
## [24] "BSgenome.Hsapiens.UCSC.hg18"
## [25] "BSgenome.Hsapiens.UCSC.hg19"
## [26] "BSgenome.Mmulatta.UCSC.rheMac2"
## [27] "BSgenome.Mmusculus.UCSC.mm10"
## [28] "BSgenome.Mmusculus.UCSC.mm8"
## [29] "BSgenome.Mmusculus.UCSC.mm9"
## [30] "BSgenome.Ptroglodytes.UCSC.panTro2"
## [31] "BSgenome.Ptroglodytes.UCSC.panTro3"
## [32] "BSgenome.Rnorvegicus.UCSC.rn4"
## [33] "BSgenome.Rnorvegicus.UCSC.rn5"
## [34] "BSgenome.Scerevisiae.UCSC.sacCer1"
## [35] "BSgenome.Scerevisiae.UCSC.sacCer2"
## [36] "BSgenome.Scerevisiae.UCSC.sacCer3"
## [37] "BSgenome.Tgondii.ToxoDB.7.0"

# 当前版本提供了18个物种的37个基因组包
unique(gsub("BSgenome\\.([^\\.]+).*", "\\1", ag))

##  [1] "Alyrata"       "Amellifera"    "Athaliana"     "Btaurus"
##  [5] "Celegans"      "Cfamiliaris"   "Dmelanogaster" "Drerio"
##  [9] "Ecoli"         "Gaculeatus"    "Ggallus"       "Hsapiens"
## [13] "Mmulatta"      "Mmusculus"     "Ptroglodytes"  "Rnorvegicus"
## [17] "Scerevisiae"   "Tgondii"

获取拟南芥的BSgenome数据包（BSgenome.Athaliana.TAIR.TAIR9）：

biocLite("BSgenome.Athaliana.TAIR.TAIR9")

载入BSgenome.Athaliana.TAIR.TAIR9包，查看数据信息：

library(BSgenome.Athaliana.TAIR.TAIR9)
# BSgenome.Athaliana.TAIR.TAIR9只定义了一个对象
ls("package:BSgenome.Athaliana.TAIR.TAIR9")

## [1] "Athaliana"

# 查看总体信息
Athaliana

## Arabidopsis genome
## |
## | organism: Arabidopsis thaliana (Arabidopsis)
## | provider: TAIR
## | provider version: TAIR9
## | release date: June 9, 2009
## | release name: TAIR9 Genome Release
## |
## | sequences (see '?seqnames'):
## |   Chr1  Chr2  Chr3  Chr4  Chr5  ChrM  ChrC
## |
## | (use the '$' or '[[' operator to access a given sequence)

# 不载入序列就可以获取染色体名称和长度信息
seqnames(Athaliana)

## [1] "Chr1" "Chr2" "Chr3" "Chr4" "Chr5" "ChrM" "ChrC"

seqlengths(Athaliana)

##     Chr1     Chr2     Chr3     Chr4     Chr5     ChrM     ChrC
## 30427671 19698289 23459830 18585056 26975502   366924   154478

二、BSgenome方法

BSgenome数据包内包含的序列为DNAString，DNAStringSet或MaskedDNAString对象，完全可以使用Biostrings包的方法进行处理：

# 可以使用多种方式获取某一序列
Athaliana$Chr1

##   30427671-letter "DNAString" instance
## seq: CCCTAAACCCTAAACCCTAAACCCTAAACCTCTG...TAGGGTTTAGGGTTTAGGGTTTAGGGTTTAGGG

Athaliana[[1]]

##   30427671-letter "DNAString" instance
## seq: CCCTAAACCCTAAACCCTAAACCCTAAACCTCTG...TAGGGTTTAGGGTTTAGGGTTTAGGGTTTAGGG

Athaliana[["Chr1"]]

##   30427671-letter "DNAString" instance
## seq: CCCTAAACCCTAAACCCTAAACCCTAAACCTCTG...TAGGGTTTAGGGTTTAGGGTTTAGGGTTTAGGG

下面使用sapply函数统计碱基组成和GC含量。如果不熟悉sapply函数，请参考这篇文章。

# 统计碱基组成
sapply(seqnames(Athaliana), function(x) alphabetFrequency(Athaliana[[x]]))

##      Chr1    Chr2    Chr3    Chr4    Chr5   ChrM  ChrC
## A 9709674 6315641 7484757 5940546 8621974 102464 48546
## C 5435374 3542973 4258333 3371349 4832253  82661 28496
## G 5421151 3520766 4262704 3356091 4858759  81609 27570
## T 9697113 6316348 7448059 5914038 8652238 100190 49866
## M      76       5       2       1       0      0     0
## R      36       7       4       0       0      0     0
## W     124      18       2       0       0      0     0
## S      30       3       1       0       0      0     0
## Y      82      12       2       0       0      0     0
## K      53      10       0       0       0      0     0
## V       0       0       0       0       0      0     0
## H       0       0       0       0       0      0     0
## D       0       0       0       1       0      0     0
## B       0       0       0       0       0      0     0
## N  163958    2506    5966    3030   10278      0     0
## -       0       0       0       0       0      0     0
## +       0       0       0       0       0      0     0

# 计算GC含量
sapply(seqnames(Athaliana), function(x) letterFrequency(Athaliana[[x]], letters = "GC", as.prob = TRUE))

## Chr1.G|C Chr2.G|C Chr3.G|C Chr4.G|C Chr5.G|C ChrM.G|C ChrC.G|C
##   0.3568   0.3586   0.3632   0.3620   0.3593   0.4477   0.3629

但是BSgenome包试图提供一些简便的方法来处理基因组水平的BSgenome类数据，例如类似于apply函数的bsapply函数。要使用bsapply函数得先构建BSParams类对象，用于设置以下参数：

X：将要处理的BSgenome对象
FUN：将要对X中每条染色体进行处理函数
exclude：字符向量，表示排除的染色体名称
simplify：逻辑向量，它的意义和与sapply的simplify参数一样。默认为FALSE，bsapply返回GenomeData类数据；如果设置为TRUE，函数的结果尽量使用表格（数据框）类型显示
maskList：逻辑向量，表示各染色体使用掩膜的应用状态
motifList：字符向量，对序列进行掩膜的motif
userMask：RangesList对象，每个元素将用于对响应染色体进行掩膜处理，为用户提供额外的掩膜。
invertUserMask：TRUE/FALSE，表示是否对userMask进行反转。

FUN函数的其他参数可以在bsapply函数中以有名参数的方式进行设置。下面代码的作用和前面sapply函数的应用是相同的：

op <- options()
options(digits = 4)
params <- new("BSParams", X = Athaliana, FUN = letterFrequency, simplify = TRUE)
bsapply(params, letters = "GC", as.prob = TRUE)

## Chr1.G|C Chr2.G|C Chr3.G|C Chr4.G|C Chr5.G|C ChrM.G|C ChrC.G|C
##   0.3568   0.3586   0.3632   0.3620   0.3593   0.4477   0.3629

BSParams是S4类，对象数据可以修改，方便重复使用：

params@motifList <- "N"
bsapply(params, letters = "GC", as.prob = TRUE)

## Chr1.G|C Chr2.G|C Chr3.G|C Chr4.G|C Chr5.G|C ChrM.G|C ChrC.G|C
##   0.3587   0.3586   0.3633   0.3620   0.3594   0.4477   0.3629

options(op)
params@FUN <- alphabetFrequency
bsapply(params, baseOnly = TRUE)

##          Chr1    Chr2    Chr3    Chr4    Chr5   ChrM  ChrC
## A     9709674 6315641 7484757 5940546 8621974 102464 48546
## C     5435374 3542973 4258333 3371349 4832253  82661 28496
## G     5421151 3520766 4262704 3356091 4858759  81609 27570
## T     9697113 6316348 7448059 5914038 8652238 100190 49866
## other     401      55      11       2       0      0     0

BSgenome包实现了BSgenome类对象的matchPWM，countPWM，vmatchPattern，vcountPattern，vmatchPDict和vcountPDict等操作，除继承了Biostrings中相应的方法外还增加了一些针对BSgenome类对象的参数设置（和bsapply函数参数类似）：

vmatchPattern("ACGTGKC", Athaliana, fix = "subject", exclude = c("ChrC", "ChrM"))

## GRanges with 13745 ranges and 0 metadata columns:
##           seqnames               ranges strand
##              ##       [1]     Chr1       [ 1540,  1546]      +
##       [2]     Chr1       [ 8276,  8282]      +
##       [3]     Chr1       [12017, 12023]      +
##       [4]     Chr1       [14697, 14703]      +
##       [5]     Chr1       [34225, 34231]      +
##       ...      ...                  ...    ...
##   [13741]     Chr5 [26911361, 26911367]      -
##   [13742]     Chr5 [26915247, 26915253]      -
##   [13743]     Chr5 [26917155, 26917161]      -
##   [13744]     Chr5 [26923565, 26923571]      -
##   [13745]     Chr5 [26933617, 26933623]      -
##   ---
##   seqlengths:
##        Chr1     Chr2     Chr3     Chr4     Chr5     ChrM     ChrC
##    30427671 19698289 23459830 18585056 26975502   366924   154478

BSgenome包还提供了其他一些方法，比如用户自行构建BSgenome数据包和SNP的处理等。这些不在我的关心范围之内，没去了解。