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 <- alphabetFrequencybsapply(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的处理等。这些不在我的关心范围之内，没去了解。