Diversidad jerárquica: desde genes hasta ecosistemas
GL, YR
31/10/2019
Partición de la diversidad y la diferenciación en sus componentes jerárquicos a través de métricas y escalas, desde genes hasta ecosistemas
Partitioning hierarchical diversity and differentiation across metrics and scales, from genes to ecosystems
Package: HierDpart
Version: 0.5.0
Date: 2019-09-30
Author: Xinghu Qin
Maintainer: Xinghu Qin qinxinghu@gmail.com
URL: https://github.com/xinghuq/HierDpart
Instalación del paquete:
install.packages("HierDpart")
Cargar paquetes:
library(HierDpart)
Ejercicio práctico 1
Tejones/Tlalcoyotes (Taxidea taxus) en Canadá
Fuente: Rico et al. (2016) Evolutionary Applications 9:1271-1284
Se genotiparon 20 loci microsatélites en 236 individuos pertenecientes a 8 poblaciones y 3 subespecies
| Subespecie | Población | n |
|---|---|---|
| T.t.jeffersonii | TO | 30 |
| T.t.taxus | EK | 30 |
| AB | 39 | |
| SK | 25 | |
| MB | 48 | |
| UP | 22 | |
| T.t.jacksoni | LP | 23 |
| ON | 19 | |
| 236 |
La estructura genética poblacional reportada por Rico et al. 
Cargar archivo de genotipos en formato ‘genepop’:
infile_Ttaxus <- "Taxidea taxus_genodata_f.gen"
infile_Ttaxus_5p <- "Taxidea taxus_genodata_5pob.gen"
Cálculo del perfil de diversidad genética (q=0,1,2)
Esta función obtiene diferentes indicadores de diversidad genética, donde q es el órden de la serie de Hill:
q=0 riqueza
q=1 exponencial de Shannon
q=2 medida relacionada con la heterocigosidad (^2D=1/(1-He)
q=0
Riqueza alélica
*Insensible a las frecuencias
qD(infile_Ttaxus,q=0,ncode=3)## TO EK AB SK MB UP LP ON
## Tt13 8 10 12 10 11 8 10 5
## Tt15 8 10 12 10 12 7 8 4
## Tt17 5 7 6 6 6 5 7 4
## Tt20 6 7 9 9 9 7 6 3
## Tt27 7 9 10 10 11 8 7 3
## Tt21 6 7 9 6 8 7 7 4
## Tt22 5 5 7 6 7 5 4 3
## Tt23 4 7 7 6 8 7 6 2
## Mp85 3 4 3 4 4 3 2 3
## Gg465 8 10 9 9 13 5 6 5
## Gg443 7 10 10 9 11 9 9 5
## Mvis72 7 5 7 7 7 5 5 3
## Tt3 5 10 10 9 11 7 7 3
## Tt4 7 8 9 10 10 9 10 3
## Mvi87 8 8 9 9 8 3 6 6
## Tt2 4 5 4 4 5 4 4 4
## Tt1 7 10 10 9 11 7 5 1
## Gu234 5 7 7 6 7 7 8 3
## Ma15 10 10 11 11 12 10 9 2
## Ma1 7 8 6 6 8 5 4 2
q=1
Diversidad alélica
*Pondera la contribución de cada alelo de acuerdo a su frecuencia, sin favorecer alelos comunes ni raros
qD(infile_Ttaxus,q=1,ncode=3)## TO EK AB SK MB UP LP
## Tt13 6.039466 8.890825 9.307952 8.085812 8.581448 5.924812 8.680116
## Tt15 4.998754 5.862826 10.007902 7.923405 9.803529 4.989340 5.802214
## Tt17 3.850213 3.956237 4.731108 4.853878 4.397755 2.776050 3.108485
## Tt20 4.366171 4.899730 5.515481 6.513457 5.877412 4.144644 4.486140
## Tt27 5.647554 3.992461 8.084900 8.375501 8.479766 5.407630 4.874079
## Tt21 4.248662 4.506320 6.111850 5.089336 5.452958 4.878211 6.398330
## Tt22 2.707854 3.066415 3.542006 3.234996 3.988592 3.595778 2.535242
## Tt23 2.966102 4.854691 5.369215 4.525666 4.652982 5.507868 3.718577
## Mp85 2.788529 2.403531 2.824242 2.475029 3.226789 2.179822 1.969830
## Gg465 5.265485 7.376335 6.175468 6.825248 7.695547 3.208035 3.757439
## Gg443 4.280599 6.915409 7.613655 7.957918 8.441518 7.159818 6.240405
## Mvis72 5.481393 4.287240 5.671412 5.690812 5.826004 4.281010 3.276945
## Tt3 2.707854 5.909792 7.752328 6.913233 8.731225 5.392406 5.293954
## Tt4 5.574882 5.053736 7.867970 7.371290 7.787561 7.564930 5.065880
## Mvi87 5.276371 5.523644 5.835918 6.591431 6.327923 2.600490 4.211295
## Tt2 3.073395 3.703451 3.697986 3.899349 3.979480 3.084634 3.091544
## Tt1 5.007498 6.876256 7.261615 6.642856 8.051406 4.483753 3.696949
## Gu234 3.192649 4.555026 4.738183 4.895127 4.903574 4.734119 4.410172
## Ma15 8.152394 7.838565 8.423818 9.535711 9.740789 6.795067 5.726580
## Ma1 4.099989 5.966991 3.294587 4.953947 4.668157 3.017497 2.380369
## ON
## Tt13 2.369992
## Tt15 2.867871
## Tt17 2.443617
## Tt20 2.085894
## Tt27 2.093695
## Tt21 2.785574
## Tt22 1.811170
## Tt23 1.931139
## Mp85 1.877668
## Gg465 3.542621
## Gg443 2.870793
## Mvis72 2.228290
## Tt3 2.644497
## Tt4 1.274652
## Mvi87 4.536723
## Tt2 3.480351
## Tt1 1.000000
## Gu234 1.485407
## Ma15 1.691002
## Ma1 1.129407
q=2
Número efectivo de alelos
*Favorece desproporcionadamente a los alelos más comunes
qD(infile_Ttaxus,q=2,ncode=3)## TO EK AB SK MB UP LP
## Tt13 5.172414 8.256881 8.266304 7.225434 7.349282 4.913706 7.934426
## Tt15 4.000000 3.956044 8.973451 6.613757 8.878613 4.321429 4.787330
## Tt17 3.522505 3.345725 4.072289 4.125413 3.668790 2.195011 2.489412
## Tt20 3.673469 3.773585 4.190083 5.364807 4.706844 3.545788 3.446254
## Tt27 4.918033 2.821317 7.312500 7.485030 7.480519 4.000000 3.947761
## Tt21 3.529412 3.370787 5.281250 4.562044 4.617234 3.983539 6.045714
## Tt22 2.107728 2.356021 2.582343 2.626050 3.107215 2.880952 2.275269
## Tt23 2.635432 4.176334 4.768025 3.858025 3.315108 4.460829 2.776903
## Mp85 2.650957 2.155689 2.661417 2.204586 3.090543 2.063966 1.941284
## Gg465 4.295943 5.980066 4.919932 5.482456 6.087186 2.659341 2.813830
## Gg443 3.180212 5.750799 6.500000 7.183908 6.727007 6.165605 5.186275
## Mvis72 4.932551 3.930131 5.102473 5.020080 5.189189 3.811024 2.664987
## Tt3 2.107728 4.422604 6.805369 6.009615 7.616529 4.792079 4.502128
## Tt4 5.013928 3.773585 6.730088 6.067961 7.035115 6.630137 3.327044
## Mvi87 4.346253 4.379562 4.436252 5.433962 5.685972 2.333333 3.431907
## Tt2 2.870813 3.364486 3.437288 3.799392 3.725141 2.838710 2.821333
## Tt1 4.390244 5.289308 6.145455 5.530973 6.939759 3.372822 3.391026
## Gu234 2.538787 3.621730 3.860406 4.370629 3.935098 3.781250 3.120944
## Ma15 6.921811 6.593407 7.365617 8.680556 8.844530 5.260870 4.426778
## Ma1 3.259690 5.096970 2.600000 4.325260 3.265769 2.554090 1.885918
## ON
## Tt13 1.748184
## Tt15 2.551237
## Tt17 2.194529
## Tt20 1.875325
## Tt27 1.773956
## Tt21 2.398671
## Tt22 1.520000
## Tt23 1.870466
## Mp85 1.593819
## Gg465 2.934959
## Gg443 2.117302
## Mvis72 1.885117
## Tt3 2.359477
## Tt4 1.112481
## Mvi87 3.658537
## Tt2 3.125541
## Tt1 1.000000
## Gu234 1.238422
## Ma15 1.519288
## Ma1 1.054015
Gráficar los perfiles de diversidad en los tres órdenes:
qDplot(infile_Ttaxus,q="all",ncode=3)
Pregunta:
¿Cuál es la interpretación biológica de la figura?
Desglose de la diversidad genética en sus diferentes niveles jerárquicos
Cálculo de la riqueza alélica jerárquica de 8 poblaciones estructuradas en 3 subespecies (regiones)
Primero necesitamos especificar la estructura jerárquica de los datos:
str_T <- as.matrix(read.table("str_Ttaxus.txt", header=FALSE, sep="\t"))
str_T_5p <- as.matrix(read.table("str_Ttaxus_5pob.txt", header=FALSE, sep="\t"))
Llamamos la función ‘HierAr’ incluyendo el arreglo jerárquico:
HAr_Ttaxus = HierAr(infile_Ttaxus, nreg=3, r=c(1,5,2), ncode=3)
donde:
nreg = número de agregados o poblaciones dentro de una metapoblación
r = número de subpoblaciones dentro de cada agregado o población
y obtenemos:
print(HAr_Ttaxus)## $Ar_pop
## Arpop 1 Arpop 2 Arpop 3 Arpop 4 Arpop 5 Arpop 6 Arpop 7
## Tt13 7.106791 9.340665 9.710135 8.957776 9.252234 7.454554 9.451770
## Tt15 6.579323 8.271192 10.301368 9.077574 9.970491 6.265598 7.181809
## Tt17 4.473905 5.254175 5.679437 5.911042 5.499904 4.489429 5.260870
## Tt20 5.417622 6.656556 7.160834 7.958140 7.123213 5.727273 5.937567
## Tt27 6.606242 6.516949 8.607606 9.213697 8.990081 7.452711 6.485002
## Tt21 5.415420 6.180746 6.999415 5.836028 6.449641 6.481630 6.871775
## Tt22 3.997553 4.641829 5.506494 4.986735 5.652730 4.876229 3.536232
## Tt23 3.748966 5.878906 6.161810 5.539326 6.385197 6.846815 5.648985
## Mp85 2.999601 3.254237 2.999962 3.444898 3.529600 2.681818 2.000000
## Gg465 6.472142 8.686694 7.587424 8.410748 8.764573 4.626850 5.521882
## Gg443 6.078266 8.120887 8.571820 8.705771 9.453971 8.304539 7.818263
## Mvis72 6.266372 4.917704 6.393890 6.675142 6.540881 4.970655 4.536181
## Tt3 3.997553 7.872696 8.470662 8.001806 9.283144 6.489337 6.498435
## Tt4 6.352947 6.719618 8.670859 8.807854 8.232743 8.674622 8.104801
## Mvi87 6.847538 7.011120 7.359484 8.165787 6.975828 2.999824 5.823104
## Tt2 3.499938 4.508474 3.996836 3.999964 4.309947 3.681801 3.652174
## Tt1 5.938512 8.493927 8.231047 8.152425 8.699075 6.656448 4.536232
## Gu234 4.438521 6.000907 6.030196 5.802821 6.316800 6.461855 6.998682
## Ma15 9.288381 8.930432 9.083795 10.178555 9.896401 8.688711 7.789350
## Ma1 5.762910 7.139460 4.630445 5.822845 6.690880 4.356263 3.641011
## Arpop 8
## Tt13 4.539108
## Tt15 3.788525
## Tt17 3.578947
## Tt20 2.789474
## Tt27 2.960171
## Tt21 3.788525
## Tt22 2.789473
## Tt23 2.000000
## Mp85 2.789474
## Gg465 4.919393
## Gg443 4.748686
## Mvis72 2.993362
## Tt3 3.000000
## Tt4 2.578947
## Mvi87 6.000000
## Tt2 3.999042
## Tt1 1.000000
## Gu234 2.782835
## Ma15 2.000000
## Ma1 1.789474
##
## $Hierareco
## [,1]
## Tt13 12.955723
## Tt15 13.953213
## Tt17 10.813559
## Tt20 11.999992
## Tt27 11.000000
## Tt21 9.000000
## Tt22 8.953390
## Tt23 8.951312
## Mp85 4.995844
## Gg465 12.955720
## Gg443 13.904701
## Mvis72 7.961538
## Tt3 11.953386
## Tt4 12.951312
## Mvi87 12.000000
## Tt2 5.953301
## Tt1 12.957381
## Gu234 9.000000
## Ma15 12.000000
## Ma1 10.960137
##
## $Ar_reg
## Ar_region 1 Ar_region 2 Ar_region 3
## 1 6.318192 10.2352 6.740778
##
## $Hierar_loc
## Ar_region 1 Ar_region 2 Ar_region 3
## Tt13 7.966102 11.975088 9.957985
## Tt15 7.932768 13.975591 7.695327
## Tt17 4.966667 7.975610 7.266781
## Tt20 5.966667 11.999478 5.997176
## Tt27 6.999435 11.000000 6.835873
## Tt21 5.966667 9.000000 6.980810
## Tt22 4.933333 8.975610 4.695353
## Tt23 3.999435 8.975088 5.943176
## Mp85 3.000000 4.998956 2.857143
## Gg465 7.900000 12.963190 6.959569
## Gg443 6.966102 11.951220 9.571144
## Mvis72 7.000000 7.000000 4.838210
## Tt3 4.933333 11.975610 6.835901
## Tt4 6.966667 11.975591 10.121590
## Mvi87 8.000000 11.000000 8.000000
## Tt2 3.966667 4.999990 4.000000
## Tt1 6.933333 12.981592 4.838210
## Gu234 4.966667 9.000000 8.671546
## Ma15 10.000000 11.999990 8.768957
## Ma1 7.000000 9.981312 3.980811
##
## $Ar_ovell
## Art Arr Arp
## [1,] 10.76103 7.764722 6.123887donde:
Ar_pop
Riqueza alélica por locus para cada una de las poblaciones
Hierareco
Riqueza alélica total por locus
Ar_reg
Riqueza alélica para cada agregado (e.g. región)
Hierar_loc
Riqueza alélica por locus p ara cada agregado (e.g. región)
Ar_overall
Riqueza alélica (efectiva) en los diferentes niveles jerárquicos
Ahora graficamos los resultados para poblaciones y regiones (en este caso, subespecies) usando los 20 loci en conjunto:
par(mfrow=c(1,2))
boxplot(HAr_Ttaxus$Ar_pop, ylab="Riqueza alélica")
boxplot(HAr_Ttaxus$Hierar_loc)
dev.off()
Preguntas:
¿Cuáles son los loci más y menos polimórficos? ¿Son los mismos que los reportados por Rico et al. (2016)?
¿Qué poblaciones y regiones son ‘significativamente’ más diversas?
Ahora separamos la diversidad alélica en sus componentes Alfa, Beta y Gamma, así como la diferenciación (Delta D)
Está función descompone el número de Hill q=1 en α, β, γ y Δ
HiD_Ttaxus = HierD(infile_Ttaxus, nreg=3, r=c(1,5,2), ncode=3)print(HiD_Ttaxus)## D_gamma D_alpha.2 D_alpha.1 D_beta.2 D_beta.1 Differentiation.2
## Locus 1 9.663865 7.885412 6.733692 1.225537 1.171038 0.22591250
## Locus 2 9.526010 7.784545 6.089762 1.223708 1.278300 0.22425360
## Locus 3 4.263340 3.946223 3.663383 1.080359 1.077207 0.08585754
## Locus 4 5.892978 5.217391 4.523122 1.129488 1.153493 0.13525487
## Locus 5 8.085174 6.353113 5.391976 1.272632 1.178253 0.26779838
## Locus 6 6.744947 5.534179 4.801177 1.218780 1.152671 0.21977123
## Locus 7 3.815363 3.322561 2.982781 1.148320 1.113914 0.15362293
## Locus 8 5.486055 4.526750 3.991221 1.211919 1.134177 0.21350040
## Locus 9 2.793655 2.579395 2.430767 1.083066 1.061145 0.08863678
## Locus 10 7.394562 6.193308 5.204096 1.193960 1.190083 0.19691671
## Locus 11 9.177483 7.231801 6.109864 1.269045 1.183627 0.26466337
## Locus 12 5.296168 4.722206 4.390728 1.121545 1.075495 0.12741657
## Locus 13 7.000544 6.018859 5.236444 1.163102 1.149417 0.16783021
## Locus 14 8.925683 6.341558 5.298933 1.407491 1.196761 0.37967917
## Locus 15 7.432867 6.229434 4.936175 1.193185 1.261996 0.19619545
## Locus 16 3.916731 3.669226 3.483384 1.067454 1.053351 0.07250904
## Locus 17 7.501690 5.672343 4.658163 1.322503 1.217721 0.31049615
## Locus 18 4.963431 4.520592 3.879349 1.097961 1.165296 0.10380874
## Locus 19 10.627916 8.009273 6.526065 1.326951 1.227275 0.31422621
## Locus 20 4.447186 3.715207 3.330022 1.197022 1.115670 0.19976219
## Differentiation.1
## Locus 1 0.13389830
## Locus 2 0.20822112
## Locus 3 0.06307058
## Locus 4 0.12109636
## Locus 5 0.13910705
## Locus 6 0.12049175
## Locus 7 0.09148683
## Locus 8 0.10677465
## Locus 9 0.05032971
## Locus 10 0.14757927
## Locus 11 0.14296622
## Locus 12 0.06172150
## Locus 13 0.11809424
## Locus 14 0.15232456
## Locus 15 0.19733519
## Locus 16 0.04407838
## Locus 17 0.16704835
## Locus 18 0.12972980
## Locus 19 0.17367574
## Locus 20 0.09282287donde:
D_alpha.1 = diversidad intrapoblacional
D_alpha.2 = diversidad intra-regional
D_beta.1 = diversidad interpoblacional (gamma.1/alfa.1)
D_beta.2 = diversidad inter-regional (gamma.2/alfa.2)
D.gamma = diversidad total en el ecosistema
Differentiation.1 = diferenciación entre poblaciones dentro de una región
Differentiation.2 = diferenciación entre regiones
Graficamos los resultados de la función ‘HiD’
par(mfrow=c(1,3))
boxplot(HiD_Ttaxus[,3:2], ylab="Diversidad alfa")
boxplot(HiD_Ttaxus[,5:4], ylab="Diversidad beta")
boxplot(HiD_Ttaxus[,1], ylab="Diversidad gamma")
dev.off()boxplot(HiD_Ttaxus[,7:6], ylab="Diferenciación (Delta D)")
Comparación con los resultados de Rico et al. (2016)
Ordinación reportada por Rico et al. 
Gráfica de ordinación basada en la diferenciación alélica
infile_Ttaxus <- "Taxidea taxus_genodata_f.gen"
region1=paste("region", rep(1,time=1))
region2=paste("region", rep(2,time=5))
region3=paste("region", rep(3,time=2))
level1=data.frame(matrix(data=0,nrow=8,ncol=1))
level1[1,1]=region1
level1[2:6,1]=region2
level1[7:8,1]=region3
colnames(level1)=c("region")plot_pcoa_aggregate(infile_Ttaxus, ncode=3, level=level1$region, label=TRUE)## Error in Fortran routine computing the spanning ellipsoid,
## probably collinear data## Warning in ellipsoidhull(X): algorithm possibly not converged in 5000
## iterations## Warning in chol.default(cov, pivot = TRUE): the matrix is either rank-
## deficient or indefinite
Aparecen varios errores (colinearidad, deficiencia de rangos y falta de convergencia en el algoritmo). Ajustamos la estructura, separando las poblaciones “LP” y “ON”.
infile_Ttaxus2 <- "Taxidea taxus_genodata_f.gen"
region1=paste("region", rep(1,time=1))
region2=paste("region", rep(2,time=4))
region3=paste("region", rep(3,time=2))
region4=paste("region", rep(4,time=1))
level1=data.frame(matrix(data=0,nrow=8,ncol=1))
level1[1,1]=region1
level1[2:5,1]=region2
level1[6:7,1]=region3
level1[8,1]=region4
colnames(level1)=c("region")plot_pcoa_aggregate(infile_Ttaxus2, ncode=3, level=level1$region, label=TRUE)## Error in Fortran routine computing the spanning ellipsoid,
## probably collinear data## Warning in ellipsoidhull(X): algorithm possibly not converged in 5000
## iterations
La colinearidad desaparece pero aún parece haber error de convergencia. Reajustamos la estructura:
infile_Ttaxus3 <- "Taxidea taxus_genodata_f.gen"
region1=paste("region", rep(1,time=1))
region2=paste("region", rep(2,time=4))
region3=paste("region", rep(3,time=1))
region4=paste("region", rep(4,time=1))
region5=paste("region", rep(5,time=1))
level1=data.frame(matrix(data=0,nrow=8,ncol=1))
level1[1,1]=region1
level1[2:5,1]=region2
level1[6,1]=region3
level1[7,1]=region4
level1[8,1]=region5
colnames(level1)=c("region")plot_pcoa_aggregate(infile_Ttaxus3, ncode=3, level=level1$region, label=TRUE)
No se presentan más errores. Cautionary tale: los resultados NO son erróneos. Solamente no se cuenta con las suficientes poblaciones para una adecuada jerarquización
Interprete el cluster de T.t.taxus comparando el método tradicional con la serie de Hill
Fst
Hier_Fst_Ttaxus=HierFst(infile_Ttaxus, nreg=3, r=c(1,5,2), ncode=3)print(Hier_Fst_Ttaxus)## region pop Ind
## Total 0.07260262 0.11391444 0.2143189
## region 0.00000000 0.04454597 0.1528108
## pop 0.00000000 0.00000000 0.1133124
Graficamos los resultados:
par(mfrow=c(1,3))
barplot(Hier_Fst_Ttaxus[3,1:3], xlab="Intra-poblacional", ylab="Fst", ylim=c(0,0.25))
barplot(Hier_Fst_Ttaxus[2,1:3], xlab="Intra-regional", ylim=c(0,0.25))
barplot(Hier_Fst_Ttaxus[1,1:3], xlab="Total (Ecosistema)", ylim=c(0,0.25))
Analizamos de nuevo, de acuerdo al reajuste en la estructura:
Hier_Fst_Ttaxus3=HierFst(infile_Ttaxus3, nreg=5, r=c(1,4,1,1,1), ncode=3)
print(Hier_Fst_Ttaxus3)## region pop Ind
## Total 0.08743353 0.10827960 0.2093226
## region 0.00000000 0.02284334 0.1335673
## pop 0.00000000 0.00000000 0.1133124
y graficamos:
par(mfrow=c(1,3))
barplot(Hier_Fst_Ttaxus3[3,1:3], xlab="Intra-poblacional", ylab="Fst", ylim=c(0,0.25))
barplot(Hier_Fst_Ttaxus3[2,1:3], xlab="Intra-regional", ylim=c(0,0.25))
barplot(Hier_Fst_Ttaxus3[1,1:3], xlab="Total (Ecosistema)", ylim=c(0,0.25))
dev.off()
Gráfica de diferenciación genética (Delta D) por loci
plotdiff1(infile_Ttaxus,ncode=3)
Gráfica de la matriz de correlación para Delta D por pares
DeltaDcorplot(infile_Ttaxus,ncode=3)
Punto extra a quien consiga el “debug”
Hierarchical heterozygosity (q=2)
He_Ttaxus = HierHe(infile_Ttaxus, nreg=3, r=c(1,5,2), ncode=3)
print(He_Ttaxus)
Ejercicio práctico 2
Clavel silvestre (Dianthus carthusianorum) en Alemania
Fuente: Rico et al. (2014) Molecular Ecology, 23:832-842
Se genotiparon 11 loci microsatélites en 1602 individuos pertenecientes a 59 poblaciones y X regiones
Manos a la obra…