把bam文件读入R，并且转为grange对象

假如你的Windows电脑有个bam文件，不想传输到linux服务器去使用samtools等命令行工具来探索它，就可以使用R语言！

有成熟的R包可以把bam文件读入R，比如Rsamtools，很简单的代码：

library(Rsamtools)
bamFile="alignResults.BAM"
quickBamFlagSummary(bamFile)
# https://kasperdanielhansen.github.io/genbioconductor/html/Rsamtools.html
bam <- scanBam(bamFile)
bam

值得注意的是，这里我虽然不再演示了，但是作为初学者的你，应该是知道

但是把读入的数据变成grange对象就需要一点点技巧，下面演示如何创建grange对象samtools等命令行工具有多复杂的功能和技巧， 那么这个R包就可以多复杂，如果你学习足够努力，那就发一个你比较Rsamtools和samtools命令行工具的心得笔记给我吧，我会给你惊喜的，我的邮箱是 jmzeng1314@163.com

names(bam[[1]])
tmp=as.data.frame(do.call(cbind,lapply(bam[[1]], as.character)))
tmp=tmp[tmp$flag!=4,] # 60885 probes
#  intersect() on two GRanges objects.
library(GenomicRanges)
my_seq <- with(tmp, GRanges(as.character(rname), 
                                 IRanges(as.numeric(pos)-60, as.numeric(pos)+60), 
                                 as.character(strand), 
                                 id = as.character(qname)))

得到对象如下：

关于 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

■   ■   ■

赞 (0)

相关推荐