DNA Sekwencjonowanie: Aligmentacja Smith-Watermana i Wykrywanie Zmiennych Pojedyncznych

Sekwencjonowanie produkuje fragmenty, a nie genom

Nowoczesne sekwencery o krótkich odczytach nie czytają chromosomu od początku do końca. Rozkładają wiele kopii DNA na miliony nakładających się fragmentów – "odczytów" – zwykle długości 100 do 300 baz, a następnie czytają każdy z nich niezależnie i niedokładnie, z błędami w około 0,1 do 1 procentu baz, w zależności od platformy. Przebudowa tego stosu krótkich, podatnych na błędy fragmentów w spójną sekwencję genomu lub w różnicę względem znanego wzorca stanowi w pełni problem dopasowania (alignment), a dopasowanie to jest miejscem, gdzie żyje prawdziwa informatyka.

Algorytm Smith-Watermana: znalezienie najlepszego lokalnego dopasowania

Algorytm Smith-Watermana (1981) jest metodą programowania dynamicznego, która znajduje najlepsze lokalne dopasowanie między dwiema sekwencjami – czyli obszar, w którym elementy się najbardziej zgadzają, pomijając te, które nie. Buduje się go krok po kroku, tworząc macierz ocen, gdzie każda komórka reprezentuje najlepszą możliwą ocenę dopasowania kończącego się na danym parze pozycji:

H(i,j) = max( 0, // reset – to jest sedno sprawy: lokalność H(i-1, j-1) + ocena(a[i], b[j]), // dopasowanie lub niezgodność, diagonalna H(i-1, j) - kara_za_przerwę, // usunięcie, z góry H(i, j-1) - kara_za_przerwę // wstawienie, z lewej strony ) Termin max(0, …) jest kluczowy: pozwala on na zresetowanie złej oceny do zera zamiast roznoszenia negatywnych wartości, co naturalnie prowadzi algorytm do odkrycia początku prawdziwego dopasowania. Śledź drogę wstecz od najwyżej ocenionej komórki w macierzy, zgodnie z kierunkiem, z którego pochodziła dana ocena, aż dojdziesz do zerowej wartości – wtedy odzyskasz najlepszą podsekwencję dopasowań. Oceny dopasowania/niezgodności zwykle pochodzą z macierzy substytucji (prosty +1/-1 dla DNA lub BLOSUM dla białek), a model przerwy uwzględnia wstawienia lub usunięcia, które rzeczywiście występują w DNA podczas błędów replikacji).

H(i,j) = max(
  0,                                       // reset — this is what makes it LOCAL
  H(i-1, j-1) + score(a[i], b[j]),         // match or mismatch, diagonal
  H(i-1, j)   - gapPenalty,                // deletion, from above
  H(i,   j-1) - gapPenalty                 // insertion, from the left
)

Dlaczego nie globalna wersja?

Smith-Waterman to zmodyfikowana wariantem jednej linii algorytmu Needleman-Wunsza (1970), który znajduje najlepsze globalne dopasowanie – zmuszając dwie sekwencje do ustawienia się obok siebie, karząc niezgodnym fragmentom. To właściwy narzut do porównywania dwóch pełnych sekwencji o zbliżonej długości, np. dwóch wersji tego samego genu. Jest to jednak niewłaściwy narzut do mapowania 150-bazowego odczytu przeciwko odniesieniu o 3 miliardach bazowych, gdzie odczyt odpowiada niewielkiej sekwencji, a reszta jest nieistotna – dokładnie to, co algorytm Smith-Watermana z jego lokalnym, resetowalnym wynikiem potrafi prawidłowo obsłużyć.

Od wyrównania do identyfikacji wariantu

Monoczęstotliwość zmienność (SNP) to pojedyncza zmiana podjednokropki względem odniesienia w danym położeniu. Po ułożeniu milionów odczytów, identyfikacja SNP jest decyzją statystyczną, a nie pojedynczym wyszukiwaniem: genom sekwencjonowany przy typowej głębokości klinicznej jest pokryty 20 do 30 niezależnymi odczytami w każdym położeniu, a prawdziwy SNP wykazuje konsekwentne zastąpienie w wielu z nich przy rozsądnym udziale alleli (około 50% dla jednej heterozygotycznej kopii, blisko 100% dla homozygotycznego zmiany). Z drugiej strony, błąd sekwencjonowania jest zasadniczo losowy na każdym odczycie i rzadko się powtarza w wielu niezależnych odczytach tej samej podjednokropki — co dokładnie stanowi sygnał statystyczny, którego używają oprogramowanie do identyfikacji wariantu, aby rozróżnić prawdziwą biologię od szumu maszynowego.

Skalowanie: pełna analiza genomowa w praktyce

Uruchamianie pełnego Smitha-Watermana – algorytmu o złożoności O(n·m) dynamicznego programowania – 150-bazowego odczytu przeciwko genomowi liczącemu 3 miliardy baz dla każdego z setek milionów odczytów byłoby obliczeniowo niemożliwe. Zamiast tego, rzeczywiste algorytmy (BWA, Bowtie) wykorzystują indeks FM, skompresowany indeks pełnego tekstu zbudowany na podstawie transformacji Burrowsa-Whelera referencyjnego genomu, aby natychmiast zawęzić odczyt do kilkudziesięciu potencjalnych lokalizacji w czasie zbliżonym do O(długości odczytu). Następnie pełne scoringowanie Smitha-Watermana jest wykonywane lokalnie, wokół tych kilku potencjalnych pozycji, aby uzyskać dokładne, zgapowane dopasowanie i ocenę jakości – łącząc indeks prawie idealnego dopasowania dla szybkości z dynamicznym programowaniem w celu uzyskania precyzji tam, gdzie jest to potrzebne.

Często zadawane pytania

Wypróbuj na żywo

Wszystko powyżej działa bezpośrednio w Twojej przeglądarce — otwórz DNA Sequencing i zmieniaj parametry podczas działania. Nic nie jest instalowane ani przesyłane na serwer, cały model działa w jednej karcie.