knitr::opts_chunk$set(warning = FALSE, message = FALSE)
library(vcfR)
##
## ***** *** vcfR *** *****
## This is vcfR 1.15.0
## browseVignettes('vcfR') # Documentation
## citation('vcfR') # Citation
## ***** ***** ***** *****
library(dartR)
## Loading required package: adegenet
## Loading required package: ade4
##
## /// adegenet 2.1.11 is loaded ////////////
##
## > overview: '?adegenet'
## > tutorials/doc/questions: 'adegenetWeb()'
## > bug reports/feature requests: adegenetIssues()
## Loading required package: ggplot2
## Loading required package: dplyr
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
## Loading required package: dartR.data
## **** Welcome to dartR.data [Version 1.2.2 ] ****
## Registered S3 method overwritten by 'pegas':
## method from
## print.amova ade4
## Registered S3 method overwritten by 'genetics':
## method from
## [.haplotype pegas
## **** Welcome to the legacy dartR [Version 2.9.9.5 ] ****
## ##############################
## Be aware this the legacy version of dartR, which is no longer developed. Please consider using the new version of dartR, called 'dartRverse'. As before you can find the latest version on CRAN. dartRverse has the same functionality as dartR, but is more up-to-date and maintained. To reduce dependencies, and improve the overall performance, dartRverse consists of several packacages, that are maintained seperately to reduce dependencies . To install dartRverse please download the package dartRverse and follow the instructions therein. For more information please visit: https://github.com/green-striped-gecko/dartRverse.
## ##############################
## Be aware that owing to CRAN requirements and compatibility reasons not all functions of the package may run after the basic installation, as some packages could still be missing. Hence for a most enjoyable experience we recommend to run the function
## gl.install.vanilla.dartR()
## This installs all missing and required packages for your version of dartR. In case something fails during installation please refer to this tutorial: https://github.com/green-striped-gecko/dartR/wiki/Installation-tutorial.
##
## For information how to cite dartR, please use:
## citation('dartR')
## Global verbosity is set to: 2
##
## **** Have fun using dartR! ****
##
## Attaching package: 'dartR'
## The following objects are masked from 'package:dartR.data':
##
## bandicoot.gl, possums.gl, testset.gl, testset.gs
library(ggplot2)
library(ggrepel)
sessionInfo()
## R version 4.5.0 (2025-04-11 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=English_United Kingdom.utf8
## [2] LC_CTYPE=English_United Kingdom.utf8
## [3] LC_MONETARY=English_United Kingdom.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United Kingdom.utf8
##
## time zone: Europe/Berlin
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggrepel_0.9.6 dartR_2.9.9.5 dartR.data_1.2.2 dplyr_1.2.1
## [5] ggplot2_4.0.2 adegenet_2.1.11 ade4_1.7-23 vcfR_1.15.0
##
## loaded via a namespace (and not attached):
## [1] Rdpack_2.6.6 gdsfmt_1.44.1 gridExtra_2.3
## [4] permute_0.9-10 rlang_1.2.0 magrittr_2.0.3
## [7] otel_0.2.0 compiler_4.5.0 mgcv_1.9-4
## [10] gdata_3.0.1 png_0.1-8 vctrs_0.7.3
## [13] maps_3.4.3 reshape2_1.4.5 combinat_0.0-9
## [16] stringr_1.6.0 pkgconfig_2.0.3 crayon_1.5.3
## [19] SNPRelate_1.42.0 fastmap_1.2.0 promises_1.5.0
## [22] rmarkdown_2.30 purrr_1.2.1 xfun_0.61
## [25] cachem_1.1.0 seqinr_4.2-36 jsonlite_2.0.0
## [28] later_1.4.8 terra_1.8-93 parallel_4.5.0
## [31] cluster_2.1.8.2 R6_2.6.1 gap.datasets_0.0.6
## [34] bslib_0.10.0 stringi_1.8.7 RColorBrewer_1.1-3
## [37] GGally_2.4.0 jquerylib_0.1.4 gdistance_1.6.5
## [40] iterators_1.0.14 Rcpp_1.1.0 knitr_1.51
## [43] fields_17.3 genetics_1.3.8.1.3 R.utils_2.13.0
## [46] httpuv_1.6.17 Matrix_1.7-4 splines_4.5.0
## [49] igraph_2.2.2 tidyselect_1.2.1 rstudioapi_0.18.0
## [52] dichromat_2.0-0.1 yaml_2.3.10 vegan_2.7-2
## [55] doParallel_1.0.17 codetools_0.2-20 lattice_0.22-9
## [58] tibble_3.3.1 plyr_1.8.9 shiny_1.12.1
## [61] withr_3.0.2 S7_0.2.1 evaluate_1.0.5
## [64] pinfsc50_1.3.0 mmod_1.3.3 ggstats_0.12.0
## [67] pillar_1.11.1 gap_1.15.2 foreach_1.5.2
## [70] generics_0.1.4 sp_2.2-1 scales_1.4.0
## [73] calibrate_1.7.7 gtools_3.9.5 xtable_1.8-4
## [76] glue_1.8.0 tools_4.5.0 data.table_1.18.2.1
## [79] gsubfn_0.7 mvtnorm_1.3-3 dotCall64_1.2
## [82] grid_4.5.0 tidyr_1.3.2 ape_5.8-1
## [85] rbibutils_2.4.1 nlme_3.1-168 patchwork_1.3.2
## [88] raster_3.6-32 StAMPP_1.6.3 RgoogleMaps_1.5.3
## [91] proto_1.0.0 cli_3.6.5 spam_2.11-3
## [94] viridisLite_0.4.3 gtable_0.3.6 R.methodsS3_1.8.2
## [97] sass_0.4.10 digest_0.6.37 farver_2.1.2
## [100] htmltools_0.5.8.1 R.oo_1.27.1 lifecycle_1.0.5
## [103] dismo_1.3-16 pegas_1.4 PopGenReport_3.1.3
## [106] mime_0.13 MASS_7.3-65
Data preparation
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
dat1
## ***** Object of Class vcfR *****
## 317 samples
## 5870 CHROMs
## 25,099 variants
## Object size: 71.5 Mb
## 0 percent missing data
## ***** ***** *****
## To remove clones detected before 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)
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
PCoA (Axis 1 vs Axis 2)
pc1 <- gl.pcoa(dat1.pre_gl, nfactors = 3, correction = "lingoes", FALSE)
## Starting gl.pcoa
## Processing genlight object with SNP data
## Performing a PCA, individuals as entities, loci as attributes, SNP genotype as state
## Completed: gl.pcoa
eig_vals <- pc1$eig
total_var <- sum(eig_vals)
axis1 <- 100 * eig_vals[1] / total_var
axis2 <- 100 * eig_vals[2] / total_var
axis3 <- 100 * eig_vals[3] / total_var
out1 <- as.data.frame(pc1$scores)
out1$SampleID <- rownames(pc1$scores)
rownames(out1) <- indNames(dat1.pre_gl)
out1$Country<-substr(out1$SampleID, 1, 2)
palette_10_light_bg <- c(
"#009988","#33BBEE","#EE7733","#CC3311", "#CC6677","#009988","#33BBEE","#EE7733","#CC3311","#CC6677")
ggplot(out1, aes(x = PC1, y = PC2, fill = Country, shape = Country)) +
stat_ellipse(level = 0.95, aes(fill = Country, color = Country), alpha =0.15, geom = "polygon") +
geom_point(size = 2.5, color = "black", alpha = 0.8) +
geom_text_repel(aes(label = SampleID, color = Country), size = 2, force = 1, max.overlaps = 10) +
labs(x = paste0("Principal Coordinate Axis 1 (", round(axis1, 2), "%)"),
y = paste0("Principal Coordinate Axis 2 (", round(axis2, 2), "%)")) +
theme_classic()+
scale_shape_manual(values = c(21,21,21,21,21,25,25,25,25,25)) +
scale_fill_manual(values = palette_10_light_bg)+
scale_color_manual(values = palette_10_light_bg)
PCoA (Axis 1 vs Axis 3)
ggplot(out1, aes(x = PC1, y = PC3, fill = Country, shape = Country)) +
stat_ellipse(level = 0.95, aes(fill = Country, color = Country), alpha =0.15, geom = "polygon") +
geom_point(size = 2.5, color = "black", alpha = 0.8) +
geom_text_repel(aes(label = SampleID, color = Country), size = 2, force = 1, max.overlaps = 10) +
labs(x = paste0("Principal Coordinate Axis 1 (", round(axis1, 2), "%)"),
y = paste0("Principal Coordinate Axis 3 (", round(axis3, 2), "%)")) +
theme_classic()+
scale_shape_manual(values = c(21,21,21,21,21,25,25,25,25,25)) +
scale_fill_manual(values = palette_10_light_bg)+
scale_color_manual(values = palette_10_light_bg)
