Wykrywanie rzadkich wariantów genomicznych z sekwencjonowania pulującego przy użyciu SPLINTER

Podsumowanie

Sekwencjonowanie puli DNA jest szybką i kosztowo efektywną strategią wykrywania rzadkich wariantów związanych ze złożonymi fenotypami w dużych kohortach. W niniejszym materiale opisujemy analizę obliczeniową sekwencjonowania nowej generacji pul 32 genów związanych z nowotworami przy użyciu pakietu oprogramowania SPLINTER. Metoda ta jest skalowalna i może być zastosowana do dowolnego analizowanego fenotypu.

Streszczenie

Wraz z gwałtownym postępem technologii sekwencjonowania DNA w ostatnich latach², staje się coraz bardziej oczywiste, że stopień zmienności genetycznej między dowolnymi dwiema osobnikami jest większy, niż wcześniej sądzono³. Z kolei genotypowanie oparte na mikromacierzach nie pozwoliło na zidentyfikowanie istotnego wpływu powszechnych wariantów sekwencji na zmienność fenotypową powszechnych chorób^4,5. Wspólnie obserwacje te doprowadziły do sformułowania hipotezy powszechnych chorób i rzadkich wariantów (Common Disease / Rare Variant), sugerującej, że większość „brakującej dziedziczności” w powszechnych i złożonych fenotypach wynika zamiast tego z indywidualnego profilu rzadkich lub prywatnych wariantów DNA danej osoby^6-8. Jednakże charakterystyka wpływu rzadkiej zmienności na złożone fenotypy wymaga analizy wielu dotkniętych osobników w wielu locusach genomowych, a w idealnym przypadku powinna być porównana z podobnym badaniem w kohorcie osób zdrowych. Pomimo mocy sekwencjonowania oferowanej przez dzisiejsze platformy, przeprowadzenie badania populacyjnego wielu locusów genomowych oraz późniejsza analiza obliczeniowa pozostają dla wielu badaczy barierą nie do przejścia.

Aby zaspokoić tę potrzebę, opracowaliśmy podejście oparte na sekwencjonowaniu puli^1,9 oraz nowy pakiet oprogramowania¹ do bardzo dokładnego wykrywania rzadkich wariantów z uzyskanych danych. Możliwość łączenia genomów z całych populacji dotkniętych osób oraz badania stopnia zmienności genetycznej w wielu celowanych regionach w ramach jednej biblioteki sekwencjonowania zapewnia doskonałe oszczędności kosztów i czasu w porównaniu z tradycyjną metodologią sekwencjonowania pojedynczych próbek. Przy średnim pokryciu sekwencjonowania na allel wynoszącym 25-krotność, nasz autorski algorytm, SPLINTER, wykorzystuje wewnętrzną strategię kontroli wywoływania wariantów do identyfikacji insercji, delecji i substytucji o długości do czterech par zasad z wysoką czułością i swoistością w pulach zawierających do 1 allelu mutanty na 500 osobników. W niniejszej pracy opisujemy metodę przygotowania biblioteki do sekwencjonowania puli, a następnie przedstawiamy instrukcje krok po kroku dotyczące korzystania z pakietu SPLINTER do analizy sekwencjonowania puli (http://www.ibridgenetwork.org/wustl/splinter). Przedstawiamy porównanie sekwencjonowania puli 947 osób, z których każda przeszła również badanie macierzowe całego genomu, przy sekwencjonowaniu ponad 20kb na osobę. Zgodność między genotypowaniem wariantów oznaczonych i nowych wariantów wykrytych w próbce pulowanej była doskonała. Metodę tę można łatwo skalować do dowolnej liczby loci genomowych i dowolnej liczby osobników. Poprzez wprowadzenie wewnętrznych pozytywnych i negatywnych kontroli amplikonów w proporcjach naśladujących badaną populację, algorytm może zostać skalibrowany dla optymalnej wydajności. Strategię tę można również zmodyfikować do zastosowania w przechwytywaniu przez hybrydyzację lub z wykorzystaniem kodów kreskowych specyficznych dla danej osoby, a także zastosować do sekwencjonowania naturalnie heterogenicznych próbek, takich jak DNA nowotworu.

Protokół

Metoda ta została wykorzystana w badaniach opisanych przez Vallanię FML et al. w czasopiśmie Genome Research 2010.

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

1. Połączyć znormalizowaną ilość DNA genomowego od każdej osoby w puli (pulach). Zastosowanie 0,3 ng DNA na osobę w jednej reakcji PCR pozwoli na wprowadzenie około 50 genomów diploidalnych na osobę do każdej reakcji PCR, co zwiększa prawdopodobieństwo jednorodnej amplifikacji dla każdego allelu w puli.
2. Sekwencje genomowe można pobrać z bazy NCBI (http://www.ncbi.nlm.nih.gov/) lub przeglądarki genomów UCSC Genome Browser (http://genome.ucsc.edu/index.html). Podczas pobierania sekwencji należy upewnić się, że użyto funkcji „RepeatMasker” (oznaczonej jako „N”), aby uniknąć zaprojektowania startera w regionie repetytywnym.
3. Do projektowania starterów należy użyć narzędzia internetowego Primer3 (http://frodo.wi.mit.edu/primer3/input.htm), kopiując interesujące regiony genomowe wraz z sekwencjami flankującymi (idealne są zazwyczaj amplikony o długości 600-2000 bp). Optymalne warunki projektowania starterów dla programu Primer 3 są następujące¹⁰: Minimalna długość startera = 19; Optymalna długość startera = 25; Maksymalna długość startera = 30; Minimalna Tm = 64 °C; Optymalna Tm = 70 °C; Maksymalna Tm = 74 °C; Maksymalna różnica Tm = 5 °C; Minimalna zawartość GC = 45; Maksymalna zawartość GC = 80; Liczba zwracanych wyników = 20 (wartość dowolna); Maksymalna stabilność końca 3' = 100. Należy zaprojektować startery do amplifikacji wszystkich interesujących loci genomowych. Po otrzymaniu starterów, liofilizowane zapasy można rozpuścić w 10 mM Tris, pH 7,5 + 0,1 mM EDTA do końcowego stężenia 100 uM, a następnie rozcieńczyć dodatkowo w stosunku 10:1 w ddH₂O do stężenia 10 uM.
4. Amplifikacja PCR: Ze względu na niski wskaźnik błędów (10^-7) i generowanie produktów z tępymi końcami (co jest niezbędne do późniejszego etapu ligacji), zalecamy stosowanie polimerazy DNA o wysokiej wierności (high-fidelity) do amplifikacji dużych amplikonów genomowych. Stosowaliśmy PfuUltra High-Fidelity, jednak 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, po 400 nM każdego startera, 200 μM dNTP, 1x bufor PfuUltra (lub bufor zawierający ≥ 2 mM Mg²⁺ w celu utrzymania wierności enzymatycznej) oraz 5-50 ng połączonego DNA w końcowej objętości 50 μL. Należy zastosować 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 bp / 1,5-3 minuty dla amplikonów 500-1000 bp / 3-5 minut dla amplikonów > 1 kb; 5. Powtórzyć kroki 2-4 przez 25-40 cykli; 6. 65 °C przez 10 minut; 7. inkubacja w 4 °C. W razie potrzeby wyniki PCR można zazwyczaj poprawić poprzez: 1) obniżenie temperatury przyłączania dla małych amplikonów; 2) podwyższenie temperatury przyłączania dla dużych amplikonów; 3. wydłużenie czasu elongacji dla dowolnego amplikonu.
5. Przygotowanie kontroli SPLINTER: Każdy eksperyment SPLINTER wymaga zastosowania kontroli negatywnej i pozytywnej w celu uzyskania optymalnej dokładności. Kontrola negatywna może składać się ze wszystkich pozycji zasad homozygotycznych w dowolnej, zakodowanej paskiem (bar-coded) próbce, 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 niniejszym raporcie kontrolą negatywną jest amplifikowany region o długości 1 934 bp z rdzenia wektora ssDNA M13mp18. Produkt PCR został poddany sekwencjonowaniu metodą Sangera przed użyciem, aby potwierdzić brak zmienności sekwencji w stosunku do materiału źródłowego lub w wyniku amplifikacji PCR. Kontrola pozytywna składa się z panelu wektorów pGEM-T Easy z wstawką o długości 72 bp, w której zaprojektowano specyficzne insercje, delecje i substytucje (Tabela 1). Wektory mieszamy z tłem typu dzikiego w stosunkach molowych tak, aby mutacje występowały z częstotliwością pojedynczego allelu w puli (tzn. dla puli 100 alleli częstotliwość pojedynczego allelu wyniesie 1%). Następnie przeprowadzamy amplifikację PCR mieszanej matrycy kontrolnej, wykorzystując miejsca starterów M13 PUC w pGEM-T Easy, uzyskując końcowy produkt PCR o długości 355 bp.

2. Przygotowanie i sekwencjonowanie połączonej biblioteki PCR

1. Łączenie produktów PCR: Każdy produkt PCR należy oczyścić z nadmiaru starterów. W przypadku oczyszczania na dużą skalę stosowaliśmy oczyszczanie kolumienkowe Qiagen Qiaquick lub płyty filtrujące 96-dołkowe z manifoldem próżniowym. Po oczyszczaniu każdy produkt PCR należy skwantyfikować przy użyciu standardowych technik. Wszystkie produkty PCR (w tym kontrole) należy połączyć w pulę znormalizowaną według liczby cząsteczek, ponieważ łączenie według stężenia spowodowałoby nadreprezentację małych amplikonów w stosunku do większych produktów. Stężenia przelicza się na bezwzględną liczbę cząsteczek DNA na objętość, korzystając ze wzoru: (g / μL) x (1 mol x bp / 660 g) x (1 / # bp w amplikonie) x (6 x 10²³ cząsteczek / 1 mol) = cząsteczek / μL. Następnie wyznaczamy objętość z każdej reakcji wymaganą do połączenia znormalizowanej liczby cząsteczek na amplikon. Liczba ta jest dowolna, może być dostosowana i w rzeczywistości zależy od objętości pipetowanych, które muszą być wystarczająco duże, aby zachować dokładność. Zazwyczaj łączymy 1-2 x 10¹⁰ cząsteczek każdego amplikonu.
2. Ligacja produktów PCR: Krok ten jest niezbędny do uzyskania jednorodnego pokrycia sekwencjonowania, ponieważ sonikacja małych amplikonów PCR spowodowałaby stronniczą reprezentację w stronę ich końców. Aby temu zapobiec, przed fragmentacją ligujemy połączone produkty PCR w duże konkatemery (>= 10 Kb). Polimeraza Pfu Ultra HF generuje tępe końce, co prowadzi do wydajnej ligacji (polimeraza oparta na Taq doda nawis „A” na końcu 3', co uniemożliwi ligację bez uprzedniego dopełnienia lub stępienia końców). Reakcję tę można w razie potrzeby zwiększyć 2-3 krotnie. Reakcja ligacji zawiera 10 U kinazy polinukleotydowej T4, 200 U ligazy T4, 15% w/v polietylenu, 1X bufor ligazy T4, glikol 8000 MW oraz do 2 μg połączonych produktów PCR w końcowej objętości 50 μL. Reakcje inkubuje się w temperaturze 22 °C przez 16 godzin, następnie w 65 °C przez 20 minut, a później przechowuje w 4 °C. Powodzenie tego etapu można sprawdzić, nanosząc 50 ng próbek na 1% żel agarozowy. Pomyślna ligacja objawi się obecnością prążka o wysokiej masie cząsteczkowej w ścieżce (patrz Rycina 2, ścieżka 3).
3. Fragmentacja DNA: Na tym etapie powinny powstać duże konkatemery (>10kb) produktów PCR. Stosujemy strategię losowej sonikacji z wykorzystaniem sonikatora Diagenode Bioruptor na 24 próbki, który może fragmentować te konkatemery w 25 minut (40 s „włączony” / 20 s „wyłączony” na minutę). Sonikacja jest hamowana przez lepkość wprowadzoną przez PEG, co można zniwelować poprzez rozcieńczenie próbki 10:1 w buforze Qiagen PB. Wyniki można sprawdzić na 2% żelu agarozowym (patrz Rycina 2, ścieżki 4 i 5).
4. Próbka jest gotowa do bezpośredniego włączenia do protokołu Illumina Genomic Library Sample Preparation, zaczynając od kroku „End Repair”. Dane przedstawione w niniejszej pracy pochodzą z odczytów single-end na urządzeniu Illumina Genome Analyzer IIx, jednak korzystaliśmy również z HiSeq 2000 i przeprowadzaliśmy odczyty single-end lub paired-end z porównywalnymi wynikami. Ze względu na skalę utworzonej biblioteki, stosowaliśmy również niestandardowe adaptery z kodami kreskowymi w celu multipleksowania wielu połączonych bibliotek, aby wykorzystać przepustowość platformy HiSeq (dane nie pokazane). Należy postępować zgodnie z protokołem i zaleceniami producenta dołączonymi do zestawu. Aby osiągnąć optymalną czułość i swoistość w wykrywaniu wariantów, zaleca się docelowe pokrycie rzędu 25-krotnego lub większego na allel (Rycina 3). Szacunek ten jest niezależny od wielkości puli i rodzaju wykrywanego wariantu. W razie potrzeby można połączyć wiele ścieżek i serii pomiarowych, aby osiągnąć odpowiednie pokrycie.

3. Mapowanie i analiza odczytów sekwencjonowania

1. Kompresja i formatowanie plików: Surowe pliki z odczytami sekwencjonowania powinny zostać skonwertowane do formatu SCARF lub skompresowane. Kompresja jest opcjonalna, ponieważ pozwala zaoszczędzić czas i miejsce w kolejnych etapach analizy bez utraty istotnych informacji. Można to osiągnąć za pomocą dołączonego skryptu RAPGAP_read_compressor_v2.pl, używając następującej komendy:
   ./RAPGAP_read_compressor_v2.pl [Read file] > [Compressed Read file]
   Akceptowane formaty wejściowe plików z odczytami to SCARF oraz FASTQ, zarówno w formie skompresowanej (gzipped), jak i nieskompresowanej:
   Przykład formatu SCARF:
   HWI-EAS440:7:1:0:316#0/1:NTCGATTCACTGCCCAACAACACCAGCTCCTCTCCC:DNWUQSPWWWWUVVPVVWVVVUVVUUPUUWWWWWUW
   Przykład formatu FASTQ:
   @HWI-EAS440_7_1_0_410#0/1
   NGTGGTTTCTTCTTTGGCTGGGGAGAGGAGCTGGTG
   +
   &/8888888888888888888854588767777666!
2. Mapowanie surowych odczytów: Surowe odczyty można teraz zmapować do adnotowanej referencyjnej sekwencji FASTA specyficznej dla regionów docelowych objętych reakcjami PCR, a także do kontroli pozytywnej i negatywnej. Mapowanie można przeprowadzić za pomocą dołączonego narzędzia RAPGAPHASH5d. Format wejściowy na tym etapie musi być SCARF lub skompresowany. Komenda do mapowania to:
   ./RAPGAPHASH5d [Compressed Read file] [FASTA file] [number of edits allowed] > [Aligned file]
   Liczba niedopasowań na odczyt dopuszczalna w porównaniu z sekwencją referencyjną jest parametrem definiowanym przez użytkownika. Odczyty przekraczające tę liczbę niedopasowań zostaną odrzucone. Zalecamy dopuszczenie 2 niedopasowań dla odczytów 36 bp, 4 niedopasowań dla odczytów 76 bp oraz 5 niedopasowań dla odczytów 101 bp. Dopuszczenie większej liczby niedopasowań zwiększy prawdopodobieństwo wprowadzenia nadmiernych błędów sekwencjonowania do zmapowanych danych. W miarę wydłużania się długości odczytów wartość ta może być dalej zwiększana.
3. Tagowanie zmapowanych plików z tej samej komórki przepływowej (flowcell): Na tym etapie cały zmapowany plik z odczytami powinien otrzymać unikalny identyfikator („tag”), aby zidentyfikować pliki z odczytami należące do tego samego przebiegu sekwencjonowania (np. wiele kanałów z tej samej komórki przepływowej może zostać zagregowanych i otrzymać jeden wspólny tag). Tag jest niezbędny, ponieważ każdy przebieg maszyny generuje unikalny profil błędów, który można scharakteryzować za pomocą tagu. Tag jest alfanumerycznym ciągiem znaków służącym do rozróżnienia zestawu odczytów (ze względu na problemy z parsowaniem nie należy używać znaku podkreślenia „_”). Różnych tagów należy używać dla zmapowanych plików z odczytami wygenerowanych na różnych komórkach przepływowych lub w różnych przebiegach maszyny. Tagi można dodać za pomocą dołączonego programu RAPGAP_alignment_tagger.pl, używając następującej komendy:
   ./RAPGAP_alignment_tagger.pl [Aligned file] [TAG] > [Aligned tagged file]
   Po tym kroku zmapowane pliki z tej samej biblioteki, wygenerowane na wielu różnych komórkach przepływowych, mogą zostać połączone, ponieważ ich poszczególne tagi utrzymają je w rozdzieleniu.
4. Generowanie modelu błędów: Jak wspomniano powyżej, każdy przebieg maszyny generuje unikalny profil błędów sekwencjonowania, który musi zostać scharakteryzowany w celu dokładnego wyznaczania wariantów (variant calling). Aby zamodelować te błędy dla każdego przebiegu maszyny, do każdej połączonej biblioteki próbek włącza się wewnętrzną sekwencję kontrolną, o której wiadomo, że jest pozbawiona zmienności sekwencji. Z zmapowanego pliku z tagami można wygenerować plik modelu błędów za pomocą dołączonego narzędzia EMGENERATOR4 oraz referencyjnej sekwencji kontroli negatywnej. Można wykorzystać całą sekwencję kontroli negatywnej lub alternatywnie tylko jej podzbiór, określony przez bazy na końcach 5' i 3' w danych wejściowych. Należy zawsze stosować unikalne odczyty i pseudoliczniki (pseudocounts):
   ./EMGENERATOR4 [Aligned tagged file] [negative control sequence] [Output file name] [5' most base of the negative control to be used] [3' most base of the negative control to be used] [include unique reads only? = Y] [alignment edits cutoff] [enter pseudocounts? =Y]
   Narzędzie EMGENERATOR4 wygeneruje 3 pliki, których nazwy będą składać się z parametru nazwy pliku wyjściowego oraz końcówek _0, _1 lub _2. Pliki te odpowiadają odpowiednio modelom błędów zero-, pierwszego i drugiego rzędu. Do wyznaczania wariantów programem SPLINTER należy zawsze używać modelu błędów drugiego rzędu.
5. Do wizualizacji profilu poziomu błędów przebiegu można użyć programu error_model_tabler_v4.pl , aby wygenerować wykres błędów w formacie PDF na podstawie pliku modelu błędów zero-rzędowego (Rysunek 4):
   ./error_model_tabler_v4.pl [Error model 0th order file] [output file name]
   Plik z wykresem ujawni trendy błędów specyficzne dla danego przebiegu i może zostać wykorzystany do określenia maksymalnej liczby baz odczytu do użycia w analizie, co zostanie wyjaśnione w następnej sekcji.

4. Wykrywanie rzadkich wariantów przy użyciu SPLINTER

1. Wykrywanie wariantów za pomocą SPLINTER: Pierwszym krokiem w analizie jest uruchomienie narzędzia SPLINTER na pliku z dopasowaniem sekwencji przy użyciu modelu błędów i sekwencji referencyjnej. Polecenie służące do tego celu to:
   ./SPLINTER6r [Aligned tagged file] [FASTA file] [2nd order error model file] [number of read bases to be used] [read bases or cycles to be excluded] [p-value cutoff = -1.301] [use unique reads = Y] [alignment edits cutoff] [pool size from the available options] [print out the absolute coverage per strand = Y] > [SPLINTER file]
   Liczba baz odczytu do wykorzystania jest zmienna i powinna być oceniana dla każdego uruchomienia. Zazwyczaj zalecamy wykorzystanie pierwszych 2/3 odczytu, ponieważ reprezentują one dane najwyższej jakości (na przykład pierwsze 24 bazy odczytu o długości 36bp). Pojedyncze bazy odczytu mogą zostać wykluczone z analizy, jeśli zostaną uznane za wadliwe (oddzielone przecinkiem lub N, np. 5,7,11 lub N). Wartość progowa p-value określa rygorystyczność analizy wykrywania wariantów. Zazwyczaj rozpoczynamy analizę, dopuszczając minimalną wartość progową -1.301 (co odpowiada p-value ≤ 0.05 w skali log10). Opcja wielkości puli (pool size) optymalizuje dyskryminację „sygnału od szumu” algorytmu poprzez eliminację potencjalnych wariantów o częstotliwości mniejszego allela niższej niż w przypadku pojedynczego allelu w rzeczywistej puli. Na przykład w puli 50 osobników można oczekiwać, że najniższy zaobserwowany wariant pojawi się z częstotliwością 0.01 lub 1 na 100 alleli. Zatem opcja wielkości puli powinna być ustawiona na najbliższą wartość, która jest większa niż rzeczywista liczba alleli analizowanych w eksperymencie (tzn. jeśli badanych jest 40 osób, oczekujemy 80 alleli, więc najbliższą opcją będzie wielkość puli 100). Warianty wykryte z częstotliwością <0.01 zostaną wówczas zignorowane jako szum. Plik ten zwraca wszystkie trafienia, które są statystycznie istotne w całej próbce, wraz z opisem pozycji wariantu, typu wariantu, p-value dla każdej nici DNA, częstotliwości wariantu i całkowitego pokrycia dla każdej nici DNA (Tabela 2).
2. Normalizacja pokrycia dla wykrytych wariantów: Fluktuacje pokrycia w obrębie próbki mogą generować fałszywe trafienia. Można to skorygować, stosując skrypt splinter_filter_v3.pl w następujący sposób:
   ./splinter_filter_v3.pl [SPLINTER file] [list file] [stringency] > [SPLINTER normalized file]
   gdzie plik listy jest listą trafień kontroli pozytywnej w formie pliku rozdzielonego tabulatorami.
   Pierwsze pole wskazuje interesujący nas amplikon, natomiast drugie pole wskazuje pozycję, w której występuje mutacja. N oznacza, że reszta sekwencji nie zawiera żadnej mutacji.
3. Wyznaczanie optymalnych progów p-value przy użyciu danych kontroli pozytywnej: Po normalizacji analiza kontroli pozytywnej jest niezbędna do maksymalizacji czułości i swoistości analizy konkretnej próbki. Można to osiągnąć poprzez znalezienie optymalnej wartości progowej p-value przy użyciu informacji z kontroli pozytywnej. Najprawdopodobniej początkowa wartość p-value -1.301 nie będzie wystarczająco rygorystyczna, co w takim przypadku doprowadzi do wykrycia wyników fałszywie dodatnich w kontroli pozytywnej lub negatywnej. Każda analiza SPLINTER wykaże rzeczywistą wartość p-value dla każdego wykrytego wariantu (patrz kolumny 5 i 6 w Tabeli 2), czego nie można było przewidzieć a priori. Jednak całą analizę można powtórzyć, stosując najmniej rygorystyczną wartość p-value wyświetloną w początkowym wyniku dla znanych prawdziwie dodatnich pozycji baz. Posłuży to do zachowania wszystkich prawdziwie dodatnich wyników przy jednoczesnym wykluczeniu większości, jeśli nie wszystkich, wyników fałszywie dodatnich, które zazwyczaj mają znacznie mniej istotne wartości p-value w porównaniu do prawdziwie dodatnich. Aby zautomatyzować ten proces, można użyć cutoff_tester.pl. Skrypt cutoff_tester.pl wymaga pliku wyjściowego SPLINTER oraz listy trafień kontroli pozytywnej w formie pliku rozdzielonego tabulatorami, takiego jak ten użyty do normalizacji:
   ./cutoff_tester.pl [SPLINTER filtered file] [list file]
   Wynikiem będzie lista wartości progowych, które stopniowo prowadzą do wartości optymalnej (patrz Tabela 3). Format to:
   [distance from max sensitivity and specificity] [sensitivity] [specificity] [cutoff]
   na przykład:
   7.76946294170104e-07 1 0.999118554429264 -16.1019999999967
   Ostatnia linia reprezentuje najbardziej optymalną wartość progową dla danego uruchomienia i może być zatem wykorzystana do analizy danych. Optymalnym wynikiem jest osiągnięcie czułości i swoistości na poziomie 1. W przypadku nieosiągnięcia tego wyniku 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: Ostateczną wartość progową można zastosować do danych za pomocą skryptu cutoff_cut.pl, który odfiltruje z pliku wyjściowego SPLINTER trafienia poniżej optymalnego progu,
   ./cutoff_cut.pl [SPLINTER filtered file] [cutoff] > [SPLINTER final file]
   Krok ten wygeneruje końcowy plik wyjściowy SPLINTER, który będzie zawierał SNP i Indele obecne w próbce. Należy pamiętać, że wynik dla insercji różni się nieco od wyników dla substytucji lub delecji (Tabela 2).

5. Reprezentatywne wyniki

Zgrupowaliśmy populację 947 osobników i skierowaliśmy do sekwencjonowania obszar o wielkości ponad 20 kb. Do wykrywania rzadkich wariantów zastosowaliśmy metodę SPLINTER zgodnie z naszym standardowym protokołem. Wcześniej u każdego osobnika przeprowadzono genotypowanie przy użyciu macierzy do genotypowania całego genomu. Zgodność między genotypowaniem wariantów znacznikowych a nowo wykrytymi wariantami w próbce zgrupowanej była doskonała (Rycina 6). Trzy warianty, z których dwa (rs3822343 i rs3776110) były rzadkie w populacji, zidentyfikowano de novo na podstawie wyników sekwencjonowania i potwierdzono poprzez indywidualne pirosekwencjonowanie. Częstotliwości rzadkiego allelu (MAF) w puli były zbliżone do MAF raportowanych w dbSNP build 129. Zgodność MAF między pirosekwencjonowaniem a sekwencjonowaniem zgrupowanym była doskonała (Tabela 3).

Tabela 1. Sekwencje oligonukleotydów DNA dla kontroli pozytywnej. Każda sekwencja składa się z fragmentu DNA różniącego się od referencyjnego typu dzikiego (Wild Type) dwiema substytucjami lub jedną insercją i jedną delecją. Kliknij tutaj, aby wyświetlić większy obraz.

Tabela 2. Przykład wyników z programu SPLINTER. Dwa pierwsze wiersze przedstawiają standardowe wyniki SPLINTER dla substytucji lub delecji (niebieski nagłówek). Ostatni wiersz przedstawia standardowe wyniki SPLINTER dla insercji (fioletowy nagłówek). Kliknij tutaj, aby wyświetlić powiększony obraz.

Tabela 3. Zidentyfikowano pięć znanych i trzy nowe warianty w dużych populacjach, które następnie zwalidowano za pomocą indywidualnego genotypowania. Walidacja indywidualna została przeprowadzona metodą pirosekwencjonowania (wiersze 1-3), testem TaqMan (wiersze 4-6) lub sekwencjonowaniem metodą Sangera (wiersze 7, 8). Dla szerokiego zakresu częstości alleli, w tym dla pięciu pozycji z MAF <1%, stwierdzono silną zgodność między szacowaną częstością alleli z sekwencjonowania pulowego a genotypowaniem indywidualnym. Pozycje oznaczone gwiazdką (*) zostały zaadaptowane z wcześniej raportowanych danych⁹.

Rycina 1. Przegląd sekwencjonowania połączonego DNA i analizy SPLINTER. DNA pacjentów jest łączone i amplifikowane w wybranych loci. Końcowe produkty PCR są łączone wraz z kontrolą pozytywną i negatywną w stosunkach ekimolarnych. Połączona mieszanina jest następnie sekwencjonowana, a otrzymane odczyty są mapowane do ich referencji. Zmapowane odczyty kontroli negatywnej służą do wygenerowania modelu błędów specyficznego dla danego przebiegu. SPLINTER może być następnie wykorzystany do wykrywania rzadkich SNP i indeli poprzez włączenie informacji z modelu błędów i kontroli pozytywnej. [Zaadaptowano z Vallania FLM et al, Genome Research 2010] Kliknij tutaj, aby wyświetlić większy obraz.

Rycina 2. Ligacja i sonikacja połączonych amplikonów PCR. W ramach demonstracji etapów ligacji i losowej fragmentacji w protokole przygotowania biblioteki, wektor pUC19 został poddany trawieniu enzymatycznemu do fragmentów widocznych w ścieżce 2. Fragmenty te znormalizowano pod kątem liczby cząsteczek, połączono i poddano losowej ligacji zgodnie z krokiem 1.7 opisanym powyżej. Powstałe duże konkatemery widoczne są w ścieżce 3. Zligowane konkatemery podzielono na równe części i poddano sonikacji zgodnie z opisem w kroku 1.8 powyżej. Powstałe smugi fragmentów DNA dla każdego powtórzenia technicznego przedstawiono w ścieżkach 4 i 5. Nawias zaznacza zakres wielkości wykorzystany do ekstrakcji z żelu i tworzenia biblioteki sekwencjonowania.

Rycina 3. Dokładność jako funkcja pokrycia dla pojedynczego allelu w próbce pulowanej. Dokładność szacowana jest jako pole pod krzywą (AUC) charakterystyki operacyjnej odbiornika (ROC), którego zakres wynosi od 0,5 (przypadkowa) do 1,0 (dokładność idealna). AUC przedstawiono jako funkcję pokrycia na allel dla detekcji pojedynczych alleli zmutowanych w pulach 200, 500 i 1000 alleli (A). AUC przedstawiono jako funkcję całkowitego pokrycia dla substytucji, insercji i delecji (B). [Zaadaptowano z Vallania FLM et al, Genome Research 2010].

Rycina 4. Wykres błędów przedstawia prawdopodobieństwo wbudowania błędnej zasady w danej pozycji. Profil błędów wykazuje niskie stopy błędów z tendencją wzrostową w kierunku końca 3' odczytu sekwencjonowania. Co istotne, różne nukleotydy referencyjne wykazują różne prawdopodobieństwa wystąpienia błędu (patrz np. prawdopodobieństwo wbudowania C przy G jako referencji). [Zadaptowano z Vallania FLM et al, Genome Research 2010].

Rysunek 5. Dokładność metody SPLINTER w szacowaniu częstości alleli dla pozycji o pokryciu większym niż 25-krotne na allel. W oparciu o wyniki przedstawione w Panelu A na Rysunku 3, wykazujące optymalną czułość wykrywania pojedynczego wariantu przy pokryciu ≥25-krotnym, porównanie częstości alleli w puli DNA oszacowanych przez SPLINTER z liczbą alleli zmierzoną za pomocą GWAS wykazuje bardzo wysoką korelację (r = 0,999). [Zaadaptowano z Vallania FLM et al, Genome Research 2010].

Rycina 6. Porównanie częstości alleli zmierzonych za pomocą GWAS z szacunkami SPLINTER z sekwencjonowania pulowanego 974 osobników. Do porównania wybrano 19 wspólnych pozycji między loci genotypowanych a regionami sekwencjonowanymi. Uzyskana korelacja jest bardzo wysoka (r = 0,99538). Kliknij tutaj, aby wyświetlić powiększoną rycinę.

Dyskusja

Coraz więcej dowodów wskazuje na to, że częstość występowania i odpowiedź terapeutyczna w przypadku powszechnych, złożonych fenotypów i chorób, takich jak otyłość⁸, hipercholesterolemia⁴, nadciśnienie tętnicze⁷ i inne, mogą być modyfikowane przez osobiste profile rzadkich wariantów. Identyfikacja genów i szlaków, w których warianty te gromadzą się w populacjach dotkniętych tymi schorzeniami, będzie miała istotne znaczenie diagnostyczne i terapeutyczne, jednak osobna analiza chorych osób może być zbyt czasochłonna i kosztowna. Analiza oparta na populacji oferuje bardziej efektywną metodę badania zmienności genetycznej w wielu loci.

Przedstawiamy nowatorski protokół sekwencjonowania zpulowanych próbek DNA w połączeniu z pakietem oprogramowania SPLINTER, zaprojektowany do identyfikacji tego typu wariancji genetycznej w populacjach. Demonstrujemy dokładność tej metody w identyfikacji i kwantyfikacji rzadkich alleli w dużej zpulowanej populacji 947 osobników, w tym wariantów rzadkich, które zostały zidentyfikowane de novo z sekwencjonowania puli i zwalidowane za pomocą indywidualnego pirosekwencjonowania. Nasza strategia różni się od innych protokołów głównie włączeniem kontroli pozytywnej i negatywnej do każdego eksperymentu. Pozwala to programowi SPLINTER osiągnąć znacznie wyższą dokładność i moc w porównaniu do innych podejść¹. Optymalne pokrycie 25-krotne na allel jest stałe, niezależnie od wielkości puli, co sprawia, że analiza dużych pul jest wykonalna, ponieważ wymóg ten skaluje się liniowo wraz z wielkością puli. Nasze podejście jest bardzo elastyczne i może być stosowane do każdego interesującego fenotypu, a także do próbek naturalnie heterogennych, takich jak mieszane populacje komórek i biopsje guzów. Biorąc pod uwagę rosnące zainteresowanie sekwencjonowaniem pul z dużych regionów docelowych, takich jak eksom lub genom, nasze przygotowanie bibliotek i analiza SPLINTER są kompatybilne z sekwencjonowaniem typu custom-capture oraz sekwencjonowaniem całego eksomu, jednak narzędzie do dopasowywania w pakiecie SPLINTER nie zostało zaprojektowane dla dużych sekwencji referencyjnych. W związku z tym z powodzeniem wykorzystaliśmy program do dopasowywania oparty na programowaniu dynamicznym, Novoalign, do dopasowań w całym genomie, a następnie do identyfikacji wariantów z zpulowanej próbki (Ramos et al., w trakcie recenzji). Zatem nasza strategia sekwencjonowania pul może być z powodzeniem skalowana do większych pul z coraz większą ilością sekwencji docelowej.

Materiały

Lista materiałów użytych w tym artykule

Nazwa	Firma	Numer katalogowy	Komentarze
PfuUltra High-Fidelity	Agilent	600384	1.4
Betaina	SIGMA	B2629	1.4
Wektor ssDNA M13mp18	NEB	N4040S	1.5
pGEM-T Easy	Promega	A1360	1.5
Polinukleotydowa kinaza T4	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