用python做生物信息数据分析（2-pysam）

写perl的思维，可能确实不能拿来学python。毕竟，python的裤子有很多。面向对象的语言，如果不好好穿裤子，怎么找对象？。手上要做的事情，需要解析sam，更或者bam文件。当然，如果有可能的话，还需要对SAM或者BAM进行排序！
这个事情我在java写过，but，最后效率不如htsjdk，故最后还是打包了htsjdk。吸取这个教训，使用python的时候。第一步，先用pysam。

pysam的下载与安装

此处，直接接上一个推文，从pycharm中，右击某个项目，就可以直接打开terminal

image.png

在网络畅通的情况下，使用pip安装pysam（其实我也不知道，windwos下是否可以）

pip install pysam

很遗憾，安装失败了。各种报错，不忍直视。
百度 google 搜索了一下，发现，似乎pysam不能直接安装，同时似乎有一个解法
https://pypi.org/project/pysam-win/

pip install pysam-win

OK。似乎这样就安装好了。可以在windows下愉快地使用python处理SAM/BAM文件了。

pysam的使用与目标

首先，当然是看官方文档，看看都怎么用的，https://pysam.readthedocs.io/en/latest/

从文档来看，pysam的速度应该不会慢，毕竟是**a lightweight wrapper of the htslib C-API.

**

随后，目标如下：

1. 读取SAM，BAM文件
2. BAM文件排序
3. ....没有了

读取

先打开一个文件对象 AlignmentFile

import pysam
samfile = pysam.AlignmentFile("ex1.bam", "rb")

然后就可以自由自在地获取任意region的reads...我总是觉得我似乎在什么时候用过pysam？还是....哦，感觉很像我写的Jsam...

for read in samfile.fetch('chr1', 100, 120):
     print read
samfile.close()

只是，此处的BAM不需要排序的吗？本来想试验一下，不过发现还是比较麻烦。继续看文档就知道，必须是先排序，因为

Without an index, random access via fetch() and pileup() is disabled.

如果需要遍历一个文件的话，那么直接

import pysam
samfile = pysam.AlignmentFile("Treat.20M.merged.sam")
lineCount = 0
for read in samfile:
    print(read)
    lineCount = lineCount + 1
    if lineCount > 10:
        break

即可。
大体过了文档，整体感觉是，我需要的可能只有pysam，而可能真的不需要biopython，毕竟有了IO，其他的似乎暂时对我来说并不重要。

pysam支持多种生信常见文件格式，包括SAM BAM CRAM fastq fasta vcf gtf gff ....

写出

import pysam
samfile = pysam.AlignmentFile("ex1.bam", "rb")
pairedreads = pysam.AlignmentFile("allpaired.bam", "wb", template=samfile)
for read in samfile.fetch():
     if read.is_paired:
             pairedreads.write(read)

pairedreads.close()
samfile.close()

排序

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

写一个用得到的小脚本

课题需要，所以要写一个python小脚本，需要完成以下内容

1. 读取SAM或者BAM文件（默认要求name-sorted）
2. 过滤去除mapped pos大于20个位置的，因为基本可以认为是repeat
3. 对文件进行pos-sorted
4. 课题内容，就不写出来了

import pysam
# filter sam file to remove reads mapped to repeat regions.\
samfile = pysam.AlignmentFile("Treat.20M.merged.sam")
tmpfile = pysam.AlignmentFile("dedup.sam", "w", template=samfile)
lineCount = 0
max_hit_num = 3
pre_read_id = ""
cur_read_list = list()
for read in samfile:
    cur_read_id = read.qname
    if cur_read_id == pre_read_id:
        cur_read_list.append(read)
    else:
        if len(cur_read_list) < max_hit_num:
            for cur_read in cur_read_list:
                tmpfile.write(cur_read)
        cur_read_list.clear()
        cur_read_list.append(read)
        pre_read_id = cur_read_id
    lineCount = lineCount + 1
    if lineCount > 100:
        break
if len(cur_read_list) != 0 & len(cur_read_list) < max_hit_num:
    for cur_read in cur_read_list:
        tmpfile.write(cur_read)
cur_read_list.clear()
# sort sam file
pysam.sort("-o", "dedup.sorted.sam", "dedup.sam")

补充

我是在windows下面写代码的，但是经过尝试，pysam确实几乎无法在windows下面正常配置，所以写代码的过程是痛苦的。
无法进行代码补全的面向对象码码，简直要命！
最后的解法是，只能在linux下，查看当前对象的属性...

import pysam
samfile = pysam.AlignmentFile("Treat.20M.merged.sam")
lineCount = 0
for read in samfile:
    print(dir(read))
    lineCount = lineCount + 1
    if lineCount > 10:
        break

python有一些好处，即几乎所有变量是全局的（除非在函数中）。所以对于上述代码的read，我之前大概翻过python的书，知道完全可以在python交互式操作中，使用

dir(read)
# 或者
help(read)

来查看read对象的文档。