Algorytm Hannenhallego




Definicje i rysunki

gen (ang. gene) - biologicznie to najmniejsza funkcjonalna część genomu, u nas to liczba całkowita ze znakiem. Geny oznacza się małymi literami: a, b, x, y.


chromosom (ang. chromosome) - ciąg genów, czyli ciąg liczb całkowitych. Utożsamiamy chromosom z jego odwróconą wersją. Chromosomy oznaczam przez X, Y i Z, np. X = (x1,x2, ...,xn), gdzie xi to geny


segment (ang. segment) - część chromosomu, czyli ciąg genów. Oznaczenie X = (X1,X2) oznacza, że chromosom X jest złożony z segmentów X1 i X2. Jeżeli segment S = (s1,s2,...,sn) to -S = (-sn,-sn-1,...,s1)


genom (ang. genome) - zbiór chromosomów, kolejność chromosomów w genomie nie ma znaczenia. Dwa genomy uznajemy za identyczne, gdy zawierają takie same chromosomy. Genomy oznaczam przez A i B. Zwykle n to liczba genów w genomie, a N to liczba chromosomów.


translokacja (ang. translocation) - operacja tworząca z dwóch chromosomów dwa nowe. Rozbija każdy chromosom na dwa niepuste segmenty i łączy je na nowo biorąc po jednym segmencie z każdego chromosomu. Okazuje się, że można tego dokonać na dwa sposoby (patrz animacja i rysunek 1). Niech X = (X1, X2) i Y = (Y1, Y2) będą dowolnymi chromosomami. Wynikiem translokacji prefix-prefix jest para (X1, Y2) i (Y1, X2), czyli następuje wymiana sufiksów. Wynikiem translokacji prefix-sufix jest para (X1, -Y1) i (-X2, Y2). Warto zauważyć, że translokacja prefix-prefix to to samo co translokacja prefix-sufix na parze z jednym odwróconym chromosomem (czyli parze X i -Y lub -X i Y). Stąd wniosek, że nie należy przywiązywać wagi do nazw translokacji, ale mieć świadomość iż dla zadanego podziału chromosomów na segmenty można otrzymać dokładnie dwie różne pary wynikowe. Oznaczenie ropp(X,Y,i,j) lub rops(X,Y,i,j) w zależności od tego czy to translokacja prefix-prefix czy prefix-sufix. Indeksy i oraz j oznaczają, że podział następuje w chromosomie X między i-1 i i-tym genem, a w Y między j-1 a j-tym.


Rysunek 1. Dwa rodzaje translokacji: prefix-prefix i prefix-sufix
[rysunek 1]
Zobacz rysunek

Zobacz animacje


odległość translokacyjna (ang. translocation distance) - minimalna liczba translokacji przekształcających jeden genom w drugi. Oznaczana przez d(A, B) gdzie A i B są genomami.


ewolucja genomu (ang. genome evolution) - ciąg zmian przekształcających jeden genom w drugi, w tych rozważaniach jedyną dopuszczalną zmianą jest translokacja.


geny końcowe (ang. nodal genes) - dwa skrajne geny w każdym chromosomie, czyli pierwszy i ostatni. Geny końcowe na animacji są zaznaczone innym kolorem.


graf cykliczny genomu (ang. cycle graph) - (patrz rys. 2 i 4). Każdemu genowi odpowiada para wierzchołków (xit,xih) dla xi > 0 i (xih,xit) dla -xi. Wierzczołki grafu zwykle oznacza się przez przez u, v, f, g. Czarne krawędzie łączą sąsiadów w A, a szare (na rysunkach zaznaczone linią przerywaną) sąsiadów w genomie docelowym, czyli wierzchołki z kolejnymi numerami: 1h z 2t, 2h z 3t, itd. Nie rysuje się krawędzi między xh i xt. Oznaczenie: GA.


Rysunek 2. Ewolucja (sortowanie) genomu i zmiany grafu cyklicznego, które towarzyszą translokacjom.
[rysunek 2]
Zobacz rysunek



cykl (ang. cycle) - cykl w grafie to zamknięty ciąg krawędzi, każde dwie sąsiednie krawędzie w ciągu mają wspólny wierzchołek. W grafie cyklicznym krawędzie w cyklach są na przemian czarne i szare. Ilość cykli w GA to cA. Długość cyklu to liczba krawędzi wchodzących w jego skład, w GA cykle są parzystej długości.


długi cykl - cykl o długości większej niż 2.


właściwa translokacja (ang. proper translocation) - translokacja zwiększająca liczbę cykli w GA, dla tej translokacji delta(cA) = 1, patrz rys. 3.


niewłaściwa translokacja (ang. improper translocation) - translokacja nie zmieniająca liczby cykli w GA, delta(cA) = 0, patrz rys. 3.


zła translokacja (ang. bad translocation) - translokacja zmniejszająca liczbę cykli w GA, delta(cA) = -1, patrz rys. 3.


Rysunek 3. Wpływ translokacji na ilość cykli w grafie.
[rysunek 3]
Zobacz rysunek



podpermutacja (ang. subpermutation, w skrócie SP) - część chromosomu, z której nie wychodzą na zewnątrz żadne szare krawędzie (patrz rys. 4). Inna, równoważna definicja: nieidentycznościowa permutacja segmentu chromosomu zawierającego kolejne geny z chromosomu decelowego, o dodatkowej własności, że pozostawia na swoim miejscu skrajne (tj. pierwszy i ostatni) elementy, np. (10 11 -13 -14 -12 15) lub (3 5 4 6) jeżeli w genomie docelowym istnieją segmenty (10 11 12 13 14 15) lub (3 4 5 6) na jednym z chromosomów.


minimalna podpermutacja (ang. minimal subpermutation, w skrócie minSP) - taka podpermutacja, która nie zawiera w sobie żadnych innych podpermutacji (jest elementem minimalnym w częściowym porządku podpermutacji wyznaczonym przez zawieranie), np. (1 3 5 7 6 8 2 4 9) - SP, ale nie minSP, a (5 7 6 8) - minSP liczbę minSP w GA to oznaczamy sA. oA to zmienna określająca parzystość liczby minSP, jeżeli liczba jest nieparzysta to oA = 1, gdy parzysta to oA = 0.


Rysunek 4. Popermutacje (SP) i minimalne podpermutacjie (minSP).
[rysunek 4]
Zobacz rysunek



parzysta izolacja (ang. even-isolation) - Własność genomu A:
(i) wszystkie minSP są na jednym chromosomie
(ii) sA jest parzyste
(iii) wszystkie minSP są na jednym SP
. Własność tę oznaczamy przez iA. iA = 1 gdy genom ma własność parzystej izolacji, iA = 0 wpp.


poprawna translokacja (ang. valid translocation) - translokacja dla której d(cA - sA - oA - 2*iA) = 1.



1. Wstęp i definicje

Opiszę pracę Sridhara Hannenhallego (z 1996 roku, patrz bibliografia), w której jest prezentowany wielomianowy algorytm wyznaczania odległości translokacyjnej. Algorytm jest konstruktywny, tzn. oprócz wyznaczenia odległości d(A, B) podaje także minimalny ciąg translokacji prowadzących z A do B. Hannenhalli, na potrzeby tej pracy, przyjmuje kilka założeń:

  • każdy gen występuje w genomie dokładnie jeden raz,
  • geny mają znak, czyli dane są skierowane,
  • są dozwolone tylko właściwe translokcje, czyli takie gdzie z dwóch chromosomów powstają dwa nowe - fuzje i fizje są zabronione,
  • algorytm pomija istnienie centromerów (ang. centromere). Centromer występuje w każdym prawdziwym chromosomie. W biologii dozwolone są tylko takie translokacje, które zostawiają centromery w obu wynikowych chromosomach, tutaj pomijamy to zagadnienie.

Geny końcowe to geny wysępujące na początku i na końcu każdego chromosomu. Warto zauważyć, że translokacje nie zmieniają zbioru genów końcowych - dobrze to widać na animacji, gdzie geny końcowe są zaznaczone innym kolorem. Stąd wniosek, że dany genom da się przekształcić translokacjami tylko w taki, który ma ten sam zbiór genów końcowych. Odległość translokacyjna d(A, B) jest określona tylko dla genomów wzajemnie osiągalnych, czyli tych z identycznym zbiorem genów końcowych.

Na potrzeby kombinatorycznej teorii ogległości translokacyjnej geny oznaczamy liczbami całkowitymi - różnym genom odpowiadają różne liczby, chromosomy to ciągi liczb, a genomy to zbiory chromosomów. Liczby a i -a odpowiadają temu samemu genowi, stąd jeżeli w genomie A występuje a to nie może wystąpić -a. Zatem genomy możemy przedstawiać jako: A = ((a11,a12,...,a1m1), (a21,a22,...,a2m2),..., (aN1,aN2,...,aNmN)) i B = ((b11,b12,...,b1,m1), (b21,b22,...,b2,m2),..., (bN1,bN2,...,bNmN)). Każdy gen występuje dokładnie raz zarówno w A jak i w B, ale mogą pojawiać się z różnymi znakami. Niech X = (x1,x2,...,xm) i Y = (y1,y2,...,yn) będą dwoma chromosomami. Translokację działającą na X i Y oznaczami przez ro(X,Y,i,j), 1 < i <= m, 1 < j <= n jeżeli podział na segmenty następuje w X między xi-1 a xi, a w Y między yi-1 i yi. Translokację prefix-prefix oznaczamy przez ropp, a prefix-sufix przez rops. Mamy ropp(X,Y,i,j) = ((x1,...,xi-1,yj,...,yn), (y1,..yj-1,xi,...,xm)) i rops(X,Y,i,j) = ((x1,...,xi-1,-yj-1,...-y1), (-xm,...,-xi,yj,...yn)). Korzystamy tu z tego, że odwrócenie segmentu odwraca zarówno kolejność genów jak i ich znak. Wynik translokacji ro na genomie A to A*ro.

Dla uproszczenia notacji, zamiast mówić o odległości d(A,B) między genomami A i B mówi się od d(A) przyjmując, że docelowy genom B jest ustalony. Można sobie wyobrazić, że B jest posortowany a sprowadzenie A do B to po prostu sortowanie A. Żeby B był posortowany to wystarczy geny ponumerować w kolejności ich występowania w B. Genomów docelowych może być wiele, np. ((1,2,3,4),(5,6)) lub ((1,2,3),(4,5,6)) - żadnego z nich nie da się przeprowadzić do drugiego, ponieważ mają inne geny końcowe.

Żeby zobrazować kombinatoryczną strukturę genomu wprowadza się graf cykliczny. Każdemu genowi +xi > 0 odpowiada w grafie skierowana para wierzchołków (xit,xih). Zachowuję tu oryginalne oznaczenie Hannenhallego, gdzie t oznacza tail (ogon) a h head (głowa) a genowi -xi odwrócona para: (xih,xit). W grafie są krawędzie czarne i szare. Krawędzie czarne łączą sąsiadów w isteniejącym genomie, a krawędzie szare łączą sąsiadów w genomie docelowym. Sąsiedzi to wierzchołki grafu odpowiadające genom, które są obok siebie w jakimś chromosomie. Nie rysuje się krawędzi między dwoma wierchołkami tego samego genu (tzn. między xih i xit). Zapisujemy graf cykliczny GA == GAB(V,E), V = {u: u = xt albo u = xh, gdzie x to gen z A}. Szczegóły na rysunku.

Kilka właściwości grafu cyklicznego:

  • liczba czarnych krawędzi (szarych też) w GA to n - N, gdzie n to liczba genów a N to liczba chromosomów w genomie A
  • z każdego wierzchołka wychodzi dokładnie jedna krawędź czarna i dokładnie jedna szara (dla skrajnych żadna)
  • druga własność gwarantuje jednoznaczne rozbicie grafu na cykle

Lemat 1.
Liczba cykli w GA osiąga maksimum (n - N) wtedy i tylko wtedy gdy A jest identyczny z genomem docelowym.
Dowód W grafie docelowym wszystkie cykle mają minimalną długość, czyli 2, a liczba krawędzi w całym grafie się nie zmienia.

Dla uproszczenia oznaczamy wszystkie pośrednie genomy w ewolucji z A do B jako A. Każdej zmianie w A (czyli translokacja) odpowiada zmiana w GA.

Niech ro == ro(X,Y,i,j) będzie translokacją na X = (x1,x2,...,xm) i Y = (y1,y2,...,yn). Niech f i g będą wierzchołkami-sąsiadami w GA, takie że f nalezy do {xi-1t,xi-1h} a g nalezy do {xit,xih}. Podobnie u i v są sąsiadami w Y, u nalezy do {yj-1t,yj-1h}, v nalezy do {yjt,yjh}. Warunki te oznaczają, że x oraz u są najbardziej prawymi wierzchołkami w lewych segmentach chromosomów odpowiednio X i Y, a y i v najbardziej lewymi w prawych segmentach (patrz rys. 3). Powiemy, że ro działa (ang. cuts) na czarnych krawędziach (f g) i (u v). Translokacje dzielimy na trzy kategorie: właściwe (ang. proper), niewłaściwe (ang. improper) i złe (ang. bad) w zależności od tego jak wplywają na liczbę cykli w GA.

Translokacja prefix-prefix działająca na czarne krawędzie (f g) i (u v) jest właściwa (ang. proper), gdy w GA istnieje cykl (u v ... f g ... u). Translokacja prefix-sufix jest właściwa gdy w GA jest cykl (u v ... g f ... u). Te z pozoru dziwnie wyglądające właściwości zapewnają to, że po wykonaniu właściwej translokacji jakiś cykl zostanie rozbity na dwa (patrz rys. 3). Translokacja jest niewłaściwa jeżeli czarne krawędzie na których działa należą do tego samego cyklu, ale nie jest właściwa. Nie zmienia ona ilości cykli. Pozostałe translacje nazywamy złymi - działają na czarnych krawędziach z dwu różnych cykli, więc zmniejszają liczbę cykli. Ogólnie można powiedzieć, że lepsze są translokacje rozbijające cykle na mniejsze gdyż w genomie docelowym ilość cykli jest maksymalna.

2. Dolne ograniczenie na odleglość translokacyjną

W tym punkcie zajmiemy się znalezieniem dolnego ograniczenia na odległość. Podając dolne ograniczenie określamy ile co najmniej potrzeba translokacji na przeprowadzenie genomu A w B.

Przyjmiemy oznaczenie, że delta(psi) to zmiana jakiejś wartości psi po dokonaniu translokacji. Tak więc delta(psi) = psiA*ro - psiA. Wielkością psi może być np. ilość cykli w grafie GA.

Lemat 2.
Dla dowolnej translokacji ro:
  • delta(cA) = 1, gdy ro jest właściwa
  • delta(cA) = 0, gdy ro jest niewłaściwa
  • delta(cA) = -1, gdy ro jest zła.
Dla dalszych rozważań ważna jest nierówność delta(cA) <= 1

Wiadomo, że genom docelowy ma maksymalną liczbę cykli równą n - N. Jeżeli mamy genom A o cA cyklach to liczba (n - N) - cA określa "niedobór" cykli w A, tzn. tyle cykli brakuje, żeby było ich tyle co w genomie docelowym. Lemat 2 mówi, że jedna translokacja może dodać co najwyżej jeden cykl, stąd Twierdzenie 3.

Twierdzenie 3.
Dla dowolnego genomu A, d(A) >= n - N - cA.

Mamy więc pierwsze dolne ograniczenie na d(A), wynikłe z istnienia długich (o długości > 2) cykli w A. Nasuwa się pytanie czy zachodzi równość d(A) = n - N - cA. Taka równość oznaczałaby, że w każdym kroku sortowania genomu można wykonać translokację właściwą, dodającą nowy cykl. Okazuje się, że nie zawsze jest to możliwe. W dalszej części pokażę własności genomów, a zatem grafu GA, które uniemożliwają takie sortowanie. Te nowe własności wprowadzą do dolnego ograniczenia nowe parametry, pierwszym parametrem jest n - N - cA, czyli "niedobór cykli".

Okazuje się, że trudna do posortowania jest taka część chromosomu, która zawiera geny mające znaleźć się obok siebie w genomie wynikowym, ale wymieszane w innej kolejności. To prowadzi do pojęcia podpermutacji. Jednak żeby zdefiniować podpermutację potrzeba kilku nowych oznaczeń. Niech I = (xi,xi+1,...xj) będzie segmentem w A. Niech VI będzie zbiorem wierzchołków związanych z I, tj. VI = {u: u = xkt lub u = xkh, i <= k <= j}. Przez LEFT(I) rozumiemy skrajnie lewy wierzchołek z I (czyli xih lub xit w zależności od tego czy xi jest dodatnie czy ujemne). Na rysunku 4 mamy, dla I = (2,4,3,5), LEFT(I) = 2t, RIGHT(I) = 5h. Niech IN(I) = VI - LEFT(I) - RIGHT(I), dla powyższego przykładu IN(2,4,3,5) to {2h,4t,4h,3t,3h,5t}. O krawędzi (u v) powiemy, że jest wewnątrz przedziału I jeżeli oba wierzchołki należą do IN(I). Wewnątrz (2,4,3,5) są trzy czarne i trzy szare krawędzie. Istnieją dwie równoważne definicje podpermutacji (SP):

  1. podpermutacja to taki przedział (xi,xi+1,...,xj) w chromosomie X, że istnieje w Y przedział (xi,permutacja(xi+1,...,xj-1),xj), gdzie permutacja(xi+1,...xj-1) <> (xi+1,...,xj-1) i oba z tych ciągów zawierają te same geny,
  2. podpermutacja to przedział I spełniający dwa warunki:
    (i) nie istnieje krawędź (u v), t. że u nalezy do IN(I) i v nie nalezy do IN(I)
    (ii) jest przynajmniej jeden długi (długość > 2) cykl w IN(I)

Oprócz tego definiuje się minimalną podpermutację (minSP), czyli taką podpermutację, która nie zawiera w sobie innych podpermutacji. Rozmiar podpermutacji to ilość genów w niej zawartych. Na rysunku 4 podpermutacja (1,2,4,3,5,6) ma rozmiar 6 a minimalna podpermutacja (2,4,3,5) ma rozmiar 4.

Wróćmy teraz do dolnego ograniczenia na odlgełość. Napisałem, że nie zawsze można znaleźć translokację zwiększającą liczbę cykli. Przeszkodą są właśnie podpermutacje (SP). Żeby zniszczyć podpermutację trzeba wykonać translokację, która działa na jednej z czarnych krawędzi należących do podpermutacji. Zniszczenie podpermutacji oznacza takie przekształcenie genomu po którym geny stanowiące podpermutacją przestaną ją tworzyć w genomie wynikowym. Każda translokacja działa na dwóch krawędziach z różnych chromosomów, więc translokacja użyta do zniszczenia podpermutacji połączy dwa cykle w jeden - jest złą translokacją. Nie da się tego uniknąć, gdyż każdy z cykli podpermutacji jest zawarty wewnątrz jednego chromosomu. Można jednak wziąć dwie czarne krawędzie należące do różnych podpermutacji w różnych chromosomach. Do szacowania odległości używa się minimalnych podpermutacji (minSP), gdyż minSP mają taką własność, że zniszczenie minSP niszczy wszystkie podpermutacje je zawierające. sA oznacza liczbę minSP w GA.

Lemat 4.
Dla każdej translokacji delta(cA - sA) <= 1

Dowód. Lemat mówi, że nie da się naraz zwiększać liczby cykli i zmniejszać liczby minSP. Wynika to bezpośrednio z tego, że translokacje niszczące minSP sklejają cykle.

W genomie docelowym nie ma minSP (sA = 0), więc nowa dolna granica to d(A) >= n - N - cA + sA. Mamy już drugi parametr utrudniający sortowanie: sA. I znowu pojawia się pytanie: czy zawsze można wykonać translokację taką, że delta(cA - sA) = 1. Jeżeli tak by było i sA > 0, to za każdym razem trzeba by niszczyć dwa minSP (zniszczenie jednego to delta(cA - sA) = 0). Nie zawsze się to da zrobić - liczba minSP może być nieparzysta. Niech oA = 1, gdy liczba minSP jest nieparzysta, oA = 0 gdy parzysta.

Lemat 5.
Dla każdej translokacji delta(cA - sA - oA) <= 1

Dowód. Jeżeli uda nam się zlikwidować własność nieparzystości (delta(-oA) > 0) to wtedy zniszczymy tylko jedno minSP, więc delta(cA - sA) = 0. W drugą stronę, gdy delta(cA - sA) = 1 (niszczę albo 2 minSP albo żadnego i rozbijam dwa cykle) to delta(-oA) = 0.

Na koniec pozostaje jeszcze jedna własność "utrudniająca": parzysta izolacja. Oto jej definicja. Genom A na własność parzystej izolacji (iA = 1), gdy:
(i) wszystkie minSP są na jednym chromosomie
(ii) sA jest parzyste
(iii) wszystkie minSP są na jednym SP
Dlaczego ta dziwna własność jest ważna? Bo gdy genom ma tę własność to, żeby go postortować trzeba wykonać translokację dla której delta(cA - sA - oA) <= -1. Genom docelowy nie ma własności parzystej izolacji (iA = 0), więc żeby posortować A trzeba w krórymś kroku tę własność zniszczyć. Rozważmy translokację, która ją niszczy - udowodnię, że zachodzi dla niej delta(cA - sA - oA) <= -1. Można próbować niszczyć parzystą izolację na kilka sposobów:

  • utowrzyć nowe minSP. Gdy delta(sA) = 1 to delta(oA) = 1, a gdy delta(sA) = 2 to delta(sA) = 0 - ogólnie delta(sA + oA) = 2. Dodając do tego to, że zawsze delta (sA) <= 1 mamy delta(cA - sA - oA) <= delta(sA) - delta(sA + oA) <= -1
  • spróbować rozrzucić minSP z jednego na dwa chromosomy, ale ponieważ są wewnątrz jednej podpermutacji (własność (iii)) a z podpermutacji nie wychodzą na zewnątrz szare krawędzie to taka translokacja będzie zła Wtedy delta(sA) = delta(oA) = 0, delta(cA) = -1, delta(cA - sA - oA) = -1
  • zniszczyć jedno minSP (dwóch się nie da): delta(sA) = delta(cA) = -1, delta(oA) = 1 i w końcu delta(cA - sA - oA) = -1.

To prowadzi do kolejnego ograniczenia na pojedynczą translokację:

Lemat 6.
Dla dowolnej translokacji zachodzi delta(cA - sA - oA - 2*iA) <= 1, iA = 1 dla patrzystej izolacji, 0 wpp.

Dowód. Wynika z poprzednich rozważań o parzystej izolacji.

Gdy dodamy do tego to, że w genomie docelowym oA = iA = 0, to otrzymamy kolejne, lepsze bo większe ograniczenie dolne:

Twierdzenie 7.
Dla dowolnego genomu A, d(A) >= n - N - cA + sA + oA + 2*iA.

Wzór należy rozumieć w ten sposób: potrzebna jest jedna translokacja na rozbicie każdego długiego cyklu (n - N - cA), dodatkowo jedna translokacja na zniszczenie minSP (sA), jedna na usunięcie niepatrzystości minSP (oA) i wreszcie dwie na usunięcie parzystej izolacji (2*iA).

3. Ograniczenie górne

Okazuje się, że podane w twierdzeniu 7 ograniczenie dolne może być osiągnięte dla każdego genomu. Żeby to udowodnić trzeba pokazać, że dla każdego genomu można wykonąć poprawną (ang. proper) translokację. Translokacja jest poprawna, gdy delta(cA - sA - oA - 2*iA) = 1. Możemy dla dowolnego genomu zdefiniować funkcję f(A) = n - N - cA + sA + oA + 2*iA. Dla genomu docelowego (i tylko dla niego) f(A) = f(B) = 0, bo cB = n - N, sB = oB = iB = 0. Po wykonaniu poprawnej translokacji na genomie A d(f(A)) = -1, więc po f(A) krokach genom A można sprowadzić do docelowego. Hannenhalli dowodzi istnienia poprawnej translokacji dla każdego genomu wskazując algorytm jej znajdywania. (Szczegóły w oryginalnej Hannenhallego, dość nużący dowód z dużą liczbą lematów i rozważanych przypadków).

Twierdzenie 14.
Dla dowolnego genomu A
  d(A) = n - N - cA + sA + oA + 2*iA.

Algorytm sortujący wygląda następująco:

dopóki genom A jest różny od genomu docelowego
    jeżeli jest właściwa translokacja w A
        zrób właściwą translokację
    w przeciwnym przypadku
        zrób złą, ale poprawną translokację
    koniec jeżeli
    A := A*ro
koniec dopóki

Algorytm Hannenhallego działa w czasie O(n3), gdzie n to liczba genów w A.


Karol Bieńkowski, 2001