O menor alelo é o menos frequente alelo na população em uma variável local. Primeiramente, descrevemos duas abordagens principais para estimar a frequência alélica menor (MAF) em um determinado local no genoma. A primeira abordagem envolve inferir genótipos individuais e tratar esses genótipos inferidos como sendo completamente precisos ao estimar o MAF. Examinamos então o desempenho de um quadro de probabilidade que leva diretamente em conta a incerteza na atribuição de genótipos. Ao longo do nosso trabalho, assumimos que todos os sites segregantes são bialélicos.

A estimativa da MAF a partir dos genótipos chamados

uma forma de estimar a MAF a partir dos dados de sequenciação da próxima geração é primeiro chamar um genótipo para cada indivíduo usando dados de sequenciação, e depois usar esses genótipos como se fossem os verdadeiros. Esta foi a abordagem tradicionalmente usada para dados de genótipo e sequenciamento de Sanger. Não está claro o quão bem irá funcionar quando aplicado aos dados de sequenciação da próxima geração.

pode utilizar-se uma abordagem de probabilidade máxima para inferir o genótipo para cada indivíduo a partir dos dados de sequenciação da próxima geração. Em cada site j, para cada indivíduo i, a probabilidade para cada um dos três genótipos possíveis (assumindo que conhecemos o alelo menor) é dada como:

, onde D i,j é o observado de dados de sequenciamento em indivíduo i no local j, g, i , j ∈ {0, 1, 2} é o menor número de alelos contidas no genótipo de cada indivíduo, e e controle para a seqüência de erros e de leitura da base de dados de qualidades, respectivamente. Os dados de sequenciamento observados para cada indivíduo podem ser considerados como o alinhamento de leituras no site j levando em conta os escores de qualidade de leitura. Isto é representado como a probabilidade do genótipo e é encontrado no genótipo probability file (GLF), que é produzido em muitos programas que analisam dados de sequenciação de próxima geração, como SOAPsnp e MAQ .para atribuir um genótipo a um determinado indivíduo, a probabilidade de cada um dos três genótipos possíveis pode ser calculada para o indivíduo. O genótipo com maior probabilidade pode então ser atribuído. No entanto, os investigadores preferem frequentemente um critério de chamada mais rigoroso e não atribuirão um genótipo a um indivíduo, a menos que o genótipo mais provável seja substancialmente mais provável do que o segundo mais provável. Aqui os três genótipos possíveis são ordenados pelo seu likelihoods: , onde g(k)corresponde ao genótipo com a th maior probabilidade. Com um determinado limiar f, pode-se chamar o genótipo g(1) se . Caso contrário, um genótipo não é chamado e o genótipo do indivíduo é considerado inexistente. Um valor limiar comum de f é 1, indicando que o genótipo mais provável é pelo menos 10 vezes mais provável do que o segundo mais provável. Note que este tipo de filtragem pode resultar em maior confiança para o genótipo “chamado”, mas também resulta em mais dados em falta.

estimador máximo de probabilidade da frequência do alelo

em vez de estimar a MAF a partir dos genótipos chamados, um método de probabilidade máxima (ML) introduzido por Kim et al. (see also Lynch for a similar approach) directly estimates MAF and takes genoty into account. Especificamente, dado um alelo menor, a probabilidade de observar os dados de sequência em cada indivíduo i é obtida somando as probabilidades correspondentes aos três genótipos possíveis.suponha que os três genótipos de likelihoods definidos na equação 1 estão disponíveis. Utilizando a mesma notação acima, D j E p j sejam os dados de sequenciação observados no site j e no MAF correspondente, respectivamente. A probabilidade do genótipo dado que a frequência alélica menor pode ser calculada assumindo o equilíbrio de Hardy-Weinberg (HWE). Em seguida, supondo independência entre os indivíduos, a probabilidade de o MAF neste locus é um produto de todos a probabilidade calculada através de todos os N indivíduos:

O ML estimativa de p j pode ser calculado directamente, maximizando a probabilidade de um restrito espaço de parâmetros usando o Broyden-Fletcher-Goldfarb-Shanno (BFGS) método ou usando a expectativa de maximização (EM) o algoritmo . Ao usar o algoritmo EM, a expectativa posterior de um genótipo é computada para cada indivíduo, e a média desses posteriores é atualizada repetidamente. A nossa implementação do BFGS foi mais rápida que o algoritmo EM. Por exemplo, para obter estimativas de 100.000 sites, BFGS levou ~16 segundos, mas EM levou ~100 segundos. No entanto, a diferença de velocidade pode ser específica para a implementação. No nosso caso, para ambos os métodos, paramos de atualizar os parâmetros quando o aumento da probabilidade foi inferior a 0,001.

estimador de probabilidade máxima com Alelo menor incerto

na prática, muitas vezes o segundo nucleótido mais comum entre os indivíduos pode ser usado como Alelo menor. No entanto, para SNPs raros (por exemplo, MAF < 1%), é difícil determinar qual alelo é o alelo menor, uma vez que todos os quatro nucleótidos podem aparecer em algumas leituras devido a erros de sequenciação. Para lidar com esta situação, descrevemos agora um quadro de probabilidade que tem em conta a incerteza na determinação do alelo menor.suponha que para o site j conhecemos o principal Alelo M. Note que decidir qual de dois alelos comuns é provável que seja o principal não é importante, uma vez que estamos principalmente preocupados com a estimativa das frequências em SNPs raros. Além disso, para alelos com frequências intermediárias (cerca de 50%), a distinção entre alelos maiores e menores é menos importante. Atribuir os outros três nucleótidos Não principais m1, m2 e m3. A probabilidade introduzida na equação 2 assume um alelo maior fixo M e um alelo menor fixo m. Portanto, para permitir a incerteza na designação do alelo menor, a função de probabilidade pode ser modificada como:

além disso, supondo-se que nenhum dos três possíveis pequenas alelos é igualmente provável, obtém-se:

, onde . Uma vez que pode ser muito pequeno com grandes conjuntos de dados (por exemplo, com muitos indivíduos), é útil calcular a probabilidade na escala log. Ordene os três log-likelihoods condicionais a (l (1), l(2), l(3)), onde l(1) é o maior. Em seguida,

G-teste, usando o chamado de genótipos para a associação de mapeamento

Em estudos de associação, SNPs, mostrando diferenças significativas em alelo frequência entre os casos e controles são disse a ser associado com o fenótipo de interesse. O mapeamento de associação pode ser realizado utilizando dados de estudos de sequenciação de próxima geração. Primeiro discutimos abordagens que exigem chamar genótipos individuais e, em seguida, realizar um teste para associação usando os chamados genótipos. Nesta abordagem, um genótipo é primeiramente chamado para cada indivíduo. Os genótipos podem ser filtrados ou não filtrados. Assumindo a independência entre indivíduos e HWE, uma tabela de contingência 2 × 2 pode ser construída contando o número de alelos maiores e menores em ambos os casos e controles. Isso leva ao conhecido teste de razão de verossimilhança para a independência, o G-teste:

onde S k,h é a freqüência observada em uma célula, e e k,h é a frequência esperado sob a hipótese nula, em que o alelo frequência é a mesma entre casos e controles. O conhecido teste de Qui-quadrado de Pearson é assintoticamente equivalente ao teste de G -. Se a tabela é gerada a partir de genótipos verdadeiros, então a estatística G segue assintoticamente uma distribuição qui-quadrado com 1 grau de liberdade (χ2(1)). No entanto, em nossos estudos, nós construímos a estatística G usando genótipos “chamados”, assim HWE pode não se manter devido ao excesso e sub-chamado de heterozigotos. Além disso, a construção da Estatística do ensaio, contando os genótipos “denominados” em vez dos genótipos “observados”, provavelmente introduz uma variabilidade adicional. Portanto, a teoria estatística pode já não ser válida. Note – se que quando um genótipo não é chamado para um determinado indivíduo, os dados são considerados omissos e não estão incluídos na tabela 2 × 2.o quadro de probabilidade permite a incerteza nos genótipos e ensaios em cada local j se a frequência do alelo é a mesma entre os casos e os controlos. Formalmente, calculamos a probabilidade das hipóteses H O: p j, 1 = p j, 2 (=p j, 0) e H A : p j, 1 ≠ p j, 2, em que p j ,1 e p j ,2 são os MAF nos casos e controlos, respectivamente.

assumindo que são conhecidos alelos menores (m) e maiores (m), a probabilidade da frequência do alelo menor pode ser calculada como descrito na equação 2, e a estatística da razão de probabilidade do ensaio é calculada como:

, onde e são os dados observados para casos e controles, respectivamente, e e são os MLEs do MAFs em casos e controles, respectivamente.

Se o alelo menor for desconhecido, a probabilidade sob a hipótese nula é calculada como na equação 3, e a estatística da GTAI é modificada como:

, onde D j é a dados observados para ambos os casos e controles, e é o alelo frequência sob a hipótese nula. Outras notações são as mesmas que na equação 6.

estimar a MAF em dados simulados

comparamos as estimativas da frequência do alelo em dados simulados usando genotipos verdadeiros (True), chamados genotipos sem qualquer filtragem (Call NF), chamados genotipos com filtragem (f = 1; Call F), e o método de máxima probabilidade (ML). Para SNPs raros, o tipo Alelo menor muitas vezes não é aparente. Ao chamar genótipos, o segundo nucleótido mais comum é assumido como sendo o alelo menor. O método ML incorpora directamente incerteza na determinação do alelo menor e, salvo indicação em contrário, são apresentados resultados utilizando o método desconhecido do alelo menor (equação 3). Note – se que o método conhecido do alelo ML menor é semelhante ao método conhecido do alelo ML menor, mas o primeiro é melhor para SNPs muito raros (arquivo adicional 1).

primeiro avaliamos o quão bem as diferentes abordagens foram capazes de estimar o MAF em 200 indivíduos através de uma gama de profundidades sequenciadoras para 1.000 SNPs com um verdadeiro MAF de 5%. A figura 1 mostra as fichas de caixa das distribuições dos MAF estimados utilizando as quatro abordagens diferentes. Como esperado, para dados de cobertura mais elevados, como uma profundidade individual de 12×, todos os métodos funcionam bem como quando os genótipos são conhecidos com certeza (verdadeiro). No entanto, quando a profundidade diminui, as estimativas da MAF obtidas pela primeira chamada de genótipos tornam-se tendenciosas. Por exemplo, a mediana MAF estimada utilizando o método Call F é de 5,3% a 6× cobertura e é de 12,5% a 2×. A razão para o viés ascendente é que se torna mais difícil chamar heterozigotos já que os verdadeiros heterozigotos muitas vezes se parecem com erros de sequenciação. Portanto, mais heterozigotos do que homozigotos menores tendem a ter genótipos em falta. No entanto, o viés geral nas estimativas de MAF de genótipos chamados nem sempre está em uma direção (dados não mostrados). Curiosamente, o viés parece ser pior para o método Call F do que o método Call NF. Este padrão pode parecer contra-intuitivo uma vez que filtrar as chamadas do genótipo parece diminuir a probabilidade de chamar um erro de sequenciação de um heterozigoto. No entanto, o método Call F também resulta em uma maior quantidade de dados em falta, uma vez que muitos homozigotes para o alelo principal não será chamado devido a erros de sequenciação. Assim, neste caso, chamar genótipos sem filtragem parece ser a melhor estratégia do que filtrar genótipos ao tentar estimar a MAF.

os resultados são dramaticamente diferentes para o novo método ML. Este método fornece estimativas imparciais do MAF (mediana de ~4,9%) em toda uma gama de profundidades. Mesmo a 2×, as estimativas mostram apenas uma variância ligeiramente maior do que as baseadas nos genótipos verdadeiros.

também comparámos a média estimada de erro ao quadrado (MSE; expectação () das diferentes estimativas do MAF numa gama de profundidades sequenciadoras (Figura 2). O método ML tem um MSE inferior aos métodos de chamada com 50 ou 200 indivíduos. In particular, the MSE computed based on the Call F method is much higher than those from the other methods especially when the depth decreases. O MSE das estimativas do MAF com base nos genótipos verdadeiros reflecte o limite inferior do MSE e não é constante nas profundidades devido à variância da amostragem e ao tamanho finito da amostra. Usando 50 indivíduos, o MSE aproxima-se de 0.0005 com profundidade crescente e ao usar um tamanho de amostra de 200 indivíduos, ele se aproxima de 0.0013 com profundidade crescente.

globalmente, o novo método de ML apresenta métodos de activação do genótipo mais eficazes.

estimar uma distribuição de Map a partir de dados simulados

em seguida, examinar como as diferentes abordagens de estimativa realizadas na estimativa da proporção de PNS em diferentes frequências na população (semelhante ao espectro de Frequência local, mas com base na frequência do alelo da população em vez da frequência da amostra). Aqui nós simulamos 20.000 SNPs, onde a distribuição dos MFS verdadeiros seguiu a distribuição estacionária padrão para uma população efetiva de 10 mil (ver Métodos). Note que na prática, no entanto, é muito difícil distinguir um SNP muito raro de um erro sequenciador. Portanto, para fins de comparação com dados reais, descartamos SNPs com MAF estimado em menos de 2%. A figura 3 mostra a proporção de PNS que caem em cada barra de frequências diferente, excluindo os PNS com uma estimativa MAF<2%.

Conforme o esperado, com uma elevada profundidade de cobertura, tais como 10× por indivíduo, todos os métodos fornecem estimado MAF distribuições que são semelhantes para a distribuição esperada com base na verdade genótipos (Figura 3). Com uma menor profundidade de cobertura, como menos de 4× por indivíduo, as distribuições de Map obtidas por métodos de chamada de genótipos divergem significativamente da distribuição prevista de MAP Com base em genótipos verdadeiros (Figura 3). Em particular, estes métodos sobrestimam a proporção de SNPs de baixa frequência. Por exemplo, a proporção esperada de SNPs no segundo bin (estimativa de MAF entre 2-4%) é de 18%. A proporção correspondente baseada no método de Call NF a uma profundidade de 4× é de 26%, o que é 1,4 vezes superior ao esperado. A superestimação da proporção de SNPs de baixa frequência ocorre devido à confusão de erros de sequenciação com verdadeiros heterozigotos, o que resulta em superlotação de genótipos heterozigóticos. A magnitude desta inflação difere em diferentes cortes de filtragem, mas um corte maior não necessariamente aumenta ou diminui a inflação.

A imagem é completamente diferente para o método ML. A distribuição estimada de MAF obtida a partir do novo método de ML segue de perto a distribuição verdadeira, mesmo com profundidades de cobertura pouco profundas. Aqui não há quase nenhum excesso de SNPs de baixa frequência. A uma profundidade de 4×, a proporção de SNPs no segundo bin do histograma é de 18,4%, o que é muito próximo da proporção esperada (18%). Assim, estimativas mais confiáveis do espectro de frequências podem ser feitas a partir de dados de baixa cobertura usando nossa abordagem de probabilidade do que usando as abordagens de chamada de genótipo.

mapeamento de associação em dados simulados

comparamos o desempenho dos métodos que tratam os genótipos inferidos como genótipos verdadeiros em testes de associação (utilizando um teste G) com o nosso teste da relação de probabilidade (GRT) que explica a incerteza nos genótipos. Examinamos a distribuição da estatística de teste sob a hipótese nula de nenhuma diferença de frequência alélica entre casos e controles. Também comparamos o poder das diferentes abordagens.

com tamanhos de amostra razoavelmente grandes, a teoria assintótica Padrão sugere que, sob a hipótese nula, tanto a estatística G-estatística quanto a LRT seguem uma distribuição qui-quadrada com um grau de liberdade (χ2(1)). Portanto, comparamos a distribuição nula da estatística G calculada com base em métodos de chamada, bem como a estatística da GTRT com a distribuição χ2(1) usando QQ-plots (Figura 4). Nós simulamos 5.000 SNPs em uma variedade de profundidades sequenciadas em 500 casos e controles onde o MAF usado para simular genotipos foi de 5% em ambos os casos e controles. A distribuição da estatística G calculada usando os genótipos verdadeiros mostra uma correspondência muito boa com uma distribuição χ2(1). No entanto, a distribuição da estatística G calculada com base nos chamados genótipos diverge substancialmente de uma distribuição χ2(1). Chamar genótipos e, em seguida, tratar esses genótipos como sendo precisos produz um grande excesso de sinais falsos positivos se os valores de p são computados usando uma distribuição χ2(1). Por exemplo, a uma profundidade de 2×, 11% dos PNS tinham um valor p inferior a 5%, em comparação com os 5% esperados. O efeito é causado por um aumento na variância, devido à superlotação homozigota como heterozigotos, no teste allelic usado aqui para detectar associação. Testes genotípicos, como o teste de tendência Armitage, que são robustos a desvios do equilíbrio Hardy-Weinberg, não mostram um aumento semelhante na taxa de falsos positivos (arquivo adicional 2). Consistente com esta observação, filtrar os genótipos chamados resulta em uma diminuição na fração de testes significativos ao usar o teste G, embora a filtragem não resolve completamente o problema. Por outro lado, a estatística da GTAI mostra apenas um ligeiro afastamento de uma distribuição χ2(1) Para 2× ou 5× profundidades de cobertura.

we also generated receiver operating characteristic (ROC) curves for each of the different association tests. Estas curvas mostram a potência do teste em diferentes taxas de falsos positivos. Uma vez que as distribuições de algumas das estatísticas de teste não seguem a distribuição χ2(1) sob a hipótese nula, para fazer uma comparação justa, obtivemos o valor crítico para cada taxa falsa positiva com base na distribuição empírica nula. A potência é calculada como a fração de loci de doença simulada que tem uma estatística que excede o valor crítico. No geral, achamos que a GTRT tem um desempenho melhor do que o G-test baseado em qualquer método de chamada de genótipo (Figura 5). Por exemplo, em um 5% taxa de falso-positivo e com uma sequência de profundidade de 5×, a capacidade para detectar uma doença locus com uma MAF de 1% e um risco relativo (RR) de 2 a 51% com o LRT, mas a potência cai para 33%, utilizando o método de chamada sem filtragem e 34% usando o método de chamada com a filtragem. Em particular, a baixa profundidade, o teste G aplicado aos chamados genótipos com filtragem tem um desempenho muito fraco (à esquerda A maior parte da coluna na Figura 5). Se compararmos o poder da GTAI com o teste de tendência do Armitage usando genótipos, descobrimos que a GTAI também tem maior poder do que o teste de tendência do Armitage (arquivo adicional 3). Isto sugere que se se deseja usar genótipos chamados, filtrando-os com base na confiança de chamada pode resultar em uma perda de energia.

aplicação a dados reais

analisámos 200 exomas de controlos para um estudo de associação de doenças que foram sequenciados utilizando tecnologia de iluminação a uma profundidade individual de 8× . Usámos o genótipo likelihoods gerado pelo programa” SOAPsnp ” para a nossa inferência. Para mais detalhes, veja métodos.em primeiro lugar, exploramos a precisão das estimativas da MAF a partir de dados sequenciais de próxima geração para 50 SNPs, comparando-os com os MAF estimados a partir de dados de genótipo Sequenom. Tanto as estimativas que utilizam o método ML como o método de chamada do genótipo sem filtragem estão altamente correlacionadas com as estimativas feitas a partir dos dados do genótipo Sequenom (ou seja, uma pequena diferença padronizada entre as duas estimativas na Figura 6). No entanto, as estimativas baseadas na chamada de genótipos com filtragem mostram pouca correspondência com as frequências estimadas a partir dos dados do genótipo Sequenom, especialmente quando a profundidade de sequenciamento é baixa. Curiosamente, há um SNP em que o MAF estimado a partir dos dados de resequenciação é muito diferente da estimativa obtida a partir dos dados do genótipo Sequenom, embora a profundidade de sequenciamento seja muito alta (14×). Especificamente, a MAF estimada a partir dos dados do genótipo de Sequenom é de 22, 5%, mas é de 17, 2% quando estimada utilizando a abordagem de ML. O exame individual mostra que em muitos indivíduos, o genótipo altamente suportado baseado nos dados de sequenciação difere dos genótipos de Sequenom. Dado que este SNP é coberto por muitas leituras nestes indivíduos e que as bases de leitura observadas têm pontuações de alta qualidade (>Q20), é provável que a diferença seja devido a erros de genotipagem Sequenom. Note – se que há um par de SNPs em que os Maps estimados do genótipo calling approach sem filtragem parecem corresponder melhor aos Maps estimados a partir da genotipagem Sequenom do que as estimativas da abordagem ML fazem. Por exemplo, em um SNP, o MAF estimado é de 25, 7% a partir dos dados do genótipo Sequenom, 25.9% do método de chamada do genótipo sem filtragem e 27, 2% do método ML. No entanto, a inspeção individual revela que existem alguns indivíduos para os quais o chamado genótipo dos dados de sequenciação difere do genótipo Sequenom. Nestes casos, os erros nos chamados genótipos foram cancelados, dando a aparência de uma melhor correspondência com os dados do genótipo Sequenom. Portanto, para estes SNPs, é difícil dizer qual método funciona melhor.

examinámos em seguida a distribuição de MAF calculada utilizando várias abordagens numa gama de profundidades sequenciais a partir dos nossos dados de sequenciação exome da próxima geração (Figura 7). Nós descartamos SNPs com estimativa MAF < 2% uma vez que é difícil distinguir estes SNPs de baixa frequência de erros de sequenciação neste conjunto de dados. Nós removemos ainda locais em que havia uma diferença significativa (valor p inferior a 10-5 usando um rank-sum-test ) na pontuação de qualidade das bases de leitura entre os alelos menores e maiores. Estes locais são susceptíveis de ser SNPs artificiais que podem ocorrer devido a mapeamento incorreto ou desvios desconhecidos introduzidos durante o procedimento experimental. Depois, classificamos cada local em contentores com base na profundidade da cobertura. O número de SNPs em cada caixa é mostrado na Tabela 1. Quando a profundidade média é inferior a 9×, as distribuições dos Map estimados com base nos métodos de chamada de genótipos são muito diferentes das do método de ML. Especificamente, as abordagens de chamada do genótipo dão origem a um grande excesso de SNPs de baixa frequência (MAF entre 2% e 4%). Este padrão espelha o que foi visto em nossos estudos de simulação (Figura 3). Além disso, para os métodos de chamada do genótipo, a distribuição de frequência do alelo muda dramaticamente como sequenciando mudanças de profundidade. Portanto, como discutido anteriormente, quando a profundidade não é muito alta, os métodos de chamada genotipagem são susceptíveis de incluir um monte de falsos SNPs que são erros de sequenciamento. Estes erros aparecem como um excesso de SNPs de baixa frequência na distribuição de frequência. A distribuição baseada no método ML é mais estável nas profundidades, mas continua a haver um excesso de SNPs com baixa frequência alélica com profundidade inferior a 9×em comparação com a proporção de SNPs de baixa frequência com profundidades superiores.

finalmente, usámos estes dados exome-resequenciadores para simular um estudo de associação de controlo de casos. Para examinar a distribuição das estatísticas de teste da Associação sob a hipótese nula, atribuímos aleatoriamente 100 indivíduos a um grupo de casos e os outros 100 ao grupo de controle. Para todos os SNPs do cromossoma 2 com estimativas de MAF > 2% (com base no método do alelo ml menor desconhecido), testámos diferenças de frequência entre casos e controlos computando a estatística de G usando genótipos chamados, tanto com e sem filtragem, como a estatística da GT. A figura 8 mostra as parcelas QQ comparando as distribuições das estatísticas de ensaio com a distribuição padrão χ2(1). Conforme observado em estudos de simulação, a distribuição nula da estatística G calculada ao chamar genótipos sem filtragem difere substancialmente da distribuição χ2(1). No entanto, a distribuição nula da estatística da GTL segue de perto a distribuição χ2(1). O Fator de inflação é 1.01, implicando que a estatística da GTL tem um bom desempenho quando aplicada a dados reais.