BAM 文件：使用 pysam 获取特定位置的所有读取答案

【问题标题】：BAM file: getting all reads on certain position with pysamBAM 文件：使用 pysam 获取特定位置的所有读取
【发布时间】：2018-06-30 17:20:10
【问题描述】：

我有一个 BAM 文件，在某个位置读取 520817（如 IGV 中所示）。但是，当我使用 pysam 获取特定位置上的读取名称和相关核苷酸时，到目前为止我没有得到那个数量（仅获得大约 7000 个读取）。我想只有当那个位置的核苷酸与参考基因组不同时，我才会得到读数。有没有解决方法，所以我得到了所有的读数？我是从生物信息学开始的……所以请告诉我你还需要什么来帮助我！

非常感谢！

这是我的代码：

import pysam
import csv
import sys

#---Get a table with in the first column: read-ID; second column: SNP-location; third column: nucleotide---#
mybam = pysam.AlignmentFile("file.bam", "rb")
w = csv.writer(open("snp.csv", "wb"), delimiter=",")
w.writerow(["Read", "Loc", "Nucl"])
for pileupcolumn in mybam.pileup('chr6', 29911198,29911199):
    print ("\ncoverage at base %s = %s" %
           (pileupcolumn.pos, pileupcolumn.n))
    for pileupread in pileupcolumn.pileups:
        if not pileupread.is_del:
            if pileupcolumn.pos == 29911198:
                w.writerow((pileupread.alignment.query_name, 29911198, pileupread.alignment.query_sequence[pileupread.query_position]))             
                print ('\tbase in read %s = %s' % (pileupread.alignment.query_name, pileupread.alignment.query_sequence[pileupread.query_position]))

mybam.close()

【问题讨论】：

您确定 520817 读取在一个位置吗？听起来相当高。您在 IGV 中从哪里获得价值？我刚刚将值ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG00096/… 与 pysam 和 IGV 进行了比较，它们是相同的。
嗨，Maximilian，我查看了 IGV 的报道轨道。几乎所有的读取都是重叠的，并且具有相似的开始和结束位置。当我获取时，我得到正确的估计。但是，在做堆积时，这不起作用，我得到的阅读量要少得多......
你能在一些公共数据上试试你的代码吗？我之前评论中的 BAM 文件？使用相同的数据进行故障排除要容易得多。

标签： python bioinformatics pysam

【解决方案1】：

检查 IGV 选项 View-->Preference-->Alignment，一些“过滤 xxxx”选项（重复、次要对齐、低质量）可能会改变输出。

通常 pysam 不会使用 BAM_FUNMAP、BAM_FSECONDARY、BAM_FQCFAIL、BAM_FDUP 标志堆积读取，因此请确保您的 IGV 视图选项与 pysam.AlignmentFile.pileup 中的选项相同。否则它们可能会产生不同的输出。

【讨论】：