Доступ к файлу Bam в определенном месте с помощью Pysam

Question

У меня есть номер хромосомы и местоположение (chr1 и местоположение 1599812). Я хочу использовать модуль pysam в python для доступа к файлу 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)

Обратите внимание, что это работает, потому что, как указано в глоссарии API pysam, 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))

Answer 2

Если вы установите одинаковые начало и конец, наложение будет относиться только к этой конкретной позиции. Например (чистые samtools):

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

Показывает 102 читает, покрывая положение 808957 для хромосомы 1.