寫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 搜索了一下，發(fā)現(xiàn)，似乎pysam不能直接安裝，同時似乎有一個解法
https://pypi.org/project/pysam-win/

pip install pysam-win

OK。似乎這樣就安裝好了?？梢栽趙indows下愉快地使用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不需要排序的嗎？本來想試驗一下，不過發(fā)現(xiàn)還是比較麻煩。繼續(xù)看文檔就知道，必須是先排序，因為

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有一些好處，即幾乎所有變量是全局的（除非在函數(shù)中）。所以對于上述代碼的read，我之前大概翻過python的書，知道完全可以在python交互式操作中，使用

dir(read)
# 或者
help(read)

來查看read對象的文檔。