Identifying A- and P-site locations on ribosome-protected mRNA fragments using Integer Programming

Integer Programming Algorithm

W analizie danych Ribo-Seq, fragmenty mRNA są początkowo wyrównywane do transkryptomu referencyjnego, a ich lokalizacja jest raportowana w odniesieniu do ich 5′ końca. Oznacza to, że jeden fragment wniesie jeden odczyt, który jest zgłaszany na współrzędnej genomu, do której wyrównany jest 5′ koniec nukleotydu tego fragmentu (Rys. 1A). W danych Ribo-Seq obserwuje się fragmenty o różnej długości, które mogą wynikać z niepełnego trawienia RNA i ze stochastycznej natury rozszczepiania mRNA przez RNazę użytą w eksperymencie (Rys. 2, Supplementary Fig. S1). Głównym wyzwaniem w ilościowej analizie danych Ribo-Seq jest zidentyfikowanie na podstawie tych odczytów Ribo-Seq, gdzie znajdowały się miejsca A i P w czasie trawienia. Jest to nietrywialne, ponieważ niekompletne trawienie i stochastyczne rozszczepianie może wystąpić na obu końcach fragmentu. Na przykład, trawienie mRNA prowadzące do powstania fragmentu o wielkości 29 nt może zachodzić na różne sposoby, z których dwa zilustrowano na Rys. 1B. Wielkością, którą musimy dokładnie oszacować, jest liczba nukleotydów oddzielających kodon w miejscu A od 5′ końca fragmentu, którą nazywamy przesunięciem i oznaczamy Δ. Znajomość Δ określa pozycję miejsca A, jak również miejsca P, ponieważ miejsce P zawsze będzie na Δ minus 3 nt.

Rysunek 1

Lokalizacja miejsca A może być określona jako przesunięcie od 5′ końca fragmentów chronionych przez rybosom. (A) Schematyczne przedstawienie rybosomu translacyjnego (rysunek górny) i przesunięcia ∆ między odczytami Ribo-Seq mapowanymi w odniesieniu do 5′ końca śladów i wyśrodkowanymi na miejscu A (niebieskie słupki). Rybosom jest pokazany jako chroniący 28 nt fragment z jego 5′ końcem w ramce odczytu 0, jak określono od kodonu startowego ATG genu. Wskazane są miejsca E-, P- i A w obrębie rybosomu. Odczyty są następnie przesunięte z końca 5′ do miejsca A o wartość przesunięcia ∆. (B) Stochastyczne trawienie nukleazami może prowadzić do powstania różnych fragmentów. Dwa najbardziej prawdopodobne warianty 29 nt śladu z końcem 5′ w ramce 1 są pokazane z ich granicami odwzorowanymi przez przerywane linie wyrównujące do genomu, co może skutkować przesunięciami odpowiednio o 15 nt (góra) i 18 nt (dół). (C) Aby zilustrować zastosowanie algorytmu programowania całkowitoliczbowego, rozważmy hipotetyczny transkrypt o długości 60 nt. Pierwszy panel pokazuje profil rybosomu pochodzący z odczytów przypisanych do 5′ końca fragmentów o rozmiarze 33 w ramce 0. Kodon startu i kodon stopu są zaznaczone, podczas gdy reszta regionu CDS jest pokolorowana na jasnobrzoskwiniowo. Algorytm przesuwa ten profil rybosomu o 3 nt i oblicza funkcję celu \(\,T({{rm{Delta }}|i,S,F)\). Zakres przesunięcia to przesunięcie Δ. Podane są wartości \(\,T({\rm{Delta }}|i,S,F)\) dla Δ = 12, 15, 18, 21 ntów. W tym przykładzie średnia liczba odczytów na kodon wynosi 7,85. Różnica między dwoma najlepszymi przesunięciami, 18 (T = 222) i 15 (T = 215), jest mniejsza niż średnia. W związku z tym sprawdzamy kryteria drugorzędne (Wyniki). Przesunięcie 18 spełnia kryteria, że liczba odczytów w kodonie startowym jest mniejsza niż jedna piąta średniej liczby odczytów w drugim, trzecim i czwartym kodonie, a także, że liczba odczytów w drugim kodonie jest większa od liczby odczytów w trzecim kodonie. Stąd, Δ = 18 nt jest optymalnym offsetem dla tego transkryptu.

Rysunek 2

rozkład wielkości fragmentów mRNA dla S. cerevisiae Ribo-Seq dataset od Pop i współpracowników (A) i Pooled dataset (B).

Nasze rozwiązanie tego problemu opiera się na biologicznym fakcie, że dla transkryptów kanonicznych, bez translacji upstream, miejsce A aktywnie translujących rybosomów musi znajdować się pomiędzy drugim kodonem a kodonem stop CDS17. Dlatego optymalną wartością offsetu Δ dla fragmentów o określonym rozmiarze (S) i ramce odczytu (F) jest ta, która maksymalizuje całkowitą liczbę odczytów \(T({{rm{Delta }}|i,S,F)\) pomiędzy tymi kodonami dla każdego genu i, na który mapowane są fragmenty. Rozmiar fragmentu mRNA S mierzony jest w nukleotydach, a ramka F ma wartości 0, 1 lub 2 określone przez kodon startowy genu ATG i odpowiada ramce, w której znajduje się 5′ końcowy nukleotyd fragmentu (Rys. 1A). Ramka 5′ końca F jest wynikiem trawienia RNazą i jest różna od ramki odczytu rybosomu, który zwykle dokonuje translacji in-frame (ramka 0 miejsca A). Innymi słowy, dla każdej kombinacji (S, F) przesuwamy profil odczytu wyrównanego do 5′ o 3 nukleotydy na raz (aby zachować ramkę odczytu F) aż do zidentyfikowania wartości ∆, która maksymalizuje ilość odczytów pomiędzy drugim i końcowym kodonem (Rys. 1C, patrz następny podrozdział). Procedura ta jest przeprowadzana systematycznie dla każdego rozmiaru fragmentu S i ramki odczytu F oddzielnie, ponieważ każdy z nich może mieć (i stwierdzamy, że niektóre mają) inne optymalne ∆.

Podczas identyfikacji wartości Δ′ dla każdego genu w naszym zbiorze danych, minimalizujemy również występowanie wyników fałszywie pozytywnych poprzez zapewnienie, że najwyższy wynik, ∆,T({{rm{Delta }}^{prime} |i,S,F)\), jest znacząco wyższy niż następny najwyższy wynik, \(T({{rm{Delta }}^{prime} |i,S,F)\), który występuje przy innym przesunięciu Δ″. Jeśli różnica między dwoma najlepszymi wynikami jest mniejsza niż średnia liczba odczytów na kodon, stosujemy następujące dodatkowe kryteria selekcji. Aby wybrać między Δ′ i Δ″, wybieramy ten, który daje liczbę odczytów w kodonie startowym mniejszą o co najmniej jedną piątą od średniej liczby odczytów w drugim, trzecim i czwartym kodonie. Ponadto wymagamy, aby drugi kodon miał większą liczbę odczytów niż kodon trzeci. Biologiczne podstawy tych dodatkowych kryteriów są takie, że prawdziwe przesunięcie (tj. rzeczywista lokalizacja miejsca A) nie może być zlokalizowane w miejscu kodonu startowego, oraz że liczba odczytów przy drugim kodonie powinna być średnio wyższa niż przy trzecim, ze względu na wkład z etapu inicjacji translacji, podczas którego rybosom gromadzi się na mRNA z kodonem startowym w miejscu P. Poniżej pokazujemy, że wyniki naszej metody są odporne na zmiany tych progów.

Ilustracja procedury optymalizacji programowania całkowego

Rozkład wielkości fragmentów i ramek fragmentów chronionych przez rybosom (Rys. 2) w S. cerevisiae nie zależy od genu (Supplementary Fig. S2), a zatem wartości offsetu Δ również nie powinny zależeć od genu. Zatem lokalizacja miejsca A, względem 5′ końca fragmentu o rozmiarze S i ramce F, odpowiada najbardziej prawdopodobnej wartości offsetu dla wszystkich genów w zbiorze danych.

Lokalizacje miejsc A w danych S. cerevisiae Ribo-Seq zależą od rozmiaru fragmentu i ramki

W pierwszej kolejności zastosowaliśmy metodę programowania całkowego do danych Ribo-Seq z S. cerevisiae opublikowanych przez Popa i współpracowników19. Dla każdej kombinacji S i F zidentyfikowaliśmy najpierw te geny, które w odpowiadającym im profilu rybosomowym mają średnio co najmniej 1 odczyt na kodon. Liczba genów spełniających to kryterium jest podana w Tabeli Dodatkowej S1. Następnie zastosowaliśmy metodę programowania całkowitoliczbowego do tego podzbioru genów. Otrzymane rozkłady wartości Δ dla różnych kombinacji długości fragmentu i ramki pokazane są na Rys. 3A. Pokazujemy tylko wyniki dla fragmentów o rozmiarach od 27 do 33 nt, ponieważ ponad 90% odczytów mapuje w tym zakresie (Rys. 2A). Najbardziej prawdopodobna wartość offsetu dla wszystkich rozmiarów fragmentów od 20 do 35 nt jest przedstawiona w tabeli offsetów (Supplementary Table S2).

Rysunek 3

Rozkład wartości offsetów z algorytmu Integer Programming zastosowanego do transkryptów z S. cerevisiae. Dane wykreślone w (A) pochodzą z zestawu danych Pop, a (B) z zestawu danych Pooled. Rozkłady są wykreślone jako funkcja wartości offsetu i dla fragmentów o rozmiarach od 27 do 33 nt, są pokazane, od lewej do prawej, dla ramek 0, 1 i 2. Dla danego rozmiaru fragmentu i ramki, lokalizacja miejsca A znajduje się przy najbardziej prawdopodobnej wartości Δ w rozkładzie, pod warunkiem, że przesunięcie występuje dla więcej niż 70% genów (linie przerywane w panelach). Słupki błędów reprezentują 95% przedziały ufności obliczone przy użyciu Bootstrappingu. Rozmiary próbek podano w Tabeli Uzupełniającej S1.

Widzimy, że optymalna wartość Δ – czyli lokalizacja miejsca A – zmienia się dla różnych kombinacji S i F, z najbardziej prawdopodobnymi wartościami albo na 15 albo 18 nt. Tak więc, lokalizacja miejsca A zależy od S i F. W większości przypadków istnieje jeden dominujący pik dla danej pary wartości S i F. Na przykład, dla fragmentów o wielkości 27-30 nt w ramce 0, ponad 70% ich pergenowo zoptymalizowanych wartości Δ znajduje się 15 nt od 5′ końca tych fragmentów. Podobne wyniki uzyskano dla innych kombinacji, takich jak rozmiary 30, 31 i 32 nt w ramce 1 oraz 28 do 32 nt w ramce 2, gdzie zoptymalizowane wartości Δ wynoszą 18 nt. Tak więc, w całym transkryptomie, pozycja kodonu miejsca A na tych fragmentach jest jednoznacznie zidentyfikowana.

Istnieją jednak kombinacje S i F, które mają niejednoznaczne lokalizacje miejsca A w oparciu o te rozkłady. Na przykład, dla fragmentów o wielkości 27 nt w ramce 1, 47% zoptymalizowanych pod kątem genu wartości Δ znajduje się na 15 nt, podczas gdy 30% na 18 nt. Podobne wyniki obserwuje się dla fragmentów 28 i 29 nt w ramce 1, oraz 31 i 32 nt w ramce 0. Tak więc, dla tych kombinacji S i F istnieje podobne prawdopodobieństwo, że miejsce A znajduje się przy tym czy innym kodonie, a zatem wydaje się, że nie możemy jednoznacznie zidentyfikować lokalizacji miejsca A.

Wyższe pokrycie prowadzi do większej liczby unikalnych offsetów

Postawiliśmy hipotezę, że niejednoznaczność w identyfikacji miejsca A dla poszczególnych kombinacji S i F może być spowodowana niskim pokryciem (tj, słabą statystyką próbkowania). Aby przetestować tę hipotezę, połączyliśmy odczyty z różnych opublikowanych zbiorów danych Ribo-Seq w jeden zbiór danych z konsekwentnie wyższym pokryciem i większą liczbą genów, które spełniają nasze kryteria selekcji (Tabela uzupełniająca S1). Zastosowanie naszej metody do tego zbioru danych Pooled daje unikalne przesunięcia dla większej liczby kombinacji S i F w porównaniu do oryginalnego zbioru danych Pop (Rys. 3B i Tabela uzupełniająca S2), co jest zgodne z naszą hipotezą. Na przykład, dla fragmentów o rozmiarze 27 i ramce 1, mamy teraz unikalne przesunięcie 15 nt z 72% zoptymalizowanych dla genów wartości Δ na 15 nt (Rys. 3B). Jednakże, nadal widzimy niejednoznaczność obecną dla pewnych kombinacji (S, F).

Zastosowaliśmy dodatkową strategię w celu zwiększenia pokrycia poprzez ograniczenie naszej analizy do genów z większą średnią ilością odczytów na kodon. Jeśli hipoteza jest poprawna, to powinniśmy zaobserwować statystycznie istotny trend wzrostu najbardziej prawdopodobnej wartości Δ wraz ze wzrostem głębokości odczytu. Zastosowaliśmy tę analizę do zbioru danych Pooled i stwierdziliśmy, że niektóre początkowo niejednoznaczne kombinacje S i F stają się jednoznaczne wraz ze wzrostem pokrycia. Na przykład, przy średnio 1 odczycie na kodon, kombinacje (S, F) (25, 0), (27, 2) i (30, 1) są niejednoznaczne, ponieważ znajdują się poniżej naszego progu 70%. Jednakże widzimy statystycznie istotny trend (nachylenie = 0,5, p = 3,94 × 10-6) dla fragmentów (25, 0), że przesunięcie o 15 nt staje się bardziej prawdopodobne wraz ze wzrostem pokrycia, ostatecznie przekraczając próg 70% (Rys. 4A). Podobnie, dla (27, 2) (nachylenie = 0,58, p = 5,77 × 10-5) i (30, 1) (nachylenie = 0,25, p = 0,009) występuje trend w kierunku offsetu 18 nt, przy czym ponad 70% genów ma taki offset przy najwyższym pokryciu (Rys. 4B,C). Stąd, dla tych fragmentów, rosnące pokrycie jednoznacznie identyfikuje Δ′, a tym samym lokalizację miejsca A. Dla kilku kombinacji (S, F), jak (32, 0), niejednoznaczność nie jest rozwiązana nawet przy bardzo wysokim pokryciu (Rys. 4D), co, jak spekulujemy, może wynikać z nieodłącznych cech trawienia nukleazą, które są równie prawdopodobne dla więcej niż jednego offsetu.

Ryc. 4

Zwiększające się pokrycie identyfikuje lokalizacje miejsca A dla kombinacji S i F, które początkowo były niejednoznaczne. Wykreślono procent transkryptów o określonej wartości Δ dla różnych kombinacji S i F z zestawu danych Pooled dla S. cerevisiae. W każdym panelu wykreślono wiele rozkładów odpowiadających transkryptom o rosnącym pokryciu, wskazanym przez legendę na dole. Na przykład, rozkłady w kolorze niebieskim i czerwonym pochodzą z transkryptów, w których na każdy kodon przypada średnio co najmniej 1 lub 2 odczyty. Obserwujemy, że lokalizacja miejsca A ma tendencję do 15 nt dla S = 25, F = 0 (A) oraz do 18 nt dla S = 27, F = 2 (B) i S = 30, F = 1 (C). Dla S = 32, F = 0 (D), nie ma tendencji nawet przy wyższym pokryciu. Zauważ, że dla S = 27, F = 2 (panel B), istnieje mniej niż 10 genów ze średnią większą niż 50 odczytów na kodon, a zatem nie uwzględniamy punktu danych poza średnią większą niż 45 odczytów na kodon (patrz Metody). Słupki błędów reprezentują 95% przedziały ufności obliczone przy użyciu Bootstrapping.

Tak więc, wystarczająco wysokie pokrycie daje optymalną tabelę offsetów przedstawioną w Tabeli 1, gdzie offset jest najbardziej prawdopodobną lokalizacją miejsca A w stosunku do 5′ końca fragmentów mRNA wygenerowanych w S. cerevisiae.

Tabela 1 Lokalizacje miejsc A (przesunięcia nukleotydów względem 5′ końca) wyznaczone przez zastosowanie algorytmu Integer Programming do zbioru danych Pooled w S. cerevisiae są przedstawione jako funkcja rozmiaru fragmentu i ramki.

Spójność pomiędzy różnymi zestawami danych

Dane Ribo-Seq są wrażliwe na protokoły eksperymentalne, które mogą wprowadzać błędy w trawieniu i ligacji fragmentów chronionych przez rybosom. Łączenie zbiorów danych oferuje korzyść w postaci większego pokrycia, ale może maskować błędy specyficzne dla poszczególnych zbiorów danych. Aby określić, czy nasze unikalne przesunięcia (Tabela 1) są zgodne z wynikami z poszczególnych zestawów danych, zastosowaliśmy algorytm programowania całkowitoliczbowego do każdego indywidualnego zestawu danych. Większość z tych zbiorów danych ma niskie pokrycie, co skutkuje mniejszą liczbą genów spełniających nasze kryteria filtrowania (Dodatkowy plik S1). Dla każdego unikalnego przesunięcia w Tabeli 1, klasyfikujemy je jako zgodne z indywidualnym zestawem danych pod warunkiem, że najbardziej prawdopodobne przesunięcie z indywidualnego zestawu danych (nawet jeśli nie osiąga progu 70% z powodu ograniczeń w głębokości pokrycia) jest takie samo jak w Tabeli 1. Stwierdziliśmy, że zdecydowana większość unikalnych offsetów (22 z 24) z Tabeli 1 jest zgodna w 75% lub więcej poszczególnych zestawów danych (statystyki w Tabeli S3). Tylko dwie kombinacje (S, F) wykazują częste niespójności. Kombinacje (S, F) (27, 1) i (27, 2) są niespójne w 33% lub więcej poszczególnych zestawów danych (Tabela uzupełniająca S3). Sugeruje to, że badacze, którzy chcą zminimalizować liczbę fałszywych pozytywów, powinni odrzucić te kombinacje (S, F) podczas tworzenia profili rybosomów miejsca A.

Odporność tabeli offsetów na zmienność progów

Algorytm programowania całkowego wykorzystuje dwa progi do identyfikacji unikalnych offsetów. Jeden z nich zakłada, że 70% genów wykazuje najbardziej prawdopodobne przesunięcie, drugi, zaprojektowany w celu zminimalizowania fałszywych pozytywów wynikających z szumu próbkowania w danych Ribo-Seq, zakłada, że odczyty w pierwszym kodonie są mniejsze niż jedna piąta średnich odczytów w drugim, trzecim i czwartym kodonie. Chociaż istnieją dobre powody do wprowadzenia tych kryteriów progowych, dokładne wartości tych progów są arbitralne. Dlatego sprawdziliśmy, czy zmiana tych progów zmienia wyniki przedstawione w tabeli 1. Zmieniliśmy pierwszy próg na 60% i 80%, a następnie ponownie obliczyliśmy tabelę offsetów. Zgłaszamy, czy unikalne przesunięcie zmieniło się przez dodanie „R” lub „S” (odpowiednio dla solidnego i wrażliwego) obok zgłoszonego przesunięcia w Uzupełniającej Tabeli S3. Stwierdzamy, że dwie trzecie unikalnych kombinacji (S, F) nie uległo zmianie (Tabela uzupełniająca S3). Kombinacje (S, F) (25, 0), (25, 2), (27, 0), (27, 1), (28, 1), (31, 0), (33, 0) i (33, 2) stają się niejednoznaczne, gdy zwiększyliśmy próg do 80%.

Zmienialiśmy drugi, wspomniany wyżej próg z jednej piątej do jednej i w dół do jednej dziesiątej i stwierdzamy, że wszystkie unikalne (S, F) kombinacje z wyjątkiem (25, 2), (33, 0), (33, 2) i (34, 1) pozostają niezmienione (raportowane jako „R” w Uzupełniającej Tabeli S3). Tak więc, podsumowując, w zdecydowanej większości przypadków, unikalne przesunięcia raportowane w Tabeli 1 zależą w bardzo niewielkim stopniu od konkretnych wartości tych progów.

Testowanie algorytmu programowania całkowego przeciwko sztucznym danym Ribo-Seq

Aby przetestować poprawność i odporność naszego podejścia, wygenerowaliśmy zestaw danych symulowanych miejsc zajmowanych przez rybosomy w 4,487 transkryptach S. cerevisiae i zapytaliśmy, czy nasza metoda może dokładnie określić lokalizacje miejsc A. Sztuczne odczyty Ribo-Seq zostały wygenerowane z tych miejsc przy założeniu rozkładu Poissonowskiego w ich wartościach (S, F) z wykorzystaniem losowych długości śladów podobnych do tych, które znaleziono w eksperymentach (patrz Metody i Uzupełniające Rys. S3A, B). Zbadaliśmy zdolność naszej metody do poprawnego określenia prawdziwych lokalizacji miejsc A dla czterech różnych zestawów wstępnie zdefiniowanych wartości przesunięć (patrz Metody). Algorytm programowania całkowitoliczbowego został następnie zastosowany do wynikowych sztucznych danych Ribo-Seq. Stwierdziliśmy, że tabela offsetów wygenerowana na podstawie algorytmu odtwarza użyte offsety wejściowe (Supplementary Fig. S3C i Supplementary Table S4). Procedura ta została powtórzona dla różnych rozkładów długości odczytów, jak również z różnymi offsetami wejściowymi i stwierdzamy, że tabele offsetów wygenerowane przez nasz algorytm odtwarzają tabele offsetów wejściowych w ponad 93% wszystkich kombinacji (S, F) (Supplementary Fig. S3B,C i Supplementary File S2). Metoda identyfikuje niewielką liczbę niejednoznacznych offsetów ze względu na niskie pokrycie odczytów w ogonach rozkładów. Odkrycie to jeszcze bardziej podkreśla znaczenie pokrycia odczytów jako krytycznego czynnika w dokładnej identyfikacji miejsca A.

A-site offsets in mouse embryonic stem cells

Biologiczny fakt, że miejsce A rybosomu znajduje się tylko pomiędzy drugim i końcowym kodonem nie jest ograniczony do S. cerevisiae i dlatego algorytm programowania całkowego powinien być stosowany do danych Ribo-Seq z dowolnego organizmu. Dlatego też zastosowaliśmy naszą metodę do zbioru danych Pooled Ribo-Seq z mysich embrionalnych komórek macierzystych (mESCs). Wynikowa tabela przesunięć miejsc A wykazywała niejednoznaczne przesunięcia we wszystkich, z wyjątkiem trzech (S, F) kombinacji (Tabela Uzupełniająca S5). W mESCs występuje powszechna elongacja translacji, która zachodzi poza granicami anotowanych regionów CDS w upstreamowych otwartych ramkach odczytu (uORFs)20. Wzbogacenie chronionych przez rybosom fragmentów z tych translacyjnych uORF-ów może utrudnić naszemu algorytmowi znalezienie unikalnych offsetów, ponieważ mogą one wnosić odczyty wokół kodonu startowego kanonicznych, anotowanych CDS-ów. Dlatego postawiliśmy hipotezę, że jeśli zastosujemy nasz algorytm tylko do transkryptów pozbawionych uORF-ów i posiadających pojedyncze miejsce inicjacji, to nasz algorytm powinien zidentyfikować więcej unikalnych offsetów. Ingolia i współpracownicy11 doświadczalnie zidentyfikowali dla dobrze translokowanych transkryptów mESCs liczbę miejsc inicjacji oraz obecność uORF-ów. Dlatego wybraliśmy te geny, które mają tylko jedno miejsce inicjacji translacji w pobliżu anotowanego kodonu startu i dodatkowo ograniczyliśmy naszą analizę do transkryptów z jedną izoformą, ponieważ wiele izoform może mieć różne miejsca terminacji.

Zastosowanie algorytmu Integer Programming do tego zestawu genów zwiększa liczbę unikalnych offsetów z 3 do 13 (S, F) kombinacji (Tabela Uzupełniająca S6). Zastosowanie tych samych testów odporności i spójności jak w przypadku S. cerevisiae ujawnia, że 77% unikalnych offsetów jest odpornych na zmienność progową, a podobny odsetek jest spójny w obu indywidualnych zestawach danych użytych do stworzenia danych zbiorczych (Tabela Uzupełniająca S6). Zatem unikalne przesunięcia, które podajemy dla mESCs, są odporne i spójne w zdecydowanej większości zbiorów danych. Wynik ten wskazuje również, że skuteczna identyfikacja lokalizacji miejsc A wymaga analizy tylko tych transkryptów, które nie zawierają uORF-ów.

Integer Programming does not yield unique offsets for E. coli

Aby sprawdzić, jak szeroko możemy zastosować nasz algorytm, zastosowaliśmy go do danych Pooled Ribo-Seq z prokariotycznego organizmu E. coli. Liczba genów spełniających nasze kryteria filtrowania jest podana w Tabeli S7. MNaza, nukleaza stosowana w protokole E. coli Ribo-Seq, trawi mRNA w sposób stronniczy – preferując trawienie od końca 5′ w stosunku do końca 3′21,22. Dlatego, podobnie jak w innych badaniach21,22,23, zastosowaliśmy nasz algorytm w taki sposób, że zidentyfikowaliśmy lokalizację miejsca A jako przesunięcie od końca 3′ zamiast od końca 5′. Policystronowe mRNA (tj. transkrypty zawierające wiele CDS) mogą powodować problemy dla naszego algorytmu, ponieważ blisko rozmieszczone odczyty na granicach sąsiadujących CDS są oceniane dla różnych offsetów w obu CDS. Aby uniknąć niedokładnych wyników, ograniczyliśmy naszą analizę do 1915 transkryptów monocystronicznych, które nie posiadają żadnego innego transkryptu w odległości 40 nt upstream lub downstream od CDS. W oparciu o nasze doświadczenie w analizie zbioru danych mESCs, odfiltrowaliśmy transkrypty z wieloma miejscami inicjacji translacji, jak również transkrypty, których przypisane miejsca inicjacji zostały zakwestionowane. Nakahigashi i współpracownicy24 wykorzystali tetracyklinę jako inhibitor translacji do zidentyfikowania 92 transkryptów w E. coli z różnymi miejscami inicjacji w stosunku do anotacji referencyjnej. Wykluczamy te transkrypty również z naszej analizy. Jednakże, dla tego zbiorczego zestawu danych o wysokim pokryciu, znajdujemy niejednoznaczne przesunięcia dla wszystkich kombinacji (S, F) (Tabela uzupełniająca S5). Analiza meta-genu znormalizowanej gęstości rybosomów w CDS i 30 nt regionie upstream i downstream ujawnia sygnatury translacji poza granicami CDS (Supplementary Fig. S4), zwłaszcza wyższe niż przeciętne wzbogacenie odczytów kilka nukleotydów przed kodonem startu. Spekulujemy, że sparowanie zasad sekwencji Shine-Dalgarno (SD) z komplementarną sekwencją anty-SD w 16S rRNA25 chroni te kilka nukleotydów przed kodonem startu przed trawieniem rybonukleazą i dlatego powoduje wzbogacenie odczytów Ribo-Seq. Ponieważ te „pseudo” chronione przez rybosom fragmenty nie mogą być odróżnione od rzeczywistych chronionych przez rybosom fragmentów zawierających kodon z miejscem A rybosomu, nasz algorytm jest ograniczony w zastosowaniu do tych danych.

Reprodukcja znanych motywów PPX i XPP, które prowadzą do spowolnienia translacji

W S. cerevisiae26 i E. coli21,27 pewne motywy polipeptydowe PPX i XPP (w których X odpowiada dowolnemu z 20 aminokwasów) mogą przeciągać rybosomy, gdy trzecia reszta znajduje się w miejscu A. Czynniki elongacji eIF5A (w S. cerevisiae) i EF-P (w E. coli) pomagają złagodzić zahamowanie wywołane przez niektóre motywy, ale nie inne26. Nawet w mESCs, Ingolia i współpracownicy11 wykryli PPD i PPE jako silne motywy pauzujące. Dlatego sprawdziliśmy, czy nasze podejście może odtworzyć znane motywy przeciągania. Zrobiliśmy to, obliczając znormalizowaną gęstość odczytu przy różnych wystąpieniach motywu PPX i XPP.

W S. cerevisiae zaobserwowaliśmy duże gęstości rybosomów przy PPG, PPD, PPE i PPN (ryc. 5A), z których wszystkie zostały sklasyfikowane jako silne motywy przeciągania w S. cerevisiae26, a także w E. coli27. W przeciwieństwie do tego, w PPP średnio nie występuje przeciąganie, co jest zgodne z innymi badaniami26. Jest to najprawdopodobniej spowodowane działaniem eIF5A. Dla motywów XPP najsilniejsze zahamowanie zaobserwowano dla motywów GPP i DPP, co jest zgodne z wynikami uzyskanymi w S. cerevisiae i E. coli (Rys. 5B). W mESCs najsilniejsze zahamowanie obserwujemy przy PPE i PPD, co jest zgodne z wynikami uzyskanymi przez Ingolia i współpracowników11 (Supplementary Fig. S5A). Dla motywów XPP zaobserwowaliśmy bardzo słabe zahamowanie tylko dla DPP (Supplementary Fig. S5B). Tak więc, nasze podejście do mapowania miejsca A na śladach rybosomów umożliwia dokładne wykrywanie ustalonych pauz w translacji w poszczególnych motywach PPX i XPP nascent polipeptydów.

Rysunek 5

Kilka motywów PPX i XPP prowadzi do przeciągnięcia rybosomalnego w S. cerevisiae. Medianę znormalizowanej gęstości rybosomów otrzymano dla wszystkich przypadków (A) motywów PPX i (B) XPP, w których X odpowiada jednemu z 20 naturalnie występujących aminokwasów. Używając testu permutacyjnego, określamy, czy mediana gęstości rybosomów jest statystycznie istotna, czy też występuje przypadkowo. Statystycznie istotne motywy są zaznaczone na ciemnoczerwono. Analiza ta została przeprowadzona na zbiorze danych Pop dla transkryptów, w których co najmniej 50% pozycji kodonowych ma zmapowane odczyty. Słupki błędów są 95% przedziałami ufności dla mediany uzyskanej za pomocą Bootstrappingu.

W badaniu danych Ribo-Seq z komórek ssaków28 zaobserwowano niezależną od sekwencji pauzę translacji, gdy 5. kodon transkryptu znajduje się w miejscu P. Ta pauza po inicjacji została również zaobserwowana w badaniu in vitro syntezy polifenylalaniny, gdzie zaobserwowano zahamowanie, gdy 4. kodon znajdował się w miejscu P29. Z profilami miejsca A uzyskanymi przy użyciu naszych tabel przesunięć dla S. cerevisiae i mESCs; obserwujemy również te zdarzenia pauzy, gdy zarówno 4. i 5. kodon znajdują się w miejscu P (Supplementary Fig. S6).

Większa dokładność lokalizacji miejsca A niż inne metody

Nie ma niezależnej metody eksperymentalnej do weryfikacji dokładności zidentyfikowanych lokalizacji miejsca A przy użyciu naszej metody lub jakiejkolwiek innej metody4,5,6,7,8,9,10,12,30,31,32,33,34,35. Twierdzimy, że dobrze ugruntowana pauza rybosomu przy poszczególnych motywach sekwencji PPX jest najlepszym dostępnym środkiem do różnicowania dokładności istniejących metod. Powodem tego jest to, że te motywy przeciągania zostały zidentyfikowane w E. coli36,37 i S. cerevisiae38 poprzez ortogonalne metody eksperymentalne (w tym badania enzymologiczne i drukowanie toe), a dokładna lokalizacja miejsca A podczas takiego spowolnienia jest znana jako w kodonie kodującym trzecią resztę motywu 36. Zatem najdokładniejszą metodą identyfikacji miejsca A będzie ta, która najczęściej przypisuje większą gęstość rybosomów do X przy każdym wystąpieniu motywu PPX.

Zastosowaliśmy ten test do najsilniej hamujących motywów PPX, tj. PPG w S. cerevisiae i PPE w mESCs. W S. cerevisiae metoda programowania całkowitego daje największą gęstość rybosomów przy glicynowym kodonie motywu PPG, gdy jest stosowana zarówno do zbioru danych Pooled (Rys. 6A), jak i Pop (Supplementary Fig. S7A). Badając każde wystąpienie PPG w naszym zbiorze danych genów, stwierdzamy, że w większości przypadków nasza metoda przypisuje większą gęstość rybosomów do glicyny niż każda inna metoda, gdy jest stosowana zarówno do Pooled (Fig. 6B, Wilcoxon signed-rank test (n = 224), P < 0.0005 dla wszystkich metod z wyjątkiem Hussmanna (P = 0,164)) i Pop datasets (Supplementary Fig. S7B, Wilcoxon signed-rank test (n = 35), P < 10-5 dla wszystkich metod z wyjątkiem Hussmanna (P = 0,026) i Ribodeblur (P = 0,01)). Te same analizy zastosowane do mESCs w motywach PPE pokazują, że nasza metoda przewyższa pozostałe dziewięć metod (Fig. 6C,D) z naszą metodą przypisującą większą gęstość rybosomów przy kwasie glutaminowym dla co najmniej 85% motywów PPE w naszym zbiorze danych w porównaniu do wszystkich innych metod (Fig. 6D, Wilcoxon signed-rank test (n = 104), P < 10-15 dla wszystkich metod). Tak więc, dla S. cerevisiae i mESCs nasza metoda programowania całkowitoliczbowego jest dokładniejsza niż inne metody w identyfikacji miejsca A na fragmentach chronionych przez rybosom.

Rysunek 6

Algorytm programowania integralnego prawidłowo przypisuje większą gęstość rybosomów niż inne metody do glicyny w motywach PPG w S. cerevisiae i do kwasu glutaminowego w motywach PPE w mESCs. (A) Znormalizowana gęstość rybosomów uzyskana przy użyciu różnych metod stosowanych do identyfikacji miejsca A jest pokazana dla przypadku motywu PPG w genie YLR375W z G w pozycji kodonu 303 w zbiorze danych Pooled z S. cerevisiae (legenda wskazuje metodę, a pełne szczegóły dla każdej metody można znaleźć w sekcji Metody). (B) Frakcja instancji PPG (n = 224), w których metoda programowania całkowitoliczbowego daje większą gęstość rybosomów przy glicynie w porównaniu do każdej innej metody. Kodowanie kolorami jest takie samo jak w legendzie w panelu (A). Nasza metoda radzi sobie lepiej, jeśli przypisuje większą gęstość rybosomów w więcej niż połowie przypadków (linia pozioma w panelu B). Metoda programowania integralnego radzi sobie lepiej niż wszystkie inne metody (P < 0,0005), z wyjątkiem Hussmanna, który nie różni się statystycznie (P = 0,164). (C) Znormalizowana gęstość rybosomów jest pokazana dla instancji motywu PPE w genie uc007zma.1 z E w pozycji kodonu 127 w Pooled dataset of mouse ESCs (patrz Legenda i tekst główny dla szczegółów na temat metod). (D) Frakcja przypadków PPE, w których metoda programowania całkowitego daje większą gęstość rybosomów w kwasie glutaminowym w porównaniu do każdej innej metody. Kodowanie kolorami jest takie samo jak w legendzie panelu (C). Metoda programowania integralnego wypada lepiej niż wszystkie inne metody (P < 10-15) w dokładnym przypisywaniu gęstości rybosomów do kwasu glutaminowego w motywach PPE (n = 104). Dla analiz przedstawionych w (B) i (D), dwustronne wartości p zostały obliczone przy użyciu testu rang Wilcoxona. Słupki błędów reprezentują 95% przedział ufności wokół mediany obliczonej przy użyciu Bootstrappingu.

Duża liczba czynników molekularnych wpływa na szybkość translacji kodonów i gęstość rybosomów wzdłuż transkryptów39. Jednym z czynników jest stężenie tRNA kognatycznego, ponieważ kodony dekodowane przez tRNA kognatowe o wyższym stężeniu powinny mieć średnio niższą gęstość rybosomów15,16,40. Dlatego też, jako dodatkowy test jakościowy, oczekujemy, że najdokładniejsza metoda miejsca A da największą antykorelację pomiędzy gęstością rybosomów w kodonie a stężeniem tRNA kognatycznego. Ten test jest tylko jakościowy, ponieważ na korelację między gęstością rybosomu w kodonie a stężeniem poznawczego tRNA mogą mieć wpływ inne czynniki, w tym wykorzystanie kodonu i ponowne wykorzystanie naładowanych tRNA w pobliżu rybosomu41,42. Używając obfitości tRNA oszacowanych wcześniej na podstawie eksperymentów RNA-Seq na S. cerevisiae16, stwierdziliśmy, że nasza metoda programowania integralnego daje największą antykorelację w porównaniu z jedenastoma innymi rozważanymi metodami (Tabela uzupełniająca S8), co dodatkowo wspiera dokładność naszej metody. Nie byliśmy w stanie przeprowadzić tego testu na mESCs, ponieważ pomiary stężenia tRNA nie zostały odnotowane w literaturze.

.

Dodaj komentarz

Twój adres e-mail nie zostanie opublikowany.