- algorytm klasyfikacji sekwencji
- tworzenie bazy danych
- struktura bazy danych i algorytm wyszukiwania
- konstruowanie symulowanych metagenomów
- ocena dokładności i szybkości
- zmniejszono rozmiar bazy danych
- wykorzystanie genomów projektowych
- eksperymenty wykluczające Klad
- klasyfikacja danych mikrobiomu ludzkiego
- dostępność oprogramowania i danych
algorytm klasyfikacji sekwencji
aby sklasyfikować sekwencję s DNA, zbieramy wszystkie k-Mery w tej sekwencji w zestaw, oznaczony jako K(S). Następnie mapujemy każdy k-mer w K (S), używając algorytmu opisanego poniżej, do taksonu LCA wszystkich genomów, które zawierają k-mer. Te taksony LCA i ich przodkowie w drzewie taksonomicznym tworzą to, co nazywamy drzewem klasyfikacyjnym, przycinanym podzbiorem używanym do klasyfikacji S. każdy węzeł w drzewie klasyfikacyjnym jest ważony liczbą k-merów w K(S), które odwzorowały takson związany z tym węzłem. Następnie każda ścieżka od korzenia do liścia (RTL) w drzewie klasyfikacji jest punktowana przez obliczenie sumy wszystkich wag węzłów wzdłuż ścieżki. Maksymalna ścieżka punktacji RTL w drzewie klasyfikacji jest ścieżką klasyfikacyjną, A S jest przypisana Etykieta odpowiadająca jego liściu (jeśli istnieje wiele ścieżek o maksymalnej punktacji, wybierany jest LCA wszystkich liści tych ścieżek). Algorytm ten, zilustrowany na rysunku 1, pozwala Krakenowi uznać każdy k-mer w sekwencji za oddzielny dowód, a następnie próbować rozwiązać wszelkie sprzeczne dowody, jeśli to konieczne. Zauważ, że dla odpowiedniego wyboru k, większość k-Mer będzie mapować unikalnie do jednego gatunku, znacznie upraszczając proces klasyfikacji. Sekwencje, dla których żaden z k-merów w K(S) nie znajduje się w żadnym genomie, nie są klasyfikowane przez ten algorytm.
zastosowanie punktacji ścieżki RTL w drzewie klasyfikacji jest konieczne w świetle nieuniknionych różnic między sekwencjami, które mają być sklasyfikowane, a sekwencjami obecnymi w dowolnej bibliotece genomów. Takie różnice mogą, nawet dla dużych wartości k, skutkować K-mer, który jest obecny w bibliotece, ale związany z gatunkiem dalekim od prawdziwego gatunku źródłowego. Punktując różne ścieżki RTL w drzewie klasyfikacji, możemy zrekompensować te różnice i poprawnie klasyfikować sekwencje, nawet jeśli niewielka mniejszość k-merów w sekwencji wskazuje, że sekwencji należy przypisać nieprawidłową Etykietę taksonomiczną.
tworzenie bazy danych
wydajna implementacja algorytmu klasyfikacji Krakena wymaga, aby mapowanie k-mers do taksonów było wykonywane przez odpytywanie wstępnie obliczonej bazy danych. Kraken tworzy tę bazę danych w wieloetapowym procesie, rozpoczynając od wyboru biblioteki sekwencji genomowych. Kraken zawiera domyślną bibliotekę opartą na ukończonych genomach drobnoustrojów w bazie danych RefSeq Narodowego Centrum Informacji biotechnologicznej (NCBI), ale biblioteka może być dostosowana do potrzeb indywidualnych użytkowników .
Po wybraniu biblioteki używamy wielowątkowego licznika K-mer Jellyfish, aby utworzyć bazę danych zawierającą wszystkie odrębne 31-mer w bibliotece. Po ukończeniu bazy danych, 4-bajtowe przestrzenie używane do przechowywania liczników k-mer w pliku bazy danych są używane przez Krakena do przechowywania taksonomicznych numerów IDENTYFIKACYJNYCH wartości LCA k-mers. Po utworzeniu bazy danych przez meduzę sekwencje genomowe w bibliotece są przetwarzane pojedynczo. Dla każdej sekwencji takson z nią związany służy do ustawiania zapisanych wartości LCA wszystkich k-merów w sekwencji. Gdy sekwencje są przetwarzane, jeśli k-mer z sekwencji miał wcześniej ustawioną wartość LCA, to LCA przechowywanej wartości i taksonu bieżącej sekwencji jest obliczany i LCA jest przechowywany dla k-mer. Informacje o taksonach pochodzą z bazy taksonomicznej NCBI.
struktura bazy danych i algorytm wyszukiwania
ponieważ Kraken bardzo często używa k-mer jako zapytania do bazy danych natychmiast po zapytaniu sąsiednich K-mer, a ponieważ sąsiednie k-Mery mają znaczną ilość sekwencji, wykorzystujemy koncepcję minimizera, aby zgrupować podobne k-Mery razem. Aby wyjaśnić nasze zastosowanie tego pojęcia, definiujemy tutaj kanoniczną reprezentację sekwencji DNA S jako leksykograficznie mniejszą od S i odwrotną dopełniacz S. Aby określić minimalizator K-mer długości M, rozważamy kanoniczną reprezentację wszystkich M-merów w k-mer i wybieramy najmniejszy leksykograficznie z tych M-Mer jako minimalizator k-mer. W praktyce sąsiadujące k-Mery często mają ten sam minimizer.
w bazie danych Krakena wszystkie k-Mery z tym samym minimizerem są przechowywane kolejno i sortowane w porządku leksykograficznym ich kanonicznych reprezentacji. Zapytanie o k-mer R może być następnie przetwarzane przez wyszukanie w indeksie pozycji w bazie danych, w których przechowywane byłyby K-Mery z minimalizerem R, a następnie przeprowadzenie wyszukiwania binarnego w tym regionie (Rysunek 5). Ponieważ sąsiednie k-Mery często mają ten sam minimizer, zakres wyszukiwania jest często taki sam między dwoma kolejnymi zapytaniami, a wyszukiwanie w pierwszym zapytaniu często może przynieść dane do pamięci podręcznej procesora, które będą używane w drugim zapytaniu. Pozwalając dostępom do pamięci w kolejnych zapytaniach na dostęp do danych w pamięci podręcznej procesora zamiast pamięci RAM, strategia ta sprawia, że kolejne zapytania są znacznie szybsze niż w innym przypadku. Indeks zawierający przesunięcia każdej grupy k-M w bazie danych wymaga 8 × 4M bajtów. Domyślnie Kraken używa minimizerów 15-bp, ale użytkownik może modyfikować tę wartość; na przykład podczas tworzenia MiniKraken użyliśmy minimizerów 13-bp, aby zapewnić całkowity rozmiar bazy danych poniżej 4 GB.
struktura bazy danych Kraken. Każdy K-mer, który ma zostać zapytany o bazę danych, ma określony podłańcuch, który jest jego minimalizatorem. Aby wyszukać k-mer w bazie danych, bada się pozycje w bazie danych, które zawierają k-Mer z tym samym minimalizatorem. Pozycje te są szybko znajdowane przez zbadanie tablicy przesunięć minimalizatora dla początkowych pozycji rekordów za pomocą minimalizatora K-mer (pomarańczowy) i następnego możliwego minimalizatora (niebieski). W zakresie rekordów związanych z danym minimizerem rekordy są sortowane według leksykograficznego porządku ich k-merów, co pozwala na wypełnienie zapytania za pomocą wyszukiwania binarnego w tym zakresie.
wdrażając Krakena, dokonaliśmy dalszych optymalizacji struktury i algorytmu wyszukiwania opisanego powyżej. Po pierwsze, jak zauważył Roberts et al. , proste leksykograficzne porządkowanie M-merów może skutkować wypaczonym rozkładem minimalizatorów, które nad-reprezentują M-Mery o niskiej złożoności. W Krakenie takie odchylenie stworzyłoby wiele dużych zakresów wyszukiwania, co wymagałoby więcej czasu na wyszukiwanie. Aby stworzyć bardziej równomierną dystrybucję minimalizatorów (a tym samym przyspieszyć wyszukiwanie), używamy operacji exclusive-or (XOR) do przełączania połowy bitów kanonicznej reprezentacji każdego m-Mera przed porównaniem M-merów ze sobą za pomocą uporządkowania leksykograficznego. Ta operacja XOR skutecznie koduje standardowe zamówienia i zapobiega dużemu nastawieniu do minimalizatorów o niskiej złożoności.
korzystamy również z faktu, że zakres wyszukiwania jest często taki sam między zapytaniami, aby zapytania Krakena były szybsze. Zamiast obliczać minimalizator za każdym razem, gdy wykonujemy zapytanie, najpierw przeszukujemy poprzedni zakres. Jeśli nasz zapytany k-mer znajdzie się w tym zakresie, zapytanie może natychmiast powrócić. Jeśli k-mer nie zostanie znaleziony, to minimizer jest obliczany; jeśli minimizer K-mer jest taki sam jak ostatni zapytany k-mer, to zapytanie nie powiedzie się, ponieważ przestrzeń wyszukiwania minimizera została pokazana, że nie ma k-mer. Tylko jeśli minimizer się zmienił, Kraken musi dostosować zakres wyszukiwania i ponownie wyszukać k-mer.
konstruowanie symulowanych metagenomów
metagenomy HiSeq i MiSeq zostały zbudowane przy użyciu 20 zestawów bakteryjnych odczytów z całego genomu. Odczyty te zostały znalezione w ramach projektu GAGE-B lub w archiwum odczytu sekwencji NCBI. Każdy metagenom zawiera sekwencje z dziesięciu genomów (dodatkowy plik 1: Tabela S1). Zarówno dla 10 000, jak i 10 milionów odczytywanych próbek każdego z tych metagenomów, 10% ich sekwencji wybrano z każdego z dziesięciu zestawów danych genomu (tj. każdy Genom miał taką samą obfitość sekwencji). Wszystkie sekwencje zostały przycięte w celu usunięcia niskiej jakości baz i sekwencji adapterów.
skład tych dwóch metagenomów stwarza pewne wyzwania dla naszych klasyfikatorów. Na przykład Pelosinus fermentans, znaleziony w naszym metagenomie HiSeq, nie może być poprawnie zidentyfikowany na poziomie rodzaju przez Krakena (lub którykolwiek z innych wcześniej opisanych klasyfikatorów), ponieważ w bazie RefSeq complete genomes database nie ma genomów Pelosinus; jednak w bazie Kraken-GB jest siedem takich genomów, w tym sześć szczepów P. fermentans. Podobnie, w metagenomie MiSeq, Proteus vulgaris jest często błędnie klasyfikowany na poziomie rodzaju, ponieważ jedynym genomem Proteusa w bazie danych Krakena jest pojedynczy Genom Proteus mirabilis. W bazie danych Kraken-GB znajduje się jeszcze pięć genomów Proteusa, co pozwala Kraken-GB lepiej sklasyfikować odczyty z tego rodzaju. Ponadto Metagenom MiSeq zawiera pięć genomów z rodziny Enterobacteriaceae (Citrobacter, Enterobacter, Klebsiella, Proteus i Salmonella). Wysokie podobieństwo sekwencji między rodzajami w tej rodzinie może utrudnić rozróżnianie rodzajów dla każdego klasyfikatora.
metagenom simBA-5 został stworzony przez symulowanie odczytów z zestawu kompletnych genomów bakteryjnych i archaealnych w RefSeq. Replikony z tych genomów były używane, jeśli były związane z taksonem, który miał wpis związany z rangą rodzaju, co skutkowało zestawem replikonów z 607 rodzajów. Następnie użyliśmy symulatora odczytu Masona z jego modelem Illumina do wytworzenia 10 milionów odczytów 100 bp z tych genomów. Najpierw stworzyliśmy symulowane genomy dla każdego gatunku, wykorzystując wskaźnik SNP wynoszący 0,1% i wskaźnik indel wynoszący 0,1% (oba parametry domyślne), z których wygenerowaliśmy odczyty. W przypadku symulowanych odczytów pomnożyliśmy domyślne współczynniki niedopasowania i indel o pięć, co skutkuje średnim współczynnikiem niedopasowania wynoszącym 2% (od 1% na początku odczytu do 6% na końcu) i wskaźnikiem indel wynoszącym 1% (0,5% prawdopodobieństwo wprowadzenia i 0,5% prawdopodobieństwo usunięcia). Dla metagenomu simBA-5 zestaw odczytu 10,000 został wygenerowany z losowej próbki zestawu odczytu 10 milionów.
ocena dokładności i szybkości
zdecydowaliśmy się zmierzyć dokładność przede wszystkim na poziomie rodzaju, który był najniższym poziomem, dla którego mogliśmy łatwo określić informacje taksonomiczne dla prognoz PhymmBL i NBC w sposób zautomatyzowany. (Wynika to ze sposobu, w jaki PhymmBL i NBC zgłaszają swoje wyniki). Ponieważ niektóre genomy nie mają pozycji taksonomicznych we wszystkich siedmiu Rang (gatunek, Rodzaj, Rodzina, kolejność, Klasa, filum i Królestwo), zdefiniowaliśmy wrażliwość na poziomie rodzaju jako A/B, gdzie A jest liczbą odczytów z przypisanym rodzajem, które zostały poprawnie sklasyfikowane w tej randze, A B jest całkowitą liczbą odczytów z dowolnym przypisanym rodzajem. Podobnie zdefiniowaliśmy czułość dla innych Rang taksonomicznych.
ponieważ Kraken może klasyfikować odczyt na poziomie powyżej gatunku, pomiar jego precyzji wymaga od nas określenia wpływu na precyzję przypisania właściwego rodzaju (na przykład), a nie przypisania gatunku w ogóle. Z tego powodu, zdefiniowaliśmy dokładność na poziomie rangi jako C / (D + E), gdzie C jest liczbą odczytów oznaczonych na lub poniżej WŁAŚCIWEGO taksonu w mierzonej randze, D jest liczbą odczytów oznaczonych na lub poniżej mierzonej rangi, A E jest liczbą odczytów nieprawidłowo oznaczonych powyżej mierzonej rangi. Na przykład, biorąc pod uwagę odczyt R, który powinien być oznaczony jako Escherichia coli, etykietowanie R jako E. coli, E. fergusonii lub Escherichia poprawiłoby precyzję na poziomie rodzaju. Etykieta Enterobacteriaceae (prawidłowa rodzina) lub Proteobacteria (prawidłowy filum) nie ma wpływu na dokładność poziomu rodzaju. Oznaczenie dla R Bacillus (nieprawidłowy rodzaj) lub Firmicutes (nieprawidłowy rodzaj) zmniejszyłoby dokładność poziomu rodzaju.
oceniając dokładność Phymmbla, zgodnie z poradami jego twórców, wybraliśmy próg pewności rodzaju dla naszych porównań. Z zestawu danych simmc (simulated medium complexity) wybraliśmy 3333 odczyty, obejmujące 31 różnych rodzajów. Aby symulować krótkie odczyty z danych sekwencji Sangera w zestawie simMC, wybraliśmy ostatnie 100 bp z każdego odczytu. Następnie porównaliśmy phymmbl z odczytami 100 bp i oceniliśmy czułość i precyzję prognoz PhymmBL z progami ufności rodzaju od 0 do 1, w krokach co 0,05. Okazało się, że próg 0,65 dawał najwyższy wynik F (harmoniczną średnią czułości i precyzji), przy czym 0,60 i 0,70 również miały wyniki F w granicach 0.5 punktów procentowych maksimum (dodatkowy plik 1: Tabela S2). Dlatego w naszych porównaniach wykorzystaliśmy próg zaufania 0,65%. Chociaż wybór progu zależy od indywidualnych potrzeb użytkownika, a więc jest do pewnego stopnia arbitralny, próg wybrany w ten sposób zapewnia bardziej właściwe porównanie do selektywnego klasyfikatora, takiego jak Kraken, niż żaden próg w ogóle.
wyniki czasu i dokładności podczas używania Megablast jako klasyfikatora uzyskano z danych dziennika wyprodukowanych przez PhymmBL, ponieważ PhymmBL używa Megablast do swojego kroku wyrównania. Przypisując Etykietę taksonomiczną do odczytu za pomocą Megablast, użyliśmy taksonu związanego z pierwszym zgłoszonym wyrównaniem. Megablast został uruchomiony z domyślnymi opcjami.
szybkość była oceniana przy użyciu jednowątkowego działania każdego programu (z wyjątkiem NBC). PhymmBL został zmieniony tak, że jego wywołanie do programu blastn używało jednego wątku zamiast dwóch. NBC zostało uruchomione z 36 jednoczesnymi procesami działającymi na disjoint zestawach genomów w swojej bibliotece genomowej, a całkowity czas klasyfikacji został określony przez zsumowanie czasów dekompresji i punktacji dla każdego genomu. Zegary ścienne były rejestrowane dla wszystkich klasyfikatorów. Porównując Krakena z innymi klasyfikatorami, użyliśmy BLAST+ 2.2.27, PhymmBL 4.0, NBC 1.1 i Metaflan 1.7.6. Wszystkie zostały uruchomione na tym samym komputerze, z 48 procesorami AMD Opteron 6172 2,1 GHz i 252 GB PAMIĘCI RAM, z systemem Red Hat Enterprise Linux 5. Zestawy danych używane do oceny prędkości miały po 10 000 odczytów dla wszystkich programów innych niż Kraken (i jego warianty) i Metaflan, które wykorzystywały 10 000 000 odczytów zestawów danych. Wyższe liczby odczytu były używane z tymi szybszymi programami, aby zminimalizować efekt początkowych i końcowych operacji, które mają miejsce podczas wykonywania programów.
chociaż Kraken jest jedynym programem, który wyraźnie wykonuje operacje, aby upewnić się, że jego dane znajdują się w pamięci fizycznej przed klasyfikacją, chcieliśmy mieć pewność, że wszystkie programy zostały ocenione w podobny sposób. Oceniając szybkość, dla każdego programu odczytujemy wszystkie pliki bazy danych (np. Pliki IMM i bazy danych BLAST Dla PhymmBL, listy częstotliwości K-mer dla NBC i indeks Bowtie dla Metaflan) do pamięci trzy razy przed uruchomieniem programu, w celu umieszczenia zawartości bazy danych w pamięci podręcznej systemu operacyjnego (która jest przechowywana w pamięci fizycznej).
zmniejszono rozmiar bazy danych
aby wygenerować bazę danych o pojemności 4 GB dla naszych wyników MiniKraken, usunęliśmy pierwsze 18 z każdego bloku 19 rekordów w standardowej bazie danych Kraken. Wybrano współczynnik zmniejszania wynoszący 19, ponieważ był to najmniejszy współczynnik całkowity, który zmniejszyłby rozmiar do mniej niż 4 GB, rozmiar, który można łatwo zmieścić w pamięci wielu popularnych komputerów osobistych. Dla użytkowników, którzy mają więcej dostępnej pamięci RAM, Kraken pozwala na użycie mniejszego współczynnika zmniejszania, co daje zwiększoną czułość.
wykorzystanie genomów projektowych
podczas konstruowania bazy danych Kraken-GB zauważyliśmy, że na końcach znajduje się kilka styków ze znanymi sekwencjami adapterów. W późniejszych testach odkryliśmy również, że niektóre sekwencje w próbkach z dużą ilością sekwencji ludzkiej były konsekwentnie błędnie klasyfikowane przez tę bazę danych, co doprowadziło nas do wniosku, że zanieczyszczenie prawdopodobnie było obecne w genomach projektowych. Próbując przeciwdziałać temu skażeniu, usunęliśmy z bazy danych te k – Mery ze znanych sekwencji adapterów, a także pierwsze i ostatnie 20 k-Mery z każdego z ciągów projektowych. Chociaż poprawiło to klasyfikację, nie wyeliminowało problemu błędnej klasyfikacji. Z tego powodu uważamy, że jeśli genomy projektowe są używane w bazie danych Kraken, należy zastosować bardzo rygorystyczne środki w celu usunięcia sekwencji zanieczyszczeń z biblioteki genomowej.
eksperymenty wykluczające Klad
podczas ponownej analizy zestawu danych simBA-5 dla naszych eksperymentów wykluczających Klad, niektóre odczyty nie były używane dla niektórych par mierzonych i wykluczonych Rang. Jeśli pochodzenie read ’ a nie posiadało wpisu taksonomicznego w jednej z mierzonych lub wykluczonych Rang, nie było używane do tego konkretnego eksperymentu.
ponadto, odczyt nie został użyty w eksperymencie, chyba że co najmniej dwa inne taksony reprezentowane w naszej bazie danych (oprócz wykluczonego kladu) w randze wykluczonej podzieliły takson kladu pochodzenia w randze mierzonej. Na przykład, odczyt z rodzaju G nie byłby użyty w eksperymencie mierzącym dokładność w randze klasy i wykluczającym rangę rodzaju, chyba że Klasa domowa G miała co najmniej dwa inne rodzaje z genomami w bibliotece genomowej Krakena. Bez tego etapu filtrowania, gdyby rodzaj był wykluczony, gdy był jedynym rodzajem w swojej klasie, Kraken nie mógłby nazwać właściwej klasy, ponieważ wszystkie wpisy w bazie danych z tej klasy również zostałyby wykluczone. Jest to to samo podejście podjęte w podobnych eksperymentach, które zostały wykorzystane do oceny PhymmBL .
klasyfikacja danych mikrobiomu ludzkiego
sklasyfikowaliśmy dane Human Microbiome Project przy użyciu bazy danych Kraken wykonanej z kompletnych genomów bakteryjnych, archaealnych i wirusowych RefSeq, wraz z ludzkim genomem GRCh37. Pobraliśmy sekwencje trzech akcesji (SRS019120, SRS014468 i SRS015055) z archiwum odczytu sekwencji NCBI, a każda akcesja miała dwa biegi. Wszystkie odczyty zostały przycięte, aby usunąć niskiej jakości bazy i sekwencje adapterów. Korona służyła do generowania wszystkich działek podziału taksonomicznego.
ponieważ wszystkie sekwencje były sparowanymi odczytami, połączyliśmy je ze sobą, łącząc pary z sekwencją 'NNNNN’ między nimi. Kraken ignoruje k – Mery z niejednoznacznymi nukleotydami, więc K-Mery, które obejmują te ” N ” znaki, nie wpływają na klasyfikację. Ta operacja pozwoliła Krakenowi sklasyfikować parę odczytów jako pojedynczą jednostkę, zamiast klasyfikować partnerów osobno.
