寫perl的思維，可能確實(shí)不能拿來學(xué)python。畢竟，python的褲子有很多。面向?qū)ο蟮恼Z言，如果不好好穿褲子，怎么找對(duì)象？。手上要做的事情，需要解析sam，更或者bam文件。當(dāng)然，如果有可能的話，還需要對(duì)SAM或者BAM進(jìn)行排序！
這個(gè)事情我在java寫過，but，最后效率不如htsjdk，故最后還是打包了htsjdk。吸取這個(gè)教訓(xùn)，使用python的時(shí)候。第一步，先用pysam。

pysam的下載與安裝

此處，直接接上一個(gè)推文，從pycharm中，右擊某個(gè)項(xiàng)目，就可以直接打開terminal

image.png

在網(wǎng)絡(luò)暢通的情況下，使用pip安裝pysam（其實(shí)我也不知道，windwos下是否可以）

pip install pysam

很遺憾，安裝失敗了。各種報(bào)錯(cuò)，不忍直視。
百度 google 搜索了一下，發(fā)現(xiàn)，似乎pysam不能直接安裝，同時(shí)似乎有一個(gè)解法
https://pypi.org/project/pysam-win/

pip install pysam-win

OK。似乎這樣就安裝好了?？梢栽趙indows下愉快地使用python處理SAM/BAM文件了。

pysam的使用與目標(biāo)

首先，當(dāng)然是看官方文檔，看看都怎么用的，https://pysam.readthedocs.io/en/latest/

從文檔來看，pysam的速度應(yīng)該不會(huì)慢，畢竟是**a lightweight wrapper of the htslib C-API.

**

隨后，目標(biāo)如下：

1. 讀取SAM，BAM文件
2. BAM文件排序
3. ....沒有了

讀取

先打開一個(gè)文件對(duì)象 AlignmentFile

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

然后就可以自由自在地獲取任意region的reads...我總是覺得我似乎在什么時(shí)候用過pysam？還是....哦，感覺很像我寫的Jsam...

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

只是，此處的BAM不需要排序的嗎？本來想試驗(yàn)一下，不過發(fā)現(xiàn)還是比較麻煩。繼續(xù)看文檔就知道，必須是先排序，因?yàn)?/p>

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

如果需要遍歷一個(gè)文件的話，那么直接

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，其他的似乎暫時(shí)對(duì)我來說并不重要。

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")

寫一個(gè)用得到的小腳本

課題需要，所以要寫一個(gè)python小腳本，需要完成以下內(nèi)容

1. 讀取SAM或者BAM文件（默認(rèn)要求name-sorted）
2. 過濾去除mapped pos大于20個(gè)位置的，因?yàn)榛究梢哉J(rèn)為是repeat
3. 對(duì)文件進(jìn)行pos-sorted
4. 課題內(nèi)容，就不寫出來了

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")

補(bǔ)充

我是在windows下面寫代碼的，但是經(jīng)過嘗試，pysam確實(shí)幾乎無法在windows下面正常配置，所以寫代碼的過程是痛苦的。
無法進(jìn)行代碼補(bǔ)全的面向?qū)ο蟠a碼，簡直要命！
最后的解法是，只能在linux下，查看當(dāng)前對(duì)象的屬性...

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ù)中）。所以對(duì)于上述代碼的read，我之前大概翻過python的書，知道完全可以在python交互式操作中，使用

dir(read)
# 或者
help(read)

來查看read對(duì)象的文檔。