在上一篇的文章里我详细介绍了 BAM(SAM/CRAM)的格式和一些需要注意的细节, 还说了该如何使用 samtools 在命令行中对其进行操作. 但是很多时候这些操作是不能满足我们的实际需要的, 比如统计比对率, 计算在某个比对质量值之上的 read 有多少, 或者计算 PE 比对的插入片段长度分布, 甚至需要你根据实际情况编写一个新的变异检测算法等. 这个时候往往难以直接通过 samtools 来实现[注] , 而是需要编写专门的程序进行计算. 因此, 在这一篇文章里我们就一起来学习应该如何在程序中借助 Pysam 来处理 BAM 文件.
[注] 关于统计比对率其实是可以通过 samtools stats 计算获得的. 不过我们这篇文章不是为了争辩 samtools 能做什么, 不能做什么, 而是要跟大家讨论该如何编写程序处理 BAM.
不过, 在开始之前我想稍微再补充一下上一节中提到的 CRAM-- 我习惯将其称为 BAM 的高压缩格式, 因为它和 BAM/SAM 的格式基本相同, 但有四点我们需要注意一下:
CRAM 的高压缩是通过借助参考序列和对其他信息的进一步编码来实现的, 它相比于 BAM 有着更高的压缩率, 能够节省 30%-50% 的空间;
CRAM 目前的 IO 效率没有 BAM 高(压得密嘛), 约慢 30%, 但在不断进步, 现在已经更新到了 3.x 版本了;
CRAM 和 BAM 可以通过 samtools 或者 picard 方便地实现互转;
CRAM 一定会取代 BAM, 这话并不是我说的, 而是 bwa/samtools 的作者 lh3 说的.
什么是 Pysam
Pysam 是一个专门用来处理 (BAM/CRAM/SAM) 比对数据和变异数据 (VCF 和 BCF) 的 Python 包. 它的核心是 htslib-- 一个高通量数据处理 API(来自 samtools 和 bwa 的核心, 基于 C 语言), 开发者们用 Python 对它直接进行轻量级包装, 因此能够在 Python 中方便地进行调用, 并且保证了它与原生 C-API 功能上的高度一致.
为什么是 Pysam
因为 Pysam 可以说是最为官方的版本, 有比较固定的开发者在维护, 它的稳定性和可靠性都很高. 虽然还有一些其它的包同样能够处理 BAM 但其实它们大多绕不开对 htslib 的使用, 但却没有 pysam 周全. 而且 Pysam 还集成了 tabix 的接口, 所以除了比对数据之外, 还能够用于处理所有用 tabix 构建过索引的文件, 总之就是全且可靠.
如果是文本格式的 sam 的话, 其实也可以直接将其当作普通文本文件来处理, 不需借助任何程序包(这在早期的数据分析中经常看到这种操作), 只是要麻烦很多(必须自己在程序中处理所有细节, 包括解析 FLAG 和 CIGAR 信息, 以前我也干过不少类似的事情), 甚至我还看到有人直接在程序中调用 samtools view 把 BAM 转换成 SAM 之后再处理的... 这样的做法实在不推荐.
所以, 只要你用的是 Python, 那么 Pysam 真的是目前看来比较好的选择. 当然如果你用 C/C++ 那么直接用 htslib 或者 bamtools, 如果是 Java, 那么直接使用 htsjdk--htslib 的 java 版本.
如何使用 Pysam
首先, 要为我们的 Python 环境安装这个包, 如果已安装过的话可以忽略这一步. 有两个方法, pip 和 bioconda, 都比较简单, 我们这里以 pip--Python 的包管理工具来进行:
$ pip install pysam
安装完成之后我们就可以在 Python 程序中调用 pysam 了.
读取 BAM/CRAM/SAM 文件
Pysam 中的函数有很多, 但是重要的读取函数主要有:
AlignmentFile: 读取 BAM/CRAM/SAM 文件
VariantFile: 读取变异数据(VCF 或者 BCF)
TabixFile: 读取由 tabix 索引的文件;
FastaFile: 读取 fasta 序列文件;
FastqFile: 读取 fastq 测序序列文件;
等以上几个, 其中尤以 AlignmentFile 和 VariantFile 为核心. 需要我们注意到的地方是, Pysam 中的有些函数由于历史原因存在重复, 比如名字上只有大小写的差异, 但功能却是一样的(比如下图的 TabixFile), 有些则只是简化了函数名, 这些情况用的时候留个心眼就行了.
另外, 这篇文章的目的是介绍如何处理比对文件, 所以我打算只介绍 AlignmentFile.
读取比对文件前, 我建议先使用 samtools index 为比对文件构建好索引. 当然如果是 SAM 文件就不必了 -- 它是文本文件, 索引的作用是让我们可以对文件进行随机读取, 而不必总是从头开始.
下面我先用一个例子作为引子:
- import pysam
- bf = pysam.AlignmentFile('in.bam', 'rb')
我在这个例子里面, 先在程序中导入 pysam 包, 然后调用 AlignmentFile 函数读取'in.bam'文件, 并把句柄赋值给了 bf,bf 其实是一个迭代器 --Python 中的术语, 意思就是适合在 for 循环中进行遍历的对象.
这样我们就是可以通过 bf 获取这份比对文件中的内容了. 比如我们想把 in.bam 中每一条 read 的比对位置(包含染色体编号和位置信息), 比对质量值和插入片段长度输出(我们的 in.bam 来自 PE 测序数据的结果), 那么可以这样做:
- import pysam
- bf = pysam.AlignmentFile('in.bam', 'rb')
- for r in bf:
- print r.reference_name, r.pos, r.mapq, r.isize
是不是很简单! 打开 in.bam 文件之后, 用 for 循环对其从头到尾地遍历, 并把每个值都赋给 r,r 在这里代表的就是比对的 read 信息, 它是一个对象(在 Pysam 由 AlignedSegment 定义), 通过它就可以获取所有的比对信息, 比如上面例子中:
r.reference_name 代表 read 比对到的参考序列染色体 id;
r.pos 代表 read 比对的位置;
r.mapq 代表 read 的比对质量值;
r.isize 代表 PE read 直接的插入片段长度, 有时也称 Fragment 长度;
这里例子的结果如下:
- chrM 160 50 235
- chrM 161 30 -283
- chrM 314 60 -207
- ...
另外, 由于 bf 是一个迭代器, 我们其实还可以用 bf.next()一个一个地对其进行访问, 而不必在 for 循环中遍历, 这在一些特殊的情况下, 这个做法是非常有用且方便的.
当然, 上面这个例子其实非常简单, 实际上 r 变量中还有很多其它关于比对的信息, 下面这个截图, 就是变量中能够获取到的所有比对相关的信息, 有好几十个.
眼尖的同学可能也发现了, 这里面存在一些名字类似的变量, 如: r.mapping_quality 和 r.mapq, 它们其实都是比对质量值. 类似的也还有几个, 这都是上面我提到的历史原因所致, 不过这种多余变量随着 Pysam 的维护也正在逐步变少.
此外, Pysam 中的位点坐标体系是 0-base(意思是染色体的起始位置是从 0 而不是 1 开始算的)而不是 1-base, 所以上面的输出的 160, 其实真实位置应该要 + 1, 也就是 161.
还有, 上文我也说过, AlignmentFile 除了能够读 / 写 BAM 之外, 还同样能够读 / 写 CRAM 和 SAM. 区别就在于函数中的第二个参数, 比如上面例子中的字符'b'就是用于明确指定 BAM 文件,'r'字符代表 "只读" 模式 (read 首字母). 如果要打开 CRAM 文件, 只需要把 b 换成 c(代表 CRAM) 就行了, 如下:
- import pysam
- cf = pysam.AlignmentFile('in.cram', 'rc')
那么, 如果是 SAM 文件呢? 去掉 b 或 c 即可:
- import pysam
- sf = pysam.AlignmentFile('in.sam', 'r')
读取特定比对区域内的数据
有时候我们并不需要遍历整一份 BAM 文件, 我们可能只想获得区中的某一个区域(比如 chrM 中 301-310 中的信息), 那么这个时候可以用 Alignmen 模块中的 fetch 函数:
- import pysam
- bf = AlignmentFile('in.bam', 'rb')
- for r in bf.fetch('chrM', 300, 310):
- print r
- bf.close()
通过 fetch 函数就可以定位特定区域了, 非常方便. 不过, 这个时候输入文件 in.bam 就必须要有索引, 不然无法实现这种读取操作. 最后用完了, 要记得关闭文件(bf.close()).
来个稍微难一点的例子
问题: 如何输出覆盖在某个位置上, 比对质量值大于 30 的所有碱基?
这个问题包含两个部分:
固定的某个位置(我们这里还是用 chrM 301 这个位置)
read 比对质量值必须是大于 30.
如何做呢? 这个时候我们要用 AlignmentFile 模块的另一个函数 --pileups 来协助解决, 代码如下:
- import pysam
- bf = pysam.AlignmentFile("in.bam", "rb" )
- for pileupcolumn in bf.pileup("chrM", 300, 301):
- for read in [al for al in pileupcolumn.pileups if al.alignment.mapq>30]:
- if not read.is_del and not read.is_refskip:
- if read.alignment.pos + 1 == 301:
- print read.alignment.reference_name,\
- read.alignment.pos + 1,\
- read.alignment.query_sequence[read.query_position]
- bf.close()
这段代码看起来虽然简单, 但其实包含了很多信息. 总的来说, 就是通过 pileup 获取了所有覆盖到该位置的 read, 并将其存到 pileupcolumn 中. 然后, 对 pileupcolumn 调用 pileups(注意多了一个 s)获得一条 read 中每个比对位置的信息 (一条 read 那么长, 并非只覆盖了一个位置), 然后通过判断语句留下覆盖到目标位点(301) 的碱基. 代码中的 read.alignment 是 Pysam 中 AlignedSegment 对象, 它包含的内容和上述其它例子中的 r 是一样的. read.alignment.pos + 1 还是 0-base 的原因. 最后结果如下:
- chrM 301 A
- chrM 301 A
- chrM 301 A
- chrM 301 C
- chrM 301 C
- chrM 301 C
- chrM 301 C
- chrM 301 C
- chrM 301 C
- chrM 301 C
- ...
创建 BAM/CRAM/SAM 文件
最后这个例子, 我想告诉大家该如何用 Pysam 输出 BAM/CRAM/SAM 格式, 具体还是看代码吧, 这里想输出结果是 BAM 文件, 所以输出模式是 "wb", 例子中我们只输出一条比对结果作为说明.
- import pysam
- header = {'HD': {'VN': '1.0'},
- 'SQ': [{'LN': 1575, 'SN': 'chr1'},
- {'LN': 1584, 'SN': 'chr2'}]
- }
- tmpfilename = "out.bam"
- with pysam.AlignmentFile(tmpfilename, "wb", header=header) as outf:
- a = pysam.AlignedSegment() # 定义一个 AlignedSegment 对象用于存储比对信息
- a.query_name = "read_28833_29006_6945"
- a.query_sequence="AGCTTAGCTAGCTACCTATATCTTGGTCTTGGCCG"
- a.flag = 99
- a.reference_id = 0
- a.reference_start = 32
- a.mapping_quality = 20
- a.cigar = ((0,10), (2,1), (0,25))
- a.next_reference_id = 0
- a.next_reference_start=199
- a.template_length=167
- a.query_qualities = pysam.qualitystring_to_array("<<<<<<<<<<<<<<<<<<<<<:<9/,&,22;;<<<")
- a.tags = (("NM", 1),
- ("RG", "L1"))
- outf.write(a)
小结
我写这篇文章的目的主要有两个: 第一, 充实上一篇文章中关于如何操作 BAM 的内容; 第二, 介绍 Pysam 这一个值得使用的包给大家. 另外, 我上面列举的例子其实都比较偏于基础操作, 这可能和我自身对认知的看法有关. 我一直认为, 只有真正理解并灵活地应用基础操作, 才可以灵活地解决一切复杂的问题.
而且, 上面几个例子中用到的模块和函数其实都是比较常用的, 所以我比较推荐优先掌握它们. 这些例子里面用到的数据我已放到 GitHub 了, 感兴趣的同学可以在公众号后台回复 "WGS" 即可获得, 后续也会陆续有其它的代码和数据可供参考.
最后, Pysam 的内容其实还有很多, 我所介绍的也仅在对于比对数据的处理, 其它很多的模块和函数, 包括对 Fasta,Fastq,VCF,BCF 和 Tabix 文件的处理, 我就不进行一一介绍了, 建议大家在使用的时候多看看它的完整文档 http://pysam.readthedocs.io/en/latest/API.html .
来源: http://www.tuicool.com/articles/go/VFnuuub