Глобальное выравнивание: Нидлмана-Вунша

Алгоритм Нидлмана-Вунша находит оптимальное глобальное выравнивание двух последовательностей. Это жемчужина динамического программирования в биологии.

Глобальное выравнивание сопоставляет последовательности по всей длине от начала до конца. Алгоритм Нидлмана-Вунша находит выравнивание с максимальным счётом за O(n·m) методом динамического программирования.

Это, возможно, самый важный алгоритм курса. Он показывает, как динамическое программирование превращает экспоненциальный перебор всех выравниваний в быстрый расчёт по таблице. Разберём по косточкам: матрица, рекуррентность, штрафы, трассировка.

Идея: разбить задачу на подзадачи

Заведём матрицу dp, где dp[i][j] — лучший счёт выравнивания первых i букв строки A с первыми j буквами строки B. В каждой клетке три варианта: буквы выровнены друг с другом (диагональ), буква A против пропуска (вверх), буква B против пропуска (влево). Берём максимум.

dp[i][j] = max(
   dp[i-1][j-1] + s(A[i], B[j]),   # выровнять буквы (совпадение/замена)
   dp[i-1][j]   + gap,             # пропуск в B (буква A против '-')
   dp[i][j-1]   + gap              # пропуск в A (буква B против '-')
)

Строим матрицу счётов

Первая строка и столбец — накопленные штрафы за пропуски (выровнять начало с пустотой). Дальше заполняем по рекуррентности.

def nw_matrix(a, b, match=1, mismatch=-1, gap=-1):
    n, m = len(a), len(b)
    dp = [[0] * (m + 1) for _ in range(n + 1)]
    for i in range(1, n + 1):
        dp[i][0] = i * gap
    for j in range(1, m + 1):
        dp[0][j] = j * gap
    for i in range(1, n + 1):
        for j in range(1, m + 1):
            s = match if a[i-1] == b[j-1] else mismatch
            dp[i][j] = max(dp[i-1][j-1] + s, dp[i-1][j] + gap, dp[i][j-1] + gap)
    return dp

a, b = 'GATTA', 'GAATA'
dp = nw_matrix(a, b)
labels = '-' + a
print('     -  ' + '  '.join(b))
for i, row in enumerate(dp):
    print(labels[i] + ' ' + ' '.join(f'{v:3d}' for v in row))
print('Итоговый счёт:', dp[-1][-1])

Вывод:

     -  G  A  A  T  A
-   0  -1  -2  -3  -4  -5
G  -1   1   0  -1  -2  -3
A  -2   0   2   1   0  -1
T  -3  -1   1   1   2   1
T  -4  -2   0   0   2   1
A  -5  -3  -1   1   1   3
Итоговый счёт: 3

Число в правом нижнем углу (3) — счёт лучшего выравнивания. Каждая клетка заполнялась за O(1), вся матрица — за O(n·m). Так перебор всех экспоненциально многих выравниваний свёлся к заполнению таблицы.

Трассировка: восстановить само выравнивание

Матрица даёт счёт, но нам нужно и выравнивание. Идём из правого нижнего угла назад, выбирая, откуда пришёл максимум, и собираем строки с пропусками.

def needleman_wunsch(a, b, match=1, mismatch=-1, gap=-1):
    n, m = len(a), len(b)
    dp = [[0] * (m + 1) for _ in range(n + 1)]
    for i in range(1, n + 1):
        dp[i][0] = i * gap
    for j in range(1, m + 1):
        dp[0][j] = j * gap
    for i in range(1, n + 1):
        for j in range(1, m + 1):
            s = match if a[i-1] == b[j-1] else mismatch
            dp[i][j] = max(dp[i-1][j-1] + s, dp[i-1][j] + gap, dp[i][j-1] + gap)
    # трассировка
    al_a, al_b = '', ''
    i, j = n, m
    while i > 0 or j > 0:
        s = match if (i > 0 and j > 0 and a[i-1] == b[j-1]) else mismatch
        if i > 0 and j > 0 and dp[i][j] == dp[i-1][j-1] + s:
            al_a = a[i-1] + al_a; al_b = b[j-1] + al_b; i -= 1; j -= 1
        elif i > 0 and dp[i][j] == dp[i-1][j] + gap:
            al_a = a[i-1] + al_a; al_b = '-' + al_b; i -= 1
        else:
            al_a = '-' + al_a; al_b = b[j-1] + al_b; j -= 1
    return dp[n][m], al_a, al_b

score, x, y = needleman_wunsch('GCATGCU', 'GATTACA')
print('Счёт:', score)
print(x)
print(''.join('|' if (p == q and p != '-') else ' ' for p, q in zip(x, y)))
print(y)

Вывод:

Счёт: 0
GCA-TGCU
| | | | 
G-ATTACA

Классический учебный пример (его приводил сам Нидлман): счёт 0, в выравнивании видны совпадения (|), замены и пропуски (-). Дефис означает вставку/удаление: чтобы выровнять, в одну строку вставлен пропуск.

Как работает под капотом: смысл штрафов

Параметры match, mismatch, gap кодируют биологию. Награда за совпадение тянет к идентичности; штраф за замену — за мутации; штраф за пропуск — за вставки/удаления. Меняя их, мы управляем характером выравнивания: большой штраф за пропуск даёт компактные выравнивания без дыр. Для белков вместо простого match/mismatch используют матрицы замен (BLOSUM/PAM) — следующий урок. Сложность O(n·m) делает алгоритм точным, но тяжёлым для длинных геномов — отсюда эвристики вроде BLAST.

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

Забыть инициализацию первой строки/столбца. Без накопленных штрафов за пропуски выравнивание неверное.
Перепутать индексы. a[i-1], b[j-1] — матрица на 1 больше длины строк (нулевая строка/столбец).
Применять глобальное выравнивание к участкам разной длины. Если похож лишь фрагмент — нужно локальное (Смит-Ватерман).

Итог

Нидлман-Вунш — глобальное выравнивание по всей длине через динамическое программирование.
dp[i][j] — лучший счёт первых i и j букв; рекуррентность из трёх вариантов (диагональ, вверх, влево).
Матрица даёт счёт за O(n·m); трассировка из угла восстанавливает само выравнивание.
Штрафы кодируют биологию мутаций и indel; для длинных последовательностей нужны эвристики.