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

2014/01/26评论4,479

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

原文来自：http://blog.csdn.net/u014801157/article/details/24372467

一、BSgenome和BSgenome数据包

二、BSgenome方法

发表评论