miércoles, 8 de junio de 2011

Mas allá de Hardy-Weinberg: Coeficiente de desequilibrio, una perspectiva Bayesiana en R!

Recientemente Helena Brentani y Colaboradores de la Universidad de Sao Paulo y de Brasilia propusieron un nuevo coeficiente de medida de desviación del equilibrio de Hardy-Weinberg, basado en métodos Bayesianos y modificado del coeficiente propuesto por Pereira & Rogatko en 1984, en el articulo de Brentani et al. 2011, ellos introducen dos rutinas para calcular el coeficiente de desequilibrio planteado en el articulo y el coeficiente de asociación de Yule en 1912.
En el presente post, explicare con tanto detalle como pueda la forma de seguir la rutina al hacer el calculo del coeficiente de desequilibrio, la rutina utilizada puede ser revisada en la pagina de uno de los autores (Carlos Alberto Pereira).











>set.seed(1234) 


>m<-10000000


# Para empezar, tenemos que generar un numero de permutaciones o el tamaño (m) de la muestra de Monte Carlo, que sea reproducible, y utilizamos "set.seed", para rotular (1234), la operación de generación de números al azar, de este modo, al volver a generar los números al azar y anteponer "set.seed (1234), se generaran los mismos, números (mas información).



>a1<-0.5    
>a2<-0.5    
>a3<-0.5    



# Escojo los parámetros a priori (los cuales en el articulo se sugiere que sean escogidos según la regla de Jeffreys de 1931; (a1, a2, a3=1/2).



>n1<-4    
>n2<-18    
>n3<-94

# Incluyo los valores de las frecuencias absolutas de los genotipos del gen observado; n1=AA (homocigoto dominante), n2=Aa (Heterocigoto) y 
n3=aa (Homocigoto recesivo).

>p1<-rgamma(m,a1+n1,2)
>p2<-rgamma(m,a2+n2,2)
>p3<-rgamma(m,a3+n3,2)

# Genero mis valores a posteriori, dados mis datos (n) y los parámetros
 (a) que haya incluido, utilizando el numero de permutaciones (m) que 
había incluido, utilizando la distribución gamma, que es una función de distribución de probabilidades, con una función de densidad determinada (mas información).
>theta<-4*p1*p3/(p2^2)
# Calculo el valor de theta de acuerdo a la siguiente formula, pero aplicada para valores de genotipos en donde pi12 y pi21 corresponden a los individuos heterocigotos, por eso son cambiados por pi2="p2": 
>lambda<- (theta^0.5 - 1)/(theta^0.5 + 1)
# Finalmente calculamos el coeficiente de desequilibrio (lambda) a partir del valor de theta que obtuvimos (para mas infomracion, ver el articulo), segun la formula:
# La unión de estas dos formulas y su adaptación a la genética de poblaciones resultaría en la formula general, planteada por los autores, como formula general de calculo de coeficiente de desequilibrio:
>lambda<-c(-1,lambda,1)
# Concateno, todos los valores de lambda que obtuve en todas las permutaciones que genere al azar y las pongo entre -1 y 1 (los limites de correlación del coeficiente).
>freq<-hist(lambda,breaks=200,main="Posterior of lambda",probability=T,ylab="Density")

# Genero un histograma de frecuencias, para realizar una interpretación de mis resultados mas gráfica, en el que se divide el histograma en 200 columnas y se rotulan, los ejes y el gráfico.
# En el diagrama podemos observar que este se encuentra por fuera del equilibrio HW, ya que la densidad máxima de la distribución posterior se haya muy alejada del valor de Lambda=0=Equilibrio HW.
>count_pos<-sum(freq$count[1:100])
>position_pos<-200
>while (position_pos>=101){
   if (freq$count[position_pos]<=freq$count[100])count_pos<-count_pos+freq$count[position_pos]
   position_pos<-position_pos-1
}

# Defino cuando mi funcion tiene un comportamiento de asociación positivo (Sensu Brentani et al, 2011), utilizando mis datos de frecuencias de histogramas!!!

>count_neg<-sum(freq$count[101:200])
>position_neg<-1
>while (position_neg<=100){
   if (freq$count[position_neg]<=freq$count[101])count_neg<-count_neg+freq$count[position_neg]
   position_neg<-position_neg+1
}

# Y realizo lo mismo para definir sì es negativa la asociación o la desviación de HW, mas exactamente!!

>evalue_pos<-(count_pos-2)/m
>evalue_neg<-(count_neg-2)/m
# Defino mis valores de probabilidad "e" (que es el área bajo la curva del histograma de densidad posterior que se encuentra (o se pasa) mas allá del valor de Lambda=0), se utiliza "evalue_pos" sí esta el histograma hacia los valores positivos de lambda y "evalue_neg" sí esta el histograma en los valores negativos de lambda.
>alpha<-0.05
# (1-alpha) es la credibilidad, en este caso del 95%.
>n.HPD<-m*alpha 
>LI<-1
>LS<-200
count.HPD<-freq$count[LI]+freq$count[LS]
while (count.HPD < n.HPD){
   if (freq$count[LI]<=freq$count[LS])
      {LI<-LI+1
       count.HPD<-count.HPD+freq$count[LI]
      }
   if (freq$count[LI]> freq$count[LS])
      {LS<-LS-1
       count.HPD<-count.HPD+freq$count[LS]
      }
}
# Para definir el intervalo fidedigno o de "confianza", dados los valores de frecuencia.
>mean(lambda) 
# Via promedio, obtengo el valor del coeficiente de desequilibrio, estimado por el metodod Bayesiano.
>(LI-100)/100  
>(LS-100)/100

# Obtenemos los limites fidedignos mas bajo (LI-100)/100  y  mas alto (LS-100)/100.
>evalue_pos  
>evalue_neg  
# Evalúa la probabilidad "evalue_pos" sí lambda es positivo, "evalue_neg" sì lambda es negativo.
Eso es todo, y espero les haya gustado este post, ya que es bastante actual y tiene varias cosas interesantes!!!

2 comentarios: