Bioconductor注释专题：BSgenome（85）

---

基于Biostrings的基因组数据包

Bioconductor项目提供了包含给定生物体全基因组序列的数据包。这些软件包称为基于Biostrings的基因组数据包，因为它们包含的序列存储在Biostrings软件包中定义的一些基本容器中，如DNAString，DNAStringSet或MaskedDNAString容器。不管它们包含的特定序列数据如何，所有的基于生物基因组的基因组数据包非常相似，可以以一致且简单的方式进行操作。为了正常工作，他们都需要使用BSgenome软件包。不同于基于Biostrings的基因组数据包，这个软件包与是一个提供支持它们所需基础的软件包（这就是为什么基于Biostrings的基因组数据包也称为BSgenome数据包）。 BSgenome软件包本身需要Biostrings软件包。

安装BSgenome包（如果没有安装）

1. source("https://bioconductor.org/biocLite.R")

2. biocLite("BSgenome.Hsapiens.UCSC.hg19")

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

1. library(BSgenome)

2. (ag <- available.genomes())

3. unique(gsub("BSgenome\\.([^\\.]+).*", "\\1", ag))

通过运行结果可得知，当前版本提供了24个物种基因组包

1. [1] "Alyrata"       "Amellifera"    "Athaliana"    

2. [4] "Btaurus"       "Celegans"      "Cfamiliaris"  

3. [7] "Dmelanogaster" "Drerio"        "Ecoli"        

4. [10] "Gaculeatus"    "Ggallus"       "Hsapiens"    

5. [13] "Mfascicularis" "Mfuro"         "Mmulatta"    

6. [16] "Mmusculus"     "Osativa"       "Ptroglodytes"

7. [19] "Rnorvegicus"   "Scerevisiae"   "Sscrofa"      

8. [22] "Tgondii"       "Tguttata"      "Vvinifera"  

获取Ecoli 的BSgenome数据包,bing载入BSgenome.Ecolide包，查看数据信息：

1. biocLite("BSgenome.Ecoli.NCBI.20080805")

2. library(BSgenome.Ecoli.NCBI.20080805)

结果如下:

1. [1] "BSgenome.Ecoli.NCBI.20080805"

2. [2] "Ecoli"

查看数据包的概况：

E. coli genome:

1. # organism: Escherichia coli (E. coli)

2. # provider: NCBI

3. # provider version: 2008/08/05 （发行日期）

4. # release date: NA

5. # release name: NA

6. # 13 sequences: （包含的序列）

7. #   NC_008253 NC_008563 NC_010468

8. #   NC_004431 NC_009801 NC_009800

9. #   NC_002655 NC_002695 NC_010498

10. #   NC_007946 NC_010473 NC_000913

11. #   AC_000091                    

12. # (use 'seqnames()' to see all the

13. # sequence names, use the '$' or '[['

14. # operator to access a given sequence)

获取每个染色体的名字还有其对应的长度：

1. > seqnames(Ecoli)

2. [1] "NC_008253" "NC_008563" "NC_010468"

3. [4] "NC_004431" "NC_009801" "NC_009800"

4. [7] "NC_002655" "NC_002695" "NC_010498"

5. [10] "NC_007946" "NC_010473" "NC_000913"

6. [13] "AC_000091"

7. > seqlengths(Ecoli)

8. NC_008253 NC_008563 NC_010468 NC_004431

9.  4938920   5082025   4746218   5231428

10. NC_009801 NC_009800 NC_002655 NC_002695

11.  4979619   4643538   5528445   5498450

12. NC_010498 NC_007946 NC_010473 NC_000913

13.  5068389   5065741   4686137   4639675

14. AC_000091

15.  4646332

查看其中一条染色体

1. Ecoli$NC_008253

2.  4938920-letter "DNAString" instance

3. seq: AGCTTTTCATTCTGAC...TTAGTAAGTGATTTTC

查看NC_008253序列中GC的数量

1. > letterFrequency(Ecoli$NC_008253,"GC")

2.    G|C

3. 2495020

查看NC_008253序列中GC的含量

1. > letterFrequency(Ecoli$NC_008253,"GC",as.prob=TRUE)

2.      G|C

3. 0.5051752

结合sapply函数统计碱基组成和GC含量

1. #统计碱基组成

2. > sapply(seqnames(Ecoli), function(x) alphabetFrequency(Ecoli[[x]]))

3.  NC_008253 NC_008563 NC_010468

4. A   1222723   1256126   1164516

5. C   1251581   1285309   1204681

6. G   1243439   1283517   1209555

7. T   1221177   1256945   1167466

8. M         0        11         0

9. R         0        34         0

10. W         0        25         0

11. S         0        17         0

12. Y         0        24         0

13. K         0        16         0

14. V         0         0         0

15. H         0         0         0

16. D         0         0         0

17. B         0         0         0

18. N         0         1         0

19. -         0         0         0

20. +         0         0         0

21. .         0         0         0

1. 统计GC含量

2. > sapply(seqnames(Ecoli), function(x) letterFrequency(Ecoli[[x]], letters = "GC",

3. +                                                         as.prob = TRUE))

4. NC_008253.G|C NC_008563.G|C NC_010468.G|C NC_004431.G|C NC_009801.G|C

5.    0.5051752     0.5054729     0.5086652     0.5047480     0.5062162

6. NC_009800.G|C NC_002655.G|C NC_002695.G|C NC_010498.G|C NC_007946.G|C

7.    0.5081961     0.5038297     0.5053706     0.5049668     0.5060419

8. NC_010473.G|C NC_000913.G|C AC_000091.G|C

9.    0.5078129     0.5078970     0.5079958

学习心得

虽然我对R也不是特别熟悉，跟着教程一个一个代码来跑，然后逐步理解也不难，希望大家也可以像我一样做到。最后也给大家推荐Coursera上讲解BSgenome专题的一个视频，帮助大家进一步理解这个包。文章链接：https://www.coursera.org/learn/bioconductor/lecture/E57YU/biostrings

---

猜你喜欢

还有更多文章，请移步公众号阅读

▼ 如果你生信基本技能已经入门，需要提高自己，请关注下面的生信技能树，看我们是如何完善生信技能，成为一个生信全栈工程师。

▼ 如果你是初学者，请关注下面的生信菜鸟团，了解生信基础名词，概念，扎实的打好基础，争取早日入门。