Laboratorio 3: Estructura genética

Author

Sofía Carvajal Rojas, 2019; Mauricio Rodríguez Bardía, 2020, 2021; Randall Hidalgo Sánchez, 2022; Alejandra Serna Sánchez, 2023, 2024

Laboratorio Estructura genética: AMOVA y DAPC

Introducción

Otra manera para estimar la estructura génica en poblaciones con estructura jerárquica es mediante el Análisis Molecular de Varianza (AMOVA). Este análisis divide la varianza molecular total en sus componentes entre grupos de poblaciones o regiones; entre las subpoblaciones dentro de los grupos y dentro de las poblaciones. Este es un método que permite estimar la diferenciación de las subpoblaciones a partir de datos moleculares y realizar pruebas de hipótesis basado en esta diferenciación. La estructura jerárquica implica que las poblaciones están estructuradas en diferentes estratos que son cada vez más inclusivos. Por ejemplo, los individuos están agrupados en sub-poblaciones, y las sub-poblaciones a su vez, están agrupadas en regiones. El AMOVA, permite determinar la diferencia en frecuencias alélicas entre sub-poblaciones dentro de una región, y también la diferenciación en frecuencias alélicas entre regiones. De esta manera, el AMOVA permite determinar cuál factor es más importante estructurando la diversidad genética. Su cálculo se basa en obtener distintas estadísticas Φ (phi).

Cuadro 1: Descripción de estadísticas Φ calculadas por AMOVA y su nivel de jerarquía.

Estadística Φ	Nivel de jerarquía
Φ_CT	Estructura entre grupos o regiones dentro de la metapoblación
Φ_ST	Estructura entre subpoblaciones dentro de la metapoblación
Φ_SC	Estructura entre subpoblaciones dentro de grupos

Existen métodos para visualizar gráficamente y estudiar la relación genética entre las poblaciones en estudio, como lo es el Análisis discriminante de componentes principales, o DAPC por sus siglas en inglés.

Objetivo general

El objetivo de esta práctica es aprender a obtener un AMOVA utilizando R, así como implementar métodos de visualización gráfica de relación genética utilizando DAPC.

1. AMOVA

El AMOVA se calcula con el comando poppr.amova (x, hier= , missing= ) de la librería poppr. Se debe indicar el objeto genind, la jerarquía o los estratos (hier=) como una fórmula jerárquica, y el método que se utilizará para corregir la presencia de valores perdidos. El método recomendado es “ignore” pues no reemplaza ni remueve los datos perdidos. Con el comando missingno(x, type=) se detallan otras opciones para lidiar con los valores perdidos. Usando help("missingno") se especifican las opciones. En este caso, se usará la opción "genotype" el cual elimina los genotipos que tengan al menos un valor perdido en alguno de sus loci.

En el primer AMOVA que realizaremos, la estructura jerárquica o la manera en que los estratos están anidados será la siguiente: los sitios están anidados dentro de las regiones en las que se encuentran. De esta forma el análisis de AMOVA nos permitirá determinar si hay diferencias en frecuencias alélicas entre sitios, entre regiones o en ambas. De la misma forma, el AMOVA determinará si hay diferencias genéticas entre los sitios dentro de las regiones. Esto puede ocurrir, si no hay un efecto de región sobre la estructura génica. La manera de indicar anidamiento de sitios dentro de región es con la fórmula: ~region/sitio. Los resultados del AMOVA se guardarán en el objeto magen.amova para análisis posteriores. Al AMOVA hay que indicarle qué hacer con los valores perdidos. Existe la opción de eliminar los loci con valores perdidos, que en este caso eliminaría a todos. También se puede eliminar lo “genotipos” o individuos con valores perdidos, y por último, se puede escoger ponerles el genotipo promedio o más común. Nosotros vamos a quitar los individuos con valores perdidos con la opción missing = "genotype".

#Recuerden establecer el directorio de trabajo donde están sus datos
setwd("/Users/asernas/Desktop/2024/ASISTENCIAS/Genética\ de\ poblaciones\ 2024-1/LAB2-Estructura_AMOVA-DAPC/")

library(poppr) 
library(hierfstat) 
library(pegas) 
library(tidyr) 
library(mmod) 
library(data.table)
library(MASS)
library(ggplot2)

data.maiz <-read.csv("maiz.csv", header = TRUE, sep = ";")
alelos.maiz <- alleles2loci(data.maiz, ploidy = 2, rownames = 1, population = 2)
magen <- loci2genind(alelos.maiz, ploidy = 2, na.alleles = "NA")
estratos<-read.csv("MaGene_Strata.csv",sep=";",header = TRUE)
otros <- list(estratos$REGION, estratos$SITIO, estratos$COLOR)
names(otros)<-c("region", "sitio", "color") #cambiar nombres
other(magen)<- otros #asignar el objeto al objeto genind magen 
strata(magen)<-data.frame(magen@other[1:3]) #las columnas 1 a la 3 contienen los estratos de interés debe ser un data.frame
nameStrata(magen) <- ~region/sitio/color 
magen

/// GENIND OBJECT /////////

 // 224 individuals; 20 loci; 351 alleles; size: 396.3 Kb

 // Basic content
   @tab:  224 x 351 matrix of allele counts
   @loc.n.all: number of alleles per locus (range: 4-37)
   @loc.fac: locus factor for the 351 columns of @tab
   @all.names: list of allele names for each locus
   @ploidy: ploidy of each individual  (range: 2-2)
   @type:  codom
   @call: adegenet::df2genind(X = as.matrix(x[, attr(x, "locicol"), drop = FALSE]), 
    sep = "/", pop = pop, ploidy = ploidy)

 // Optional content
   @pop: population of each individual (group size range: 81-143)
   @strata: a data frame with 3 columns ( region, sitio, color )
   @other: a list containing: region  sitio  color

magen.amova <- poppr.amova(magen, hier= ~region/sitio, missing = "genotype", within = F)
magen.amova

$call
ade4::amova(samples = xtab, distances = xdist, structures = xstruct)

$results
                              Df    Sum Sq  Mean Sq
Between region                 1  25.28252 25.28252
Between samples Within region  8 131.44638 16.43080
Within samples                66 672.97504 10.19659
Total                         75 829.70395 11.06272

$componentsofcovariance
                                               Sigma          %
Variations  Between region                 0.1435103   1.275519
Variations  Between samples Within region  0.9110310   8.097238
Variations  Within samples                10.1965915  90.627244
Total variations                          11.2511328 100.000000

$statphi
                          Phi
Phi-samples-total  0.09372756
Phi-samples-region 0.08201854
Phi-region-total   0.01275519

De acuerdo con el resultado del AMOVA anterior, el valor Phi-region-total es equivalente a Φ_CT (Ver Cuadro 1), por lo que la variación genética total que se debe a diferencias entre regiones es de Φ_CT =0.0127. Es decir, aproximadamente un 1.27% de la variación genética total, se da porque hay individuos en regiones distintas. El valor Phi-samples-region es equivalente a Φ_SC (Ver Cuadro 1), la variación entre los sitios dentro de las regiones es de Φ_SC = 0.0820. Finalmente, la estructuración que vemos entre sub-poblaciones, ignorando la región, es moderada con un valor de Φ_ST =0.0937.

Comparando los resultados de las estadísticas Φ se podría afirmar que el principal factor que afecta las diferencias en diversidad genética son los sitios de muestreo, no tanto así las diferencias entre regiones. Para determinar si los estadísticos Φ difieren estadísticamente de 0, se utiliza el comando randtest(xtest) de la librería “ade4”. Este comando utiliza permutaciones de los individuos entre sitios y regiones para comprobar si los valores phi difieren estadísticamente de cero. Al comando se le indica sobre cuál objeto de clase AMOVA debe operar y el número de permutaciones que se desean realizar. Este comando puede tomar mucho tiempo para calcular, por lo que utilizaremos el argumento nrepet = 99, para limitar las permutaciones a 99. En análisis más robustos, este valor debería ser superior a 999. Como resultado se produce un objeto de tipo randtest que también puede ser graficado como un histograma de los valores simulados y la posición del valor observado.

#|cache: true
magen.amova.sig <- randtest(magen.amova, nrepet = 99)
magen.amova.sig

class: krandtest lightkrandtest 
Monte-Carlo tests
Call: randtest.amova(xtest = magen.amova, nrepet = 99)

Number of tests:   3 

Adjustment method for multiple comparisons:   none 
Permutation number:   99 
                        Test        Obs     Std.Obs   Alter Pvalue
1  Variations within samples 10.1965915 -12.8530806    less   0.01
2 Variations between samples  0.9110310   9.6609063 greater   0.01
3  Variations between region  0.1435103   0.8518311 greater   0.24

De acuerdo con la prueba de hipótesis y los resultados anteriores, no se encontró una estructuración genética a nivel de regiones (Φ_CT =0.0127, p=0.24), siendo el mayor efecto sobre la estructura la estratificación por sitio (Φ_ST =0.0937, p=0.01), teniendo efecto también entre los sitios dentro de regiones (Φ_SC =0.0820, p=0.01). Este último resultado es lógico, ya que las diferencias en frecuencias alélicas se dan entre sub-poblaciones, por lo que no importa si se considera la región (Φ_SC) o no (Φ_ST), en ambos casos, la estructuración espacial es significativa.

A continuación, realizaremos otro ejemplo con los mismos datos pero con otra estructura jerárquica: color de mazorca, anidada dentro del sitio. Para ello, ejecutaremos el siguiente comando:

magen.amova2 <- poppr.amova(magen,hier= ~sitio/color, missing = "genotype", within=F)


Found 8529 missing values.

148 genotypes contained missing values greater than 5%

Removing 148 genotypes: BE12.372, BE12.381, BE12.544, BE12.547,
BE13.353, BE13.354, BE13.356, BE13.361, BE13.363, BE13.390, BE13.552,
BE15.528, BE15.529, BE15.530, BE15.531, BE15.533, BE15.536, BE15.541,
BE16.516, BE16.517, BE16.527, BE2.251, BE2.253, BE2.254, BE2.26,
BE2.263, BE2.264, BE2.265, BE2.301, BE2.304, BE5.275, BE5.276, BE5.277,
BE5.278, BE5.279, BE5.280, BE5.320, BE5.321, BE5.495, BE7.271, BE7.272,
BE7.273, BE7.28, BE7.30, BE9.465, BE9.466, BE9.467, BE9.468, CN11.100,
CN11.305, CN11.307, CN11.453, CN11.456, CN11.457, CN11.96, CN23.158,
CN23.159, CN23.161, CN23.162, CN23.163, CN23.164, CN24.316, CN24.486,
CN24.487, CN24.488, CN25.166, CN25.168, CN25.169, CN25.290, CN25.458,
CN25.459, CN25.461, CN27.176, CN27.177, CN27.179, CN27.286, CN27.287,
CN27.288, CN27.289, CN27.330, CN27.331, CN28.182, CN28.183, CN28.434,
CN28.435, CN28.436, CN28.437, CN28.438, CN28.439, CN28.440, CN28.441,
CN28.443, CN28.445, CN28.446, CN28.563, CN29.469, CN29.470, CN29.472,
CN29.474, CN29.475, CN29.477, CN29.479, CN29.480, CN29.481, CN29.482,
CN35.217, CN35.295, CN35.296, CN35.299, CN35.300, CN35.350, CN35.368,
CN35.52, CN39.329, CN39.344, CN39.345, CN39.348, CN39.448, CN39.449,
CN39.450, CN39.451, CN4.309, CN4.358, CN4.359, CN4.360, CN4.417,
CN8.322, CN8.323, CN8.325, CN8.406, CN8.407, CN8.408, CN8.409, CN8.91,
CN8.92, CN8.93, CN8.94, CC1.31, CC1.32, CC1.33, CC1.34, CC1.351, CC1.4,
CC1.490, CN44.463, CN44.464, CN44.506, CN44.510

magen.amova2

$call
ade4::amova(samples = xtab, distances = xdist, structures = xstruct)

$results
                             Df   Sum Sq   Mean Sq
Between sitio                 9 156.7289 17.414323
Between samples Within sitio  9 121.7471 13.527457
Within samples               57 551.2279  9.670665
Total                        75 829.7039 11.062719

$componentsofcovariance
                                              Sigma         %
Variations  Between sitio                 0.3264098   2.91727
Variations  Between samples Within sitio  1.1918047  10.65169
Variations  Within samples                9.6706654  86.43104
Total variations                         11.1888799 100.00000

$statphi
                        Phi
Phi-samples-total 0.1356896
Phi-samples-sitio 0.1097176
Phi-sitio-total   0.0291727

Ahora verificaremos la significancia de cada estadístico, con el siguiente comando:

#|cache: true
magen.amova.sig2 <- randtest(magen.amova2, nrepet = 99) 
magen.amova.sig2

class: krandtest lightkrandtest 
Monte-Carlo tests
Call: randtest.amova(xtest = magen.amova2, nrepet = 99)

Number of tests:   3 

Adjustment method for multiple comparisons:   none 
Permutation number:   99 
                        Test       Obs     Std.Obs   Alter Pvalue
1  Variations within samples 9.6706654 -11.7261533    less   0.01
2 Variations between samples 1.1918047   5.4939683 greater   0.01
3   Variations between sitio 0.3264098   0.6183935 greater   0.30

De acuerdo a los dos comandos anteriores, podemos afirmar que el mayor efecto de estructuración para los datos de maíz se obtiene a nivel de diferencias entre los distintos colores de mazorca (Φ_ST =0.136, p=0.01), obteniéndose efecto también entre colores dentro de un mismo sitio, pero con menor magnitud en la estructura (Φ_SC =0.110, p=0.01). Algo curioso es que la significancia a nivel de sitio se pierde al realizar el análisis incluyendo el color (Φ_CT =0.029, p=0.30), por lo que podemos afirmar que la mayor estructuración es causada por el color de la mazorca. Esto se debe, probablemente, a que en que cada sitio usualmente sólo se planta una accesión, es decir sólo se planta un único color de mazorca. De esta manera, las diferencias entre sitios, realmente se deben a diferencias entre colores. Esto es un claro ejemplo de cómo el diseño experimental y la manera en que se colectan los datos, puede tener un impacto importante sobre los datos y las conclusiones que de ellos se deriven.

2. DAPC

Este análisis nos permite organizar los individuos en un espacio de componentes principales bajo grupos definidos a priori. En este caso, los grupos definidos a priori serán aquellos definidos por la población de muestreo o los diferentes estratos definidos por el investigador. El análisis de DAPC (Discriminant analysis of principal componentes) particiona la varianza en componentes entre grupos y dentro de grupos, con el fin de maximizar la discriminación entre grupos. Los datos son primeramente transformados en sus componentes principales y luego los clusters formados son identificados usando un análisis de discriminantes. Para realizar el análisis se usa el comando dapc(x), de la librería “adegenet”.

Lo primero que se debe hacer es definir los componentes principales. Normalmente se selecciona el menor número de componentes principales que describa correctamente la variación genética. Si no se tiene mucha información, se escogen N/3 componentes principales, donde N es el número de datos y se coloca en el comando como argumento n.pca=N/3. Además, se guardan los 3 primeros discriminantes, con el argumento n.da=3. Cuando se tiene un set de datos desconocido, es una buena práctica visualizar las gráficas de PCA y DA, para luego asignar la cantidad de componentes principales y discriminantes. Esto se logra con solo quitar el argumento n.pca y n.da de la función. Hecho esto, el programa solicita ingresar los valores de cada uno en la consola una vez vista la gráfica.

Primero definimos que el análisis lo realizaremos a nivel de sitios. Seguido realizamos el DAPC con el comando dapc, le indicamos el número de PCA’s y el número de funciones discriminantes. El comando scatter realiza la figura de dispersión. Las opciones son las siguientes clabel indica que no se adicione cuadros sobre el centroide de cada grupo que indique el nombre del cluster. La opción legend = T adiciona una leyenda. La opción cleg = 0.5 indica el tamaño de la leyenda, col indica los colores para cada grupo. La opción pch indica la forma de los puntos, por lo que asigné una lista del 1 al 10 ya que tenemos 10 grupos. La opción posi.leg="bottomleft" indica la posición de la leyenda. Existen otras opciones para cambiar el diseño del gráfico, pero no las veremos en este curso:

setPop(magen) <- ~sitio #haremos la determinación a nivel de sitios 
dpcmag <- dapc(magen, n.pca = 74, n.da = 3)
miscolores<- RColorBrewer::brewer.pal(12, "Paired")
scatter(dpcmag, legend = T, cleg = 0.5, col = miscolores, pch = c(1:10), cellipse = 0, cex = 0.8, solid = 1, clabel = 0, posi.leg = "bottomleft")

Según lo anterior, solo la población Concepción y Boruca se pueden diferenciar genéticamente entre sí, mientras que el resto forman un conglomerado, siendo similares genéticamente.

Para más información de DAPC, así como de otras funciones, visitar:

http://adegenet.r-forge.r-project.org/files/tutorial-dapc.pdf