Comparação bayesiana de modelos com uma aplicação para o equilíbrio de Hardy-Weinberg usando o coeficiente de desequilíbrio

Bayesian comparison of models with an application to the Hardy-Weinberg equilibrium using the disequilibrium coefficient

ARTIGOS CIENTÍFICOS

BIOLOGIA

Comparação bayesiana de modelos com uma aplicação para o equilíbrio de Hardy-Weinberg usando o coeficiente de desequilíbrio

Bayesian comparison of models with an application to the Hardy-Weinberg equilibrium using the disequilibrium coefficient

Ricardo Luis dos Reis^I; Joel Augusto Muniz^I,1; Fabyano Fonseca e Silva^II; Thelma Sáfadi^I; Luiz Henrique de Aquino^I

^IPrograma de Pós-graduação em Estatística e Experimentação Agropecuária, Universidade Federal de Lavras (UFLA), Campus Universitário, CP 3037, 37200-000, Lavras, MG, Brasil. E-mail: [e-mail]

^IIDepartamento de Informática, Universidade Federal de Viçosa (UFV), Viçosa, MG, Brasil

RESUMO

O equilíbrio de Hardy-Weinberg é um dos principais assuntos estudados pela Genética de populações. Neste contexto, o presente trabalho aborda a análise e a comparação bayesiana de modelos utilizando o coeficiente de desequilíbrio (D_A). Para isso, realizou-se um estudo de simulação no qual as seguintes distribuições a priori foram consideradas: Dirichlet (modelo 1); beta - função degrau uniforme (modelo 2); uniforme - função degrau uniforme (modelo 3); e as prioris independentes uniformes (modelo 4). Exemplos de aplicação a dados reais de grupos raciais também são apresentados e discutidos. As amostras das distribuições marginais a posteriori para os parâmetros de interesse foram obtidas mediante o algoritmo Metropolis-Hastings, o qual foi implementado no software livre R. A convergência das cadeias geradas por este algoritmo foi monitorada pelos critérios de Geweke e Gelman & Rubin, os quais estão implementados no pacote BOA do R. Quanto às comparações entre os modelos, efetuadas por meio do fator de Bayes, observa-se que, para os dados simulados, o modelo 4 é o mais indicado para os casos de D_A=0,146, D_A=0,02 e D_A=-0,02 com n=200; o modelo 2 é o mais indicado para D_A=-0,02 e n=50 e o modelo 3 é o mais indicado para D_A=-0,02 e n=1000. Para os dados reais, em cada caso analisado, nota-se uma grande diferenciação na escolha de modelos, em que apenas o modelo 1 não é recomendado.

Palavras-chave: fator de Bayes, Genética de populações, simulação de dados.

ABSTRACT

One of the main subjects studied by population genetics is the Hardy-Weinberg equilibrium. In this context, this paper addresses the analysis and comparison of bayesian models used in its evaluation by the coefficient of disequilibrium. For this, it was carried out a simulation study in which the following prior distributions were considered: Dirichlet (model 1), beta - uniform step function (model 2), uniform - uniform step function (model 3) and independent uniform priors (model 4). Examples of application to real data for racial groups are presented and discussed. Samples from the marginal posterior distributions for parameters of interest were obtained by Metropolis-Hastings algorithm, which was implemented in the software R. The convergence of the chains generated by this algorithm was monitored by criteria of Geweke and Gelman & Rubin, which are implemented in the BOA package R. Regarding comparisons between models, performed using the Bayes factor, it was observed that model 4 is the most suitable for the cases of D_A=0.146, D_A=0.02 and D_A=-0.02 with n=200, the model 2 is the most suitable for D_A=-0.02 with n=50 and the model 3 is the most suitable for D_A=-0.02 and n=1000. For real data, in each case examined, there is a large difference in choice of models, where model 1 is the only one not recommended.

Key words: Bayes factor, population genetics, simulation data.

INTRODUÇÃO

Estudos referentes ao modo como os genes estão distribuídos nos indivíduos assumem grande importância para a obtenção de informações úteis ao estabelecimento de estratégias mais seguras na coleta e conservação da variabilidade genética. Em 1908, o matemático inglês Godfrey Harold Hardy e o médico alemão Wilhelm Weinberg chegaram independentemente e quase que simultaneamente à Lei do Equilíbrio de Hardy-Weinberg. Estes pesquisadores perceberam que, se não existissem fatores evolutivos atuando sobre uma população, as proporções alélicas permaneceriam inalteradas, e as proporções genotípicas atingiriam um equilíbrio estável, mostrando a mesma relação constante entre si ao longo do tempo. Portanto, considerando um gene com dois tipos de alelos, A e B, e definindo as proporções alélicas p_A e p_B=1-p_A, pode-se determinar as proporções genotípicas da população, como: P_AA = p²_A (proporção do genótipo homozigoto AA), P_AB=2p_Ap_B (proporção do genótipo heterozigoto AB) e P_BB=P²_B(proporção do genótipo homozigoto BB).

Um dos assuntos mais pesquisados na área da Genética de populações refere-se ao estudo das violações à Lei de Hardy-Weinberg. Neste caso, um dos parâmetros mais utilizados na avaliação desse desvio é o coeficiente de desequilíbrio DA (HERNÁNDEZ & WEIR, 1989), em que este expressa a relação entre as proporções alélicas e o coeficiente de endogamia de uma população. Portanto, o coeficiente de desequilíbrio mede estas discrepâncias entre as proporções genotípicas sob cruzamentos aleatórios e estas sob cruzamentos endogâmicos na população (WEIR, 1996). A endogamia é definida como um sistema em que os acasalamentos se dão entre indivíduos aparentados, ou seja, relacionados pela ascendência, afetando diretamente a diversidade genética da população (MUNIZ et al., 2008; MUNIZ et al., 2010).

As proporções genotípicas homozigotas e heterozigotas para o caso de dois alelos sob a hipótese de violação do modelo de Hardy-Weinberg são definidas por (HERNÁNDEZ & WEIR, 1989):

em que os limites de D_A são dados por ,que dependem das proporções alélicas.

Utilizando este modelo, AYRES & BALDING (1998) e ARMBORST (2005) aplicaram o método bayesiano na estimação dos parâmetros e concluíram que este possibilitou a incorporação dos efeitos de incerteza relativa aos parâmetros nuisance (parâmetros pelos quais não se tem interesse direto), isto é, as proporções alélicas. Assim, na comparação com os métodos frequentistas, o método bayesiano, utilizando distribuições a priori uniformes, apresentou os melhores resultados. SHOEMAKER et al. (1998) descreveram uma metodologia bayesiana para estudar o equilíbrio de Hardy-Weinberg considerando dois parâmetros, o coeficiente de desequilíbrio e o coeficiente de endogamia. Estes autores usaram três distribuições a priori para cada parâmetro (Dirichlet, beta - função degrau uniforme e uniforme - função degrau uniforme), mas não consideraram nenhuma forma específica de comparação entre estas.

O trabalho desenvolvido por REIS et al. (2008) utilizou a metodologia bayesiana na estimação do coeficiente de endogamia e da taxa de fecundação cruzada de uma população diploide por meio do modelo aleatório de COCKERHAM (1969) para frequências alélicas. Este trabalho propiciou resultados condizentes, validados pelo estudo de simulação de dados adotado. REIS et al. (2009) descreveram um método bayesiano para estudar o equilíbrio de Hardy-Weinberg através do coeficiente de endogamia. Neste trabalho, os autores analisaram vários modelos e concluíram que o melhor modelo é aquele que utiliza distribuições a priori Dirichlet.

Tendo em vista os aspectos apresentados, o objetivo do presente trabalho foi reescrever Hardy-Weinberg por meio do coeficiente de desequilíbrio. Objetivou-se também testar o método por meio de estudos de simulação de dados e aplicá-lo a um conjunto de dados reais.

MATERIAL E MÉTODOS

Atualmente, a inferência bayesiana é um método bastante utilizado em pesquisas genéticas para estimação do parâmetro DA. Esta é definida através do Teorema de Bayes, que associa um modelo relacionado aos dados (função de verossimilhança) com a distribuição a priori dos parâmetros, que são considerados aleatórios e, a partir daí, resume essas informações através da distribuição condicional dos parâmetros sobre os dados observados, a distribuição a posteriori (GELMAN et al., 2000). Portanto, a estimação do parâmetro é realizada a partir da distribuição a posteriori, e esta informação pode ser resumida pela média, moda, mediana ou pelos intervalos de credibilidade (PAULINO et al., 2003).

Neste trabalho, as distribuições a priori utilizadas foram as Dirichlet, com hiperparâmetros inteiros definidas como em que .

De (1),

A distribuição a priori conjunta beta-função degrau uniforme foi obtida por ,em que a distribuição a priori para , foi condicionada por uma distribuição beta com hiperparâmetros α e β, e a distribuição condicional a priori para D_A, dado p_A, foi determinada por uma função degrau uniforme sob cada um dos intervalos apresentados em SHOEMAKER et al. (1998). Dessa forma, a distribuição a priori conjunta é dada por:

A distribuição a priori conjunta uniforme - função degrau uniforme - é obtida por diferenciando da distribuição a priori conjunta beta - função degrau uniforme - apenas pela distribuição a priori para , condicionada por uma distribuição uniforme. Dessa forma, a distribuição a priori conjunta é dada por: . Considerando a independência entre os parâmetros e a falta de informação a priori, optou-se também pela utilização de uma distribuição uniforme para cada um dos parâmetros. Portanto, a distribuição a priori conjunta é dada por , em que

Para a definição da função de verossimilhança, considere que n₁, n₂ e n₃ representem as quantidades observadas de genótipos AA, AB e BB, respectivamente, em uma amostra de tamanho n=n₁+n₂+n₃. Para esses dados, que apresentam uma distribuição multinomial com parâmetros n, p_A e D_A, a função de verossimilhança é dada por . De acordo com (1), tem-se: .

A distribuição conjunta a posteriori encontrada deve ser integrada em relação a todos os outros parâmetros que a constituem, obtendo-se, assim, a distribuição marginal a posteriori de um parâmetro θ (PAULINO et al, 2003). Esta integração geralmente não é analítica, necessitando de algoritmos iterativos especializados denominados algoritmos MCMC (Markov Chain Monte Carlo), dentre os quais, destaca-se o algoritmo de Metropolis-Hastings. Este algoritmo gera um valor de uma distribuição auxiliar ou candidata e este valor é aceito com uma dada probabilidade (METROPOLIS et al., 1953; HASTINGS, 1970), caso contrário é rejeitado e um novo valor é amostrado. Para avaliar a convergência da cadeia de valores gerada pelo Metropolis-Hastings, são utilizados os critérios de GEWEKE (1992) e GELMAN & RUBIN (1992).

Muitas distribuições a priori são associadas aos parâmetros ligados ao desequilíbrio de Hardy-Weinberg, e a sua adequabilidade é avaliada por uma das grandes áreas da inferência bayesiana, a análise de sensitividade ou robustez, a qual se caracteriza pela comparação de distribuições a priori através de avaliadores de qualidade, como o fator de Bayes (KASS & RAFTERY, 1995). Este é definido por: em que representam as verossimilhanças marginais de cada modelo. Em algumas situações, as quantidades podem ser calculadas analiticamente, mas, em geral, os métodos MCMC são usados para obter soluções aproximadas. Uma interpretação para o fator de Bayes (FB) é dada em JEFFREYS (1961), em que valores de demonstram evidência a favor de M_j, demonstram evidência muito fraca a favor de M_i, demonstram evidência fraca a favor de M_i, demonstram evidência forte a favor de M_i e demonstram evidência muito forte a favor de M_i.

Um estudo de simulação foi realizado no intuito de avaliar o método empregado por cada uma das distribuições a priori implementadas. Assim, a partir de (1), vários cenários foram abordados e estes diferiram pelo tamanho da amostra (n=50; 200; 1000) e pela intensidade do parâmetro analisado, em que se considerou um valor próximo ao limite inferior do parâmetro (-0,02), um valor positivo próximo do EHW (0,02) e outro com alta endogamia (0,146) (ARMBORST, 2005), totalizando 9 cenários. Foram simuladas m=100 amostras para cada distribuição a posteriori, em que se estimaram os valores pontuais e por intervalo para cada um dos 9 cenários propostos. Utilizaram-se também os dados FBI e Cellmark retirados do trabalho de SHOEMAKER et al. (1998), que se referem às proporções genotípicas de três grupos raciais de imigrantes dos Estados Unidos (afro-americanos, caucasianos e hispânicos) localizados em três locos diferentes (D7S8, LDLR e GYPA).

A implementação dos algoritmos computacionais foi realizada no software livre R (R DEVELOPMENT CORE TEAM, 2009). Vale ressaltar que: (a) utilizou-se como função candidata a distribuição uniforme no intervalo entre o limite inferior e superior de cada parâmetro; (b) em relação aos hiperparâmetros das distribuições beta e Dirichlet, foi utilizado o valor 2, pois, neste caso, as distribuições cobriam todo o espaço paramétrico das proporções alélicas e genotípicas, respectivamente; (c) para se alcançar uma taxa de aceitação (número de vezes em que o parâmetro foi aceito ao longo das iterações) entre 20 a 50% (GILKS et al., 1996) foram testados valores adequados para o erro da proporção alélica e do coeficiente de desequilíbrio até atingirem valores satisfatórios para a convergência; (d) o procedimento proposto por NOGUEIRA et al. (2004) sugeriu 50000 iterações, em que se descartaram as 10000 iniciais (burn-in) e considerou-se um espaçamento de tamanho 40 (thin), obtendo uma amostra de tamanho 1000; (e) o critério de Geweke apresentou o valor p sempre maior que o nível de significância pré-fixado (5%), e o critério de Gelman & Rubin estimou valores de próximos a 1.

RESULTADOS E DISCUSSÃO

A partir da função de verossimilhança e das distribuições a priori adotadas, encontram-se as distribuições conjuntas a posteriori para os modelos 1, 2, 3 e 4 dadas, respectivamente, por:

As distribuições condicionais completas a posteriori para p_A, e D_A, , apresentam a mesma expressão e correspondem à distribuição conjunta a posteriori dada em (2) para a distribuição a priori Dirichlet. Para as distribuições a priori beta - função degrau uniforme e uniforme - função degrau uniforme, as distribuições condicionais completas a posteriori para p_A são dadas pelas distribuições conjuntas a posteriori (3) e (4), respectivamente. Para o caso do parâmetro D_A, estas apresentam a mesma forma e são dadas por:

As distribuições condicionais completas a posteriori para as distribuições a priori uniformes independentes são: para D_A, é

Os resultados referentes aos valores estimados (média, mediana, moda, desvio padrão e intervalo de credibilidade) para o parâmetro D_A, via simulação de dados, são apresentados na tabela 1. Observa-se que, exceto em alguns cenários para D_A =0,146, o intervalo de credibilidade contém o valor real adotado para o parâmetro. Nota-se também que os modelos apresentaram estimativas (média, mediana e moda) de D_Abastante próximas entre si, principalmente para os cenários de n=200 e n=1000, evidenciando distribuições do tipo simétricas. Vale ressaltar que estas estimativas estão bastante próximas do valor real, demonstrando que o método utilizado pode ser adaptado às diferentes situações encontradas.

Os resultados de comparação dos modelos, por meio do fator de Bayes, para os dados simulados e reais, são apresentados na tabela 2. Para os dados simulados, o valor real esteve presente em pelo menos 95% das vezes nos intervalos de credibilidade, demonstrando que o processo de simulação foi eficiente e propiciando uma melhor análise dos modelos, pois considerou vários cenários possíveis. Nestas análises observou-se que o modelo 4 é o mais indicado nos casos de D_A =0,146, D_A =0,02 e D_A =-0,02 com n=200, obtendo muitas vezes evidências fraca e muito fraca a seu favor, já o modelo 2 é o mais indicado para o caso de D_A =-0,02 com n=50, e o modelo 3 é o mais indicado considerando o parâmetro D_A =-0,02 com n=1000. Nota-se uma grande diferenciação na escolha de modelos para o parâmetro D_A, em que apenas o modelo 1 não é indicado, o que pode ser explicado, principalmente, pela distribuição a priori conjugada utilizada. Para o caso de D_A =-0,02, as diferentes escolhas dos modelos podem ser justificadas pela proximidade desse valor ao equilíbrio de Hardy-Weinberg, ou seja, D_A = 0.

Para os dados reais, observou-se que os modelos 2, 3 e 4, na maioria das vezes, apresentaram evidências muito forte e forte a seu favor quando comparados com o modelo 1 e, portanto, em nenhum dos casos, este modelo pode ser indicado, o que pode ser explicado pela utilização de uma distribuição a priori conjugada. Quando comparados entre si, os resultados foram os mais diversos possíveis, especificando um modelo para cada caso analisado, o que pode ser observado na tabela 2. Resultados semelhantes foram encontrados por SHOEMAKER et al. (1998), quando comparadas as prioris Dirichlet, beta - função degrau uniforme e uniforme - função degrau uniforme.

CONCLUSÃO

O modelo 1 (Dirichlet) é o único não indicado tanto para o processo de simulação como para cada grupo dos dados reais. A metodologia bayesiana mostrou-se eficiente no estudo do modelo de Hardy-Weinberg, sendo avaliada e comprovada pelo estudo de simulação, apresentando estimativas bem próximas ao valor real.

Recebido para publicação 26.05.10

Aprovado em 19.02.11

Devolvido pelo autor 29.04.11

CR-3632

• ARMBORST, T. Métodos para medir o desequilíbrio de Hardy-Weinberg através de medidas de endocruzamento 2005. 187f. Dissertação (Mestrado em Estatística) - Curso de Pós-graduação em Estatística, Universidade Federal de Minas Gerais, Belo Horizonte, MG.
• AYRES, K.L.; BALDING, D.J. Measuring departures from Hardy-Weinberg: a Markov chain Monte Carlo method for estimating the inbreeding coefficient. Heredity, v.80, p.769-777, 1998. Disponível em: <http://www.nature.com/hdy/journal/v80/n6/full/6883600a.html>. Acesso em: 14 set. 2009. doi:10.1046/j.1365-2540.1998.00360.x.
• COCKERHAM, C.C. Variance of gene frequencies. Evolution, v.23, p.72-84, 1969. Disponível em: <http://www.jstor.org/stable/2406485>. Acesso em: 12 jan. 2009.
• GELMAN, A.; RUBIN, D.B. Inference from iterative simulation using multiple sequences. Statistical Science, v.7, p.457-511, 1992. Disponível em: <http://projecteuclid.org/DPubS?service=UI&version=1.0&verb=Display&handle=euclid.ss/1177011136>. Acesso em: 12 out. 2009. doi: 10.1214/ss/1177011136.
• GELMAN, A. et al. Bayesian data analysis USA: Chapman & Hall/CRC, 2000. 526p.
• GEWEKE, J. Evaluating the accuracy of sampling-based approaches to calculating posterior moments. In: BERNARDO, J.M. et al. Bayesian statistics New York: Oxford University, 1992. p.625-631.
• GILKS, W.R. et al. Markov chain Monte Carlo in practice London: Chapman & Hall, 1996. 481p.
• HASTINGS, W.K. Monte Carlo sampling methods using Markov chains and their applications. Biometrika, v.57, p.97-109, 1970. Disponível em: <http://biomet.oxfordjournals.org/content/57/1/97.abstract>. Acesso em: 12 ago. 2009. doi: 10.1093/biomet/57.1.97.
• HERNÁNDEZ, J.L.; WEIR, B.S. A disequilibrium approach to Hardy-Weinberg testing. Biometrics, v.45, p.53-70, 1989. Disponível em: <http://www.jstor.org/stable/2532034>. Acesso em: 23 out. 2009.
• JEFFREYS, H. Theory of probability UK: Claredon, 1961. 325p.
• KASS, R.E.; RAFTERY, A.E. Bayes factors and model uncertainty. Journal of the American Statistical Association, v.90, p.773-795, 1995. Disponível em: <http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.143.835>. Acesso em: 13 out. 2009. doi: 10.1.1.143.835.
• METROPOLIS, N. et al. Equations of state calculations by fast computing machines. Journal of Chemical Physics, v.21, p.1087-1092, 1953. Disponível em: <http://m.jcp.aip.org/jcpsa6/v21/i6/p1087_s1>. Acesso em: 11 jun. 2009. doi: 10.1063/1.1699114.
• MUNIZ, J.A. et al. Métodos de estimação do coeficiente de endogamia em uma população diplóide com alelos múltiplos. Ciência e Agrotecnologia, v.32, p.93-102, 2008. Disponível em: <http://www.scielo.br/scielo.php?script=sci_arttext&pid=S1413-70542008000100014>. Acesso em: 23 jul. 2009. doi: 10.1590/S1413-70542008000100014.
• MUNIZ, J.A. et al. Comparação entre métodos de estimação do coeficiente de endogamia com dados de frequências alélicas em uma população diplóide. Ciência e Agrotecnologia, v.34, p.43-54, 2010. Disponível em: <http://www.scielo.br/scielo.php?pid=S1413-70542010000100005&script=sci_arttext >. Acesso em: 15 out. 2009. doi: 10.1590/S1413-70542010000100005.
• NOGUEIRA, D.A. et al. Avaliação de critérios de convergência para o método de Monte Carlo via cadeias de Markov. Revista Brasileira de Estatística, v.65, p.59-88, 2004.
• PAULINO, C.D. et al. Estatística Bayesiana Lisboa: Fundação Calouste Gulbenkian, 2003. 446p.
• R Development Core Team. R: a language and environment for statistical computing Vienna. Disponível em: <http://www.R-project.org>. Acesso em: 12 mar. 2009.
• REIS, R.L. et al. Inferência bayesiana na análise genética de populações diplóides: estimação do coeficiente de endogamia e da taxa de fecundação cruzada. Ciência Rural, v.38, n.5, p.1258-1265, 2008. Disponível em: <http://www.scielo.br/scielo.php?script=sci_arttext&pid=S0103-84782008000500009&lng=es&nrm=iso&tlng=es>. Acesso em: 18 out. 2009. doi: 10.1590/S0103-84782008000500009.
• REIS, R.L. et al. Abordagem bayesiana da sensitividade de modelos para o coeficiente de endogamia. Ciência Rural, v.39, n.6, p.1752-1759, 2009. Disponível em: <http://www.scielo.br/scielo.php?pid=S0103-84782009000600018&script=sci_arttext>. Acesso em: 17 dez. 2009. doi: 10.1590/S0103-84782009000600018.
• SHOEMAKER, J.S. et al. A Bayesian characterization of Hardy-Weinberg disequilibrium. Genetics, v.149, n.4, p.2079-2088, 1998. Disponível em: <http://www.genetics.org/cgi/content/full/149/4/2079>. Acesso em: 24 jun. 2009.
• WEIR, B.S. Genetic data analysis II. Methods for discrete population genetic data Sunderland: Sinauer Associates, 1996. 445p.

1

Autor para correspondência.

Datas de Publicação

• Publicação nesta coleção
  06 Jun 2011
• Data do Fascículo
  Maio 2011

Histórico

• Recebido
  26 Maio 2010
• Aceito
  19 Fev 2011