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