Optymalizacja podkładów PCR ukierunkowanych na bakteryjne 16S ribosomal RNA Gen / BMC Bioinformatyka

ograniczenia problemowe

jak wspomniano w poprzednim akapicie, optymalna para primer-set-powinna jednocześnie maksymalizować wydajność i zasięg oraz minimalizować dopasowanie-bias. Poniżej opiszemy jak ilościowo zakodowaliśmy te ograniczenia.

wydajność

idealne pary primer-set-powinny spełniać kilka ograniczeń, mających na celu poprawę wydajności i specyficzności PCR ., Jednak jednoczesne spełnianie wszystkich ograniczeń jest często niepraktyczne i większość najnowocześniejszych podkładów narusza jedno lub więcej ograniczeń . Dlatego postanowiliśmy wprowadzić efektywność jako wynik optymalizacji, kodując wiele ograniczeń jako funkcje oceny rozmytej. Dokładniej, nasz wynik efektywności zdefiniowaliśmy jako sumę dziesięciu terminów: siedem terminów rozmytych związanych z ograniczeniami wydajności pojedynczego podkładu, uśrednionych dla wszystkich podkładów w parach primer-set-oraz trzy terminy dotyczące wydajności PAR primer-set – jako całości., Ponieważ wszystkie terminy mają się różnić od 0 do 1, wynik optymalizacji waha się od 0 (minimalna wydajność) do 10 (maksymalna wydajność).

Ogólnie rzecz biorąc, nasz wynik rozmyty liczy 1 dla każdego ograniczenia, które jest idealnie spełnione, lub, alternatywnie, wartość między 0 a 1 w zależności od tego, jak blisko jest podkład do limitu ograniczeń. Jako przykład rozważmy temperaturę topnienia podkładu, Tm. Tm powinien być większy lub równy 52 stopni w idealnym podkładu, ale 51 jest nadal tolerowalny, choć nie idealny., W tym przypadku nasza funkcja fuzzy scoring przypisuje 1 do temperatur 52 stopni lub większych, 0 do temperatur 50 stopni lub mniejszych i uwzględnia liniową funkcję zwiększania między 50 a 52 stopniami. Każdy termin jest dokładnie opisany w tym, co następuje.

7 terminów oceny pojedynczego podkładu to:

temperatura topnienia: Temperatura topnienia TM podkładu oblicza się ze wzoru najbliższego sąsiada . Termin wynik jest 1 jeśli TM ≥ 52, 0 jeśli TM ≤ 50 i (Tm-50) / 2 Jeśli 50 < Tm < 52.,
GC-content: GC-content jest ułamkiem fGC par zasad w sekwencji startera równym G (guanina) lub C (cytozyna). Termin wynik wynosi 1, jeśli 0, 5 ≤ fGC ≤ 0, 7, 0, jeśli fGC > 0, 7 lub fGC < 0, 4 i (0, 5 – fGC)/0, 1, jeśli 0, 4 ≤ fGC < 0, 5.
3′-end stability – score term 1: dwa terminy score są zdefiniowane dla 3 ' – end stability. Pierwszym terminem jest 0, jeśli trzy ostatnie podstawy podkładu składają się w całości z As (adeniny) i Ts, (tyminy) oraz 1 w przeciwnym razie.,
3 ' – end stability-score term 2: drugi termin score wynosi 0, jeśli ostatnie 5 baz zawiera więcej niż 3 Cs lub Gs, a 1 w przeciwnym razie.
Homopolimery: homopolimer jest sekwencją identycznych nukleotydów. Termin score wynosi 1, jeśli nie ma homopolimerów dłuższych niż 4 nt, 0,5, jeśli nie ma homopolimerów dłuższych niż 5 nt i 0, jeśli w sekwencji znajduje się co najmniej homopolimer dłuższy niż 5 nt.
SELF-dimers: obecność regionów samo-dopełniających się między parami identycznych primerów może prowadzić do generowania SELF-dimers., Biorąc pod uwagę maksymalną liczbę meczów w wyrównaniu bez przerw między podkładem z jego odwrotnym uzupełnieniem, maxM, termin wynik wynosi 1, jeśli maxM ≤ 8, 0, jeśli MAXM ≥ 11 i (11 – maxM)/3, Jeśli 8 < maxM < 11.
spinki do włosów: spinka może być utworzona w obecności samo-komplementarności w sekwencji podkładu, szczególnie na jej 3′-końcu., Termin score wynosi 0, jeśli dla co najmniej jednego wyrównania bez przerwy między podkładem a odwrotnym dopełnieniem jego 3′-końca, zarówno ostatni nukleotyd, jak i 3 lub więcej z 4 poprzedzających nukleotydów pasują do siebie, a 1 w przeciwnym razie.

3 terminy oceny par Grunt-zestaw są zdefiniowane w następujący sposób:

Zakres temperatury topnienia: Zakres temperatury topnienia ΔTm pary Grunt-zestaw jest obliczany jako maksimum minus minimum temperatur topnienia wszystkich gruntujących w zestawie pary., Termin wynik jest 1 jeśli ΔTm ≤ 3, 0 jeśli ΔTm ≥ 5 i (5 – ΔTm)/2 Jeśli 3 < ΔTm < 5.
Dimery: bierzemy pod uwagę maksymalną liczbę dopasowań maxM we wszystkich możliwych wyrównaniach między wszystkimi możliwymi kombinacjami podkładów do przodu i do tyłu z pary primer-set -. Termin wynik jest 1 jeśli maxM ≤ 8, 0 jeśli maxM ≥ 11 i (11-maxM) / 3 Jeśli 8 <maxM< 11.,
Zakres długości Amplikonów: ze względu na znane zmniejszenie wydajności PCR ze wzrostem długości amplikonów, chcemy, aby długości generowanych amplikonów znajdowały się w wąskim zakresie. Szczególnie chcemy uniknąć amplikonów znacznie krótszych od docelowej długości, ponieważ byłyby nadmiernie wzmocnione w stosunku do pozostałych. Chcemy jednak być w stanie tolerować niewielką część wartości odstających, aby uniknąć karania potencjalnie cennych par primer-set-z powodu kilku rzadkich sekwencji., Biorąc pod uwagę reprezentatywny zbiór sekwencji bakteryjnych 16S, zwany od teraz „zestawem odniesienia”, rozważamy różnicę Δamplen między medianą i pierwszym percentylem długości amplikonów we wszystkich możliwych amplikonach, utworzoną przez dopasowanie wszystkich kombinacji początkowych i wstecznych z pary zestawu z zestawem odniesienia. Termin wynik wynosi 1, jeśli Δamplen ≤ 50 nukleotydów, 0, jeśli Δamplen ≥ 100 i (100 – Δamplen) / 50, Jeśli 50 < Δamplen < 100.,

wybór kryteriów punktacji i domyślnego progu opiera się na wcześniejszej literaturze . Jednak zarówno progi, jak i przedziały tolerancji rozmytej mogą być ustawione przez użytkownika inaczej niż domyślne i zgodnie z jego eksperymentalnymi potrzebami poprzez określenie żądanych wartości jako parametrów wejściowych podczas wywoływania narzędzia wiersza poleceń.

pokrycie

wynik pokrycia jest zdefiniowany jako liczba sekwencji 16S dopasowanych przez co najmniej jeden element., Biorąc pod uwagę sekwencje podkładu i bakterii 16S, definiujemy nasiona ostatnich 5 nukleotydów na 3 ' – końcu podkładu i rozważamy sekwencję 16S dopasowaną do podkładu, jeśli istnieje region sekwencji 16S, który pasuje i) ziarno podkładu dokładnie; I ii) pozostała część podkładu z maksymalnie 2 niedopasowaniami . Sekwencja 16S z zestawu odniesienia jest uważana za pokrytą parą Grunt-zestaw, jeśli co najmniej jeden Grunt przedni i jeden Grunt Odwrotny w parze Grunt-zestaw pasują do sekwencji., Ponieważ wydajność PCR maleje wraz z długością amplikonu, nakładamy kolejne ograniczenie: biorąc pod uwagę parę primer-set – i zestaw odniesienia sekwencji 16S, szacujemy docelową długość amplikonu jako medianę długości wszystkich amplikonów uzyskanych przez dopasowanie wszystkich kombinacji podkładów przednich i odwrotnych z pary primer-set-z zestawem odniesienia. Następnie uznajemy, że nie obejmuje wszystkich sekwencji referencyjnych 16S, których długość amplikonu różni się od długości docelowej o więcej niż 100 nukleotydów (dłuższych lub krótszych).,

Matching-bias

biorąc pod uwagę zestaw referencyjny sekwencji 16S i parę primer-set -, trzeci wynik optymalizacji mierzy zmienność liczby kombinacji pierwszych i odwrotnych primerów pasujących do każdej sekwencji referencyjnej 16S. Zmienność zasięgu spowodowana dopasowaniem powinna być zminimalizowana lub przynajmniej uwzględniona, gdy badanie ma na celu oszacowanie względnej obfitości różnych gatunków bakterii, ze względu na odchylenie amplifikacji w stosunku do gatunków objętych większą liczbą kombinacji podkładów do przodu i do tyłu., Jako miarę dopasowania biasu wykorzystujemy współczynnik zmienności pokrycia w sekwencjach docelowych, obliczany jako odchylenie standardowe ponad średnią liczby kombinacji pasujących do każdej sekwencji.

Zestaw referencyjny sekwencji 16S, przygotowanie i adnotacja

aby zoptymalizować trzy powyższe wyniki, opieramy się na reprezentatywnym zestawie bakteryjnych sekwencji 16S wyodrębnionych z publicznej bazy sekwencji 16s, GreenGenes ., Baza sekwencji GreenGenes 16S jest zorganizowana w operacyjne jednostki taksonomiczne (OTUs), które są zagnieżdżonymi klastrami sekwencji w bazie danych, zorganizowanymi na różnych poziomach podobieństwa między klastrami. Dla każdego poziomu podobieństwa, Sekwencja odniesienia jest powiązana z każdym klastrem, maksymalnie podobna do wszystkich innych sekwencji w tym samym klastrze . Zbiór sekwencji referencyjnych można zatem uznać za reprezentatywny podzbiór całej bazy sekwencji, stając się coraz dokładniejszym dla zwiększenia poziomów podobieństwa między klastrami (a tym samym liczby sekwencji referencyjnych)., Jako dobry kompromis między reprezentatywnością a złożonością wybraliśmy poziom podobieństwa między klastrami wynoszący 85%, odpowiadający zestawowi 5088 sekwencji reprezentatywnych do oceny kryteriów optymalizacji.

chociaż jest bardzo wrażliwa na opisywanie domen bakterii i archeonów, Taksonomia Greengenesa nie jest zaprojektowana do rozróżniania sekwencji należących do eukariotów lub wirusów., Z tego powodu zdecydowaliśmy się na ponowną adnotację sekwencji bakteryjnych 16S, wykorzystując oryginalną taksonomię NCBI, aby dokładnie zidentyfikować, wśród sekwencji reprezentatywnych, tylko te należące do domeny bakterii. Ponieważ brak informacji o domenie w adnotacji NCBI dla około 20% sekwencji, zaprojektowaliśmy procedurę ad hoc w celu identyfikacji sekwencji bakteryjnych wśród nich. Procedura jest szczegółowo opisana w materiałach dodatkowych(patrz plik dodatkowy 1)., Zachowawczo zdecydowaliśmy się wziąć pod uwagę tylko sekwencje opatrzone adnotacją jako bakterie zarówno w naszej opatrzonej, opartej na NCBI adnotacji, jak i w oryginalnej adnotacji GreenGenes. W rezultacie powstał zestaw 4573 reprezentatywnych sekwencji 16S należących do domeny bakterii.,

algorytm optymalizacji

ponieważ problem optymalnego wyboru podkładów wymaga jednoczesnej optymalizacji różnych konkurencyjnych wyników, można go rzucić jako problem optymalizacji wielu celów, gdzie przestrzeń wyszukiwania jest zbiorem wszystkich możliwych par podkładów i funkcja punktacji, lub kryterium optymalizacji, może być zdefiniowana tak, aby zmaksymalizować wydajność i zasięg oraz zminimalizować dopasowanie-bias., Gdy zachodzi potrzeba jednoczesnej optymalizacji więcej niż jednego kryterium, a cele, które mają być optymalizowane, są sprzeczne, zwykle nie interesuje nas jedno rozwiązanie, ale raczej zestaw rozwiązań Pareto optymalnych, tj. zestaw rozwiązań, dla których żaden z celów nie może być ulepszony bez poświęcenia co najmniej jednego innego celu ., Wynikiem optymalizacji wielu celów nie jest już unikalna optymalna para primer-set -, jak w optymalizacji jednego celu, ale raczej zbiór par primer-set -, które nie są gorsze niż jakakolwiek inna para primer-set-i ściśle lepsze zgodnie z co najmniej jednym z kryteriów. Dokładniej, dla problemu optymalizacji Tri-objective polegającego na maksymalizacji wyników optymalizacji efektywności (E) i pokrycia (C) oraz zminimalizowaniu wyniku matching-bias (M), jak zdefiniowano w poprzedniej sekcji, kandydujące pary primer-set-są oceniane zgodnie z wektorem funkcji celu f = (f E ; F C ; fM)., Biorąc pod uwagę dwie pary P I p', mówimy, że p dominuje p' (p ≺ p') wtedy i tylko wtedy, gdy f (p) ≠ f (p'), Fe (P) ≥ Fe (p'), fC (P) ≥ fC (p') i fM (P) ≤ fM (p'). Jeśli p 'nie istnieje tak, że p' ≺ p, to Pareta-zbiór-p nazywa się Pareto-optymalną. W tym kontekście celem optymalnego wyboru primerów jest wyznaczenie (lub przybliżenie) zbioru wszystkich par Pareto-optymalnych primer-set -, których obraz w przestrzeni trójobiektywowej nazywany jest frontem Pareto .

aby wyszukać optymalny front Pareto opieramy się na dwufazowym iterowanym podejściu best improvement local search zaproponowanym przez Dubois-Lacoste et al., i skutecznie wykorzystywane w Sambo i in. i Borrotti i in. dla optymalnego wieloobiektywnego projektowania eksperymentów.

Wyszukiwanie lokalne zaczyna się od początkowego rozwiązania i iteracyjnie je udoskonala, stosując małe zmiany lokalne i oceniając za każdym razem ich wpływ na jakość rozwiązania; zatrzymuje się, gdy żadne dalsze zmiany lokalne nie mogą poprawić rozwiązania. Proces jest powtarzany z kilku różnych punktów początkowych i powraca najlepsze rozwiązanie, jakie kiedykolwiek znaleziono, jako przybliżenie nieznanego optimum ., Powszechnym rozszerzeniem wyszukiwania lokalnego do przypadku multi-objective jest rozpoczęcie od zestawu początkowych rozwiązań Pareto, próbkowanie jednego rozwiązania z przodu, Optymalizacja za pomocą wyszukiwania lokalnego przypadkowa skalaryzacja problemu, tj. liniowa kombinacja wyników optymalizacji z wagami próbkowanymi równomiernie losowo z jednostki simplex, aktualizacja frontu Pareto i iteracja do momentu spełnienia warunku zakończenia .,

procedura MULTI-OBJECTIVE-SEARCH, której pseudo-kod jest podany poniżej, otrzymuje jako wejście żądany zakres długości amplikonów (rangeamplen), reprezentatywny zestaw sekwencji 16S (repset), początkowy zestaw (prawdopodobnie zdegenerowanych) par primerów (init) i liczbę restartów (NRE). Procedura rozpoczyna się od wybrania z init wszystkich możliwych par podkładu o żądanej długości amplikonu, długości podkładu (od 17 do 21 nukleotydów) i domeny docelowej (bakteryjnej lub uniwersalnej).,

zdegenerowane pary primer są konwertowane na nie-zdegenerowane pary primer-set-i dodawane do archiwum. Procedura następnie iteruje czasy nrest, za każdym razem próbkowanie losowego primer-set-pair pstart z przodu Pareto i losowego wektora α względnych wag dla wyników optymalizacji, z ciężarami próbkowanymi równomiernie z jednostki simplex; procedura, następnie, rozwiązuje skalaryzację problemu wielu celów, tj. problem jednego celu, w którym kombinację liniową trzech celów z względnymi wag α jest zmaksymalizowana i dodaje wynik do archiwum., W tym celu wyniki efektywności, zasięgu i dopasowania bias są znormalizowane do maksimum, tak że każdy znormalizowany wynik mieści się w zakresie od 0 do 1, a matching-bias jest redefiniowany jako 1-matching-bias, dzięki czemu można go zmaksymalizować jak inne wyniki., długość amplikonu w zakresie

2 Dodaj do archiwum odpowiednie Nie zdegenerowane pary primer-set-

3 dla r = 1 do nrest

4 PF = Pareto-FRONT(archiwum)

5 próbka PStart z pf

6 próbka α z 3, z Σi ai = 1

7 P = Local-SEARCH(PStart , α , repset)

8 Dodaj p do archiwum

9 wróć do archiwum

optymalizacja pojedynczego celu jest uzyskiwana przy użyciu algorytmu wyszukiwania lokalnego best improvement : zaczynając od początkowej pary primer-set-pair, algorytm wyszukiwania lokalnego cykli przez primery pary set-i, dla każdego primer, skanuje jego sąsiedztwo, i.,e. zbiór wszystkich możliwych lokalnych perturbacji podkładu. Perturbacje lokalne polegają na wszystkich możliwych przewrotach jednego nukleotydu (ocena trzech innych możliwych zasad) oraz wszystkich możliwych dodatkach i odjęciach jednego nukleotydu na końcach., Wyszukiwanie w przestrzeni rozwiązań odbywa się z najlepszą metodą wyszukiwania lokalnego: po wygenerowaniu całej dzielnicy, jak wyjaśniono powyżej, algorytm wybiera najlepszą perturbację sąsiada, rozpoczyna od niej generowanie następnej dzielnicy i powtarza aż do osiągnięcia rozwiązania, dla którego nie można znaleźć lepszej perturbacji sąsiada. Procedura kończy się, gdy nie można zastosować dalszych lokalnych ulepszeń do żadnego podkładu w parze primer-set -., Funkcja WEIGHTED-SCORE oblicza trzy wyniki optymalizacji z pary primer-set-i zestawu reprezentatywnego, mnoży wyniki przez względne wagi α i Zwraca sumę wyników.

opracowaliśmy narzędzie programowe wdrażające nasze podejście i wydaliśmy je na licencji GNU General Public License jako narzędzie programowe mopo16S (Multi-Objective Primer Optimization for 16S experiments) pod adresemhttp://sysbiobig.dei.unipd.it/?q=Software#mopo16S., mopo16S jest zaimplementowany jako wielowątkowe narzędzie wiersza poleceń C++; narzędzie programowe opiera się na wydajnych algorytmach i strukturach danych z biblioteki SeqAn i wykorzystuje bibliotekę openMP do wielowątkowości.,

4 dla i = 1 do |pcurr|

5 pri = i-ten podkład pcurr

6 dla pnew = pcurr ze wszystkimi możliwymi dodatkami i usunięciami podstawy na końcach i zamiennikami podstawy pri

7 scorenew = ważony wynik(pnew , α , repset)

8 jeśli scorenew >scorebest

9 pbest = pnew

10 scorebest = scorenew

11 pcurr = pbest

12 Return pcurr

najnowocześniejsze pary primer jako początkowe rozwiązania

wybraliśmy bazę online database probebase jako źródło kandydujących par primer-Set-do użycia jako początkowe rozwiązania przez mopo16s., Baza danych zawiera ponad 500 par (prawdopodobnie zdegenerowanych) podkładów i raporty dla każdego podkładu jego sekwencji, nici i pozycji, w której pasuje do referencyjnego genu 16S Escherichia coli, oraz domeny docelowej, dla której jest przeznaczony (są to bakterie, Archeea lub uniwersalne).,

biorąc pod uwagę pożądany zakres dla docelowej długości amplikonu jako wejście mopo16S, wybraliśmy wszystkie pary podkładów z bazy danych sonbase spełniające wszystkie następujące właściwości:

Długość amplikonu w pożądanym zakresie;
długość obu podkładów większa lub równa 17 nt i mniejsza lub równa 21 nt;
bakterie lub uniwersalna domena docelowa obu podkładów.,

ponieważ naszym podejściem jest praca z zestawami nie-zdegenerowanych primerów, w przypadku degeneracji w primer forward lub reverse, zastępujemy zdegenerowany primer odpowiednim zestawem Nie-zdegenerowanych primerów, otrzymanym przez przypisanie wszystkich możliwych kombinacji wartości do zdegenerowanych nukleotydów w primer. Przykład tej procedury przedstawiono w tabeli 1.

obliczyliśmy trzy wyniki dla każdej pary primer-set-i zidentyfikowaliśmy, wśród nich, pary primer-set-tworzące początkowy przód Pareto.