SEKWENCJONOWANIE I ASEMBLACJA DNA(2).pdf

Tom 58 2009

Numer 1–2 (282–283)

Strony 17–28

M arta K asprzaK 1,2 , a leKsandra Ś wiercz 1,2

1 Instytut Informatyki,

Politechnika Poznańska

Piotrowo 2, 60-965 Poznań

2 Instytut Chemii Bioorganicznej PAN

Noskowskiego 12/14, 61-704 Poznań

e-mail: mkasprzak@cs.put.poznan.pl

aswiercz@cs.put.poznan.pl

SEKWENCJONOWANIE I ASEMBLACJA DNA — PODEJŚCIA, MODELE GRAFOWE,

ALGORYTMY

WPROWADZENIE

Kwas deoksyrybonukleinowy (DNA) jest

cząsteczką kodującą informację genetyczną

organizmu, złożoną z dwóch nici połączo-

nych ze sobą za pomocą wiązań wodoro-

wych (tzw. helisa DNA). Każda nić to łań-

cuch nukleotydów, których kolejność stano-

wi tę informację, zapisywanych symbolicznie

jako ‘A’, ‘C’, ‘G’, ‘T’. ‘A’ oznacza nukleotyd z

zasadą azotową adeniną, ‘C’ z cytozyną, ‘G’ z

guaniną, ‘T’ z tyminą. Nukleotyd A z jednej

nici helisy łączy się za pomocą wiązania wo-

dorowego z nukleotydem T położonym na

przeciwko niego w drugiej nici, natomiast

nukleotyd G łączy się z C. Ta zasada wystę-

powania w parach A z T oraz G z C nazywa-

na jest zasadą komplementarności. Dzięki tej

właściwości, znając fragment jednej nici, mo-

żemy odtworzyć komplementarną do niego

sekwencję nukleotydów w drugiej nici. Krót-

ki fragment jednoniciowego DNA nazywany

jest oligonukleotydem.

Odczytywanie kolejności nukleotydów w

sekwencji DNA (czyli rozpoznanie informacji

genetycznej organizmu) odbywa się w kilku

etapach. W pierwszym etapie, sekwencjono-

waniu, odczytywane są sekwencje DNA o

długości zazwyczaj kilkuset nukleotydów. W

kolejnym etapie, asemblacji, sekwencje te są

łączone w dłuższe odcinki, w wyniku czego

otrzymujemy sekwencje o długości nawet

do miliona nukleotydów. Dla krótkich geno-

mów, np. wirusów lub bakterii, te dwa etapy

często wystarczą, aby odczytać cały genom.

Dla dłuższych, potrzebny jest trzeci etap, ma-

powanie, polegający na właściwym uszerego-

waniu względem siebie zasemblowanych se-

kwencji. W tym celu stosuje się inne niż w

poprzednich etapach podejścia biochemicz-

ne, np. z użyciem enzymów restrykcyjnych

(mapowanie restrykcyjne, patrz np. B łaże -

wicz i współaut. 2001).

Etap sekwencjonowania może być reali-

zowany na kilka sposobów. Najbardziej po-

pularną metodą do niedawna była elektrofo-

reza żelowa (M axaM i G ilBert 1977, s anGer

i współaut. 1977), która generuje sekwencje

o długości kilkuset nukleotydów. Metoda ta

ogranicza się do etapu laboratoryjnego (eks-

peryment biochemiczny) i nie wymaga etapu

algorytmicznego (przetwarzanie danych), co

z jednej strony czyni ją prostą i przystępną,

z drugiej strony jednak nieodporną na błędy

eksperymentalne.

Bardziej zaawansowanym technologicznie

i koncepcyjnie podejściem jest sekwencjono-

wanie przez hybrydyzację (SBH), w wyniku

którego otrzymujemy sekwencje o długości

do tysiąca nukleotydów, ale w tym przypad-

ku potrzebne są metody algorytmiczne do

przetwarzania danych eksperymentalnych

(s outhern 1988). Na mikromacierz DNA na-

nosi się bibliotekę oligonukleotydów, czyli

zbiór krótkich, jednoniciowych fragmentów

łańcucha DNA (np. wszystkie 4 l oligonukle-

M arta K asprzaK , a leKsandra Ś wiercz

otydy o zadanej długości l ). Podczas ekspe-

rymentu hybrydyzacyjnego badana sekwen-

cja DNA przykleja się do komplementarnych

oligonukleotydów na mikromacierzy. W jego

wyniku znajdowany jest zbiór S (spektrum)

fragmentów, które zawarte są w oryginal-

nej sekwencji (szerszy opis eksperymentu z

mikromacierzą DNA zawarty jest np. w (s o -

uthern i współaut. 1992)). Celem części obli-

czeniowej SBH jest rekonstrukcja oryginalnej

sekwencji o znanej długości n z elementów

spektrum. Rekonstrukcja polega na takim

uszeregowaniu wszystkich elementów spek-

trum, aby kolejne oligonukleotydy w rozwią-

zaniu nakładały się na siebie z przesunięciem

równym 1 (przy założeniu braku błędów eks-

perymentalnych, patrz Ryc. 1).

rodzajem błędów negatywnych są powtórze-

nia fragmentów sekwencji DNA, które są co

najmniej tak długie jak oligonukleotydy użyte

w eksperymencie. Jeśli oligonukleotyd wystę-

puje kilkakrotnie w sekwencji, w spektrum

pojawia się on tylko raz. Odtworzenie se-

kwencji oryginalnej w przypadku powtórzeń

jest trudniejsze ze względu na wiele możliwo-

ści optymalnego (z obliczeniowego punktu

widzenia) dopasowania fragmentów. W celu

zminimalizowania liczby błędów pochodzą-

cych z eksperymentu zaproponowane zostały

nowe podejścia, które w trakcie hybrydyzacji

zamiast oligonukleotydów o równej długości

wykorzystują oligonukleotydy, dla których

temperatura zajścia idealnej hybrydyzacji jest

taka sama, albo oligonukleotydy zawierające

zasady uniwersalne mogące łączyć się z do-

wolną zasadą azotową. Algorytmy do rekon-

strukcji oryginalnej sekwencji dla każdego z

tych podejść (z oligonukleotydami o równej

długości, o równej temperaturze oraz zawie-

rającymi zasady uniwersalne) zostaną przed-

stawione w części „Sekwencjonowanie”.

Od niedawna stosuje się wysoce zautoma-

tyzowane podejścia do sekwencjonowania,

które w stosunkowo krótkim czasie generują

miliony sekwencji o długości kilkudziesięciu

(technologie Solexa i SOLID) lub kilkuset

nukleotydów (metoda 454). Każdy z tych

systemów posiada inną strategię generowa-

nia wysokiej jakości danych. Metoda 454

firmy Roche (M arGulies i współaut. 2005)

polega na sekwencjonowaniu przez syntezę

fragmentów DNA, które przyczepiane są do

małych koralików, a następnie klonowane w

emulsji wodno-oleistej. Kolejne nukleotydy są

syntetyzowane do fragmentów DNA w wyni-

ku szeregu reakcji chemicznych. W jednym

kroku może się przykleić nawet kilka nukle-

otydów tego samego typu. Podczas sekwen-

cjonowania metodą Solexa firmy Illumina

(B ennett 2004) sekwencja jest przyczepiana

do powierzchni płytki. Do sekwencji kolejno

syntetyzowane są specjalnie zaprojektowane,

znakowane kolorem nukleotydy, które każdo-

razowo kończą syntezę DNA, a następnie (po

odczytaniu koloru) za pomocą enzymów są

odblokowywane. W metodzie SOLID firmy

Applied Biosystems (F u i współaut. 2008),

zamiast pojedynczych nukleotydów, przycze-

piane są znakowane kolorem oligonukleoty-

dy, w których znane są pierwsze dwa nukle-

otydy. W kolejnych cyklach znajdowane są

następne dwu-nukleotydy tak, że w sumie

otrzymujemy podwójne pokrycie każdej po-

zycji w sekwencji.

Ryc. 1. Rekonstrukcja sekwencji AGGACG ze

spektrum S= {ACG, AGG, GAC, GGA}.

Kolejne oligonukleotydy nakładają się na siebie z

przesunięciem równym 1. Pokrycie każdej pozycji w

sekwencji wyznacza liczba oligonukleotydów, które

nakładają się w tym miejscu na sekwencję, np. po-

krycie czwartego nukleotydu w sekwencji wynosi 3.

Znalezienie rozwiązania (sekwencji) jest

problemem łatwym obliczeniowo, o ile w

spektrum nie pojawią się błędy (p evzner

1989). Jednakże w rzeczywistym ekspery-

mencie pojawiają się błędy i mogą być one

dwojakiego rodzaju: błędy negatywne — je-

śli brakuje jakiegoś oligonukleotydu w spek-

trum, pomimo że występuje on w badanej

sekwencji — oraz błędy pozytywne, gdy w

spektrum znajduje się nadmiarowy oligonu-

kleotyd, który nie występuje w sekwencji. W

przypadku, gdy w spektrum występują błędy,

podczas rekonstrukcji badanej sekwencji za-

kłada się, że przesunięcie pomiędzy kolejny-

mi elementami spektrum może być większe

niż 1 (co odpowiada błędom negatywnym),

a niektóre oligonukleotydy nie znajdą się w

ogóle w rozwiązaniu (błędy pozytywne). Pro-

blem ten jest wówczas obliczeniowo trudny

(B łażewicz i K asprzaK 2003). Szczególnym

Sekwencjonowanie i asemblacja dNA — podejścia, modele grafowe, algorytmy

Drugi etap odczytywania sekwencji DNA,

asemblacja, polega na połączeniu fragmen-

tów DNA pochodzących z etapu sekwencjo-

nowania w dłuższe odcinki. Spośród wielu

algorytmów, które rozwiązują problem asem-

blacji, kilka z nich zostanie omówionych

tutaj ze względu na interesujące podejścia

grafowe. Są wśród nich zarówno algorytmy

skonstruowane dla dłuższych fragmentów

wejściowych (klasyczna asemblacja), jak i dla

krótszych, pochodzących z nowych metod

sekwencjonowania. Wszystkie one zostaną

przedstawione w części „Asemblacja”.

SEKWENCJONOWANIE

Pierwszy algorytm do rekonstrukcji ory-

ginalnej sekwencji na podstawie wyników

eksperymentu hybrydyzacyjnego z mikroma-

cierzą DNA został zaproponowany w (B ains

i s Mith 1988). Autorzy założyli brak błędów

w spektrum. Algorytm buduje drzewo, w

którym wierzchołkami są elementy spek-

trum, natomiast łuki łączą wierzchołki, jeśli

l -1 ostatnich liter poprzednika nakłada się

na l -1 pierwsze litery następnika. Rozwiąza-

niem jest sekwencja odpowiadająca ścieżce

w tym drzewie od korzenia do liścia zawiera-

jącej wszystkie elementy spektrum dokładnie

jeden raz (Przykład 1). Jako korzeń drzewa

wybierany jest ten element, który jest począt-

kiem badanej sekwencji. Jeśli natomiast nie

jest on znany, algorytm musi skonstruować

| S | drzew z każdym kolejnym elementem

spektrum jako korzeniem.

l ysov i współaut. (1988) zauważyli, że

problem SBH bez błędów w spektrum moż-

na sprowadzić do znanego problemu poszu-

kiwania ścieżki Hamiltona w pewnym grafie.

Skierowany graf H konstruowany jest w na-

stępujący sposób. Każdy wierzchołek grafu

odpowiada innemu elementowi spektrum.

Łuk ( u , v ) łączy wierzchołek u z wierzchoł-

kiem v jeśli l -1 ostatnich liter etykiety wierz-

chołka u nakłada się na l -1 pierwsze litery

etykiety v . W takim grafie poszukiwana jest

ścieżka przechodząca przez wszystkie wierz-

chołki dokładnie jeden raz (ścieżka Hamilto-

na) (Przykład 1).

ce (Ryc. 2). Tylko dolna ścieżka przechodzi

przez wszystkie elementy spektrum. Odczy-

tując kolejne etykiety można zrekonstruować

sekwencję oryginalną.

GAT

AGG GGA

GAC

ACG CGA GAT

Ryc. 2. Drzewo dla metody zaproponowanej

przez B ainsa i s Mitha (1988).

Metoda l ysova i współaut. (1988) utwo-

rzy na podstawie tego samego spektrum graf,

który zaprezentowany jest na Ryc. 3. W gra-

fie istnieje dokładnie jedna ścieżka Hamilto-

na, która odpowiada badanej sekwencji.

AGG GGA

GAC

GAT

CGA

ACG

Ryc. 3. Graf dla metody zaproponowanej przez

l ysova i współaut. (1988).

Powyższe metody działają tylko dla ide-

alnego, bezbłędnego spektrum, jednakże ich

złożoność obliczeniowa jest wykładnicza

(tzn. liczbę elementarnych operacji algoryt-

mu można wyrazić funkcją wykładniczą, w

której rozmiar instancji rozwiązywanego pro-

blemu jest wykładnikiem potęgi). Algorytm

o wielomianowej złożoności czasowej (tzn.

w którym liczba operacji jest wyrażona funk-

cją wielomianową) dla problemu SBH został

zaproponowany przez p evznera (1989), co

dowodzi, że problem należy do klasy pro-

blemów łatwych obliczeniowo. Algorytm ten

szuka ścieżki przechodzącej przez wszystkie

łuki w grafie skierowanym dokładnie raz (tj.

PRZYKŁAD 1

Załóżmy, że dla badanej sekwencji nu-

kleotydów AGGACGAT eksperyment hybry-

dyzacji przebiegł bez błędów, w wyniku

czego otrzymano spektrum: S = {ACG, AGG,

CGA, GAC, GAT, GGA}. Długość badanej se-

kwencji n = 8, długość oligonukleotydów l

= 3, natomiast | S | = n–l+ 1. Dla uproszczenia

przykładu załóżmy, że znany jest pierwszy

oligonukleotyd (AGG). Metoda Bainsa i Smi-

tha tworzy drzewo dodając kolejne elementy

o ile nie znajdują się już w bieżącej ścież-

M arta K asprzaK , a leKsandra Ś wiercz

ścieżki Eulera). Tym razem elementy spek-

trum odpowiadają łukom, a każdy z łuków

wychodzi z wierzchołka, który jest zaetykie-

towany l -1-literowym prefiksem etykiety łuku

i wchodzi do wierzchołka zaetykietowanego

l -1 - literowym sufiksem (Przykład 2).

Transformacja grafu, w którym poszuki-

wana jest ścieżka Hamiltona w graf, w któ-

rym poszukiwana jest ścieżka Eulera, zmienia

złożoność obliczeniową problemu. Klasa gra-

fów etykietowalnych, dla których taka trans-

formacja jest możliwa, została szeroko omó-

wiona przez B łażewicza i współaut. (1999c).

Grafy budowane na podstawie spektrum w

metodzie l ysova i współaut. (1988), nazwa-

ne grafami DNA, należą do klasy grafów ety-

kietowalnych. Grafy skonstruowane na pod-

stawie metody l ysova i współaut. (1988) są

grafami liniowymi grafów p evznera (1989).

Dla takiej pary grafów poszukiwanie ścieżek

Hamiltona i Eulera jest równoważne (patrz

B łażewicz i współaut. 1999c).

Kolejny algorytm zaprezentowany przez

p evznera (1989) dopuszczał błędy negatyw-

ne w spektrum. Złożoność czasowa algo-

rytmu jest wielomianowa, jednakże nie w

każdym przypadku algorytm potrafi znaleźć

rozwiązanie, więc nie może być traktowany

jako algorytm dokładny. W tej metodzie naj-

pierw wyznaczana jest liczba brakujących oli-

gonukleotydów, równa n-l +1-| S |, a następnie

oligonukleotydy te są znajdowane poprzez

transformację problemu sekwencjonowania

do problemu poszukiwania przepływu w

sieci zbudowanej na podstawie grafu dwu-

dzielnego K m,m . Graf dwudzielny zbudowany

jest z wierzchołków grafu Pevznera, które

posiadają różną liczbę łuków wchodzących

i wychodzących. Jeśli różnica ta jest większa

niż 1 dla pewnego wierzchołka, liczba jego

wystąpień jest odpowiednio zwiększana. Po

lewej stronie grafu dwudzielnego umiesz-

czane są wierzchołki z większą liczbą łuków

wchodzących, a po prawej z większą liczbą

łuków wychodzących. Liczba wierzchołków

z lewej strony równa jest liczbie wierzchoł-

ków z prawej strony i wynosi m . Z każdego

wierzchołka z lewej strony wychodzi łuk do

każdego wierzchołka z prawej strony. Koszt

łuku jest równy najmniejszemu przesunięciu

względem siebie etykiet wierzchołków minus

1. Zatem koszt jest równy 1, jeśli nałożenie

etykiet jest równe l- 2, natomiast jeśli etykie-

ty w ogóle się na siebie nie nakładają koszt

równy jest l -1 (koszt łuku oznacza ile wierz-

chołków/oligonukleotydów brakowałoby w

grafie Pevznera, gdybyśmy chcieli dany łuk

wykorzystać). Dodatkowo do grafu dwudziel-

nego dodawane jest źródło s , z którego wy-

chodzą łuki do każdego wierzchołka z lewej

strony, oraz ujście t , do którego dochodzą

łuki od każdego wierzchołka z prawej strony.

W tak skonstruowanej sieci, gdzie wszystkie

łuki mają pojemność równą 1, poszukiwany

jest przepływ o wartości m- 1 i o minimal-

nym koszcie. Jeśli koszt okazałby się równy

liczbie brakujących oligonukleotydów (czyli

n-l +1-| S |), wówczas graf Pevznera zostałby

uzupełniony o brakujące łuki (lub ścieżki) i

można by w nim poszukiwać ścieżki Eulera

(Przykład 2).

PRZYKŁAD 2

Rozważmy tę samą sekwencję jak w Przy-

kładzie 1. Spektrum idealne S = {ACG, AGG,

CGA, GAC, GAT, GGA}. Graf Pevznera został

przedstawiony na Ryc. 4.

GGA

GAC

AGG

GAT

ACG

CGA

Ryc. 4. Graf dla metody p evznera (1989).

Ścieżka Eulera w tym grafie odpowiada

badanej sekwencji AGGACGAT. Załóżmy te-

raz, że podczas eksperymentu hybrydyzacji

wystąpiły błędy i w spektrum brakuje ele-

mentu CGA (błąd negatywny). Aby utworzyć

kompletny graf Pevznera należy wyznaczyć

liczbę brakujących oligonukleotydów n-l+ 1 -

|S| = 1 (przepływu o takim koszcie będzie-

my szukać) a następnie skonstruować sieć z

wierzchołków o różnej liczbie łuków wcho-

dzących i wychodzących w grafie z Ryc. 4

pozbawionego łuku CGA (Ryc. 5).

Ryc. 5. Sieć dla metody p evzner a (1989).

Koszty na łukach są równe przesunięciom

pomiędzy etykietami wierzchołków przy zało-

Sekwencjonowanie i asemblacja dNA — podejścia, modele grafowe, algorytmy

żeniu dokładnego nałożenia prefiksu lewego

wierzchołka z sufiksem prawego wierzchołka.

W takiej sieci szukamy przepływu o wartości

m -1=1 i o koszcie równym 1. Istnieje tylko

jeden taki przepływ zawierający wierzchołki

(CG, GA), co odpowiada brakującemu oligo-

nukleotydowi CGA. Po dodaniu tego łuku do

grafu Pevznera możemy poszukiwać ścieżki

Eulera.

Pierwszy algorytm dokładny, który do-

puszcza występowanie błędów zarówno po-

zytywnych jak i negatywnych w spektrum i

nie wymaga żadnej dodatkowej informacji

o elementach spektrum został zapropono-

wany przez B łażewicza i współaut. (1999a).

Problem sekwencjonowania został sformuło-

wany jako wariant problemu komiwojażera

(ang. selective traveling salesman problem).

Tworzony jest graf skierowany pełny, w któ-

rym wierzchołki są zaetykietowane oligo-

nukleotydami a koszt łuku pomiędzy dwo-

ma wierzchołkami jest równy minimalnemu

przesunięciu odpowiednich etykiet przy za-

łożeniu ich dokładnego nałożenia. Z każdym

wierzchołkiem skojarzony jest zysk o warto-

ści 1. W takim grafie poszukiwana jest ścież-

ka o największym sumarycznym zysku i kosz-

cie nie przekraczającym n-l , która jest równo-

ważna sekwencji o długości nie większej niż

n utworzonej z maksymalnej liczby oligonu-

kleotydów ze spektrum (Przykład 3).

nie większym niż n-l =5. Jako rezultat otrzy-

mujemy dwie ścieżki, które przechodzą przez

5 wierzchołków i o koszcie równym 5. Od-

czytując etykiety wierzchołków otrzymamy

rozwiązania: AGGACGAT oraz AGGACACG.

Problem SBH w przypadku, gdy wystę-

pują błędy, nawet błędy tylko jednego typu

(pozytywne lub negatywne), jest problemem

silnie NP-trudnym (tzn. prawdopodobnie nie

może dla niego zostać skonstruowany algo-

rytm o wielomianowej złożoności czasowej).

Zatem algorytm dokładny dla większych in-

stancji problemu może nie zakończyć dzia-

łania w sensownym czasie. W efekcie więk-

szość zaproponowanych dla tego problemu

algorytmów to heurystyki, które w znacznej

mierze przyspieszają przeszukanie przestrze-

ni rozwiązań, jednocześnie nie gwarantując,

że znalezione rozwiązanie jest optymalne.

K ruGlyaK (1998) zaproponował wielo-

stopniowe podejście do SBH. Kolejno wyko-

nywane są eksperymenty z oligonukleotyda-

mi o zwiększanej długości. Pierwszy ekspe-

ryment przeprowadzany jest dla kompletnej

biblioteki oligonukleotydów o długości l . W

kolejnych eksperymentach biblioteka oligo-

nukleotydów składa się z połączonych ze

sobą oligonukleotydów, które znalazły się w

spektrum po poprzednim eksperymencie,

przy założeniu jak największego możliwego

nałożenia. Ustawienie progu dla minimalnej

liczby nakładających się nukleotydów na war-

tość większą niż 0 (0 oznacza konkatenację

pary oligonukleotydów) zmniejsza znacząco

rozmiar biblioteki generowanej dla następne-

go eksperymentu. Zaletą tego podejścia jest

redukcja błędów wynikających z powtórzeń

ciągu nukleotydów w sekwencji oryginalnej,

jednakże inne błędy negatywne oraz błędy

pozytywne będą propagowane w kolejnych

krokach tej metody.

Interaktywne podejście do SBH zapropo-

nowali p han i s Kiena (2001). Zwykły ekspe-

ryment hybrydyzacji wzbogacony jest o do-

datkowe eksperymenty, które mają na celu

rozwianie wszelkich niejednoznaczności

podczas rekonstrukcji oryginalnej sekwen-

cji. Metoda ta konstruuje graf, taki sam jak

w przypadku metody l ysova i współaut.

(1988), gdzie oligonukleotydy są etykietami

wierzchołków a łuki łączą wierzchołki, któ-

rych etykiety przesunięte są względem siebie

o jedną pozycję. Seria dodatkowych zapytań

(eksperymentów biochemicznych) z oligonu-

kleotydami o zwiększanej długości pozwala

na wyeliminowanie niepotrzebnych wierz-

PRZYKŁAD 3

Załóżmy, że dla sekwencji AGGACGAT spek-

trum z błędami negatywnymi i pozytywnymi

wygląda następująco: S = {ACG, AGG, CAC,

GAC, GAT, GGA}. W spektrum, oprócz błędu

negatywnego CGA, pojawił się jeden błąd po-

zytywny CAC. Graf pełny skonstruowany dla tej

metody (B łażewicz i współaut. 1999a) zawierał-

by łuki o kosztach przedstawionych Tabela 1.

Zysk za odwiedzenie każdego wierzchoł-

ka jest równy 1. Algorytm szuka ścieżki prze-

chodzącej co najwyżej jeden raz przez każdy

wierzchołek, o maksymalnym zysku i koszcie,

Tabela 1. Koszty łuków dla metody wg B łaże -

wicza i współaut. (1999a).

ACG AGG CAC GAC GAT GGA

ACG – 3 3 2 2 2

AGG 3 – 3 2 2 1

CAC 1 3 – 3 3 3

GAC 1 3 2 – 3 3

GAT 3 3 3 3 – 3

GGA 2 2 3 1 1 –

Plik z chomika:

Inne pliki z tego folderu:

Inne foldery tego chomika: