在R里面对坐标进行映射
比如把自己制作好的bam文件的坐标，跟提取自gtf文件的坐标信息对应起来，使用GenomicRanges包自带的函数即可。
ann1 <- data.frame(
GeneID=c("gene1","gene1","gene2","gene2"),
Chr="chr_dummy",
Start=c(100,1000,3000,5000),
End=c(500,1800,4000,5500),
Strand=c("+","+","-","-"),
stringsAsFactors=FALSE)
ann1
ann1 <- with(ann1, GRanges(as.character(Chr),
IRanges(as.numeric(Start), as.numeric(End)),
as.character(Strand),
id = as.character(GeneID)))
ann1
ann2 <- data.frame(
peakID=c("peak1","peak1","peak2","peak2"),
Chr="chr_dummy",
Start=c(50,1500,3100,5400),
End=c(400,1900,4700,5500),
Strand=c("+","+","-","-"),
stringsAsFactors=FALSE)
ann2
ann2 <- with(ann2, GRanges(as.character(Chr),
IRanges(as.numeric(Start) , as.numeric(End)),
as.character(Strand),
id = as.character(peakID)))
ann2
gr3 = intersect(ann1,ann2)
gr3library(ChIPpeakAnno)
o = findOverlaps(ann1,ann2)
str(o)
lo=cbind(as.data.frame(ann1[queryHits(o)]),
as.data.frame(ann2[subjectHits(o)]))
head(lo)
gr3这里的重点其实是grange对象和intersect及findOverlaps函数的使用。
关于 GRANGE对象
三年前我在生信菜鸟团博客就多次强调过这个重点了，在R里面处理生物信息学数据是躲不过这个定义的，有点类似于各式各样的生物信息学文件格式，是一个标准。
对这个grange对象也会有很多很多的方法，假设有一个grange对象命名为exon_txdb，来自于代码
library("TxDb.Hsapiens.UCSC.hg19.knownGene")
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
exon_txdb=exons(txdb)
genes_txdb=genes(txdb)那么操作它的函数有：
seqnames(exon_txdb)返回一个class ‘Rle’ [package “S4Vectors”] with 4 slots，有93个染色体信息，以及每条染色体上面有多少个外显子信息
ranges(exon_txdb)返回外显子的起始终止位点，长度，以及其它信息，也是一个对象class ‘IRanges’ [package “IRanges”] with 6 slots
strand(exon_txdb)返回外显子的正负链信息，要么在正链要么在负链
mcols(exon_txdb)返回exon的id编号，1到27750个
seqlengths(exon_txdb)返回每条染色体的长度信息
names
length
GRanges对象还有很多其它类型的操作，非常好玩的，split,shift,resize,flank,reduce,gaps,disjoin,coverage
其它求交集并集和都可以用,union,intersect,setdiff,pintersect,psetdiff
关于 FINDOVERLAPS函数
本来应该是ChIPpeakAnno包带有的一个非常实用的peaks分析小工具，在我的GitHub很早以前关于ChIP-seq分析流程代码分析里面有提到。
https://github.com/jmzeng1314/NGS-pipeline
评论专区