Определите длину последовательности, используя СИГАРУ

Question

Чтобы дать вам немного контекста: я пытаюсь преобразовать файл sam в bam

samtools view -bT reference.fasta sequences.sam > sequences.bam

который выходит со следующей ошибкой

[E::sam_parse1] CIGAR and query sequence are of different length
[W::sam_read1] parse error at line 102
[main_samview] truncated file

и оскорбительная строка выглядит так:

SRR808297.2571281       99      gi|309056|gb|L20934.1|MSQMTCG   747     80      101M    =       790     142     TTGGTATAAAATTTAATAATCCCTTATTAATTAATAAACTTCGGCTTCCTATTCGTTCATAAGAAATATTAGCTAAACAAAATAAACCAGAAGAACAT      @@CFDDFD?HFDHIGEGGIEEJIIJJIIJIGIDGIGDCHJJCHIGIJIJIIJJGIGHIGICHIICGAHDGEGGGGACGHHGEEEFDC@=?CACC>CCC      NM:i:2  MD:Z:98A1A

Моя последовательность длиной 98 символов, но возможная ошибка при создании файла sam, о котором было сообщено в CIGAR, 101. Я могу позволить себе роскошь потерять пару чтений, и в данный момент у меня нет доступа к исходному коду, который создал файлы sam, поэтому нет возможности выследить ошибку и повторно выполнить выравнивание. Другими словами, мне нужно прагматичное решение, чтобы двигаться дальше (пока). Поэтому я разработал скрипт на python, который подсчитывает длину моей строки нуклеотидов, сравнивает ее с тем, что зарегистрировано в CIGAR, и сохраняет "нормальные" строки в новом файле.

#!/usr/bin/python
import itertools
import cigar

with open('myfile.sam', 'r') as f:
    for line in itertools.islice(f,3,None): #Loop through the file and skip the first three lines
            cigar=line.split("\t")[5]
            cigarlength = len(Cigar(cigar)) #Use module Cigar to obtain the length reported in the CIGAR string
            seqlength = len(line.split("\t")[9])

            if (cigarlength == seqlength):
                    ...Preserve the line in a new file...

Как видите, для перевода СИГАРЫ в целое число, показывающее длину, я использую модуль СИГАР. Если честно, я немного настороженно отношусь к его поведению. Этот модуль, похоже, неправильно рассчитал длину в очень очевидных случаях. Есть ли другой модуль или более явная стратегия для перевода СИГАРЫ в длину последовательности?

Sidenote: Интересно, что, мягко говоря, эта проблема широко освещалась, но в интернете не найдено прагматического решения. Смотрите ссылки ниже:

https://github.com/COMBINE-lab/RapMap/issues/9
http://seqanswers.com/forums/showthread.php?t=67253
http://seqanswers.com/forums/showthread.php?t=21120
https://groups.google.com/forum/#!msg/snap-user/FoDsGeNBDE0/nRFq-GhlAQAJ

Answer 1

Спецификация SAM предлагает нам эту таблицу операций CIGAR, которая указывает, какие из них "потребляют" запрос или ссылку, вместе с явными инструкциями о том, как вычислить длину последовательности из строки CIGAR:

                                                             Consumes  Consumes
Op  BAM Description                                             query  reference
M   0   alignment match (can be a sequence match or mismatch)   yes   yes
I   1   insertion to the reference                              yes   no
D   2   deletion from the reference                             no    yes
N   3   skipped region from the reference                       no    yes
S   4   soft clipping (clipped sequences present in SEQ)        yes   no
H   5   hard clipping (clipped sequences NOT present in SEQ)    no    no
P   6   padding (silent deletion from padded reference)         no    no
=   7   sequence match                                          yes   yes
X   8   sequence mismatch                                       yes   yes

• "Потребляет запрос" и "потребляет ссылку" указывает, вызывает ли операция CIGAR выравнивание по последовательности запроса и эталонной последовательности соответственно.

...

• Сумма длин операций M/I/S/=/X должна равняться длине SEQ.

Это позволяет нам тривиально вычислить длину последовательности из ее CIGAR, сложив длины всех операций "потребляет запрос" в CIGAR. Это именно то, что происходит в модуле сигары (см. https://github.com/brentp/cigar/blob/754cfed348364d390ec1aa40c951362ca1041f7a/cigar.py), поэтому я не знаю, почему OP здесь посчитал реализацию этого модуля был неправ.

Если мы извлечем соответствующий код из (уже очень короткого) модуля сигары, у нас останется что-то вроде этого в виде короткой реализации на Python операции суммирования, описанной в приведенной выше цитате:

from itertools import groupby

def query_len(cigar_string):
    """
    Given a CIGAR string, return the number of bases consumed from the
    query sequence.
    """
    read_consuming_ops = ("M", "I", "S", "=", "X")
    result = 0
    cig_iter = groupby(cigar_string, lambda chr: chr.isdigit())
    for _, length_digits in cig_iter:
        length = int(''.join(length_digits))
        op = next(next(cig_iter)[1])
        if op in read_consuming_ops:
            result += length
    return result

Answer 2

Я подозреваю, что причина, по которой нет инструмента для решения этой проблемы, заключается в том, что нет общего решения, кроме повторного выравнивания с использованием программного обеспечения, которое не демонстрирует эту проблему. В вашем примере последовательность запросов идеально выравнивается по ссылке, и поэтому в этом случае строка CIGAR не очень интересна (всего одна Mоперация atch с префиксом общей длины запроса). В этом случае исправление просто требует изменения 101M в 98M,

Тем не менее, для более сложных строк CIGAR (например, те, которые включают Insertions, Dэлементы или другие операции), вы не сможете узнать, какая часть строки CIGAR слишком длинная. Если вы вычтете из неправильной части строки CIGAR, у вас останется неправильное чтение, что, вероятно, хуже для вашего последующего анализа, чем просто полное считывание.

Тем не менее, если это оказывается тривиальным, чтобы сделать это правильно (возможно, ваша процедура неправильного выравнивания всегда добавляет дополнительные базы к первой или последней операции CIGAR), то вам нужно знать, как правильно рассчитать длину запроса в соответствии с Строка CIGAR, чтобы вы знали, что из нее вычесть.

samtools рассчитывает это с помощью htslib функция bam_cigar2qlen.

Другие функции, которые bam_cigar2qlen вызовы определены в sam.h, включая полезный комментарий, показывающий таблицу истинности, для которой операции используют последовательность запросов против ссылочной последовательности.

Короче говоря, чтобы вычислить длину запроса строки CIGAR так, как это делает samtools (действительно htslib), вы должны добавить данную длину для операций CIGAR. M, I, S, =, или же X и игнорировать длину операций CIGAR для любых других операций.

Текущая версия модуля сигары Python, похоже, использует тот же набор операций и алгоритм для расчета длины запроса (что len(Cigar(cigar)) вернулся бы) выглядит прямо для меня. Что заставляет вас думать, что это не дает правильных результатов?

Похоже, вы должны иметь возможность использовать модуль Python для сигары, чтобы жестко обрезать левый или правый конец, используя mask_left или же mask_right метод с mask="H",