使用 Pysam 在特定位置访问 Bam 文件答案

【问题标题】：Accessing Bam file at particular location using Pysam使用 Pysam 在特定位置访问 Bam 文件
【发布时间】：2019-11-19 05:55:43
【问题描述】：

我有一个给定的染色体编号和位置（chr1 和位置 1599812）。我想使用 python 的 pysam 模块来访问 bam 文件以获取仅特定区域 chr1 和位置 1599812 的读取数字信息。我尝试使用 pileup() 但它需要一系列位置，而在我的情况下我只想要一个特定的位置，而不是这样的范围。

【问题讨论】：

尝试添加更多标签，以便可能能够帮助您找到问题的人更容易找到您的问题
你的 bam 被索引了吗？
您是否尝试将开始和结束索引设置为相同的坐标（警告基于 0，而不是基于残基）

标签： python bioinformatics python-module biopython pysam

【解决方案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 中所述， pysam 使用半开区间，因此范围 [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))

【讨论】：

【解决方案2】：

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

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

显示 102 个读数，覆盖 1 号染色体的 808957 位置。

【讨论】：