python - BAM 文件 : getting all reads on certain position with pysam

Question

我有一个 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()

Answer 1

检查 IGV 选项 View -->首选项-->对齐，某些“过滤 xxxx”选项(重复、辅助对齐、低质量)可能会更改输出。

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