坐标注释最简单的生物学应用就是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)
是不是很简单呀!
联系客服