(4)

där. Eftersom kan vara mycket liten med stora datamängder (t.ex. med många individer) är det användbart att beräkna sannolikheten i loggskalan. Beställ de tre villkorliga log-sannolikheterna till (l (1), l(2), l(3)), där l(1) är den största. Sedan

G-test med kallade genotyper för associeringskartläggning

i associeringsstudier sägs SNP som visar signifikanta skillnader i allelfrekvens mellan fall och kontroller vara associerade med fenotypen av intresse. Associeringsmappning kan utföras med hjälp av data från nästa generations sekvenseringsstudier. Vi diskuterar först tillvägagångssätt som kräver att man ringer enskilda genotyper och utför sedan ett test för association med hjälp av de kallade genotyperna. I detta tillvägagångssätt kallas en genotyp först för varje individ. Genotyperna kan filtreras eller ofiltreras. Förutsatt oberoende mellan individer och HWE, en 2 kub 2 beredskapstabell kan byggas genom att räkna antalet stora och mindre alleler i både fall och kontroller. Detta leder till det välkända sannolikhetsprovet för oberoende, G-testet:

(5)

där O k,h är den frekvens som observeras i en cell, och E k,h är den frekvens som förväntas under nollhypotesen där allelfrekvensen är densamma mellan fall och kontroller. Det välkända Pearsons chi-square-test är asymptotiskt ekvivalent med G-testet. Om tabellen genereras från sanna genotyper, följer G-statistiken asymptotiskt en chi-kvadratisk fördelning med 1 frihetsgrad (kub2(1)). Men i våra studier konstruerar vi g – statistiken med hjälp av ”kallade” genotyper, så HWE kanske inte håller på grund av över-och underkallning av heterozygoter. Dessutom konstruerar teststatistiken genom att räkna ”kallade” genotyper istället för ”observerade” genotyper sannolikt extra variation. Därför kanske den statistiska teorin inte längre är giltig. Observera att när en genotyp inte kallas för en viss individ anses data saknas och ingår inte i tabellen 2 oz 2.

sannolikhetsprov redovisning av osäkerhet i de observerade genotyperna för associeringskartläggning

istället för att ringa genotyper möjliggör sannolikhetsramen osäkerhet i genotyperna och testerna på varje plats j om allelfrekvensen är densamma mellan fall och kontroller. Formellt beräknar vi sannolikheten för hypoteserna H O: p j ,1 = P j, 2 (=p j, 0) och H A : p j, 1 kg P j ,2, där p j ,1 och p j, 2 är MAF i fall och kontroller, respektive.

förutsatt att mindre (m) och större (m) alleler är kända, sannolikheten för den mindre allelfrekvensen kan beräknas enligt beskrivningen i ekvation 2, och sannolikhetsgraden teststatistik beräknas som:

(6)

där och är observerade data för fall respektive kontroller och och är mles för Mafs i fall respektive kontroller.

om den mindre allelen är okänd beräknas sannolikheten under nollhypotesen som i ekvation 3 och LRT-statistiken modifieras som:

(7)

där D j är observerade data för både fall och kontroller, och är allelfrekvensen under nollhypotesen. Andra beteckningar är desamma som i ekvation 6.

uppskatta MAF i simulerade data

vi jämför uppskattningarna av allelfrekvens på simulerade data med hjälp av sanna genotyper (True), kallade genotyper utan filtrering (Call NF), kallade genotyper med filtrering (f = 1; Call F) och maximal sannolikhetsmetod (ML). För sällsynta SNP är den mindre alleltypen ofta inte uppenbar. När man ringer genotyper antas den näst vanligaste nukleotiden vara den mindre allelen. ML-metoden innehåller direkt osäkerhet vid bestämning av mindre allel och om inte annat anges visas resultat med den okända mindre allelmetoden (ekvation 3). Observera att den okända mindre allel ML-metoden fungerar på samma sätt som den kända mindre allel ML-metoden men den tidigare bättre för mycket sällsynta SNP: er (ytterligare fil 1).

vi utvärderade först hur väl de olika metoderna kunde uppskatta MAF hos 200 individer över ett intervall av sekvenseringsdjup för 1000 SNP med en sann MAF på 5%. Figur 1 visar boxplots av fördelningarna av uppskattade MAF: er med de fyra olika metoderna. Som förväntat, för högre täckningsdata, till exempel ett individuellt djup på 12 KB, utför alla metoder såväl som när genotyperna är kända med säkerhet (Sant). Men när djupet minskar blir uppskattningarna av MAF som erhållits genom att först ringa genotyper förspända. Till exempel är medianmaf som uppskattas med Call F-metoden 5,3% vid 6 kg täckning och är 12,5% vid 2 kg. Anledningen till den uppåtgående förspänningen är att det blir svårare att ringa heterozygoter eftersom sanna heterozygoter ofta ser ut som sekvenseringsfel. Därför tenderar fler heterozygoter än mindre homozygoter att sakna genotyper. Den övergripande bias i MAF-uppskattningar från kallade genotyper är emellertid inte alltid i en riktning (data visas inte). Intressant verkar bias vara värre för Call F-metoden än Call NF-metoden. Detta mönster kan verka kontraintuitivt eftersom filtrering av genotyp samtal verkar minska sannolikheten för att ringa ett sekvenseringsfel en heterozygot. Call F-metoden resulterar emellertid också i en större mängd saknade data eftersom många homozygoter för den stora allelen inte kommer att ringas på grund av sekvenseringsfel. Således, i detta fall, att ringa genotyper utan filtrering verkar vara den bättre strategin än att filtrera genotyper när man försöker uppskatta MAF.

uppskattningar av allelfrekvens på platser med en sann MAF på 5% för olika djup av täckning. Vid varje djup simulerades 1000 platser med 200 individer, och vid varje plats beräknas en uppskattning av allelfrekvensen med: (1) sanna genotyper (sanna); (2) kallade genotyper utan filtrering (samtal NF); (3) kallade genotyper med filtrering (samtal F); och (4) metoden för maximal sannolikhet (ML). För mer information om uppskattningsmetoderna, se metoder.

resultaten är dramatiskt olika för den nya ML-metoden. Denna metod ger objektiva uppskattningar av MAF (median på ~4.9%) över ett djupområde. Även vid 2 kg visar uppskattningarna endast en något större variation än de som är baserade på de sanna genotyperna.

vi jämförde också det uppskattade genomsnittliga kvadratfelet (MSE; förväntan () för de olika uppskattningarna av MAF över ett intervall av sekvenseringsdjup (Figur 2). ML-metoden har en lägre MSE än anropsmetoderna med 50 eller 200 individer. I synnerhet är MSE beräknad baserat på Call F-metoden mycket högre än de från de andra metoderna, särskilt när djupet minskar. MSE för uppskattningarna av MAF baserat på de sanna genotyperna återspeglar den nedre gränsen för MSE och är inte konstant över djup på grund av provtagningsvarians och en ändlig provstorlek. Med hjälp av 50 individer närmar sig MSE 0,0005 med ökande djup och när man använder en provstorlek på 200 individer närmar sig den 0,0013 med ökande djup.

Mean squred error (MSE; förväntat) av fyra olika typer av allelfrekvensestimatorer för olika provstorlekar (vänster och höger panel) och djup av täckning (X-axel). Vid varje djup beräknades MSE från allelfrekvensestimat som gjordes med fyra olika metoder: True, Call NF, Call F och ML (för detaljer om metoderna, se bildtexten i Figur 1).

sammantaget utför den nya ML-metoden genotyp anropande metoder.

uppskatta en fördelning av MAF från simulerade data

Vi undersöker nästa hur de olika uppskattningsmetoderna utfördes vid uppskattning av andelen SNP vid olika frekvenser i befolkningen (liknande platsfrekvensspektret men baserat på populationsallelfrekvens istället för provfrekvens). Här simulerade vi 20 000 SNP där distributionen av de sanna MAF: erna följde standard stationär distribution för en effektiv befolkningsstorlek på 10 000 (se metoder). Observera att det i praktiken är mycket svårt att skilja en mycket sällsynt SNP från ett sekvenseringsfel. Därför, för jämförelse ändamål med verkliga data, kasserade vi SNP med beräknad MAF mindre än 2%. Figur 3 visar andelen SNP som faller i varje frekvensbehållare efter att ha exkluderat dessa SNP med beräknad MAF<2%.

fördelning av Allelfrekvenser av SNPs simulerade förutsatt att standard stationär distribution av allelfrekvenser. Vid varje djup (varje panel) simulerades 20 000 SNP, och för varje SNP erhölls uppskattningar av MAF med fyra olika metoder (se bildtexten i Figur 1). För varje metod (varje färg) används endast platser med uppskattade allelfrekvenser > 2% för att generera varje histogram (x-axel).

som förväntat, med ett högt täckningsdjup, såsom 10 kg per individ, ger alla metoder uppskattade MAF-distributioner som liknar den förväntade fördelningen baserat på de sanna genotyperna (Figur 3). Med ett grundare täckningsdjup, såsom mindre än 4 kb per individ, avviker fördelningarna av MAF som erhållits med genotyp anropande metoder avsevärt från den förväntade MAF-distributionen baserat på sanna genotyper (Figur 3). I synnerhet överskattar dessa metoder andelen lågfrekventa SNP. Till exempel är den förväntade andelen SNP i det andra facket (beräknad MAF mellan 2-4%) 18%. Den motsvarande andelen baserad på Call NF-metoden på ett djup av 4 kg är 26%, vilket är 1,4 gånger högre än väntat. Överskattningen av andelen lågfrekventa SNP uppstår på grund av förvirring av sekvenseringsfel med sanna heterozygoter, vilket resulterar i överkallning av heterozygot genotyper. Storleken på denna inflation skiljer sig åt mellan olika filtreringsavbrott, men en större cutoff ökar eller minskar inte nödvändigtvis inflationen.

bilden är helt annorlunda för ML-metoden. Den uppskattade MAF-distributionen som erhållits från den nya ML-metoden följer noggrant den verkliga fördelningen även med grunda täckningsdjup. Här finns det nästan inget överskott av lågfrekventa SNP. På ett djup av 4 kg är andelen SNP i histogrammets andra bin 18,4%, vilket ligger mycket nära den förväntade andelen (18%). Således kan mer tillförlitliga uppskattningar av frekvensspektret göras från data med låg täckning genom att använda vår sannolikhetsmetod än genom att använda genotyp anropande metoder.

Associationskartläggning i simulerade data

vi jämför prestanda för metoder som behandlar härledda genotyper som sanna genotyper i associeringstester (med hjälp av ett G-test) med vårt sannolikhetsprov (LRT) som står för osäkerhet i genotyperna. Vi undersöker fördelningen av teststatistiken under nollhypotesen om ingen allelfrekvensskillnad mellan fall och kontroller. Vi jämför också kraften i de olika metoderna.

med rimligt stora provstorlekar föreslår standard asymptotisk teori att under nollhypotesen följer både g-statistik och LRT-statistik en chi-kvadratisk fördelning med en grad av frihet (kub2(1)). Därför har vi jämfört nollfördelningen av G-statistiken beräknad baserat på anropsmetoder såväl som LRT-statistiken till distributionen av 2(1) med QQ-tomter (Figur 4). Vi simulerade 5000 SNP över en mängd sekvenseringsdjup i 500 fall och kontroller där MAF som användes för att simulera genotyper var 5% i båda fallen och kontrollerna. Fördelningen av G-statistiken som beräknas med hjälp av de sanna genotyperna visar en mycket god korrespondens med en distribution av 2(1) i form av en distribution. Fördelningen av G-statistiken som beräknas baserat på de kallade genotyperna avviker emellertid väsentligt från en distribution av 2(1) i form av en distribution. Att ringa genotyper och sedan behandla dessa genotyper som korrekta ger ett stort överskott av falskt positiva signaler om p-värdena beräknas med hjälp av en distribution av 2(1). Till exempel hade 11% av SNP: erna ett p-värde mindre än 5%, jämfört med de förväntade 5%. Effekten orsakas av en ökning av variansen på grund av överkallande homozygoter som heterozygoter, i det alleliska testet som används här för att upptäcka association. Genotypiska tester som Armitage trend test, som är robusta för avvikelser från Hardy-Weinberg jämvikt, visar inte en liknande ökning av den falska positiva frekvensen (ytterligare Fil 2). I överensstämmelse med denna observation resulterar filtrering av de kallade genotyperna i en minskning av fraktionen av signifikanta tester vid användning av G-testet, även om filtrering inte helt löser problemet. Å andra sidan visar LRT-statistiken endast en mycket liten avvikelse från en distribution av 2(1) för 2 eller 5 (5) djup av täckning.

QQ-diagram som jämför null-fördelningen av teststatistiken av intresse med en distribution av ubic2(1). Varje kolumn motsvarar en annan teststatistik: (1) g-statistik beräknad med de sanna genotyperna (True); (2) g-statistik beräknad med kallade genotyper utan filtrering (Call NF); (3) G-statistik beräknad med kallade genotyper med filtrering (samtal F); och (4) teststatistiken för sannolikhetsförhållande med okänd mindre allel (LRT). Förutsatt att 500 fall och 500 kontroller, under nollhypotesen, simulerades en uppsättning av 5000 platser med en MAF på 5% med ett sekvenseringsdjup på 2 kg (övre paneler) och 5 kg (nedre paneler). ”Inflation” – Faktorn visas i det övre vänstra hörnet av varje figur.

vi genererade också receiver operating characteristic (ROC) kurvor för var och en av de olika associeringstesterna. Dessa kurvor visar testets kraft vid olika falskt positiva hastigheter. Eftersom fördelningarna av en del av teststatistiken inte följer fördelningen av 2(1) enligt nollhypotesen, för att göra en rättvis jämförelse, fick vi det kritiska värdet för varje falsk positiv hastighet baserat på den empiriska nollfördelningen. Effekten beräknas som fraktionen av simulerade sjukdomslokaler som har en statistik som överstiger det kritiska värdet. Sammantaget finner vi att LRT presterar bättre än G-testet baserat på antingen genotyp anropande metod (Figur 5). Till exempel, vid en 5% falsk positiv hastighet och med ett sekvenseringsdjup på 5 kg, är kraften att upptäcka ett sjukdomslokus med en MAF på 1% och en relativ risk (rr) på 2 51% Med LRT, men kraften sjunker till 33% Med anropsmetoden utan filtrering och till 34% Med anropsmetoden med filtrering. I synnerhet vid lågt djup utför G-testet applicerat på kallade genotyper med filtrering mycket dåligt (vänster mest kolumn i Figur 5). Om vi jämför kraften i LRT med Armitage trend test med kallade genotyper, finner vi att LRT också har högre effekt än Armitage trend test (ytterligare fil 3). Detta tyder på att om man vill använda kallade genotyper kan filtrering av dem baserat på samtalsförtroende leda till strömförlust.

Receiver operating characteristic (ROC) kurvor av fyra tester av association. För definitionen av de fyra statistiken, se bildtexten i Figur 4. Förutsatt att 500 fall och 500 kontroller simulerades en uppsättning av 20 000 platser under null och under alternativet vid individuella sekvenseringsdjup på 2 kg, 5 kg och 10 kg (tre kolumner). Vid varje falsk positiv hastighet (x-axel) beräknades motsvarande kritiska värde med hjälp av den empiriska nollfördelningen. Den sanna positiva hastigheten (effekt; y-axel) erhölls genom att beräkna fraktionen av orsakssamband med teststatistik som överstiger det kritiska värdet.

Application to real data

vi analyserade 200 exomer från kontroller för en sjukdomsföreningsstudie som har sekvenserats med Illumina-teknik vid ett individuellt djup av 8 kg . Vi använde genotypens Sannolikhet genererad av” SOAPsnp ” – programmet för vår slutsats. För mer information, se metoder.

först undersökte vi noggrannheten i uppskattningarna av MAF från nästa generations sekvenseringsdata för 50 SNP genom att jämföra dem med de uppskattade MAF: erna från Sequenom-genotypdata. Både uppskattningarna med ML-metoden och genotyp-anropsmetoden utan filtrering är starkt korrelerade med uppskattningarna gjorda från Sequenom-genotypdata (dvs. en liten standardiserad skillnad mellan de två uppskattningarna i Figur 6). Uppskattningar baserade på genotypsamtal med filtrering visar emellertid dålig korrespondens med frekvenserna uppskattade från Sequenom-genotypdata, särskilt när sekvensdjupet är lågt. Intressant är att det finns en SNP där den uppskattade MAF från resekvenseringsdata skiljer sig mycket från uppskattningen som erhållits från Sequenom-genotypdata, även om sekvensdjupet är mycket högt (14 kg). Specifikt är den uppskattade MAF från Sequenom-genotypdata 22,5%, men är 17,2% när den uppskattas med ML-metoden. Individuell undersökning visar att hos många individer skiljer sig den starkt stödda genotypen baserad på sekvenseringsdata från Sequenomgenotyperna. Med tanke på att denna SNP omfattas av många läsningar hos dessa individer och att de observerade läsbaserna har högkvalitativa poäng (>Q20), är det troligt att skillnaden beror på Sequenom genotypfel. Observera att det finns ett par SNP: er där de uppskattade MAF: erna från genotyp-anropsmetoden utan filtrering verkar bättre motsvara MAF: erna uppskattade från Sequenom-genotypningen än uppskattningarna från ML-metoden gör. Till exempel, vid en SNP är den uppskattade MAF 25, 7% från Sequenom-genotypdata, 25.9% från genotyp-anropsmetoden utan filtrering och 27,2% från ML-metoden. Individuell inspektion avslöjar emellertid att det finns några individer för vilka den kallade genotypen från sekvenseringsdata skiljer sig från Sekvensgenotypen. I dessa fall avbröts felen i de kallade genotyperna, vilket gav upphov till bättre korrespondens med Sequenom-genotypdata. För dessa SNP är det därför svårt att säga vilken metod som fungerar bäst.

uppskattningar av allelfrekvens beräknad från 200 individer som använder nästa generations sekvenseringsdata jämfört med Sekvensgenotypdata. På varje plats användes endast individer som har både Sequenom-genotypdata och sekvenseringsdata för uppskattning av allelfrekvens. För sekvenseringsdata erhölls uppskattningar av MAF med tre olika metoder (Ring NF; Ring F; och ML). Den standardiserade skillnaden för varje uppskattning beräknades som , där och är de uppskattade MAF: erna från sekvenseringsdata respektive Sekvensgenotypdata, och n är antalet individer som används för uppskattningen. Varje webbplats är indelade i fyra fack baserat på genomsnittlig individ djupet av täckning (färg): mindre än 4×, högre än 4× men mindre än 8×, högre än 8× men mindre än 16×, och högre än 16×.

vi undersökte därefter fördelningen av MAF: er beräknade med flera tillvägagångssätt över ett intervall av sekvenseringsdjup från vår nästa generations exome-sekvenseringsdata (Figur 7). Vi kasserade SNP: er med beräknad MAF <2% eftersom det är svårt att skilja dessa mycket lågfrekventa SNP: er från sekvenseringsfel i denna dataset. Vi tog vidare bort platser där det fanns en signifikant skillnad (p-värde mindre än 10-5 med hjälp av ett rank-sum-test ) i kvalitetsresultatet för läsbaser mellan de mindre och stora allelerna. Dessa platser kommer sannolikt att vara artificiella SNP som kan uppstå på grund av felaktig kartläggning eller okända fördomar infördes under den experimentella proceduren. Sedan klassificerade vi varje webbplats i lagerplatser baserat på täckningsdjupet. Antalet SNP i varje fack visas i Tabell 1. När medeldjupet är mindre än 9 ml skiljer sig fördelningarna av uppskattade MAF: er baserade på genotyp-anropsmetoder mycket från den som baseras på ML-metoden. Specifikt ger genotyp anropande tillvägagångssätt upphov till ett stort överskott av lågfrekventa SNP (MAF mellan 2% och 4%). Detta mönster speglar vad som sågs i våra simuleringsstudier (Figur 3). Också, för genotyp anropande metoder, allel frekvensfördelningen förändras dramatiskt som sekvensering djupförändringar. Därför, som diskuterats tidigare, när djupet inte är särskilt högt, kommer genotypningssamtalsmetoderna sannolikt att innehålla många falska SNP: er som är sekvenseringsfel. Dessa fel visas som ett överskott av lågfrekventa SNP i frekvensfördelningen. Fördelningen baserad på ML-metoden är stabilare över djup, men det finns fortfarande ett överskott av SNP med låg allelfrekvens med djup mindre än 9 ml jämfört med andelen lågfrekventa SNP på större djup.

fördelning av den mindre allelfrekvensen uppskattad från exomerna hos 200 sekvenserade individer. För varje plats uppskattades den mindre allelfrekvensen med fyra olika metoder: (1) ML-metoden med okänd mindre allel, (2) ML-metoden med en känd eller fast mindre allel, (3) anropande genotyper utan filtrering (samtal NF) och (4) anropande genotyper med filtrering (samtal F). Varje webbplats klassificeras i lagerplatser baserat på täckningsdjupet. Vidare beaktas inte platser med beräknad MAF mindre än 2% i varje histogram. För antalet SNP som användes för denna analys, se Tabell 1.

slutligen använde vi denna exome-resequencing data för att simulera en fallkontrollföreningsstudie. För att undersöka fördelningen av föreningsteststatistiken under nollhypotesen tilldelade vi slumpmässigt 100 individer till en fallgrupp och de andra 100 till kontrollgruppen. För alla SNP på kromosom 2 med MAF-uppskattningar > 2% (baserat på den okända mindre allel ML-metoden) testade vi för allelfrekvensskillnader mellan fall och kontroller genom att beräkna G-statistiken med kallade genotyper både med och utan filtrering samt LRT-statistiken. Figur 8 visar QQ-tomterna som jämför utdelningarna av teststatistiken med standarddistributionen för 2(1). Som framgår av simuleringsstudier avviker null-fördelningen av G-statistiken som beräknas när genotyper anropas utan filtrering väsentligt från distributionen av kub2(1). Null-fördelningen av LRT-statistiken följer emellertid noggrant fördelningen av 2(1) i enlighet med den. Inflationsfaktorn är 1.01, vilket innebär att LRT-statistiken fungerar bra när den tillämpas på verkliga data.

QQ-plottar som jämför associeringsteststatistiken för allelfrekvensskillnader mellan 100 fall och 100 kontroller till en distribution av 2(1). Fenotyper tilldelades slumpmässigt till indivdiduals i exome resequencing dataset så att det finns 100 fall och 100 kontroller. För varje webbplats beräknades tre statistik: g-statistiken använder kallade genotyper utan filtrering( samtal NF), G-statisk med kallade genotyper med filtrering (samtal F) och LRT-statistiken. För att minimera inkludering av falska SNP: er kasseras webbplatser med ML MAF-uppskattningar mindre än 2%. För visningsändamål visas resultat från platser på kromosom 2. Observera att inflationsfaktorn visas i det övre vänstra hörnet av varje QQ-plot.