Bioconductor注释专题：BSgenome（85）

---

基于Biostrings的基因组数据包

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

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

1. source('https:///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./learn/bioconductor/lecture/E57YU/biostrings

---