Informace

H je sada všech možných diplomotypů, které jsou v souladu s daty genotypu


Jsem nový v Biology.SX. Mám statistiky a téměř nulové znalosti v genetice. Snažím se porozumět některým věcem souvisejícím s genetikou v určitém článku biostatistiky.

  1. Data haplotypu byla simulována pomocí haplotypových obrazců a frekvencí (ukázáno níže) pro 5-SNP podél oblasti citlivosti na dabetes na chromozomu 22, hlášené ve studii FUSION.

Vím, že haplotypy mohou být reprezentovány jako binární sekvence. Zajímalo by mě, proč všechno možné $2^5$ zde nejsou přítomny. (?)

Dokument to také říká

  1. Nechat $ G = (g_1, ldots, g_M) $ označují nefázovaná data genotypu pro $ M $ loci. $ mathcal {H} _G $ označují množinu všech možných diplotypů, které jsou v souladu s údaji o genotypu $ G $.

Pokud subjekt nese nejvýše jednu kopii kauzálního haplotypu '01100', patří k dominantnímu genetickému modelu, pokud nese dvě kopie tohoto haplotypu, patří do recesivního modelu.

Sada $ mathcal {H} _G $ není mi jasné. Při výpočtu pravděpodobnosti potřebuji vědět $ mathcal {H} _G $. Nějaká pomoc nebo návrhy?


Disclaimer: Stejně jako OP vím o genetice velmi málo a předpokládám, že ostatní lidé na webu mohou dát lepší odpovědi než já. Každopádně, protože otázka byla měsíce nezodpovězena, posílám svou odpověď. Snad to někdo vylepší.

Jen pozor, že „možné“ v tomto kontextu neznamená všechny haplotypy, které si dokážeme představit. Jak píše list, "$ mathcal {H} _G $ označuje množinu všech možných diplotypů, které jsou v souladu s údaji o genotypu $ G $" a to je další omezení, které umožňuje jen některé diplotypy - uvedených sedm, nikoli 32, které si dokážeme představit.

Dále si všimněte, že součet frekvencí v tabulce se rovná 1, proto v sadě nemohou být žádné další haplotypy.

Pomohlo by však, kdybyste propojili zdroj, odkud jste tato data a uvozovky získali.


Imputace genotypu pomocí transformace Poziční Burrows Wheeler

Imputace genotypu je proces predikce nepozorovaných genotypů ve vzorku jedinců pomocí referenčního panelu haplotypů. Za posledních 10 let se velikost referenčních panelů zvýšila více než 100krát. Zvětšení velikosti referenčního panelu zvyšuje přesnost markerů s nízkými frekvencemi vedlejších alel, ale představuje stále větší výpočetní problémy pro metody imputace. Zde představujeme IMPUTE5, metodu imputace genotypu, kterou lze škálovat na referenční panely s miliony vzorků. Tato metoda pokračuje v upřesňování pozorování provedeného metodou IMPUTE2, že přesnost je optimalizována použitím vlastní podmnožiny haplotypů při imputaci každého jednotlivce. Rychlé, přesné a paměti efektivní imputace dosahuje výběrem haplotypů pomocí transformace poziční Burrows Wheeler Transform (PBWT). Použitím datové struktury PBWT u genotypových markerů identifikuje IMPUTE5 lokálně nejlépe odpovídající haplotypy a dlouhé identické podle stavových segmentů. Metoda poté použije vybrané haplotypy jako kondicionační stavy v rámci modelu IMPUTE. Pomocí referenčního panelu HRC, který má ∼65 000 haplotypů, ukazujeme, že IMPUTE5 je až 30krát rychlejší než MINIMAC4 a až 3x rychlejší než BEAGLE5.1 a že využívá méně paměti než obě tyto metody. Pomocí simulovaných referenčních panelů ukazujeme, že IMPUTE5 měří sublineárně s velikostí referenčního panelu. Například udržování konstantního počtu imputovaných markerů a zvýšení velikosti referenčního panelu z 10 000 na 1 milion haplotypů vyžaduje méně než dvojnásobek času výpočtu. Jak se velikost referenčního panelu zvětšuje, je IMPUTE5 schopen využívat menší počet referenčních haplotypů, čímž snižuje výpočetní náklady.


Pozadí

Identifikace genetické varianty založená na referencích zahrnuje dva související procesy: genotypování a fázování. Genotypizace je proces určování, které genetické varianty jsou přítomny v genomu jedince. Genotyp v daném místě popisuje, zda obě chromozomální kopie nesou variantní alelu, zda ji nese pouze jedna z nich, nebo zda variantní alela není vůbec přítomna. Fázování se týká určování haplotypů jednotlivce, které se skládají z variant, které leží blízko sebe na stejném chromozomu a jsou zděděny společně. K úplnému popisu všech genetických variací v organismu je zapotřebí jak genotypizace, tak fázování. Společně se tyto dva procesy nazývají diplotypování.

Mnoho stávajících variantních analytických potrubí je navrženo pro krátká čtení sekvenování DNA [1, 2]. Ačkoli jsou krátká čtení velmi přesná na úrovni základny, mohou trpět obtížným jednoznačným zarovnáním s genomem, zejména v opakujících se nebo duplikovaných oblastech [3]. Výsledkem je, že miliony bází referenčního lidského genomu nejsou v současné době spolehlivě genotypizovány krátkými čteními, primárně ve vícemegabázových mezerách v blízkosti centromer a krátkých ramen chromozomů [4]. Zatímco krátká čtení nemohou být jedinečně mapována do těchto oblastí, dlouhá čtení mohou potenciálně zasahovat do nich nebo dokonce přes ně. Dlouhá čtení se již ukázala jako užitečná pro haplotypování založené na čtení, detekci velkých strukturálních variant a sestavování de novo [5–8]. Zde demonstrujeme užitečnost dlouhých čtení pro komplexnější genotypizaci. Vzhledem k historicky vyšším relativním nákladům a vyšším chybám sekvenování těchto technologií byla tomuto problému dosud věnována malá pozornost. Technologie sekvenování DNA s dlouhým čtením však rychle klesají na ceně a obecně se zvyšují. Mezi takové technologie patří sekvenování jedné molekuly v reálném čase (SMRT) od Pacific Biosciences (PacBio) a sekvenování nanopórů od Oxford Nanopore Technologies (ONT), obojí zde hodnotíme.

Problém genotypizace souvisí s úkolem odvozovat haplotypy z dlouho čitelných sekvenčních dat, o nichž existuje bohatá literatura a mnoho nástrojů [8–14], včetně našeho vlastního softwaru WhatsHap [15, 16]. Nejběžnější formalizací rekonstrukce haplotypu je problém s minimální opravou chyb (MEC). Problém MEC se snaží rozdělit čtení podle haplotypu tak, aby bylo nutné opravit minimální počet chyb, aby byly čtení ze stejného haplotypu navzájem konzistentní. V zásadě by tato formulace problému mohla sloužit k odvození genotypů, ale v praxi se vychází z předpokladu „všech heterozygotů“: nástroje pro rekonstrukci haplotypu obecně předpokládají, že jako vstup je uveden soubor heterozygotních poloh, který pracuje výhradně na těchto místech.

Navzdory tomuto obecnému nedostatku nástrojů byly navrženy některé metody pro genotypování pomocí dlouhého čtení. Guo a kol. [17] popisují způsob volání a nukleotidové rekonstrukce jedno nukleotidové varianty (SNV) s dlouhým čtením a rekonstrukce haplotypu, která identifikuje příklad čtení na každém místě SNV, které nejlépe odpovídá blízkým čtení, překrývající se s místem. Poté rozdělí čtení kolem místa na základě podobnosti s příkladem na sousedních místech SNV. Tato metoda však nezaručuje nalezení optimálního rozdělení čtení mezi haplotypy a autoři uvádějí poměrně vysokou míru falešných objevů (15,7%) a míru falešně negativních (11,0%) pro data PacBio z NA12878, což odpovídá skóre F1 pouze 86,6%. Kromě toho dvě skupiny v současné době vyvíjejí volající na základě učení, které ukazují, že mohou být vyladěny tak, aby fungovaly pomocí dlouhých, hlučných čtení: V nedávném předtisku Luo a kol. [18] popisují metodu, která využívá konvoluční neuronovou síť (CNN) k volání variant z dat s dlouhým čtením, které hlásí k dosažení skóre F1 mezi 94,90 a 98,52%, v závislosti na parametrizaci (při nácviku čtení dat od jednoho jednotlivce a volání variant na jiného jedince, viz tabulka 3 z [18]). Poplin a kol. [19] představují další nástroj založený na CNN, který dosahuje skóre F1 92,67% na datech PacBio (podle doplňkové tabulky 3 z [19]). Tato opatření se zdají slibná, nicméně tyto metody systematicky nevyužívají informace o propojení mezi variantami poskytovanými dlouhým čtením. Nevyužívají tedy jednu z klíčových výhod dlouhého čtení.

Pro ilustraci potenciálního přínosu používání dlouhých čtení k diplotypu napříč sousedními místy zvažte obr. 1a. Jsou ukázány tři pozice SNV, které jsou pokryty dlouhým čtením. Šedé sekvence představují skutečné sekvence haplotypu a čtení jsou zbarvena modře a červeně, kde barvy odpovídají haplotypu, ze kterého pochází příslušné čtení: červené z horní sekvence a modré z dolní. Protože se mohou vyskytnout chyby sekvenování, alely podporované čtením nejsou vždy stejné jako ty pravé v haplotypech zobrazených šedě. Když vezmeme SNV jednotlivě, bylo by rozumné nazvat první jako A/C, druhé jako T/G a třetí jako G/C, protože počet čtení podporujících každou alelu je stejný. To vede ke špatné předpovědi pro druhé SNV. Pokud bychom však věděli, z jakého haplotypu každý čtený pramení, tedy pokud bychom znali jejich barvy, pak bychom věděli, že na druhém místě SNV musí být chyby sekvenování. Jelikož čtení pocházející ze stejných haplotypů musí podporovat stejné alely a mezi haplotypovými čteními v tomto místě existují nesrovnalosti, musí být jakákoli predikce genotypu v tomto místě považována za vysoce nejistou. Používání informací o haplotypech během genotypizace proto umožňuje detekovat nejistoty a potenciálně vypočítat spolehlivější předpovědi genotypu.

Motivace a přehled diplotypování. A Šedé sekvence ilustrují haplotypy, které jsou přečtené zobrazeny červeně a modře. Červené čtení pochází z horního haplotypu, modré z dolního. Genotypizace každého SNV jednotlivě by vedla k závěru, že všichni jsou heterozygotní. Použití kontextu haplotypu odhaluje nejistotu ohledně genotypu druhého SNV. b Ve směru hodinových ručiček, počínaje vlevo nahoře: nejprve jsou sekvenční čtení zarovnána s referenčním genomem zadána jako druhá, zarovnání čtení jsou použita k nominaci kandidátských variant (červené svislé pruhy), které se vyznačují rozdíly vůči referenčnímu genomu třetímu, skrytému Markovu model (HMM) je konstruován tak, že každá kandidátská varianta dává vznik jedné „řadě“ stavů, které představují možné způsoby přiřazení každého čtení jednomu ze dvou haplotypů a také možných genotypů (podrobnosti viz část „Metody“) dále, k provedení se používá HMM diplotypovánítj. odvozujeme genotypy každé kandidátské varianty a také to, jak jsou alely přiřazeny k haplotypům

Příspěvky

V tomto článku ukazujeme, že u současných technologií dlouhého čtení lze fázovou inferenci založenou na čtení současně kombinovat s genotypizačním procesem pro SNV za účelem produkce přesných diplotypů a detekce variant v oblastech, které nelze mapovat krátkými čteními. Ukazujeme, že klíčem k tomuto závěru je detekce vazebných vztahů mezi heterozygotními místy v rámci čtení. Za tímto účelem popíšeme nový algoritmus pro přesnou předpověď diplotypů z hlučných dlouhých čtení, které se škálují do hluboce sekvenovaných lidských genomů.

Tento algoritmus pak použijeme k diplotypu jednoho jedince z projektu 1000 Genomes, NA12878, pomocí dlouhých čtení z PacBio i ONT. NA12878 byla rozsáhle sekvenována a studována a konsorcium genomu v láhvi publikovalo sady oblastí s vysokou spolehlivostí a odpovídající sadu vysoce spolehlivých variantních volání uvnitř těchto genomových oblastí [20]. Ukazujeme, že naše metoda je přesná, že ji lze použít k potvrzení variant v oblastech nejistoty a že umožňuje objevovat varianty v oblastech, které jsou nemapovatelné, pomocí technologií sekvenování krátkého čtení DNA.


H je sada všech možných diplomotypů, které jsou v souladu s údaji genotypu - biologie

Kopie zdroje snphap-1.3.1, napsaná Davidem Claytonem, vložena sem, protože stávající web se zavírá.

Nejsem správce, ale uživatel, který si chce ponechat kopii, ze které bych mohl v budoucnu instalovat. Následuje Davidova dokumentace.

Program pro odhad frekvencí velkých haplotypů SNP

Ústav lékařské genetiky Cambridge Institute for Medical Research Wellcome Trust/MRC Building Addenbrooke’s Hospital, Cambridge, CB2 2XY

Naše skupina a další nyní tento program docela využívají. Přichází však bez záruk. Rád bych věděl o jakýchkoli potížích/chybách, se kterými se lidé setkávají, a pokusí se je vyřešit, ale může to chvíli trvat.

Tento program implementuje poměrně standardní metodu pro odhad frekvencí haplotypů pomocí dat od nepříbuzných osob. Pro výpočet ML odhadů haplotypových frekvencí daných měřeními genotypu, které nespecifikují fázi, používá EM algoritmus. Algoritmus také umožňuje, aby některá měření genotypu chyběla (například v důsledku selhání PCR). Umožňuje také vícenásobnou imputaci jednotlivých haplotypů.

Známý algoritmus nejprve rozšíří data pro každý subjekt do kompletní sady (plně fázovaných) párů haplotypu v souladu s pozorovanými daty („instance haplotypu“). Algoritmus EM pak probíhá následovně:

  • Krok E: Vzhledem k aktuálním odhadům pravděpodobností haplotypu a za předpokladu Hardy-Weinbergovy rovnováhy vypočítejte pravděpodobnost každého fázovaného a úplného přiřazení genotypu pro každý subjekt. Změňte měřítko tak, aby součet činil 1,0 v každém předmětu. To jsou pak POSTERIÁRNÍ PROBABILITY přiřazení genotypů.
  • M krok :: Vypočítejte další sadu odhadů pravděpodobností haplotypu sečtením pozdějších pravděpodobností přes instance každého odlišného haplotypu. Po škálování na součet na 1,0 poskytují tyto nové odhady PRIOR PROBABILITIES.

Tento algoritmus je široce používán a není nový. U velkého počtu lokusů se však počet možných instancí haplotypu rychle stává nemožně velkým, přestože případné řešení může poskytnout znatelnou podporu pouze dosti omezené sadě haplotypů. Tento program se této obtížnosti vyhýbá tím, že začíná přizpůsobením 2-lokusových haplotypů a rozšiřuje řešení vždy o jeden lokus. Když se přidá každý nový lokus, počet instancí haplotypu, které je třeba zvážit, se nejprve rozšíří zvážením všech možných větších haplotypů. Poté, po aplikaci EM algoritmu k odhadu „předchozích“ haplotypových pravděpodobností a pozdějších pravděpodobností haplotypových instancí, jsou haplotypové instance vyřazeny dvěma způsoby:

  • Oříznutí posteriorní: Jakékoli přiřazení genotypu, jehož zadní pravděpodobnost klesne pod danou prahovou hodnotu, je odstraněno a zadní pravděpodobnosti přiřazení genotypu subjektu jsou přepočítány.
  • Předchozí oříznutí :: Všechny instance jakéhokoli haplotypu, jehož předchozí pravděpodobnost klesne pod prahovou hodnotu, budou odstraněny. Tato možnost se ve výchozím nastavení nepoužívá, protože může vést k obtížím při porovnávání pravděpodobností (viz níže).

Přidáváme po jednom lokusu až do dokončení.

Proces rušení přiřazení haplotypů v raných fázích může vést k řešením, která nejsou optimální. Například haplotyp 1.1.1 může mít nulovou odhadovanou frekvenci v analýze maximální pravděpodobnosti haplotypu se třemi lokusy, zatímco 1.1.1.x může mít nenulovou odhadovanou frekvenci v ML řešení problému čtyř lokusů. Není jasné, jak často to bude problém. Částečným řešením je zkusit zahrnout lokusy v různých řádech a zjistit, zda se získaná duše liší. Další ochranou není vyřazení haplotypů po zahrnutí každého lokusu, ale pouze každého k lokusu, i když bude penalizace jak v počítačovém čase, tak ve využití paměti, způsobené volbou velkých hodnot k. Všimněte si toho, že výběr k & gt1 získáte pouze tehdy, když bude každý EM algoritmus spuštěn z náhodného počátečního bodu (viz volba -rs).

Vzorkování (bayesovské) pozdější distribuce jednotlivých haplotypových dat se pohodlně provádí pomocí Gibbsova vzorkovače. To napodobuje EM algoritmus, ale používá spíše stochastické kroky než deterministické. Byl nazýván algoritmem IP (Imputation/Posterior Sampling). V našem případě algoritmus funguje následovně:

Krok I (nahrazuje krok E): U každého subjektu vyberte přiřazení haplotypu z možných instancí s pravděpodobností danou aktuální zadní nebo „plnou podmíněnou“ distribucí.

P-krok (nahrazuje M-krok): Ukázka frekvencí populace haplotypu z jejich úplného podmíněného Bayesian posterior. Pokud je převorem Dirichletova distribuce s konstantním df na všech možných haplotypech*, plná podmíněná zadní distribuce je také Dirichletova. Abychom získali sadu parametrů df pro tento pozdější Dirichlet, jednoduše přidáme konstantní df k počtu chromozomů aktuálně přiřazených ke každému haplotypu.

(* tj. pokud jsou relativní frekvence neznámého populačního haplotypu označeny p_1, p2, ..., p_i, ..., p_n, pak se jejich předchozí hustota považuje za úměrnou

kde d je parametr předchozího stupně volnosti)

Pokud je v předchozím kroku Dirchlet df brán jako nulový, pokud není v jednom kroku haplotyp přiřazen žádnému jednotlivci, může dojít k potížím se vzorkováním celého prostoru pomocí tohoto algoritmu, pak je plný podmíněný zadní nevhodný a haplotypu bude dána nulová pravděpodobnost v dalším kroku. Poté již nelze nikdy znovu odebrat vzorky. Pokud existuje pravděpodobnost více maxim, může se algoritmus „zaseknout“ pod jedním vrcholem. Aby se předešlo těmto obtížím, je učiněno opatření ke spuštění předchozího parametru df na relativně velké hodnotě, což dává všem haplotypům znatelnou pravděpodobnost vzorkování. Poté se předchozí parametr df v každém kroku sníží. Tento algoritmus se opakuje pro pevný počet kroků k získání jediné imputace. Předchozí parametr df se poté nastaví zpět na vysokou hodnotu, frekvence haplotypu populace se obnoví na jejich MLE a proces se opakuje, aby se získala další imputace. A tak dále.

Varování: Přestože je vícenásobná imputace pomocí IP algoritmu zavedenou technikou (viz Schafer J.L. „Analýza neúplných vícerozměrných dat“ Chapman a Hall: London, 1997), zůstává v této aplikaci přísně validována.

Je dobře známo, že pravděpodobnostní povrch tohoto problému může mít více maxim a že EM algoritmus bude konvergovat pouze k lokálnímu maximu. Poté, co byly přidány všechny lokusy a byl vypočítán konečný oříznutý seznam instancí haplotypu, může být EM algoritmus několikrát opakován z náhodných počátečních bodů za účelem hledání globálního maxima. Náhodná výchozí místa mohou být vybrána jedním ze dvou způsobů: (a) z náhodně vybraných hodnot pro předchozí pravděpodobnosti haplotypu, nebo (b) z náhodně vybraných pozdějších pravděpodobností pro každé přiřazení haplotypu. V první sadě iterací EM lze také vybrat náhodná výchozí místa a v tomto případě se použije metoda (b).

Program je vyvolán z příkazového řádku pomocí

Vstupní soubor by měl obsahovat data v předmětném pořadí s identifikátorem subjektu následovaným dvojicemi alel každého lokusu. Identifikátor subjektu nemusí být číselný, ale nesmí obsahovat „prázdné znaky“ (mezery nebo tabulátory). Alely by měly být buď kódovány 1, 2 (číselné kódování), nebo A, C, G nebo T („nukleotidové“ kódování). Chybějící data jsou označena 0 v numerickém kódování a pro kódování nukleotidů jakýmkoli znakem, který nebyl dosud uveden. Datová pole by měla být oddělena libovolným „mezerou“ (libovolným počtem mezer, tabulátorů nebo znaků nového řádku).

Ve výchozím nastavení jsou lokusy přidány ve stejném pořadí, ve kterém jsou uvedeny ve vstupním souboru, ale volitelně je lze přidat do

  1. Obrácené pořadí
  2. Náhodná objednávka
  3. V sestupném pořadí úplnost údajů
  4. V sestupném pořadí frekvence menších alel (MAF)

Zatím mám jen málo zkušeností s účinkem změny pořadí zařazení. Myšlenkou (3) je zastavit příliš mnoho šíření možných haplotypů v raných fázích procesu, kdy je k dispozici málo dat, na jejichž základě lze spolehlivě vybrat vzácné haplotypy k likvidaci. Myšlenkou (4) je soustředit se nejprve na starší haplotypy.

Výstup pravděpodobnosti protokolu z tohoto programu by měl být používán s určitou opatrností, zvláště když bylo použito předchozí oříznutí, protože pravděpodobnosti, které nezohledňují stejné podmnožiny možných haplotypů, nemusí být srovnatelné.

Volitelné příznaky umožňují nastavit následující parametry:

Více možností imputace:

Pokud je příkaz vydán bez voleb nebo argumentů, zapíše se na obrazovku stručný popis dostupných možností.

  1. Zprávy o průběhu iterace (zapsané na obrazovku). Všimněte si, že některé emulátory terminálu, které poskytují „rolování“, mohou vážně zpomalit provoz programu. V takovém případě byste měli použít buď standardní xterm bez rolování, nebo vyvolat volbu -q, která tento výstup potlačí.
  2. Soubor se seznamem nalezených haplotypů a jejich pravděpodobností (výstupní soubor-1). Seznam je v sestupném pořadí pravděpodobnosti a je také uvedena kumulativní pravděpodobnost. Kumulativní pravděpodobnost je potlačena, pokud je v platnosti možnost -ss.
  3. Soubor se seznamem přiřazení haplotypů k subjektům (výstupní soubor-2). Tento soubor obsahuje všechna přiřazení, jejichž zadní pravděpodobnost přesahuje násobek nejpravděpodobnějšího přiřazení (viz -th možnost).
  4. Soubor (pojmenovaný „snphap-warnings“), který obsahuje jakékoli varovné zprávy.

Výstupní soubory output-file-1 a output-file-2 jsou v komprimovaném a snadno čitelném formátu. Alternativně je lze uložit jako textové soubory oddělené tabulátory vhodné pro čtení do tabulkového procesoru nebo statistického programu, jako je „Stata“. Oba názvy souborů jsou volitelné a chybějící argument lze označit jediným „.“ (tečka nebo tečka). Ale protože je třeba předpokládat, že chcete NĚKTERÝ výstup, vynechání obou názvů souborů způsobí, že program bude pro výstupní soubor-1 výchozí „snphap.out“.

V režimu vícenásobné imputace se vytvoří další řada souborů. Každá imputace způsobí, že bude zapsán nový soubor (nebo dvojice souborů). Názvy souborů jsou uvedeny na příkazovém řádku, ale připojí se řetězce .001, .002, .003 ... atd.

Dodán je primitivní Makefile. Toto používá kompilátor gcc a bude nutné jej upravit, pokud má být použit jiný kompilátor C. Také budete muset upravit možnosti CMP_FLAGS a LD_FLAGS (které poskytují příznaky používané kompilátorem ve fázích kompilace a načítání)

Pro uživatele systému Microsoft Windows doporučuji použít emulační balíček „Cygwin“ Unix. Vidět

Zjistil jsem, že nastavení LD_FLAGS na -lm mi fungovalo na Linuxu i Solarisu (toto je výchozí nastavení), ale na Cygwinu jsem musel tento příznak vynechat.

Výchozí generátor jednotných náhodných čísel (UNIFORM_RANDOM) je nastaven na standardní 48bitovou funkci „drand48“ a odpovídající secí funkce (RANDOM_SEED) je „srand48“. Pro systémy, které nepodporují 48bitové funkce (včetně Cygwinu), však lze zvolit 32bitové verze:

„drand ()“ je definováno jako makro hodnotící na (0,5+rand ())/(1+RAND_MAX).

Součástí je také soubor krátkých testovacích dat. Obsahuje typizaci 100 subjektů pro 51 SNP v malé oblasti. Testování programu:

Alternativně, pokud chcete do výstupu zahrnout názvy loků,

Děkujeme Newtonovi Mortonovi a Nikolasovi Maniatisovi za užitečné komentáře a návrhy k rané předchozí verzi. Děkuji také každému, kdo upozornil na chyby v předchozích verzích.


HaplotypeCaller Sledujte

HaplotypeCaller je schopen volat SNP a indels současně prostřednictvím místní de-novo sestavy haplotypů v aktivní oblasti. Jinými slovy, kdykoli program narazí na oblast vykazující známky variací, zahodí stávající informace o mapování a kompletně znovu sestaví čtení v této oblasti. To umožňuje HaplotypeCaller být přesnější při volání regionů, které jsou tradičně obtížně volat, například když obsahují různé typy variant blízko sebe. Díky tomu je HaplotypeCaller mnohem lepší při volání indelů než volající na základě polohy, jako je UnifiedGenotyper.

V pracovním postupu GVCF používaném pro škálovatelnou variantu vyvolávající data sekvence DNA spouští HaplotypeCaller na vzorek pro generování přechodného GVCF (není použito v konečné analýze), který pak může být použit v GenotypeGVCF pro společné genotypování více vzorků ve velmi efektivní způsob. Pracovní postup GVCF umožňuje rychlé přírůstkové zpracování vzorků při jejich sjíždění ze sekvenceru a také škálování na velmi velké velikosti kohort (např. 92K exomy ExAC).

Kromě toho je HaplotypeCaller schopen zpracovat nediploidní organismy i sdružená experimentální data. Všimněte si však, že algoritmy používané pro výpočet pravděpodobnosti variant nejsou příliš vhodné pro extrémní frekvence alel (vzhledem k ploidii), takže jeho použití se nedoporučuje pro objevování somatických (rakovinných) variant. Za tímto účelem použijte místo toho Mutect2.

Nakonec je HaplotypeCaller také schopen správně zpracovat spojovací spoje, díky nimž je RNAseq výzvou pro většinu variantních volajících, za podmínky, že vstupní načtená data byla dříve zpracována podle našich doporučení, jak je zdokumentováno zde.

Jak HaplotypeCaller funguje

1. Definujte aktivní oblasti

Program určuje, které oblasti genomu potřebuje operovat (aktivní oblasti), na základě přítomnosti důkazů variability.

2. Určete haplotypy sestavením aktivní oblasti

Pro každou aktivní oblast program sestaví graf podobný De Bruijnovi, aby znovu sestavil aktivní oblast a identifikoval, jaké jsou v datech možné haplotypy. Program poté znovu zarovná každý haplotyp s referenčním haplotypem pomocí algoritmu Smith-Waterman, aby identifikoval potenciálně variantní místa.

3. Určete pravděpodobnost haplotypů s ohledem na načtená data

Pro každou aktivní oblast program provede párové zarovnání každého čtení proti každému haplotypu pomocí algoritmu PairHMM. To vytváří matici pravděpodobnosti haplotypů vzhledem k načteným datům. Tyto pravděpodobnosti jsou poté marginalizovány, aby se získaly pravděpodobnosti alel pro každé potenciálně variantní místo s ohledem na přečtená data.

4. Přiřaďte vzorové genotypy

Na každé potenciálně variantní místo program aplikuje Bayesovo pravidlo, přičemž pomocí pravděpodobností alel daných čtenými daty vypočítá pravděpodobnosti každého genotypu na vzorek s ohledem na odečtená data pozorovaná pro daný vzorek. Poté je vzorku přiřazen nejpravděpodobnější genotyp.

Vstup

Zadejte soubory bam, ze kterých chcete volat varianty

Výstup

Buď soubor VCF nebo GVCF se surovými, nefiltrovanými hovory SNP a indel. Běžné VCF musí být před použitím v navazujících analýzách filtrovány buď rekalibrací variant (Best Practice), nebo tvrdým filtrováním. Pokud používáte pracovní postup GVCF, výstupem je soubor GVCF, který je nutné nejprve spustit přes GenotypeGVCF a poté filtrovat před další analýzou.

Příklady použití

Toto jsou příklady příkazů, které ukazují, jak spustit HaplotypeCaller pro typické případy použití. Podívejte se na dokumentaci metod pro základní pracovní postup GVCF.

Jedno ukázkové volání GVCF (výstupy mezilehlé GVCF)

Jedno ukázkové volání GVCF s anotacemi specifickými pro alely

Variantní volání s bamoutem k zobrazení přeuspořádaných čtení

Upozornění

  • Dosud jsme plně netestovali interakci mezi voláním založeným na GVCF nebo voláním více vzorků a funkcemi specifickými pro RNAseq. Používejte je v kombinaci na vlastní nebezpečí.

Zvláštní poznámka k ploidii

Tento nástroj je schopen zvládnout mnoho případů nediploidního použití, přičemž požadovanou ploidii lze zadat pomocí argumentu -ploidy. Všimněte si však, že velmi vysoké ploidie (jaké se vyskytují ve velkých souhrnných experimentech) mohou způsobit problémy s výkonem, včetně nadměrné pomalosti. Na vyřešení těchto omezení pracujeme.

Další poznámky

  • Při práci s daty bez PCR nezapomeňte nastavit `-pcr_indel_model NONE` (viz argument níže).
  • Při spuštění v režimech `-ERC GVCF` nebo` -ERC BP_RESOLUTION` je prahová hodnota spolehlivosti automaticky nastavena na 0. Toto nelze přepsat příkazovým řádkem. Prah lze nastavit ručně na požadovanou úroveň v dalším kroku pracovního postupu (GenotypeGVCFs)
  • Pro urychlení analýzy doporučujeme použít seznam intervalů. Podrobnosti najdete v tomto dokumentu.

PRŮMYSLOVÁ VYUŽITELNOST

Podle předkládaného vynálezu je možné provést odhad maximální pravděpodobnosti frekvencí haplotypu v populaci, distribuce diplotypu jednotlivců (pozdější distribuce pravděpodobnosti konfigurací diplotypu) a penetrace pomocí genotypových dat a fenotypových dat, bez nutnosti definitivního určení diplotypu konfigurace jednotlivců. Pokud je použit algoritmus podle předkládaného vynálezu, lze asociaci mezi existencí haplotypu a jednoho fenotypu testovat pomocí genotypových dat a fenotypových dat získaných například v důsledku kohortové studie, klinického hodnocení nebo případová kontrolní studie.


Tato práce byla podpořena Čínskou národní přírodovědnou nadací [čísla grantu 81460169 a 8196030236] a cenou „Excellence Award“, kterou financuje Grant pro rozvoj kreativního výzkumu z First Affiliated Hospital of Guangxi Medical University. Tato práce byla umožněna mým stipendiem financovaným IACN-ISN-HKSN.

LP: návrh studie, sběr a analýza dat, sběr vzorků, vedení experimentů a příprava rukopisu RXY: návrh studie, analýza dat a revize rukopisu YHL: návrh studie, sběr dat a vzorků, revize rukopisu MQM: sběr vzorků, vedení experimentu a sběr dat a QHZ: analýza dat. Všichni autoři přečetli a schválili konečný rukopis.


Metody

Sběr osiva a půdy

V létě 2015 jsme navštívili 192 divokých populací pokrývajících původní rozšíření H. annuus, H. petiolaris, a H. argophyllusa shromáždili semena od 21 do 37 jedinců z každé populace. Semena z deseti dalších populací H. annuus byly dříve odebrány v létě 2011. Tři až pět vzorků půdy (hloubka 0 - 25 cm) bylo odebráno pomocí jádrového jádra v každé populaci z celé oblasti, ve které byla sbírána semena. Půdy se sušily na vzduchu na poli, dále se sušily při 60 ° C v laboratoři a nechaly se projít sítem 2 mm, aby se odstranily kořeny a kameny. Půdy byly poté předloženy k analýze společnosti Midwest Laboratories Inc. (Omaha, NE, USA).

Společná zahrada

V létě 2016 bylo na polní stanici Totem Plant Science Field of University of British Columbia (Vancouver, Kanada) pěstováno deset rostlin z každé ze 151 vybraných populací. Páry rostlin ze stejné populace původu byly vysety pomocí zcela randomizovaného designu. . Nejméně tři květiny z každé rostliny byly pytlovány před anthézou, aby se zabránilo opylování, a ručně kříženy s jedincem ze stejné populace původu. Fenotypová měření byla prováděna v průběhu růstu rostlin a byly sbírány listy, stonky, květenství a semena a digitálně zobrazovány za účelem získání příslušných morfometrických dat (viz doplňková tabulka 1).

Příprava a sekvenování knihovny

Whole-genome shotgun (WGS) sequencing libraries were prepared for 719 H. annuus, 488 H. petiolaris, 299 H. argophyllus individuals, and twelve additional samples from annual and perennial sunflowers (Supplementary Table 1). Genomic DNA was sheared to ∼400 bp fragments using a Covaris M220 ultrasonicator (Covaris, Woburn, Massachusetts, USA) and libraries were prepared using a protocol largely based on Rowan a kol., 2015 51 , the TruSeq DNA Sample Preparation Guide from Illumina (Illumina, San Diego, CA, USA) and Rolhand a kol., 2012 52 . In order to reduce the proportion of repetitive sequences, libraries were treated with a Duplex-Specific Nuclease (DSN Evrogen, Moscow, Russia), following the protocols reported in Shagina a kol. 2010 10 and Matvienko a kol. 2013 53 , with modifications (see Supplementary Methods for details). All libraries were sequenced at the McGill University and Génome Québec Innovation Center on HiSeq2500, HiSeq4000 and HiSeqX instruments (Illumina, San Diego, CA, USA), to produce paired end, 150 bp reads. Libraries with fewer reads were re-sequenced to increase genome coverage. After quality filtering (see below), a total of 60.7 billion read pairs were retained, equivalent to 14.5 Tbp of sequence data.

Variantní volání

The call set included the 1518 samples described above, the Sunflower Association Mapping (SAM) population (a set of cultivated H. annuus lines 54 ), and wild Helianthus samples previously sequenced for other projects 54–56 , for a total of 2392 samples (Supplementary Table 1). The additional samples were included to improve SNP calling, and to identify haploblock genotypes. Sequences were trimmed for low quality using Trimmomatic 57 (v0.36) and aligned to the H. annuus XRQv1 genome 9 using NextGenMap 58 (v0.5.3). We followed the best practices recommendations of The Genome Analysis ToolKit (GATK) 59 , and executed steps documented in GATK’s germline short variant discovery pipeline (for GATK 4.0.1.2). During genotyping, to reduce computational time and improve variant quality, genomic regions containing transposable elements were excluded 9 . Since performing joint-genotyping on the whole ensemble of samples would have been computationally impractical, genotyping was performed independently on three per-species cohorts (H. annuus, H. argophyllus a H. petiolaris).

Variant quality filtering

Genotyping produced VCF files featuring an extremely large number of variant sites (222M, 78M and 167M SNPs and indels for H. annuus, H. argophyllus a H. petiolaris). Over the called portion of the genome, this corresponds to 0.07 to 0.2 variants per bp, with 30-47% percent of variable sites being indel variation. To remove low-quality calls and produce a dataset of a more manageable size, we used GATK’s VariantRecalibrator (v4.0.1.2), which filters variants in the call set according to a machine learning model inferred from a small set of “true” variants. In the absence of an externally-validated set of known sunflower variants to use as calibration, we computed a stringently-filtered set from top-N samples with highest sequencing coverage for each species (N=67 (SAM) samples for H. annuus, and N=20 otherwise). The stringency of the algorithm in classifying true/false variants was adjusted by comparing variant sets produced for different parameter values (tranche 100.0, 99.0, 90.0, 70.0, and 50.0). For each cohort, results for tranche = 90.0 were chosen for downstream analysis, based on heuristics: the number of novel SNPs identified, and improvements to the transition/transversion ratio (towards GATK’s default target of 2.15).

Remapping sites to the HA412-HOv2 reference genome

Our initial analysis of haploblocks (see section “Population genomic detection of haploblocks”), as well as GWA/GEA results for haploblocks regions, found many instances of disconnected haploblocks and high linkage between distant parts of the genome, suggesting problems in contig ordering. We remapped genomic locations from XRQv1 9 to HA412-HOv2 11 using BWA 60 . Measures of LD using vcftools 61 showed that remapping significantly improved LD decay (Extended Data Fig. 1a) and produced more contiguous haploblocks (Extended Data Fig. 1b), supporting the accuracy of the new genome assembly and our remapping procedure. While we recognize that this approach reduces accuracy at the local scale, and would not be appropriate, for example, for determining the effects of variants on coding sequences, it produces a more accurate reflection of the genome and linkage structure.

Fylogenetická analýza

Variants were called for 20 windows of 1 Mbp, randomly selected across the genome. Indels were removed and SNP sites were filtered for <20% missing data and minor allele frequency >0.1%. All sites were then concatenated and analyzed using IQtree 62–64 with ascertainment bias and otherwise default parameters.

Mapování asociací v celém genomu

Genome-wide association analyses were performed for 86, 30 and 69 phenotypic traits in H. annuus, H. argophyllus a H. petiolaris, respectively, using the EMMAX (v07Mar2010) or the EMMAX module in EasyGWAS 65 an annotated list of candidate genes is reported in Supplementary Table 2. Inflorescence and seed traits could not be collected for H. argophyllus, since most plants of this species flowered very late in our common garden, and failed to form fully-developed inflorescences and set seeds before temperatures became too low for their survival.

Genome-environment association analyses

Twenty-four topo-climatic factors were extracted from climate data collected over a 30-year period (1961-1990) for the geographic coordinates of the population collection sites, using the software package Climate NA 66 . Soil samples from each population were also analyzed for 15 soil properties (Supplementary Table 1). The effects of each environmental variable were analyzed using BayPass 67 version 2.1. Following Gautier, 2015 67 , we employed Jeffreys’ rule 68 , and quantified the strength of associations between SNPs and variables as “strong” (10 dB ≤ BFje < 15 dB), “very strong” (15 dB ≤ BFje < 20 dB) and decisive (BFje ≥ 20 dB). An annotated list of candidates genes from GEA analyses is reported in Supplementary Table 2.

Transgenes and expression assays

The complete coding sequences (CDS) of HaFT1, HaFT2 a HaFT6 were amplified from complementary DNA (cDNA) from H. argophyllus individuals carrying the early and late haplotype for arg06.01. Two alleles of the HaFT2 CDS were identified in late-flowering H. argophyllus plants (one of them identical to the HaFT2 CDS from early-flowering individuals), differing only for two synonymous substitutions at position 285 and 288. All alleles were placed under control of the constitutive CaMV 35S promoter in pFK210 derived from pGREEN 69 . Constructs were introduced into plants by Agrobacterium tumefaciens-mediated transformation 70 . Col-0 and ft-10 seeds were obtained from the Arabidopsis Biological Resource Center. All primer sequences are reported in Supplementary Table 3.

Population genomic detection of haploblocks

The program lostruct (local PCA/population structure) was used to detect genomic regions with abnormal population structure 28 . Lostruct divides the genome into non-overlapping windows and calculates a PCA for each window. It then compares the PCAs derived from each window and calculates a similarity score. The matrix of similarity scores is then visualized using a multidimensional scaling (MDS) transformation. Lostruct analyses were performed on the H. annuus, H. argophyllus, H. petiolaris petiolaris, a H. petiolaris fallax datasets, as well as in a H. petiolaris dataset including both H. petiolaris petiolaris a H. petiolaris fallax Jednotlivci. For each dataset, lostruct was run with 100 SNP-wide windows and independently for each chromosome. Each MDS axis was then visualized by plotting the MDS score against the position of each window in the chromosome.

Many localized regions of extreme MDS values with high variation in MDS scores and sharp boundaries were detected (Fig. 4a Extended Data Fig. 4). Localized changes to population structure could occur due to selection or introgression, but both the size and discrete nature of the regions are consistent with underlying structural changes defining the boundaries and preventing recombination. For example, inversions prevent recombination between orientations and if inversion haplotypes are diverged enough, they will show up in lostruct scans 28 . Since we are interested in recombination suppression in the context of adaptation, we focused on regions that had the following features: (1) a PCA in the region should divide samples into three groups representing 0/0, 0/1 and 1/1 genotypes, (2) the middle 0/1 genotype should have higher average heterozygosity and (3) there should be high linkage disequilibrium (LD) within the region.

The combined evidence of PCA and linkage suggests that the lostruct outlier regions are characterized by long haplotypes with little or no recombination between haplotypes. We refer to these as haploblocks. To explore the haplotype structure underlying the haploblocks, sites correlated (R 2 > 0.8) with PC1 in the PCA of the haploblock were extracted as haplotype diagnostic sites and used to genotype the haploblocks. Since there is seemingly little recombination between haplotypes, this is conceptually similar to a hybrid index and we expect all samples to be consistently homozygous for one haplotypes alleles or be heterozygous at all sites (i.e. similar to an F1 hybrid). Haploblock genotypes were assigned to all samples using equation (1), where p is the proportion of haplotype 1 alleles and h is the observed heterozygosity. The haplotype structure was also visualized by plotting diagnostic SNP genotypes for each sample, with samples ordered by the proportion of alleles from haplotype 1 (e.g. Fig. 2f).

Lostruct was run in SNP datasets containing H. petiolaris petiolaris, H. petiolaris fallax, and both subspecies together. Although each dataset produced a collection of haploblocks, they were not identical. Some haploblocks were identified in one subspecies, but not the other, and some were only identified when both subspecies were analyzed together. In some cases, it was clear that haploblocks identified in both subspecies represented the same underlying haploblock because they physically overlapped and had overlapping diagnostic markers. We manually curated the list of haploblocks and merged those found in multiple datasets. We set the boundaries of these merged haploblocks to be inclusive (i.e. include windows found in either) and the diagnostic markers to be exclusive (i.e. only include sites found in both). For this merged set of haploblocks, all H. petiolaris samples were genotyped using diagnostic markers.

Design of genetic markers for haploblock screening

Diagnostic SNPs for haploblocks were extracted from filtered vcf files. The resulting markers Cleaved-Amplified Polymorphic Sequence (CAPS) or direct sequencing markers were tested on representative subsets of individuals included in the original local PCA analysis (Fig. 4a, Extended Data Fig. 4), for which the genotype at haploblocks of interest was known. Marker information are reported in Supplementary Table 3.

Sequencing coverage analysis

To detect the presence of potential deletions in the late-flowering allele of arg06.01, SNP in the haploblock region with average coverage of at least 4 across at least one of the genotypic classes were selected (in order to exclude positions with overall low mapping quality). SNP positions with coverage 0 or 1 in one genotypic class were counted as missing data for that genotypic class (Extended Data Fig. 2c).

H. annuus reference assemblies comparisons

Masked reference sequences for the H. annuus cultivars HA412-HOv2 and PSC8 11, 12 were aligned using MUMmer 71 (v4.0.0b2). The programs nucmer (parameters -b 1000 -c 200 -g 500) and dnadiff within the MUMmer package were used. Only orthologous chromosomes were aligned together because of the high similarity and known conservation of chromosome structure. The one-to-one output file was then visualized in R and only included alignments where both sequences were > 5000 bp. Inversion boundaries and sequence identity between haplotypes were further determined using Syri 72 .

Genetic maps comparisons

Fourteen genetic maps were used: the seven H. annuus genetic maps used in the creation of the XRQv1 genome 9 three newly generated H. annuus maps obtained from wild X cultivar F2 populations (E.B.M.D., M.T., G.L.O., L.H.R., in preparation) two previously published H. petiolaris genetic maps obtained from F1 crosses 50 and two newly generated H. petiolaris maps (K.H., Rose L. Andrews, G.L.O., K.L.O., L.H.R., in preparation). Whenever necessary, marker positions relative to XRQv1 were re-mapped to the HA12-HOv2 assembly (see above). Six of the previously described H. annuus maps were obtained from crosses between cultivars (the seventh one was obtained from a wild X cultivar cross) in order to determine which haploblock could be expected to segregate in the genetic maps, all of the H. annuus SAM population lines were genotyped for each H. annuus haploblock using diagnostic markers identified in wild H. annuus. Ann01.01 and ann05.01 were found to be highly polymorphic in the SAM population, while other haploblocks were fixed or nearly fixed for a single allele. For all fourteen maps, marker order was compared to physical positions in the HA412-HOv2 reference assembly, and evidence for suppressed recombination or structural variation was recorded (Extended Data Table 1).

Pairs of H. petiolaris a H. argophyllus populations that diverged for a large number of haploblocks were selected. Individuals from these populations were genotypes using haploblock diagnostic markers (see “Design of genetic markers for haploblock screening”) to identify, for each species, a pair of individuals with different genotypes at the largest possible number of haploblocks. Chromosome conformation capture sequencing 36, 73 (Hi-C) libraries were prepared by Dovetail Genomics (Scotts Valley, CA, USA) and sequenced on a single lane of HiSeq X with 150 bp paired end reads. Reads were trimmed for enzyme cut site and base quality using the tool trim in the package HOMER 74 (v4.10) and aligned to the HA412-HOv2 reference genome using NextGenMap 58 (v0.5.4). Interactions were quantified using the calls ‘makeTagDirectory - tbp 1-mapq 10’ and ‘analyzeHiC -res 1000000 -coverageNorm’ from HOMER. Hi-C data were used in two ways to identify structural changes. First, the difference between interaction matrices for samples of the same species was plotted for each haploblock region where the two samples had different genotypes. Second, the difference between interaction matrices for H. annuus (using the HiC data that were generated to scaffold the HA412-HOv2 reference assembly 11 ) and each H. petiolaris a H. argophyllus sample were plotted.

Haploblock phenotype and environment associations

Since haploblocks are large enough to affect genome wide population structure, their associations with phenotypes of environmental variables may be masked when controlling for population structure. Therefore, a version of the variant file was created with all haploblock sites removed GWA and GEA analyses were performed as before, but kinship, PCA and genetic covariance matrix were calculated using this haploblock-free variant file. Regions of high associations co-localizing with haploblock regions were identified, and haploblocks were also directly tested by coding each haploblock as a single bi-allelic locus.

To examine the relative importance of haploblocks to trait evolution and environmental adaptation, association results were compared between haploblocks and SNPs. Using SNPs as a baseline allows to control for the correlation between traits or environmental variables. To make values comparable, both SNPs and SVs with minor allele frequency ≤ 0.03 were removed. Each locus was classified as associated (p < 0.001 or BFje > 10 dB) or not to each trait. The number of traits or climate variable each locus was associated with was then counted. The proportion of loci with ≥ 1 traits/climate variables associated for SNPs and haploblocks was then compared using prop.test in R 75 (Extended Data Fig. 9b).

Haploblocks phylogenies and dating

The phylogeny of each haploblock region was estimated by Bayesian inference using BEAST 76 1.10.4 for 100 genes within the region. The dataset was partitioned, assuming unlinked substitution and clock models for the genes, and analyzed under the HKY model with 4 Gamma categories for site heterogeneity: a strict clock, a “Constant Size” tree prior with a Gamma distribution with shape parameter 10.0 and a scale parameter 0.004 for the population size. Default priors were used for the other parameters. A custom Perl script was used to combine FASTA sequences and the model parameters into XML format for BEAST input. The Markov chain Monte Carlo (MCMC) process was run for 1 million iterations and sampled every 1000 states. The convergence of chains was inspected in Tracer 77 1.7.1. In order to estimate divergence times, the resulting trees were calibrated using a mutation rate estimate of 6.9 × 10 −9 substitutions/site/year for sunflowers 78 , and visualized with R package ggtree 79 and Figtree v1.4.4 80 . Divergence times were extracted from the trees and plotted showing the 95% highest posterior density (HPD) interval based on the BEAST posterior distribution. This was repeated for 100 non-haploblock genes to estimate the species divergence times.

For the 10 Mb region on chromosome 6 controlling flowering time in H. argophyllus, the early flowering haplotype grouped with H. annuus. To determine if it is the product of an ancient haplotype that has retained polymorphism only in H. annuus or if it is introgressed from H. annuus, the phylogeny of 10 representative H. argophyllus samples homozygous for each haploblock allele, as well as 200 H. annuus samples, was inferred using IQtree. SNPs from the 10 Mb region were concatenated and the maximum likelihood tree was constructed using the GTR model with ascertainment bias correction. Branch support was estimated using ultrafast bootstrap implemented in IQtree 62–64 with 1,000 bootstrap replicates. Phylogenies of haploblock arg03.01, arg03.02 and arg06.02 were inferred using the same approach. To explore intra-specific history of the H. petiolaris haploblocks, all samples homozygous for either allele for each haploblock were selected, and phylogenies were constructed using IQtree with the same settings.


Metoda

Elektrofyziologie

Cell lines stably expressing Kv11.1-1A or Kv11.1-3.1 channels were maintained as previously described (13). Patch clamp electrophysiology recordings were undertaken as previously described (13 a summary is provided in the data supplement that accompanies the online edition of this article).

Drug block was calculated as lék / Iřízení and dose response curves were fitted with a modified Hill equation: kde lék is current recorded in the presence of drug, řízení is current recorded in control conditions, D is the drug concentration, h is the Hill coefficient, and IC50 is the half maximal inhibitory concentration of D.

Data are presented as mean and standard error of the mean. Data were analyzed using one-tailed paired t tests and analysis of variance, followed by the Tukey t test for pairwise comparison. The significance threshold was set at 0.05.

Clinical Cohort

The clinical cohort consisted of patients randomly assigned to one of five antipsychotic medications during phase 1/1A (first drug assigned) of the CATIE trial. The details of the overall design for the CATIE study, genotyping of the KCNH2 SNPs, and the participants’ demographic characteristics have been described previously (6, 9).

Because of the outpatient and parallel design of the original CATIE study, information about compliance based on drug clearance is an important factor determining symptom change during the CATIE trial. Thus, we only analyzed treatment response from subjects of European ancestry for whom we had drug clearance and genotype data (N=362). We previously showed that drug clearance data substantially improve prediction of treatment response (14). As an ancillary study to the CATIE trial, blood samples were drawn during study visits to measure antipsychotic drug concentrations. Data were collected on the amount of the last dose of medication, time the last dose was taken, and time the blood sample was drawn. This information was used with the drug concentration data to estimate drug clearance for each subject based on nonlinear mixed-effect modeling using NONMEM, version 5 (GloboMax, Ellicott City, Md.) (12). A one-compartment linear model with first-order absorption (NONMEM ADVAN5) using the first-order estimation method was used to estimate drug clearance (12).

We used estimated drug clearance instead of plasma concentrations because it is a dose-independent and time-independent measure, which allows for comparison of drug exposure across all subjects, as described in detail elsewhere (14, 15).

For this analysis, we focused on three SNPs in KCNH2—rs3800779 (SNP1), rs748693 (SNP2), and rs1036145 (SNP3)—which have been associated with increased expression of the novel Kv11.1-3.1 isoform in human postmortem brain samples (8) and overall response to treatment in the CATIE trial (9). Since the three SNPs were in moderate to strong linkage disequilibrium (see Table S1 in the online data supplement), in order to reduce multiple testing and to gain statistical power for detecting association, we constructed three SNP diplotypes to be used for testing diplotype-by-risperidone interaction on the treatment response. Haplotype construction was performed and phased diplotype was assigned using the Phase program (16). Details of genotyping and construction of diplotypes are provided in the data supplement. Diplotype was grouped into three categories according to the number of minor alleles that a diplotype contains at SNP1 and SNP3, coded “0” for no minor allele of either SNP1 or SNP3, “1” for one or two copies of minor alleles, and “2” for three or four copies of minor alleles. The distribution of diplotypes in individuals with drug clearance data was consistent with the total European ancestry sample in phase 1/1A of the CATIE trial, suggesting minimal selection bias (see Table S2 in the data supplement).

Clinical Data Analysis

In the CATIE sample, because all patients were receiving treatment and because the time and number of Positive and Negative Syndrome Scale (PANSS) evaluations in the study varied considerably among subjects, we treated the baseline PANSS rating as “before treatment” and the last rating as “after treatment” to test for genetic variant-by-risperidone treatment interaction on the treatment response. Since each subject had two measures in the analysis, we used a general linear mixed model to incorporate the relatedness between two observations within a subject (9). We did not perform a separate analysis with only those subjects who completed the trial because that subset was too small (N=39).

We performed this analysis on all subjects for whom drug clearance data were available and for whom diplotypes were assigned with good confidence (N=362), and we controlled for potential covariates of sex, age, time on medication, and whether the patient completed the 18-month trial or discontinued medication before the end of 18 months and therefore switched to phase 2 of the trial.

Individuals who were on risperidone (N=88) had a mean estimated drug clearance rate of 20.62 L/hour (SD=10.73, range=3.61–40.05). Based on tertile distribution of the clearance data range, we classified individuals into three groups: slow (N=30), intermediate (N=29), and fast metabolism (N=29) groups.

The clinical response analysis consisted of a diplotype-by-treatment interaction, as previously described (9). To test our specific hypothesis of a differential effect of diplotype on the antipsychotic response to risperidone, which was based on the differential affinity of risperidone for Kv11.1 in contrast to all other drugs in the trial, we combined all other medications as one group (see the online data supplement for more detail). Because the mean and variance of estimated drug clearance varied with different drugs, we assigned an ordinal measure of 1 to 3 according to the tertile distribution of each drug clearance to make the estimated measurements comparable between drugs. Using an ordinal measure based on the tertiles of estimated drug clearance for each drug while adjusting for the type of drugs in the same model of analysis allowed us to capture the likely nonlinear relationship between estimated drug clearance and treatment response (9). For this analysis, however, since non-risperidone medications were all combined into one group, we considered the possibility that the effect of drug clearance using tertiles may be different between drugs and consequently may affect our assessment of the overall effect of drug clearance. Therefore, we performed a leave-one-out sensitivity analysis to assess for such a potential bias.


Our primary aim is to learn haplotype-cluster models from large training sets and use them to phase samples efficiently and accurately. Here we introduce some modifications to BEAGLE so that the algorithm is better suited to this aim. Our new algorithm is called Underdog.

BEAGLE only represents haplotypes that actually appear in the training examples. However, since we would like to phase new genotype samples that do not necessarily appear in the training set, we set the transition probability for allele a at a given SNP to

(Eq. B1)

kde nA is the number of times allele A is observed in training data, and nā is the number of times the other allele is observed. This is compared with the BEAGLE formula shown in Algorithm 3. Here, γ is a positive number between 0 and 1. To illustrate the rationale for this choice of transition probability, consider the bottom state of level 2 in Figure 2.1. Instead of having only one transition (to the bottom state in level 3) with 100% probability, we add a second transition for the blue allele (also to the bottom state in level 3) that is visited with probability γ. We define all transition probabilities in the haplotype-cluster model in this way. These transition probabilities are only noticeably different from the transition probabilities in BEAGLE when one allele occurs very infrequently in the training set within a given cluster of haplotypes. With this modification, Underdog allows for genotype phase based on haplotypes that did not appear in the training set.

Although the BEAGLE haplotype-cluster models are intended to be parsimonious, building these models from hundreds of thousands of haplotypes can still yield very large models with millions of states, making it difficult to phase genotype samples in a reasonable amount of time. To address this problem, we first observe that although there is typically a large number of possible ways of phasing a sample, most of these possibilities are extremely unlikely conditioned on a specific haplotype-cluster model. In other words, most of the probability mass is typically concentrated on a small subset of paths through the HMM. To avoid considering all possible paths (which is computationally expensive), at a given level d we retain the smallest number of states such that the probability of being in one of those states is greater than 1 - ε. Even for small values of ε, this heuristic dramatically decreases the computational cost of sampling from the HMM, and computing the most likely phase using the Viterbi algorithm (Figure B1), while incurring very few additional phasing errors.

Figure B1: Relationship between choice of HMM parameter ε and average computation time for phasing a genotype sample (based on chromosome 1 only). If we set ε = 0, the average sample phasing time is 63 seconds, and the average phasing error rate is 0.93%. For choices of ε that are larger, but not too large, we achieve comparable phasing accuracy with a dramatic reduction in computational expense. Note that the computation time here does not include file input/output, nor the time taken to merge the phasing results from multiple windows.

The second modification we make to BEAGLE concerns the criterion for deciding whether two haplotype clusters (tj., nodes of the haploid Markov model) should be merged during model learning (see Algorithm 4). Since the standard method is overly confident for frequencies that are close to 0 or 1, we regularize the estimates using a symmetric beta distribution as a prior. Specifically, haplotype clusters X a y are not merged unless the following condition is satisfied for some haplotype h:

(Eq. B2)

kde nX a ny are the sizes of clusters X a y. The posterior allele frequency estimates in this formula are

(Eq. B3)

kde nX(h) a ny(h) are the numbers of haplotypes that begin with haplotype h. We set the parameters of the Beta prior (the prior counts), α and β, to 0.5. Compare this criterion to the one used in Browning (2006), (also refer to Algorithm 3), which merges two clusters unless the following relation holds for some h:

(Eq. B4)

kde pX (h) is the proportion of haplotypes in cluster X with that begin with haplotype h, a py (h) is the proportion of haplotypes in cluster y that begin with h. We evaluated the phasing accuracy of the algorithm using a few different values for constant C and settled on C = 20.

Algorithm 4 is the modified version of BEAGLE's procedure (Algorithm 3) that applies Eq. B2 to merging haplotypes during model building.

For computational efficiency, on each chromosome we estimate the genotype phase within 500-SNP windows separately. This can result in a loss of phasing accuracy at the beginning and end of each window because information outside the window is ignored, and therefore there is less information about the genotypes at the two extremities of the window. To address this problem, we learn haplotype-cluster models in overlapping windows specifically, we use 500-SNP windows in which two adjacent windows on the same chromosome overlap by 100 SNPs. Since the final phasing estimates produced in the two windows may disagree in the overlapping portion, it is not immediately clear how to combine the phasing estimates from adjacent windows. We propose a simple solution to this problem. First, we select the SNP nearest the midpoint of the overlapping portion at which the genotype is heterozygous (that is, the two allele copies are not the same). We call this the "switch-point SNP." We then join the sequences from the overlapping windows that share the same allele at this switch-point SNP. For example, in Figure B2 we join the top sequence in the left-hand window with the bottom sequence in the right-hand window because they are both estimated to carry the blue allele at the selected switch-point SNP.

Figure B2: Underdog learns haplotype-cluster models in overlapping windows. This figure illustrates how we obtain the final genotype phase from these overlapping windows.


Podívejte se na video: Genetika IV pro 3 A (Leden 2022).