Reading data
dat1 <- read.vcfR("md6ct95ms200.vcf")
## Scanning file to determine attributes.
## File attributes:
## meta lines: 10
## header_line: 11
## variant count: 25099
## column count: 326
## Meta line 10 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
## Character matrix gt rows: 25099
## Character matrix gt cols: 326
## skip: 0
## nrows: 25099
## row_num: 0
## Processed variant 1000Processed variant 2000Processed variant 3000Processed variant 4000Processed variant 5000Processed variant 6000Processed variant 7000Processed variant 8000Processed variant 9000Processed variant 10000Processed variant 11000Processed variant 12000Processed variant 13000Processed variant 14000Processed variant 15000Processed variant 16000Processed variant 17000Processed variant 18000Processed variant 19000Processed variant 20000Processed variant 21000Processed variant 22000Processed variant 23000Processed variant 24000Processed variant 25000Processed variant: 25099
## All variants processed
## remove clones detected on phylogenetic tree
clones_r <- c(
"DE1-4","DE1-7","DE1-5",
"DE2-4",
"DE3-9","DE3-10","DE3-11","DE3-5",
"DE4-1","DE4-3","DE4-6","DE4-8",
"DE5-3","DE5-2","DE5-7","DE5-9",
"DE7-12","DE7-6","DE7-7","DE7-10","DE7-2","DE7-17","DE7-18","DE7-19","DE7-20", "DE7-8",
"DE9-13", "DE9-2r",
"DE10-11","DE10-12","DE10-2","DE10-5","DE10-7","DE10-8",
"DE11-12","DE11-15","DE11-14r", "DE12-5",
"DK-4","DK-5",
"EE-7",
"FI1-7-1","FI5-2","FI5-10","FI5-8",
"LT1-9",
"LT2-3","LT2-4","LT2-5","LT2-6","LT2-7", "LT2-1r",
"LT3-11","LT3-12","LT3-5",
"LV-12","LV-8",
"PL2-4","PL2-6-2",
"PL4-14","PL4-1","PL4-3","PL4-7",
"PL5-20","PL5-21",
"PL6-10","PL6-9",
"PL7-12","PL7-5","PL7-18",
"PL8-17","PL8-8","PL9-4","PL9-8", "PL8-19-",
"RU1-4","RU2-3","RU2-5","RU2-6",
"SE3-4-2","SE3-4-2r", "SE1-4","SE1-5","SE1-6","SE1-7","SE1-9",
"UA1-13","UA1-14","UA1-15","UA1-17", "UA1-16r",
"UA2-1","UA2-12","UA2-18","UA2-3")
dat1.pre_gl <- !colnames(dat1@gt)[-1] %in% clones_r
keep <- !colnames(dat1@gt)[-1] %in% clones_r
dat1@gt <- dat1@gt[, c(TRUE, keep)]
dat1.pre_gl <- vcfR2genlight(dat1)
## Warning in vcfR2genlight(dat1): Found 245 loci with more than two alleles.
## Objects of class genlight only support loci with two alleles.
## 245 loci will be omitted from the genlight object.
ploidy(dat1.pre_gl) <- 2
dat1.pre_gl <- gl.compliance.check(dat1.pre_gl)
## Starting gl.compliance.check
## Processing genlight object with SNP data
## Warning: data include loci that are scored NA across all individuals.
## Consider filtering using gl <- gl.filter.allna(gl)
## The slot loc.all, which stores allele name for each locus, is empty.
## Creating a dummy variable (A/C) to insert in this slot.
## Checking coding of SNPs
## SNP data scored NA, 0, 1 or 2 confirmed
## Checking locus metrics and flags
## Recalculating locus metrics
## Checking for monomorphic loci
## Dataset contains monomorphic loci
## Checking for loci with all missing data
## Dataset contains loci with all missing dat
## Checking whether individual names are unique.
## Checking for individual metrics
## Warning: Creating a slot for individual metrics
## Checking for population assignments
## Population assignments not detected, individuals assigned
## to a single population labelled 'pop1'
## Spelling of coordinates checked and changed if necessary to
## lat/lon
## Completed: gl.compliance.check
dat1.pre_gl <- gl.filter.allna(dat1.pre_gl)
## Starting gl.filter.allna
## Processing genlight object with SNP data
## Identifying and removing loci and individuals scored all
## missing (NA)
## Deleting loci that are scored as all missing (NA)
## Deleted
## Deleting individuals that are scored as all missing (NA)
## Completed: gl.filter.allna
DAPC K=2
Note: DAPC randomly assigns numbers to the clusters during each
run
#grp <- find.clusters(dat1.pre_gl, max.n.clust = 20)
# with the function above (see output in Supplemetary Figure S2) I found that K=2 and 3 the best for describing the data
grp <- find.clusters(dat1.pre_gl, n.pca = 1000, n.clust = 2, choose.n.clust = FALSE)
dapc_res <- dapc(dat1.pre_gl, grp$grp, n.pca = 100, n.da = 5)
var_exp <- dapc_res$eig / sum(dapc_res$eig) * 100
var_exp
## [1] 100
dapc_res$grp
## DE1-2 DE10-1 DE10-10 DE10-3 DE11-1 DE11-10 DE11-11 DE11-14 DE11-16 DE11-17
## 2 1 2 1 1 1 1 1 1 1
## DE11-2 DE11-3 DE11-4 DE11-5 DE11-6 DE11-9 DE12-1 DE12-11 DE12-12 DE12-14
## 1 1 1 1 1 1 1 1 1 1
## DE12-2 DE12-3 DE12-4 DE12-6 DE12-7 DE12-8 DE12-9 DE2-1 DE2-2 DE2-3
## 1 1 1 1 1 1 1 2 2 2
## DE3-1 DE3-2 DE3-3 DE3-4 DE3-8 DE4-10 DE4-11 DE4-2 DE4-4 DE4-5
## 2 2 2 2 2 2 2 2 2 2
## DE5-1 DE5-10 DE5-4 DE5-6 DE6-4 DE6-5 DE6-6 DE6-7 DE7-1 DE7-11
## 2 2 2 2 1 1 2 2 1 1
## DE7-13 DE7-23 DE7-25 DE7-3 DE7-4 DE8-1 DE8-6 DE8-7 DE8-8 DE8-9
## 1 1 1 1 1 1 1 1 1 1
## DE9-1 DE9-10 DE9-11 DE9-12 DE9-15 DE9-16 DE9-2 DE9-4 DE9-6 DE9-7
## 1 1 1 1 1 1 1 1 1 1
## DE9-9 DK-8 DK-9 EE-5 EE-6 FI1-1 FI1-12 FI1-14 FI1-15 FI1-16
## 1 2 2 1 1 1 1 1 1 1
## FI1-17 FI1-4 FI1-5-1 FI1-5-2 FI1-6 FI1-7-2 FI1-8 FI2-10 FI2-13 FI2-15
## 1 1 1 1 1 1 1 1 1 1
## FI2-17 FI2-19 FI2-20 FI2-3 FI2-6 FI3-14 FI3-17 FI3-19 FI3-21 FI3-24
## 1 1 1 1 1 1 1 1 1 1
## FI3-3 FI3-5 FI3-9 FI4-1 FI4-13 FI4-16 FI4-2 FI4-22 FI4-24 FI4-27
## 1 1 1 1 1 1 1 1 1 1
## FI4-3 FI4-32 FI5-1 FI5-12 FI5-4 FI5-5 FI5-7 FI5-9 LT1-10 LT1-11
## 1 1 1 1 1 1 1 1 1 1
## LT1-4 LT1-5 LT1-6 LT2-1 LT2-2 LT3-10 LT3-14 LT3-7 LT3-9 LV-10
## 1 1 1 1 1 1 1 1 1 1
## LV-9 PL1-10 PL1-5 PL1-6 PL1-7 PL1-8 PL1-9 PL2-10 PL2-2 PL2-3
## 1 1 1 1 1 1 1 1 1 1
## PL2-5 PL2-6_1 PL4-10 PL4-11 PL4-13 PL4-15 PL4-2 PL4-5 PL4-8 PL5-15
## 1 1 1 1 1 1 1 1 1 1
## PL5-19 PL5-4 PL6-1 PL6-13 PL6-3 PL6-4 PL6-5 PL6-6 PL6-7 PL7-10
## 1 1 1 1 1 1 1 1 1 1
## PL7-13 PL7-14 PL7-17 PL7-7 PL7-8 PL7-9 PL8-1 PL8-11 PL8-12 PL8-14
## 1 1 1 1 1 1 1 1 1 1
## PL8-15 PL8-16 PL8-18 PL8-19 PL8-2 PL8-4 PL8-5 PL8-6 PL8-7 PL9-1
## 1 1 1 1 1 1 1 1 1 1
## PL9-10 PL9-13 PL9-5 PL9-6 PL9-7 RU1-1 RU1-2 RU1-3 RU1-5 RU1-6
## 1 1 1 1 1 1 1 1 1 1
## RU2-1 RU2-2 RU2-4 SE1-1 SE1-10 SE2-1 SE2-10 SE2-2 SE2-3 SE2-6
## 1 1 1 1 1 1 1 1 1 1
## SE2-8 SE3-1 SE3-10 SE3-2-1 SE3-2-2 SE3-4-1 SE3-5 SE3-8-1 SE3-9 UA1-1
## 1 1 1 1 1 1 1 1 1 1
## UA1-10 UA1-12 UA1-16 UA1-3 UA1-7 UA2-11 UA2-14 UA2-16 UA2-19 UA2-2
## 1 1 1 1 1 1 1 1 1 1
## UA2-5 UA2-7
## 1 1
## Levels: 1 2
cluster_colors <- c("#009988", "#EE7733")
scatter.dapc(dapc_res,
bg = "white",
col = cluster_colors,
pch = 19,
cstar = 0, solid = 0.75,
scree.da = TRUE,
scree.pca =TRUE,
posi.pca = "topright",
posi.da = "bottomright",
label.inds = list(air = 0.8, pch = NA, col = "gray40", cex = 0.6, pos = 3))

DAPC K=3
grp3 <- find.clusters(dat1.pre_gl, n.pca = 1000, n.clust = 3, choose.n.clust = FALSE)
dapc_res3 <- dapc(dat1.pre_gl, grp3$grp, n.pca = 100, n.da = 3)
var_exp <- dapc_res3$eig / sum(dapc_res3$eig) * 100
var_exp
## [1] 81.47256 18.52744
dapc_res3$grp
## DE1-2 DE10-1 DE10-10 DE10-3 DE11-1 DE11-10 DE11-11 DE11-14 DE11-16 DE11-17
## 1 3 1 2 3 2 2 2 2 2
## DE11-2 DE11-3 DE11-4 DE11-5 DE11-6 DE11-9 DE12-1 DE12-11 DE12-12 DE12-14
## 2 3 2 2 3 2 2 2 2 2
## DE12-2 DE12-3 DE12-4 DE12-6 DE12-7 DE12-8 DE12-9 DE2-1 DE2-2 DE2-3
## 2 3 3 2 2 3 2 1 1 1
## DE3-1 DE3-2 DE3-3 DE3-4 DE3-8 DE4-10 DE4-11 DE4-2 DE4-4 DE4-5
## 1 1 1 1 1 1 1 1 1 1
## DE5-1 DE5-10 DE5-4 DE5-6 DE6-4 DE6-5 DE6-6 DE6-7 DE7-1 DE7-11
## 1 1 1 1 3 2 1 1 2 3
## DE7-13 DE7-23 DE7-25 DE7-3 DE7-4 DE8-1 DE8-6 DE8-7 DE8-8 DE8-9
## 3 2 2 2 2 2 2 2 2 2
## DE9-1 DE9-10 DE9-11 DE9-12 DE9-15 DE9-16 DE9-2 DE9-4 DE9-6 DE9-7
## 2 2 2 3 2 2 3 2 2 2
## DE9-9 DK-8 DK-9 EE-5 EE-6 FI1-1 FI1-12 FI1-14 FI1-15 FI1-16
## 3 1 1 2 2 2 2 2 2 2
## FI1-17 FI1-4 FI1-5-1 FI1-5-2 FI1-6 FI1-7-2 FI1-8 FI2-10 FI2-13 FI2-15
## 2 2 2 2 2 2 2 2 2 2
## FI2-17 FI2-19 FI2-20 FI2-3 FI2-6 FI3-14 FI3-17 FI3-19 FI3-21 FI3-24
## 2 2 2 2 2 2 2 2 2 2
## FI3-3 FI3-5 FI3-9 FI4-1 FI4-13 FI4-16 FI4-2 FI4-22 FI4-24 FI4-27
## 2 2 2 2 3 3 2 2 2 2
## FI4-3 FI4-32 FI5-1 FI5-12 FI5-4 FI5-5 FI5-7 FI5-9 LT1-10 LT1-11
## 2 2 2 3 3 3 3 3 3 2
## LT1-4 LT1-5 LT1-6 LT2-1 LT2-2 LT3-10 LT3-14 LT3-7 LT3-9 LV-10
## 2 3 2 3 3 3 3 3 3 2
## LV-9 PL1-10 PL1-5 PL1-6 PL1-7 PL1-8 PL1-9 PL2-10 PL2-2 PL2-3
## 2 2 2 2 2 2 2 3 3 2
## PL2-5 PL2-6_1 PL4-10 PL4-11 PL4-13 PL4-15 PL4-2 PL4-5 PL4-8 PL5-15
## 3 3 2 3 3 2 3 2 2 3
## PL5-19 PL5-4 PL6-1 PL6-13 PL6-3 PL6-4 PL6-5 PL6-6 PL6-7 PL7-10
## 3 3 3 3 3 3 3 3 3 3
## PL7-13 PL7-14 PL7-17 PL7-7 PL7-8 PL7-9 PL8-1 PL8-11 PL8-12 PL8-14
## 3 3 3 3 3 3 3 3 3 3
## PL8-15 PL8-16 PL8-18 PL8-19 PL8-2 PL8-4 PL8-5 PL8-6 PL8-7 PL9-1
## 3 3 3 3 3 3 3 3 3 3
## PL9-10 PL9-13 PL9-5 PL9-6 PL9-7 RU1-1 RU1-2 RU1-3 RU1-5 RU1-6
## 3 3 3 3 3 3 2 3 3 3
## RU2-1 RU2-2 RU2-4 SE1-1 SE1-10 SE2-1 SE2-10 SE2-2 SE2-3 SE2-6
## 2 3 2 2 2 3 3 3 3 3
## SE2-8 SE3-1 SE3-10 SE3-2-1 SE3-2-2 SE3-4-1 SE3-5 SE3-8-1 SE3-9 UA1-1
## 3 2 3 2 3 2 3 3 2 3
## UA1-10 UA1-12 UA1-16 UA1-3 UA1-7 UA2-11 UA2-14 UA2-16 UA2-19 UA2-2
## 3 3 3 3 3 3 3 3 3 3
## UA2-5 UA2-7
## 3 3
## Levels: 1 2 3
cluster_colors <- c("#0077BB","#009988","#EE7733" )
scatter.dapc(dapc_res3,
bg = "white",
col = cluster_colors,
pch = 19,
cstar = 0, solid = 0.75,
scree.da = TRUE,
scree.pca =TRUE,
posi.pca = "topright",
posi.da = "bottomright",
label.inds = list(air = 0.8, pch = NA, col = "gray40", cex = 0.6, pos = 3))
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "pos" is not a graphical
## parameter

DAPC K=4
grp4 <- find.clusters(dat1.pre_gl, n.pca = 1000, n.clust = 4, choose.n.clust = FALSE)
dapc_res4 <- dapc(dat1.pre_gl, grp4$grp, n.pca = 100, n.da = 5)
dapc_res4$grp
## DE1-2 DE10-1 DE10-10 DE10-3 DE11-1 DE11-10 DE11-11 DE11-14 DE11-16 DE11-17
## 3 2 3 4 4 4 4 4 4 4
## DE11-2 DE11-3 DE11-4 DE11-5 DE11-6 DE11-9 DE12-1 DE12-11 DE12-12 DE12-14
## 4 4 4 4 2 4 4 4 4 4
## DE12-2 DE12-3 DE12-4 DE12-6 DE12-7 DE12-8 DE12-9 DE2-1 DE2-2 DE2-3
## 4 1 4 4 4 1 4 3 3 3
## DE3-1 DE3-2 DE3-3 DE3-4 DE3-8 DE4-10 DE4-11 DE4-2 DE4-4 DE4-5
## 3 3 3 3 3 3 3 3 3 3
## DE5-1 DE5-10 DE5-4 DE5-6 DE6-4 DE6-5 DE6-6 DE6-7 DE7-1 DE7-11
## 3 3 3 3 2 4 3 3 4 4
## DE7-13 DE7-23 DE7-25 DE7-3 DE7-4 DE8-1 DE8-6 DE8-7 DE8-8 DE8-9
## 4 4 4 4 4 4 4 4 4 4
## DE9-1 DE9-10 DE9-11 DE9-12 DE9-15 DE9-16 DE9-2 DE9-4 DE9-6 DE9-7
## 4 4 4 4 4 4 4 4 4 4
## DE9-9 DK-8 DK-9 EE-5 EE-6 FI1-1 FI1-12 FI1-14 FI1-15 FI1-16
## 4 3 3 4 4 4 4 4 4 4
## FI1-17 FI1-4 FI1-5-1 FI1-5-2 FI1-6 FI1-7-2 FI1-8 FI2-10 FI2-13 FI2-15
## 4 4 4 4 4 4 2 4 4 4
## FI2-17 FI2-19 FI2-20 FI2-3 FI2-6 FI3-14 FI3-17 FI3-19 FI3-21 FI3-24
## 4 4 4 4 4 4 4 4 4 4
## FI3-3 FI3-5 FI3-9 FI4-1 FI4-13 FI4-16 FI4-2 FI4-22 FI4-24 FI4-27
## 4 4 4 4 2 2 4 4 4 4
## FI4-3 FI4-32 FI5-1 FI5-12 FI5-4 FI5-5 FI5-7 FI5-9 LT1-10 LT1-11
## 4 4 2 2 2 2 2 2 2 4
## LT1-4 LT1-5 LT1-6 LT2-1 LT2-2 LT3-10 LT3-14 LT3-7 LT3-9 LV-10
## 4 2 4 2 2 2 2 2 2 4
## LV-9 PL1-10 PL1-5 PL1-6 PL1-7 PL1-8 PL1-9 PL2-10 PL2-2 PL2-3
## 4 4 4 4 4 4 4 4 2 4
## PL2-5 PL2-6_1 PL4-10 PL4-11 PL4-13 PL4-15 PL4-2 PL4-5 PL4-8 PL5-15
## 1 4 4 1 1 1 1 4 4 2
## PL5-19 PL5-4 PL6-1 PL6-13 PL6-3 PL6-4 PL6-5 PL6-6 PL6-7 PL7-10
## 2 4 2 1 2 2 1 1 2 2
## PL7-13 PL7-14 PL7-17 PL7-7 PL7-8 PL7-9 PL8-1 PL8-11 PL8-12 PL8-14
## 2 2 2 2 1 2 1 1 1 1
## PL8-15 PL8-16 PL8-18 PL8-19 PL8-2 PL8-4 PL8-5 PL8-6 PL8-7 PL9-1
## 1 1 1 1 1 1 1 1 1 1
## PL9-10 PL9-13 PL9-5 PL9-6 PL9-7 RU1-1 RU1-2 RU1-3 RU1-5 RU1-6
## 1 1 1 1 1 4 4 4 4 4
## RU2-1 RU2-2 RU2-4 SE1-1 SE1-10 SE2-1 SE2-10 SE2-2 SE2-3 SE2-6
## 4 1 4 4 4 2 2 2 2 2
## SE2-8 SE3-1 SE3-10 SE3-2-1 SE3-2-2 SE3-4-1 SE3-5 SE3-8-1 SE3-9 UA1-1
## 2 4 4 4 2 4 4 2 4 1
## UA1-10 UA1-12 UA1-16 UA1-3 UA1-7 UA2-11 UA2-14 UA2-16 UA2-19 UA2-2
## 1 1 1 1 1 1 1 1 1 1
## UA2-5 UA2-7
## 1 1
## Levels: 1 2 3 4
var_exp <- dapc_res4$eig / sum(dapc_res4$eig) * 100
var_exp
## [1] 65.88197 20.77813 13.33990
cluster_colors <- c("#009988", "#5E4FA2", "#EE7733", "#0077BB" )
scatter.dapc(dapc_res4,
bg = "white",
col = cluster_colors,
pch = 19,
cstar = 0, solid = 0.75,
scree.da = TRUE,
scree.pca =TRUE,
posi.pca = "topright",
posi.da = "bottomright",
label.inds = list(air = 0.8, pch = NA, col = "gray40", cex = 0.6, pos = 3))
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "pos" is not a graphical
## parameter
