the minor allele is the less frequent allele in the population at a variable site. Najpierw opisujemy dwa główne podejścia do oszacowania mniejszej częstotliwości alleli (MAF)w określonym miejscu w genomie. Pierwsze podejście polega na wnioskowaniu poszczególnych genotypów i traktowaniu tych genotypów jako całkowicie dokładnych przy szacowaniu MAF. Następnie badamy wydajność struktury prawdopodobieństwa, która bezpośrednio bierze pod uwagę niepewność przy przypisywaniu genotypów. W trakcie naszej pracy Zakładamy, że wszystkie strony segregujące są bialleliczne.

oszacowanie MAF na podstawie nazwanych genotypów

jednym ze sposobów oszacowania MAF na podstawie danych sekwencjonowania nowej generacji jest najpierw wywołanie genotypu dla każdej osoby za pomocą danych sekwencjonowania, a następnie użycie tych genotypów tak, jakby były one prawdziwe. Było to podejście tradycyjnie stosowane dla danych genotypowych i danych sekwencjonowania Sangera. Nie jest jasne, jak dobrze będzie on działać, gdy stosowane do danych sekwencjonowania nowej generacji.

podejście maksymalne prawdopodobieństwo może używać wnioskować genotyp dla każdego osobnika od następnej generacji sekwencjonowanie dane. W każdym miejscu j, dla każdego osobnika i, prawdopodobieństwo dla każdego z trzech możliwych genotypów (zakładając, że znamy mniejszy allel) jest podane jako:

(1)

gdzie D i,j to obserwowane dane dotyczące sekwencjonowania u osobnika i w miejscu j, g I , j ∈ {0, 1, 2} to liczba drobnych alleli zawartych w genotypie każdego osobnika, a I kontrola odpowiednio błędów sekwencjonowania i odczytu właściwości bazowych. Obserwowane dane sekwencjonowania dla każdej osoby mogą być traktowane jako wyrównanie odczytów w miejscu j, biorąc pod uwagę wyniki jakości odczytu. Jest to reprezentowane jako prawdopodobieństwo genotypu i znajduje się w pliku prawdopodobieństwa genotypu (GLF), który jest produkowany w wielu programach, które analizują sekwencjonowanie nowej generacji danych, takich jak SOAPsnp i MAQ .

aby przypisać genotyp konkretnej osobie, prawdopodobieństwo każdego z trzech możliwych genotypów można obliczyć dla danej osoby. Genotyp z największym prawdopodobieństwem można następnie przypisać. Jednak badacze często preferują bardziej rygorystyczne kryterium wywołania i nie przypisują genotypu jednostce, chyba że najbardziej prawdopodobny genotyp jest znacznie bardziej prawdopodobny niż drugi najbardziej prawdopodobny. Tutaj trzy możliwe genotypy są posortowane według ich prawdopodobieństwa: , gdzie g(k)odpowiada genotypowi z największym prawdopodobieństwem k. Z podanym progiem f można wywołać genotyp g(1), jeśli . W przeciwnym razie genotyp nie jest wywoływany, a genotyp jednostki jest uważany za brakujący. Wspólna wartość progowa F wynosi 1, co wskazuje, że najbardziej prawdopodobny genotyp jest co najmniej 10 razy bardziej prawdopodobny niż drugi najbardziej prawdopodobny. Zauważ, że tego typu filtrowanie może skutkować większym zaufaniem dla „wywołanego” genotypu, ale również skutkować większą liczbą brakujących danych.

Estymator maksymalnego prawdopodobieństwa częstości alleli

zamiast szacowania MAF z nazwanych genotypów, metoda maksymalnego prawdopodobieństwa (ML) wprowadzona przez Kim i wsp. (patrz również Lynch dla podobnego podejścia) bezpośrednio szacuje MAF i bierze pod uwagę niepewność genotypu. W szczególności, biorąc pod uwagę niewielki allel, prawdopodobieństwo obserwacji danych sekwencji u każdego osobnika i uzyskuje się przez zsumowanie prawdopodobieństwa odpowiadającego wszystkim trzem możliwym genotypom.

przypuśćmy, że dostępne są trzy podobieństwa genotypów określone w równaniu 1. Stosując tę samą notację jak powyżej, niech D j I p j będą odpowiednio obserwowanymi danymi sekwencjonowania w miejscu j i odpowiadającym mu MAF. Prawdopodobieństwo genotypu biorąc pod uwagę, że mniejsza częstotliwość alleli może być obliczona przez założenie równowagi Hardy ’ ego-Weinberga (HWE). Następnie, zakładając niezależność wśród osobników, prawdopodobieństwo MAF w tym locus jest iloczynem wszystkich podobieństw obliczonych dla wszystkich N osobników:

(2)

oszacowanie ML p j można obliczyć poprzez bezpośrednie maksymalizowanie prawdopodobieństwa ograniczonej przestrzeni parametrów za pomocą Broydena-Fletchera-Goldfarba-metoda shanno (bfgs) lub za pomocą algorytmu expectation-maximization (em). Podczas korzystania z algorytmu EM, tylne oczekiwanie genotypu jest obliczane dla każdego osobnika, a średnia tych posterów jest wielokrotnie aktualizowana. Nasza implementacja BFGS była szybsza niż algorytm em. Na przykład, aby uzyskać szacunki z witryn 100,000, BFGS zajęło ~16 sekund, ale EM zajęło ~100 sekund. Jednak różnica w szybkości może być specyficzna dla implementacji. W naszym przypadku, dla obu metod, przestaliśmy aktualizować parametry, gdy wzrost prawdopodobieństwa był mniejszy niż 0.001.

Estymator maksymalnego prawdopodobieństwa z niepewnym drobnym allelem

w praktyce często jako drobny allel można stosować drugi najczęściej występujący nukleotyd u poszczególnych osobników. Jednak w przypadku rzadkich SNP (np. MAF < 1%) trudno jest określić, który allel jest mniejszym allelem, ponieważ wszystkie cztery nukleotydy mogą pojawić się w niektórych odczytach z powodu błędów sekwencjonowania. Aby poradzić sobie z tą sytuacją, opisujemy teraz ramy prawdopodobieństwa, które uwzględniają niepewność przy określaniu podrzędnego allelu.

Załóżmy, że dla strony j znamy główny allel M. Należy zauważyć, że decydowanie, który z dwóch wspólnych alleli może być głównym, nie jest ważne, ponieważ zajmujemy się głównie szacowaniem częstotliwości w rzadkich SNP. Ponadto, w przypadku alleli o częstotliwości pośredniej (około 50%), rozróżnienie między głównym i mniejszym allelem jest mniej ważne. Przydziel pozostałe trzy inne niż główne nukleotydy M1, M2 i m3. Prawdopodobieństwo wprowadzone w równaniu 2 zakłada stały allel główny M i stały allel podrzędny M. dlatego, aby umożliwić niepewność w oznaczeniu drobnego allelu, funkcję prawdopodobieństwa można modyfikować jako:

(3)

ponadto, zakładając, że którykolwiek z trzech możliwych drobnych alleli jest jednakowo prawdopodobny, otrzymujemy:

(4)

gdzie . Ponieważ może być bardzo mała w przypadku dużych zbiorów danych (np. dla wielu osób), przydatne jest obliczenie prawdopodobieństwa w skali logów. Uporządkuj trzy logarytmy warunkowe jako (l(1), l(2), l(3)), gdzie L(1) jest największym. Następnie,

g-test z użyciem nazwanych genotypów do mapowania asocjacji

w badaniach asocjacji SNP wykazujące znaczące różnice w częstości alleli między przypadkami i grupami kontrolnymi uważa się za związane z interesującym fenotypem. Mapowanie asocjacji może być wykonywane przy użyciu danych z badań sekwencjonowania nowej generacji. Najpierw omawiamy podejścia, które wymagają wywołania poszczególnych genotypów, a następnie wykonujemy test na asocjację przy użyciu nazwanych genotypów. W tym podejściu genotyp jest najpierw wywoływany dla każdej osoby. Genotypy mogą być filtrowane lub niefiltrowane. Zakładając niezależność między jednostkami i HWE, można zbudować tabelę awaryjną 2 × 2, licząc liczbę głównych i mniejszych alleli zarówno w przypadkach, jak i w kontrolach. Prowadzi to do dobrze znanego testu współczynnika prawdopodobieństwa dla niezależności, g-test:

(5)

gdzie O k,H jest częstotliwością obserwowaną w komórce, A E k,h jest częstotliwością oczekiwaną pod hipotezą zerową, w której częstotliwość alleli jest taka sama między przypadkami i kontrolami. Dobrze znany Test chi-kwadrat Pearsona jest asymptotycznie równoważny testowi G. Jeśli tabela jest generowana z prawdziwych genotypów, to statystyka G asymptotycznie podąża za rozkładem chi-kwadrat z 1 stopniem swobody (χ2 (1)). Jednak w naszych badaniach konstruujemy g-statystykę za pomocą” zwanych ” genotypów, więc HWE może nie trzymać z powodu nadmiernego i niedostatecznego wywoływania heterozygotów. Ponadto konstruowanie statystyki testu poprzez liczenie” zwanych „genotypów zamiast” obserwowanych ” genotypów prawdopodobnie wprowadza dodatkową zmienność. Dlatego teoria statystyczna może już nie być aktualna. Zauważ, że gdy genotyp nie jest wymagany dla określonej osoby, dane są uważane za brakujące i nie są uwzględnione w tabeli 2 × 2.

Test współczynnika prawdopodobieństwa uwzględniający niepewność w obserwowanych genotypach dla mapowania asocjacji

zamiast wywoływania genotypów, ramy prawdopodobieństwa pozwalają na niepewność w genotypach i testach w każdym miejscu j, czy częstotliwość alleli jest taka sama między przypadkami i kontrolami. Formalnie obliczamy prawdopodobieństwo hipotez H O: p j, 1 = p j, 2 (=p j ,0) I H A : p j ,1 ≠ P j, 2 ,gdzie P j ,1 i p j, 2 to odpowiednio mafy w przypadkach i kontroli.

zakładając, że znane są allele drobne (m) i główne (m), prawdopodobieństwo występowania częstości alleli drobnych można obliczyć w sposób opisany w równaniu 2, a statystykę badania współczynnika prawdopodobieństwa oblicza się jako:

(6)

gdzie I są obserwowanymi danymi odpowiednio dla przypadków i kontroli, oraz I są wartościami mle MAF odpowiednio w przypadkach i kontrolach.

Jeśli mały allel jest nieznany, prawdopodobieństwo pod hipotezą zerową jest obliczana jak w równaniu 3, a statystyka LRT jest modyfikowana jako:

(7)

gdzie D j jest obserwowanymi danymi dla obu przypadków i kontroli, a jest częstotliwością alleli w hipotezie zerowej. Inne zapisy są takie same jak w równaniu 6.

Szacowanie MAF w symulowanych danych

porównujemy szacunki częstości alleli na symulowanych danych za pomocą prawdziwych genotypów (True), zwanych genotypami bez filtrowania (Call NF), nazywanymi genotypami z filtrowaniem (f = 1; Call F) i metodą maksymalnego prawdopodobieństwa (ML). W przypadku rzadkich SNP Typ alleli drobnych często nie jest widoczny. Przy określaniu genotypów przyjmuje się, że drugim najczęściej występującym nukleotydem jest allel mniejszy. Metoda ML bezpośrednio uwzględnia niepewność przy określaniu allelu podrzędnego i o ile nie zaznaczono inaczej, przedstawiono wyniki przy użyciu metody nieznanego allelu podrzędnego (równanie 3). Należy zauważyć, że metoda nieznanego allelu podrzędnego ML działa podobnie do metody znanego allelu podrzędnego ML, ale ta pierwsza jest lepsza dla bardzo rzadkich SNP (dodatkowy plik 1).

najpierw oceniliśmy, jak dobrze różne podejścia były w stanie oszacować MAF U 200 osób w zakresie głębokości sekwencjonowania dla 1000 SNPs z prawdziwym MAF wynoszącym 5%. Rysunek 1 przedstawia boxploty rozkładów szacowanych Maf przy użyciu czterech różnych podejść. Zgodnie z oczekiwaniami, w przypadku danych o większym zasięgu, takich jak indywidualna głębokość 12×, wszystkie metody działają, jak również wtedy, gdy genotypy są znane z pewnością (prawdziwe). Jednak gdy głębokość maleje, szacunki MAF uzyskane przez pierwsze wywołanie genotypów stają się stronnicze. Na przykład mediana MAF oszacowana przy użyciu metody Call F wynosi 5,3% przy pokryciu 6× i 12,5% przy pokryciu 2×. Powodem odchylenia w górę jest to, że trudniej jest nazwać heterozygotami, ponieważ prawdziwe heterozygoty często wyglądają jak błędy sekwencjonowania. Dlatego więcej heterozygotów niż mniejszych homozygotów zwykle ma brakujące genotypy. Jednak ogólne odchylenie w szacunkach MAF od nazwanych genotypów nie zawsze jest w jednym kierunku (dane nie są pokazane). Co ciekawe, błąd wydaje się być gorszy dla metody Call F niż dla metody Call NF. Ten wzorzec może wydawać się intuicyjny, ponieważ filtrowanie wywołań genotypu wydaje się zmniejszać prawdopodobieństwo nazwania błędu sekwencjonowania heterozygotą. Jednak metoda wywołania F powoduje również większą ilość brakujących danych, ponieważ wiele homozygotów dla głównego allelu nie będzie wywoływanych z powodu błędów sekwencjonowania. Tak więc w tym przypadku nazywanie genotypów bez filtrowania wydaje się być lepszą strategią niż filtrowanie genotypów podczas próby oszacowania MAF.

oszacowanie częstotliwości alleli w miejscach o rzeczywistym MAF 5% dla różnych głębokości pokrycia. Na każdej głębokości symulowano 1000 miejsc przy użyciu 200 osobników, a w każdym miejscu oszacowanie częstości alleli oblicza się za pomocą: (1) prawdziwych genotypów (True); (2) zwanych genotypami bez filtrowania (Call nF); (3) zwanych genotypami z filtrowaniem (Call F); i (4) metody maksymalnego prawdopodobieństwa (ML). Więcej szczegółów na temat metod szacowania można znaleźć w Methods.

wyniki są diametralnie różne dla nowej metody ML. Metoda ta zapewnia bezstronne szacunki MAF (mediana ~4,9%) na różnych głębokościach. Nawet przy 2× szacunki pokazują tylko nieco większą wariancję niż te oparte na prawdziwych genotypach.

porównaliśmy również szacowany średni błąd kwadratowy (MSE; oczekiwanie () różnych szacunków MAF w zakresie głębokości sekwencjonowania (Rysunek 2). Metoda ML ma niższy MSE niż metody wywołujące z 50 lub 200 osobnikami. W szczególności MSE obliczane na podstawie metody Call F jest znacznie wyższe niż te z innych metod, zwłaszcza gdy głębokość maleje. MSE szacunków MAF opartych na prawdziwych genotypach odzwierciedla dolną granicę MSE i nie jest stała na głębokościach ze względu na wariancję pobierania próbek i skończoną wielkość próby. Przy użyciu 50 osobników, MSE zbliża się do 0,0005 ze wzrostem głębokości, a przy użyciu próbki o wielkości 200 osobników, zbliża się do 0,0013 ze wzrostem głębokości.

średni błąd squred (MSE; oczekiwany) czterech różnych typów estymatory częstotliwości alleli dla różnych rozmiarów próbek (lewy i prawy panel) i głębokości pokrycia (oś X). Na każdej głębokości MSE obliczono na podstawie oszacowań częstotliwości alleli wykonanych przy użyciu czterech różnych metod: True, Call NF, Call F I ML (szczegóły metod można znaleźć w podpisie na rysunku 1).

ogólnie nowa metoda ML wykonuje metody wywołania genotypu.

Szacowanie rozkładu Maf na podstawie symulowanych danych

następnie badamy, w jaki sposób różne podejścia do szacowania przeprowadzono w szacowaniu odsetka SNP przy różnych częstotliwościach w populacji (podobne do widma częstotliwości miejsca, ale oparte na częstotliwości alleli populacji zamiast częstotliwości próbki). Tutaj symulowaliśmy 20 000 SNPs, gdzie rozkład prawdziwych Maf podążał za standardowym rozkładem stacjonarnym dla efektywnej wielkości populacji wynoszącej 10 000 (patrz metody). Należy jednak zauważyć, że w praktyce bardzo trudno jest odróżnić bardzo rzadki SNP od błędu sekwencjonowania. W związku z tym, w celu porównania z rzeczywistymi danymi, odrzuciliśmy SNP o szacowanym MAF poniżej 2%. Rysunek 3 pokazuje udział SNP wchodzących w każdy inny zakres częstotliwości po wykluczeniu tych SNP z szacowanym MAF<2%.

rozkład częstotliwości alleli SNPs symulowany przy założeniu standardowego stacjonarnego rozkładu częstotliwości alleli. Na każdej głębokości (każdy panel) symulowano 20 000 SNP, a dla każdego SNP oszacowania MAF uzyskano przy użyciu czterech różnych metod (patrz podpis na rysunku 1). Następnie, dla każdej metody (każdego koloru), do wygenerowania każdego histogramu (Oś x) używane są tylko strony o szacowanych częstotliwościach alleli > 2%.

zgodnie z oczekiwaniami, przy dużej głębokości pokrycia, na przykład 10× na osobę, wszystkie metody zapewniają szacunkowe rozkłady MAF, które są podobne do oczekiwanego rozkładu w oparciu o prawdziwe genotypy (rycina 3). Przy płytszej głębokości pokrycia, takiej jak mniej niż 4× na osobnika, rozkłady Maf uzyskane metodami wywoływania genotypu znacznie odbiegają od oczekiwanego rozkładu MAF opartego na prawdziwych genotypach(ryc. 3). W szczególności metody te zawyżają udział SNP o niskiej częstotliwości. Na przykład oczekiwany odsetek SNP w drugim pojemniku (szacowany MAF między 2-4%) wynosi 18%. Odpowiedni odsetek oparty na metodzie Call NF na głębokości 4× wynosi 26%, co jest 1,4-krotnie wyższe niż oczekiwano. Nadmierne oszacowanie proporcji SNP niskiej częstotliwości występuje z powodu mylenia błędów sekwencjonowania z prawdziwymi heterozygotami, co skutkuje nadmiernym obciążeniem heterozygotycznych genotypów. Wielkość tej inflacji różni się w zależności od odcięć filtrujących, ale większe odcięcie niekoniecznie zwiększa lub zmniejsza inflację.

obraz jest zupełnie inny dla metody ML. Szacowany rozkład MAF uzyskany za pomocą nowej metody ML jest ściśle zgodny z rzeczywistym rozkładem nawet przy niewielkich głębokościach pokrycia. Tutaj prawie nie ma nadmiaru SNP o niskiej częstotliwości. Na głębokości 4× odsetek SNPs w drugim pojemniku histogramu wynosi 18,4%, co jest bardzo zbliżone do oczekiwanej proporcji (18%). Tak więc, bardziej wiarygodne szacunki widma częstotliwości mogą być wykonane z danych o niskim zasięgu przy użyciu naszego podejścia prawdopodobieństwa niż przy użyciu podejścia wywołującego genotyp.

mapowanie Asocjacji w symulowanych danych

porównujemy skuteczność metod, które traktują wnioskowane genotypy jako prawdziwe genotypy w testach asocjacji (przy użyciu testu G) z naszym testem prawdopodobieństwa (LRT), który odpowiada za niepewność genotypów. Badamy rozkład testu-statystyki pod hipotezą zerową braku różnicy częstotliwości alleli między przypadkami i kontrolami. Porównujemy również siłę różnych podejść.

przy dość dużych rozmiarach próbek, standardowa teoria asymptotyczna sugeruje, że pod hipotezą zerową zarówno statystyka G, jak i statystyka LRT podążają za rozkładem chi-kwadrat z jednym stopniem swobody (χ2 (1)). Dlatego porównaliśmy rozkład null g-statistic obliczony na podstawie metod wywołujących, jak również statystykę LRT do rozkładu χ2(1) za pomocą wykresów QQ (Rysunek 4). Symulowaliśmy 5000 SNP na różnych głębokościach sekwencjonowania w 500 przypadkach i kontrolach, w których MAF używany do symulacji genotypów wynosił 5% w obu przypadkach i kontrolach. Rozkład g-statystyka obliczona przy użyciu prawdziwych genotypów wykazuje bardzo dobrą zgodność z rozkładem χ2(1). Jednak rozkład g-statystycznej obliczonej na podstawie nazwanych genotypów zasadniczo odbiega od rozkładu χ2 (1). Wywołanie genotypów, a następnie traktowanie tych genotypów jako dokładnych, daje ogromny nadmiar fałszywie dodatnich sygnałów, jeśli wartości p są obliczane przy użyciu rozkładu χ2(1). Na przykład na głębokości 2× 11% SNP miało wartość p mniejszą niż 5%, w porównaniu do oczekiwanych 5%. Efekt jest spowodowany zwiększeniem wariancji, spowodowanym nadmiernym homozygotami jako heterozygotami, w teście allelicznym stosowanym tutaj do wykrywania asocjacji. Badania genotypowe, takie jak Armitage trend test, które są odporne na odchylenia od równowagi Hardy ’ ego-Weinberga, nie wykazują podobnego wzrostu wskaźnika fałszywie dodatniego (dodatkowy plik 2). Zgodnie z tą obserwacją, filtrowanie nazwanych genotypów powoduje zmniejszenie ułamka istotnych testów przy użyciu testu G, chociaż filtrowanie nie rozwiązuje całkowicie problemu. Z drugiej strony, statystyka LRT pokazuje tylko bardzo niewielkie odejście od rozkładu χ2 (1) dla 2× lub 5× głębokości pokrycia.

QQ-wykresy porównujące zerowy rozkład interesującej statystyki badania z rozkładem χ2(1). Każda kolumna odpowiada innej statystyce badania: (1) g-statystyka obliczona przy użyciu prawdziwych genotypów (True); (2) g-statystyka obliczona przy użyciu nazwanych genotypów bez filtrowania (Call NF); (3) g-statystyka obliczona przy użyciu nazwanych genotypów z filtrowaniem (Call F); oraz (4) statystyka badania współczynnika prawdopodobieństwa z nieznanym drobnym allelem (LRT). Zakładając 500 przypadków i 500 kontroli, pod hipotezą zerową, Zestaw 5000 miejsc symulowano z MAF 5% z sekwencjonowania głębokości 2× (górne panele) i 5× (dolne panele). Współczynnik „inflacji” jest pokazany w lewym górnym rogu każdej cyfry.

wygenerowaliśmy również krzywe charakterystyki operacyjnej odbiornika (Roc) dla każdego z różnych testów asocjacji. Krzywe te pokazują moc testu przy różnych wskaźnikach fałszywie dodatnich. Ponieważ rozkłady niektórych statystyk testowych nie podążają za rozkładem χ2 (1) pod hipotezą zerową, aby dokonać uczciwego porównania, otrzymaliśmy wartość krytyczną dla każdego fałszywie dodatniego wskaźnika w oparciu o empiryczny rozkład zerowy. Moc jest obliczana jako ułamek symulowanych loci choroby, które mają statystykę przekraczającą wartość krytyczną. Ogólnie rzecz biorąc, stwierdzamy, że LRT działa lepiej niż g-test na podstawie obu metod wywołania genotypu (Rysunek 5). Na przykład, przy 5% fałszywie dodatnim wskaźniku i przy głębokości sekwencjonowania 5×, moc wykrywania locus choroby z MAF 1% i względnym ryzykiem (RR) 2 wynosi 51% W przypadku LRT, ale moc spada do 33% przy użyciu metody wywołania bez filtrowania i do 34% przy użyciu metody wywołania z filtrowaniem. W szczególności, na małej głębokości, Test G stosowany do nazwanych genotypów z filtrowaniem wypada bardzo słabo (lewa większość kolumny na fig.5). Jeśli porównamy moc LRT z testem trendu Armitage 'a przy użyciu nazwanych genotypów, okaże się, że LRT ma również wyższą moc niż test trendu Armitage’ a (dodatkowy plik 3). Sugeruje to, że jeśli ktoś chce używać nazwanych genotypów, filtrowanie ich na podstawie zaufania do połączenia może spowodować utratę mocy.

krzywe charakterystyki operacyjnej odbiornika (ROC) czterech testów asocjacji. Definicja czterech statystyk znajduje się w podpisie na rysunku 4. Zakładając 500 przypadków i 500 kontroli, zestaw 20 000 miejsc symulowano pod null i pod alternatywą na indywidualnych głębokościach sekwencjonowania 2×, 5× i 10× (trzy kolumny). Przy każdej fałszywie dodatniej szybkości (oś x) odpowiadającą wartość krytyczną obliczano przy użyciu empirycznego rozkładu zerowego. Rzeczywistą wartość dodatnią (moc; oś y) uzyskano przez obliczenie ułamka miejsc sprawczych ze statystykami testowymi, które przekraczają wartość krytyczną.

zastosowanie do rzeczywistych danych

przeanalizowaliśmy 200 egzomów z kontroli do badania związku chorobowego, które zostały zsekwencjonowane przy użyciu technologii Illumina na indywidualnej głębokości 8× . Użyliśmy genotypu likelihoods generowanego przez program „SOAPsnp” do naszego wnioskowania. Aby uzyskać więcej informacji, zobacz metody.

najpierw zbadaliśmy dokładność szacunków MAF z danych sekwencjonowania nowej generacji dla 50 SNPs, porównując je z szacowanymi MAF z danych genotypu Sequenomu. Zarówno szacunki przy użyciu metody ML, jak i metoda wywoływania genotypu bez filtrowania są silnie skorelowane z szacunkami dokonanymi na podstawie danych o genotypie Sequenomu (tj. mała standaryzowana różnica między tymi dwoma szacunkami na fig. 6). Jednak szacunki oparte na wywołaniu genotypu z filtrowaniem wykazują słabą zgodność z częstotliwościami szacowanymi na podstawie danych genotypu Sekwenom, zwłaszcza gdy głębokość sekwencjonowania jest niska. Co ciekawe, istnieje jeden SNP, w którym szacowany MAF z danych resekwencjonowania jest bardzo różny od oszacowania uzyskanego z danych genotypu Sekwenom, mimo że głębokość sekwencjonowania jest bardzo wysoka (14×). W szczególności, szacowany MAF z danych o genotypie Sequenomu wynosi 22,5%, ale 17,2%, gdy szacowany jest przy użyciu metody ML. Badanie Indywidualne pokazuje, że u wielu osób wysoce wspierany genotyp oparty na danych sekwencjonujących różni się od genotypów Sekwenomowych. Biorąc pod uwagę, że ten SNP jest objęty wieloma odczytami u tych osób i że obserwowane bazy odczytu mają wysoką jakość wyników (>Q20), jest prawdopodobne, że różnica wynika z błędów genotypowania sekwencji. Zauważ, że istnieje kilka SNP, w których szacowane MAF z podejścia wywołującego genotyp bez filtrowania wydają się lepiej odpowiadać Maf oszacowanym z Sekwenom genotypowania niż szacunki z podejścia ML. Na przykład przy jednym SNP szacowany MAF wynosi 25,7% z danych dotyczących genotypu Sequenomu, 25.9% od metody wywołania genotypu bez filtrowania i 27,2% od metody ML. Jednak indywidualna kontrola ujawnia, że istnieje kilka osób, dla których zwany genotyp z danych sekwencjonowania różni się od genotypu Sekwenom. W takich przypadkach błędy w nazwanych genotypach anulowane, dając wygląd lepszej korespondencji z danymi o genotypie sekwencji. Dlatego w przypadku tych SNP trudno jest określić, która metoda działa najlepiej.

szacunki częstości alleli obliczone od 200 osób przy użyciu danych sekwencjonowania nowej generacji w porównaniu z danymi genotypu Sekwenomów. W każdym miejscu, tylko osoby, które mają zarówno dane Sekwenom genotypu i dane sekwencjonowania zostały wykorzystane do oszacowania częstości alleli. Dla danych sekwencjonowania, oszacowania MAF uzyskano przy użyciu trzech różnych metod (Call NF; Call F; I ML). Znormalizowaną różnicę dla każdego oszacowania obliczono jako , gdzie I są szacunkowymi wartościami MAF z danych sekwencjonowania i genotypu Sekwenomów, odpowiednio, a n jest liczbą osobników użytych do oszacowania. Każde miejsce jest klasyfikowane do jednego z czterech pojemników na podstawie średniej indywidualnej głębokości pokrycia (koloru): mniejszej niż 4×, wyższej niż 4×, ale mniejszej niż 8×, wyższej niż 8×, ale mniejszej niż 16× i wyższej niż 16×.

następnie zbadaliśmy rozkład Maf obliczony przy użyciu kilku podejść w zakresie głębokości sekwencjonowania z naszych danych sekwencjonowania exome nowej generacji (Rysunek 7). Odrzuciliśmy SNP z szacowanym MAF <2%, ponieważ trudno jest odróżnić te SNP o bardzo niskiej częstotliwości od błędów sekwencjonowania w tym zbiorze danych. Następnie usunęliśmy miejsca, w których wystąpiła znacząca różnica (wartość p mniejsza niż 10-5 za pomocą testu rank-sum ) w wyniku jakości odczytywanych baz pomiędzy allelami podrzędnymi i głównymi. Miejsca te mogą być sztucznymi SNP, które mogą wystąpić z powodu nieprawidłowego mapowania lub nieznanych uprzedzeń wprowadzonych podczas procedury eksperymentalnej. Następnie sklasyfikowaliśmy każdą witrynę w pojemnikach na podstawie głębokości pokrycia. Liczbę SNP w każdym pojemniku przedstawiono w tabeli 1. Gdy średnia głębokość jest mniejsza niż 9×, rozkłady szacunkowych MAF w oparciu o metody wywołania genotypu są bardzo różne od tych opartych na metodzie ML. W szczególności, metody wywołania genotypu powodują duży nadmiar SNP o niskiej częstotliwości (MAF między 2% a 4%). Wzór ten odzwierciedla to, co widzieliśmy w naszych badaniach symulacyjnych (Rysunek 3). Ponadto, dla metod wywoływania genotypu, rozkład częstotliwości alleli zmienia się dramatycznie wraz ze zmianami głębokości sekwencjonowania. Dlatego, jak wspomniano wcześniej, gdy głębokość nie jest bardzo wysoka, metody wywoływania genotypowania prawdopodobnie zawierają wiele fałszywych SNP, które są błędami sekwencjonowania. Błędy te pojawiają się jako nadmiar SNP o niskiej częstotliwości w rozkładzie częstotliwości. Rozkład oparty na metodzie ML jest bardziej stabilny na głębokościach, ale nadal występuje nadmiar SNP o niskiej częstotliwości alleli o głębokości mniejszej niż 9×w porównaniu z odsetkiem SNP o niskiej częstotliwości na większych głębokościach.

rozkład częstości występowania drobnego allelu oszacowanej na podstawie eksomów 200 zsekwencjonowanych osobników. Dla każdego miejsca częstość występowania drobnych alleli oszacowano za pomocą czterech różnych metod: (1) Metoda ML z nieznanym drobnym allelem, (2) metoda ML ze znanym lub stałym drobnym allelem, (3) wywołanie genotypów bez filtrowania (wywołanie NF) i (4) wywołanie genotypów z filtrowaniem (wywołanie F). Każda witryna jest klasyfikowana w pojemnikach na podstawie głębokości pokrycia. Ponadto w każdym histogramie nie uwzględnia się miejsc o szacowanym MAF poniżej 2%. Liczba SNP, które zostały wykorzystane do tej analizy, znajduje się w tabeli 1.

W końcu użyliśmy tych danych resekwencjonowania exome do symulacji badania skojarzenia case-control. Aby zbadać rozkład statystyk testowych asocjacji pod hipotezą zerową, losowo przypisaliśmy 100 osób do grupy przypadków, a pozostałe 100 do grupy kontrolnej. Dla wszystkich SNP na chromosomie 2 z szacunkami MAF > 2% (w oparciu o nieznaną metodę ml), przetestowaliśmy różnice częstotliwości alleli między przypadkami i kontrolami, obliczając statystykę G przy użyciu nazwanych genotypów, zarówno z filtrowaniem, jak i bez, oraz statystykę LRT. Rysunek 8 przedstawia wykresy QQ porównujące rozkłady statystyki testu ze standardowym rozkładem χ2 (1). Jak widać w badaniach symulacyjnych, rozkład zerowy g-statystyki obliczonej przy wywoływaniu genotypów bez filtrowania znacznie odbiega od rozkładu χ2(1). Jednak rozkład zerowy statystyki LRT jest ściśle zgodny z rozkładem χ2(1). Współczynnik inflacji wynosi 1,01, co oznacza, że statystyka LRT działa dobrze, gdy jest stosowana do rzeczywistych danych.

QQ-wykresy porównujące statystyki testu asocjacji dla różnic częstotliwości alleli między 100 przypadkami i 100 kontrolami do rozkładu χ2(1). Fenotypy zostały losowo przypisane indivdiduals w exome resekwencjonowania zestaw danych tak, że istnieje 100 przypadków i 100 kontroli. Dla każdej witryny obliczono trzy statystyki: statystyka G wykorzystująca nazwane genotypy bez filtrowania (Call NF), statystyka g wykorzystująca nazwane genotypy z filtrowaniem (Call F) oraz statystyka LRT. Aby zminimalizować włączenie fałszywych SNP, miejsca z ml MAF szacowane na mniej niż 2% są odrzucane. Do celów wyświetlania pokazywane są wyniki z miejsc na chromosomie 2. Należy pamiętać, że współczynnik inflacji jest pokazany w lewym górnym rogu każdego wykresu QQ.