L’allele minore è l’allele meno frequente nella popolazione in un sito variabile. Per prima cosa descriviamo due approcci principali per stimare la frequenza dell’allele minore (MAF) in un particolare sito nel genoma. Il primo approccio consiste nell’inferire i singoli genotipi e nel trattare quei genotipi inferiti come completamente accurati quando si stima la MAF. Esaminiamo quindi le prestazioni di un quadro di verosimiglianza che prende direttamente in considerazione l’incertezza nell’assegnazione dei genotipi. Durante il nostro lavoro, assumiamo che tutti i siti di segregazione siano biallelici.

Stima della MAF da genotipi chiamati

Un modo per stimare la MAF dai dati di sequenziamento di prossima generazione è quello di chiamare prima un genotipo per ogni individuo usando i dati di sequenziamento, e quindi utilizzare quei genotipi come se fossero quelli veri. Questo è stato l’approccio tradizionalmente utilizzato per i dati di genotipo e Sanger sequenziamento dei dati. Non è chiaro quanto bene si esibirà quando applicato ai dati di sequenziamento di prossima generazione.

Un approccio di massima verosimiglianza può essere utilizzato per dedurre il genotipo per ogni individuo dai dati di sequenziamento di prossima generazione. In ogni sito j, per ogni individuo i, la probabilità per ciascuno dei tre genotipi possibili (supponendo che conosciamo l’allele minore) è data come:

(1)

D,i, j è la cosa osservata dati di sequenziamento in singolo sono al sito j, g , i, j ∈ {0, 1, 2} è il numero di alleli minori contenute nel genotipo di ogni individuo, e e controllo per errori di sequenziamento e di lettura di base di qualità, rispettivamente. I dati di sequenziamento osservati per ogni individuo possono essere pensati come l’allineamento delle letture al sito j tenendo conto dei punteggi di qualità di lettura. Questo è rappresentato come la probabilità di genotipo e si trova nel file di probabilità di genotipo (GLF) che viene prodotto in molti programmi che analizzano i dati di sequenziamento di nuova generazione, come SOAPsnp e MAQ .

Per assegnare un genotipo a un particolare individuo, la probabilità di ciascuno dei tre possibili genotipi può essere calcolata per l’individuo. Il genotipo con la più alta probabilità può quindi essere assegnato. Tuttavia, i ricercatori spesso preferiscono un criterio di chiamata più rigoroso e non assegneranno un genotipo a un individuo a meno che il genotipo più probabile non sia sostanzialmente più probabile del secondo più probabile. Qui i tre genotipi possibili sono ordinati in base alle loro probabilità: , dove g(k)corrisponde al genotipo con la k esima probabilità maggiore. Con una data soglia f, si può chiamare il genotipo g(1) se . Altrimenti, un genotipo non viene chiamato e il genotipo dell’individuo è considerato mancante. Un valore soglia comune di f è 1, indicando che il genotipo più probabile è almeno 10 volte più probabile del secondo più probabile. Si noti che questo tipo di filtraggio può comportare una maggiore confidenza per il genotipo “chiamato”, ma si traduce anche in più dati mancanti.

Stimatore di massima verosimiglianza della frequenza allele

Invece di stimare la MAF dai genotipi chiamati, un metodo di massima verosimiglianza (ML) introdotto da Kim et al. (vedi anche Lynch per un approccio simile) stima direttamente le MAF e tiene conto dell’incertezza del genotipo. In particolare, dato un allele minore, la probabilità di osservare i dati della sequenza ad ogni individuo i è ottenuta sommando le probabilità corrispondenti a tutti e tre i genotipi possibili.

Supponiamo che le tre probabilità di genotipo definite nell’equazione 1 siano disponibili. Usando la stessa notazione di cui sopra, sia DJ e pj i dati di sequenziamento osservati nel sito j e la corrispondente MAF, rispettivamente. La probabilità genotipo dato che minore frequenza allele può essere calcolata assumendo equilibrio Hardy-Weinberg (HWE). Quindi, supponendo che l’indipendenza tra gli individui, la probabilità che la MAF in questo locus è un prodotto di tutte le verosimiglianze calcolata su tutti gli N individui:

(2)

La stima ML di p j può essere calcolata direttamente massimizzare la probabilità per un limitato spazio di parametro utilizzando il Broyden-Fletcher-Goldfarb-Shanno (BFGS) il metodo o utilizzando il expectation-maximization (EM) algoritmo . Quando si utilizza l’algoritmo EM, l’aspettativa posteriore di un genotipo viene calcolata per ogni individuo e la media di quei posteriori viene ripetutamente aggiornata. La nostra implementazione di BFGS è stata più veloce dell’algoritmo EM. Ad esempio, per ottenere stime da 100.000 siti, BFGS ha richiesto ~16 secondi ma EM ha richiesto ~100 secondi. Tuttavia, la differenza di velocità può essere specifica per l’implementazione. Nel nostro caso, per entrambi i metodi, abbiamo interrotto l’aggiornamento dei parametri quando l’aumento della probabilità era inferiore a 0,001.

Stimatore di massima verosimiglianza con allele minore incerto

In pratica, spesso il secondo nucleotide più comune tra gli individui può essere usato come allele minore. Tuttavia, per SNP rari (ad esempio, MAF < 1%), è difficile determinare quale allele sia l’allele minore, poiché tutti e quattro i nucleotidi possono apparire in alcune letture a causa di errori di sequenziamento. Per affrontare questa situazione, descriviamo ora un quadro di verosimiglianza che tiene conto dell’incertezza nella determinazione dell’allele minore.

Supponiamo che per il sito j si conosca l’allele maggiore M. Si noti che decidere quale dei due alleli comuni è probabile che sia il maggiore non è importante poiché ci occupiamo principalmente di stimare le frequenze a SNP rari. Inoltre, per gli alleli con frequenze intermedie (circa il 50%), la distinzione tra allele maggiore e minore è meno importante. Assegnare gli altri tre nucleotidi non principali m1, m2 e m3. La probabilità introdotta nell’equazione 2 assume un allele maggiore fisso M e un allele minore fisso m. Pertanto, per consentire l’incertezza nella designazione dell’allele minore, la funzione di verosimiglianza può essere modificata come:

(3)

Inoltre, supponendo che uno dei tre possibili alleli minori è altrettanto probabile, otteniamo:

(4)

dove . Poiché può essere molto piccolo con grandi set di dati (ad esempio, con molti individui), è utile calcolare la probabilità nella scala di log. Ordina le tre probabilità di log condizionali come(l(1), l(2), l(3)), dove l (1) è il più grande. Quindi,

G-test usando i genotipi chiamati per la mappatura delle associazioni

Negli studi di associazione, gli SNP che mostrano differenze significative nella frequenza degli alleli tra casi e controlli sono associati al fenotipo di interesse. La mappatura delle associazioni può essere eseguita utilizzando i dati provenienti da studi di sequenziamento di nuova generazione. Per prima cosa discutiamo gli approcci che richiedono la chiamata di singoli genotipi e quindi eseguiamo un test per l’associazione utilizzando i genotipi chiamati. In questo approccio, viene prima chiamato un genotipo per ogni individuo. I genotipi possono essere filtrati o non filtrati. Assumendo l’indipendenza tra individui e HWE, una tabella di contingenza 2 × 2 può essere costruita contando il numero di alleli maggiori e minori sia nei casi che nei controlli. Questo porta al ben noto test del rapporto di verosimiglianza per l’indipendenza, il G-test:

(5)

dove O k,h è la frequenza osservata in una cella, e E k,h è la frequenza prevista nell’ipotesi in cui l’allele frequenza è la stessa tra casi e controlli. Il noto test chi-quadrato di Pearson è asintoticamente equivalente al test G. Se la tabella è generata da veri genotipi, allora la statistica G segue asintoticamente una distribuzione chi-quadrata con 1 grado di libertà (χ2(1)). Tuttavia, nei nostri studi, costruiamo la statistica G usando genotipi” chiamati”, quindi HWE potrebbe non reggere a causa di over – e under-calling di eterozigoti. Inoltre, la costruzione della statistica del test contando i genotipi” chiamati “anziché i genotipi “osservati” probabilmente introduce una maggiore variabilità. Pertanto, la teoria statistica potrebbe non essere più valida. Si noti che quando un genotipo non viene chiamato per un determinato individuo, i dati sono considerati mancanti e non sono inclusi nella tabella 2 × 2.

Test del rapporto di verosimiglianza contabilizzazione dell’incertezza nei genotipi osservati per la mappatura delle associazioni

Invece di chiamare i genotipi, il quadro di verosimiglianza consente l’incertezza nei genotipi e nei test in ogni sito j se la frequenza dell’allele è la stessa tra casi e controlli. Formalmente, calcoliamo la probabilità delle ipotesi H O : pj ,1 = pj ,2(= pj ,0) e H A : pj, 1 p pj, 2, dove pj ,1 e pj, 2 sono rispettivamente le MAF nei casi e nei controlli.

Supponendo che gli alleli minori (m) e maggiori (M) siano noti, la probabilità della frequenza dell’allele minore può essere calcolata come descritto nell’equazione 2 e la statistica del test del rapporto di verosimiglianza è calcolata come:

(6)

dove e sono i dati osservati per casi e controlli, rispettivamente, e e sono i MLEs di MAFs in casi e controlli, rispettivamente.

Se l’allele minore è sconosciuto, la probabilità sotto l’ipotesi nulla viene calcolata come nell’equazione 3 e la statistica LRT viene modificata come:

(7)

dove dj è il dato osservato per entrambi i casi e controlli, e è la frequenza allele sotto l’ipotesi nulla. Altre notazioni sono le stesse dell’equazione 6.

Stima della MAF nei dati simulati

Confrontiamo le stime della frequenza degli alleli sui dati simulati utilizzando genotipi veri (True), chiamati genotipi senza filtraggio (Chiamata NF), chiamati genotipi con filtraggio (f = 1; Chiamata F) e il metodo di massima verosimiglianza (ML). Per gli SNP rari, il tipo di allele minore spesso non è evidente. Quando si chiamano genotipi, si presume che il secondo nucleotide più comune sia l’allele minore. Il metodo ML incorpora direttamente l’incertezza nella determinazione dell’allele minore e, se non diversamente indicato, vengono mostrati i risultati utilizzando il metodo dell’allele minore sconosciuto (Equazione 3). Si noti che il metodo sconosciuto minor allele ML esegue in modo simile al metodo noto minor allele ML, ma il primo migliore per SNPS molto rari (file aggiuntivo 1).

Per prima cosa abbiamo valutato quanto bene i diversi approcci fossero in grado di stimare la MAF in 200 individui in una gamma di profondità di sequenziamento per 1.000 SNPS con una MAF vera del 5%. La figura 1 mostra i riquadri delle distribuzioni delle MAF stimate utilizzando i quattro diversi approcci. Come previsto, per dati di copertura più elevati, come una profondità individuale di 12×, tutti i metodi funzionano così come quando i genotipi sono noti con certezza (True). Tuttavia, quando la profondità diminuisce, le stime della MAF ottenute chiamando prima i genotipi diventano distorte. Ad esempio, la MAF mediana stimata utilizzando il metodo Call F è 5.3% a 6× copertura ed è 12.5% a 2×. La ragione del pregiudizio verso l’alto è che diventa più difficile chiamare eterozigoti poiché i veri eterozigoti spesso sembrano errori di sequenziamento. Pertanto, più eterozigoti rispetto agli omozigoti minori tendono ad avere genotipi mancanti. Tuttavia, il pregiudizio generale nelle stime MAF dai genotipi chiamati non è sempre in una direzione (dati non mostrati). È interessante notare che il bias sembra essere peggiore per il metodo Call F rispetto al metodo Call NF. Questo modello può sembrare contro-intuitivo poiché filtrare le chiamate genotipo sembrerebbe diminuire la probabilità di chiamare un errore di sequenziamento un eterozigote. Tuttavia, il metodo Call F si traduce anche in una maggiore quantità di dati mancanti poiché molti omozigoti per l’allele principale non verranno chiamati a causa di errori di sequenziamento. Pertanto, in questo caso, chiamare genotipi senza filtraggio sembra essere la strategia migliore rispetto al filtraggio dei genotipi quando si cerca di stimare la MAF.

Stime della frequenza degli alleli in siti con una vera MAF del 5% per diverse profondità di copertura. Ad ogni profondità, 1.000 siti sono stati simulati utilizzando 200 individui, e in ogni sito, una stima della frequenza allele è calcolata utilizzando: (1) genotipi veri (True); (2) chiamati genotipi senza filtraggio (Chiamata NF); (3) chiamati genotipi con filtraggio (Chiamata F); e (4) il metodo di massima verosimiglianza (ML). Per ulteriori dettagli sui metodi di stima, vedere Metodi.

I risultati sono drammaticamente diversi per il nuovo metodo ML. Questo metodo fornisce stime imparziali della MAF (mediana di ~4,9%) su una gamma di profondità. Anche a 2×, le stime mostrano solo una varianza leggermente più grande di quelle basate sui veri genotipi.

Abbiamo anche confrontato l’errore medio quadrato stimato (MSE; Expectation () delle diverse stime della MAF in un intervallo di profondità di sequenziamento (Figura 2). Il metodo ML ha un MSE inferiore rispetto ai metodi di chiamata con 50 o 200 individui. In particolare, l’MSE calcolato in base al metodo Call F è molto più alto di quelli degli altri metodi, specialmente quando la profondità diminuisce. L’MSE delle stime della MAF basate sui genotipi reali riflette il limite inferiore dell’MSE e non è costante attraverso le profondità a causa della varianza di campionamento e di una dimensione del campione finita. Utilizzando 50 individui, l’MSE si avvicina a 0,0005 con profondità crescente e quando si utilizza una dimensione del campione di 200 individui, si avvicina a 0,0013 con profondità crescente.

Media squred error (MSE; Atteso ) di quattro diversi tipi di frequenza allele estimatori per le diverse dimensioni del campione (a sinistra e a destra del pannello) e la profondità della copertura (asse x). Ad ogni profondità, MSE è stato calcolato dalle stime di frequenza allele effettuate utilizzando quattro diversi metodi: True, Chiamata NF, Chiamata F, e ML (per i dettagli dei metodi, vedere la didascalia della Figura 1).

Nel complesso, il nuovo metodo ML esegue i metodi di chiamata genotipo.

Stima di una distribuzione di MAFs da dati simulati

Esaminiamo successivamente come i diversi approcci di stima si sono comportati nella stima della proporzione di SNP a frequenze diverse nella popolazione (simile allo spettro di frequenza del sito ma basato sulla frequenza dell’allele di popolazione anziché sulla frequenza del campione). Qui abbiamo simulato 20.000 SNP in cui la distribuzione delle MAF vere seguiva la distribuzione stazionaria standard per una dimensione effettiva della popolazione di 10.000 (vedi Metodi). Si noti che in pratica, tuttavia, è molto difficile distinguere un SNP molto raro da un errore di sequenziamento. Pertanto, a scopo di confronto con i dati reali, abbiamo scartato SNP con MAF stimato inferiore al 2%. La figura 3 mostra la percentuale di SNP che cadono in ogni bin di frequenza diversa dopo aver escluso quegli SNP con MAF stimato < 2%.

la Distribuzione delle frequenze alleliche di SNPs simulato assumendo standard stazionario distribuzione delle frequenze alleliche. Ad ogni profondità (ogni pannello), sono stati simulati 20.000 SNP e per ogni SNP, le stime della MAF sono state ottenute utilizzando quattro metodi diversi (vedere la didascalia della Figura 1). Quindi, per ogni metodo (ogni colore), vengono utilizzati solo siti con frequenze allele stimate > 2% per generare ciascun istogramma (asse x).

Come previsto, con un’elevata profondità di copertura, ad esempio 10× per individuo, tutti i metodi forniscono distribuzioni MAF stimate simili alla distribuzione prevista basata sui veri genotipi (Figura 3). Con una profondità di copertura inferiore, ad esempio inferiore a 4× per individuo, le distribuzioni di MAF ottenute con i metodi di chiamata del genotipo si discostano significativamente dalla distribuzione MAF prevista basata su genotipi veri (Figura 3). In particolare, questi metodi sovrastimano la percentuale di SNP a bassa frequenza. Ad esempio, la percentuale prevista di SNP nel secondo bin (MAF stimato tra il 2-4%) è del 18%. La proporzione corrispondente basata sul metodo Call NF a una profondità di 4× è del 26%, che è 1,4 volte superiore al previsto. La sovrastima della proporzione di SNP a bassa frequenza si verifica a causa della confusione degli errori di sequenziamento con veri eterozigoti, che si traduce in overcalling genotipi eterozigoti. L’entità di questa inflazione differisce tra diversi tagli di filtraggio, ma un taglio più grande non aumenta o diminuisce necessariamente l’inflazione.

L’immagine è completamente diversa per il metodo ML. La distribuzione stimata di MAF ottenuta dal nuovo metodo ML segue da vicino la vera distribuzione anche con basse profondità di copertura. Qui non c’è quasi nessun eccesso di SNP a bassa frequenza. A una profondità di 4×, la proporzione di SNP nel secondo contenitore dell’istogramma è del 18,4%, che è molto vicina alla proporzione prevista (18%). Pertanto, stime più affidabili dello spettro di frequenza possono essere fatte da dati a bassa copertura utilizzando il nostro approccio di verosimiglianza piuttosto che utilizzando gli approcci di chiamata del genotipo.

Mappatura dell’associazione nei dati simulati

Confrontiamo le prestazioni dei metodi che trattano i genotipi inferiti come veri genotipi nei test di associazione (usando un G-test) con il nostro test del rapporto di verosimiglianza (LRT) che spiega l’incertezza nei genotipi. Esaminiamo la distribuzione del test-statistica sotto l’ipotesi nulla di nessuna differenza di frequenza allele tra casi e controlli. Confrontiamo anche la potenza dei diversi approcci.

Con campioni di dimensioni ragionevolmente grandi, la teoria asintotica standard suggerisce che sotto l’ipotesi nulla sia la statistica G che la statistica LRT seguono una distribuzione chi-quadrata con un grado di libertà (χ2(1)). Pertanto, abbiamo confrontato la distribuzione null della statistica G calcolata in base ai metodi di chiamata e la statistica LRT con la distribuzione χ2 (1) utilizzando QQ-plots (Figura 4). Abbiamo simulato 5.000 SNP in una varietà di profondità di sequenziamento in 500 casi e controlli in cui la MAF utilizzata per simulare i genotipi era del 5% in entrambi i casi e nei controlli. La distribuzione della statistica G calcolata utilizzando i veri genotipi mostra un’ottima corrispondenza con una distribuzione χ2(1). Tuttavia, la distribuzione della statistica G calcolata sulla base dei genotipi chiamati si discosta sostanzialmente da una distribuzione χ2 (1). Chiamare i genotipi e quindi trattare quei genotipi come accurati produce un vasto eccesso di segnali falsi positivi se i valori p sono calcolati usando una distribuzione χ2(1). Ad esempio, a una profondità di 2×, l ‘ 11% degli SNP aveva un valore p inferiore al 5%, rispetto al previsto 5%. L’effetto è causato da un aumento della varianza, dovuto al superamento degli omozigoti come eterozigoti, nel test allelico qui utilizzato per rilevare l’associazione. I test genotipici come Armitage trend test, che sono robusti alle deviazioni dall’equilibrio Hardy-Weinberg, non mostrano un aumento simile nel tasso di falsi positivi (file aggiuntivo 2). Coerentemente con questa osservazione, il filtraggio dei genotipi chiamati si traduce in una diminuzione della frazione di test significativi quando si utilizza il G-test, sebbene il filtraggio non risolva completamente il problema. D’altra parte, la statistica LRT mostra solo un leggero allontanamento da una distribuzione χ2(1) per profondità di copertura 2× o 5×.

QQ-grafici che confrontano la distribuzione nulla della statistica di test di interesse con una distribuzione χ2(1). Ogni colonna corrisponde a una diversa statistica di test: (1) G-statistica calcolata utilizzando i genotipi veri (True); (2) G-statistica calcolata utilizzando i genotipi chiamati senza filtraggio (Chiamata NF); (3) G-statistica calcolata usando genotipi chiamati con filtraggio (chiamata F); e (4) la statistica del test del rapporto di verosimiglianza con allele minore sconosciuto (LRT). Supponendo 500 casi e 500 controlli, sotto l’ipotesi nulla, un insieme di 5.000 siti sono stati simulati con una MAF del 5% con una profondità di sequenziamento di 2× (pannelli superiori) e 5× (pannelli inferiori). Il fattore “Inflazione” è mostrato nell’angolo in alto a sinistra di ogni figura.

Abbiamo anche generato curve ROC (Receiver Operating Characteristic) per ciascuno dei diversi test di associazione. Queste curve mostrano la potenza del test a diversi tassi di falsi positivi. Poiché le distribuzioni di alcune delle statistiche di test non seguono la distribuzione χ2(1) sotto l’ipotesi nulla, per fare un confronto equo, abbiamo ottenuto il valore critico per ogni tasso falso positivo basato sulla distribuzione nulla empirica. La potenza viene calcolata come la frazione di loci di malattia simulati che hanno una statistica superiore al valore critico. Nel complesso, troviamo che LRT si comporta meglio del G-test basato su entrambi i metodi di chiamata genotipo (Figura 5). Ad esempio, con un tasso di falsi positivi del 5% e con una profondità di sequenziamento di 5×, la potenza per rilevare un locus di malattia con una MAF dell ‘ 1% e un rischio relativo (RR) di 2 è del 51% con LRT, ma la potenza scende al 33% usando il metodo di chiamata senza filtraggio e al 34% usando il metodo di chiamata con filtraggio. In particolare, a bassa profondità, il G-test applicato ai genotipi chiamati con filtraggio si comporta molto male (colonna più a sinistra in Figura 5). Se confrontiamo la potenza dell’LRT con l’Armitage trend test usando i genotipi chiamati, scopriamo che l’LRT ha anche una potenza superiore rispetto all’Armitage trend test (file aggiuntivo 3). Ciò suggerisce che se si desidera utilizzare i genotipi chiamati, filtrarli in base alla fiducia delle chiamate può comportare una perdita di potenza.

Curve caratteristiche di funzionamento del ricevitore (ROC) di quattro prove di associazione. Per la definizione delle quattro statistiche, vedere la didascalia della Figura 4. Supponendo 500 casi e 500 controlli, un insieme di 20.000 siti sono stati simulati sotto il null e sotto l’alternativa a singole profondità di sequenziamento di 2×, 5× e 10× (tre colonne). Ad ogni tasso di falsi positivi (asse x), il valore critico corrispondente è stato calcolato utilizzando la distribuzione nulla empirica. Il vero tasso positivo (potenza; asse y) è stato ottenuto calcolando la frazione dei siti causali con statistiche di test che superano il valore critico.

Applicazione a dati reali

Abbiamo analizzato 200 exomi dai controlli per uno studio di associazione di malattie che sono stati sequenziati utilizzando la tecnologia Illumina a una profondità individuale di 8× . Abbiamo usato le probabilità di genotipo generate dal programma “SOAPsnp” per la nostra inferenza. Per ulteriori dettagli, vedere Metodi.

In primo luogo, abbiamo esplorato l’accuratezza delle stime della MAF dai dati di sequenziamento di prossima generazione per 50 SNPS confrontandoli con le MAF stimate dai dati del genotipo Sequenom. Sia le stime utilizzando il metodo ML e il metodo di chiamata genotipo senza filtraggio sono altamente correlati con le stime fatte dai dati di genotipo Sequenom (cioè, una piccola differenza standardizzata tra le due stime in Figura 6). Tuttavia, le stime basate sulla chiamata al genotipo con filtraggio mostrano una scarsa corrispondenza con le frequenze stimate dai dati del genotipo Sequenom, specialmente quando la profondità di sequenziamento è bassa. È interessante notare che esiste un SNP in cui la MAF stimata dai dati di resequencing è molto diversa dalla stima ottenuta dai dati del genotipo Sequenom, anche se la profondità di sequenziamento è molto alta (14×). In particolare, la MAF stimata dai dati del genotipo Sequenom è del 22,5%, ma è del 17,2% se stimata utilizzando l’approccio ML. L’esame individuale mostra che in molti individui, il genotipo altamente supportato basato sui dati di sequenziamento differisce dai genotipi Sequenom. Dato che questo SNP è coperto da molte letture in questi individui e che le basi di lettura osservate hanno punteggi di alta qualità (>Q20), è probabile che la differenza sia dovuta a errori di genotipizzazione Sequenom. Si noti che ci sono un paio di SNP in cui le MAF stimate dall’approccio di chiamata genotipo senza filtraggio sembrano corrispondere meglio alle MAF stimate dalla genotipizzazione Sequenom rispetto alle stime dell’approccio ML. Ad esempio, a un SNP la MAF stimata è del 25,7% dai dati del genotipo Sequenom, 25.9% dal metodo di chiamata genotipo senza filtraggio e 27,2% dal metodo ML. Tuttavia, l’ispezione individuale rivela che ci sono alcuni individui per i quali il genotipo chiamato dai dati di sequenziamento differisce dal genotipo Sequenom. In questi casi, gli errori nei genotipi chiamati annullati, dando l’apparenza di una migliore corrispondenza con i dati del genotipo Sequenom. Pertanto, per questi SNP, è difficile dire quale metodo esegue meglio.

Stime della frequenza degli alleli calcolate da 200 individui utilizzando i dati di sequenziamento di prossima generazione rispetto ai dati del genotipo Sequenom. In ogni sito, solo gli individui che hanno sia i dati di genotipo Sequenom che i dati di sequenziamento sono stati utilizzati per la stima della frequenza dell’allele. Per i dati di sequenziamento, le stime di MAF sono state ottenute utilizzando tre diversi metodi (Chiamata NF; Chiamata F; e ML). Il standardizzate differenza per ogni stima è stata calcolata come , dove e sono stimati MAFs dai dati di sequenziamento e di Sequenom dati di genotipo, rispettivamente, e n è il numero di individui utilizzati per la stima. Ogni sito è classificato in uno dei quattro bidoni in base alla profondità media individuale di copertura (colore): inferiore a 4×, superiore a 4× ma inferiore a 8×, superiore a 8× ma inferiore a 16× e superiore a 16×.

Abbiamo quindi esaminato la distribuzione di MAFs calcolata utilizzando diversi approcci in una gamma di profondità di sequenziamento dai nostri dati di sequenziamento exome di prossima generazione (Figura 7). Abbiamo scartato gli SNP con MAF stimato < 2% poiché è difficile distinguere questi SNP a bassissima frequenza dagli errori di sequenziamento in questo set di dati. Abbiamo inoltre rimosso i siti in cui vi era una differenza significativa (valore p inferiore a 10-5 utilizzando un rank-sum-test ) nel punteggio di qualità delle basi di lettura tra gli alleli minori e maggiori. È probabile che questi siti siano SNP artificiali che possono verificarsi a causa di una mappatura errata o di pregiudizi sconosciuti introdotti durante la procedura sperimentale. Poi abbiamo classificato ogni sito in bidoni in base alla profondità di copertura. Il numero di SNP in ciascun contenitore è mostrato nella Tabella 1. Quando la profondità media è inferiore a 9×, le distribuzioni delle MAF stimate basate sui metodi di chiamata del genotipo sono molto diverse da quella basata sul metodo ML. In particolare, gli approcci di chiamata del genotipo danno origine a un grande eccesso di SNP a bassa frequenza (MAF tra il 2% e il 4%). Questo modello rispecchia ciò che è stato visto nei nostri studi di simulazione (Figura 3). Inoltre, per i metodi di chiamata del genotipo, la distribuzione di frequenza dell’allele cambia drammaticamente mentre la profondità di sequenziamento cambia. Pertanto, come discusso in precedenza, quando la profondità non è molto alta, i metodi di chiamata genotipizzazione sono suscettibili di includere un sacco di falsi SNP che sono errori di sequenziamento. Questi errori appaiono come un eccesso di SNP a bassa frequenza nella distribuzione di frequenza. La distribuzione basata sul metodo ML è più stabile attraverso profondità, ma c’è ancora un eccesso di SNPS con bassa frequenza allele con profondità inferiore a 9×rispetto alla proporzione di SNPS a bassa frequenza a profondità maggiori.

Distribuzione della frequenza dell’allele minore stimata dagli esomi di 200 individui sequenziati. Per ogni sito, la frequenza dell’allele minore è stata stimata utilizzando quattro metodi diversi: (1) il metodo ML con sconosciuto allele minore, (2) il metodo ML con un noto o fisso allele minore, (3) chiamata genotipi senza filtraggio (Chiamata NF), e (4) chiamata genotipi con filtraggio (Chiamata F). Ogni sito è classificato in bidoni in base alla profondità di copertura. Inoltre, in ogni istogramma, i siti con MAF stimato inferiore al 2% non sono considerati. Per il numero di SNP utilizzati per questa analisi, vedere Tabella 1.

Infine, abbiamo usato questi dati exome-resequencing per simulare uno studio di associazione case-control. Per esaminare la distribuzione delle statistiche dei test di associazione sotto l’ipotesi nulla, abbiamo assegnato casualmente 100 individui a un gruppo di casi e gli altri 100 al gruppo di controllo. Per tutti gli SNP sul cromosoma 2 con stime MAF > 2% (basato sul metodo sconosciuto ML allele minore), abbiamo testato per le differenze di frequenza allele tra i casi e controlli calcolando la G-statistica utilizzando genotipi chiamati sia con e senza filtraggio, nonché la statistica LRT. La figura 8 mostra i grafici QQ che confrontano le distribuzioni delle statistiche di prova con la distribuzione standard χ2(1). Come visto negli studi di simulazione, la distribuzione nulla della statistica G calcolata quando si chiamano genotipi senza filtraggio si discosta sostanzialmente dalla distribuzione χ2(1). Tuttavia, la distribuzione nulla della statistica LRT segue da vicino la distribuzione χ2(1). Il fattore di inflazione è 1.01, il che implica che la statistica LRT si comporta bene quando applicata ai dati reali.

QQ-plot confrontando il test di associazione statistica per l’allele differenze di frequenza tra 100 casi e 100 controlli di un χ2(1) distribuzione. I fenotipi sono stati assegnati in modo casuale a indivdiduals nel set di dati di resequenziamento dell’esoma in modo tale che ci siano 100 casi e 100 controlli. Per ogni sito sono state calcolate tre statistiche: la G-statistica utilizzando genotipi chiamati senza filtraggio (Chiamata NF), la G-statica utilizzando genotipi chiamati con filtraggio (Chiamata F), e la statistica LRT. Per ridurre al minimo l’inclusione di falsi SNP, i siti con stime ML MAF inferiori al 2% vengono scartati. Ai fini della visualizzazione, vengono mostrati i risultati dei siti sul cromosoma 2. Si noti che il fattore di inflazione è mostrato nell’angolo in alto a sinistra di ogni QQ-plot.