Поиск в ДНК — удивительно важная задача биоинформатики. Некоторые генетические цепи, например, используют нити РНК для передачи сигналов через клетку. Было бы очень плохо, если бы нити РНК, уже присутствующие в клетке, могли мешать этим цепям, поэтому важно обнаруживать перекрытия при разработке схемы.

Вот тут-то и появляется мой алгоритм. Мой алгоритм основан на этом алгоритме поиска строк, но с некоторыми модификациями, чтобы заставить его работать с ДНК. Мой алгоритм имеет сложность в наихудшем случае O(N) и сложность в наилучшем случае O(M), когда присутствует искомая последовательность ДНК, и сложность в наилучшем случае O(N/M), когда нити нет, где M — длина подстроки.

Шаги алгоритма

Предварительная обработка и настройка

Прежде всего, входные данные должны быть обработаны для алгоритма. Последовательности ДНК преобразуются в байты, где A равно 0b10, T равно 0b11, C равно 0b00, а G равно 0b01.

Затем мы создаем кеш, который содержит подпоследовательности длиной 4 п.н., присутствующие в искомой последовательности. Алгоритм создает массив из 256 «сегментов», к которым мы можем получить доступ с помощью битового сдвига. Алгоритм перебирает последовательность, которую мы ищем, и добавляет подпоследовательность в этом месте в массив, используя следующий формат.

cache[(sequence[i] << 6) + (sequence[i + 1] << 4) + ...] = true

Этот кеш позволит нам обнаруживать наличие фрагментов размером 4 бита за постоянное время. Этот шаг предварительной обработки — самая медленная часть алгоритма, но мы можем кэшировать результаты в двоичном файле, который можно очень быстро загрузить во время выполнения.

Основной алгоритм

Алгоритм начинается с выравнивания подпоследовательности (то, что мы ищем) с более длинной последовательностью, например:

Seq1:   AAAAATTCGCGAT
SubSeq: ATTCG

Затем алгоритм входит в основной цикл. В этом цикле алгоритм использует кеш, чтобы проверить, присутствует ли в последовательности текущий фрагмент из 4 бит.

Seq1:   A[AAAA]TTCGCGAT
SubSeq: ATTCG

Если фрагмент из 4 битов отсутствует, алгоритм продвигается вперед на длину подпоследовательности минус 1.

Seq1:   AAAAATTCGCGAT
SubSeq:     ATTCG

Затем алгоритм повторяет этот процесс, пока не достигнет конца или не найдет совпадение. Если он находит совпадение, он перемещается назад, пока не найдет место совпадения в последовательности. Вот пример:

      Found match↓
Seq1:   AAAAATTCGCGAT
SubSeq:    TCGCGAT
Shift forwards by 1:
      Found match↓
Seq1:   AAAAATTCGCGAT
SubSeq:     TCGCGAT
Shift forwards by 1:
      Found match↓
Seq1:   AAAAATTCGCGAT
SubSeq:      TCGCGAT
Shift forwards by 1, sequences overlap:
      Found match↓
Seq1:   AAAAATTCGCGAT
SubSeq:       TCGCGAT

Теперь, когда алгоритм нашел совпадение, он движется назад по подпоследовательности и сравнивает две последовательности, пока либо не найдет несоответствие, либо не достигнет конца подстроки:

Mismatch example:
             Matches↓
Seq1:   AAAAATTCGCGAT
SubSeq:       TCGCCAT
Shift backwards by 1:
            Matches↓
Seq1:   AAAAATTCGCGAT
SubSeq:       TCGCCAT
Shift backwards by 1, mismatch detected:
     Doesn't match↓
Seq1:   AAAAATTCGCGAT
SubSeq:       TCGCCAT

Если алгоритм достигает конца подстроки, то точное совпадение найдено и функция возвращает значение. Если он находит несоответствие, он продвигается на расстояние несоответствия от текущего индекса:

Doesn't match↓
Seq1:   AAAAATTCGCGAT
SubSeq:       TCGCCAT
Advances:
Seq1:   AAAAATTCGCGAT
SubSeq:          TCGCCAT

Как только алгоритм достигает конца последовательности ДНК, он возвращается и завершается.

Алгоритм есть! Я опубликовал его на GitHub под лицензией MIT, так что не стесняйтесь включать его в свои собственные проекты. Вот ссылка! Я провел очень неформальный бенчмаркинг с помощью Visual Studio, и этому алгоритму требуется менее 1 мс, чтобы найти уникальную цепь РНК, отсутствующую в транскриптоме человека.

Привет 👋

Большое спасибо за то, что прочитали эту статью. Я надеюсь, что это было не слишком сложно, и я надеюсь, что это будет полезно для кого-то. Вполне возможно, что я не первый, кто это делает, я просто собрал это за выходные.

Если вы хотите быть в курсе других моих проектов, подобных этому, подпишитесь на мою рассылку и следуйте за мной в Twitter.