---
title: "Equilibrio Hardy-Weinberg"
author: "Enrique Lessa"
date: "09/04/2019"
output:
  word_document: default
  pdf_document: default
  html_document: default
---

**Equilibrio Hardy-Weinberg**

**Introducción**

Vamos a usar R, un paquete de programas de uso libre para estadística y un gran conjunto de aplicaciones relacionadas (análisis de datos, gráficas, presentación de resultados).  Como interfase usamos RStudio, ampliamente utilizado, y para producir los archivos de salida usamos Rmarkdown.  En Rmarkdown, intercalamos texto común (incluyendo ecuaciones en Latex) con bloques ("chunks") de código (las líneas de código en R propiamente dichas).  Si producimos una salida (en nuestro caso en html), esta intercalará texto, código y resultados del código.

El primer bloque de código activa "knitr" y "markdown", necesarios para las salidas en html.

```{r global_options, include=FALSE}
knitr::opts_chunk$set(warning=FALSE, message=FALSE) #solo para que no salten mensajes de que tal o cual paquete está desarrollado con otra version de R y queden en el htm
```


```{r setup}
library("knitr")
knitr::opts_chunk$set(echo = TRUE)
library("rmarkdown")
library("lattice")
#no se si era necesario pero activé readbitmap y saqué las comillas a RColorbrewer
library("readbitmap")
# algunos detalles (algo oscuros) para usar en lattice
library("RColorBrewer")
 MisColores = brewer.pal(6, "Accent")
  my.settings = list(col = MisColores[], superpose.polygon=list(col=MisColores[1:3]), strip.background=list(col=MisColores[6]))
```

**Distribución binomial y frecuencias alélicas**

Ya realizamos una exploración de la distribución binomial, de la que solamente repetimos la introducción.

Si conocemos la frecuencia real del alelo A1 *p=f(A1)* en la población, podemos aplicar la binomial para calcular la probabilidad de observar *i* alelos de tipo A en una muestra de tamaño *n*. Como es lógico, dicha probabilidad depende de la frecuencia del alelo y del tamaño de la muestra.  En concreto, toma la siguiente forma:

$P(i) = {n \choose i}p^{i}(1-p)^{n-i}$,                        *i = 0, 1, ..., n*

Notamos que:

$p^{i}$ es la probabilidad de muestrear *n* veces el alelo A1.
$(1-p)^{n-i}$ es la probabilidad de muestrear *n-i* veces el alelo A1.
${n\choose i}$ es el número de formas de obtener los resultados anteriores.

Este es un ejercicio deductivo: dada una frecuencia alélica conocida (*p*) y un tamaño de muestra, deducimos con qué probabilidad podemos obtener todos los resultados posibles, desde *i=0* hasta *i=n*.

**Aplicación al modelo de Hardy-Weinberg**

Por ejemplo, tomamos una muestra de *n=50* alelos, y queremos saber cual es la probabilidad de que 20 de ellos sean de tipo A1 (*i=20*). 

Recordemos que, en el modelo de Hardy-Weinberg, consideramos una población de reproducción sexuada de tamaño infinito, sin inmigración, sin mutación, con frecuencias alélicas idénticas en los dos sexos, en la que los genotipos de una generación son combinaciones al azar de los alelos de la generación precedente. La oferta de alelos está dada por los gametos, en los que las frecuencias alélicas son idénticas en los dos tipos de gametos.  

Para el caso del modelo original (con dos tipos de alelos), podemos aplicar la distribución binomial para formar pares de alelos, tomados al azar en base a la frecuencia de las clases alélicas en la población.  En cada par, puede haber 0, 1 o 2 copias del alelo de referencia A1, lo que corresponde a los genotipos A2A2, A1A2, y A1A1, respectivamente.

```{r Parámetros}
# Definimos aquí los parámetros para las secciones siguientes:
x1 = 50 # número de genotipos por muestra
n1 = 2  # tamaño de la muestra en cada réplica (en este caso, los dos alelos que forman un genotipo diploide)
p1 = 0.6 # frecuencia del alelo A1 en la población [Nota: el código actual, provisorio, funciona mejor lejos del 0 y del 1; ver #C más abajo]
y = 4 # número de muestras (número de veces que repetiremos el muestreo al azar en la sección #F)

```



```{r Hardy-Weinberg usando la binomial}

#A. Frecuencias esperadas relativas y absolutas:
  Esperadas_relativas = c((1-p1)**2, 2*p1*(1-p1), p1**2) # frecuencias relativas esperadas de A2A2, A1A2, y A1A1
  Esperadas = Esperadas_relativas * x1 # frecuencias absolutas de los tres genotipos 
  print(Esperadas)

#B. Una muestra al azar de genotipos obtenida con la binomial (función *rbinom*)
  Muestra_1 = rbinom(x1,n1,p1)
  print(Muestra_1)
  Observadas_1 = table(Muestra_1) # tabla de frecuencias en la muestra
  print(Observadas_1)

#C. Combinamos las frecuencias observadas y esperadas en una misma tabla
  # [Nota: falla si, por azar, el muestreo no incluye los tres genotipos]
  Esp = c(rep(0, Esperadas[1]), rep(1, Esperadas[2]), rep(2, Esperadas[3]))
  par(mfrow=c(1,2))
  hist(Esp, breaks = 10, col=rgb(1,1,0,0.7), main = "E: amarillo; O: azul", xlab = "0:A2A2, 1:A1A2, 2:A1A1")
  par(new = TRUE)
  hist(Muestra_1, breaks = 10, col=rgb(0,1,1,0.4), main ="", xlab = "", axes=FALSE)
 
# Frec. esperadas en amarillo y observadas en azul, con la superposición en verde por transparencia)  
```

**Usando la función "sample"**

En el siguiente bloque repetimos el muestreo del anterior con un código ligeramente distinto:  
  - definimos los objetos (genotipos) a ser muestreados;
  - definimos las probabilidades correspondientes usando las frecuencias "Esperadas_relativas" definidas;
  - en el bloque anterior, seguimos el orden natural de los genotipos (0, 1, y 2 copias de A); ahora usamos el orden más usual en genética, A1A1, A1A2, A2A2, usando esos genotipos como objetos de muestra en la función *sample*;
   - (notar la opción "replace = TRUE", que corresponde al muestreo con reposición [indica que las probabilidades se mantienen constantes: no cambian a medida que vamos muestreando los distintos genotipos].
  

```{r Hardy-Weinberg_función_sample, fig.width=5, fig.height=5}

#D. Una muestra al azar, usando la función *sample* de manera equivalente a *rbinom*.
  Genotipos = c("A1A1", "A1A2", "A2A2") # creamos la lista de genotipos de interés
  Esperadas_relativas_2 = c(p1**2, 2*p1*(1-p1), (1-p1)**2)
  Muestra_2 = sample(Genotipos,  x1, replace = TRUE, Esperadas_relativas_2) # obtenemos una muestra de tamaño x1, muestreando al azar con reposición los genotipos, cada uno con una probabilidad definida más arriba (Esperadas_relativas_2)
  Observadas_2 = table(Muestra_2) # tabla de frecuencias en la muestra
  print(Observadas_2)
  
  Esperadas_2 = Esperadas_relativas_2 * x1 
  
#E. Combinamos las frecuencias observadas y esperadas en una misma tabla
  Resumen_2 = rbind(Esperadas_2, Observadas_2)
  print(Resumen_2)
 # probemos crear un resumen sin depender del formato generado en table
  Resumen_3 = matrix(data = NA, nrow = y+2 , ncol= 3)
  colnames(Resumen_3) = c("A1A1", "A1A2", "A2A2")
  rownames(Resumen_3) = c("Esp.", letters[1:(y+1)]) # la primera fila toma el nombre "Esp", y las siguientes se etiquetan 
                                                    # con letras hasta el valor definido en "y", usado más abajo para 
                                                    # replicar los muestreos
  Resumen_3[1,] = Esperadas_2
  Resumen_3[2,] = Observadas_2
  
#F. COMBINANDO MÚLTIPLES MUESTRAS ADICIONALES ()
  Resumen_y = (Resumen_3) # comenzamos incorporando las observadas y esperadas generadas más arriba
  for (i in seq(1:y)){
    Muestra_y = sample(Genotipos,  x1, replace = TRUE, Esperadas_relativas_2) # obtenemos una muestra de tamaño x1, muestreando al azar con reposición los genotipos, cada uno con una probabilidad definida más arriba (Esperadas_relativas)
    Observadas_y = table(Muestra_y) 
     # if((length(Observadas_y) =3)
        Resumen_y[i+2,] = Observadas_y
    }
    print(Resumen_y)
    
# Pruebas con gráficas básicas de "lattice"
   Fig_1 = barchart(Resumen_y[1:(y+2),], 
         xlab = "Frecuencias absolutas",  
         main = "Frecuencias genotípicas",  
         auto.key= list(space = "top", columns = 1, points = FALSE, 
         rectangles = TRUE, title = "Genotipos", cex.title = 1),
         panel = lattice.getOption("panel.barchart"),
         default.prepanel = lattice.getOption("prepanel.default.barchart"),
         box.ratio = 2,
         par.settings = my.settings    )
  plot(Fig_1)

```


**Incorporando la endocría**

Muchas poblaciones reales se apartan de la panmixia, de modo que los apareamientos ocurren con una mayor probabilidad entre individuos emparentados.  Aunque con frecuencia la reproducción entre individuos fuertemente emparentados se evitan, efectos de vecindario, organización social, y otros tienen a hacer que los apareamientos entre individuos con un parentesco mayor que el promedio de la población sean más frecuentes de lo esperado por azar.

Para incorporar de manera sencilla este fenómeno general (sin obligarnos por ello a estudiar pedigrís de manera directa), introducimos *F*, el coeficiente de endocría (o endogamia) de la población.  Se trata de un único parámetro adicional que procura capturar el efecto neto de la endocría sobre las frecuencias genotípicas esperadas.

Si los alelos se aparean con sus símiles (A1 con A1, y A2 con A2) con probabilidad *F*, entonces aumentarán las frecuencias de homocigotas A1A1 y A2A2 a expensas de las frecuencias de heterocigotas (A1A2).  

Planteamos ahora lo razonado más arriba de manera explícita:

(1)	Una fracción *1-F* de las combinaciones gaméticas se realizan al azar, produciendo los siguientes resultados parciales:
  $frec.(A1A1) = p^2 (1-F)$
  $frec.(A1A2) = 2pq (1-F)$			
  $frec.(A2A2) = q^2 (1-F)$
  
(2)	Por otra parte, la restante fracción *F* de las combinaciones gaméticas corresponden a la endocría, es decir que no son al azar, sino que combinan gaméticos idénticos.  Por tanto, estas combinaciones estarán en proporción a las frecuencias de los alelos, y serán *pF* y *qF* para los genotipos AA y aa, respectivamente.  En combinación con el resultado anterior, obtenemos:

  $frec.(A1A1) = p^2 (1-F) + pF$
  $frec.(A1A2) = 2pq (1-F)$
  $frec.(A2A2) = q^2 (1-F) + qF$
  
Sumando las frecuencias, verificamos el resultado:

$p^2 (1-F) + 2pq (1-F) + q^2 (1-F) + pF + qF = (1-F) (p^2 + 2pq + q^2) + F (p + q) = 1$

Notamos, de paso, que la expresión de frecuencias esperadas obtenida en (2) es también válida cuando F = 0, en cuyo caso se simplifica a la expresión en (1).


