使用Python怎么对BAM进行处理-创新互联
本篇文章为大家展示了使用Python怎么对BAM进行处理,内容简明扼要并且容易理解,绝对能使你眼前一亮,通过这篇文章的详细介绍希望你能有所收获。
成都创新互联是专业的富拉尔基网站建设公司,富拉尔基接单;提供网站制作、网站设计,网页设计,网站设计,建网站,PHP网站建设等专业做网站服务;采用PHP框架,可快速的进行富拉尔基网站开发网页制作和功能扩展;专业做搜索引擎喜爱的网站,专业的做网站团队,希望更多企业前来合作!什么是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)
上述内容就是使用Python怎么对BAM进行处理,你们学到知识或技能了吗?如果还想学到更多技能或者丰富自己的知识储备,欢迎关注创新互联行业资讯频道。
本文标题:使用Python怎么对BAM进行处理-创新互联
网页URL:http://scyanting.com/article/dooogp.html