我正在尝试使用pysam.view()从BAM文件中过滤出某些对齐方式。我面临的问题是如何在过滤器中包含多个区域。
pysam.view()模拟samtools view命令,该命令允许一个人输入由空格字符分隔的多个区域,例如:
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是否支持多个区域,一个区域如何指定此列表?我已经搜索了与此相关的文档,但未找到任何内容。
您的问题的简短答案是您要使用的格式是
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中的功能进行进一步的工作。