如何使用 pysam.view 模拟 samtools 视图的所有功能答案

【问题标题】：How to use pysam.view to emulate all functions of samtools view如何使用 pysam.view 模拟 samtools 视图的所有功能
【发布时间】：2016-06-06 13:59:02
【问题描述】：

我正在尝试使用 pysam.view() 从 BAM 文件中过滤掉某些对齐。我面临的问题是如何在过滤器中包含多个区域。

pysam.view() 模拟 samtools 视图命令，该命令允许进入由空格字符分隔的多个区域，例如：

samtools view opts bamfile chr1:2010000-20200000 chr2:2010000-20200000

但是对应的pysam.view调用：

pysam.view(ops, bamfile, '1:2010000-20200000 2:2010000-20200000')

不起作用。它不返回任何对齐方式。我很确定问题在于如何指定区域列表，因为以下命令可以正常工作：

pysam.view(ops, bamfile, '1:2010000-20200000')

并返回对齐方式。

我的问题是：pysam.view 是否支持多个区域以及如何指定此列表？我已经搜索了有关此的文档，但没有找到任何东西。

【问题讨论】：

我应该补充一点，我已经尝试了最明显的指定列表的方法：使用逗号、分号、句点、制表符作为分隔符并将区域放在列表和元组中。

标签： samtools pysam

【解决方案1】：

对您的问题的简短回答是您使用的格式是

pysam.view(ops, bamfile, '1:2010000-20200000','2:2010000-20200000')

（另请注意，表示您的每个区域结束的数字比开头大约 10 倍 - 看来您可能打算改为 2010000-2020000。）

我已经使用以下代码对其进行了测试：

import pysam

my_bam_file = '/path/to/my/bam_file.bam'
alignments1 = pysam.view(my_bam_file, '1:2010000-4000000')
alignments2 = pysam.view(my_bam_file, '1:5000000-6000000')
alignments3 = pysam.view(my_bam_file, '1:2010000-4000000', '1:5000000-6000000')

print(len(alignments1) + len(alignments2) == len(alignments3))

[Output:] True

但是，这种提取对齐方式的效率不是很高，因为您得到的输出是一个大的str，而不是单独的对齐方式。要获得单独对齐的list，请使用以下代码：

import pysam

my_bam_file = '/path/to/my/bam_file.bam'
imported = pysam.AlignmentFile(my_bam_file, mode = 'rb')
regions = ('1:2010000-20200000','2:2010000-20200000')
alignments = []
for region in regions:
    bam = imported.fetch(region = region, until_eof = True)
    alignments.extend([alignment for alignment in bam])

alignment 的每个元素最终都会成为一个pysam.AlignedSegment 对象，您可以使用pysam API 中的函数进一步处理它。

【讨论】：