在上一篇的文章里我详细介绍了 BAM(SAM/CRAM)的格式和一些需要注意的细节,还说了该如何使用 samtools 在命令行中对其进行操作。但是很多时候这些操作是不能满足我们的实际需要的,比如统计比对率、计算在某个比对质量值之上的 read 有多少,或者计算 PE 比对的插入片段长度分布,甚至需要你根据实际情况编写一个新的变异检测算法等。这个时候往往难以直接通过 samtools 来实现【注】,而是需要编写专门的程序进行计算。因此,在这一篇文章里我们就一起来学习应该如何在程序中借助 Pysam 来处理 BAM 文件。
【注】关于统计比对率其实是可以通过 samtools stats 计算获得的。不过我们这篇文章不是为了争辩 samtools 能做什么,不能做什么,而是要跟大家讨论该如何编写程序处理 BAM。
不过,在开始之前我想稍微再补充一下上一节中提到的 CRAM——我习惯将其称为 BAM 的高压缩格式,因为它和 BAM/SAM 的格式基本相同,但有四点我们需要注意一下:
Pysam 是一个专门用来处理(BAM/CRAM/SAM)比对数据和变异数据(VCF 和 BCF)的 Python 包。它的核心是 htslib——一个高通量数据处理 API(来自 samtools 和 bwa 的核心,基于 C 语言),开发者们用 Python 对它直接进行轻量级包装,因此能够在 Python 中方便地进行调用,并且保证了它与原生 C-API 功能上的高度一致。
因为 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 了。
Pysam 中的函数有很多,但是重要的读取函数主要有:
等以上几个,其中尤以 AlignmentFile 和 VariantFile 为核心。需要我们注意到的地方是,Pysam 中的有些函数由于历史原因存在重复,比如名字上只有大小写的差异,但功能却是一样的(比如下图的 TabixFile),有些则只是简化了函数名,这些情况用的时候留个心眼就行了。
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 定义),通过它就可以获取所有的比对信息,比如上面例子中:
这里例子的结果如下:
- 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()。
这个问题包含两个部分:
如何做呢?这个时候我们要用 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
- ...
最后这个例子,我想告诉大家该如何用 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://www.jianshu.com/p/9ac58493dce4