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 Vallania FML et al. w Genome Research 2010.

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

1. Połączyć znormalizowaną ilość DNA genomowego od każdej osoby w danej puli (puliach). 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, co zwiększa prawdopodobieństwo jednorodnej amplifikacji 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). Przy pobieraniu sekwencji należy upewnić się, że użyto narzędzia „RepeatMasker” (oznaczonego jako „N”), aby uniknąć projektowania startera w regionie repetytywnym.
3. Do projektowania starterów należy wykorzystać narzędzie webowe Primer3 (http://frodo.wi.mit.edu/primer3/input.htm), kopiując i wklejając interesujące regiony genomowe wraz z sekwencjami flankującymi (zazwyczaj idealne są 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. 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 stężenia końcowego 100 uM, a następnie przeprowadzić dodatkowe rozcieńczenie 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) oraz generowanie produktów o tępych końcach (co jest niezbędne w kolejnym etapie ligacji), do amplifikacji dużych amplikonów genomowych zaleca się stosowanie polimerazy DNA o wysokiej wierności. Użyto polimerazy PfuUltra High-Fidelity, jednak enzymy o podobnych charakterystykach (takie jak Phusion) powinny zapewnić porównywalne wyniki. Każda reakcja PCR zawiera końcowe stężenia: 2,5 U polimerazy PfuUltra High-Fidelity, 1 M betainy, 400 nM każdego startera, 200 μM dNTP, 1x bufor PfuUltra (lub bufor zawierający ≥ 2 mM Mg²⁺ w celu zachowania wierności enzymatycznej) oraz 5-50 ng DNA z puli 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. Przechowywanie 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 próbce z 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 niniejszym raporcie kontrolą negatywną jest amplifikowany region o długości 1 934 bp z rdzenia wektora ssDNA M13mp18. Produkt PCR został zsekwencjonowany metodą Sangera przed użyciem, aby potwierdzić brak wariacji sekwencji względem materiału wyjściowego lub wynikającej z amplifikacji PCR. Kontrola pozytywna składa się z panelu wektorów pGEM-T Easy z klonowanym wstawkiem o długości 72 bp, w którym 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 zmieszanego szablonu kontrolnego przy użyciu miejsc wiązania starterów M13 PUC w pGEM-T Easy, generują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. Użyliśmy oczyszczania kolumienkowego Qiagen Qiaquick lub płytek filtrujących 96-dołkowych z manifoldem próżniowym do oczyszczania na dużą skalę. Po oczyszczeniu każdy produkt PCR należy poddać ilościowemu oznaczeniu przy użyciu standardowych technik. Wszystkie produkty PCR (łącznie z kontrolami) 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 względem większych produktów. Stężenia przelicza się na bezwzględną liczbę cząsteczek DNA w jednostce objętości, 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 zależy w rzeczywistości 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 sekwencjonowaniem, ponieważ sonikacja małych amplikonów PCR spowodowałaby stronniczą reprezentację w kierunku ich końców. Aby temu zapobiec, łączymy pulę produktów PCR w duże konkatemery (>= 10 Kb) przed fragmentacją. Polimeraza Pfu Ultra HF generuje tępe końce, co prowadzi do wydajnej ligacji (polimeraza oparta na Taq doda zwis „A” na końcu 3', co uniemożliwi ligację bez wcześniejszego dopełnienia lub stępienia końców). Reakcję tę można zwiększyć 2-3 krotnie, jeśli jest to konieczne. Reakcja ligacji zawiera 10 U kinazy polinukleotydowej T4, 200 U ligazy T4, 15% w/v glikolu polietylenowego (PEG), 1X bufor ligazy T4, glikol 8000 MW, 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ć, nakładając 50 ng próbek na 1% żel agarozowy. Prawidłowa ligacja objawi się obecnością prążka o wysokiej masie cząsteczkowej w odpowiedniej ścieżce (patrz Rycina 2, ścieżka 3).
3. Fragmentacja DNA: Na tym etapie powinny być obecne duże konkatemery (>10kb) produktów PCR. Stosujemy strategię randomowej sonikacji przy użyciu sonikatora Diagenode Bioruptor na 24 próbki, który może sfragmentować 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 przygotowania biblioteki genomicznej Illumina, rozpoczynając od etapu „End Repair” (naprawa końców). Przedstawione tutaj dane pochodzą z odczytów jednokierunkowych (single-end) wykonanych na analizatorze Illumina Genome Analyzer IIx, jednak korzystaliśmy również z HiSeq 2000 i przeprowadzaliśmy odczyty jednokierunkowe lub parzyste (paired-end) z porównywalnymi wynikami. Biorąc pod uwagę skalę stworzonej biblioteki, stosowaliśmy również dedykowane adaptery z kodami kreskowymi, aby zmultipleksować wiele połączonych bibliotek w celu wykorzystania przepustowości platformy HiSeq (dane nie pokazano). Należy postępować zgodnie z protokołem i zaleceniami producenta dołączonymi do zestawu. Aby uzyskać optymalną czułość i swoistość wykrywania 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. Dopasowanie i analiza odczytów sekwencjonowania

1. Kompresja i formatowanie plików: Surowe pliki z odczytami sekwencjonowania powinny zostać przekonwertowane do formatu SCARF lub skompresowane. Kompresja jest opcjonalna, ponieważ oszczędza czas i miejsce w kolejnych etapach analizy bez utraty istotnych informacji. Osiąga się to za pomocą dołączonego skryptu RAPGAP_read_compressor_v2.pl i następującego polecenia:
   ./RAPGAP_read_compressor_v2.pl [Plik z odczytami] > [Skompresowany plik z odczytami]
   Akceptowane formaty wejściowe plików z odczytami to SCARF oraz FASTQ, w wersji skompresowanej (gzipped) lub 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 sekwencji referencyjnej FASTA specyficznej dla regionów docelowych uwzględnionych w reakcjach PCR, a także do kontroli pozytywnej i negatywnej. Mapowanie można przeprowadzić za pomocą dołączonego narzędzia do wyrównywania RAPGAPHASH5d. Format wejściowy na tym etapie musi być SCARF lub skompresowany. Polecenie mapowania brzmi:
   ./RAPGAPHASH5d [Skompresowany plik z odczytami] [Plik FASTA] [liczba dopuszczalnych edycji] > [Plik zmapowany]
   Liczba niedopasowań (mismatches) na odczyt dopuszczalna w stosunku do sekwencji referencyjnej jest parametrem definiowanym przez użytkownika. Odczyty przekraczające tę liczbę niedopasowań zostaną odrzucone. Zaleca się 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 nadmiarowych błędów sekwencjonowania do zmapowanych danych. Wraz ze wzrostem długości odczytów wartość ta może być odpowiednio 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 należące do tego samego cyklu sekwencjonowania (np. wiele linii z tej samej komórki przepływowej może zostać zagregowanych i otrzymać jeden wspólny tag). Tag jest niezbędny, ponieważ każdy cykl pracy urządzenia 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 odróżnienia zestawu odczytów (ze względu na problemy z parsowaniem nie należy używać znaku podkreślenia „_”). Różne tagi powinny być stosowane dla zmapowanych plików z odczytami wygenerowanych na różnych komórkach przepływowych lub w różnych cyklach pracy urządzenia. Tagi można dodać za pomocą dołączonego skryptu RAPGAP_alignment_tagger.pl i następującego polecenia:
   ./RAPGAP_alignment_tagger.pl [Plik zmapowany] [TAG] > [Zmapowany plik otagowany]
   Po tym kroku zmapowane pliki z tej samej biblioteki, wygenerowane na wielu różnych komórkach przepływowych, mogą zostać połączone, ponieważ ich odpowiednie tagi pozwolą zachować ich rozdzielność.
4. Generowanie modelu błędów: Jak wspomniano powyżej, każdy cykl pracy urządzenia generuje unikalny profil błędów sekwencjonowania, który musi zostać scharakteryzowany w celu dokładnego wywoływania wariantów (variant calling). Aby zamodelować te błędy dla każdego cyklu, do każdej puli bibliotek próbek włącza się wewnętrzną sekwencję kontrolną, o której wiadomo, że jest pozbawiona wariacji sekwencji. Z zmapowanego pliku otagowanego można wygenerować plik modelu błędów za pomocą dołączonego narzędzia EMGENERATOR4 wraz z referencyjną sekwencją kontroli negatywnej. Można wykorzystać całą sekwencję kontroli negatywnej lub alternatywnie tylko jej podzbiór, określając najbardziej 5' i 3' zasady w wejściu. Należy zawsze stosować unikalne odczyty i pseudoliczniki (pseudocounts):
   ./EMGENERATOR4 [Zmapowany plik otagowany] [sekwencja kontroli negatywnej] [Nazwa pliku wyjściowego] [najbardziej 5' zasada kontroli negatywnej do użycia] [najbardziej 3' zasada kontroli negatywnej do użycia] [uwzględnić tylko unikalne odczyty? = Y] [próg edycji mapowania] [wprowadzić pseudoliczniki? =Y]
   Narzędzie EMGENERATOR4 wygeneruje 3 pliki o nazwie określonej w parametrze pliku wyjściowego z dodatkiem końcówek _0, _1 lub _2. Pliki te odpowiadają odpowiednio modelowi błędów 0., 1. i 2. rzędu. Do wywoływania wariantów za pomocą SPLINTER należy zawsze stosować model błędów 2. rzędu.
5. W celu wizualizacji profilu poziomu błędów danego cyklu można użyć skryptu error_model_tabler_v4.pl do wygenerowania wykresu błędów w formacie PDF na podstawie pliku modelu błędów 0. rzędu (Rysunek 4):
   ./error_model_tabler_v4.pl [Plik modelu błędów 0. rzędu] [nazwa pliku wyjściowego]
   Plik z wykresem ujawni specyficzne dla danego cyklu trendy błędów i może służyć do określenia maksymalnej liczby zasad odczytu do wykorzystania w analizie, co zostanie omówione w następnej sekcji.

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

1. Wykrywanie wariantów za pomocą SPLINTER: Pierwszym krokiem analizy jest uruchomienie narzędzia SPLINTER na pliku z dopasowanymi sekwencjami, przy użyciu modelu błędów i sekwencji referencyjnej. Komenda służąca 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 wykorzystywanych zasad odczytu jest zmienna i powinna zostać oceniona dla każdego uruchomienia. Zazwyczaj zalecamy wykorzystanie pierwszych 2/3 odczytu, ponieważ reprezentują one dane najwyższej jakości (na przykład pierwsze 24 zasady odczytu w odczycie o długości 36bp). Pojedyncze zasady odczytu mogą zostać wykluczone z analizy, jeśli zostaną uznane za wadliwe (oddzielone przecinkiem lub N, np. 5,7,11 lub N). Próg p-value określa rygorystyczność analizy wykrywania wariantów. Zazwyczaj rozpoczynamy analizę, przyjmując minimalny próg -1.301 (co odpowiada p-value ≤ 0.05 w skali log10). Opcja wielkości puli (pool size) optymalizuje dyskryminację sygnału od szumu w algorytmie poprzez eliminację potencjalnych wariantów o częstotliwości mniejszego allelu niższej niż częstotliwość pojedynczego allelu w rzeczywistej puli. Na przykład w puli 50 osobników najniższy obserwowany wariant może wystąpić z częstotliwością 0.01 lub 1 na 100 alleli. Zatem opcja wielkości puli powinna zostać 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, spodziewamy się 80 alleli, więc najbliższą opcją będzie wielkość puli 100). Warianty wykryte z częstotliwościami <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 oraz 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 (list file) jest listą trafień kontroli pozytywnej w formie pliku rozdzielanego 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 danej próbki. Można to osiągnąć poprzez znalezienie optymalnego progu p-value z wykorzystaniem informacji z kontroli pozytywnej. Najprawdopodobniej początkowe p-value wynoszące -1.301 nie będzie wystarczająco rygorystyczne, co może skutkować wykryciem wyników fałszywie dodatnich w kontroli pozytywnej lub negatywnej. Każda analiza SPLINTER wykaże rzeczywiste 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że całą analizę można powtórzyć, stosując najmniej rygorystyczne p-value wyświetlone w początkowym wyniku dla znanych prawdziwie dodatnich pozycji zasad. Służy to zachowaniu wszystkich wyników prawdziwie dodatnich 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 wyników prawdziwie dodatnich. Aby zautomatyzować ten proces, można użyć cutoff_tester.pl. cutoff_tester.pl wymaga pliku wyjściowego SPLINTER oraz listy trafień kontroli pozytywnej w formie pliku rozdzielanego tabulatorami, takiego jak ten użyty do normalizacji:
   ./cutoff_tester.pl [SPLINTER filtered file] [list file]
   Wynikiem będzie lista progów, 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 optymalny próg dla danego uruchomienia i może zatem zostać wykorzystana 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 zasad odczytu, aż do uzyskania najbardziej optymalnych warunków.
4. Końcowe filtrowanie wariantów: Ostateczny próg 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]
   Ten krok wygeneruje końcowy plik wyjściowy SPLINTER, który będzie zawierał SNP i Indele obecne w próbce. Należy pamiętać, że format wyniku dla insercji jest nieco inny niż dla substytucji lub delecji (Tabela 2).

5. Reprezentatywne wyniki

Zgrupowaliśmy populację 947 osób i wytypowaliśmy ponad 20 kb do sekwencjonowania. Zastosowaliśmy metodę SPLINTER do wykrywania rzadkich wariantów zgodnie z naszym standardowym protokołem. Każda osoba została wcześniej poddana genotypowaniu za pomocą macierzy genomowej. Zgodność między genotypowaniem wariantów znacznikowych a nowymi wariantami zidentyfikowanymi w próbce zgrupowanej była doskonała (Rycina 6). Trzy warianty, z których dwa (rs3822343 i rs3776110) były rzadkie w populacji, zostały zidentyfikowane de novo na podstawie wyników sekwencjonowania i zwalidowane poprzez indywidualne pirosekwencjonowanie. Częstości alleli mniejszościowych (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 referencji typu dzikiego (Wild Type) dwiema substytucjami lub jedną insercją i jedną delecją. Kliknij tutaj, aby wyświetlić powiększony obraz.

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

Tabela 3. Z dużych populacji zidentyfikowano pięć znanych i trzy nowe warianty, które następnie zwalidowano poprzez genotypowanie indywidualne. Walidację indywidualną przeprowadzono za pomocą pirosekwencjonowania (wiersze 1-3), testu TaqMan (wiersze 4-6) lub sekwencjonowania metodą Sangera (wiersze 7, 8). Dla szerokiego zakresu częstości alleli, w tym dla pięciu pozycji z MAF <1%, stwierdzono wysoką zgodność między szacowaniem częstości alleli z sekwencjonowania zagregowanego a genotypowaniem indywidualnym. Pozycje oznaczone gwiazdką (*) zostały zaadaptowane z wcześniej opublikowanych danych⁹.

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

Rycina 2. Ligacja i sonikacja połączonych amplikonów PCR. W ramach demonstracji etapów ligacji i randomowej fragmentacji w protokole przygotowania biblioteki, wektor pUC19 został poddany trawieniu enzymatycznemu do fragmentów widocznych w ścieżce 2. Fragmenty te znormalizowano pod względem liczby cząsteczek, połączono i poddano randomowej ligacji zgodnie z krokiem 1.7 opisanym powyżej. Powstałe duże konkatamery przedstawiono w ścieżce 3. Zligowane konkatamery 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 sekwencyjnej.

Rycina 3. Dokładność w funkcji pokrycia dla pojedynczego allelu w próbce połączonej. Dokładność szacowana jest jako pole pod krzywą (AUC) charakterystyki operacyjnej odbiornika (ROC), gdzie wartości wahają się od 0,5 (przypadek losowy) do 1,0 (dokładność idealna). AUC przedstawiono w funkcji pokrycia na allel dla detekcji pojedynczych alleli zmutowanych w pulach zawierających 200, 500 i 1000 alleli (A). AUC przedstawiono w funkcji całkowitego pokrycia dla substytucji, insercji i delecji (B). [Zadaptowano z Vallania FLM et al, Genome Research 2010].

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

Rycina 5. Dokładność metody SPLINTER w szacowaniu częstotliwości alleli dla pozycji o pokryciu przekraczającym 25-krotność na allel. W oparciu o wyniki przedstawione w Panelu A na Rycinie 3, wykazujące optymalną czułość wykrywania pojedynczych wariantów przy pokryciu ≥25-krotnym, porównanie częstotliwości alleli w puli DNA oszacowanych przez SPLINTER z liczbą alleli zmierzoną w GWAS wykazuje bardzo wysoką korelację (r = 0,999). [Adaptowano z Vallania FLM et al, Genome Research 2010].

Rycina 6. Porównanie częstości alleli zmierzonych za pomocą GWAS z szacunkami SPLINTER uzyskanymi z sekwencjonowania pulowego 974 osobników. Do porównania wybrano 19 wspólnych pozycji pomiędzy locusy genotypowane a regionami sekwencjonowanymi. Wynikowa 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