Поиск в ДНК — удивительно важная задача биоинформатики. Некоторые генетические цепи, например, используют нити РНК для передачи сигналов через клетку. Было бы очень плохо, если бы нити РНК, уже присутствующие в клетке, могли мешать этим цепям, поэтому важно обнаруживать перекрытия при разработке схемы.
Вот тут-то и появляется мой алгоритм. Мой алгоритм основан на этом алгоритме поиска строк, но с некоторыми модификациями, чтобы заставить его работать с ДНК. Мой алгоритм имеет сложность в наихудшем случае 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.