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