Как использовать pysam.view для эмуляции всех функций представления samtools

Question

Я пытаюсь использовать 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 несколько регионов и как один из них определяет этот список? Я искал документацию по этому вопросу, но ничего не нашел.

Answer 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.