在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)
gr3
library(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

■   ■   ■

赞 (0)

相关推荐