在R里面对坐标进行基因组区域注释

坐标注释最简单的生物学应用就是peaks区域的注释，通常我们可以使用linux的各种软件加上gtf等格式的基因组注释信息来完成，在R里面当然也是可以轻松完成的啦！

假设有如下格式的坐标：

> head(pos)
    chr     start       end
1 chr10 100505299 100505300
2 chr10 100505299 100505300
3 chr10 104125494 104125495
4 chr10  11320827  11320828
5 chr10 118691247 118691248
6 chr10 119123605 119123606

这里可以使用大名鼎鼎的Y书开发的ChIPseeker包，加上人类的注释信息包TxDb.Hsapiens.UCSC.hg38.knownGene来进行注释，示例代码如下：

pos=data.frame(chr=str_split(dat$id,':',simplify = T)[,1],
                  start=as.numeric(str_split(dat$id,':',simplify = T)[,2]) )
pos$end=pos$start+1 
pos_anno=as.data.frame(peakAnno)
require(ChIPseeker)
library(org.Hs.eg.db)
library(org.Mm.eg.db)
library(GenomicRanges)
peak <- GRanges(seqnames=Rle(pos[,1]),
                ranges=IRanges(pos[,2], pos[,3]), strand=rep(c("*"), nrow(pos)))
peak
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
txdb=TxDb.Hsapiens.UCSC.hg38.knownGene
peakAnno <- annotatePeak(peak, tssRegion=c(-3000, 3000),
                         TxDb=txdb, annoDb="org.Hs.eg.db")
pos_anno=as.data.frame(peakAnno)

是不是很简单呀！

猜你喜欢： 生信基础知识100讲

独家福利

如果需要组装自己的服务器；代办生物信息学服务器

如果需要帮忙下载海外数据(GEO/TCGA/GTEx等等)，点我？

如果需要线下辅导及培训，看招学徒

如果需要个人电脑：个人计算机推荐

如果需要置办生物信息学书籍，看：生信人必备书单

如果需要实习岗位：实习职位发布

如果需要售后：点我