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] maps_3.4.3 scatterpie_0.2.6 rnaturalearthdata_1.0.0
## [4] rnaturalearth_1.2.0 ggforce_0.5.0 sf_1.0-24
## [7] ggrepel_0.9.6 reshape2_1.4.5 lubridate_1.9.5
## [10] forcats_1.0.1 stringr_1.6.0 dplyr_1.2.1
## [13] purrr_1.2.1 readr_2.1.6 tidyr_1.3.2
## [16] tibble_3.3.1 ggplot2_4.0.2 tidyverse_2.0.0
##
## loaded via a namespace (and not attached):
## [1] gtable_0.3.6 xfun_0.61 bslib_0.10.0 tzdb_0.5.0
## [5] vctrs_0.7.3 tools_4.5.0 generics_0.1.4 yulab.utils_0.2.4
## [9] proxy_0.4-29 pkgconfig_2.0.3 KernSmooth_2.23-26 RColorBrewer_1.1-3
## [13] S7_0.2.1 lifecycle_1.0.5 compiler_4.5.0 farver_2.1.2
## [17] ggfun_0.2.1 htmltools_0.5.8.1 class_7.3-23 sass_0.4.10
## [21] yaml_2.3.10 pillar_1.11.1 jquerylib_0.1.4 MASS_7.3-65
## [25] classInt_0.4-11 cachem_1.1.0 tidyselect_1.2.1 digest_0.6.37
## [29] stringi_1.8.7 polyclip_1.10-7 fastmap_1.2.0 grid_4.5.0
## [33] cli_3.6.5 magrittr_2.0.3 dichromat_2.0-0.1 e1071_1.7-17
## [37] withr_3.0.2 scales_1.4.0 rappdirs_0.3.4 timechange_0.4.0
## [41] rmarkdown_2.30 otel_0.2.0 hms_1.1.4 evaluate_1.0.5
## [45] knitr_1.51 rlang_1.2.0 Rcpp_1.1.0 glue_1.8.0
## [49] DBI_1.3.0 tweenr_2.0.3 rstudioapi_0.18.0 jsonlite_2.0.0
## [53] R6_2.6.1 plyr_1.8.9 fs_2.1.0 units_1.0-0
Determinung number of K
cv <- read.table("admixture/cv_errors.txt", header=FALSE)
colnames(cv) <- c("K", "CV")
cv_df <- data.frame(K = cv$K,CV = cv$CV)
ggplot(cv_df, aes(x = K, y = CV)) +
geom_line(linewidth = 0.8, color = "#009988") +
geom_point(size = 3, shape = 21, fill = "white", colour = "#009988", stroke = 1) +
scale_x_continuous(breaks = cv_df$K) +
labs(x = "K (number of clusters)",
y = "CV error") +
theme_minimal(base_size = 14) +
theme(
panel.grid.minor = element_blank(),
plot.title = element_text(hjust = 0.5, face = "bold"),
axis.title = element_text(face = "bold"))
qpath <- "admixture/"
famfile <- "admixture/md6ct95ms200_numeric.fam"
samples <- read.table(famfile, stringsAsFactors = FALSE)
sample_ids <- samples$V2 # V2 = individual IDs
readQ <- function(K) {
qfile <- paste0(qpath, "md6ct95ms200_numeric.", K, ".Q")
q <- read.table(qfile, stringsAsFactors = FALSE)
colnames(q) <- paste0("Pop", 1:K)
q$SampleID <- sample_ids
q$K <- K
q }
Ks <- c(2,3,4)
Qlist <- lapply(Ks, readQ)
Qdata <- bind_rows(Qlist)
# Reshape to long format for ggplot
Qlong <- melt(Qdata, id.vars = c("SampleID", "K"),
variable.name = "Cluster", value.name = "Ancestry")
length(unique(Qlong$SampleID))
## [1] 311
# To sort samples consistently (optional: here by SampleID)
Qlong$SampleID <- factor(Qlong$SampleID, levels = unique(Qlong$SampleID))
palette <- c("#E69F00","#009E73","#0072B2" ,"#F0E442")
ggplot(Qlong, aes(x = SampleID, y = Ancestry, fill = Cluster)) +
geom_bar(stat = "identity", width = 1) +
facet_wrap(~K, ncol = 1, scales = "free_x") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 95, size = 4)) +
labs(title = "ADMIXTURE ancestry proportions",
x = "Individuals",
y = "Ancestry proportion")+
scale_fill_manual(values = palette)
## Warning: Removed 933 rows containing missing values or values outside the scale range
## (`geom_bar()`).
Prepare coordinates
Coordinates_pops <- read_delim("Coordinates_pops.txt", delim = "\t", escape_double = FALSE, trim_ws = TRUE)
## Rows: 310 Columns: 3
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (1): SampleID
## dbl (2): Latitude, Longitude
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Coordinates_pops$Population <- sub("-.*", "", Coordinates_pops$SampleID)
Qwide <- Qdata %>% left_join(Coordinates_pops, by = c("SampleID"))
Qlong <- Qdata %>%
pivot_longer(cols = matches("^Pop\\d+$"),
names_to = "Cluster",
values_to = "Ancestry")
Qlong <- Qlong %>%
left_join(Coordinates_pops, by = "SampleID")
Qlong <- Qlong %>%
filter(!is.na(Longitude))
pop_means_long <- Qlong %>%
group_by(Population, K, Cluster) %>%
summarise(Ancestry = mean(Ancestry, na.rm = TRUE))
## `summarise()` has regrouped the output.
## ℹ Summaries were computed grouped by Population, K, and Cluster.
## ℹ Output is grouped by Population and K.
## ℹ Use `summarise(.groups = "drop_last")` to silence this message.
## ℹ Use `summarise(.by = c(Population, K, Cluster))` for per-operation grouping
## (`?dplyr::dplyr_by`) instead.
pop_coords <- Coordinates_pops %>%
group_by(Population) %>%
summarise(Longitude = mean(Longitude, na.rm = TRUE),
Latitude = mean(Latitude, na.rm = TRUE))
pop_means_wide <- pop_means_long %>%
pivot_wider(names_from = Cluster, values_from = Ancestry) %>%
left_join(pop_coords, by = "Population")
head(pop_means_wide)
## # A tibble: 6 × 8
## # Groups: Population, K [6]
## Population K Pop1 Pop2 Pop3 Pop4 Longitude Latitude
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 DE1 2 1.000 0.00001 NaN NaN 7.89 51.8
## 2 DE1 3 0.00001 0.00001 1.000 NaN 7.89 51.8
## 3 DE1 4 0.00001 1.000 0.00001 0.00001 7.89 51.8
## 4 DE10 2 0.577 0.423 NaN NaN 9.74 48.5
## 5 DE10 3 0.190 0.289 0.520 NaN 9.74 48.5
## 6 DE10 4 0.112 0.414 0.0286 0.445 9.74 48.5
Admixture on the map K = 2
The pie locations are slightly jittered for visualisation
purpose
Kplot <- 2
cols <- paste0("Pop", 1:Kplot) # required columns for that K
# Keep only relevant columns and K = 2
plot_data <- pop_means_wide %>%
dplyr::filter(K == Kplot) %>%
dplyr::select(any_of(c("Population","Pop1", "Pop2","K", "Longitude", "Latitude")), all_of(cols))
# Reproject both map and data
europe <- ne_countries(continent = "Europe", scale = "medium", returnclass = "sf")
europe_proj <- st_transform(europe, 3035)
plot_data_sf <- st_as_sf(plot_data, coords = c("Longitude", "Latitude"), crs = 4326)
plot_data_proj <- st_transform(plot_data_sf, 3035)
# Extract projected coordinates as plain data frame
plot_data_proj_df <- cbind(st_drop_geometry(plot_data_proj), st_coordinates(plot_data_proj))
# Make sure Pop columns are numeric
plot_data_proj_df <- plot_data_proj_df %>%
mutate(across(all_of(cols), as.numeric))
palette <- c("#E69F00","#009E73")
ggplot(data = europe_proj) +
geom_sf(fill = "grey80", color = "grey60", linewidth = 0.8) +
geom_scatterpie(
data = plot_data_proj_df,
aes(x = X, y = Y, r = 40000),
cols = cols) +
geom_text_repel(
data = plot_data_proj_df,
aes(x = X, y = Y, label = Population),
size = 3, fontface = "bold",
segment.color = NA,
max.overlaps = 20, box.padding = 0.5, # push label away from itself
point.padding = 0.2, # push label away from pie center
min.segment.length = 0) +
coord_sf(xlim = c(3900000, 5700000), ylim = c(2500000, 4500000), expand = FALSE) +
scale_fill_manual(values = palette[1:length(cols)]) +
theme_minimal()
Admixture K = 3 on the map
The pie locations are slightly jittered for visualisation
purpose
Kplot <- 3
cols <- paste0("Pop", 1:Kplot) # required columns for that K
# Keep only relevant columns and K = 3
plot_data <- pop_means_wide %>%
dplyr::filter(K == Kplot) %>%
dplyr::select(any_of(c("Population","Pop1", "Pop2","Pop3","K", "Longitude", "Latitude")), all_of(cols))
# Reproject both map and data
europe <- ne_countries(continent = "Europe", scale = "medium", returnclass = "sf")
europe_proj <- st_transform(europe, 3035)
plot_data_sf <- st_as_sf(plot_data, coords = c("Longitude", "Latitude"), crs = 4326)
plot_data_proj <- st_transform(plot_data_sf, 3035)
# Extract projected coordinates as plain data frame
plot_data_proj_df <- cbind(st_drop_geometry(plot_data_proj), st_coordinates(plot_data_proj))
# Make sure Pop columns are numeric
plot_data_proj_df <- plot_data_proj_df %>%
mutate(across(all_of(cols), as.numeric))
palette <- c("#0072B2","#009E73","#E69F00")
ggplot(data = europe_proj) +
geom_sf(fill = "grey80", color = "grey60", linewidth = 0.8) +
geom_scatterpie(
data = plot_data_proj_df,
aes(x = X, y = Y, r = 40000),
cols = cols) +
geom_text_repel(
data = plot_data_proj_df,
aes(x = X, y = Y, label = Population),
size = 3, fontface = "bold",
segment.color = NA,
max.overlaps = 20, box.padding = 0.5, # push label away from itself
point.padding = 0.2, # push label away from pie center
min.segment.length = 0) +
coord_sf(xlim = c(3900000, 5700000), ylim = c(2500000, 4500000), expand = FALSE) +
scale_fill_manual(values = palette[1:length(cols)]) +
theme_minimal()
Admixture K = 4 on the map
The pie locations are slightly jittered for visualisation
purpose
Kplot <- 4
cols <- paste0("Pop", 1:Kplot) # required columns for that K
# Keep only relevant columns and K = 4
plot_data <- pop_means_wide %>%
dplyr::filter(K == Kplot) %>%
dplyr::select(any_of(c("Population","Pop1", "Pop2","Pop3","Pop4","K", "Longitude", "Latitude")), all_of(cols))
# Reproject both map and data
europe <- ne_countries(continent = "Europe", scale = "medium", returnclass = "sf")
europe_proj <- st_transform(europe, 3035)
plot_data_sf <- st_as_sf(plot_data, coords = c("Longitude", "Latitude"), crs = 4326)
plot_data_proj <- st_transform(plot_data_sf, 3035)
# Extract projected coordinates as plain data frame
plot_data_proj_df <- cbind(st_drop_geometry(plot_data_proj), st_coordinates(plot_data_proj))
# Make sure Pop columns are numeric
plot_data_proj_df <- plot_data_proj_df %>%
mutate(across(all_of(cols), as.numeric))
palette <- c("#009E73","#E69F00","#0072B2" ,"#5E4FA2")
ggplot(data = europe_proj) +
geom_sf(fill = "grey80", color = "grey60", linewidth = 0.8) +
geom_scatterpie(
data = plot_data_proj_df,
aes(x = X, y = Y, r = 40000),
cols = cols) +
geom_text_repel(
data = plot_data_proj_df,
aes(x = X, y = Y, label = Population),
size = 3, fontface = "bold",
segment.color = NA,
max.overlaps = 20, box.padding = 0.5, # push label away from itself
point.padding = 0.2, # push label away from pie center
min.segment.length = 0) +
coord_sf(xlim = c(3900000, 5700000), ylim = c(2500000, 4500000), expand = FALSE) +
scale_fill_manual(values = palette[1:length(cols)]) +
theme_minimal()
