python - 使用 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)

Note that this works because, as noted in the pysam API glossary, pysam uses half-open intervals, so the range [1599812, 1599813) will include exactly one base-pair.

运行上面的代码会给你这样的结果:

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

这是一个数组元组，分别包含覆盖该基因组位置的读数中 A、C、G 和 T 碱基的数量。如果您只是对映射到该特定基因组位置总数的读取数量感兴趣，则可以对该元组求和:

import numpy as np

print(np.sum(coverage))