Wykrywanie rzadkich wariantów genomowych z sekwencjonowania zbiorczego za pomocą SPLINTER

Summary

Sekwencjonowanie DNA w puli to szybka i opłacalna strategia wykrywania rzadkich wariantów związanych ze złożonymi fenotypami w dużych kohortach. W tym miejscu opisujemy analizę obliczeniową zbiorczego sekwencjonowania nowej generacji 32 genów związanych z rakiem przy użyciu pakietu oprogramowania SPLINTER. Ta metoda jest skalowalna i ma zastosowanie do każdego fenotypu będącego przedmiotem zainteresowania.

Abstract

W miarę jak technologia sekwencjonowania DNA znacznie się rozwinęła w ostatnich latach², staje się coraz bardziej oczywiste, że ilość zmienności genetycznej między dowolnymi dwoma osobnikami jest większa niż wcześniej sądzono3. W przeciwieństwie do tego, genotypowanie oparte na tablicach nie zidentyfikowało znaczącego udziału wspólnych wariantów sekwencji w zmienności fenotypowej powszechnej choroby^4,5. Podsumowując, obserwacje te doprowadziły do ewolucji hipotezy powszechnej choroby / rzadkiego wariantu, sugerującej, że większość "brakującej odziedziczalności" w powszechnych i złożonych fenotypach wynika z osobistego profilu danej osoby w zakresie rzadkich lub prywatnych wariantów DNA^6-8. Jednak scharakteryzowanie, w jaki sposób rzadka zmienność wpływa na złożone fenotypy, wymaga analizy wielu osób dotkniętych chorobą w wielu loci genomu i jest idealnie porównane z podobnym badaniem w niedotkniętej chorobą kohorcie. Pomimo możliwości sekwencjonowania oferowanych przez dzisiejsze platformy, populacyjne badanie wielu loci genomowych i wymagana następnie analiza obliczeniowa pozostają zaporowe dla wielu badaczy.

Aby zaspokoić tę potrzebę, opracowaliśmy metodę sekwencjonowania zbiorczego^1,9 oraz nowy pakiet oprogramowania¹ do bardzo dokładnego wykrywania rzadkich wariantów na podstawie uzyskanych danych. Możliwość łączenia genomów z całych populacji osób dotkniętych chorobą i badania stopnia zmienności genetycznej w wielu docelowych regionach w jednej bibliotece sekwencjonowania zapewnia doskonałą oszczędność kosztów i czasu w porównaniu z tradycyjną metodologią sekwencjonowania pojedynczej próby. Przy średnim pokryciu sekwencjonowania na allel wynoszącym 25-krotnie, nasz niestandardowy algorytm, SPLINTER, wykorzystuje wewnętrzną strategię kontroli wywoływania wariantów do wywoływania insercji, delecji i podstawień o długości do czterech par zasad o wysokiej czułości i swoistości z pul do 1 zmutowanego allelu na 500 osób. W tym miejscu opisujemy metodę przygotowania biblioteki sekwencjonowania w puli, a następnie instrukcje krok po kroku, jak używać pakietu SPLINTER do analizy sekwencjonowania w puli (http://www.ibridgenetwork.org/wustl/splinter). Pokazujemy porównanie sekwencjonowania zbiorczego 947 osób, z których wszystkie również przeszły macierz całego genomu, przy ponad 20 kb sekwencjonowania na osobę. Zgodność między genotypowaniem oznaczonych i nowych wariantów wywołanych w próbie zbiorczej była doskonała. Metodę tę można łatwo przeskalować do dowolnej liczby loci genomowych i dowolnej liczby osobników. Poprzez włączenie wewnętrznych dodatnich i ujemnych kontroli amplikonów w proporcjach naśladujących badaną populację, algorytm można skalibrować w celu uzyskania optymalnej wydajności. Strategia ta może być również modyfikowana do użytku z wychwytywaniem hybrydyzacji lub indywidualnymi kodami kreskowymi i może być stosowana do sekwencjonowania naturalnie heterogenicznych próbek, takich jak DNA guza.

Protocol

Ta metoda została użyta w badaniach opisanych w Vallania FML et al. Badania genomu 2010.

1. Łączenie próbek i przechwytywanie metodą PCR docelowych loci genomowych

1. Połącz znormalizowaną ilość genomowego DNA od każdej osoby w swojej puli (pulach). Użycie 0,3 ng DNA na osobę na reakcję PCR włączy około 50 diploidalnych genomów na osobę do każdej reakcji PCR, co zwiększa prawdopodobieństwo jednolitej amplifikacji na allel w puli.
2. Sekwencje genomowe można uzyskać z NCBI (http://www.ncbi.nlm.nih.gov/) lub UCSC Genome Browser (http://genome.ucsc.edu/index.html). Upewnij się, że używasz "RepeatMasker" (oznaczonego jako "N") podczas uzyskiwania sekwencji, aby uniknąć projektowania podkładu w powtarzającym się regionie.
3. Użyj internetowego narzędzia Primer3 (http://frodo.wi.mit.edu/primer3/input.htm) do projektowania starterów poprzez wycinanie i wklejanie interesujących regionów genomu oraz niektórych sekwencji flankujących (amplikony 600-2000 pz są zazwyczaj idealne). Optymalne warunki projektowania podkładu dla stosowania podkładu 3 to¹⁰ : Minimalny rozmiar podkładu = 19; Optymalna wielkość podkładu = 25; Maksymalny rozmiar podkładu = 30; Minimum Tm = 64 °C; Optimum Tm = 70 °C; Maksymalne Tm = 74 °C; Maksymalna różnica Tm = 5 °C; Minimalna zawartość GC = 45; Maksymalna zawartość GC = 80; Liczba do zwrócenia = 20 (jest to dowolne); Maksymalna stabilność końcowa 3' = 100. Zaprojektuj startery, aby amplifikować wszystkie interesujące nas loci genomowe. Po otrzymaniu starterów, liofilizowane masy można rozcieńczyć w 10 mM Tris, pH 7,5 + 0,1 mM EDTA do końcowego stężenia 100 uM, a następnie dodatkowo rozcieńczyć 10:1 w ddH₂O do 10 uM.
4. Amplifikacja PCR: Zalecamy stosowanie polimerazy DNA o wysokiej wierności do amplifikacji dużych amplikonów genomowych ze względu na niski poziom błędów (^10-7) i generowanie produktów końcowych (jest to konieczne w przypadku dalszego etapu ligacji). Użyliśmy PfuUltra High-Fidelity, ale enzymy o podobnych właściwościach (takie jak Phusion) powinny zapewnić porównywalne wyniki. Każda reakcja PCR zawiera końcowe stężenie 2,5 U polimerazy PfuUltra High-Fidelity, 1 M Betainy, 400 nM każdy starter, 200 μM dNTP, 1x bufor PfuUltra (lub bufor zawierający ≥ 2 mM Mg²⁺ w celu utrzymania wierności enzymatycznej), 5-50 ng połączonego DNA w końcowej objętości 50 μL. Stosować następujące warunki PCR: 1. 93-95 °C przez 2 minuty; 2. 93-95 °C przez 30 sekund; 3. 58-60 °C przez 30 sekund; 4. 65-70 °C przez 60-90 sekund dla amplikonów 250-500 pz / 1,5-3 minuty dla amplikonów 500-1000 pz / 3-5 minut dla amplikonów > 1 kb; 5. Powtórz kroki 2-4 przez 25-40 cykli; 6. 65 °C przez 10 minut; 7. Trzymaj w 4 °C. W razie potrzeby wyniki PCR można zazwyczaj poprawić poprzez: 1) obniżenie temperatury wyżarzania dla małych amplikonów; 2) podniesienie temperatury wyżarzania dla dużych amplikonów; 3. Wydłużenie czasu przedłużenia dla dowolnego wzmacniacza.
5. Przygotowanie kontroli SPLINTER: Każdy eksperyment SPLINTER wymaga obecności kontroli ujemnej i pozytywnej w celu uzyskania optymalnej dokładności. Kontrola ujemna może składać się ze wszystkich homozygotycznych pozycji podstawowych w dowolnej pojedynczej próbce oznaczonej kodem kreskowym, która została wcześniej zsekwencjonowana (np. próbka HapMap). Kontrola pozytywna składałaby się wówczas z mieszaniny dwóch lub więcej takich próbek. W tym raporcie kontrola ujemna to amplifikowany region o długości 1,934 pz ze szkieletu wektora ssDNA M13mp18. Produkt PCR został zsekwencjonowany metodą Sangera przed jego użyciem w celu potwierdzenia, że nie ma żadnych zmian sekwencji w stosunku do materiału źródłowego lub amplifikacji PCR. Kontrola pozytywna składa się z panelu wektorów pGEM-T Easy ze sklonowaną wstawką 72 pz, zaprojektowaną z określonymi insercjami, delecjami, podstawieniami (Tabela 1). Mieszamy wektory ze sobą na tle typu dzikiego w proporcjach molowych tak, aby mutacje były obecne z częstotliwością pojedynczego allelu w puli (tj. dla puli 100 alleli częstość pojedynczego allelu wyniesie 1%). Następnie PCR amplifikujemy mieszaną matrycę kontrolną przy użyciu miejsc startera M13 PUC w pGEM-T Easy, generując końcowy produkt PCR o długości 355 pz.

2. Przygotowanie i sekwencjonowanie biblioteki PCR w puli

1. Grupowanie produktów PCR: Każdy produkt PCR powinien zostać oczyszczony z nadmiaru starterów. Użyliśmy oczyszczania kolumny Qiagen Qiaquick lub 96-dołkowych płyt filtracyjnych z kolektorem próżniowym do czyszczenia na dużą skalę. Po oczyszczeniu każdy produkt PCR powinien zostać oznaczony ilościowo przy użyciu standardowych technik. Połącz każdy produkt PCR (w tym kontrole) w pulę znormalizowaną przez liczbę cząsteczek, ponieważ łączenie według stężenia spowoduje nadreprezentację małych amplikonów w porównaniu z większymi produktami. Stężenia przelicza się na bezwzględną liczbę cząsteczek DNA na objętość za pomocą wzoru: (g / μL) x (1 mol x bp / 660 g) x (1 / # bp w amplikonie) x (6 x 10²³ cząsteczki / 1 mol) = cząsteczki / μL. Następnie określamy objętość z każdej reakcji wymaganą do zebrania znormalizowanej liczby cząsteczek na amplikon. Liczba ta jest dowolna, można ją dostosować i tak naprawdę zależy od objętości pipetowania wystarczająco dużej, aby zachować dokładność. Zazwyczaj łączymy 1-2 x 10¹⁰ cząsteczek każdego amplikonu.
2. Podwiązanie produktów PCR: Ten krok jest niezbędny do uzyskania jednolitego pokrycia sekwencjonowania, ponieważ sonikacja małych amplikonów PCR będzie stronnicza w kierunku ich reprezentacji. Aby temu zaradzić, przed fragmentacją podwiązujemy połączone produkty PCR do dużych konkamerów (> = 10 Kb). Polimeraza Pfu Ultra HF generuje końce, co prowadzi do wydajnej ligacji (polimeraza na bazie Taq doda zwis 3p "A", który nie pozwoli na ligację bez wcześniejszego wypełnienia lub stępienia). W razie potrzeby reakcję tę można zwiększyć 2-3 krotnie. Reakcja ligacji zawiera 10 kinazy polinukleotydowej U T4, 200 ligazy U T4, 15% w/v polietylenu, bufor ligazy 1X T4, glikol 8000 MW, do 2 μg zbiorczych produktów PCR w końcowej objętości 50 μL. Reakcje inkubuje się w temperaturze 22 °C przez 16 godzin, a następnie w temperaturze 65 °C przez 20 minut, a następnie utrzymuje w temperaturze 4 °C. Powodzenie tego kroku można sprawdzić, ładując 50 ng próbek do 1% żelu agarozowego. Pomyślne podwiązanie spowoduje, że na torze pojawi się pasmo o dużej masie cząsteczkowej (patrz rysunek 2, tor 3).
3. Fragmentacja DNA: W tym momencie powinieneś mieć duże konkamery (>10kb) produktów PCR. Mamy losową strategię sonikacji przy użyciu 24-próbkowego sonikatora Diagenode Bioruptor, który może rozdrobnić te konkamery w ciągu 25 minut (40 sekund "włączenia" / 20 sekund "wyłączenia" na minutę). Sonikacja jest hamowana przez lepkość wprowadzoną przez PEG, więc można to przezwyciężyć poprzez rozcieńczenie próbki w stosunku 10:1 w buforze Qiagen PB. Wyniki można sprawdzić na 2% żelu agarozowym (patrz Rysunek 2, pasy 4 i 5).
4. Próbka jest gotowa do włączenia bezpośrednio do protokołu przygotowania próbki biblioteki Illumina Genomic Library, zaczynając od kroku "Naprawa końcowa". Przedstawione tutaj dane pochodzą z odczytów pojedynczych końcówek w analizatorze genomu Illumina IIx, ale użyliśmy HiSeq 2000 i przeprowadziliśmy odczyty pojedyncze lub sparowane z porównywalnymi wynikami. Biorąc pod uwagę skalę utworzonej biblioteki, użyliśmy również niestandardowych adapterów z kodami kreskowymi w celu multipleksowania wielu bibliotek w puli, aby dostosować się do przepustowości dostarczanej przez platformę HiSeq (dane nie pokazane). Postępuj zgodnie z protokołem producenta i zaleceniami dołączonymi do zestawu. W celu osiągnięcia optymalnej czułości i swoistości do wykrywania wariantów zaleca się 25-krotne lub więcej pokrycie celu na allel (ryc. 3). Oszacowanie to jest niezależne od wielkości puli i rodzaju wariantu, który ma zostać wykryty. W razie potrzeby można połączyć wiele pasów ruchu i biegów, aby uzyskać odpowiedni zasięg.

3. Sekwencjonowanie odczytuje wyrównanie i analizę

1. Kompresja i formatowanie plików: Surowe sekwencjonowane pliki odczytu powinny być przekonwertowane do formatu SCARF lub skompresowane. Kompresja jest opcjonalna, ponieważ oszczędza czas i miejsce na kolejne etapy analizy bez utraty istotnych informacji. Osiąga się to za pomocą dołączonego skryptu RAPGAP_read_compressor_v2.pl za pomocą następującego polecenia:
   ./RAPGAP_read_compressor_v2.pl [Odczyt pliku] > [Skompresowany plik odczytu]
   Akceptowane formaty wejściowe odczytywanych plików to SCARF i FASTQ, spakowane lub nieskompresowane:
   Przykładowy format szalika:
   HWI-EAS440:7:1:0:316#0/1:NTCGATTCACTGCCCAACAACACCAGCTCCTCTCCC:DNWUQSPWWWWUVVPVVVVVVUVVUUPUUWWWWWUW
   Przykładowy format FASTQ:
   @HWI-EAS440_7_1_0_410#0/1
   NGTGGTTTCTTTTTTGGGGGGAGAGGAGCTGGTG
   +
   &/8888888888888888888854588767777666!
2. Wyrównanie odczytu surowego: Odczyty surowe można teraz wyrównać do sekwencji referencyjnej FASTA z adnotacjami specyficznej dla regionów docelowych zawartych w reakcjach PCR, a także kontroli dodatniej i ujemnej. Wyrównanie można wykonać za pomocą dołączonego narzędzia do wyrównywania RAPGAPHASH5d. Format wejściowy w tym momencie musi być SCARF lub skompresowany. Polecenie wyrównania to:
   ./RAPGAPHASH5d [Skompresowany plik do odczytu] [Plik FASTA] [dozwolona liczba edycji] > [Wyrównany plik]
   Liczba niezgodności na odczyt, które są dozwolone w porównaniu z sekwencją odniesienia, jest parametrem zdefiniowanym przez użytkownika. Odczyty, które mają nadmierną liczbę niezgodności, zostaną odrzucone. Zalecamy zezwolenie na 2 niezgodności dla odczytów 36 pz, 4 niezgodności dla odczytów 76 pz i 5 niezgodności dla odczytów 101 pz. Zezwolenie na większą liczbę niezgodności zwiększy prawdopodobieństwo dopuszczenia nadmiaru błędów sekwencjonowania do wyrównanych danych. Ponieważ długość odczytu staje się coraz dłuższa, wartość ta może być jeszcze większa.
3. Znakowanie wyrównanych plików z tej samej komórki przepływu: W tym momencie cały wyrównany plik odczytu powinien otrzymać unikalny identyfikator ("znacznik") w celu zidentyfikowania odczytanych plików należących do tego samego przebiegu sekwencjonowania (tj. wiele linii z tej samej komórki przepływowej można zagregować i nadać jej jeden znacznik). Znacznik jest niezbędny, ponieważ każde uruchomienie maszyny generuje unikatowy profil błędu, który można scharakteryzować za pomocą tagu. Tag to alfanumeryczny ciąg znaków używany do rozróżniania zestawu odczytów (znak podkreślenia "_" nie powinien być używany do analizowania problemów). Różne znaczniki powinny być używane dla wyrównanych plików odczytu generowanych w różnych komórkach przepływu lub przebiegach maszyn. Tagi można dodawać za pomocą dołączonego RAPGAP_alignment_tagger.pl za pomocą następującego polecenia:
   ./RAPGAP_alignment_tagger.pl [Wyrównany plik] [TAG] > [Wyrównany plik oznakowany]
   Po tym momencie wyrównane pliki z tej samej biblioteki wygenerowane na wielu różnych komórkach przepływu mogą być łączone ze sobą, ponieważ ich odpowiednie znaczniki będą je rozdzielać.
4. Generowanie modelu błędów: Jak wspomniano powyżej, każdy przebieg maszyny generuje unikalny profil błędu sekwencjonowania, który należy scharakteryzować w celu dokładnego wywołania wariantu. Aby modelować te błędy dla każdego uruchomienia maszyny, w każdej bibliotece próbek w puli znajduje się wewnętrzna sekwencja kontrolna, o której wiadomo, że jest pozbawiona zmienności sekwencji. Z wyrównanego pliku ze znacznikami można wygenerować plik modelu błędu za pomocą dołączonego narzędzia EMGENERATOR4 z sekwencją odniesienia do kontroli ujemnej. Można użyć całej sekwencji kontroli ujemnej lub alternatywnie tylko jej podzbioru, określonego przez większość zasad 5' i 3' na wejściu. Zawsze należy używać unikalnych odczytów i pseudoliczb:
   ./EMGENERATOR4 [Wyrównany plik oznaczony] [sekwencja kontroli ujemnej] [Nazwa pliku wyjściowego] [5' najbardziej podstawowa kontroli ujemnej, która ma być użyta] [3' najbardziej podstawowa kontroli ujemnej, która ma być użyta] [uwzględnij tylko unikalne odczyty? = Y] [odcięcie edycji wyrównania] [wprowadź pseudoliczby? =Y]
   Narzędzie EMGENERATOR4 wygeneruje 3 pliki o nazwie nazwanej jako parametr nazwy pliku wyjściowego, po którym następuje _0, _1 lub _2. Pliki te odpowiadają odpowiednio modelowi błędów 0., 1. i 2. rzędu. W przypadku wywoływania wariantów za pomocą SPLINTER należy zawsze używać modelu błędu 2. rzędu.
5. W celu wizualizacji profilu współczynnika błędów przebiegu można użyć error_model_tabler_v4.pl do wygenerowania wykresu błędów PDF w pliku modelu błędu 0. rzędu (Rysunek 4):
   ./error_model_tabler_v4.pl [Błąd modelu pliku 0. rzędu] [nazwa pliku wyjściowego]
   Plik wykresu ujawni trendy błędów specyficzne dla przebiegu i może być używany do wywnioskowania maksymalnej liczby baz odczytu, które mają być używane do analizy, co wyjaśniono w następnej sekcji.

4. Wykrywanie rzadkich wariantów za pomocą SPLINTER

1. Wywołanie wariantu przez SPLINTER: Pierwszym krokiem w analizie jest uruchomienie narzędzia SPLINTER na wyrównanym pliku przy użyciu modelu błędów i sekwencji referencyjnej. Polecenie, aby to zrobić, to:
   ./SPLINTER6r [Wyrównany plik otagowany] [Plik FASTA] [plik modelu błędu 2. rzędu] [liczba odczytanych baz do użycia] [odczytaj bazy lub cykle do wykluczenia] [odcięcie wartości p = -1,301] [użyj unikalnych odczytów = Y] [odcięcie edycji wyrównania] [rozmiar puli z dostępnych opcji] [wydrukuj bezwzględne pokrycie na pasmo = Y] > [Plik SPLINTER]
   Liczba baz odczytu, które mają być używane, jest różna i powinna być oceniana zgodnie z każdym przebiegiem. Generalnie zalecamy używanie pierwszych 2/3 odczytu, ponieważ reprezentują one dane najwyższej jakości (na przykład pierwsze 24 podstawy odczytu o długości 36 bp). Pojedyncze bazy odczytu mogą zostać wyłączone z analizy, jeśli zostaną uznane za wadliwe (oddzielone przecinkiem lub N, np. 5,7,11 lub N). Wartość odcięcia wartości p określa, jak rygorystyczna będzie analiza wywołań wariantów. Zwykle rozpoczynamy analizę od zezwolenia na minimalną wartość odcięcia -1,301 (co odpowiada wartości p ≤ 0,05 w skali log10). Opcja wielkości puli optymalizuje dyskryminację algorytmów "sygnał do szumu" poprzez eliminację potencjalnych wariantów z częstością występowania alleli mniejszych niż w przypadku pojedynczego allelu w rzeczywistej puli. Na przykład w puli 50 osobników najniższego obserwowanego wariantu można spodziewać się przy częstości 0,01 lub 1 na 100 alleli. W związku z tym opcja wielkości puli powinna być ustawiona na najbliższą wartość, która jest >większa niż rzeczywista liczba alleli analizowanych w eksperymencie (tj. jeśli badanych jest 40 osób, spodziewamy się 80 alleli, więc najbliższą opcją byłaby wielkość puli 100). Warianty wywoływane na częstotliwościach <0,01 będą następnie ignorowane jako szum. Plik ten zwraca wszystkie trafienia, które są statystycznie istotne w całej próbce, wraz z opisem pozycji wariantu, typem wariantu, wartością p na nić DNA, częstością wariantu i całkowitym pokryciem na nić DNA (tabela 2).
2. Normalizowanie pokrycia dla wywoływanych wariantów: wahania pokrycia w próbce mogą generować fałszywe trafienia. Można to poprawić, stosując skrypt splinter_filter_v3.pl w następujący sposób:
   ./splinter_filter_v3.pl [plik SPLINTER] [plik listy] [ścisłość] > [plik znormalizowany SPLINTER]
   gdzie plik listy jest listą trafień kontroli dodatniej w postaci pliku rozdzielanego tabulatorami.
   Pierwsze pole wskazuje interesujący nas amplikon, podczas gdy drugie pole wskazuje pozycję, w której mutacja jest obecna. N wskazuje, że reszta sekwencji nie zawiera żadnej mutacji.
3. Określenie optymalnych progów wartości p na podstawie danych kontroli dodatniej: Po normalizacji analiza kontroli dodatniej jest niezbędna do maksymalizacji czułości i swoistości analizy danej próbki. Można to osiągnąć, znajdując optymalną wartość odcięcia wartości p na podstawie informacji z kontroli dodatniej. Najprawdopodobniej początkowa wartość p -1,301 nie będzie wystarczająco rygorystyczna, co jeśli tak, spowoduje wywołanie wyników fałszywie dodatnich z kontroli dodatniej lub ujemnej. Każda analiza SPLINTER pokaże rzeczywistą wartość p dla każdego wywoływanego wariantu (patrz kolumny 5 i 6 w tabeli 2), której nie można było przewidzieć a priori. Jednak całą analizę można powtórzyć, stosując najmniej rygorystyczną wartość p wyświetlaną na początkowym wyjściu dla znanych prawdziwie dodatnich pozycji bazowych. Będzie to służyć do zachowania wszystkich wyników prawdziwie dodatnich przy jednoczesnym wykluczeniu większości, jeśli nie wszystkich, wyników fałszywie dodatnich i zazwyczaj mają one znacznie mniej istotne wartości p w porównaniu z wartościami prawdziwie dodatnimi. Aby zautomatyzować ten proces, można użyć cutoff_tester.pl. cutoff_tester.pl wymaga pliku wyjściowego SPLINTER i listy trafień kontroli pozytywnej w postaci pliku rozdzielanego tabulatorami, takiego jak ten używany do normalizacji:
   ./cutoff_tester.pl [Filtrowany plik SPLINTER] [plik listy]
   Wynikowym wynikiem będzie lista punktów odcięcia, które stopniowo osiągają optymalny (patrz tabela 3). Format to :
   [Odległość od maksymalnej czułości i swoistości] [Czułość] [specyfika] [odcięcie]
   Na przykład :
   7.76946294170104e-07 1 0.999118554429264 -16.1019999999967
   Ostatni wiersz reprezentuje najbardziej optymalną wartość odcięcia dla przebiegu i dlatego może być używany do analizy danych. Optymalnym wynikiem jest osiągnięcie czułości i swoistości na poziomie 1. W przypadku, gdy wynik ten nie zostanie osiągnięty, analizę SPLINTER można powtórzyć, zmieniając liczbę włączonych baz odczytu, aż do osiągnięcia najbardziej optymalnych warunków.
4. Końcowe filtrowanie wariantów: Końcowe odcięcie można zastosować do danych za pomocą skryptu cutoff_cut.pl, który przefiltruje plik wyjściowy SPLINTER z trafień poniżej optymalnego odcięcia,
   ./cutoff_cut.pl [filtrowany plik SPLINTER] [odcięcie] > [końcowy plik SPLINTER]
   Ten krok spowoduje wygenerowanie końcowego pliku wyjściowego >SPLINTER, który będzie zawierał SNP i Indels obecne w próbce. Należy pamiętać, że dane wyjściowe dla wstawień są nieco inne niż dla podstawień lub usunięć (Tabela 2).

5. Reprezentatywne wyniki

Zebraliśmy populację 947 osobników i skupiliśmy się na ponad 20 kb do sekwencjonowania. Zastosowaliśmy SPLINTER do wykrywania rzadkich wariantów zgodnie z naszym standardowym protokołem. Każda osoba została wcześniej poddana genotypowaniu za pomocą genotypowania całego genomu. Zgodność między genotypowaniem oznaczonych i nowych wariantów wywołanych w próbie zbiorczej była doskonała (ryc. 6). Trzy warianty, z których dwa (rs3822343 i rs3776110) były rzadkie w populacji, nazwano de novo na podstawie wyników sekwencjonowania i zwalidowano przez indywidualne pirosekwencjonowanie. Częstości występowania alleli drugorzędnych (MAF) w puli były podobne do MAF zgłoszonych w kompilacji DBSNP 129. Zgodność MAF między pirosekwencjonowaniem a sekwencjonowaniem zbiorczym była doskonała (Tabela 3).

Tabela 1. Sekwencje oligonukleotydowe DNA dla kontroli pozytywnej. Każda sekwencja składa się z fragmentu DNA różniącego się od odniesienia typu dzikiego dwoma podstawieniami lub jedną insercją i jedną delecją. Kliknij tutaj, aby zobaczyć większy obraz.

Tabela 2. Przykład wyjścia SPLINTER. Pierwsze dwa wiersze reprezentują standardowe dane wyjściowe SPLINTER do podstawienia lub usunięcia (niebieski nagłówek). Ostatni wiersz reprezentuje standardowe dane wyjściowe SPLINTER dla wstawienia (fioletowy nagłówek). Kliknij tutaj, aby wyświetlić większy obraz.

Tabela 3. Z dużych populacji zidentyfikowano pięć znanych i trzy nowe warianty, które zweryfikowano za pomocą indywidualnego genotypowania. Walidację indywidualną przeprowadzono za pomocą pirosekwencjonowania (rzędy 1-3), testu TaqMan (rzędy 4-6) lub sekwencjonowania Sangera (rzędy 7,8). Dla szerokiego zakresu częstości występowania alleli i w tym pięciu pozycji z MAF <1%, zgodność między estymacją częstości alleli sekwencjonowania zbiorczego a indywidualnym genotypowaniem była silna. Pozycje oznaczone gwiazdką (*) są dostosowane do wcześniej zgłoszonych danych⁹.

Rysunek 1. Omówienie sekwencjonowania zbiorczego DNA i analizy SPLINTER. DNA pacjenta jest łączone i amplifikowane w wybranych loci. Końcowe produkty PCR są łączone razem z kontrolą dodatnią i ujemną w proporcjach równomolowych. Mieszanka w puli jest następnie sekwencjonowana, a wynikowe odczyty są mapowane z powrotem do ich odwołania. Zamapowane odczyty kontroli ujemnej są używane do generowania modelu błędów specyficznego dla przebiegu. SPLINTER może być następnie wykorzystany do wykrywania rzadkich SNP i indeli poprzez włączenie informacji z modelu błędu i kontroli pozytywnej. [Na podstawie Vallania FLM i in., Badania nad genomem 2010] Kliknij tutaj, aby zobaczyć większy obraz.

Rysunek 2. Zbiorcza ligacja amplikonów PCR i sonikacja. Jako demonstracja etapów ligacji i losowej fragmentacji w protokole przygotowania biblioteki, wektor pUC19 został enzymatycznie strawiony do fragmentów pokazanych na torze 2. Fragmenty te zostały znormalizowane przez liczbę cząsteczek, połączone i losowo ligowane zgodnie z krokiem 1.7 powyżej. Powstałe w ten sposób duże konkatamery są pokazane na torze 3. Podwiązane konkatamery zostały równo podzielone i poddane sonikacji, jak opisano w kroku 1.8 powyżej. Otrzymany rozmaz fragmentów DNA dla każdej kontrpróby technicznej pokazano na torach 4 i 5. Nawias podkreśla zakres rozmiarów używany do ekstrakcji żelu i tworzenia biblioteki sekwencjonowania.

Rysunek 3. Dokładność w funkcji pokrycia pojedynczego allelu w próbce zbiorczej. Dokładność jest szacowana jako obszar pod krzywą (AUC) krzywej operatora odbiornika (ROC), która waha się od 0,5 (losowo) do 1,0 (dokładność doskonała). AUC jest wykreślane jako funkcja pokrycia na allel w celu wykrycia pojedynczych zmutowanych alleli w pulach 200, 500 i 1000 alleli (A). AUC jest wykreślane jako całkowite pokrycie funkcji dla podstawień, insercji i delecji (B). [Na podstawie Vallania FLM i in., Genome Research 2010].

Rysunek 4. Wykres błędu pokazuje prawdopodobieństwo uwzględnienia błędnej podstawy w danej pozycji. Profil błędu pokazuje niskie wskaźniki błędów z tendencją wzrostową w kierunku końca 3' odczytu sekwencjonowania. Warto zauważyć, że różne nukleotydy referencyjne wykazują różne prawdopodobieństwa błędu (patrz na przykład prawdopodobieństwo włączenia C z podaniem G jako odniesienia). [Na podstawie Vallania FLM i in., Genome Research 2010].

Rysunek 5. Dokładność SPLINTER w szacowaniu częstości alleli dla pozycji, które miały ponad 25-krotne pokrycie na allel. Na podstawie wyników w Panelu A, Rysunek 3 przedstawiający optymalną czułość do wykrywania pojedynczego wariantu z ≥25-krotnym pokryciem, porównanie częstości alleli puli DNA oszacowanych przez SPLINTER z liczbą alleli mierzoną przez GWAS daje bardzo wysoką korelację (r = 0,999). [Na podstawie Vallania FLM i in., Genome Research 2010].

Rysunek 6. Porównanie częstości występowania alleli mierzonych za pomocą GWAS w porównaniu z szacunkami SPLINTER z sekwencjonowania zbiorczego 974 osób. Dla porównania stwierdzono 19 wspólnych pozycji między genotypowanymi loci a regionami sekwencji. Uzyskana korelacja jest bardzo wysoka (r = 0,99538). Kliknij tutaj, aby zobaczyć większą postać.

Discussion

Istnieje coraz więcej dowodów na to, że częstość występowania i odpowiedź terapeutyczna powszechnych, złożonych fenotypów i chorób, takich jak otyłość8, hipercholesterolemia4, nadciśnienie^tętnicze7 i inne, może być moderowana przez osobiste profile rzadkiej zmienności. Identyfikacja genów i szlaków, w których te warianty agregują się w dotkniętych populacjach, będzie miała głębokie implikacje diagnostyczne i terapeutyczne, ale osobna analiza osób dotkniętych chorobą może być zbyt czasochłonna i kosztowna. Analiza populacyjna oferuje bardziej efektywną metodę badania zmienności genetycznej w wielu loci.

Przedstawiamy nowatorski protokół sekwencjonowania zbiorczego DNA w połączeniu z pakietem oprogramowania SPLINTER zaprojektowanym w celu identyfikacji tego typu zmienności genetycznej w populacjach. Wykazujemy dokładność tej metody w identyfikacji i ilościowym określaniu pomniejszych alleli w dużej puli 947 osobników, w tym rzadkich wariantów, które zostały nazwane de novo z sekwencjonowania zbiorczego i zwalidowane przez indywidualne pirosekwencjonowanie. Nasza strategia różni się od innych protokołów przede wszystkim tym, że w każdym eksperymencie włącza się kontrolę pozytywną i negatywną. Dzięki temu SPLINTER może osiągnąć znacznie wyższą dokładność i moc w porównaniu z innymi podejściami¹. Optymalne pokrycie 25-krotnego allelu jest ustalane niezależnie od wielkości puli, dzięki czemu analiza dużych pul jest wykonalna, ponieważ wymóg ten skaluje się tylko liniowo wraz z wielkością puli. Nasze podejście jest bardzo elastyczne i może być stosowane do każdego fenotypu, ale także do próbek, które są naturalnie niejednorodne, takich jak mieszane populacje komórek i biopsje guza. Biorąc pod uwagę stale rosnące zainteresowanie sekwencjonowaniem zbiorczym z dużych regionów docelowych, takich jak egzom lub genom, nasze przygotowanie biblioteki i analiza SPLINTER są kompatybilne z niestandardowym przechwytywaniem i sekwencjonowaniem całego eksomu, ale narzędzie do wyrównywania w pakiecie SPLINTER nie zostało zaprojektowane dla dużych sekwencji referencyjnych. Dlatego z powodzeniem wykorzystaliśmy wyrównywacz do programowania dynamicznego, Novoalign, do dopasowania całego genomu, a następnie wywołania wariantu z próbki zbiorczej (Ramos i wsp., przedłożone). W ten sposób nasza strategia sekwencjonowania w puli może być z powodzeniem skalowana do większych pul przy rosnącej ilości sekwencji docelowej.

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
PfuUltra High-Fidelity	Agilent	600384	1.4
Betaina	SIGMA	B2629	1.4
M13mp18 wektor ssDNA	NEB	N4040S	1.5
pGEM-T Easy	Promega	A1360	1.5
T4 Kinaza polinukleotydowa	NEB	M0201S	2.2
Ligaza T4	NEB	M0202S	2.2
Glikol polietylenowy 8000 MW	SIGMA	P5413	2.2
Sonikator Bioruptor	Diagenode	UCD-200-TS	2.3