威尼斯5848vip电子游戏(wns认证网站)-BinG百科

市场动态MARKET DYNAMICS

实用生信小技巧：利用pysam高效处理bam文件

发布时间：2021-4-19 11:18:04阅读次数：次分享到：

★

★ ★ ★

★

做生信的小伙伴一定对bam文件很熟悉，这种文件是由测序数据和参考基因组比对得到的。文件主要分为两个部分：一部分为头文件，以@符号开头；一部分为主文件，记录了详细的比对信息，每行共11列。详细的bam文件说明可以参考：

https://www.cnblogs.com/raisok/p/10917769.html。这里就不再赘述啦。

由于bam文件是二进制的压缩文件，正常情况下我们无法用普通编辑器打开查看，一般情况下我们借助samtools软件来对文件进行处理。但很多小伙伴一定遇到过这样的问题：如果我想对bam文件进行个性化的处理，需要怎么做呢？常规的做法是利用samtools将bam文件先转化成可读写的sam文件，再进行后续的处理。这样做一方面牺牲了大量的时间和存储，另一方面在编程时需要调用shell命令十分繁琐，今天介绍的pysam可以说是一条龙服务完美解决问题。

安装pysam

pysam是python的一个库，安装好python以后可以利用python自带的pip安装pysam。在后续使用中请注意，samtools软件的坐标是1-base，而pysam对于基因组坐标的处理是0-base：

pip install pysam

读入bam文件

pysam可读入的文件格式有bam、sam、cram，以test前缀的文件作为示例，读入方法如下，读入后生成对象bf，后续对文件的操作都基于对bf的访问：

import pysam

bf = pysam.AlignmentFile("test.bam", "rb") #对bam文件的读入

bf = pysam.AlignmentFile("test.sam", "r") #对sam文件的读入

bf = pysam.AlignmentFile("test.cram", "rc") #对cram文件的读入

检查bam文件是否存在index：

bf.check_index() #返回True则存在索引，返回False则索引不存在

pysam封装了samtools，因此可以避开命令行调用直接对bam文件进行sort：

pysam.sort("-o", "output.bam", "test.bam")

想要提取bam文件的头文件，可以这样操作：

head=bf.header #返回一个迭代器，可以利用循环来访问

想要提取bam文件的前10行比对结果，可以这样操作：

out=bf.head(10) #该命令会以迭代器的形式返回bam文件的前10行

如果想提取基因组某个区域的比对结果，比如现在我想提取bam文件中比对到血红蛋白基因上的结果，可以这样操作：

hemoglobin_region=bf.fetch(contig="chr11",start=5246694,stop=5250625) #结果同样是以迭代器的形式返回

想要计算落在某个区间的reads数目，可以这样操作：

readsNumber=bf.count(contig="chr1",start=1,stop=1000000) #返回reads数目值

想要计算某个区间的每个碱基上的覆盖深度，可以这样操作：

baseCov=bf.count_coverage(contig="chr1",start=1,stop=1000000) #返回一个数组，数组的长度和区间长度一致，表示每个碱基的覆盖深度

对参考基因组的信息进行统计，可以这样操作：

bf.get_index_statistics() #返回一个元组，每个元素记录一条染色体比对上的reads数目和未比对上的reads数目

bf.lengths #返回一个元组，每个元素记录一条染色体长度

bf.references #以元组的形式返回参考基因组每条染色体的名称

bf.nreferences #返回参考基因组的染色体条数

对比对信息进行统计，可以这样操作：

bf.mapped #返回比对上的reads数目

bf.unmapped #返回未比对上的reads数目

对于bf对象可以利用for循环进行访问，每个循环是一个read，对于read的各种处理操作汇总如下：

for read1 in bf:

read2=bf.mate(read1) #获取与read1配对的read2的比对信息

read1.cigarstring #字符串形式返回read1的cigar值

read1.is_duplicate #返回布尔型变量，判断read1是否为PCR重复

read1.is_paired #返回布尔型变量，判断read1是否为成对reads

read1.is_read1 #返回布尔型变量，判断read1是否为左端reads

read1.is_read2 #返回布尔型变量，判断read1是否为右端reads

read1.is_reverse #返回布尔型变量，判断read1是否为反向比对

read1.is_secondary #返回布尔型变量，判断read1是否为次级比对

read1.is_unmapped #返回布尔型变量，判断read1是否比对上参考基因组

read1.is_qcfail #返回布尔型变量，判断read1是否qc质控失败

read1.mapping_quality #返回数值型变量，代表read1的比对质量

read1.mate_is_reverse #返回布尔型变量，判断read1的配对read2是否为反向比对

read1.get_blocks() #返回数值型元组，表示read1在参考基因组比对的起始坐标和终止坐标

read1.get_overlap(5246694,5250625) #返回数值，表示read1在基因组上的坐标和输入参数坐标的重叠区间长度

read1.get_reference_sequence() #字符串形式返回read1所在参考基因组位置的参考序列

小结：

pysam库的使用可以大大减少在脚本中对bash命令行的调用，一方面美化代码便于后期维护更新，另一方面可以个性化调取reads信息。本篇代码涵盖了pysam的常用方法，其他参数详解可见官网https://pysam.readthedocs.io/en/latest/usage.html。

上一条： Promoter capture Hi-C解析染色质三维结构及增强子活性变化在人类细胞多能性状态转变过程中的重要作用

下一条： 5848vip威尼斯电子游戏全转录组线上培训课开始报名啦！！！

威尼斯5848vip电子游戏(wns认证网站)-BinG百科

农学科研

医学临检

市场与支持

关于5848vip威尼斯电子游戏

加入我们

实用生信小技巧：利用pysam高效处理bam文件