Поиск подстроки и мотивов

Поиск короткой последовательности-мотива внутри длинной — одна из самых частых операций. Сделаем её правильно.

Мотив — короткая последовательность с биологическим смыслом (сайт связывания белка, рестрикционный сайт, сигнал старта), которую ищут в геноме.

Биологи постоянно спрашивают: «где в этом гене сайт связывания транскрипционного фактора?», «сколько раз встречается этот мотив?». Это поиск подстроки — но с нюансами, из-за которых наивный str.find может соврать. Разберёмся.

Все вхождения, а не первое

seq.find(motif) даёт только первую позицию. Нам нужны все. Идём в цикле, сдвигая старт поиска.

def find_all(seq, motif):
    positions = []
    i = seq.find(motif)
    while i != -1:
        positions.append(i)
        i = seq.find(motif, i + 1)
    return positions

dna = 'ATGATGCATGATG'
print('Позиции ATG:', find_all(dna, 'ATG'))

Вывод:

Позиции ATG: [0, 3, 7, 10]

Перекрывающиеся вхождения — частая ловушка

Сдвиг на i + 1 важен: он находит и перекрывающиеся вхождения. Метод str.count и split их пропускают. Для мотива ATA в ATATA это критично.

def find_all(seq, motif):
    positions = []
    i = seq.find(motif)
    while i != -1:
        positions.append(i)
        i = seq.find(motif, i + 1)
    return positions

seq = 'ATATATA'
print('count (не видит перекрытий):', seq.count('ATA'))
print('find_all (видит перекрытия):', find_all(seq, 'ATA'))

Вывод:

count (не видит перекрытий): 2
find_all (видит перекрытия): [0, 2, 4]

В ATATATA мотив ATA встречается три раза (позиции 0, 2, 4), но они перекрываются, и count насчитал лишь 2. В биологии перекрывающиеся сайты реальны, поэтому почти всегда нужен поиск со сдвигом на 1.

Сложность и большие геномы

Наивный поиск в худшем случае $O(nm)$, где $n$ — длина текста, $m$ — длина мотива. Для одного короткого мотива в геноме это приемлемо. Но если мотивов миллионы (например, картирование ридов), наивный подход умрёт — там применяют индексы (хеши k-меров, суффиксные структуры), о которых речь в разделе про k-меры. Питоновский find внутри использует быстрый алгоритм, так что для учебных и многих практических задач его достаточно.

Полезно прочувствовать масштаб. Геном бактерии E. coli — около 4,6 миллиона букв. Один проход поиска короткого мотива по нему — это несколько миллионов операций, доли секунды. Но представьте задачу картирования: у вас 30 миллионов ридов, и каждый надо найти в геноме. Наивно это 30 000 000 поисков по 4 600 000 букв — десятки триллионов операций, что нереально. Именно поэтому для массового поиска геном предварительно индексируют один раз, а потом каждый запрос отвечается почти мгновенно. Разница между «искать каждый раз заново» и «построить индекс однажды» — водораздел между учебным кодом и инструментом, который реально гоняют на данных секвенатора.

Как работает под капотом: почему не регуляркой всегда

Можно искать и через re, но для точного фиксированного мотива это избыточно и иногда медленнее. Регулярки незаменимы, когда мотив неоднозначен (например, «G, потом A или C, потом T») — этим займёмся в следующем уроке. Для точного совпадения связка find в цикле проще и предсказуемее.

Текст:  A T A T A T A
        ^ATA
            ^ATA   (перекрытие!)
                ^ATA

Что вообще считается «мотивом» в биологии и почему его поиск так важен — стоит проговорить. Мотив обычно отражает место, к которому что-то прикрепляется или что-то происходит: сайт посадки белка-регулятора, точка разреза фермента, сигнал начала или конца гена, участок, узнаваемый системой защиты клетки. Найдя все вхождения мотива, биолог отвечает на содержательные вопросы: активен ли ген (есть ли перед ним нужные регуляторные сайты), уязвима ли ДНК для конкретного фермента, не вирусная ли это последовательность. То есть за сухим поиском подстроки стоит биологическая интерпретация, и именно она превращает позиции [0, 3, 7, 10] в осмысленный вывод. Хороший биоинформатик никогда не останавливается на «нашёл столько-то совпадений» — он спрашивает, что эти совпадения значат.

Частые ошибки

• Брать только find без цикла. Найдёте лишь первое вхождение.
• Сдвигать на длину мотива. Тогда пропустите перекрывающиеся вхождения; сдвигайте на 1.
• Использовать count для подсчёта мотивов. Он не считает перекрытия.

Итог

• Для всех вхождений мотива ищите в цикле, сдвигая старт на 1.
• Сдвиг на 1 ловит перекрывающиеся вхождения, которые count и split теряют.
• Наивный поиск $O(nm)$ годится для одного мотива; для массового поиска нужны индексы.
• Точный мотив ищут find-ом, неоднозначный — регулярками.