Rozdział 16

SIECI NEURONOWE W BIOINFORMATYCE

Rafał ADAMCZAK¹

Wstęp

W ostatnich latach nastąpił rozwój technologii, który pozwolił na szybkie sekwencjonowanie genomów wielu organizmów. Dzięki temu powstają bazy danych sekwencji białkowych, których analiza odkrywa ogromną różnorodność białek i zróżnicowaną funkcję jaką pełnią w systemach biologicznych. Można wyróżnić między innymi: funkcje katalityczną, transportową, magazynową, immunologiczną, regulatorową. Niestety zwiększenie liczby znanych sekwencji białkowych nie prowadzi w bezpośredni sposób do lepszego zrozumienia ich roli w organizmie. Do zrozumienia funkcji białka, w jaki sposób i z jakimi innymi białkami może oddziaływać, pomocna jest znajomość struktury trzeciorzędowej, czyli ułożenia w przestrzeni 3D wszystkich atomów występujących w białku.

W chwili obecnej znanych jest około 6 milionów sekwencji białkowych (wielkości bazy NR – Non Redundant, baza dostępna jest na stronie NCBI – National Center for Biotechnology Information) ale w bazie PDB (Protein Data Bank) [] , w której znajdują się białka o znanej strukturze trzeciorzędowej, występuje tylko około 71 tysięcy struktur, co stanowi mniej niż 1% wszystkich znanych sekwencji. Doświadczalne metody wyznaczania struktur białkowych wykorzystują takie techniki jak: krystalografia rentgenowska czy jądrowy rezonans magnetyczny (NMR – nuclear magnetic resonance). Cały proces odkrywania struktury trzeciorzędowej jest bardzo kosztowny i czasochłonny oraz często dla wielu białek jak na razie technicznie niemożliwy. Nie wydaję się aby w najbliższych czasie sytuacja miała się poprawić tzn. aby było możliwe uzyskanie za pomocą tych metod struktur trzeciorzędowych dla wszystkich znanych białek. Co więcej, liczba znanych sekwencji podwaja się mniej więcej co 18 miesięcy, natomiast znanych struktur co 3 lata. Stąd różnica pomiędzy liczbą znanych sekwencji, a odkrytych struktur będzie się powiększała.

Nie jest możliwe uzyskanie doświadczalnymi metodami struktur trzeciorzędowych dla wszystkich białek, dlatego konieczne jest stosowanie i rozwijanie metod in silico. Zgodnie z tzw. dogmatem Anfinsena [] informacja o strukturze natywnej białka zakodowana jest w sekwencji aminokwasów oraz środowisku, w którym białko się znajduje. Teoretycznie stosując więc metody ab initio, w których minimalizuje się energie swobodną zdefiniowaną w oparciu o sekwencje białka i fizyczny model oddziaływań pomiędzy atomami białka i atomami otaczającego środowiska, powinno uzyskać się strukturę trójwymiarową białka. Niestety w praktyce takie podejście nie prowadzi do dobrych rezultatów. Można wyróżnić dwie główne przyczyn tej porażki, pierwsza z nich związana jest z dużą liczba atomów i wysokim poziom złożoności oddziaływań występujących w białku, co powoduje, że do minimalizacji funkcji energetycznej konieczne są ogromne zasoby mocy obliczeniowej, które nie są i w najbliższym czasie nie będą dostępne. Druga natomiast związana jest z tym, że różnica energetyczna pomiędzy postacią zwiniętą i rozwiniętą białka może by niewielka (nawet rzędu 1 kcal/mol), co przy drobnych nieścisłościach w definicji energii może prowadzić do niepoprawnej struktury trzeciorzędowej. Innymi słowy, niedokładności związane z oceną energii oddziaływań uzyskiwane podczas eksperymentów oraz ograniczony czas obliczeniowy powodują, że niemożliwe staje się określanie struktur białkowych tylko przy użyciu podstawowych zasad fizyki. Dlatego też stworzenie metod pozwalających przewidywać struktury białek stało się jednym z najważniejszych wyzwań bioinformatyki.

Tradycyjnie metody przewidywań struktury trzeciorzędowej dzielone są, na podstawie podobieństwa analizowanej sekwencji białkowej do białek o znanej strukturze, na trzy grupy: homologiczne, rozpoznawania zwoju i ab initio.

Pierwszym krokiem jest analiza sekwencji, która ma na celu identyfikację sekwencji podobnych do zadanej. Identyfikacja sekwencji składa się z dwóch kroków: identyfikacji domen oraz porównania sekwencji z bazami motywów i profili sekwencyjnych charakterystycznych dla znanych domen. Porównanie sekwencji polega na jej nałożeniu (w przypadku białka wielodomenowego – poszczególnych fragmentów sekwencji reprezentujących domeny) na sekwencje znajdujące się w bazie danych zawierającej białka o znanych strukturach PDB. Jeżeli uda się znaleźć białko homologiczne posiadające znaną strukturę wówczas struktura takiego białka staje się szablonem wykorzystywanym następnie do modelowania. Może się zdarzyć, że istnieje więcej niż jeden szablon, w takim przypadku konieczne jest dokonanie superpozycji struktur, która pozwoli wykryć obszary zachowane strukturalnie. Zazwyczaj obszary te związane są ze strukturami drugorzędowymi lub też obszarami istotnymi dla zachowania funkcji białka.

W przypadku, gdy nie daje się znaleźć w bazie PDB homologa poszukuje się tzw. odległego ewolucyjnie homologa czy też typu zwoju (architektury) białka. Występuje wiele technik poszukiwania zwoju, można je podzielić na dwie grupy: metody przewlekania, czyli nakładania sekwencji na strukturę, która może być reprezentowana np. przez tzw. zmienne środowiskowe charakteryzujące poszczególne aminokwasy [] oraz metody bazujące tylko na podobieństwie sekwencji z uwzględnieniem informacji ewolucyjnej poprzez wykorzystanie tzw. profili []. Profil generowany jest przez nałożenie danej sekwencji na inne podobne sekwencje znajdujące się w bazie danych sekwencji nieredundantnych np. NR. Profil może być wygenerowany zarówno dla szablonu czyli dla sekwencji ze znaną strukturą trzeciorzędową, wówczas można dokonać nałożenia sekwencji na profil, jak i dla targetu (czyli sekwencji dla której poszukiwana jest struktura trzeciorzędowa) i wówczas dokonuje się nałożenie profilu na profil. Zastosowanie profili znacznie poprawia czułość w wyszukiwaniu odległych homologów.

Po znalezieniu zwoju konieczne jest optymalne nałożenie sekwencji na strukturę, gdyż dokładność nałożenia otrzymanego przez opisane metody może nie być wystarczająca. Do poprawy otrzymanego nałożenia można wykorzystać porównanie przewidywanej struktury drugorzędowej czy też przewidywanych wartości dostępności aminokwasów dla cząsteczek solwenta z wygenerowanymi wartościami otrzymanymi z rozpatrywanej struktury. Obszary sekwencji, które nie zostały nałożone na szablon są modelowane w postaci pętli i są jednocześnie zazwyczaj największym źródłem błędu w modelu.

Gdy nie daje się znaleźć struktury homologicznej ani odpowiedniego zwoju stosuje się metodę ab initio, w której poszukiwanie struktury trzeciorzędowej opiera się na minimalizowaniu energii swobodnej układu. Ze względu na wielkość układu (tysiące atomów) proces minimalizacji jest długotrwały i może utykać w licznych minimach lokalnych. Dlatego też nawet w przypadku metod ab initio wykorzystuje się metody pośrednie takie jak np. przewidywania struktur drugorzędowych czy kontaktów, które pozwalają narzucić więzy na układ i w ten sposób ograniczyć przeszukiwaną przestrzeń konformacyjną.

Niektóre etapy przewidywania struktury trzeciorzędowej nie dają się w łatwo przedstawić w postaci algorytmu. W związku z czym stosuje się różne techniki uczenia maszynowego, które wykorzystywane są najczęściej na etapie przewidywań struktur drugorzędowych [] [] [], dostępności aminokwasów dla cząsteczek rozpuszczalnika [] [], oszacowania jakości modelu [] [] [], rozpoznawania typu zwoju [].

W tym rozdziale przedstawiono rozwój metod przewidywania struktur drugorzędowych i przewidywania relatywnej solwatacji aminokwasów, począwszy od najprostszych korzystających z prostych statystyk po metody bardziej złożone bazujące na zaawansowanych metodach uczenia maszynowego (najczęściej są to sieci neuronowe). Przedstawiono również dwa serwery SABLE (Solvent AccessiBiLitiEs, dostępny pod adresem sable.cchmc.org) i MINNOU (Membrane protein IdeNtificatioN withOUt explicit use of hydropathy profiles and alignments, minnou.cchmc.org) przeznaczone do przewidywania solwatacji aminokwasów oraz struktur drugorzędowych, w powstaniu których autor miał swój znaczący wkład.

Struktury drugorzędowe

Struktura drugorzędowa wyprowadzana jest ze struktury trzeciorzędowej białka na podstawie różnych kryteriów takich jak: wiązania wodorowe, uproszczone kryteria odległości dla separacji donora i akceptora, wzorce wiązań wodorowych w połączeniu z analizą kątów torsyjnych jak również na podstawie obserwacji tzw. kryteria wizualne []. Z powodu różnych definicji struktur drugorzędowych stosowanych przez programy możliwe jest, po zastosowaniu tych programów, uzyskanie różnych przypisań dla danego aminokwasu. Największe różnice występują dla obszarów brzegowych beta kartek oraz alfa helis. W chwili obecnej do przygotowania zbioru treningowego dla metod przewidywania struktur drugorzędowych najczęściej korzysta się z programu DSSP [].

DSSP

Struktury drugorzędowe są formalnie zdefiniowane za pomocą wiązań wodorowych obserwowanych pomiędzy grupą aminową i karboksylową dwóch aminokwasów. Niestety nie istnieje jedna poprawna definicja wiązania wodorowego, często definicje wiązania są dostosowywane do określonego celu. W przypadku DSSP wiązanie wodorowe definiowane jest poprzez funkcję energetyczną zdefiniowaną poniżej:

	(1)

gdzie q₁=0.42e, q₂=0.2e, r(AB) odległość pomiędzy atomami A i B, f czynnik stały równy 332. Ładunek umieszczony jest w miejscach gdzie znajdują się atomy C, O (ładunek +q₁, -q₁) oraz N, H (-q₂, +q₂). Wiązanie wodorowe ma energie na poziomie -3 kcal/mol jednak DSSP jest znacznie mniej restrykcyjne i definiuje wiązanie wodorowe dla energii < -0.5 kcal/mol.

Struktury drugorzędowe są definiowane przez DSSP na podstawie występowania pewnych wzorców wiązań wodorowych. Najprostszym wzorcem jest skręt_n, w którym występuje pojedyncze wiązanie wodorowe pomiędzy grupą karboksylową i-tego aminokwasu i grupą aminową i+n-tego aminokwasu, gdzie n=3, 4, 5. Minimalna helisa zdefiniowana jest przez 2 następujące po sobie skręty_n np. alfa helisa₄ występująca pomiędzy aminokwasami i i i+3 wymaga dwóch skrętów₄ dla aminokwasów i-1 i i. DSSP wyróżnia jeszcze: helisy₃ i helisy₅. Poniżej przedstawiona jest definicja dla omówionych helis.

Dłuższe alfa helisy zdefiniowane są jako nałożenie helis minimalnych. Przez DSSP helisy₃ oznaczone są za pomocą litery G, helisy₄ H, helisy₅ I. Helisy₃ lub helisy₅, które zostały zmniejszone poniżej minimalnej długości na skutek nakładania się struktur zostają przez DSSP oznaczone przez T. Jednak w przewidywaniach struktur drugorzędowych zazwyczaj nie dokonuje się rozróżnień helis w związku z czym wszystkie typy helis G, H, I, T oznaczane są przez jeden symbol.

Dwa r aminokwasy i i j tworzą równoległy (pB) lub antyrównoległy (aB) mostek, gdy występują wiązania wodorowe (Hb) pomiędzy odpowiednimi aminokwasami :

Mostki występujące obok siebie tworzą struktury nazywane beta drabiną i beta kartką. Beta drabiną nazywa się jeden lub kilka po sobie następujących beta mostków jednego typu. Zbiór jednej lub więcej beta drabin, których aminokwasy łączą się ze sobą nazywane są beta kartkami. DSSP zaznacza beta drabinę i beta kartkę tym samym symbolem - E, natomiast pojedynczy mostek symbolem B. Wszystkie obszary, którym nie zostały przypisane powyżej wspomniane struktury traktowane są jako obszary nieustrukturyzowane, nazywane pętlami.

Przewidywanie struktur drugorzędowych

Budowa systemu do przewidywania struktur drugorzędowych z sekwencji aminokwasów wymaga zbioru białek ze znanymi strukturami drugorzędowymi. Niestety określenie struktur drugorzędowych nawet przy ustalonych kryteriach definiujących struktury drugorzędowe nie może być jednoznaczne co wynika z tego, że białka w swym naturalnym środowisku ciągle się poruszają i zmienia swoją strukturę. Zmienność struktury drugorzędowej można zaobserwować w strukturach białkowych otrzymywanych metodą NMR, w tej metodzie generowanych jest zazwyczaj od 10 do 20 modeli dla danego białka, jak również poprzez analizę struktur przedstawicieli rodzin białkowych występujących np. w bazie danych Pfam []. W związku z tym nie jest możliwe zbudowanie predyktora, który będzie przewidywał ze 100% dokładnością struktury drugorzędowe wszystkich białek. Szacuje się, że limitem dokładności przewidywań jest 88% []. W chwili obecnej najlepsze metody osiągają wyniki na poziomie 80%. Przedstawione dokładności dotyczą miary Q₃, zdefiniowanej poniżej, uśrednionej po przewidywaniach uzyskanych dla wielu sekwencji białkowych.

Ocena dokładności przewidywań

Najczęściej używaną miarą oceny przewidywań struktur drugorzędowych jest miara Q₃:

	(2)

gdzie H_p, E_p i L_p oznaczają liczbę poprawnie przewidzianych struktur drugorzędowych przypisanych aminokwasom występującym odpowiednio w: helisach, beta kartkach i pętlach, N liczba wszystkich aminokwasów, dla których dokonane zostało przewidywanie. Niestety miara Q₃ nie zawsze poprawnie oddaje jakość przewidywań. Dla przykładu rozpatrzmy mioglobinę 1m6c, której struktura drugorzędowa przedstawiona jest na poniższym rysunku (rysunek wygenerowany przy pomocy serwera polyview.cchmc.org [])

Jak widać na 153 aminokwasy tylko 35 nie należy do alfa helisy, oznacza to, że w przypadku predyktora, który przewidzi występowanie alfa helisy dla każdego aminokwasu otrzymamy Q₃=77% co jest dość dobrym rezultatem. Tymczasem przydatność takiego przewidywania w przewidywaniach struktur trzeciorzędowych jest znikoma. Dlatego też obok miary Q₃ stosuje się również inne miary jak np. miarę SOV (Segment Overlap) [], która wyliczana jest na podstawie segmentów struktur drugorzędowych, a nie poszczególnych aminokwasów. Poniższy rysunek przedstawia struktury drugorzędowe z zaznaczonymi segmentami reprezentującymi alfa helisy i pętle.

Obserwowane CCCCCCCCCCCCCCHHHHHHHHCCCCCCCHHHHHHHHHHHCCCC

Przewidywane CCCCCCCCCCCCCCCCCCCCCCCCCHHHHHHHCCCCCHHHCCCC

Niech s₁ oznacza segment z obszaru obserwowanego a s₂ z przewidywanego, S(i) oznacza zbiór wszystkich par nakładających się segmentów (s₁,s₂) dla i-tej struktury drugorzędowej, S'(i) zbiór wszystkich segmentów s₁, dla których nie ma nakładających się segmentów w s₂. Miarę SOV dla wybranej struktury drugorzędowej definiuje poniższy wzór

	(3)

gdzie minov(s₁,s₂) oznacza długość nakładania się segmentu s₁ i s₂, maxov(s₁,s₂) całkowitą długość obszaru zajmowanego przez nakładające się segmenty, len(s₁) długość segmentu s₁, N i δ(s₁,s₂) zdefiniowane są poniżej

	(4)
	(5)

Miara SOV przyjmuje wartości od 0 do 1, 0 informuje o braku jakiegokolwiek nałożenia segmentów, 1 idealne nałożenie. Dla przypadku przewidywań mioglobiny opisanego powyżej miara SOV daje wynik 0.12.

Publikowane wyniki metod przewidywań struktur drugorzędowych trudno jest obiektywnie porównać z kilku powodów:

1. Metody testowane są na różnych zbiorach danych.

2. Wszystkie najnowsze metody korzystają z wielosekwencyjnego dopasowania, którego jakość silnie zależny od zastosowanej bazy danych sekwencji. W związku z tym nawet jeśli do testowania zostaną użyte dokładnie te same sekwencje białkowe, a metody opublikowane zostały w różnym czasie, to z uwagi na wygenerowanie MSA (Multiple Sequnece Alignment) na różnych bazach danych (powstałych w różnym czasie) metody nie mogą być obiektywnie porównane.

3. Metody korzystają z innych definicji struktur drugorzędowych.

W celu uniknięcia przedstawionych powyżej problemów stworzony został serwer EVA (Evaluation of Automatic protein structure prediction, niestety serwer ten w chwili obecnej jest nieaktywny, jednak ma zostać ponownie uruchomiony w 2011 roku) [], który do zarejestrowanych serwerów webowych wysyła w tym samym czasie sekwencje białek dodawanych do PDB. Serwer EVA wysyła tylko te sekwencje, które nie mają homologów wśród istniejących sekwencji w PDB, a więc nie mogły być użyte do utworzenia zbioru treningowego. Do wygenerowania struktur drugorzędowych ze struktury trzeciorzędowej serwer EVA wykorzystuje program DSSP. W chwili obecnej zarejestrowanych jest 13 serwerów do przewidywań struktur drugorzędowych, rezultaty porównań można znaleźć na stronie http://www.pdg.cnb.uam.es/eva/doc/intro_sec.html.

Przewidywanie struktur drugorzędowych białek globularnych

2.2.2.1. Pierwsza generacja metod przewidywania struktur drugorzędowych

Pierwsze generacja oparta jest na statystyce występowania pojedynczych aminokwasów w odpowiednich strukturach drugorzędowych. Jedną z najstarszych i jednocześnie najprostszych metod przewidywania struktur drugorzędowych jest metoda opisana przez Chou-Fasman-a []. Idea tej metody polega na przypisaniu każdemu aminokwasowi jego skłonności do pojawiania się w odpowiednich strukturach drugorzędowych (alfa helisach – p(alfa), beta kartach – p(beta), zagięciach – p(zagięcia) i pętlach). Skłonność ta wyznaczona została na podstawie znanych struktur białkowych. Przypisanie poszczególnym aminokwasom struktur drugorzędowych odbywa się na podstawie skłonności aminokwasów do występowania w strukturach drugorzędowych z uwzględnieniem sąsiadujących aminokwasów. Najpierw skanowana jest zadana sekwencja w poszukiwaniu regionów, dla których 4 z 6 po sobie następujących aminokwasów ma skłonność p(alfa) powyżej zadanego progu. Następnie takie regiony są rozbudowywane tak długo, jak długo średnia skłonność p(alfa) z 4 następujących po sobie aminokwasów daje wartość powyżej zadanego progu. Podobnie postępuje się z przewidywaniem beta kartek przy czym poszukiwane regiony zamiast 4 z 6 mają 3 z 5 aminokwasów. Przypisania aminokwasom klasę zagięcia następuje, gdy dla 4 kolejno po sobie następujących aminokwasów spełniona jest jedna z poniższych reguł:

• Otrzymuje się średnią p(zagięcia) większą od zadanego progu lub średnia p(zagięcia) jest większa od p(alfa) i p(beta).

• Iloczyn częstotliwości wystąpień dla tych 4 aminokwasów jest powyżej zadanego progu.

Wszystkie aminokwasy, którym nie zostały przypisane żadne struktury drugorzędowe zostają zaklasyfikowane jako pętle. Dokładność przewidywań tej metody przedstawiona w pracy była na poziomie 70-80%, jednak statystyki wyliczone zostały na podstawie tylko 15 sekwencji. Statystyki zostały przeliczone w późniejszych latach 1978 dla 29 i w 1989 dla 64 białek []. W chwili obecnej, gdy znana jest znacznie większa ilość struktur białkowych, dokładność tej metody jest na poziomie Q₃=50%.

2.2.2.2. Druga generacja metod przewidywania struktur drugorzędowych

Druga generacja metod oparta jest na statystykach wyliczonych na podstawie segmentów aminokwasów. Metody wchodzące w skład tej generacji bazują na: teorii informacji, własnościach fizykochemicznych oraz metodach uczenia maszynowego: sieciach neuronowych, metodach minimalno-odległościowych.

Jedną z najbardziej znanych metod tej generacji jest GOR []. Pierwsze wersje algorytmu GOR (GOR I i GOR II) przewidywały 4 struktury drugorzędowe. Oprócz standardowych 3 struktur dodatkowo wyróżnione zostały skręty występujące w pętlach, jednak od wersji III [] stosowano już 3 stanową predykcję. Predykcja w metodzie GOR polega na wyliczeniu ilości informacji dla zadanego aminokwasu z uwzględnieniem kontekstu, czyli 8 otaczających aminokwasów, co daje okna 17 aminokwasów ().

Rys. 2 Kolorem jasno szarym zaznaczone jest okno aminokwasów wokół pozycji wyróżnionej kolorem ciemno szarym

Ilość informacji wyliczana jest przy użyciu prawdopodobieństw wyznaczonych dla każdego typu struktury drugorzędowej. Z uwagi na małe zbiory danych dostępne w tych czasach zastosowano przybliżoną postać informacji zdefiniowaną jako suma informacji pochodząca od poszczególnych aminokwasów w oknie lub pochodząca od par aminokwasów w oknie, co przedstawiono poniższymi wzorami:

	(7)

lub

	(8)

gdzie S_j struktura drugorzędowa j-tego aminokwasu,, reprezentuje struktury drugorzędowe, do których j-ty aminokwas nie należy, R typ aminokwasów na odpowiedniej pozycji w oknie. Informacja I dla pojedynczego aminokwasu zdefiniowana jest w następujący sposób

	(9)

gdzie P(S_j=X|R_j) to prawdopodobieństwo warunkowe tego, że dla zadanego typu aminokwasu na pozycji j otrzymamy strukturę drugorzędową X, natomiast P(S_j=X) prawdopodobieństwo wystąpienia struktury drugorzędowej X.

Wzór (7) reprezentuje prostszą wersję aproksymacji informacji, która pozwala wyliczać ilość informacji dla zadanej struktury drugorzędowej reprezentowanej przez aminokwas j zależnej od typu aminokwasu i jego położenia w oknie.

Do wyliczenia informacji I ze wzoru (7) konieczne było wygenerowanie macierzy oceniających o wymiarach 17x20, w których kolumny reprezentują prawdopodobieństwa występowania danego aminokwasu na określonej pozycji. Macierz taka generowana jest dla każdej z klas (struktury drugorzędowej) osobno, na podstawie istniejących struktur białkowych. Przypisanie struktury drugorzędowej dla danego aminokwasu dokonywane jest na podstawie wartości informacji I, wygrywa ta klasa, dla której wartość ta jest największa. Jakość przewidywań dla GOR III była na poziomie 63%. W GOR próbowano również wykorzystać postać informacji przedstawioną we wzorze (20), jednak baza danych struktur trzeciorzędowych w tamtych czasach była na tyle mała, że wyznaczanie prawdopodobieństw dla par aminokwasów na podstawie częstotliwości wystąpień było obarczone dużym błędem.

2.2.2.3. Trzecia generacja metod przewidywania struktur drugorzędowych

W metodach 3 generacji głównym atrybutem było zastosowanie informacji ewolucyjnej zawartej w wielosekwencyjnym dopasowaniu (MSA - Multiple Sequence Alignment) oraz zastosowanie zaawansowanych metod uczenia maszynowego. MSA uzyskuje się poprzez nałożenie analizowanej sekwencji na inne sekwencje znajdujące się w bazie danych (np. NR) co umożliwia określenie zmienności aminokwasów na odpowiednich pozycjach w nałożeniu. Dzięki temu można wykryć obszary, w których zmienność aminokwasów jest niewielka, co może świadczyć o identycznych lub bardzo podobnych obszarach struktury trzeciorzędowej (drugorzędowej) białek.

Pierwszy algorytm tej generacji został opisany w [], opierał się na algorytmie GOR I, do którego dołączono informacje zawartą w dopasowaniu wielosekwencyjnym. Rezultatem było osiągnięcie dokładności na poziomie 66%, czyli poprawa o 9% w stosunku do GOR I. Od tego momentu wszystkie najlepsze metody korzystają z informacji ewolucyjnej.

Najbardziej znaną metodą 3 generacji, która jako pierwsza przekroczyła dokładność 70% jest PHD (Profile network from HeiDelberg) []. W PHD zastosowano dwa nowe elementy: zbalansowanie zbioru treningowego oraz specjalną architekturę predyktora.

Zbiór treningowy

Rozkład klas, czyli struktur drugorzędowych, w typowym zbiorze treningowym nie jest zbalansowany, gdyż w większości białek znacznie więcej jest aminokwasów związanych z pętlami niż alfa helisami czy beta kartkami. W wyniku treningu na takim zbiorze otrzymuje się przewidywania, które odpowiadają rozkładowi zbioru treningowego czyli beta kartki są bardzo słabo przewidywane. Jednym z elementów wprowadzonych przez PHD było zbalansowanie zbioru treningowego w taki sposób, żeby wszystkie klasy były równo licznie reprezentowane.

Zbiór treningowy został zbudowany w oparciu o przedstawicieli rodzin białkowych występujących w bazie danych HSSP [] (homology-derived structures of proteins). Pierwszym krokiem do utworzenia wektora treningowego było zbudowanie MSA w oparciu o bazę danych SWISS-PROT [], a następnie wygenerowanie profilu dla otrzymanego MSA.

Architektura PHD

PHD składa się z 3 poziomów: sieci neuronowej dokonującej odwzorowania sekwencji na strukturę, sieci neuronowej dokonującej odwzorowania struktury na strukturę oraz komitetu predyktorów, schemat przedstawiony został na . Na pierwszy poziomie trenowane są sieci, których wyniki przekazywane są do sieci drugiego poziomu. Każda sieć pierwszego poziomu ma swój odpowiednik (jedną sieć) na drugim poziomie, do którego przekazuje wyniki. Ostatnim etapem jest zebranie wyników ze wszystkich sieci poziomu drugiego i na ich podstawie dokonywana jest końcowa klasyfikacja.

Rys. 3. Schemat architektury PHD: sekwencja na strukturę, struktura na strukturę oraz komitet predyktorów

Odwzorowanie sekwencji na strukturę

Do odwzorowania sekwencji na strukturę zastosowana została sieć neuronowa z jedną warstwą ukrytą. Każdy j-ty aminokwas sekwencji na wejściu sieci reprezentowany jest w postaci okna 13 aminokwasów j-6...j...j+6. Zastosowanie informacji ewolucyjnej otrzymanej z MSA pozwala przekształcić każdy aminokwas w oknie w wektor 20 liczb będących odzwierciedleniem częstotliwości z jaką dany aminokwas jest w MSA zastępowany przez inne aminokwasy. Oprócz informacji z MSA każdy aminokwas w oknie charakteryzowany jest przez jeszcze jedną cechę pozwalającą zakodować informację zawiązaną z końcami białka (N-koniec, C-koniec). W wyniku użycia przedstawionej przestrzeni cech konieczne było zastosowanie sieci neuronowej z 273 węzłami wejściowymi. Na wyjściu sieci, dla każdego okna, wprowadzana była struktura drugorzędowa przypisana do środkowego aminokwasu w oknie. Struktura drugorzędowa zakodowana została w postaci wektora binarnego o długości 3 i tak H->1 0 0 (helisa), E-> 0 1 0 (beta kartka), C-> 0 0 1(pętla).

Reprezentacja aminokwasów w postaci okna powoduje, że każda sekwencja o długości l zamieniana jest na l wektorów uczących. Do utworzenia zbioru treningowego wykorzystano 130 sekwencji białkowych, z których uzyskano około 25000 wektorów.

Odwzorowanie struktury na strukturę

W odwzorowaniu sekwencji na strukturę dokonuje się klasyfikacji okien (segmentów), które są od siebie niezależne. Jednakże w przypadku struktur drugorzędowych pojawienie się alfa helisy na jakiejś pozycji implikuje duże prawdopodobieństwo wystąpienia tej samej struktury w bliskim otoczeniu. To oznacza, że pomiędzy segmentami z pierwszego odwzorowania powinna występować korelacja. Korelację tą wykorzystuje się w drugim odwzorowaniu użytym w PHD nazywanym odwzorowaniem struktury na strukturę. Podobnie jak w poprzednim etapie stosuje się sieć neuronową. Wejście sieci składa się z 51 węzłów wejściowych, do których przekazywane są wyniki przewidywań uzyskanych z odwzorowania sekwencji na strukturę. Wektor treningowy zbudowany jest na podstawie przewidywań struktur drugorzędowych dla okna 17 aminokwasów. Każdy aminokwas w oknie reprezentowany jest za pomocą trzech liczb otrzymanych z wyjścia pierwszej sieci, natomiast żądane wyjście danego wektora treningowego reprezentuje rzeczywistą strukturę drugorzędową środkowego aminokwasu w oknie zakodowaną w ten sam sposób co w odwzorowaniu sekwencji na strukturę.

Komitet sieci neuronowych

Trenowanie sieci neuronowych związane jest z minimalizacją funkcji kosztu, którą zazwyczaj jest funkcja błędu średniokwadratowego przedstawiona poniżej

	(10)

gdzie d_ij wartość i-tego wektora treningowego żądana na j-tym wyjściu sieci, f_j(x_i) wartość uzyskiwana na j-tym wyjściu sieci dla i-tego wektora treningowego, W parametry adaptacyjne, których wartość jest zmieniana podczas uczenia w taki sposób aby zminimalizować funkcję błędu. Do minimalizacji powyższej funkcji zazwyczaj stosuje się metody gradientowe, skutkiem czego uczenie kończy się najczęściej w jakimś minimum lokalnym. Ponieważ startuje się z losowych wartości parametrów W to powtarzanie treningu może prowadzić do innych minimów lokalnych i w końcowym efekcie do trochę innych przewidywań.

Kombinacja sieci neuronowych, tzw. komitet, wytrenowanych na różnych podzbiorach danych treningowych lub przy użyciu różnych architektur czy algorytmów uczenia zazwyczaj daje lepsze rezultaty niż pojedynczy klasyfikator. Najprostszym sposobem połączenia wyników różnych sieci neuronowych jest policzenie średniej arytmetycznej poszczególnych wyjść wszystkich sieci. W metodzie PHD wytrenowano 12 sieci, których wyniki zostały zsumowane. Zastosowanie komitetu pozwoliło na poprawę dokładności o 2% w stosunku do wyników pojedynczej sieci z drugiego poziomu.

Metody 3 generacji oceniane przez serwer EVA

Zgodnie z oceną meta serwera EVA jedną z najlepszych metod do przewidywania struktur drugorzędowych jest PSIPRED []. Bardzo dobre wyniki osiągnięte przez PSIPRED są skutkiem zastosowania macierzy PSSM wygenerowanej przy pomocy programu PSI-BLAST []. PSI-BLAST jako jeden z pierwszych programów zastosował iteracyjne nałożenie co znacząco poprawiło (w stosunku do innych metod) czułość w wykrywaniu sekwencji homologicznych. Do wygenerowania PSSM użyta została nieredundantna baza sekwencji, którą utworzono wybierając z istniejących publicznych baz danych sekwencje nie identyczne, które następnie zostały przefiltrowane w celu usunięcia obszarów o niskiej złożoności (zmienności) oraz obszarów transbłonowych. Architektura PSIPRED jest zbliżona do PHD, występuje odwzorowanie sekwencji na strukturę i struktury na strukturę, brak jest natomiast komitetu predyktorów. W obu odwzorowaniach zastosowano sieci neuronowe, w których algorytmem uczenia była wsteczna propagacja z momentem. Architektura sieci neuronowej w odwzorowaniu sekwencji na strukturę była następującą 315-75-3. Na wejście sieci wprowadzono fragment macierzy PSSM otrzymany z przesuwającego się okna 15 aminokwasów ().

Oprócz 300 węzłów wejściowych reprezentujących macierz PSSM, występuje jeszcze 1 węzeł, dla każdego aminokwasów w oknie, mający na celu określenie czy dany aminokwas znajduje się na jednym z końców (N lub C) białka. Do odwzorowania struktury na strukturę również zastosowano sieć neuronową, natomiast sekwencja reprezentowana jest przy użyciu przesuwającego się okna o tej samej długości co w pierwszym odwzorowaniu. Każdy aminokwas w oknie reprezentowany jest przez 3 wartości przewidywań struktur drugorzędowych z pierwszego etapu oraz znacznik N lub C końca białka co w sumie daje 60 węzłów warstwy wejściowej. W warstwie ukrytej użyto również 60 węzłów.

W pracy [] opisującej system Jnet przeanalizowano wpływ różnych metod nakładania wielosekwencyjnego oraz wyliczania profili z tych nałożeń na przewidywanie struktur drugorzędowych. Do przewidywania struktur drugorzędowych użyto sieci neuronowej i schematu predyktora identycznego do PHD czyli składającego się z trzech etapów. Do utworzenia nałożeń zastosowano PSIBLAST, CLUSTALW [] oraz AMPS [] natomiast profile wyznaczane były na trzy różne sposoby:

• Jako częstotliwości każdego z aminokwasów występujących na odpowiednich pozycjach w dopasowaniu wielosekwencyjnym tak jak to było zrobione w PHD.

• Każdy aminokwas w kolumnie wielosekwencyjnego dopasowania oceniany jest na podstawie macierzy podstawień BLOSUM62, wartości podstawień w danej kolumnie są uśredniane.

• Profile wygenerowane za pomocą HMMER [] .

• Profile PSSM wygenerowane przez PSIBLAST.

Wyniki przewidywań struktur drugorzędowych bardzo mocno zależą od użytych profili, różnice w dokładności przewidywań sięgały 6%, najgorsze rezultaty otrzymano przy reprezentacji wielosekwencyjnego nałożenia za pomocą BLOSUM62, najlepsze, gdy uśredniono profile z PSIBLAST oraz HMMER. Porównanie wpływu różnych profili na przewidywania struktur drugorzędowych pojawiło się również w pracy opisującej systemy SSpro-1.0 i SSpro-2.0 []. Zastosowane profile generowane były przez BLAST i PSI_BLAST, profile z PSI-BLAST poprawiły wyniki przewidywań o 1.5%.

W pracy [] zaprezentowany został predyktor SABLE, który podobnie jak PHD składa się z 3 poziomów, na 2 pierwszych poziomach wytrenowanych zostało 9 sieci. Przewidywania struktur drugorzędowych na poziomie trzecim uzyskane zostały poprzez zsumowanie wartości wyjść o wszystkich sieci drugiego poziomu i wybranie klasy, dla której wartość sumy jest największa. Poniższy wzór przedstawia wyznaczenie klasy C dla i-tego wektora x_i

	(11)

Schemat architektury sieci neuronowej dla predyktora na poziomie pierwszym przedstawiony jest na , jest to sieć rekurencyjna typu Elman z dwiema warstwami ukrytymi, w każdej z nich występuje 30 węzłów, oraz warstwą wyjściową zawierającą 9 węzłów tzw. rozszerzone wyjście reprezentujące struktury drugorzędowe dla trzech aminokwasów z sekwencji, które wcześniej zaproponowane zostało w pracy []. Optymalna architektura sieci została dobrana przy użyciu CV (walidacji krzyżowej). Wszystkie sieci trenowane były przy użyciu algorytmu RP (Resilient Propagation) [], zastosowanie tego algorytmu pozwoliło znacząco przyspieszyć czas uczenia w porównaniu do algorytmu standardowej wstecznej propagacji, jednakże w przypadku części sieci zamiast minimalizować standardową funkcję błędu przedstawioną we wzorze (10) zastosowana została zmodyfikowana funkcja postaci

	(12)
	(13)

gdzie i jest indeksem wektora treningowego natomiast trzy człony wewnątrz sumowania związane są z rozszerzonym wyjściem sieci i odpowiadają błędowi przedstawionemu we wzorze 25 dla i tego wektora treningowego reprezentującego klasę C_i oraz jego dwóch sąsiadów i-1 i i+1

Celem modyfikacji funkcji błędu jest zbalansowanie wpływu wektorów z poszczególnych klas na modyfikowaną funkcję. W PHD uzyskano to poprzez modyfikację zbioru treningowego w taki sposób aby wszystkie klasy były równo liczne w zbiorze treningowym. W SABLE dla każdej struktury drugorzędowej wprowadzono wagę α, której wartości dla poszczególnych klas były następujące: C – 0.3, H – 0.5 i E – 1.0.

W sieci Elman wewnętrzna reprezentacja problemu, wzbudzenia węzłów warstwy ukrytej, jest przekazywane ponownie na wejście odpowiedniej warstwy, co powoduje wykorzystanie stanów poprzednich w wyliczaniu stanów bieżących, dzięki czemu oprócz informacji znajdującej sie w pojedynczym oknie użyta zostaje również informacja pochodząca z okna poprzedniego.

Zastosowane zostało okno długości 11 aminokwasów. Podobnie jak we wcześniej opisanych metodach każdy aminokwas w oknie reprezentowany jest przez odpowiedni wiersz w macierzy PSSM. Poza tym każda pozycja w oknie charakteryzowana jest przez entropie na tej pozycji (wyliczanej na podstawie MSA stworzonego dla tej sekwencji). Dodatkowo zastosowano następujące cechy: poziom cystein w oknie, dla każdej pozycji w oknie uśredniona hydrofobowość i uśredniony rozmiar aminokwasów obserwowanych w MSA oraz dla każdej z trzech środkowych pozycji w oknie wprowadzono wektor binarny zawierający 5 bitów oznaczających przynależność danego aminokwasu do jednej z następujących grup: {A, E, L}, {V, I}, {S, N}, {P}, {G}. W efekcie zastosowania przedstawionej reprezentacji okna wymiar wektora wejściowego wynosił 269.

Zbiór treningowy utworzony został na podstawie reprezentantów rodzin występujących w bazie Pfam. Wybranych zostało ponad 800 sekwencji, które przy pomocy przesuwającego się okna zostały zamienione na około 200000 wektorów.

Przestrzeń cech na etapie drugim składa się z przewidywań sieci z etapu pierwszego oraz dodatkowo wprowadzone zostały cechy reprezentujące przewidywania relatywnej dostępności aminokwasów dla cząsteczek rozpuszczalnika otrzymane również przy pomocy serwera SABLE. Chociaż przestrzeń cech użyta do przewidywania relatywnej dostępności aminokwasów(RSA), metoda opisana w 3.2, jest identyczna do przestrzeni cech użytej do przewidywania SS, to użycie RSA poprawiło wyniki o 1%.

Znacząca poprawa generalizacji na trzecim poziomie jest możliwa tylko wtedy, gdy wyniki sieci z drugiego poziomu są jak najmniej skorelowane. W tym celu dla sieci z pierwszego poziomu stworzone zostały dwa zbiory treningowe składające się z tych samych okien tylko mające inne przestrzenie cech. W jednym zbiorze okna opisane były zbiorem cech przedstawionym powyżej, w drugim przestrzeń cech została okrojona do 232 wymiarów, w której występowała tylko informacja ewolucyjna, entropia oraz poziom cystein w oknie.

Opisaną przestrzeń cech o wymiarze 269 zastosowano również w pracy [], w której zbadano wpływ redukcji wymiarowości macierzy PSSM na przewidywania struktur drugorzędowych. Co prawda zredukowanie macierzy PSSM nie umożliwiło poprawy przewidywań, ale pozwoliło wprowadzić nową strategię uczenia sieci neuronowych.

Idea redukcji macierzy PSSM opiera się na traktowaniu kolumn macierzy (zdefiniowanych na ) jako pojedynczych cech. Liczba cech jest więc niezależna od wybranej długości okna i równa jest 20. Do selekcji cech wykorzystana została macierz PSSM pochodząca z połączenia wszystkich macierzy PSSM otrzymanych dla poszczególnych sekwencji białkowych, tak więc jej rozmiar wynosił 20 kolumn i około 200000 wierszy. Na przedstawiona jest analiza korelacji pomiędzy poszczególnymi kolumnami macierzy PSSM. Występuje szereg kolumn, które są bardzo silnie skorelowane np. kolumny reprezentujące podstawienie przez walinę (V) i izoleucynę (I) mają korelację 0.82, izoleucyna (I) i leucyna (L) 0.74, co nie powinno dziwić, gdyż aminokwasy te charakteryzują się dużym podobieństwem w budowie. Usunięcia korelacji w cechach można dokonać stosując PCA (Principal Component Analysis), przy czym konieczna jest reprezentacja PCA w pełnej postaci 220 na około 200000 wektorów co wykonane zostało w pracach [] []. Zastosowanie PCA poprawia co prawda wyniki na pierwszym poziomie predyktora (odwzorowanie sekwencji na strukturę), ale nie daje lepszych rezultatów w końcowej fazie przewidywań. Poza tym istotną wadą jest brak możliwości zrozumienia, która część macierzy PSSM jest najistotniejsza w przewidywaniach struktur drugorzędowych, gdyż potrzebna jest duża liczba wektorów własnych (około 60). Dlatego potraktowanie kolumn macierzy PSSM jako cech, czyli każda cecha składa się z wektora o wymiarze równym długości stosowanego okna, znacznie upraszcza interpretację wyników.

Do selekcji cech użyte zostały dwie proste i co najważniejsze szybkie (ze względu na liczbę wierszy) metody: t-statystyka oraz wzmocnienie informacji zdefiniowane w następujący sposób:

	(14)

gdzie μ_plus oznacza wartość średnią danej cechy wyznaczoną dla wszystkich elementów zbioru treningowego reprezentujących klasę plus, μ_minus wartość średnią dla klasy minus, σ odchylenie standardowe i n liczba elementów występujących w poszczególnych klasach.

Dla danej cechy F wzmocnienie informacji I(F) wyliczane jest przy użyciu entropii H(C) wyznaczonej dla rozkładu klas C występującego w danych, oraz entropii warunkowej H(C|F)

	(15)

Wartości wzmocnienia informacji oraz t-statystyki wyznaczone zostały dla wszystkich kolumn ciągu treningowego reprezentujących macierz PSSM (a więc 220 kolumn), a następnie wartości reprezentujące podstawienia określonego aminokwasu na poszczególnych pozycjach w oknie zostały zsumowane reprezentując w ten sposób miarę istotności kolumny reprezentującej podstawiany aminokwas. Na przedstawione zostały profile reprezentujące wartości wymienionych wcześniej miar dla poszczególnych aminokwasów (kolumn macierzy podstawień). Kształty profili są do siebie bardzo zbliżone szczególnie w obszarach o największych i najmniejszych wartościach. Najmniej istotnymi cechami są kolumny podstawieniowe następujących aminokwasów: tryptofan, cysteina, tyrozyna.

Na etapie odwzorowania sekwencji na strukturę, w oparciu o przedstawione oceny istotności cech, wykonana została selekcja cech. Stopniowo odrzucano cechy o najniższym stopniu istotności, po każdej eliminacji oceniana była generalizacja przy użyciu walidacji krzyżowej. Jako predyktora użyto prostej sieci neuronowej z jedną warstwą ukrytą. Najlepszą generalizację zarówno w CV jak i dla zewnętrznych zbiorów testowych uzyskano, gdy pozostawiono zaledwie 5 kolumn w PSSM odpowiadających podstawieniom przez prolinę (P), alaninę (A), walinę (V), kwas glutaminowy (E) oraz leucynę (L).

Niestety końcowy predyktor czyli komitet sieci neuronowych dawał znacząco gorsze wyniki od SABLE. Metoda selekcji kolumn macierzy PSSM okazała się jednak bardzo pomocna w budowaniu komitetu sieci neuronowych, gdyż odpowiedni dobór kolumn umożliwił wytrenowanie sieci, które są mniej skorelowane ze sobą niż w przypadku wygenerowaniu różnych sieci tylko w oparciu o walidację krzyżową i różne architektury sieci. Umożliwiło to zastosowanie jednego typu sieci (jednej architektury), co znacznie upraszcza budowanie końcowego predyktora. Pokazano również, że warunkiem otrzymanie dobrego wyniku końcowego jest budowanie zbiorów treningowych w oparciu o pięć najważniejszych kolumn macierzy PSSM wymienionych wcześniej, a w szczególności dotyczy to kolumny reprezentującej podstawienie przez prolinę (P). Wyniki otrzymane z takiego predyktora są właściwie identyczne z SABLE.

W ostatnich latach najlepszą dokładność odnotowaną przez serwer EVA osiąga Porter []. Powstał na bazie SSpro, który zbudowany jest przy użyciu sieci BRNN (Bidirectional Recurrent Neural Network) przedstawionej na . W architekturze tej można wyróżnić trzy składniki: środkowy zestaw węzłów T, do którego napływają sygnały z lokalnego okna (tak jak w tradycyjnej architekturze), oraz dwa boczne reprezentujące tzw. kontekst.

Rys. 8. Architektura sieci BRNN

Kontekst składa się z informacji otrzymanych w poprzednim przewidywaniu oraz następnym. Sygnał lokalny podawany do węzłów wejściowych T reprezentowany jest przez przesuwające się okno w macierzy PSSM.

Poprawa wyników przewidywań otrzymana w Porterze jest skutkiem kilku zmian w stosunku do SSpro:

• Zwiększenia zbioru treningowego (użyto 2171 białek mających w sumie 344653 aminokwasów).

• Likwidacja redundancji występującej w bazie NR użytej do utworzenia MSA.

• Zastosowania sieci z rekurencją przez co, jak uważają autorzy, uwzględnione zostają oddziaływania daleko zasięgowe, gdyż informacja na podstawie, której następuje klasyfikacja nie opiera się tylko na lokalnym oknie, lecz uwzględnia również sygnał pochodzący z dalszych obszarów sekwencji.

• Zwiększona została liczba aminokwasów w PSSM, zamiast standardowych 20 zastosowano 25: B (kwas asparaginowy lub asparagina), U (selenocysteina), X (nieznany), Z (kwas glutaminowy lub glutamina) i przerwa.

• Wprowadzenie mapowania struktury na strukturę.

Porter w nowszej wersji (Porter_H) [] dokonuje przewidywań struktur drugorzędowych w oparciu o szablony struktur trzeciorzędowych co znacznie poprawia jego przewidywania. Przestrzeń cech wykorzystywana przez Porter_H jest rozszerzeniem przestrzeni cech użytej w Porterze o 10 cech strukturalnych otrzymanych ze znalezionych szablonów. Pierwszych 8 tworzy wektor binarny, w którym zakodowane są uśrednione po wszystkich szablonach struktury drugorzędowe otrzymane z DSSP, natomiast pozostałe 2 reprezentują stopień pokrycia sekwencji przez szablon, dla której przewidywana jest struktura drugorzędowa oraz stopień identyczności aminokwasów w sekwencji i szablonach. Wszystkie cechy uśredniane są zgodnie z następującym schematem

	(16)

gdzie P liczba szablonów, c_p,j uśredniana cecha rozpatrywana dla j-tego aminokwasu oraz w_p jest wagą przypisaną p-temu szablonowi, id_p stopień identyczności sekwencyjnej pomiędzy sekwencją zadaną i sekwencją szablonu, q_p indeks jakości struktury krystalicznej.

2.2.2.4. Wyniki

Do porównania wybrane zostały 4 najlepsze metody przewidywania struktur drugorzędowych ocenione za pomocą serwera EVA. W każdej z tych metod zastosowane zostały sieci neuronowe, z tym, że w przypadku SAM-T99 [] do wykrywania homologów zastosowane zostały ukryte łańcuchy markowa. Ponieważ metody zarejestrowane zostały w serwerze EVA w różnym czasie porównywanie dokonywane jest parami na podstawie tylko tych sekwencji białkowych, dla których obie metody wygenerowały przewidywania. W poniższych tabelach przedstawione zostały średnie różnice porównywanych metod dla wybranej miary wraz z odchyleniem standardowym oraz wielkość próbki użyta w porównaniu (liczba sekwencji białkowych). Z poniższych tabel widać, że zastosowanie jednej miary jest mało miarodajne. Porównanie wyników pomiędzy PSiPred a SAM-T99 przy użyciu miary Q₃ jest bardzo wysokie natomiast przy użyciu miary SOV już nie.

Tabela 1 Porównanie parami 4 najlepszych metod (na podstawie oceny serwera EVA- http://www.pdg.cnb.uam.es/eva/sec/pairwise.html) przy użyciu miary Q₃

	PSipred	SAM-T99	PROFSec	SABLE
PSipred		-0.01±7.93 [ 277]	0.97 ± 8.87 [ 323]	3.16±5.77 [ 30]
SAM-T99	0.01±7.93 [ 277]		0.53 ± 8.63 [ 274]	0.58± 8.92 [ 30]
PROFSec	-0.97±8.87 [ 323]	0.53 ± 8.63 [ 274]		1.00±13.89 [ 26]
SABLE	-3.16 ± 5.77 [ 30]	-0.58 ± 8.92 [ 30]	-1.00 ± 13.89 [ 26]

Tabela 2 Porównanie parami 4 najlepszych metod (na podstawie oceny serwera EVA) przy użyciu miary SOV

	PSipred	PROFSec	SABLE	SAM-T99
PSipred		0.85 ± 13.48 [ 323]	0.88±12.22[72]	1.06±12.66 [277]
PROFSec	-0.85±13.48[323]		-1.78±15.54[80]	0.16±13.01[274]
SABLE	-0.88 ± 12.22 [ 72]	1.78±15.54[80]		0.20±10.94 [84]
SAM-T99	-1.06±12.66[277]	-0.16±13.01[274]	-0.20±10.94[84]

Białka transbłonowe

Rys. 9. Schemat reprezentujący różne struktury białek transbłonowych: pojedyncza alfa helisa przechodząca przez błonę komórkową, zespół kilku alfa helisy przechodzących przez błonę oraz zespół beta kartek tworzących strukturę beczki przechodzącą przez błonę

Wszystkie komórki otoczone są przez błonę komórkową zbudowaną głównie z lipidów. Lipidy tworzą dwuwarstwową strukturę, która jest hydrofilowa po zewnętrznych stronach i hydrofobowa pomiędzy warstwami. Cząsteczki białek transbłonowych są połączone z warstwą lipidów, w najprostszej postaci w cząsteczkach tych można wyróżnić trzy obszary: jeden obszar hydrofobowy zanurzony w membranie i 2 hydrofilowe znajdujące się na zewnątrz i wewnątrz komórki. W sytuacjach bardziej złożonych może występować większa ilość obszarów hydrofobowych. W większości przypadków segmenty białka wchodzące w błonę komórkową tworzę strukturę alfa helisy, zdarzają się jednak przypadki beta kartek.

Błona komórkowa pełni rolę bariery przepuszczalnej, której zadaniem jest chronić zawartość komórki przed środowiskiem zewnętrznym. Białka transbłonowe pełnią przede wszystkim rolę receptorów, czyli umożliwiają kontaktowanie się komórki ze środowiskiem zewnętrznym oraz rolę transportową.

W chwili obecnej w bazie danych białek transbłonowych MPtopo [] występuje tylko 185 struktur trzeciorzędowych z czego 150 to białka alfa helikalne a 35 białka, w których obszary transbłonowe tworzą beta kartki. Jest to mniej niż 1% liczby wszystkich dostępnych struktur trzeciorzędowych białek. Z drugiej jednak strony szacuje się, że białka transbłonowe stanowią od 20-30% liczby wszystkich białek w genomach. Tak niewielka liczba dostępnych struktur białek transbłonowych wynika z ogromnych trudności zastosowania technik eksperymentalnych, dlatego też bardzo istotne jest zastosowanie metod obliczeniowych pozwalających wyznaczyć struktury takich białek.

W białkach transbłonowych właściwie wyróżnia się tylko dwa stany drugorzędowe: alfa helisy i pętle lub beta kartki i pętle. Z tego powodu trenuje się niezależne predyktory do przewidywania alfa helis i beta kartek. Jako miarę dokładności predyktora stosuje się miarę Q₂, która jest podobna do Q₃ z tym, że występują tylko dwa stany.

Przewidywanie struktur drugorzędowych białek transbłonowych

Można wyróżnić trzy grupy metod przewidywania struktur drugorzędowych w białkach transbłonowych oparte na: skalach hydrofobowych i innych własnościach fizyko-chemicznych aminokwasów, obserwowanych preferencjach aminokwasów, profilach ewolucyjnych.

Pierwsza grupa metod wynika z obserwacji różnic na poziomie sekwencji pomiędzy białkami transbłonowymi i globularnymi. Białka transbłonowe zawierają dużą ilość aminokwasów hydrofobowych takich jak: leucyna, izoleucyna czy walina. Aminokwasy hydrofobowe występują w segmentach o długości od 20 do 30, są to te części białka, które zanurzone są w membranie komórkowej. W związku z tym jedną z możliwości wykrywania białek transbłonowych jest analiza rozkładu aminokwasów hydrofobowych występujących w strukturze pierwszorzędowej.

Kluczowym elementem w tej grupie metod jest zastosowanie odpowiedniej skali hydrofobowości aminokwasów. W literaturze występują dziesiątki różnych skal hydrofobowości, jednak w jednym z pierwszych algorytmów KD (Kyte-Doolittle) [] do przewidywania obszarów transbłonowych zastosowano skalę hydrofobowości KD zaproponowaną przez autorów. Do zidentyfikowania regionów transbłonowych zastosowano metodę okna, o długości 19 aminokwasów, przesuwającego się wzdłuż sekwencji. Wartość hydrofobowości aminokwasów w oknie jest sumowana, a gdy suma ta przekracza wartość progu T=1.6, obszar zdefiniowany przez okno klasyfikowany jest jako transbłonowy.

Metoda TopPred [] przewiduje tzw. pełną topologie białek transbłonowych czyli położenie helis transbłonowych oraz określenie po której stronie membrany znajdują się pętle łączące te helisy. W tym celu wykonywana jest analiza hydrofobowa, na podstawie której generowane są możliwe topologie oceniane następnie przez tzw. regułę dodatnie-wewnątrz. Najpierw za pomocą okna wyszukiwane są obszary o wysokiej hydrofobowości, przy użyciu skali GES (Goldman-Engelman-Steitz)[], następnie stosowana jest reguła stwierdzająca, że dodatnio naładowane aminokwasy (np. arginina czy lizyna) częściej występują po wewnętrznej stronie błony.

W metodzie SOSUI [] oprócz indeksu hydrofobowości KD zastosowano również indeks amfifilowy, reprezentujący amfifilowość polarnych łańcuchów bocznych. Indeks amfifilowy ma szczególne znaczenie dla określenia początków i końców helis transbłonowych

Chociaż metody oparte na skalach hydrofobowości są bardzo użyteczne to jednak mają jedną wadę, nie są w stanie odróżnić obszarów transbłonowych od globularnych zawierających dużą ilość aminokwasów hydrofobowych. W algorytmie PRED-TMR [], który można zaliczyć do drugiej grupy algorytmów, oprócz standardowych cech związanych z hydrofobowością wprowadzono cechy aminokwasów, których zadaniem jest charakteryzowanie początków i końców helis transbłonowych. Do tej grupy algorytmów można również zaliczyć metodę SPLIT [], w której użyto preferencji aminokwasów, wyliczonych na podstawie średnich wybranych własności (np. skali zagrzebania aminokwasu, czyli stopnia powierzchni aminokwasu występującego na powierzchni białka [], czy też skali hydrofobowości) 4 otaczających aminokwasów (okno o długości 8, wybrana pozycja jest pomijana) w 7 wyróżnionych strukturach drugorzędowych: N-koniec helisy, C-koniec helisy, środek helisy, krótkie helisy (o długości poniżej 5), beta kartki, skręty i obszary nieustrukturyzowane.

Kolejnym usprawnieniem było wprowadzenie w metodzie TMpred [] dodatkowych cech reprezentujących skłonności aminokwasów wyliczone dla obszarów transbłonowych i globularnych z uwzględnieniem aminokwasów sąsiadujących (podobnie jak w metodzie Chou-Fasmana stosowanej w białkach globularnych).

Trzecia grupa predyktorów bazuje na informacji ewolucyjnej. Publikowane rezultaty porównawcze [] pokazują, że ta grupa daje najlepsze wyniki. Jedną z najlepszych metod jest PHDhtm []. PHDhtm jest bardzo podobny do wcześniej opisanego systemu PHD przeznaczonego do przewidywania struktur drugorzędowych białek globularnych. Do wytrenowania PHDhtml wykorzystane zostały dwa zbiory treningowe zbalansowany i nie zbalansowany podobnie jak w PHD. Również architektura obu systemów jest identyczna czyli składa się z 3 poziomów: odwzorowania sekwencji na strukturę, struktury na strukturę i komitetu predyktorów. W odwzorowaniu sekwencji na strukturę sekwencja reprezentowana jest w postaci okna o długości 13 aminokwasów, różnice między tymi dwoma systemami występują na poziomie przestrzeni cech. W PHDhtml w odwzorowaniu sekwencji na strukturę każde okno reprezentowane jest przez grupę cech odzwierciedlającą informację lokalną dla tego okna oraz globalną dla całej sekwencji. Dla każdej pozycji w oknie wyliczana jest częstotliwość występowania 20 aminokwasów w przyrównaniu wielosekwencyjnym otrzymanym dla zadanej sekwencji, liczba wstawień i usunięć w przyrównaniu i waga zachowania charakterystyczna dla danej pozycji w oknie. Informacja globalna natomiast charakteryzowana była przez kompozycję aminokwasów w sekwencji, długość sekwencji oraz odległość pierwszego aminokwasu w oknie od N-końca białka i odległość ostatniego aminokwasu w oknie od C-końca. W drugim etapie czyli odwzorowaniu struktury na strukturę podobnie jak w pierwszym można wyróżnić dwa typy cech: lokalne i globalne. Cechy lokalne czyli odpowiedzi sieci dokonującej odwzorowania sekwencji na strukturę, każdy aminokwas w oknie reprezentowane jest przez dwie wartości otrzymane z wyjść tej sieci i wagę reprezentującą poziom konserwacji, oraz cechy globalne, które są identyczne do cech globalnych z odwzorowania sekwencji na strukturę. Wytrenowano 4 sieci, jedną na zbiorze zbalansowanym i jedną na nie zbalansowanym, reprezentujące odwzorowanie sekwencji na strukturę i odpowiadające im sieci odwzorowania struktury na strukturę. Wynik komitetu 2 sieci wyznaczony został poprzez wyliczenie średniej arytmetycznej dla poszczególnych węzłów wyjściowych, natomiast końcowa klasyfikacja reprezentowana jest przez węzeł wyjściowy o większej wartości. Ostatnim etapem przewidywań było zastosowanie filtrów, które reprezentują wiedzę o długościach helis transbłonowych. Helisy zbyt długie (>35) są dzielone na dwie o równych długościach lub jeżeli długość helisy L=22*n gdzie n=3,4... dzielona jest na n helis o długościach L/n, zbyt krótkie są usuwane lub rozszerzane tak aby ich długość była większa lub równa 17.

Znacząca poprawa wyników w porównaniu do metod opisanych wcześniej uzyskana została w metodzie MINNOU [], gdzie reprezentacja informacji ewolucyjnej zawarta została w postaci profilu przewidywanej relatywnej dostępności aminokwasów dla cząsteczek rozpuszczalnika (RSA – Relative Solvent Accessibility). Profil RSA składał się z ciągłych wartości RSA przewidzianych dla każdego aminokwasu w sekwencji przy użyciu serwera SABLE. Dla przedstawienia powiązania przewidywań RSA ze strukturami drugorzędowymi w białkach transbłonowych wykonano . Obszary zaznaczone żółtym kolorem odnoszą się do aminokwasów, które zgodnie z informacjami znajdującymi się w bazie MPtopo tworzą helisy transbłonowe. Pasek odcieni szarości poniżej sekwencji aminokwasów jest graficzną reprezentacją stopnia dostępności aminokwasów dla środowiska zewnętrznego, im ciemniejszy kolor tym aminokwas jest bardziej zagrzebany. Jak można zauważyć RSA dla aminokwasów z helis transbłonowych jest w większości przypadków bardzo mały (aminokwasy są zagrzebane), stąd widać, że wartości RSA są bardzo dobrą cechą do przewidywania obszarów transbłonowych. Ponieważ serwer SABLE optymalizowany był dla białek globularnych, w zbiorze treningowym nie występowały białka transbłonowe, to aminokwasy w helisach transbłonowych, które zanurzone są w lipidach, są przez SABLE traktowane jako zagrzebane, gdyż nie mają kontaktu z cząsteczkami wody, chociaż z punktu widzenia struktury białka zagrzebane nie są.

Architektura predyktora MINNOU jest identyczna z PHD. Na każdym etapie (odwzorowanie sekwencji na strukturę oraz struktury na strukturę) wytrenowanych zostało, przy użyciu algorytmu Rprop, wiele sieci neuronowych o różnych architekturach, których końcowa decyzja wyznaczana jest na podstawie prostego głosowania większościowego. Zmienność architektur uzyskano zmieniając szerokość okna przesuwanego wzdłuż sekwencji (od 8-18), co powodowało zmienną liczbę węzłów wejściowych, oraz zmieniając liczbie węzłów ukrytych. W warstwie wyjściowej natomiast występowały dwa neurony reprezentujące odpowiednie klasy.

W odwzorowaniu sekwencji na strukturę każdy aminokwas w oknie reprezentowany był przez 5 liczb: przewidywaną ciągłą wartość RSA, współczynnik zaufania dla przewidzianej wartości, otrzymane z serwera SABLE, prawdopodobieństwa przewidywania każdej ze struktur drugorzędowych również otrzymane z serwera SABLE. W odwzorowaniu struktury na strukturę wektor treningowy składał się z przewidywań predyktora pierwszego etapu, które przedstawione zostały jako uśrednione wartości, z wszystkich sieci biorących udział w konsensusie pierwszego etapu. Oprócz tych dwóch wartości każdy z aminokwasów reprezentowany jest przez te same cechy co na pierwszym poziomie predyktora MINNOU.

Jak pokazano w pracy [] serwer MINNOU uzyskuje dokładność 88% (miara Q₂), oszacowanie otrzymane w wyniku wykonania walidacji krzyżowej, natomiast dokładność otrzymywana przy użyciu profili hydrofobowych jest na poziomie 84%. Co więcej użycie dodatkowych cech (oprócz RSA i SS), które w poprzednich serwerach były bardzo pomocne: hydrofobowość oraz profil PSSM pogarszało (w niewielkim stopniu) końcowe rezultaty. MINNOU porównany został również z innymi metodami przy użyciu meta serwera TMH [] (idea podobna do serwera EVA). Metoda MINNOU, przy zastosowaniu miary Q₂, osiągnęła najlepsze wyniki (89% pozostałe metody na poziomie 80%).

Dostępność białka dla cząsteczek rozpuszczalnika

Pojęcie stopnia dostępności białka dla cząsteczek rozpuszczalnika zostało wprowadzone w pracy [] i oznacza powierzchnię białka opisaną przez wszystkie możliwe położenie cząsteczki wody będące w kontakcie (zdefiniowanym przez oddziaływanie Van der Waalsa’a) z atomami białka. Powierzchnia białka, a w szczególności stopnień dostępności aminokwasów dla cząsteczek rozpuszczalnika (SA - Solvent Accessibility), jest znaczącym uproszczeniem szczegółowej budowy trzeciorzędowej. Mimo to ma poważną zaletę – daje się dużo łatwiej przewidzieć w porównaniu do struktury trzeciorzędowej. Co więcej, mimo znacznego zubożenia informacji o budowie białka, przewidywane wartości SA są bardzo przydatne w przewidywaniu interfejsów białkowych [], struktur drugorzędowych [], [] oraz struktur trzeciorzędowych białek.

W przewidywaniach zamiast dostępności aminokwasów dla cząsteczek rozpuszczalnika dla wygody przewiduje się tzw. relatywną dostępność, która reprezentuje procentową dostępność danego aminokwasu, jest więc niezależna od jego wielkości. Relatywną dostępność dla cząsteczek rozpuszczalnika (RSA – Relative Solvent Accessibility) i-tego aminokwasu w białku definiuje się jako stosunek powierzchni tego aminokwasu, wystawionej dla rozpuszczalnika (SA_i), do maksymalnej możliwej powierzchni danego aminokwasu (MSAA_i) wystawionej dla cząsteczek rozpuszczalnika:

	(18)

RSA_i przyjmuje wartości od 0% do 100%, gdzie 0% odpowiada aminokwasowi w pełni zagrzebanemu, natomiast 100% w pełni dostępnemu dla cząsteczek rozpuszczalnika. Do wyliczenia SA_i konieczna jest znajomość struktury trzeciorzędowej. Istnieje wiele programów pozwalających wyliczyć SA_i jak np. DSSP, NACCESS [], jednak do najczęściej stosowanych należy DSSP. Wartości MSAA dla poszczególnych aminokwasów są wyliczone dla trójpeptydu Gly-X-Gly i można je znaleźć w pracach [], []. Niestety wartości przedstawione w tych pracach nie są identyczne w związku z tym wartości RSA dla różnych standardów MSAA mogą być różne co utrudnia porównanie dokładności metod przewidujących RSA.

Ocena dokładności przewidywań

Miarą dokładności przewidywań RSA jest MAE (Mean Absolute Error) oraz RMSD (Root Mean Squere Error), które zdefiniowane są poniżej

	(19)
	(20)

gdzie N liczba wektorów w zbiorze, d_i żądana wartość dla i tego wektora, o_i otrzymana wartość z predyktora. Ponieważ jednak duża część metod sprowadza problem przewidywania RSA do problemu klasyfikacyjnego to używane są również miary dokładności przewidywań opisane w przewidywaniach struktur drugorzędowych.

Przewidywania RSA

Istniejące metody do przewidywań RSA najczęściej dokonują podziału wartości RSA na kilka klas np. dwie: reprezentujące aminokwasy zagrzebane i na powierzchni, trzy: zagrzebane, częściowo na powierzchni i całkowicie na powierzchni jak również 10 klas. Niestety progi oddzielające różne klasy wybierane są w sposób arbitralny i w różnych pracach można spotkać zastosowane różne progi. Podział na klasy ma podłoże historyczne, i wywodzi się z analizy korelacji hydrofobowości z powierzchnią aminokwasu dostępną dla rozpuszczalnika [], []. Na rozpowszechnienie się podziału RSA na klasy miały też wpływ sukcesy w przewidywaniach struktur drugorzędowych, które z natury mają charakter dyskretny.

Najprostszym klasyfikatorem jest klasyfikator bazowy, którym zazwyczaj jest klasyfikator większościowy, jednak w pracy [] zaproponowano inną metodę wyznaczania poziomu odniesienia dla przewidywań RSA, która definiuje minimalny poziom jakości przewidywań. Algorytm ten opiera się na statystykach występowania aminokwasów w poszczególnych klasach RSA wyznaczonych na podstawie bazy białek, dla których znane są struktury trzeciorzędowe. Zastosowano podziały na klasy, których definicje można było najczęściej spotkać w literaturze. W związku z tym dwie klasy definiowane były przy użyciu następujących progów RSA: 9, 16, 20, 23, trzy klasy: 5-40, 9-36, 9-64 oraz dziesięć klas: 1-2-9-16-25-36-49-64-81, 0-3-8-15-23-32-42-53-67. Poziom dokładności uzyskanych przewidywań silnie zależy oczywiście od liczby klas oraz zastosowanych progów. W przypadku dwu klasowym waha się pomiędzy 68-71%, trzy klasowym między 50-61%, i dziesięcioklasowym 19-23%.

Bardziej zaawansowane metody, podobnie jak w przypadku przewidywań struktur drugorzędowych, wykorzystują metody uczenia maszynowego, dla których wektory treningowe zbudowane są w oparciu o kontekst zdefiniowany otaczającymi aminokwasami. Przykładem takiej metody jest Netasa [], system oparty na sieci neuronowej, w którym każdy aminokwas w sekwencji białkowej opisany jest przez okno o długości 17 aminokwasów. Każdy aminokwas w oknie zakodowany jest w postaci 21 wymiarowego wektora binarnego, 20 wymiarów reprezentujących poszczególne aminokwasy i wymiar 21, który koduje tzw. aminokwas nieznany (niestandardowy). Wektor binarny ma wartość jeden tylko na jednej pozycji reprezentującej dany aminokwas. Netasa przewiduje RSA dla problemu dwu klasowego, klasy zdefiniowane są przy użyciu następujących progów: 0, 5, 10, 25, 50. Dla każdego progu konieczne było wytrenowanie nowej sieci. Chociaż porównanie wyników Netasy z wcześniej omawianą metodą nie jest proste, gdyż nie są stosowane te same progi do definicji klas, to jednak łatwo zauważyć, że pomimo zastosowania okna i sieci neuronowej, nie ma znaczących różnic w przewidywaniach. Podobnie jest przy zastosowaniu innych metod, których predykcja opiera się na informacji zawartej w pojedynczej sekwencji, jak np. klasyfikator bayesowski [], czy też metoda stosująca teorię informacji []. Metoda oparta na teorii informacji jest właściwie identyczna z metodą GOR stosowaną do przewidywań struktur drugorzędowych, z tym, że do wyliczenia informacji zamiast wzoru (7) zastosowano wzór (8). Co prawda wyniki podane w pracy są znacząco lepsze (o 5%) jednak autorzy wyliczali RSA na podstawie własnego programu, sami potwierdzają, że różnice w RSA wyliczone na podstawie DSSP i ich metody są na poziomie 5-6%.

Znacząca poprawa wyników pojawiła się dopiero przy wprowadzeniu do przestrzeni cech informacji ewolucyjnej. Jak wykazano w pracy [] konserwacja RSA dla struktur homologicznych jest znacznie mniejsza w porównaniu do konserwacji struktur drugorzędowych. Najlepiej zachowane są aminokwasy tworzące rdzeń białka, a więc te, która mają niewielką wartość RSA, są zagrzebane, najgorzej te, które są na powierzchni. W przypadku trzy klasowym różnica w konserwacji aminokwasów między strukturami drugorzędowymi a RSA sięga 15%, co oczywiści skutkuje znacznie niższą poprawą jakości przewidywań RSA w stosunku do przewidywań SS.

PHDacc [] jest jedną z pierwszych metod wykorzystujących informację ewolucyjną. Przestrzeń cech w PHDacc jest bardzo podobna do przestrzeni cech zastosowanej w PHD. Podobnie jest z architekturą predyktora, która składa się z 2 warstw. Pierwsza warstwa zbudowana jest z sieci neuronowych wytrenowanych na dwóch różnych zbiorach treningowych wygenerowanych przy użyciu okien o długości 9 i 13 aminokwasów. Zastosowane zostały sieci neuronowe z jedną warstwą ukrytą oraz warstwą wyjściową zawierającą tą samą liczbę węzłów co liczba klas w rozpatrywanym problemie. Sieci trenowane były standardowym algorytmem wstecznej propagacji tak długo, aż osiągnięty został zadany próg dokładności. Przestrzeń cech dla tej warstwy składa się z informacji ewolucyjnej zakodowanej w identyczny sposób jak w PHD, oprócz tego występują cechy reprezentujące kompozycje aminokwasów (częstotliwość aminokwasów) w rozpatrywanej sekwencji, długość sekwencji, odległość okna od N i C końców białka.

Drugą warstwą jest komitet utworzony z wytrenowanych sieci, wyjściem komitetu jest średnia arytmetyczna ze wszystkich sieci. Dodatkowo wprowadzony został filtr polegający na uśrednieniu wartości przewidywań sąsiednich klas. PHDacc został wytrenowany dla problemu trzy klasowego (progi 16-36) oraz dziesięcioklasowego (progi 1-2-9-16-25-36-49-64-81). Zastosowanie informacji ewolucyjnej pozwoliło poprawić wyniki klasyfikacji o 3% w porównaniu do metod, w których informacja ta nie została użyta. Wykazano również, że długość okna kontekstowego, w przeciwieństwie do przewidywań struktur drugorzędowych, ma niewielki wpływ na jakość przewidywań RSA.

Pollastri [] zastosował swoją sieć BRNN, która odniosła duży sukces w przewidywaniach struktur drugorzędowych (PORTER), do przewidywania RSA. Na wejście sieci BRNN wprowadzona została macierz PSSM. Dla porównania użyto macierzy PSSM wygenerowanej przez BLAST oraz PSI-BLAST. Wygenerowanych zostało 20 predyktorów dla 20 problemów 2 klasowych zdefiniowanych przy progach: 0, 5, 10, ..., 95. Zaobserwowano poprawę w jakości przewidywań dzięki zastosowaniu PSSM z PSI-BLAST w stosunku do PSSM z BLAST w najlepszym przypadku o 0.8% czyli o połowę mniej niż w przewidywaniach SS. Oczywiście poprawa w przewidywaniach silnie zależy od zastosowanego progu definiującego klasę, najwyższa była dla progu 15%, natomiast powyżej progu 55% praktycznie nie zaobserwowano różnic.

Traktowanie problem przewidywania RSA jako problemu klasyfikacji wymusza podejmowanie nienaturalnych arbitralnych wyborów związanych z doborem liczby klas oraz progów definiujących te klasy. Również interpretacja wyników nie jest prosta, gdyż w klasyfikacji nie bierze się pod uwagę odległości od granic decyzyjnych, w związku z tym aminokwasy, dla których RSA wynosi 25.1 i 24.9, przy progu definiującym klasy zagrzebane i na powierzchni o wartości 25, należą do różnych klas. Co więcej, analiza średnich różnic RSA dla aminokwasów w strukturach homologicznych, które są nałożone w bazie danych Pfam, przedstawiona w pracy [] i zaprezentowana na pokazuje, że są one znaczne i zależne od regionów RSA. Zgodnie z wykresem przy RSA na poziomie 25, zróżnicowanie mierzone w RMSE wynosi od 20 do 15 w zależności od stopnia identyczności sekwencji, w związku z tym wymuszanie dwóch różnych klas w powyższym przykładzie jest nie fizyczne.

Rys. 12. Średnie RMSE pomiędzy aminokwasami strukturalnie równoważnymi dla różnych regionów SA. Pary aminokwasów zostały wzięte z bazy danych Pfam, która zawiera nałożenia wielosekwencyjne struktur homologicznych. Poziomo zróżnicowania RSA zależy od stopnia homologii, co reprezentują trzy krzywe na wykresie. Linia ciągła uzyskana ze struktur nałożonych w bazie danych Pfam reprezentuje odległe homologi, linia przerywana uzyskana dla struktur, których podobieństwo sekwencyjne wynosiło co najmniej 20%, nakładanych parami przy pomocy programu BLAST oraz linia kropkowana nałożenie struktur, których podobieństwo sekwencyjne wynosiło co najmniej 50%

Rozwiązanie tego problemu czyli wyjście poza klasyfikację w przewidywaniach RSA pojawiło się dość późno bo dopiero w 2003 roku w pracy []. Przestrzeń cech zastosowana w systemie opisanym w tej pracy jest identyczna z przestrzenią użytą przez system Netasa, czyli nie użyto informacji ewolucyjnej. System, o nazwie RVP-net, oparty jest na sieci neuronowej składającą się z jednej warstwy ukrytej zawierającej 2 węzły oraz warstwy wyjściowej, w której również użyto dwóch węzłów. Przewidywana wartość RSA wyliczana jest na podstawie wzbudzeń dwóch węzłów wyjściowych, jeżeli przez y₁ i y₂ oznaczyć wartości wzbudzeń węzłów wyjściowych to RSA wyliczane jest w następujący sposób

	(21)

Przewidywania RVP-net były oceniane na trzech niezależnych zbiorach danych, średnia dokładność jest na poziomie MAE=19.0, współczynnik korelacji 0.48. Jakość przewidywań zależna jest od obserwowanej wartości RSA im mniejsza wartość RSA tym lepsze przewidywania, najgorsze przewidywania są dla aminokwasów o dużym RSA czyli aminokwasów na powierzchni, dla których błąd sięga nawet 70%.

Znacząca poprawa przewidywań dla ciągłych wartości RSA zaprezentowana została w pracy [], w której uzyskano średnią wartość MAE=15.6 natomiast współczynnik korelacji 0.64. W celu poprawy wyników przewidywań RSA zastosowano informacje ewolucyjną. Przestrzeń cech zastosowana do wytrenowania predyktora jest identyczna z przestrzenią cech opisaną dla predyktora SABLE przeznaczonego do przewidywań struktur drugorzędowych. Jako predyktora użyto sieci neuronowej typu Elman posiadającej tylko jedno wyjście. Dokonana jednak została istotna zmiana w funkcji kosztu minimalizowanej podczas treningu. Zamiast standardowej funkcji błędu, używanej w sieciach neuronowych, zastosowano modyfikację uwzględniającą zmienność wartości RSA przedstawioną na .

	(22)

gdzie α(o_i) są to wagi zależne od obserwowanych RSA uwzględniające poziom wariancji RSA, y_i przewidywana wartość RSA i o_i wartość obserwowana ( żądana na wyjściu sieci). Wartości wag wyznaczone zostały dla 10 przedziałów RSA ([0,10] (10,20], (20,30] …), jako 1/RMSE*100, Modyfikacja funkcji błędu powoduje, że błędy dla małych wartości RSA mają większe znaczenie niż dla dużych.

Zastosowanie sieci neuronowych w problemach klasyfikacyjnych pozwala w łatwy sposób wyznaczyć współczynniki zaufania do przewidywań, poprzez podzielenie wartości wzbudzenia węzła wyjściowego przez sumę wzbudzeń wszystkich węzłów wyjściowych. W przypadku aproksymacji takie podejście nie jest możliwe, gdyż występuje tylko jeden węzeł wyjściowy, poza tym ocena błędu zależna jest od rzeczywistych wartości RSA. W związku z tym w celu wyprowadzenia współczynników zaufania stworzony został pomocniczy problem klasyfikacyjny, którego zadaniem było wyznaczenie współczynnika zaufania dla danego przewidywania. Wyróżnione zostały dwie klasy: przewidywanie poprawne i niepoprawne. Przewidywanie uznawane jest jako poprawne, gdy różnica pomiędzy wartością przewidzianą RSA a obserwowaną była mniejsza niż średnia różnica otrzymywana dla równoważnych aminokwasów z sekwencji pochodzących z tej samej rodziny pomnożona przez współczynnik 1.5. Dla przykładu, jeśli wartość obserwowana RSA wynosi 1%, natomiast przewidywana 9%, to takie przewidywanie traktowane jest jako poprawne, gdyż średnia różnica dla RSA z przedziału [0,10] wynosi 8.1. Utworzony został zbiór treningowy składający się z tego samego zbioru wektorów, które zostały użyte do wytrenowania RSA z tym, że każdemu wektorowi została przypisana klasa reprezentująca poprawność przewidywań tego wektora przez predyktor RSA. Tak utworzony zbiór treningowy posłużył do wytrenowania kolejnych sieci, których architektura była następująca 278-50-2, 278-30-2, 278-20-2. Współczynnik zaufania do przewidywań RSA wyliczony jest na podstawie wartości węzłów wyjściowych sieci zgodnie z poniższym wzorem

	(23)

gdzie sumowanie przebiega po wszystkich wytrenowanych sieciach, o^corr wartość wzbudzenia węzła wyjściowego reprezentującego klasę poprawnego przewidywania RSA i o^inc wartość wzbudzenia węzła reprezentującego klasę niepoprawnego przewidywania RSA.

Podsumowanie

Przewidywanie struktur drugorzędowych oraz wartości RSA na podstawie sekwencji aminokwasów ma bardzo istotne znaczenie w przewidywaniu struktur trzeciorzędowych białek. W chwili obecnej większość metod przewidywania struktur trzeciorzędowych korzysta z tych przewidywań. Jednak historia rozwoju metod przewidywania struktur drugorzędowych pokazuje, że utworzenie metody uzyskującej najlepsze rezultaty nie jest łatwe, gdyż wymaga perfekcyjnego wykonania wszystkich etapów budowy predyktora tj.: doboru przestrzeni cech, budowy zbioru treningowego, wyboru metody uczenia maszynowego, procedury uczenia. Różnice pomiędzy istniejącymi metodami wynikają ze sposobu wprowadzania wiedzy o problemie, która może być przekazywana w postaci konstruowanej przestrzeni cech zbioru treningowego, wyboru wektorów do zbioru treningowego lub poprzez odpowiednią modyfikację algorytmu uczącego. W przypadku najlepszych metod różnice występują prawie na wszystkich etapach budowy predyktora, są jednak również podobieństwa: przestrzeń cech składa się głównie z informacji ewolucyjnej zapisanej w macierzy PSSM, najczęściej wybieranym algorytmem uczenia maszynowego są sztuczne sieci neuronowe, gdyż doskonale nadają się do problemów, w których występują duże zbiory treningowe. Zastosowanie sieci neuronowych ma jednak również liczne wady, do których niewątpliwie należy zaliczyć trudności w zrozumieniu na jakiej podstawie dokonywane są takie a nie inne decyzje, lub też inaczej mówiąc, które cechy mają najistotniejsze znacznie w przewidywaniach. W ramach badań nad tym problemem przedstawiona została selekcji cech macierzy PSSM, która pozwoliła określić elementy macierzy PSSM mające największy wpływ na przewidywaniach.

Przewidywane wartości RSA otrzymane z serwera SABLE wykorzystane zostały w budowie przestrzeni cech dla predyktora przeznaczonego do przewidywań struktur drugorzędowych białek globularnych (predyktor ten stał się częścią serwera SABLE) oraz dla serwera MINNOU (minnou.cchmc.org) przeznaczonego do przewidywań struktur drugorzędowych białek transbłonowych. W metodzie MINNOU zamiast informacji ewolucyjnej reprezentowanej przez macierz PSSM, stosowanej przez wszystkie inne metody, użyto przewidywanych wartości RSA oraz SS otrzymanych z serwera SABLE. To pozwoliło znacząco zmniejszyć przestrzeń cech i poprawić wyniki przewidywań.

O dużym zainteresowaniu tego typu metodami może świadczyć fakt, że serwer SABLE od chwili stworzenia do dnia dzisiejszego został użyty około 180000 razy natomiast serwer MINNOU 45000.

Bibliografia

[1] Adamczak R.: Dimensionality reduction of pssm matrix and its influence on secondary structure and relative solvent accessibility predictions, in World Academy Of Science, Engineering And Technology, vol. 58, October 2009, pp. 657–664.

[2]  Adamczak R., Porollo A., Meller J. : Accurate prediction of solvent accessibility using neural networks-based regression., Proteins, 56 (2004), pp. 753–767.

[3] Adamczak R., Porollo A., Meller J.: Combining prediction of secondary structure and solvent accessibility in proteins., Proteins, 59 (2005), pp. 467–475.

[4] Ahmad S., Gromiha M. M.: Netasa: neural network based prediction of solvent accessibility., Bioinformatics, 18 (2002), pp. 819–824.

[5] Ahmad S., Gromiha M. M., Sarai A., Real value prediction of solvent accessibility from amino acid sequence., Proteins, 50 (2003), pp. 629–635.

[6] Altschul S. F., Madden T. L., Schoffer A. A., Zhang J., Zhang Z., Miller W., Lipman D. J.: Gapped blast and psi-blast: a new generation of protein database search programs., Nucleic Acids Res, 25 (1997), pp. 3389–3402.

[7] Anfinsen C. B., Haber E., Sela M., White F. H.: The kinetics of formation of native ribonuclease during oxidation of the reduced polypeptide chain., Proc Natl Acad Sci U S A, 47 (1961), pp. 1309–1314.

[8] Bairoch A., Apweiler R.: The swiss-prot protein sequence data bank and its supplement trembl in 1999., Nucleic Acids Res, 27 (1999), pp. 49–54.

[9] Barton G. J.: Protein multiple sequence alignment and flexible pattern matching., Methods Enzymol, 183 (1990), pp. 403–428.

[10] Benkert P., Tosatto S. C. E., Schomburg D.: Qmean: A comprehensive scoring function for model quality assessment., Proteins, 71 (2008), pp. 261–277.

[11] Bernstein F. C., Koetzle T. F., Williams G. J., Meyer E. F., Brice M. D., Rodgers J. R., Kennard O., Shimanouchi T., Tasumi M.: The protein data bank: a computer-based archival file for macromolecular structures., J Mol Biol, 112 (1977), pp. 535–542.

[12] Bowie J. U., Lathy R., Eisenberg D., A method to identify protein sequences that fold into a known three-dimensional structure., Science, 253 (1991), pp. 164–170.

[13] Cao B., Porollo A., Adamczak R., Jarrell M., Meller J.: Enhanced recognition of protein transmembrane domains with prediction-based structural profiles., Bioinformatics, 22 (2006), pp. 303–309.

[14] Chen H., Gu F., Huang Z.: Improved chou-fasman method for protein secondary structure prediction., BMC Bioinformatics, 7 Suppl 4 (2006), p. S14.

[15] Chothia C.: The nature of the accessible and buried surfaces in proteins., J Mol Biol, 105 (1976), pp. 1–12.

[16] Chou P. Y., Fasman G. D.: Prediction of protein conformation., Biochemistry, 13 (1974), pp. 222–245.

[17] Cuff J. A., Barton G. J.: Application of multiple sequence alignment profiles to improve protein secondary structure prediction., Proteins, 40 (2000), pp. 502–511.

[18] Eddy S. R., Profile hidden markov models., Bioinformatics, 14 (1998), pp. 755–763.

[19] Engelman D. M., Steitz T. A., Goldman A.: Identifying nonpolar transbilayer helices in amino acid sequences of membrane proteins., Annu Rev Biophys Biophys Chem, 15 (1986), pp. 321–353.

[20] Eyrich V. A., Marti-Renom M. A., Przybylski D., Madhusudhan M. S., Fiser A., Pazos F., Valencia A., Sali A., Rost B.: Eva: continuous automatic evaluation of protein structure prediction servers., Bioinformatics, 17 (2001), pp. 1242–1243.

[21] Finn R. D., Mistry J., Tate J., Coggill P., Heger A., Pollington J. E., Gavin O. L., Gunasekaran P., Ceric G., Forslund K., Holm L.: Sonnhammer E. L. L., Eddy S. R., Bateman A., The pfam protein families database., Nucleic Acids Res, (2009).

[22] Frishman D., Argos P.: Knowledge-based protein secondary structure assignment., Proteins, 23 (1995), pp. 566–579.

[23] Gibrat J. F., Garnier J., Robson B.: Further developments of protein secondary structure prediction using information theory. new parameters and consideration of residue pairs., J Mol Biol, 198 (1987), pp. 425–443.

[24] Hirokawa T., Boon-Chieng S., Mitaku S.: Sosui: classification and secondary structure prediction system for membrane proteins., Bioinformatics, 14 (1998), pp. 378–379.

[25] Hofmann K., Stoffel W.: Tmbase - a database of membrane spanning proteins segments., Biol. Chem. Hoppe-Seyler, 374 (1993), p. 166.

[26] Hubbard S. J., Thornton M.: Naccess, tech. report, Department of Biochemistry and Molecular Biology, University College: London, 1993.

[27] Janin J.: Surface and inside volumes in globular proteins., Nature, 277 (1979), pp. 491–492.

[28] Jayasinghe S., Hristova K., White S. H.: Mptopo: A database of membrane protein topology., Protein Sci, 10 (2001), pp. 455–458.

[29] Jones D.T.: Genthreader: an efficient and reliable protein fold recognition method for genomic sequences., J Mol Biol, 287 (1999), pp. 797–815.

[30] Jones D.T.: Protein secondary structure prediction based on position-specific scoring matrices., J Mol Biol, 292 (1999), pp. 195–202.

[31] Juretić D., Lee B., Trinajstić N., Williams R. W.: Conformational preference functions for predicting helices in membrane proteins., Biopolymers, 33 (1993), pp. 255–273.

[32] Kabsch W., Sander C.: Dictionary of protein secondary structure: pattern recognition of hydrogen-bonded and geometrical features., Biopolymers, 22 (1983), pp. 2577–2637.

[33] Karplus K., Barrett C.: Hughey R., Hidden markov models for detecting remote protein homologies., Bioinformatics, 14 (1998), pp. 846–856.

[34] Kernytsky A., Rost B.: Static benchmarking of membrane helix predictions., Nucleic Acids Res, 31 (2003), pp. 3642–3644.

[35] Kyte J., Doolittle R. F.: A simple method for displaying the hydropathic character of a protein., J Mol Biol, 157 (1982), pp. 105–132.

[36] Lundstrom J., Rychlewski L., Bujnicki J., Elofsson A.: Pcons: a neural-network-based consensus predictor that improves fold recognition., Protein Sci, 10 (2001), pp. 2354–2362.

[37] Melo J. C. B., Cavalcanti G. D. C., Guimaraes K. S.: Pca feature extraction for protein structure prediction, in Proc. International Joint Conference on Neural Networks, vol. 4, 20–24 July 2003, pp. 2952–2957.

[38] Naderi-Manesh H., Sadeghi M., Arab S., Movahedi A. A. M.: Prediction of protein surface accessibility with information theory., Proteins, 42 (2001), pp. 452–459.

[39] Pasquier C., Promponas V. J., Palaios G. A., Hamodrakas J. S., Hamodrakas S. J.: A novel method for predicting transmembrane segments in proteins based on a statistical analysis of the swissprot database: the pred-tmr algorithm., Protein Eng, 12 (1999), pp. 381–385.

[40] Pawlowski M., Gajda M. J., Matlak R., Bujnicki J. M.: Metamqap: a meta-server for the quality assessment of protein models., BMC Bioinformatics, 9 (2008), p. 403.

[41] Petersen T. N., Lundegaard C., Nielsen M., Bohr H., Bohr J., Brunak S., Gippert G. P., Lund O.: Prediction of protein secondary structure at 80% accuracy., Proteins, 41 (2000), pp. 17–20.

[42] Pollastri G., Baldi P., Fariselli P., Casadio R.: Prediction of coordination number and relative solvent accessibility in proteins., Proteins, 47 (2002), pp. 142–153.

[43] Pollastri G., Martin A. J. M., Mooney C., Vullo A.: Accurate prediction of protein secondary structure and solvent accessibility by consensus combiners of sequence and structure information., BMC Bioinformatics, 8 (2007), p. 201.

[44] Pollastri G., McLysaght A.: Porter: a new, accurate server for protein secondary structure prediction., Bioinformatics, 21 (2005), pp. 1719–1720.

[45] Pollastri G., Przybylski D., Rost B., Baldi P.: Improving the prediction of protein secondary structure in three and eight classes using recurrent neural networks and profiles., Proteins, 47 (2002), pp. 228–235.

[46] Porollo A., Meller J.: Prediction-based fingerprints of protein-protein interactions., Proteins, 66 (2007), pp. 630–645.

[47] Porollo A. A. , Adamczak R., Meller J., Polyview: a flexible visualization tool for structural and functional annotations of proteins., Bioinformatics, 20 (2004), pp. 2460–2462.

[48] Richardson C. J., Barlow D. J.: The bottom line for prediction of residue solvent accessibility., Protein Eng, 12 (1999), pp. 1051–1054.

[49] Riedmiller M., Braun H.: Rprop- a fast adaptive learning algorithm, tech. report, Proc. of ISCIS VII, Universitat, 1992.

[50] Rose G. D., Geselowitz A. R., Lesser G. J., Lee R. H., Zehfus M. H.: Hydrophobicity of amino acid residues in globular proteins., Science, 229 (1985), pp. 834–838.

[51] Rost B., Review: protein secondary structure prediction continues to rise., J Struct Biol, 134 (2001), pp. 204–218.

[52] Rost B., Casadio R., Fariselli F., Sander C.: Transmembrane helices predicted at 95% accuracy., Protein Sci, 4 (1995), pp. 521–533.

[53] Rost B., Sander C.: Prediction of protein secondary structure at better than 70% accuracy., J Mol Biol, 232 (1993), pp. 584–599.

[54] Rost B., Sander C.: Conservation and prediction of solvent accessibility in protein families., Proteins, 20 (1994), pp. 216–226.

[55] Sander C., Schneider R.: Database of homology-derived protein structures and the structural meaning of sequence alignment., Proteins, 9 (1991), pp. 56–68.

[56] Simas G. M., Botelho S. S. C., Grando N., Colares R. G.: Dimensional reduction in the protein secondary structure prediction — nonlinear method improvements, in Innovations in Hybrid Intelligent Systems, vol. 44/2008 of Advances in Soft Computing, Springer Berlin / Heidelberg, 2008, pp. 425–432.

[57] Thompson J. D., Higgins D. G., Gibson T. J.: Clustal w: improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice., Nucleic Acids Res, 22 (1994), pp. 4673–4680.

[58] Thompson M. J., Goldstein R. A.: Predicting solvent accessibility: higher accuracy using bayesian statistics and optimized residue substitution classes., Proteins, 25 (1996), pp. 38–47.

[59] von Heijne G.: Membrane protein structure prediction. hydrophobicity analysis and the positive-inside rule., J Mol Biol, 225 (1992), pp. 487–494.

[60] Zemla A., Venclovas C., Fidelis K., Rost B.: A modified definition of sov, a segment-based measure for protein secondary structure prediction assessment., Proteins, 34 (1999), pp. 220–223.

[61] Zvelebil M. J., Barton G. J., Taylor W. R., Sternberg M. J., Prediction of protein secondary structure and active sites using the alignment of homologous sequences., J Mol Biol, 195 (1987), pp. 957–961.

---

1. Katedra Informatyki Stosowanej, Uniwersytet Mikołaja Kopernika, ul. Grudziądzka 5, 87-100 Toruń; [e-mail]↩