practica_1
LC
2023-05-24
En este documento encontrará dos prácticas del sitio web https://popgen.nescent.org/index.html para el análisis de datos SNPs. Las dos prácticas son
- Estadística de base
- Diferenciación de poblaciones.
#PRÁCTICA nº 1 Cálculo de estadísticas genéticas poblacionales básicas a partir de datos SNP
Introducción En esta viñeta, calculará estadísticas genéticas poblacionales básicas a partir de datos SNP utilizando paquetes R. Estas estadísticas sirven como análisis exploratorio y requieren trabajar a nivel poblacional. Calcularemos:
Diversidad genética, Prueba Hardy Weinberg Fis y Fst global.
El conjunto de datos utilizado para estos análisis corresponde al pino lodgepole (Pinus contorta, Pinaceae), una especie vegetal Para que los cálculos sean más rápidos, trabajamos sólo con un subconjunto del conjunto de datos completo.
Flujo de trabajo para datos SNP
Importación de datos Los datos se almacenan en un archivo de texto (genotype=AA..). Importaremos el conjunto de datos a R como un marco de datos y, a continuación, convertiremos el archivo de datos SNP en un objeto “genind”.
El conjunto de datos “Master_Pinus_data_genotype.txt” puede descargarse aquí.
El archivo de texto es una matriz de (550 filas x 3086 columnas). Contiene 4 columnas adicionales: la primera columna es la etiqueta de los individuos, las otras tres son la descripción de la región, todas las demás columnas son para los genotipos como (AA o AT…).
Al importar los datos en R, el archivo de datos debe estar en su directorio de trabajo, o ajustar la ruta en la invocación read.table() a continuación en consecuencia.
library("adegenet")## Loading required package: ade4##
## /// adegenet 2.1.10 is loaded ////////////
##
## > overview: '?adegenet'
## > tutorials/doc/questions: 'adegenetWeb()'
## > bug reports/feature requests: adegenetIssues()library("hierfstat")##
## Attaching package: 'hierfstat'## The following objects are masked from 'package:adegenet':
##
## Hs, read.fstatlibrary("pegas")## Loading required package: ape##
## Attaching package: 'ape'## The following objects are masked from 'package:hierfstat':
##
## pcoa, varcomp## Registered S3 method overwritten by 'pegas':
## method from
## print.amova ade4##
## Attaching package: 'pegas'## The following object is masked from 'package:ape':
##
## mst## The following object is masked from 'package:ade4':
##
## amova#1 Cargar los datos.
Mydata <- read.table("D:/Gene_avan/Genetica de poblaciones/Practica/Master_Pinus_data_genotype.txt", header = TRUE)
dim(Mydata)## [1] 550 3086Para trabajar con los datos, necesitamos convertir el objeto R devuelto por read.table() en un objeto “genind”. Para ello, creamos una matriz con sólo genotipos, y mantenemos sólo un subconjunto de los 11 primeros loci SNP (para hacer los cálculos más rápidos). A continuación, el resultado puede convertirse en un objeto “genind” (para el paquete adegent). A continuación, el objeto “genind” puede convertirse fácilmente en un objeto “hierfstat” (paquete hierfstat).
#Creación del objeto
locus <- Mydata[, -c(1, 2, 3, 4, 17:3086)]
colnames(locus) <- gsub("\\.", "_", colnames(locus)) # locus names can't have "."
ind <- as.character(Mydata$tree_id) # labels of the individuals
population <- as.character(Mydata$state) # labels of the populations
Mydata1 <- df2genind(locus, ploidy = 2, ind.names = ind, pop = population, sep = "")
Mydata1## /// GENIND OBJECT /////////
##
## // 550 individuals; 12 loci; 24 alleles; size: 102.8 Kb
##
## // Basic content
## @tab: 550 x 24 matrix of allele counts
## @loc.n.all: number of alleles per locus (range: 2-2)
## @loc.fac: locus factor for the 24 columns of @tab
## @all.names: list of allele names for each locus
## @ploidy: ploidy of each individual (range: 2-2)
## @type: codom
## @call: df2genind(X = locus, sep = "", ind.names = ind, pop = population,
## ploidy = 2)
##
## // Optional content
## @pop: population of each individual (group size range: 4-177)Identificación de número de alelos
nAll(Mydata1)## X0_10037_01_257 X0_10040_02_394 X0_10044_01_392 X0_10048_01_60 X0_10051_02_166
## 2 2 2 2 2
## X0_10054_01_402 X0_10067_03_111 X0_10079_02_168 X0_10112_01_169 X0_10113_01_119
## 2 2 2 2 2
## X0_10116_01_165 X0_10151_01_86
## 2 2#Creando el objeto “Hierfstat”
Mydata2 <- genind2hierfstat(Mydata1) # Create hierfstat objectDiversidad genetica esperada
Diversidad genética (heterocigosidad observada y esperada)
Con adegenet:
div <- summary(Mydata1)
div##
## // Number of individuals: 550
## // Group sizes: 10 73 177 58 169 7 10 24 9 9 4
## // Number of alleles per locus: 2 2 2 2 2 2 2 2 2 2 2 2
## // Number of alleles per group: 22 24 24 23 24 20 19 23 21 19 18
## // Percentage of missing data: 1.62 %
## // Observed heterozygosity: 0.42 0.49 0.46 0.4 0.12 0.47 0.07 0.08 0.11 0.46 0.04 0.14
## // Expected heterozygosity: 0.41 0.5 0.5 0.45 0.13 0.5 0.07 0.08 0.11 0.44 0.04 0.16names(div)## [1] "n" "n.by.pop" "loc.n.all" "pop.n.all" "NA.perc" "Hobs"
## [7] "Hexp"Gráfico de la heteocigocidad observada por locus
plot(div$Hobs, xlab="Loci number", ylab="Observed Heterozygosity",
main="Observed heterozygosity per locus")
Gráfico de la heteocigocidad esperada por locus
plot(div$Hobs,div$Hexp, xlab="Hobs", ylab="Hexp",
main="Expected heterozygosity as a function of observed heterozygosity per locus")
#Barlett test
bartlett.test(list(div$Hexp, div$Hobs)) # a test : H0: Hexp = Hobs##
## Bartlett test of homogeneity of variances
##
## data: list(div$Hexp, div$Hobs)
## Bartlett's K-squared = 0.011457, df = 1, p-value = 0.9148No se rechaza la hipotesis nula.
Observamos que la heterocigosidad varía entre loci. No observamos diferencias entre la heterocigosidad esperada y la observada.
Estadística básica con hierfstat. Las poblaciones son estados. La función basic.stats() proporciona la heterocigosidad observada (Ho), las diversidades génicas medias dentro de la población (Hs), Fis y Fst.
basicstat <- basic.stats(Mydata2, diploid = TRUE, digits = 2)
names(basicstat) ## [1] "n.ind.samp" "pop.freq" "Ho" "Hs" "Fis"
## [6] "perloc" "overall"La función boot.ppfis() proporciona el intervalo de confianza para Fis. La función indpca() realiza un PCA en la matriz centrada de las frecuencias alélicas de los individuos.
boot.ppfis(Mydata2) ## $call
## boot.ppfis(dat = Mydata2)
##
## $fis.ci
## ll hl
## 1 -0.0265 0.2753
## 2 -0.0528 0.0369
## 3 0.0156 0.1011
## 4 -0.0606 0.0650
## 5 -0.0397 0.0788
## 6 -0.4200 0.1737
## 7 -0.2240 0.1363
## 8 -0.2157 -0.0674
## 9 -0.1464 0.3218
## 10 0.0371 0.3514
## 11 -0.4348 -0.1490Revisión
x <- indpca(Mydata2)
plot(x, cex = 0.7)
#Pruebas del equilibrio de Hardy-Weinberg Utilizamos la función hw.test() del paquete pegas.
hw.test(Mydata1, B = 1000)## chi^2 df Pr(chi^2 >) Pr.exact
## X0_10037_01_257 0.75465490 1 0.385006451 0.466
## X0_10040_02_394 0.16628227 1 0.683437231 0.750
## X0_10044_01_392 2.96053195 1 0.085319867 0.090
## X0_10048_01_60 6.39076800 1 0.011471539 0.008
## X0_10051_02_166 1.12035389 1 0.289842265 0.306
## X0_10054_01_402 1.54582494 1 0.213752853 0.223
## X0_10067_03_111 0.04868128 1 0.825374035 0.566
## X0_10079_02_168 0.01579313 1 0.899992594 0.590
## X0_10112_01_169 0.41126025 1 0.521330544 1.000
## X0_10113_01_119 0.74248768 1 0.388865151 0.430
## X0_10116_01_165 1.87098306 1 0.171362531 0.238
## X0_10151_01_86 8.94232481 1 0.002786379 0.007Obtenemos para cada locus una prueba de significación de la hipótesis nula: H0 el locus está en equilibrio HW en la población. Todos los locus menos uno están en equilibrio HW.
Conclusión ¿Qué hemos aprendido hoy? En esta viñeta, hemos aprendido a explorar los patrones de diversidad genética y a estimar el nivel de diferenciación genética en una población. Además, tiene una idea de las posibles violaciones del conjunto de datos al modelo nulo de Wright-Fischer.
¿Cuál es el siguiente paso? Puede que ahora desee profundizar en la diferenciación de poblaciones (véase Cálculo de la diferenciación genética y métodos de agrupación a partir de datos SNP).