R语言实现基因组的注释

BiocManager::install("cellbaseR")
BiocManager::install("Gviz")
BiocManager::install("VariantAnnotation")

###创建基础库
library(cellbaseR)
cb <-CellBaseR()

##获取基础数据
res <-getMeta(object=cb, resource="species")

##获取基因信息
genes <-c("TP73","TET1")
res <-getGene(object = cb, ids = genes, resource = "info")

>colnames(res)
 [1] "id"          "name"        "biotype"     "status"    
 [5] "chromosome"  "start"       "end"         "strand"    
 [9] "source"      "description""transcripts" "annotation"

###抽取详细信息
transcripts<- res$transcripts[[1]]
str(transcripts,1)

##获取区域的信息。Resource需要自己进行根据需要进行设置包括genes, snps,clincally relevant variants, proteins等。
res <-getRegion(object=cb,ids="17:1000000-1005000",resource="clinical")
res <-getRegion(object=cb,ids="17:1000000-1005000",resource="conservation")
res <-getRegion(object=cb,ids="17:1000000-1005000",resource="regulatory")
str(res,1)

##获取突变数据
res2 <-getVariant(object=cb, ids="1:169549811:A:G",resource="annotation")

##获取公共突变信息
cbParam <-CellBaseParam(feature=c("TP73","TET1"), limit=100)
res <-getClinical(object=cb,param=cbParam)

##获取SNP的信息
res <-getSnp(object=cb, ids="rs6025", resource="info")

##获取蛋白质的信息
res <-getProtein(object=cb, ids="O15350", resource="info")

##获取转录因子信息
param <- CellBaseParam(limit= 12)
res <-getTf(object=cb, ids="CTCF", resource="tfbs", param=param)

##获取基因在不同数据库中的信息
res <-getXref(object=cb, ids="ENST00000373644", resource="xref")

##基于Gviz实现基因组的可视化
test <-createGeneModel(object = cb, region = "17:1500000-1550000")
if(require("Gviz")){
  testTrack <- Gviz::GeneRegionTrack(test)
  Gviz::plotTracks(testTrack,transcriptAnnotation='symbol')
}

##注释VCF文件
library(VariantAnnotation)
fl <-system.file("extdata","hapmap_exome_chr22_200.vcf.gz",package = "cellbaseR" )
res <-AnnotateVcf(object=cb, file=fl, BPPARAM = bpparam(workers=2),batch_size = 100)