Индексация генома по k-мерам

Чтобы искать в геноме много раз, его индексируют по k-мерам — как поисковик индексирует веб. Запрос отвечается почти мгновенно.

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

В уроке про поиск мотивов мы выяснили: искать каждый запрос наивно по огромному геному нереально. Решение — заплатить один раз за построение индекса, а потом отвечать на запросы за доли секунды. Именно так устроены выравниватели ридов вроде Bowtie и BWA (хоть и с более хитрыми структурами).

Строим k-мерный индекс

Проходим геном один раз и для каждого k-мера запоминаем его позицию. Получаем словарь «k-мер → список позиций».

def build_index(seq, k):
    index = {}
    for i in range(len(seq) - k + 1):
        kmer = seq[i:i+k]
        index.setdefault(kmer, []).append(i)
    return index

genome = 'ATGATGGCAGCA'
idx = build_index(genome, 3)
print('ATG встречается в позициях:', idx['ATG'])
print('GCA встречается в позициях:', idx['GCA'])
print('Уникальных 3-меров:', len(idx))

Вывод:

ATG встречается в позициях: [0, 3]
GCA встречается в позициях: [6, 9]
Уникальных 3-меров: 8

Поиск запроса через индекс

Чтобы найти, где в геноме встречается запрос, берём его первый k-мер, мгновенно получаем кандидатов-позиции из индекса, а затем проверяем полное совпадение только в этих местах. Это «фильтр + проверка».

def build_index(seq, k):
    index = {}
    for i in range(len(seq) - k + 1):
        index.setdefault(seq[i:i+k], []).append(i)
    return index

def query(genome, idx, k, pattern):
    seed = pattern[:k]
    hits = []
    for pos in idx.get(seed, []):
        if genome[pos:pos+len(pattern)] == pattern:
            hits.append(pos)
    return hits

genome = 'ATGGGATGCATT'
idx = build_index(genome, 3)
print('Поиск ATGCA:', query(genome, idx, 3, 'ATGCA'))

Вывод:

Поиск ATGCA: [5]

Seed ATG встречается в этом геноме дважды (позиции 0 и 5), но полное слово ATGCA подтвердилось лишь в позиции 5. Мы проверили только две позиции-кандидата вместо всего генома — на больших данных это разница между секундами и часами.

Почему это быстро

Построение индекса — O(n) (один проход). Зато каждый запрос — это поиск в словаре O(1) плюс проверка лишь нескольких кандидатов. Если запросов много (миллионы ридов), амортизированная стоимость крошечная. Сравните: наивно — O(запросов · n), с индексом — O(n + запросов · кандидатов).

Без индекса:  каждый запрос сканирует ВЕСЬ геном
С индексом:   запрос -> словарь -> 2-3 позиции -> проверка
              (геном просканирован ОДИН раз при постройке)

Как работает под капотом: цена памяти

За скорость платят памятью: индекс хранит позиции всех k-меров и может быть в разы больше самого генома. Поэтому реальные инструменты используют экономные структуры: FM-индекс на основе преобразования Барроуза-Уилера (BWT) сжимает индекс почти до размера генома, сохраняя быстрый поиск. Идея k-мерного индекса — упрощённая модель того же принципа: предобработка ради быстрых запросов. Понимая её, вы понимаете суть всех быстрых выравнивателей.

Полезно соотнести наш индекс с тремя классическими структурами, которые вы могли встречать в курсе структур данных. Хеш-таблица (наш словарь) проста и быстра в среднем, но прожорлива по памяти. Суффиксный массив хранит отсортированные позиции всех суффиксов и позволяет искать любую подстроку бинарным поиском, занимая меньше места. FM-индекс идёт ещё дальше — сжимает данные так, что индекс сопоставим по размеру с самим геномом, при этом поддерживая поиск. Все три решают одну задачу — быстрый повторный поиск в фиксированном тексте — но по-разному балансируют память и скорость. Наш k-мерный словарь — самый наглядный вход в эту тему: он показывает идею «предобработай один раз, ищи мгновенно», а какую именно структуру выбрать, зависит от того, сколько у вас памяти и насколько длинные запросы.

Стоит проговорить и связь с приближённым поиском из предыдущего раздела. Точный индекс находит только точные совпадения seed, но риды содержат ошибки и мутации. Трюк, которым пользуются выравниватели: рид содержит несколько seed-ов, и хотя бы один из них при небольшом числе ошибок совпадёт точно (тот самый «принцип ящиков»). Найдя точное совпадение хотя бы одного seed, инструмент уже знает примерное место в геноме и там делает дорогое точное выравнивание с допуском ошибок (его мы разберём в следующем разделе). Получается двухступенчатая схема: быстрый k-мерный фильтр грубо находит район, медленное выравнивание уточняет. Эта же идея — основа BLAST. Так k-мерный индекс оказывается не отдельным трюком, а первой ступенью почти всех быстрых методов сравнения последовательностей.

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

Строить индекс на каждый запрос. Смысл в том, чтобы построить один раз и переиспользовать.
Забыть проверку полного совпадения. Seed-попадание — лишь кандидат; полное слово надо подтвердить.
Недооценить память. На больших геномах простой словарь не влезет — нужны компактные индексы.

Итог

k-мерный индекс — словарь «k-мер → позиции», строится за один проход O(n).
Запрос ищет seed в индексе за O(1) и проверяет лишь нескольких кандидатов.
Это превращает многократный поиск из неподъёмного в мгновенный.
Реальные выравниватели используют сжатые индексы (BWT/FM) ради экономии памяти.