minor alela je méně časté alely v populaci na proměnnou stránky. Nejprve popíšeme dva hlavní přístupy k odhadu frekvence malých alel (MAF) v určitém místě v genomu. První přístup zahrnuje odvození jednotlivých genotypů a léčbu těchto odvozených genotypů jako zcela přesných při odhadu MAF. Poté zkoumáme výkon rámce pravděpodobnosti, který přímo bere v úvahu nejistotu při přiřazování genotypů. Během naší práce, předpokládáme, že všechna segregační místa jsou biallelic.

Odhad MAF z tzv. genotypy

Jeden způsob, jak odhadnout MAF z next-generation sekvenování dat je nejprve zavolat genotyp pro každé jednotlivé použití sekvenčních dat, a pak použít tyto genotypy, jako kdyby oni jsou ty skutečné. To byl přístup, který se tradičně používá pro data genotypu a data Sangerova sekvenování. Není jasné, jak dobře bude fungovat při aplikaci na sekvenční data nové generace.

přístup s maximální pravděpodobností lze použít k odvození genotypu pro každého jednotlivce z údajů o sekvenování nové generace. Na každém místě j, pro každého jednotlivce i, je pravděpodobnost pro každý ze tří možných genotypů (za předpokladu, že známe menší alelu) dána jako:

kde D i,j je pozorován sekvenčních dat v jednotlivých i v místě j, g, i , j ∈ {0, 1, 2} je číslo menší alely obsažené v genotypu každého jedince, a ovládání pro sekvenování chyby a přečtěte si základní vlastnosti, respektive. Pozorovaná data sekvenování pro každého jednotlivce lze považovat za zarovnání čtení na místě j s přihlédnutím ke skóre kvality čtení. Toto je reprezentováno jako pravděpodobnost genotypu a Nachází se v souboru pravděpodobnosti genotypu (GLF), který je produkován v mnoha programech, které analyzují sekvenční data nové generace, jako jsou SOAPsnp a MAQ.

pro přiřazení genotypu konkrétnímu jednotlivci lze vypočítat pravděpodobnost každého ze tří možných genotypů pro jednotlivce. Genotyp s nejvyšší pravděpodobností pak může být přiřazen. Nicméně, vědci často dávají přednost přísnějšímu kritériu volání a nepřidělí genotyp jednotlivci, pokud nejpravděpodobnější genotyp není podstatně pravděpodobnější než druhý nejpravděpodobnější. Zde jsou tři možné genotypy seřazeny podle jejich pravděpodobnosti: , kde g(k)odpovídá genotypu s největší pravděpodobností k. S daným prahem f lze nazvat genotyp g (1), pokud . V opačném případě není genotyp nazýván a genotyp jedince je považován za chybějící. Společná prahová hodnota f je 1, což naznačuje, že nejpravděpodobnější genotyp je nejméně 10krát pravděpodobnější než druhý nejpravděpodobnější. Všimněte si, že tento typ filtrování může mít za následek vyšší důvěru pro“ nazývaný “ genotyp,ale také to má za následek více chybějících dat.

Maximální věrohodnosti odhad frekvenci alel

Namísto odhadování MAF z tzv. genotypů, maximální pravděpodobnost (ML) metoda zavedená Kim et al. (viz také Lynch pro podobný přístup) přímo odhaduje MAF a bere v úvahu nejistotu genotypu. Konkrétně, vzhledem k tomu, minor alela, pravděpodobnost pozorování posloupnosti dat na každé jednotlivé já se získá sečtením přes pravděpodobností odpovídající na všechny tři možné genotypy.

Předpokládejme, že jsou k dispozici tři pravděpodobnosti genotypu definované v rovnici 1. Při použití stejného zápisu jako výše nechť D j a p j jsou pozorovaná data sekvenování v místě j a odpovídající MAF. Pravděpodobnost genotypu vzhledem k tomu, že frekvence menší alely může být vypočtena za předpokladu Hardy-Weinbergovy rovnováhy (HWE). Pak, za předpokladu nezávislosti mezi jednotlivci, pravděpodobnost MAF na tomto lokusu je produktem všech pravděpodobností vypočítanou přes všech N jedinců:

ML odhad p j lze vypočítat buď přímo maximalizace pravděpodobnosti pro omezený parametr prostoru pomocí Broyden-Fletcher-Goldfarb-Shanno (BFGS) metoda, nebo pomocí expectation-maximization (EM) algoritmus . Při použití algoritmu EM se pro každého jednotlivce vypočítá zadní očekávání genotypu a průměr těchto posteriorů se opakovaně aktualizuje. Naše implementace BFGS byla rychlejší než algoritmus EM. Například pro získání odhadů ze 100 000 stránek bfgs trvalo ~16 sekund, ale EM trvalo ~100 sekund. Rozdíl v rychlosti však může být specifický pro implementaci. V našem případě jsme pro obě metody zastavili aktualizaci parametrů, když zvýšení pravděpodobnosti bylo menší než 0,001.

maximální odhad pravděpodobnosti s nejistou malou alelou

v praxi lze často jako menší alelu použít druhý nejběžnější nukleotid napříč jednotlivci. MAF < 1%) je však těžké určit, která alela je menší alela, protože všechny čtyři nukleotidy se mohou objevit v některých čteních kvůli chybám sekvenování. Abychom se vypořádali s touto situací, nyní popisujeme rámec pravděpodobnosti, který bere v úvahu nejistotu při určování menší alely.

Předpokládejme, že pro místo j známe hlavní alelu m. Všimněte si, že rozhodování o tom, která ze dvou společných alel bude pravděpodobně hlavní, není důležité, protože se většinou zabýváme odhadem frekvencí u vzácných SNP. Dále je pro alely s mezilehlými frekvencemi (kolem 50%) rozdíl mezi hlavní a vedlejší alelou méně důležitý. Přiřaďte další tři nemajetkové nukleotidy m1, m2 a m3. Pravděpodobnost zavedená v rovnici 2 předpokládá pevnou hlavní alelu M a pevnou menší alelu m. Proto, umožnit nejistotu při určování menší alely, funkci pravděpodobnosti lze upravit jako:

Dále, za předpokladu, že každá ze tří možných menší alely je stejně pravděpodobný, dostáváme:

kde . Jelikož může být velmi malý, s velkými soubory dat (např., s mnoha jednotlivci), to je užitečné pro výpočet pravděpodobnosti v log měřítku. Objednejte si tři podmíněné log-pravděpodobností jako k (l(1), l(2), jsem(3)), kde l(1) je největší. Pak,

G-test pomocí tzv. genotypů pro asociační mapování

V asociační studie, Snp ukazuje významné rozdíly ve frekvenci alel mezi případy a kontroly jsou řekl, aby byl spojen s fenotypem zájmu. Asociační mapování lze provést pomocí dat ze sekvenačních studií nové generace. Nejprve diskutujeme přístupy, které vyžadují volání jednotlivých genotypů, a poté provedeme test asociace pomocí tzv. V tomto přístupu je pro každého jednotlivce nejprve povolán genotyp. Genotypy mohou být filtrovány nebo nefiltrovány. Za předpokladu nezávislosti mezi jednotlivci a HWE, kontingenční tabulku 2 × 2 lze sestavit počítáním počtu hlavních a vedlejších alel v případech i kontrolách. To vede k dobře známo, že pravděpodobnost, že poměr test pro nezávislost, G-test:

O, k,h je frekvence pozorovány v buňce, a E k,h je frekvence očekává, že za nulové hypotézy, v nichž se alely frekvence je stejná mezi případy a kontrolami. Známý Pearsonův chí-kvadrát test je asymptoticky ekvivalentní g-testu. Pokud je tabulka generována ze skutečných genotypů, pak g-statistika asymptoticky následuje rozdělení chí-kvadrát s 1 stupněm volnosti(χ2 (1)). Nicméně, v našich studiích, konstruujeme g-statistiku pomocí“ nazývaných “ genotypů, takže HWE nemusí držet kvůli nadměrnému a nedostatečnému volání heterozygotů. Dále, konstrukce statistiky testu počítáním“ nazývaných „genotypů namísto“ pozorovaných “ genotypů pravděpodobně přináší další variabilitu. Statistická teorie tedy již nemusí být platná. Všimněte si, že pokud není genotyp povolán pro určitého jedince, data se považují za chybějící a nejsou zahrnuta v tabulce 2 × 2.

Pravděpodobnost, že poměr test účetnictví pro nejistotu v pozorované genotypy pro asociační mapování

Místo volání genotypů, pravděpodobnost rámec umožňuje pro nejistotu v genotypy a testy na každém místě j zda alela frekvence je stejná mezi případy a kontrolami. Formálně vypočítáme pravděpodobnost hypotéz H O: p j, 1 = p j, 2 (= p j, 0 )a H A : p j, 1 ≠ p j, 2, kde p j, 1 a p j, 2 jsou MAF v případech a kontrolách.

za předpokladu, že jsou známy menší (m) a hlavní (M) alely, lze vypočítat pravděpodobnost frekvence menší alely, jak je popsáno v rovnici 2, a statistika testu poměru pravděpodobnosti se vypočítá jako:

kde jsou pozorovány data pro případy a kontroly, respektive, a jsou MLEs z MAFs v případů a kontrol, resp.

není-li menší alela známa, je pravděpodobnost podle nulové hypotézy vypočtena jako v rovnici 3 a statistika LRT je upravena jako:

kde D j je pozorovaná data pro oba případy a kontroly, a je alela frekvence pod nulovou hypotézu. Ostatní zápisy jsou stejné jako v rovnici 6.

Odhad MAF v simulované údaje

porovnat odhady frekvenci alel na simulovaných dat pomocí true genotypů (Pravda), tzv. genotypy bez jakékoliv filtrace (Call NF), tzv. genotypů s filtraci (f = 1; F) a metoda maximální věrohodnosti (ML). U vzácných SNP není typ menší alely často zřejmý. Při volání genotypů se předpokládá, že druhým nejčastějším nukleotidem je menší alela. Metoda ML přímo zahrnuje nejistotu při určování minor alely a pokud není uvedeno jinak, jsou uvedeny výsledky pomocí metody neznámé minor alely (rovnice 3). Všimněte si, že neznámá metoda ml minor alely funguje podobně jako známá metoda ml minor alely, ale první z nich je lepší pro velmi vzácné SNP (další soubor 1).

nejprve jsme vyhodnotili, jak dobře byly různé přístupy schopny odhadnout MAF u 200 jedinců v rozsahu hloubek sekvenování pro 1 000 SNP se skutečným MAF 5%. Figura 1 ukazuje boxplots distribucí odhadovaných MAF pomocí čtyř různých přístupů. Jak se očekávalo, pro vyšší údaje o pokrytí, jako je individuální hloubka 12×, fungují všechny metody stejně dobře, jako když jsou genotypy s jistotou známy (pravdivé). Když se však hloubka sníží, odhady MAF získané prvním voláním genotypů se stanou předpojatými. Například medián MAF odhadovaný metodou volání F je 5,3% při 6× pokrytí a 12,5% při 2×. Důvodem vzestupného zkreslení je to, že je těžší volat heterozygoty, protože skutečné heterozygoty často vypadají jako chyby sekvenování. Proto více heterozygotů než menších homozygotů má tendenci mít chybějící genotypy. Celkové zkreslení odhadů MAF z tzv. genotypů však není vždy v jednom směru (údaje nejsou zobrazeny). Zajímavé je, že zkreslení se zdá být horší pro metodu volání F než metodu volání NF. Tento vzor se může zdát pult-intuitivní, protože filtrování genotypu hovory by se zdát, snížit pravděpodobnost volání sekvenování chyba heterozygot. Metoda volání F však také vede k většímu množství chybějících dat, protože mnoho homozygotů pro hlavní alelu nebude voláno kvůli chybám sekvenování. Tím pádem, v tomto případě, volání genotypů bez filtrování se zdá být lepší strategií než filtrování genotypů při pokusu o odhad MAF.

výsledky se u nové metody ML dramaticky liší. Tato metoda poskytuje nestranné odhady MAF (medián ~4,9%) v celém rozsahu hloubek. I při 2×, odhady ukazují jen o něco větší rozptyl než odhady založené na skutečných genotypech.

porovnávali jsme také odhadovanou střední kvadratickou chybu (MSE; Expectation () různých odhadů MAF v celém rozsahu hloubek sekvenování (Obrázek 2). Metoda ML má nižší MSE než metody volání s 50 nebo 200 jedinci. Zejména, MSE vypočtená na základě metody volání F je mnohem vyšší než u jiných metod, zejména když hloubka klesá. MSE odhadů MAF na základě skutečných genotypů odráží spodní hranici MSE a není konstantní v hloubkách kvůli rozptylu odběru vzorků a konečné velikosti vzorku. Pomocí 50ti jedinců, MSE přístupy 0.0005 s rostoucí hloubkou a při použití vzorku velikosti 200 osob, přístupy 0.0013 s rostoucí hloubkou.

celkově nová metoda ML provádí metody volání genotypu.

Odhad distribuce MAFs od simulované údaje

dále Jsme zkoumat, jak různé odhadu přístupy provádí v odhadování podílu Snp na různé frekvence v populaci (podobné stránky frekvenční spektrum, ale na základě populační frekvenci alel místo vzorkovací frekvence). Tady jsme simulovali 20,000 Snp, kde distribuce pravda MAFs následoval standardní stacionární rozdělení pro efektivní velikost populace 10 000 (viz Metody). Všimněte si, že v praxi je však velmi obtížné rozlišit velmi vzácný SNP od chyby sekvenování. Proto jsme pro srovnání s reálnými daty vyřadili SNP s odhadovaným MAF méně než 2%. Obrázek 3 ukazuje podíl SNP spadajících do každého zásobníku s různou frekvencí po vyloučení těchto SNP s odhadovaným MAF<2%.

Jak se očekávalo, s vysokou hloubku pokrytí, jako jsou 10× za jednotlivé, všechny metody poskytují odhaduje MAF distribucí, které jsou podobné očekává, že distribuce založená na skutečné genotypů (Obrázek 3). S mělčí hloubkou pokrytí, například menší než 4× na jednotlivce, se distribuce MAF získaných metodami volání genotypů významně odchylují od očekávané distribuce MAF založené na skutečných genotypech (obrázek 3). Zejména tyto metody nadměrně odhadují podíl nízkofrekvenčních SNP. Například očekávaný podíl SNP ve druhém zásobníku (odhadovaný MAF mezi 2-4%) je 18%. Odpovídající podíl založený na metodě Call NF v hloubce 4× je 26%, což je 1,4 krát vyšší, než se očekávalo. Přes-odhad podílu low-frekvence Snp se vyskytuje v důsledku zmatku sekvenování chyby s pravou heterozygoti, což má za následek overcalling heterozygotní genotypy. Velikost této inflace se liší v různých filtrování v kraťáskách, ale větší cutoff nemusí nutně zvýšit nebo snížit inflaci.

obrázek je u metody ML zcela odlišný. Odhadovaná distribuce MAF získaná z nové metody ML úzce sleduje skutečnou distribuci i při mělkých hloubkách pokrytí. Zde není téměř žádný přebytek nízkofrekvenčních SNP. V hloubce 4× je podíl SNP ve druhém zásobníku histogramu 18,4%, což je velmi blízko očekávanému podílu (18%). Tím pádem, spolehlivější odhady frekvenčního spektra lze provést z dat s nízkým pokrytím pomocí našeho přístupu pravděpodobnosti než pomocí přístupů volání genotypu.

Sdružení mapování v simulované údaje

porovnat výkonnost metody, které léčí odvodit genotypy jako pravda genotypů v testech asociace (pomocí G-test), aby naše pravděpodobnost ratio test (LRT), že účty pro nejistotu v genotypů. Zkoumáme distribuci test-statistiky podle nulové hypotézy o rozdílu frekvence alel mezi případy a kontrolami. Porovnáváme také sílu různých přístupů.

S přiměřeně velký vzorek velikosti, standardní asymptotické teorie naznačuje, že za nulové hypotézy jak G-statistika a LRT statistika podle chi-kvadrát rozdělení s jedním stupněm volnosti (χ2(1)). Proto jsme porovnali nulovou distribuci g-statistiky vypočtenou na základě metod volání a statistiky LRT s distribucí χ2 (1) pomocí QQ-grafů (obrázek 4). Simulovali jsme 5 000 SNP v různých hloubkách sekvenování v 500 případech a kontrolách, kde MAF použitý k simulaci genotypů byl 5% v obou případech a kontrolách. Distribuce g-statistiky vypočtená pomocí skutečných genotypů vykazuje velmi dobrou korespondenci s distribucí χ2 (1). Distribuce g-statistiky vypočtená na základě nazývaných genotypů se však podstatně odchyluje od distribuce χ2(1). Volání genotypů a pak léčení těchto genotypů jako je přesné produkuje obrovský přebytek falešně pozitivní signály, pokud p-hodnoty jsou počítány pomocí χ2(1) distribuce. Například v hloubce 2× mělo 11% SNP hodnotu p menší než 5% ve srovnání s očekávanými 5%. Tento efekt je způsoben zvýšením rozptylu, vzhledem k overcalling homozygotů jako heterozygoti, v allelickou test zde použit pro detekci sdružení. Genotypové testy, jako je Armitage trend zkoušky, které jsou robustní na odchylky od Hardy-Weinberg equilibrium, nevykazují podobné zvýšení falešné pozitivity (Další soubor 2). V souladu s tímto pozorováním vede filtrování nazývaných genotypů ke snížení podílu významných testů při použití G-testu, i když filtrování problém zcela nevyřeší. Na druhé straně statistika LRT ukazuje pouze velmi mírný odklon od distribuce χ2 (1) pro 2× nebo 5× hloubky pokrytí.

Jsme také generovány receiver operating characteristic (ROC) křivky pro každou z různých asociací. Tyto křivky ukazují sílu testu při různých falešně pozitivních rychlostech. Protože distribuce některých statistik testů nesledují distribuci χ2 (1) podle nulové hypotézy, pro spravedlivé srovnání jsme získali kritickou hodnotu pro každou falešně pozitivní míru na základě empirického nulového rozdělení. Výkon se počítá jako zlomek simulovaných lokusů onemocnění, které mají statistiku překračující kritickou hodnotu. Celkově zjistíme, že LRT funguje lépe než g-test založený na obou metodách volání genotypu (obrázek 5). Například, při 5% falešně pozitivním záchytu a s sekvenční hloubky 5×, moc odhalit onemocnění, locus s MAF 1% a relativní riziko (RR) 2 je 51% s LRT, ale výkon klesne na 33% pomocí volání metody bez filtrace a 34% pomocí volání metody s filtrací. Zejména v malé hloubce funguje G-test aplikovaný na tzv. genotypy s filtrací velmi špatně (vlevo většina sloupců na obrázku 5). Pokud porovnáme sílu LRT s Armitage trend testem pomocí tzv. genotypů, zjistíme, že LRT má také vyšší výkon než Armitage trend test (další soubor 3). To naznačuje, že pokud si člověk přeje používat tzv. genotypy, jejich filtrování na základě důvěry volání může mít za následek ztrátu energie.

Aplikace na reálná data

analyzovali Jsme 200 exomes z kontroly onemocnění, asociační studie, které byly sekvenovány pomocí Illumina technologie na jednotlivé-jednotlivé hloubky 8× . Pro náš závěr jsme použili pravděpodobnosti genotypu generované programem „SOAPsnp“. Další podrobnosti viz metody.

Nejprve jsme zkoumali přesnost odhadů MAF z dat sekvenování nové generace pro 50 SNP jejich porovnáním s odhadovanými MAF z údajů o genotypu Sekvenomu. Oba odhady pomocí ML způsob a genotypu volá metodu bez filtrování jsou vysoce korelovány s odhady z Sequenom genotyp dat (tj. malé standardizovaný rozdíl mezi oběma odhady v Obrázku 6). Nicméně, odhady založené na volání genotypu s filtrováním ukazují špatnou korespondenci s frekvencemi odhadovanými z údajů genotypu Sekvenomu, zvláště když je hloubka sekvenování nízká. Zajímavé je, že existuje jeden SNP, kde se odhadovaný MAF z údajů o resekvencování velmi liší od odhadu získaného z údajů o genotypu Sekvenomu, i když hloubka sekvenování je velmi vysoká (14×). Konkrétně, odhad MAF z Sequenom genotyp dat je 22.5%, ale je 17.2% při odhadované pomocí ML přístupu. Individuální vyšetření ukazuje, že u mnoha jedinců se vysoce podporovaný genotyp založený na sekvenačních datech liší od genotypů Sekvenomu. Vzhledem k tomu, že tento SNP je pokryt mnoha čte v těchto jednotlivců a že pozorované číst základny mají vysoké skóre kvality (>Q20), je pravděpodobné, že rozdíl je způsoben Sequenom genotypizace chyby. Všimněte si, že existuje několik Snp v jejichž odhadovaná MAFs z genotypu volá přístup, bez filtrování zdá lépe odpovídat MAFs odhadnout z Sequenom genotypizace než odhady z ML přístupu. Například u jednoho SNP je odhadovaný MAF 25,7% z údajů o genotypu Sequenomu, 25.9% z metody volání genotypu bez filtrování a 27,2% z metody ML. Nicméně, individuální kontrola odhaluje, že existuje několik jedinců, u nichž se nazývaný genotyp ze sekvenačních dat liší od genotypu Sequenomu. V těchto případech, chyby v tzv. genotypů zrušen, dávat vzhled lepší korespondence s Sequenom genotyp data. Proto je pro tyto SNP těžké říci, která metoda funguje nejlépe.

dále Jsme zkoumali distribuci MAFs, které se vypočítají pomocí několika přístupů v celé řadě sekvenování hloubkách od naší příští generace exome sequencing data (Obrázek 7). Jsme zlikvidovat Snp se odhaduje MAF <2%, protože je obtížné rozlišit tyto velmi nízkou frekvencí Snp od sekvenování chyby v tomto souboru. Dále jsme odstranili místa, ve kterých byl významný rozdíl (hodnota p menší než 10-5 pomocí testu rank-sum) V skóre kvality čtecích základen mezi menšími a hlavními alelami. Tyto stránky jsou pravděpodobně umělé SNP, ke kterým může dojít v důsledku nesprávného mapování nebo neznámých předsudků zavedených během experimentálního postupu. Pak jsme klasifikovali každé místo do popelnic na základě hloubky pokrytí. Počet SNP v každém zásobníku je uveden v tabulce 1. Pokud je průměrná hloubka menší než 9×, distribuce odhadovaných MAF založených na metodách volání genotypu se velmi liší od distribuce založené na metodě ML. Konkrétně přístupy vyvolávající genotyp vedou k velkému přebytku nízkofrekvenčních SNP (MAF mezi 2% a 4%). Tento vzor odráží to, co bylo vidět v našich simulačních studiích (obrázek 3). Taky, pro metody volání genotypu, distribuce frekvence alely se dramaticky mění se změnami hloubky sekvenování. Proto, jak bylo uvedeno výše, když hloubka není příliš vysoká, metody volání genotypizace pravděpodobně zahrnují mnoho falešných SNP, které jsou sekvenačními chybami. Tyto chyby se objevují jako přebytek nízkofrekvenčních SNP ve frekvenční distribuci. Distribuce založená na metodě ML je stabilnější napříč hloubkami, ale stále existuje přebytek SNP s nízkou frekvencí alely s hloubkou menší než 9×ve srovnání s podílem nízkofrekvenčních SNP ve větších hloubkách.

Konečně, jsme použili tento exome-resekvencování data pro simulaci case-control studie sdružení. Zkoumat distribuci asociace testovací statistiky za nulové hypotézy, jsme náhodně přiřazeno 100 osob, pro případ, skupinu a ostatní 100 s kontrolní skupinou. Pro všechny Modifikace na chromozomu 2 s MAF odhadů > 2% (na základě neznámé minor alela ML metoda), testovali jsme pro frekvenci alel rozdíly mezi případy a kontrolami ze strany computing G-statistika pomocí tzv. genotypů, a to jak s a bez filtrování, stejně jako LRT statistika. Obrázek 8 ukazuje grafy QQ porovnávající rozdělení statistik testů se standardním rozdělením χ2 (1). Jak je vidět v simulačních studiích, nulová distribuce g-statistiky vypočtená při volání genotypů bez filtrování se podstatně odchyluje od distribuce χ2 (1). Nulová distribuce statistiky LRT však úzce sleduje rozdělení χ2 (1). Inflační faktor je 1,01, což znamená, že statistika LRT funguje dobře při aplikaci na reálná data.