W niniejszej pracy, przedstawiamy zintegrowane podejście łączące analizę statystyczną i metodę klasyfikacji hiper-boxów opartą na MILP do wczesnego przewidywania działania leków skierowanych na Ache, BZR, COX-2, DHFR_TG, DHFR_RL, DHFR_PC i wreszcie Cytochrom P450 C17.
Podejście zastosowane w niniejszej pracy składa się z pięciu głównych kroków. W pierwszym kroku, struktury molekularne kandydatów na leki są budowane i optymalizowane przez Marvin Sketch. Następnie, deskryptory molekularne tych kandydatów na leki są uzyskiwane przy użyciu serwera internetowego E-Dragon . Drugi krok polega na zbudowaniu modelu regresji przy użyciu PLS, co pozwoli na wyselekcjonowanie najbardziej znaczących deskryptorów. Następnie kandydaci na leki są klasyfikowani na podstawie najbardziej znaczących deskryptorów, które zostały uzyskane w poprzednim kroku, przy użyciu metody hiper-boxów opartej na MILP. Ta pierwotna klasyfikacja może skutkować relatywnie niższą dokładnością klasyfikacji ze względu na istnienie kilku nieistotnych deskryptorów w modelu, dlatego też przeprowadzana jest analiza testowania istotności w celu określenia nieistotnych deskryptorów, które mogą zakłócić naszą dokładność klasyfikacji w czwartym kroku. Jeśli w modelu występują deskryptory nieistotne, zastępujemy je deskryptorami bardziej istotnymi, a następnie wracamy do kroku trzeciego, w którym ponownie klasyfikujemy działania leków za pomocą nowego modelu uzyskanego w kroku piątym. Po przeprowadzeniu testów istotności, jeśli wszystkie deskryptory są znaczące, budujemy model z najbardziej znaczącymi i podajemy wyniki klasyfikacji.
Wykorzystujemy algorytm iteracyjny, dzięki czemu niektóre kroki mogą być powtarzane, gdy testy istotności dają niezadowalające wyniki dla wybranych deskryptorów danego modelu. Mniej znaczące deskryptory są zastępowane bardziej znaczącymi, co wpływa na ostateczną klasyfikację leków w każdej iteracji, a tym samym zwiększa sukces badania. Zarys naszej metody przedstawiono na rysunku 1.
Zarys podejścia do klasyfikacji.
Zestawy danych
Zastosowaliśmy nasz algorytm do szeroko znanych zestawów danych QSAR dostępnych w literaturze. Reduktaza Dihydrofolianu (DHFR), Acetylocholinoesteraza (AchE), Receptor Benzodiazepinowy (BZR) oraz inhibitory cyklooksygenazy-2 (COX-2) zostały wykorzystane do klasyfikacji. Przedstawiamy również nowy zestaw danych inhibitorów Cytochromu P450 C17, które zaczerpnęliśmy z literatury i obliczyliśmy ich struktury 3D.
Siedem zestawów danych zostało wykorzystanych do walidacji naszej metodologii poprzez zastosowanie algorytmu na tych dużych i znanych zestawach danych oraz porównanie naszej dokładności klasyfikacji na tych zestawach danych z innymi szeroko stosowanymi klasyfikatorami dostępnymi w pakiecie eksploracji danych WEKA. Reprezentatywne związki z każdego zbioru danych są pokazane na Rysunku 2. Eksperymentalne wartości IC50 dla zestawu inhibitorów reduktazy dihydrofolianowej (DHFR) zostały obliczone i podane dla enzymu DHFR z trzech różnych gatunków: P. carinii (PC), T. gondii (TG) i wątroby szczura (RL), gdzie aktywność inhibitorów DHFR wobec enzymów z różnych gatunków jest różna. Dlatego też w naszej pracy aktywność inhibitorów wobec enzymów tych trzech gatunków dla inhibitorów DHFR badana jest oddzielnie. Zestaw 397 inhibitorów reduktazy dihydrofolianowej (DHFR) zastosowano dla DHFR P. carinii z wartościami IC50 od 0,31 nM do 3700 μM, zestaw 378 inhibitorów zastosowano dla DHFR T. gondii z wartościami od 0,88 nM do 392 μM, a 397 inhibitorów zastosowano dla DHFR wątroby szczura z wartościami od 0,156 nM do 7470 μM. Zastosowano zestaw 111 inhibitorów acetylocholinoesterazy (AchE) z obliczonymi doświadczalnie wartościami IC50, podawanymi w zakresie od 0,3 nM do 100 μM . Zestaw danych dotyczących inhibitorów receptora benzodiazepinowego (BZR) obejmował 163 inhibitory, których wartości IC50 obliczono doświadczalnie w zakresie od 1,2 nM do 5 μM. W zestawie 322 cząsteczek inhibitorów cyklooksygenazy-2 (COX2) uzyskano takie wartości IC50 od 1 nM do 100 μM . Zestawy QSAR użyte w tej pracy zostały również wykorzystane w badaniu porównawczym metod QSAR przeprowadzonym przez Sutherlanda i in. Porównaliśmy również wartości R2 naszych modeli deskryptorów 3D, które zostały obliczone przez przebieg Minitab PLS w pierwszej fazie naszego algorytmu, z podanymi wartościami R2 przez Sutherlanda i in. dla kilku modeli PLS na tych samych zestawach danych.
Reprezentatywne związki z każdego z danych QSAR.
Budowanie struktury i otrzymywanie modelu deskryptorowego
Jak przedstawiono powyżej, w naszym badaniu pierwszym krokiem jest znalezienie deskryptorów molekularnych dla kandydatów na leki. Dlatego też, Marvin Sketch został użyty do obliczenia struktur molekularnych każdego kandydata na lek, które powinny być zbudowane poprzez zbudowanie ich struktury i zoptymalizowanie ich energii poprzez minimalizację w celu określenia ich potwierdzenia w przestrzeni 3-D. Następnie, zoptymalizowane struktury 3-D są ładowane do E-Dragon i deskryptory molekularne są obliczane przy użyciu serwera internetowego.
E-Dragon proponuje wiele bloków deskryptorowych, z których każdy zawiera parametry opisujące charakterystykę cząsteczek, a te, które zostały wykorzystane w niniejszej pracy można wymienić następująco: deskryptory konstytucyjne (48), deskryptory topologiczne (119), indeksy połączalności (33), indeksy informacyjne (47), indeksy addytywności krawędzi (107), indeksy ładunku topologicznego (21), deskryptory geometryczne (74), deskryptory 3D-MoRSE (160), liczebności grup funkcyjnych (154), fragmenty atomocentryczne (120), właściwości molekularne (29). Łączna liczba deskryptorów branych pod uwagę przy budowie modelu deskryptorowego QSAR wynosi zatem 912. Do analizy regresji wybrano PLS, ponieważ liczba instancji jest znacznie mniejsza niż liczba atrybutów (deskryptorów) przy użyciu programu MINITAB. Jak wspomnieliśmy wcześniej, PLS jest szeroko stosowany do opracowywania modeli QSAR poprzez redukcję liczby atrybutów w zbiorze deskryptorów do niewielkiej liczby atrybutów skorelowanych ze zdefiniowaną właściwością podlegającą modelowaniu, którą w naszym badaniu są eksperymentalne wartości IC50.
Budowa modelu z PLS w celu wyboru najbardziej informacyjnych deskryptorów
Głównym celem analizy regresji jest wyznaczenie modelu przewidującego aktywność (IC50) kandydatów na leki pod względem deskryptorów. PLS może być określona jako metoda MLR blisko spokrewniona z regresją składowych głównych. Zasadniczo, przeprowadzając badanie PLS możemy przewidywać zestaw zmiennych zależnych Y na podstawie zestawu zmiennych niezależnych X za pomocą programu MINITAB, który podawał nam przebiegi PLS automatycznie w oparciu o ustalony przez nas górny limit liczby najbardziej znaczących deskryptorów. Każdy przebieg PLS dostarcza modelu liniowego zmiennej zależnej (wartości IC50) w odniesieniu do zmiennych niezależnych (najbardziej znaczących deskryptorów). W tym momencie budowany jest odpowiedni model i wyznaczane są najbardziej znaczące deskryptory. Kolejnym krokiem będzie wstępna klasyfikacja leków na podstawie deskryptorów. Wybór znaczących deskryptorów przez pierwsze przebiegi PLS może nie być najbardziej efektywny w klasyfikacji. Dlatego też przeprowadzamy testy istotności na wybranych deskryptorach za pomocą analizy regresji, aby zwiększyć dokładność klasyfikacji.
Klasyfikacja kandydatów na leki metodą hiper-boxów opartą na MILP
Trzeci krok poświęcony jest klasyfikacji leków; stosujemy metodę hiper-boxów opartą na MILP wykorzystując wybrane deskryptory z poprzedniego kroku.
Celem w problemach klasyfikacji danych jest przypisanie punktów danych, które są opisane pewną liczbą atrybutów, do predefiniowanych klas. The strength of hyper-boxes classification method is from its ability to use more than one hyper-box when defining a class as shown in Figure 3, and this ability prevents overlapping in the classes, which would not be prevented if the classes were defined with a single hyper-box only.
Schematic representation of multi-class data classification using hyper-boxes.
The data classification problem is solved in two steps: training step and testing step. W kroku szkoleniowym granice klas są tworzone poprzez konstrukcję hiper-boxów, natomiast skuteczność skonstruowanych klas jest testowana w kroku testowym.
Problem MILP dla klasyfikacji jest skonstruowany w taki sposób, że funkcją celu jest minimalizacja błędnych klasyfikacji w zbiorze danych przy minimalnej liczbie hiper-boxów w kroku szkoleniowym. Minimalizacja liczby hiper-boxów, czyli eliminacja zbędnego użycia hiper-boxów, jest wymuszona przez karanie istnienia pudełka małym skalarem w funkcji celu. W części szkoleniowej górna i dolna granica każdego hiper-boxu jest również obliczana na podstawie punktów danych zawartych w tym hiper-boxie.
W etapie testowania punkty danych są przypisywane do klas poprzez obliczenie odległości pomiędzy punktem danych a każdym z pudełek i określenie pudełka, które jest najbliżej punktu danych. W końcu, oryginalne i przypisane klasy testowych punktów danych są porównywane, a skuteczność klasyfikacji jest uzyskiwana poprzez poprawnie sklasyfikowane instancje.
Rozwiązanie zaproponowanego problemu MILP do optymalności jest obliczeniowo trudne dla dużych zbiorów danych ze względu na dużą liczbę zmiennych binarnych. W związku z tym, opracowano trójstopniową metodę dekompozycji do uzyskiwania optymalnych rozwiązań problemów klasyfikacji dużych danych. W pierwszym etapie identyfikowane są przypadki trudne do sklasyfikowania, które określamy mianem przetwarzania wstępnego. Ponadto, dla każdej klasy wyznaczane są nasiona w celu poprawy efektywności obliczeniowej. Przy większym nacisku na te obserwacje, w drugim etapie otrzymywane jest rozwiązanie problemu przy użyciu zmodyfikowanego modelu. Wreszcie, w trzecim etapie przeprowadzane są ostateczne przypisania i eliminacje przecięć.
W niniejszej pracy zastosowaliśmy opisaną powyżej metodę do klasyfikacji aktywności molekuł leków dla rozważanych zbiorów danych. Przeprowadzamy 10-krotną walidację krzyżową przy wyborze zbiorów treningowych i testowych, gdzie losowo dzielimy zbiory danych na 10 podprób o równej liczbie członków. Z tych 10 podpróbek 9 z nich jest łączonych i używanych jako zbiór treningowy, a pozostała 1 podpróba jest używana jako zbiór testowy. Następnie klasyfikacja jest wykonywana 10 razy z każdą z 10 podpróbek użytych dokładnie raz jako zbiór testowy. Na koniec, dokładność klasyfikacji jest raportowana jako średnia z tych 10 klasyfikacji.
Klasyfikujemy każdego z kandydatów na leki w zbiorze testowym jako posiadającego niską lub wysoką wartość IC50. W tym iteracyjnym badaniu, ten etap klasyfikacji jest wykonywany kilkakrotnie: najpierw z początkowym zestawem deskryptorów, a następnie przy użyciu ulepszonego zestawu deskryptorów uzyskanych z analizy istotności.
Analiza istotności
W czwartym kroku wykonywane są testy istotności. Po przeprowadzeniu PLS możliwe jest uznanie deskryptora za istotny, podczas gdy w rzeczywistości taki nie jest i problem ten jest rozwiązywany poprzez przeprowadzenie testów istotności po klasyfikacji pierwotnej. Główna idea testu istotności jest następująca: Jeśli Z jest całym zbiorem kandydatów na leki, załóżmy, że po klasyfikacji jest on podzielony na dwie klasy, A i B. Aby klasyfikacja była udana, wariancja wartości deskryptorów powinna być mniejsza w obrębie klas A i B niż jest dla całej populacji, Z.
Równanie podane poniżej w Eq. 2.1 przedstawia rozkład F.
gdzie, S i j 2 MathType@MTEF@5@5@+=feaagaart1ev2aqatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaem4uam1aa0baaSqaaiabdMgaPjabdQgaQbqaaiabikdaYaaaaaa@30DC@ jest wariancją próby wartości dla deskryptora i dla zestawu leków j, ν = n-1 i η = m-1 są stopniami swobody, n jest liczbą wartości deskryptora i dla zbioru leków j, a m jest liczbą wartości deskryptora i dla zbioru leków k.
Testowanie hipotezy odbywa się za pomocą hipotezy zerowej S i j 2 = S i k 2 MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaem4uam1aa0baaSqaaiabdMgaPjabdQgaQbqaaiabikdaYaaakiabg2da9iabdofatnaaDaaaleaacqWGPbqAcqWGRbWAaeaacqaIYaGmaaaaaa@36F3@ , co sugeruje, że wariancja całego zbioru kandydatów na leki jest równa wariancji leków z tej samej klasy. Ponieważ wariancja całego zbioru leków powinna być większa niż wariancja w obrębie klasy, definiujemy naszą hipotezę alternatywną jako: H a = S i j 2 ≻ S i k 2 MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemisaG0aaSbaaSqaaiabdggaHbqabaGccqGH9aqpcqWGtbWudaqhaaWcbaGaemyAaKMaemOAaOgabaGaeGOmaidaaOGaeS4EIyMaem4uam1aa0baaSqaaiabdMgaPjabdUgaRbqaaiabikdaYaaaaaa@3B21@ , gdzie j jest członkiem całego zbioru danych, a k jest członkiem klasy. Zauważmy, że wartość p fvη w bieżącym modelu powinna być mniejsza niż wartość p fvη w poprzednim modelu, aby przyjąć hipotezę alternatywną.
Budujemy nowy model klasyfikacyjny
Ten ostatni krok jest wykonywany, gdy stwierdzimy, że istnieją przeszacowane deskryptory w modelu podczas kroku czwartego.
W związku z tym, całkowita liczba 3 modeli jest konstruowana poprzez analizę regresji, wybierając odpowiednio 7, 10 i 15 deskryptorów jako zmienne reprezentatywne dla każdego modelu, a analiza istotności jest stosowana do wszystkich deskryptorów w tych 3 modelach. Jeśli stwierdzimy istnienie nieistotnej zmiennej w jednym z tych modeli, zastępujemy je tymi, które są istotne w pozostałych modelach. Udowodniono, że taka korekta poprawia naszą dokładność klasyfikacji. Kiedy zastępujemy mniej znaczące zmienne, pozostałe 880 deskryptorów, które zostały wyeliminowane podczas analizy PLS, są ignorowane, ponieważ te 7, 10 i 15 atrybutów zostało wybranych przez analizę regresji PLS i mają udowodnioną siłę w opisywaniu wartości IC50. Głównym celem analizy regresji PLS w rzeczywistości jest wyeliminowanie statystycznie nieistotnych cech i dostarczenie nam najbardziej znaczącej przestrzeni próbki do dalszej pracy.
Wyniki uzyskane przez naszą metodę są porównywane ze wszystkimi 63 metodami klasyfikacji dostępnymi w WEKA, a 16 najlepszych klasyfikatorów WEKA raportuje wyniki uzyskane przez nasz algorytm w Tabeli 3, wraz z odpowiadającą im dokładnością klasyfikacji. Atrybuty użyte w klasyfikatorach WEKA są tymi samymi deskryptorami, które zostały znalezione po testach istotności, a 10-krotna walidacja krzyżowa została zastosowana do każdego klasyfikatora, w tym naszej metody klasyfikacji.
WEKA jest potężnym narzędziem eksploracji danych do wykorzystania w celach porównawczych, ponieważ zawiera wszystkie powszechnie znane algorytmy uczenia maszynowego wśród swoich 63 klasyfikatorów. Sukces tych istniejących algorytmów uczenia maszynowego w binarnej klasyfikacji aktywnych i nieaktywnych związków na podstawie ich wartości deskryptorów został również wcześniej opisany. Poniżej znajduje się krótki przegląd najlepiej działających metod klasyfikacji danych dostępnych w WEKA. Sieć bayesowskaB = <N, A, Φ > jest skierowanym grafem acyklicznym <N, A> z warunkowym rozkładem prawdopodobieństwa dołączonym do każdego węzła, reprezentowanym zbiorczo przez Φ. Każdy węzeł n ∈ N reprezentuje atrybut zbioru danych, a każdy łuk a ∈ A pomiędzy węzłami reprezentuje zależność probabilistyczną. Klasyfikator Naive Bayes zakłada, że wszystkie zmienne są od siebie niezależne, gdzie węzeł klasyfikacyjny jest reprezentowany jako węzeł nadrzędny wszystkich innych węzłów. Naive Bayes Simple używa rozkładu normalnego do modelowania atrybutów i obsługuje atrybuty numeryczne używając dyskretyzacji nadzorczej, podczas gdy Naive Bayes Updateable jest wersją przyrostową, która przetwarza jedną instancję na raz i używa estymatora jądra zamiast dyskretyzacji.
Klasyfikator logistyczny buduje dwuklasowy model regresji logistycznej. Jest to model regresji statystycznej, gdzie regresja logistyczna zakłada, że stosunek logarytmów prawdopodobieństwa rozkładów klas jest liniowy w obserwacjach. Prosty klasyfikator logistyczny buduje liniowe modele regresji logistycznej w oparciu o jeden atrybut. Model ten jest uogólnionym modelem zwykłego modelu regresji najmniejszych kwadratów. Perceptron wielowarstwowy jest siecią neuronową, która wykorzystuje propagację wsteczną. Perceptron, który jest elementem przetwarzającym, oblicza pojedyncze wyjście, nieliniową funkcję aktywacji liniowej kombinacji wielu wejść, której parametry są uczone w fazie treningu. SMO (sequential minimal optimization), zwana również WEKA SVM (support vector machine), jest metodą trenowania klasyfikatora wektora wsparcia przy użyciu wielomianowych jąder poprzez rozbicie dużego problemu optymalizacji programowania kwadratowego na mniejsze problemy optymalizacji QP.
IB1 jest wymieniony jako leniwy klasyfikator, w tym sensie, że przechowuje instancje treningowe i nie wykonuje żadnej pracy aż do czasu klasyfikacji. IB1 jest uczniem opartym na instancji. Znajduje instancję treningową najbliższą w odległości euklidesowej do danej instancji testowej. IBk to klasyfikator k-najbliższych sąsiadów, który wykorzystuje ten sam pomysł.
Logit Boost wykorzystuje addytywną regresję logistyczną. Algorytm może zostać przyspieszony poprzez przypisanie wagom określonego progu. Multi Class Classifier wykorzystuje cztery różne metody klasyfikacji dwuklasowej do rozwiązywania problemów wieloklasowych. Threshold Selector, który jest meta-uczniem, optymalizuje miarę F poprzez wybór progu prawdopodobieństwa na wyjściu klasyfikatora.
Random forest i LMT są metodami drzew decyzyjnych. Random Forest generuje drzewa losowe poprzez zbieranie zespołów drzew losowych, podczas gdy LMT buduje drzewa modelu logistycznego i używa walidacji krzyżowej do określenia liczby iteracji podczas dopasowywania funkcji regresji logistycznej w każdym węźle. OneR (jedna reguła) buduje jednopoziomowe drzewo decyzyjne i uczy się reguły na podstawie każdego atrybutu i wybiera regułę o najmniejszym współczynniku błędu jako jedyną regułę.