R专题：Rsamtools的简单探索

Rsamtools安装bam文件导入简单的操作

Rsamtools是一个用来读取bam文件并进行相应操作的R包，它是很多R生信包的基础，如ShortRead（辅助创建AlignedRead 对象，用于质量评价和read操作）和GenomicAlignments （辅助创建GAlignments 对象，常用于RNA-seq）。

Rsamtools安装

Rsamtools是一个bioconductor工具，安装方法如下：

## try http:// if https:// URLs are not supported
source("https://bioconductor.org/biocLite.R")
biocLite("Rsamtools")

安装Rsamtools时会有很多依赖包会被安装，虽然这个过程是自动完成的，但是由于依赖包过多，有时候还是容易有一些问题（其他包同样道理），所以在安装后导入Rsamtools时，一定要看一遍警告和错误信息，如果有错误信息，要逐一排查。我当时安装完Rsamtools后，导入它并没有工作，看提示是需要需要RCurl，安装后RCurl才可以使用（install.packages("RCurl"))。

bam文件导入

简单来说，bam文件的导入支持以下3种情况：

1. 可以直接导入bam全文件；

2. 逐次从bam文件中读取特定数目的read；

3. 读入特定基因组范围的特定内容；

直接导入bam文件

## 以Rsamtools包安装目录下自带的ex1.bam为示例文件

#1. BamFile创建文件对象
library(Rsamtools)
# 找到ex1.bam的目录全名
bamPath <- system.file("extdata", "ex1.bam", package="Rsamtools")
# BamFile创建文件对象，此时还不能直接查看bam文件的文本内容
bamFile <- BamFile(bamPath)
bamFile # 会给出文件对象的简单信息，如图1所示

#2. scanBam读取文本内容，aln就是整个bam文件内容，list格式

aln <- scanBam(bamFile)# scanBam读取BAM文件的文本内容，其结果为一个list文件
length(aln) # 长度为1，因为sanBam可以读取多个bam文件，并合成一个list文件
class(aln) # list
# 取出bam文件的内容，并查看bam文件各个列下的第一个read信息
aln<-aln[[1]] 
names(aln) # bam文件的各列名，"qname"  "flag"   "rname"  "strand" "pos" 等共13个
lapply(aln, function(x) x[1])

图1

bamFile是一个bam文件对象，如上图，它有多个属性，如isopen、yieldSize、obeyQname、asMates等等，他们含义为：

isopen：是否为打开状态，用于逐次读取read（下述）；
yieldSize：逐次读取的read数目，需要文件处于打开状态；
obeyQname:Qname是read名，此参数意味bam文件是否安装read名排序；
asMates：是否为mate pair；
qnamePrefixEnd：配合asMates使用，用于标识和去除read名中的前缀；
qnameSuffixStart：配合asMates使用，用于标识和去除read名中的后缀；

aln是一个列表，其列表元素的名称基本上是bam文件的列名，重要的名称如下：

qname：读段的名称；
rname：染色体/序列/ contig的名称；
strand：比对的链；
pos：比对的最左边部分的坐标；
qwidth：读段的长度；
mapq：对齐的映射质量；
seq：对齐的实际顺序；
qual：对齐的质量字符串；
cigar：CIGAR字符串；
flag：flag；

逐次导入read

##每次读取bam文件中一定数量的内容，本例中依次读取一个read

open(bamFile) # 文件对象需处于打开状态
yieldSize(bamFile) <- 1 #每次读入一个read
scanBam(bamFile)[[1]]$seq
scanBam(bamFile)[[1]]$seq # 重复使用本条命令，可以依次读取每一条read的序列
close(bamFile) # 关闭文件
yieldSize(bamFile) <- NA # 取消逐次读取数量

如果bam文件对象没有处于打开状态，则逐次读取失效。如果已经读取到最后一个read，则需要先关闭文件对象，再重新打开。

读取特定基因组范围的特定内容

## ScanBamParam有两个参数which和what，分别用于指定读取的read范围以及相应的read内容（即序列名、参考序列、flag等等）

# 先指定which内容，本例的which代表取出seq1中的两个区域：read1-3、read2-4（不包括3与4）
whichs <- GRanges(seqnames = "seq1",
              ranges = IRanges(start = c(1, 2), end = c(3,4)))
# what指定需要显示的bam内容
whats <- c("qname","rname", "strand", "pos", "qwidth", "seq")
params <- ScanBamParam(which = whichs, what = whats)
aln2 <- scanBam(bamFile, param = params)
names(aln2) # "seq1:1-3" "seq1:2-4"
aln2[[1]]$qname # 显示seq1:1-3的两个read名

读取bam的特定区域及特定内容只需要使用ScanBamParam构造一个参数对象，然后使用scanBam读取文件时，将参数传给scanBam即可。

假如查看了部分内容后，想要查看bam的全部列名，可以使用what = scanBamWhat()，则会重置上一次选取的字段范围，如params <- ScanBamParam(which = whichs, what = scanBamWhat()) 。

简单的操作

# bamFile是使用BamFile创建的bam文件对象

seqinfo(bamFile) # 参考序列信息
countBam(bamFile) # bam文件的read数及核酸数等信息
idxstatsBam(bamFile) # 序列比对情况
scanBamHeader(bamFile) # head信息

sortBam #排序，具体可以在R中查看帮助（?sortBam）
indexBam #索引，具体可以在R中查看帮助（?indexBam）
mergeBam #合并，具体可以在R中查看帮助（?mergeBam）

参考资料

1. Rsamtools的bioconductor主页http://www.bioconductor.org/packages/release/bioc/html/Rsamtools.html

2. Rsamtools参考手册http://www.bioconductor.org/packages/release/bioc/manuals/Rsamtools/man/Rsamtools.pdf

3. Rsamtools的生信技能树论坛介绍http://www.biotrainee.com/forum.php?mod=viewthread&tid=1820&extra=page%3D1%26filter%3Dtypeid%26typeid%3D52

---

猜你喜欢

还有更多文章，请移步公众号阅读

▼ 如果你生信基本技能已经入门，需要提高自己，请关注下面的生信技能树，看我们是如何完善生信技能，成为一个生信全栈工程师。

▼ 如果你是初学者，请关注下面的生信菜鸟团，了解生信基础名词，概念，扎实的打好基础，争取早日入门。