使用Pysam访问特定位置的Bam文件

Question

我有一个给定的染色体数目和位置(chr1和位置1599812)。我想使用python的pysam模块访问bam文件，以仅获得该特定区域chr1和位置1599812的读取编号信息。我尝试过使用pileup()，但它需要一个位置范围，而在我的例子中，我只想要一个特定的位置，而不是一个这样的范围。

Answer 1

我不认为pileup()是你想要的--根据pysam API的说法，这个函数返回“一个关于基因组位置的迭代器”，特别是，“所有与区域重叠的读取都会被返回。返回的第一个碱基将是第一个读取的第一个碱基，而不一定是查询中使用的区域的第一个碱基。”

你是说你想要获得“读数信息”--也就是那个特定位置的读数，对吗？为此，count_coverage()应该完成这项工作。在你的例子中，我认为这段代码应该会给你你想要的答案：

import pysam

my_bam_file = '/path/to/your/bam_file.bam'
imported = pysam.AlignmentFile(my_bam_file, mode = 'rb')  # 'rb' ~ read bam
coverage = imported.count_coverage(
                  contig = '1',     # Chromosome ID; also might be "chr1" or similar 
                  start = 1599812,
                  stop = 1599813,
                  )
print(coverage)

请注意，这是可行的，因为正如pysam API glossary中所指出的，

使用半开区间，因此范围[1599812,1599813]将恰好包括一个碱基对。

运行上面的代码会得到类似下面这样的结果：

> (array('L', [0]), array('L', [0]), array('L', [0]), array('L', [0]))

其是分别包含覆盖该基因组位置的读物中的A、C、G和T碱基的数目的阵列的元组。如果您只对映射到此特定基因组位置的读取总数感兴趣，则可以对此元组求和：

import numpy as np

print(np.sum(coverage))

Answer 2

如果您设置相同的开始和结束，则堆积将仅引用该特定位置。例如(纯samtools)：

$ samtools mpileup -r chr1:808957-808957 YourFile.bam
chr1    808957  N   102 READSTRING READQUALITYSTRING

显示了覆盖1号染色体位置808957的第102个阅读片段。