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 = ((a₁₁,a₁₂,...,a_1m1), (a₂₁,a₂₂,...,a_2m2),..., (a_N1,a_N2,...,a_NmN)) i B = ((b₁₁,b₁₂,...,b_1,m1), (b₂₁,b₂₂,...,b_2,m2),..., (b_N1,b_N2,...,b_NmN)). 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 = (x₁,x₂,...,x_m) i Y = (y₁,y₂,...,y_n) będą dwoma chromosomami. Translokację działającą na X i Y oznaczami przez (X,Y,i,j), 1 < i <= m, 1 < j <= n jeżeli podział na segmenty następuje w X między x_i-1 a x_i, a w Y między y_i-1 i y_i. Translokację prefix-prefix oznaczamy przez _pp, a prefix-sufix przez _ps. Mamy _pp(X,Y,i,j) = ((x₁,...,x_i-1,y_j,...,y_n), (y₁,..y_j-1,x_i,...,x_m)) i _ps(X,Y,i,j) = ((x₁,...,x_i-1,-y_j-1,...-y₁), (-x_m,...,-x_i,y_j,...y_n)). Korzystamy tu z tego, że odwrócenie segmentu odwraca zarówno kolejność genów jak i ich znak. Wynik translokacji na genomie A to A*.

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 +x_i > 0 odpowiada w grafie skierowana para wierzchołków (x_i^t,x_i^h). Zachowuję tu oryginalne oznaczenie Hannenhallego, gdzie t oznacza tail (ogon) a h head (głowa) a genowi -x_i odwrócona para: (x_i^h,x_i^t). 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 x_i^h i x_i^t). Zapisujemy graf cykliczny G_A == G_AB(V,E), V = {u: u = x^t albo u = x^h, gdzie x to gen z A}. Szczegóły na rysunku.

Kilka właściwości grafu cyklicznego:

• liczba czarnych krawędzi (szarych też) w G_A 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 G_A.

Niech == (X,Y,i,j) będzie translokacją na X = (x₁,x₂,...,x_m) i Y = (y₁,y₂,...,y_n). Niech f i g będą wierzchołkami-sąsiadami w G_A, takie że f {x_i-1^t,x_i-1^h} a g {x_i^t,x_i^h}. Podobnie u i v są sąsiadami w Y, u {y_j-1^t,y_j-1h}, v {y_j^t,y_j^h}. 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 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 G_A.

Translokacja prefix-prefix działająca na czarne krawędzie (f g) i (u v) jest właściwa (ang. proper), gdy w G_A istnieje cykl (u v ... f g ... u). Translokacja prefix-sufix jest właściwa gdy w G_A 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ą

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

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

Wiadomo, że genom docelowy ma maksymalną liczbę cykli równą n - N. Jeżeli mamy genom A o c_A cyklach to liczba (n - N) - c_A 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 - c_A. 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 G_A, które uniemożliwają takie sortowanie. Te nowe własności wprowadzą do dolnego ograniczenia nowe parametry, pierwszym parametrem jest n - N - c_A, 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 = (x_i,x_i+1,...x_j) będzie segmentem w A. Niech V_I będzie zbiorem wierzchołków związanych z I, tj. V_I = {u: u = x_k^t lub u = x_k^h, i <= k <= j}. Przez LEFT(I) rozumiemy skrajnie lewy wierzchołek z I (czyli x_i^h lub x_i^t w zależności od tego czy x_i jest dodatnie czy ujemne). Na rysunku 4 mamy, dla I = (2,4,3,5), LEFT(I) = 2^t, RIGHT(I) = 5^h. Niech IN(I) = V_I - LEFT(I) - RIGHT(I), dla powyższego przykładu IN(2,4,3,5) to {2^h,4^t,4^h,3^t,3^h,5^t}. 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ł (x_i,x_i+1,...,x_j) w chromosomie X, że istnieje w Y przedział (x_i,permutacja(x_i+1,...,x_j-1),x_j), gdzie permutacja(x_i+1,...x_j-1) <> (x_i+1,...,x_j-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 IN(I) i v 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. s_A oznacza liczbę minSP w G_A.

Lemat 4. Dla każdej translokacji (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 (s_A = 0), więc nowa dolna granica to d(A) >= n - N - c_A + s_A. Mamy już drugi parametr utrudniający sortowanie: s_A. I znowu pojawia się pytanie: czy zawsze można wykonać translokację taką, że (c_A - s_A) = 1. Jeżeli tak by było i s_A > 0, to za każdym razem trzeba by niszczyć dwa minSP (zniszczenie jednego to (c_A - s_A) = 0). Nie zawsze się to da zrobić - liczba minSP może być nieparzysta. Niech o_A = 1, gdy liczba minSP jest nieparzysta, o_A = 0 gdy parzysta.

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

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

Na koniec pozostaje jeszcze jedna własność "utrudniająca": parzysta izolacja. Oto jej definicja. Genom A na własność parzystej izolacji (i_A = 1), gdy:
(i) wszystkie minSP są na jednym chromosomie
(ii) s_A 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 (c_A - s_A - o_A) <= -1. Genom docelowy nie ma własności parzystej izolacji (i_A = 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 (c_A - s_A - o_A) <= -1. Można próbować niszczyć parzystą izolację na kilka sposobów:

• utowrzyć nowe minSP. Gdy (s_A) = 1 to (o_A) = 1, a gdy (s_A) = 2 to (s_A) = 0 - ogólnie (s_A + o_A) = 2. Dodając do tego to, że zawsze (s_A) <= 1 mamy (c_A - s_A - o_A) <= (s_A) - (s_A + o_A) <= -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 (s_A) = (o_A) = 0, (c_A) = -1, (c_A - s_A - o_A) = -1
• zniszczyć jedno minSP (dwóch się nie da): (s_A) = (c_A) = -1, (o_A) = 1 i w końcu (c_A - s_A - o_A) = -1.

To prowadzi do kolejnego ograniczenia na pojedynczą translokację:

Lemat 6. Dla dowolnej translokacji zachodzi (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 o_A = i_A = 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 - c_A), dodatkowo jedna translokacja na zniszczenie minSP (s_A), jedna na usunięcie niepatrzystości minSP (o_A) i wreszcie dwie na usunięcie parzystej izolacji (2*i_A).

3. Ograniczenie górne

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*
koniec dopóki

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