【生信进阶练习1000days】day9-BSgenome和AnnotationHub

最新推荐文章于 2024-12-24 09:06:01 发布

原创最新推荐文章于 2024-12-24 09:06:01 发布 · 1.3k 阅读

1 ·

CC 4.0 BY-SA版权

2.3 生信进阶练习1000days 专栏收录该内容

14 篇文章

订阅专栏

博客介绍了BSgenome包和AnnotationHub。BSgenome包包含物种/基因组版本序列信息，可获取序列数据。AnnotationHub是R包和通道，能不安装注释包就下载和查询注释对象，还介绍了查询方法、相关练习及使用时的报错解决参考网址。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

文章目录

1. BSgenome packages

BSgenome包包含了给定物种/基因组版本的序列信息。我们可以通过使用 available.genomes函数，获取可用的BSgenome系列包


library(BSgenome)
head(available.genomes())
#> [1] "BSgenome.Alyrata.JGI.v1"                 "BSgenome.Amellifera.BeeBase.assembly4"  
#> [3] "BSgenome.Amellifera.UCSC.apiMel2"        "BSgenome.Amellifera.UCSC.apiMel2.masked"
#> [5] "BSgenome.Athaliana.TAIR.04232008"        "BSgenome.Athaliana.TAIR.TAIR9"

加载并检查BSgenome包

library(BSgenome.Hsapiens.UCSC.hg19)
Hsapiens

使用 getSeq() 函数获取序列数据，或者通过传入GRanges对象来获取某个区域的数据

getSeq(Hsapiens, "chr1")
#>   249250621-letter "DNAString" instance
#> seq: NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN...NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
getSeq(Hsapiens, gns["5467",])
#>   A DNAStringSet instance of length 1
#>     width seq                                                                                       names               
#> [1] 85634 GCGGAGCGTGTGACGCTGCGGCCGCCGCGGACCTGGGGATTAA...ACTTTAAATAAATGGGAATTAAATATTTAAGAGCTGACTGGAA 5467

Biostrings包包含处理这些 *StringSet 对象的大部分代码—可以参阅Biostrings小片段和帮助页面以获得更多信息。

BSgenome 练习 Get the sequences for all transcripts of the TP53 gene

2. AnnotationHub

AnnotationHub是一个R包，也可以视作一个通道，我们可以在不安装相应注释包的情况下，通过它下载和查询多种不同的注释对象（annotation objects）

library(AnnotationHub)
hub <- AnnotationHub()
#> snapshotDate(): 2018-06-27
hub

2.1 查询AnnotationHub

想在AnnotationHub上找到“正确”的资源，就像使用谷歌一样，我们需要使用一个经过良好配置的查询检索方法，才能获得自己想要的信息。
在AnnotationHub中进行高效的查询，需要注意根据以下规则分类检索：

Data provider
Data class
Species
Data source

names(mcols(hub))
#>  [1] "title"              "dataprovider"       "species"            "taxonomyid"         "genome"            
#>  [6] "description"        "coordinate_1_based" "maintainer"         "rdatadateadded"     "preparerclass"     
#> [11] "tags"               "rdataclass"         "rdatapath"          "sourceurl"          "sourcetype"

2.2 AnnotationHub Data providers, Data classes, Species， Data sources

## AnnotationHub Data providers
unique(hub$dataprovider)

## AnnotationHub Data classes
unique(hub$rdataclass)

## AnnotationHub Species
head(unique(hub$species))
length(unique(hub$species))

## Data sources
unique(hub$sourcetype)

2.3 AnnotationHub query

## AnnotationHub query
qry <- query(hub, c("granges","homo sapiens","ensembl"))
qry
qry$sourceurl

## Selecting AnnotationHub resource
whatIwant <- qry[["AH50377"]]
GRCh38TxDb <- makeTxDbFromGRanges(whatIwant)
GRCh38TxDb

2.4 AnnotationHub exercises

How many resources are on AnnotationHub for Atlantic salmon (Salmo salar)?
Get the most recent Ensembl build for domesticated dog (Canis familiaris) and make a TxDb