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

01 安装pysam

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

pip install pysam

02 读入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文件的读入

03

检查bam文件是否存在index：

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

04

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

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

05

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

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

06

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

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

07

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

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

08

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

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

09

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

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

10

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

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

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

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

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

11

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

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

12

对于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。