在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)是不是很简单呀！
本文参与 腾讯云自媒体分享计划，分享自微信公众号。
原始发表：2018-09-10，如有侵权请联系 cloudcommunity@tencent.com 删除
评论
登录后参与评论
推荐阅读