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

健明 2021-07-14

展开全文

坐标注释最简单的生物学应用就是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)