practica_2
LC
2023-05-24
#PRÁCTICA nº 2: Cálculo de la diferenciación genética y métodos de agrupación a partir de datos SNP
Introducción En esta viñeta, discutiremos cómo evaluar la estructura genética poblacional a partir de datos SNP a nivel poblacional. Estimaremos Fst por población, Fst por pares AMOVA (Fst jerárquica ). Finalmente evaluaremos la estructura genética a nivel individual asumiendo que no conocemos las poblaciones usando un análisis multivariante.
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':
##
## amovaImportar datos
Mydata <- read.table("D:/Gene_avan/Genetica de poblaciones/Practica/Master_Pinus_data_genotype.txt", header = T, check.names = FALSE)
dim(Mydata)## [1] 550 3086ind <- as.character(Mydata$tree_id) # use later with adegenet (individual labels)
population <- as.character(Mydata$state) # use later with adegenet (population labels)
county <- Mydata$county
dim(Mydata) # 550 individuals x 3082 SNPs## [1] 550 3086Conversión de datos Para convertir Mydata en un objeto “genind” (adegenet), la entrada sólo debe contener genotipos. Disminuimos el número de SNPs para que los cálculos sean más rápidos y mantenemos sólo 20 SNPs en el locus objeto. A continuación, convertimos Mydata1 en un objeto “hierfstat” (Mydata2).
locus <- Mydata[, 5:24]
Mydata1 <- df2genind(locus, ploidy = 2, ind.names = ind, pop = population, sep = "")
Mydata1## /// GENIND OBJECT /////////
##
## // 550 individuals; 20 loci; 40 alleles; size: 141.5 Kb
##
## // Basic content
## @tab: 550 x 40 matrix of allele counts
## @loc.n.all: number of alleles per locus (range: 2-2)
## @loc.fac: locus factor for the 40 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)Conversión final
Mydata2 <- genind2hierfstat(Mydata1) Heterocigosidad observada y esperada: Fst Estos estadísticos proceden del paquete hierfstat.
basic.stats(Mydata1) # Fst following Nei (1987) on genind object## $perloc
## Ho Hs Ht Dst Htp Dstp Fst Fstp
## X0.10037.01.257 0.4986 0.4079 0.4259 0.0180 0.4277 0.0198 0.0422 0.0462
## X0.10040.02.394 0.4866 0.4971 0.4968 -0.0003 0.4968 -0.0003 -0.0005 -0.0006
## X0.10044.01.392 0.3638 0.4232 0.4931 0.0699 0.5000 0.0768 0.1417 0.1537
## X0.10048.01.60 0.4261 0.4626 0.4953 0.0327 0.4986 0.0359 0.0660 0.0721
## X0.10051.02.166 0.0596 0.0613 0.0634 0.0020 0.0636 0.0023 0.0323 0.0354
## X0.10054.01.402 0.4584 0.4481 0.4761 0.0280 0.4789 0.0308 0.0588 0.0643
## X0.10067.03.111 0.0879 0.0853 0.0853 0.0000 0.0853 -0.0001 -0.0005 -0.0006
## X0.10079.02.168 0.0833 0.0808 0.0830 0.0022 0.0832 0.0024 0.0263 0.0288
## X0.10112.01.169 0.0764 0.0766 0.0760 -0.0006 0.0760 -0.0007 -0.0079 -0.0088
## X0.10113.01.119 0.4436 0.4331 0.4294 -0.0037 0.4290 -0.0041 -0.0087 -0.0095
## X0.10116.01.165 0.0399 0.0407 0.0402 -0.0004 0.0402 -0.0005 -0.0105 -0.0115
## X0.10151.01.86 0.1063 0.1113 0.1125 0.0012 0.1126 0.0013 0.0106 0.0116
## X0.10162.01.255 0.0863 0.0879 0.0862 -0.0018 0.0860 -0.0020 -0.0207 -0.0228
## X0.10207.01.280 0.2875 0.3521 0.3613 0.0092 0.3622 0.0101 0.0254 0.0279
## X0.10210.01.41 0.1081 0.1108 0.1089 -0.0019 0.1087 -0.0021 -0.0174 -0.0194
## X0.10219.01.433 0.3618 0.3648 0.3656 0.0009 0.3657 0.0009 0.0023 0.0026
## X0.1022.02.173 0.4346 0.4267 0.4490 0.0222 0.4512 0.0245 0.0495 0.0542
## X0.10240.01.410 0.4576 0.4530 0.5006 0.0477 0.5054 0.0524 0.0952 0.1037
## X0.10262.01.558 0.2256 0.2292 0.2305 0.0013 0.2306 0.0014 0.0055 0.0060
## X0.10266.01.426 0.0677 0.0720 0.0716 -0.0004 0.0715 -0.0005 -0.0059 -0.0064
## Fis Dest
## X0.10037.01.257 -0.2224 0.0334
## X0.10040.02.394 0.0210 -0.0006
## X0.10044.01.392 0.1403 0.1332
## X0.10048.01.60 0.0790 0.0669
## X0.10051.02.166 0.0286 0.0024
## X0.10054.01.402 -0.0229 0.0558
## X0.10067.03.111 -0.0293 -0.0001
## X0.10079.02.168 -0.0308 0.0026
## X0.10112.01.169 0.0036 -0.0007
## X0.10113.01.119 -0.0243 -0.0072
## X0.10116.01.165 0.0193 -0.0005
## X0.10151.01.86 0.0443 0.0015
## X0.10162.01.255 0.0192 -0.0021
## X0.10207.01.280 0.1834 0.0156
## X0.10210.01.41 0.0245 -0.0024
## X0.10219.01.433 0.0082 0.0015
## X0.1022.02.173 -0.0185 0.0427
## X0.10240.01.410 -0.0102 0.0958
## X0.10262.01.558 0.0160 0.0018
## X0.10266.01.426 0.0591 -0.0005
##
## $overall
## Ho Hs Ht Dst Htp Dstp Fst Fstp Fis Dest
## 0.2580 0.2612 0.2725 0.0113 0.2737 0.0124 0.0415 0.0454 0.0124 0.0168Estimación de Weir y Cockerham
wc(Mydata1) # Weir and Cockerham's estimate## $FST
## [1] 0.02300324
##
## $FIS
## [1] 0.03090781Pruebas Fst jerárquicas (=AMOVA para conjunto de datos SNP) La función varcomp.glob() produce un Fst Jerárquico (=AMOVA para SNPs o marcadores bialélicos) Es posible hacer permutaciones en los diferentes niveles: La función test.g() comprueba el efecto de la población en la diferenciación genética. Los individuos se permutan aleatoriamente entre estados. Los estados influyen en la diferenciación genética a un nivel del 5%. Con la función test.between(), los condados se permutan entre estados. Los estados influyen significativamente en la estructuración genética.
loci <- Mydata2[, -1] # Remove the population column
varcomp.glob(levels = data.frame(population, county), loci, diploid = TRUE) ## $loc
## [,1] [,2] [,3] [,4]
## X0.10037.01.257 4.631785e-03 -0.0075801286 -0.0110876895 0.42329020
## X0.10040.02.394 -3.927184e-06 0.0057909958 0.0035980465 0.48979592
## X0.10044.01.392 1.276810e-02 0.0039247749 0.0247870243 0.46198830
## X0.10048.01.60 6.490717e-03 0.0189152357 0.0255163684 0.39741220
## X0.10051.02.166 2.977513e-03 0.0037478140 -0.0001317105 0.11970534
## X0.10054.01.402 1.806575e-02 -0.0002977277 0.0134753262 0.47329650
## X0.10067.03.111 9.490247e-04 -0.0022733109 0.0022479222 0.07339450
## X0.10079.02.168 1.359482e-03 0.0015693361 -0.0020878338 0.07692308
## X0.10112.01.169 5.570862e-04 -0.0011207466 -0.0022195075 0.11151737
## X0.10113.01.119 -3.038733e-03 0.0190331181 -0.0321043585 0.45871560
## X0.10116.01.165 -4.102906e-04 0.0010983394 0.0019041211 0.04204753
## X0.10151.01.86 1.180089e-03 0.0038425248 0.0163965166 0.14180479
## X0.10162.01.255 2.762257e-04 0.0025535687 -0.0007521869 0.09778598
## X0.10207.01.280 2.419336e-02 0.0016539601 0.0286306732 0.38051471
## X0.10210.01.41 -1.206136e-03 0.0040434594 0.0069081606 0.09775967
## X0.10219.01.433 2.608653e-03 0.0035115484 0.0048812171 0.30755064
## X0.1022.02.173 7.258406e-03 0.0006657100 -0.0222437933 0.41198502
## X0.10240.01.410 3.603309e-02 0.0212763776 0.0035015937 0.41263941
## X0.10262.01.558 5.048435e-04 -0.0006787502 0.0310839391 0.15201465
## X0.10266.01.426 -7.102129e-04 0.0020016810 0.0053850735 0.11151737
##
## $overall
## population county Ind Error
## 0.11448482 0.08167778 0.09768890 5.24165877
##
## $F
## population county Ind
## Total 0.02068189 0.03543713 0.05308481
## population 0.00000000 0.01506685 0.03308722
## county 0.00000000 0.00000000 0.01829604test.g(loci, level = population) ## $g.star
## [1] 205.0137 261.2203 240.0722 204.4959 217.7883 232.8068 225.4453 219.3582
## [9] 228.3732 244.3867 247.1552 233.3499 205.4381 206.9693 265.6322 194.8710
## [17] 244.1548 239.5045 204.1266 228.7004 229.0208 206.4280 197.9791 227.7793
## [25] 201.9512 264.5635 213.7352 227.8742 235.4000 203.4867 208.4111 223.6889
## [33] 225.0723 259.9158 222.9178 222.1141 185.9082 247.8932 203.5019 210.8898
## [41] 241.0161 226.1597 202.5719 228.1248 264.4577 246.7554 263.5900 189.8385
## [49] 236.2845 218.3503 218.9117 252.6959 231.3135 227.5882 216.1272 246.3573
## [57] 251.3572 247.2888 237.5621 182.3552 203.8249 184.6740 220.7069 225.4622
## [65] 249.3244 202.1369 239.8745 212.6512 197.9028 206.8366 251.7288 183.8597
## [73] 219.2770 231.4649 228.2586 221.6990 202.8121 215.1941 192.5580 178.9972
## [81] 206.8197 214.8628 191.6519 230.3396 238.4507 208.7480 226.9017 217.6387
## [89] 221.2079 220.1301 209.7764 254.0793 208.0785 259.4231 218.2337 212.3895
## [97] 229.5584 240.9149 221.9940 378.7654
##
## $p.val
## [1] 0.01test.between(loci, test.lev = population, rand.unit = county, nperm = 100) ## $g.star
## [1] 226.3518 216.3278 240.2997 234.4399 259.3469 207.5079 263.4204 275.9181
## [9] 280.7377 265.6555 260.3884 252.8274 234.4215 288.0257 229.2636 258.8305
## [17] 252.5587 273.5986 240.7157 268.7078 272.8535 232.9686 276.0974 305.5461
## [25] 262.9456 218.3979 294.4636 265.8197 262.4904 254.2350 266.8151 249.1796
## [33] 242.5928 291.3419 265.8925 298.4618 240.2064 225.4393 297.7730 263.5586
## [41] 294.0938 267.5754 236.3846 295.2488 255.8374 251.7269 248.8596 246.8049
## [49] 294.0309 235.1101 259.8851 235.1804 270.5308 276.7476 269.6473 269.8785
## [57] 211.6253 242.6638 271.4075 269.9564 215.6159 238.1041 279.8214 294.0705
## [65] 240.7484 233.9395 310.6488 279.7521 245.6695 253.1752 215.4911 270.7321
## [73] 271.6947 233.7996 263.3020 275.6045 256.5686 266.6917 242.6900 296.5988
## [81] 279.0743 267.6146 289.4335 225.1005 280.2719 228.8582 228.6465 257.9564
## [89] 261.4010 243.5127 286.9465 225.0518 273.9091 262.8850 299.5269 280.7104
## [97] 268.1913 215.2070 261.4161 378.7654
##
## $p.val
## [1] 0.01Fst por pares
genet.dist(Mydata1, method = "WC84")## georgia virginia northcarolina southcarolina
## virginia 0.014916721
## northcarolina 0.035531517 0.012476717
## southcarolina 0.029831181 0.009878778 0.006809538
## alabama 0.046469226 0.027303304 0.009066562 0.007107402
## oklahoma 0.092718257 0.125242248 0.116089566 0.107662986
## arkansas 0.028847476 0.065029859 0.057810705 0.052707863
## florida 0.001379763 0.027697122 0.033602798 0.019836955
## mississippi 0.018627828 0.034121040 0.030048967 0.030695068
## texas 0.046361901 0.068423280 0.046659978 0.052491099
## louisiana 0.133824385 0.092337839 0.043365828 0.064729942
## alabama oklahoma arkansas florida mississippi
## virginia
## northcarolina
## southcarolina
## alabama
## oklahoma 0.104844418
## arkansas 0.052033165 0.022743098
## florida 0.034205805 0.080741693 0.031222010
## mississippi 0.019629941 0.061834123 -0.011045367 0.008707686
## texas 0.038153327 0.033860376 0.013982640 0.025059824 -0.025535193
## louisiana 0.031028123 0.070815392 0.035722729 0.096253248 0.030209963
## texas
## virginia
## northcarolina
## southcarolina
## alabama
## oklahoma
## arkansas
## florida
## mississippi
## texas
## louisiana 0.036298136Agrupación no supervisada No conocemos las poblaciones que buscamos. Como recomienda T. Jombart, con la función find.clusters() utilizamos el máximo número posible de ejes PCA que aquí es 20. Ver tutorial detallado sobre este método para más información (https://github.com/thibautjombart/adegenet/raw/master/tutorials /tutorial-basics.pdf) En este ejemplo, hemos utilizado choose.n.clust = FALSE pero es bueno utilizar la opción TRUE y entonces usted será capaz de elegir el número de clusters
# using Kmeans and DAPC in adegenet
set.seed(20160308) # Setting a seed for a consistent result
grp <- find.clusters(Mydata1, max.n.clust = 10, n.pca = 20, choose.n.clust = FALSE)
names(grp)## [1] "Kstat" "stat" "grp" "size"grp$grp## 1066 2040 4004 4005 4018 6009 6013 7033 7056 7069 7088
## 1 1 2 1 1 1 2 2 2 2 1
## 7105 8001 8061 8068 8076 8120 8195 8203 8222 8223 8231
## 2 1 1 1 1 1 1 1 2 1 1
## 8237 8301 8302 8303 8304 8305 8307 8308 8309 8310 8313
## 2 1 1 1 1 2 2 1 1 1 1
## 8314 8316 8317 8318 8319 8320 8323 8324 8327 8328 8329
## 1 1 1 2 1 1 1 1 1 2 1
## 8330 8332 8333 8334 8335 8336 8337 8338 8339 8340 8342
## 2 1 1 2 1 2 2 2 1 1 1
## 8343 8344 8345 8346 8347 8349 8350 8351 8352 8353 8354
## 1 1 1 1 2 1 1 1 1 2 1
## 8355 8356 8357 8358 8364 8365 8366 8367 8368 8369 8370
## 1 1 1 1 1 1 1 2 1 1 1
## 8371 8372 8373 8374 8375 8376 8377 8378 8379 8381 8382
## 1 1 1 1 2 2 1 1 1 2 1
## 8383 8384 8385 8386 8387 8388 8389 8390 8391 8392 8393
## 2 2 2 1 1 2 1 1 2 2 2
## 8395 8400 8402 8403 8559 8565 8567 8568 8569 8570 8571
## 2 1 2 2 1 2 2 2 1 1 1
## 8572 8573 8574 8601 8602 8603 8604 8606 8607 8608 8609
## 2 1 1 1 1 2 2 2 1 1 1
## 8610 8611 8613 8615 8616 8618 8619 8620 8621 8622 8624
## 2 1 2 2 1 2 2 1 2 1 2
## 8626 8628 8629 8630 8631 8633 8634 8635 8636 8637 8638
## 2 1 1 1 1 1 2 1 2 1 1
## 8639 8640 8641 8642 8643 8644 8645 8646 8647 8648 8651
## 2 2 2 1 2 1 2 1 2 1 1
## 8652 8653 8654 8655 8656 8657 8658 8659 8660 8661 8662
## 1 1 1 2 1 2 1 1 1 1 2
## 8663 8664 8667 8669 8670 8671 8672 8673 8674 8675 8676
## 1 2 1 1 2 1 2 2 2 1 1
## 8677 8678 8679 8680 8683 8685 8686 8687 8688 8689 8691
## 1 2 2 1 2 2 1 2 1 1 2
## 8692 8693 8694 8695 8696 8697 8698 8699 8700 8701 8702
## 2 2 1 2 2 2 2 1 2 1 1
## 8703 8704 8705 8706 9003 9006 9015 10005 11010 11503 11532
## 2 1 2 2 1 1 2 1 1 1 2
## 12008 12012 14010 14015 22212 68087 68088 68090 68095 68130 68131
## 2 1 2 1 1 1 1 2 2 1 2
## 68133 68134 68135 105A 108A 109B 110B 112C 115B 117B 118B
## 1 1 1 2 2 1 1 1 1 1 2
## 11A 120A 121C 127A 128A 131B 132B 136C 138A 139B 140B
## 1 1 1 2 2 2 2 2 1 2 2
## 141A 142B 144C 145A 146C 147A 149B 150A 151A 152B 153B
## 2 2 2 2 1 2 2 1 2 2 2
## 154C 155B 156C 157A 158B 15A 162A 166B 16A 171A 173A
## 2 2 1 2 1 1 1 2 1 1 2
## 174A 17C 188A 189A 190A 191A 19A 205B 20A 212A 213A
## 2 2 1 2 2 2 2 2 2 1 1
## 217C 219A 21A 220A 224A 226C 227C 234B 235A 238A 23A
## 2 1 2 1 1 2 2 1 2 2 2
## 245B 248A 250C 253A 254C 257B 258A 260B 262A 264C 265A
## 1 2 1 1 2 1 2 2 1 1 1
## 268A 269A 270C 271A 272B 275A 276B 277A 27A 281A 282B
## 1 1 2 2 1 2 1 2 1 1 2
## 283C 285C 286B 287C 288C 289A 290C 291C 292C 298B 299C
## 1 1 2 1 1 2 2 1 1 2 1
## 300C 302A 303C 305A 306A 307A 311A 316B 320C 322A 323B
## 2 2 1 2 1 2 1 2 2 1 2
## 324A 326A 327A 328B 329B 32A 330A 331B 332C 334A 335A
## 1 2 1 2 2 2 2 2 2 2 1
## 336A 339B 340A 341C 346A 349B 34A 351A 353C 355A 35A
## 2 2 1 1 2 1 2 2 2 1 1
## 360B 361B 362C 365B 366A 368B 369A 36A 370C 371B 372A
## 2 2 1 1 1 2 2 1 2 1 1
## 373B 375A 377B 378B 379B 382A 383C 384A 385B 387A 388B
## 1 2 2 2 2 2 2 2 2 1 1
## 389A 390B 391C 392A 393C 395A 397A 398C 400C 407A 408C
## 1 1 2 2 2 2 2 1 1 2 2
## 409B 410A 411C 412C 414A 415A 416B 417A 418A 419B 41C
## 2 1 2 2 2 2 1 2 1 1 2
## 420B 421A 422B 423C 424B 425B 426B 427C 428C 429B 42A
## 1 2 1 1 1 1 2 2 1 1 2
## 430B 431B 433B 434B 435C 436A 43B 441C 442C 443C 448C
## 2 2 1 2 2 2 1 2 1 1 1
## 449A 450B 451B 459A 461A 463A 469C 470A 471B 481A 483A
## 2 1 1 2 2 1 1 2 2 1 1
## 484A 485A 486B 487C 489B 48B 492C 493A 496B 498B 499A
## 1 1 1 1 2 2 1 1 2 1 2
## 49A 500B 501A 502C 514A 515A 519B 520B 526A 527B 528C
## 2 1 2 2 2 1 2 1 2 1 2
## 52A 531A 532A 533A 534A 535C 536A 539A 53C 540A 541A
## 2 2 2 2 1 1 2 2 2 2 2
## 542C 543C 544B 545A 546C 548C 549A 54C 551C 552A 553B
## 2 2 1 2 1 2 2 1 1 1 2
## 554A 555C 556A 557B 558A 559A 55A 561A 562A 563A 564B
## 1 2 2 2 2 2 2 2 2 2 2
## 565C 566B 568B 570A 571A 572C 573C 574B 576A 577B 578B
## 1 2 1 1 2 2 2 2 2 2 2
## 579A 57A 580A 581C 600A 601A 603A 605B 606B 60A 612C
## 1 2 2 1 1 2 1 1 1 1 1
## 613A 618A 619A 61A 620C 621A 633B 634C 635A 636C 637B
## 1 2 2 1 1 1 1 2 1 1 1
## 63B 644B 645A 646A 66A 67C 69A 73B 77B 7A 86C
## 1 1 1 1 2 2 2 1 2 1 1
## 89C 90C 92A 93C 94C 97B 98A 99C 9A CRO108 CRO120
## 1 2 2 2 2 1 1 2 2 2 1
## CRO121 CRO133 DF3364 FM406 FM417 FM428 FM442 FM445 S4PT6 SH13 SH7
## 1 1 1 1 1 2 1 2 1 2 1
## Levels: 1 2El procedimiento de medias K detectó 4 grupos. Utilizaremos este número de grupo en el análisis discriminante (función dapc()). En su propio conjunto de datos, tendrá que dedicar más tiempo a estimar el número de grupos.
dapc1 <- dapc(Mydata1, grp$grp, n.pca = 20, n.da = 6)
scatter(dapc1) # plot of the group
Está claro que un subconjunto de 20 SNP no tiene una señal lo
suficientemente fuerte como para separar las muestras en grupos
distintos. ¿Qué pasaría si utilizáramos más SNPs?
Conclusiones ¿Qué hemos aprendido hoy? En esta viñeta, aprendimos a calcular Fst en poblaciones existentes y a investigar el efecto de la estructura de la población en la diferenciación genética a partir del análisis Fst jerárquico (como AMOVA en el caso de SNP). También ejecutamos un análisis multivariante para investigar la estructura genética de los datos a nivel individual asumiendo que no hay estructura poblacional.
Traducción realizada con la versión gratuita del traductor www.DeepL.com/Translator