Population label key: The code labels used in the analysis below correspond to the following manuscript labels: AUT = AUTO, MAN = NON-AUTO, NEW = NON-AUTO-FIELD.

1. R libraries and software

library(tidyverse)
library(here)
library(qvalue)
library(ggvenn)
library(ggrepel)
library(pcadapt)

2. Helper functions

# Load custom ggplot theme
source(here("notebooks", "helpers", "my_theme2.R"))

# Chromosome color palette
chrom_colors <- c("#CCF6D6", "#F6E1CC", "#CCD8F6")

# Format genomic positions as Mb
mb_format <- function(x) {
  paste0(format(x / 1e6, big.mark = "", scientific = FALSE), "Mb")
}

# Run pcadapt, build results data frame, mark BH outliers, save SNP list
run_pcadapt <- function(bed_file, bim_file, K, out_snp_file) {
  x <- read.pcadapt(here("output", "outflank", bed_file), type = "bed")
  res <- pcadapt(x, method = "mahalanobis", min.maf = 0.1,
                 LD.clumping = NULL, tol = 1e-04, K = K)
  snps <- read_delim(
    here("output", "outflank", bim_file),
    col_names = FALSE, show_col_types = FALSE, col_types = "ccidcc"
  )
  colnames(snps) <- c("Scaffold", "SNP", "Cm", "Position", "Allele1", "Allele2")
  df <- data.frame(
    SNP        = snps$SNP,
    Chromosome = snps$Scaffold,
    Position   = snps$Position,
    pvalues    = res$pvalues,
    stat       = res$stat
  ) |> drop_na()
  padj <- p.adjust(df$pvalues, method = "BH")
  df$highlight <- padj < 0.05
  outlier_snps <- df$SNP[df$highlight]
  write.table(outlier_snps,
    file = here("output", "pcadapt", out_snp_file),
    row.names = FALSE, quote = FALSE, col.names = FALSE, sep = "\n"
  )
  list(res = res, df = df)
}

# Manhattan plot for pcadapt results
plot_manhattan <- function(df, title_label, save_path) {
  color_vector <- chrom_colors
  names(color_vector) <- unique(df$Chromosome)
  min_pval_snps <- df |>
    filter(highlight) |>
    group_by(Chromosome) |>
    slice(which.min(pvalues))
  p <- ggplot(df, aes(x = Position, y = -log10(pvalues))) +
    geom_point(aes(color = Chromosome),
               data = subset(df, !highlight), size = .5) +
    geom_point(color = "magenta",
               data = subset(df, highlight), size = .5) +
    scale_color_manual(values = color_vector, guide = "none") +
    labs(
      title   = title_label,
      x       = "Position",
      y       = "-log10(p-value)",
      caption = "Magenta: significant after Benjamini-Hochberg adjustment (alpha = 0.05)."
    ) +
    facet_wrap(~Chromosome, scales = "free_x") +
    scale_x_continuous(labels = mb_format) +
    my_theme() +
    theme(
      panel.spacing  = unit(.2, "lines"),
      plot.margin    = unit(c(1, 1, 2, 2), "lines"),
      plot.caption   = element_text(face = "italic", color = "#574E4E")
    ) +
    geom_text_repel(
      data = min_pval_snps, aes(label = SNP),
      size = 3, nudge_y = 0.01, segment.color = NA, max.overlaps = Inf
    )
  ggsave(here("output", "pcadapt", "figures", save_path),
         width = 8, height = 5, units = "in")
  p
}

# Outlier method comparison: q-values, BH, Bonferroni
compare_thresholds <- function(res) {
  qval   <- qvalue(res$pvalues)$qvalues
  cat("q-values (alpha=0.05):    ", sum(qval < 0.05), "\n")
  cat("BH adjusted (alpha=0.05): ",
      sum(p.adjust(res$pvalues, method = "BH") < 0.05), "\n")
  cat("Bonferroni (alpha=0.05):  ",
      sum(p.adjust(res$pvalues, method = "bonferroni") < 0.05), "\n")
}

# Venn diagram and shared-SNP file for outFlank vs pcadapt
venn_outflank_pcadapt <- function(outflank_file, pcadapt_file, venn_pdf, common_txt,
                                   fill_colors = c("steelblue", "darkorange")) {
  of_snps <- read.table(here("output", "outflank", outflank_file),
                        stringsAsFactors = FALSE) |> drop_na()
  pc_snps <- read.table(here("output", "pcadapt", pcadapt_file),
                        stringsAsFactors = FALSE)
  lst <- list(outFlank = of_snps$V1, pcadapt = pc_snps$V1)
  vd  <- ggvenn(lst, fill_color = fill_colors)
  print(vd)
  common <- Reduce(intersect, lst)
  write.table(common,
    file = here("output", "pcadapt", common_txt),
    row.names = FALSE, col.names = FALSE, quote = FALSE
  )
  ggsave(here("output", "pcadapt", "figures", venn_pdf),
         vd, height = 5, width = 5, dpi = 300)
  common
}

# SNP cluster table (10 Mb windows) as flextable + Word doc
cluster_table <- function(filtered_df, docx_file) {
  cluster_snps_df <- filtered_df |>
    mutate(Window_Start = floor(Position / 1e7) * 1e7) |>
    group_by(Chromosome, Window_Start) |>
    summarise(SNPs = list(SNP), .groups = "drop") |>
    mutate(SNPs = sapply(SNPs, paste, collapse = ", "))
  knitr::kable(cluster_snps_df)
}

3. pcadapt with all 3 populations (NON-AUTO, NON-AUTO-FIELD and AUTO)

pcadapt performs principal component analysis and tests for selection outliers based on correlations between genetic variation and the first K principal components.

K = 2 was selected based on scree plot inspection.

x_all3 <- read.pcadapt(here("output", "outflank", "autogenous.bed"), type = "bed")
plot(pcadapt(x_all3, K = 30), option = "screeplot")

out_all3 <- run_pcadapt("autogenous.bed", "autogenous.bim",
                         K = 2, out_snp_file = "SNPs_pcadapt.txt")
pcadapt_res <- out_all3$res
pcadapt_res_df <- out_all3$df
cat("BH outliers:", sum(pcadapt_res_df$highlight), "\n")
## BH outliers: 37
plot_manhattan(pcadapt_res_df,
               title_label = "pcadapt: NON-AUTO + NON-AUTO-FIELD + AUTO",
               save_path   = "pcadapt.pdf")

plot(pcadapt_res, option = "qqplot")

plot(pcadapt_res, option = "scores")

hist(pcadapt_res$pvalues, xlab = "p-values", main = NULL, breaks = 50, col = "pink")

compare_thresholds(pcadapt_res)
## q-values (alpha=0.05):     NA 
## BH adjusted (alpha=0.05):  NA 
## Bonferroni (alpha=0.05):   NA
par(mfrow = c(2, 2))
for (i in 1:2)
  plot(pcadapt_res$loadings[, i], pch = 19, cex = .3, ylab = paste0("Loadings PC", i))

3.1 Venn diagram outFlank vs pcadapt

Note: This section was not used in the final version of the manuscript.

common_SNPs <- venn_outflank_pcadapt(
  outflank_file = "SNPs_outFlank.txt",
  pcadapt_file  = "SNPs_pcadapt.txt",
  venn_pdf      = "significant_snps.pdf",
  common_txt    = "common_SNPs_pcadapt_outflank.txt"
)

filtered_df <- pcadapt_res_df |> filter(SNP %in% common_SNPs)
color_vector <- chrom_colors
names(color_vector) <- unique(pcadapt_res_df$Chromosome)
min_pval_snps <- pcadapt_res_df |> filter(highlight) |>
  group_by(Chromosome) |> slice(which.min(pvalues))

ggplot(pcadapt_res_df, aes(x = Position, y = -log10(pvalues))) +
  geom_point(aes(color = Chromosome), data = subset(pcadapt_res_df, !highlight), size = .3) +
  geom_point(color = "magenta", data = subset(filtered_df, highlight), size = .3) +
  scale_color_manual(values = color_vector, guide = "none") +
  labs(x = "Position", y = "-log10(p-value)",
       caption = "Magenta: BH-adjusted (alpha=0.05) outliers shared between pcadapt and outFlank.") +
  facet_wrap(~Chromosome, scales = "free_x") +
  scale_x_continuous(labels = mb_format) +
  my_theme() +
  theme(panel.spacing = unit(.2, "lines"),
        plot.margin = unit(c(1, 1, 2, 2), "lines"),
        plot.caption = element_text(face = "italic")) +
  geom_text_repel(data = min_pval_snps, aes(label = SNP),
                  size = 3, nudge_y = 0.01, segment.color = NA, max.overlaps = Inf)

ggsave(here("output", "pcadapt", "figures", "pcadapt_outflank.pdf"),
       width = 8, height = 5, units = "in")
cluster_table(filtered_df, "snps_pcadap_outflank.docx")
Chromosome Window_Start SNPs
1 1.4e+08 AX-583515734
2 3.0e+07 AX-584444558
2 1.8e+08 AX-585196879
2 3.2e+08 AX-584907866
2 3.8e+08 AX-579548089
2 3.9e+08 AX-579560686, AX-579584369, AX-579583040, AX-579584883
2 4.0e+08 AX-579604213, AX-579607140, AX-579619995
2 4.1e+08 AX-579630465, AX-579632196, AX-579632983, AX-579632927, AX-579661979, AX-579667335
2 5.1e+08 AX-579987164
2 5.6e+08 AX-580195612
3 0.0e+00 AX-580344997
3 1.8e+08 AX-581275893
3 1.9e+08 AX-581302901
3 2.1e+08 AX-581437212, AX-581438408, AX-581442470, AX-581461793, AX-581462648, AX-581467653
3 2.2e+08 AX-581504582, AX-581512892
3 2.3e+08 AX-581527485, AX-581534474

4. pcadapt: NON-AUTO vs AUTO (pairwise)

pcadapt performs principal component analysis and tests for selection outliers based on correlations between genetic variation and the first K principal components.

rm(list = setdiff(ls(), c("run_pcadapt", "plot_manhattan", "compare_thresholds",
                           "venn_outflank_pcadapt", "cluster_table",
                           "chrom_colors", "mb_format", "my_theme")))
gc()
##           used  (Mb) gc trigger  (Mb) max used  (Mb)
## Ncells 2508471 134.0    4241018 226.5  4189766 223.8
## Vcells 4782954  36.5   13285406 101.4 13285406 101.4
x_ma <- read.pcadapt(here("output", "outflank", "man_aut.bed"), type = "bed")
plot(pcadapt(x_ma, K = 30), option = "screeplot")

out_ma <- run_pcadapt("man_aut.bed", "man_aut.bim",
                       K = 2, out_snp_file = "man_aut_SNPs_pcadapt.txt")
pcadapt_res <- out_ma$res
pcadapt_res_df <- out_ma$df
cat("BH outliers:", sum(pcadapt_res_df$highlight), "\n")
## BH outliers: 115
plot_manhattan(pcadapt_res_df,
               title_label = "pcadapt: NON-AUTO vs AUTO",
               save_path   = "man_aut_pcadapt.pdf")

plot(pcadapt_res, option = "qqplot")

plot(pcadapt_res, option = "scores")

hist(pcadapt_res$pvalues, xlab = "p-values", main = NULL, breaks = 50, col = "pink")

compare_thresholds(pcadapt_res)
## q-values (alpha=0.05):     NA 
## BH adjusted (alpha=0.05):  NA 
## Bonferroni (alpha=0.05):   NA
par(mfrow = c(2, 2))
for (i in 1:2)
  plot(pcadapt_res$loadings[, i], pch = 19, cex = .3, ylab = paste0("Loadings PC", i))

4.1 Venn diagram outFlank vs pcadapt

Note: This section was not used in the final version of the manuscript.

common_SNPs <- venn_outflank_pcadapt(
  outflank_file = "man_aut_SNPs_outFlank.txt",
  pcadapt_file  = "man_aut_SNPs_pcadapt.txt",
  venn_pdf      = "man_aut_significant_snps.pdf",
  common_txt    = "man_aut_common_SNPs_pcadapt_outflank.txt"
)

filtered_df <- pcadapt_res_df |> filter(SNP %in% common_SNPs)
color_vector <- chrom_colors
names(color_vector) <- unique(pcadapt_res_df$Chromosome)
min_pval_snps <- pcadapt_res_df |> filter(highlight) |>
  group_by(Chromosome) |> slice(which.min(pvalues))

ggplot(pcadapt_res_df, aes(x = Position, y = -log10(pvalues))) +
  geom_point(aes(color = Chromosome), data = subset(pcadapt_res_df, !highlight), size = .3) +
  geom_point(color = "magenta", data = subset(filtered_df, highlight), size = .3) +
  scale_color_manual(values = color_vector, guide = "none") +
  labs(x = "Position", y = "-log10(p-value)",
       caption = "Magenta: BH-adjusted (alpha=0.05) outliers shared between pcadapt and outFlank.") +
  facet_wrap(~Chromosome, scales = "free_x") +
  scale_x_continuous(labels = mb_format) +
  my_theme() +
  theme(panel.spacing = unit(.2, "lines"),
        plot.margin = unit(c(1, 1, 2, 2), "lines"),
        plot.caption = element_text(face = "italic")) +
  geom_text_repel(data = min_pval_snps, aes(label = SNP),
                  size = 3, nudge_y = 0.01, segment.color = NA, max.overlaps = Inf)

ggsave(here("output", "pcadapt", "figures", "man_aut_pcadapt_outflank.pdf"),
       width = 8, height = 5, units = "in")
cluster_table(filtered_df, "man_aut_snps_pcadap_outflank.docx")
Chromosome Window_Start SNPs
1 0.0e+00 AX-583054970
1 1.0e+07 AX-583095890
1 9.0e+07 AX-583324654
1 1.0e+08 AX-583355868
1 1.2e+08 AX-583423493, AX-583426050
1 1.3e+08 AX-583467551, AX-583504302
1 1.4e+08 AX-583514148, AX-583515734, AX-583518994, AX-583517999
1 1.7e+08 AX-583631472
1 2.7e+08 AX-583889791, AX-583905588, AX-583924569
1 3.4e+08 AX-584133664
2 2.0e+07 AX-584419201, AX-584427590
2 3.0e+07 AX-584444558, AX-584453518
2 4.0e+07 AX-584473131
2 5.0e+07 AX-585010439, AX-584502184, AX-585023213
2 7.0e+07 AX-585042322
2 8.0e+07 AX-584546335
2 1.1e+08 AX-585101539
2 1.7e+08 AX-585182739
2 1.8e+08 AX-585190366, AX-585196879
2 2.0e+08 AX-582494788, AX-582519052, AX-582532563
2 2.6e+08 AX-584765066
2 3.2e+08 AX-585400879, AX-585413917
2 3.5e+08 AX-584951049
2 3.6e+08 AX-579474142, AX-579474650
2 3.7e+08 AX-579509845
2 3.8e+08 AX-579548089
2 3.9e+08 AX-579560686, AX-579580051, AX-579584369, AX-579583040, AX-579584883
2 4.0e+08 AX-579604213, AX-579603553, AX-579604110, AX-579607140, AX-579606795, AX-579619995, AX-579620627
2 4.1e+08 AX-579630462, AX-579630465, AX-579632196, AX-579632927, AX-579636106, AX-579638540, AX-579640459, AX-579646083, AX-579661979, AX-579662371
2 4.5e+08 AX-579768613, AX-579778684
2 4.6e+08 AX-579808128, AX-579816168
2 4.9e+08 AX-579909895
2 5.1e+08 AX-579955661, AX-579982534
2 5.2e+08 AX-580001880, AX-580012078, AX-580016916, AX-580026657, AX-580027537, AX-580042929
2 5.4e+08 AX-580099245
2 5.5e+08 AX-580165355
2 5.6e+08 AX-580195612
2 5.8e+08 AX-580318010, AX-580321888
3 0.0e+00 AX-580342672
3 2.0e+07 AX-580425125
3 7.0e+07 AX-580620056
3 8.0e+07 AX-580673467
3 9.0e+07 AX-580719743
3 1.3e+08 AX-580926966, AX-580930203
3 1.8e+08 AX-581275893
3 1.9e+08 AX-581298415, AX-581302704
3 2.0e+08 AX-581399096
3 2.1e+08 AX-581408938, AX-581437212, AX-581438408, AX-581442470, AX-581458697, AX-581461793, AX-581462648, AX-581467653, AX-581470736
3 2.2e+08 AX-581485809, AX-581488185
3 2.3e+08 AX-581572119
3 3.5e+08 AX-582085824
3 4.4e+08 AX-582853985
3 4.5e+08 AX-582884221
3 4.8e+08 AX-583025646

5. pcadapt: NON-AUTO-FIELD vs AUTO (pairwise)

pcadapt performs principal component analysis and tests for selection outliers based on correlations between genetic variation and the first K principal components.

rm(list = setdiff(ls(), c("run_pcadapt", "plot_manhattan", "compare_thresholds",
                           "venn_outflank_pcadapt", "cluster_table",
                           "chrom_colors", "mb_format", "my_theme")))
gc()
##           used  (Mb) gc trigger  (Mb) max used  (Mb)
## Ncells 2507940 134.0    4241018 226.5  4241018 226.5
## Vcells 4762886  36.4   16024048 122.3 16024048 122.3
x_na <- read.pcadapt(here("output", "outflank", "new_aut.bed"), type = "bed")
plot(pcadapt(x_na, K = 30), option = "screeplot")

out_na <- run_pcadapt("new_aut.bed", "new_aut.bim",
                       K = 2, out_snp_file = "new_aut_SNPs_pcadapt.txt")
pcadapt_res <- out_na$res
pcadapt_res_df <- out_na$df
cat("BH outliers:", sum(pcadapt_res_df$highlight), "\n")
## BH outliers: 67
plot_manhattan(pcadapt_res_df,
               title_label = "pcadapt: NON-AUTO-FIELD vs AUTO",
               save_path   = "new_aut_pcadapt.pdf")

plot(pcadapt_res, option = "qqplot")

plot(pcadapt_res, option = "scores")

hist(pcadapt_res$pvalues, xlab = "p-values", main = NULL, breaks = 50, col = "pink")

compare_thresholds(pcadapt_res)
## q-values (alpha=0.05):     NA 
## BH adjusted (alpha=0.05):  NA 
## Bonferroni (alpha=0.05):   NA
par(mfrow = c(2, 2))
for (i in 1:2)
  plot(pcadapt_res$loadings[, i], pch = 19, cex = .3, ylab = paste0("Loadings PC", i))

5.1 Venn diagram outFlank vs pcadapt

Note: This section was not used in the final version of the manuscript.

common_SNPs <- venn_outflank_pcadapt(
  outflank_file = "new_aut_SNPs_outFlank.txt",
  pcadapt_file  = "new_aut_SNPs_pcadapt.txt",
  venn_pdf      = "new_aut_significant_snps.pdf",
  common_txt    = "new_aut_common_SNPs_pcadapt_outflank.txt"
)

filtered_df <- pcadapt_res_df |> filter(SNP %in% common_SNPs)
color_vector <- chrom_colors
names(color_vector) <- unique(pcadapt_res_df$Chromosome)
min_pval_snps <- pcadapt_res_df |> filter(highlight) |>
  group_by(Chromosome) |> slice(which.min(pvalues))

ggplot(pcadapt_res_df, aes(x = Position, y = -log10(pvalues))) +
  geom_point(aes(color = Chromosome), data = subset(pcadapt_res_df, !highlight), size = .3) +
  geom_point(color = "magenta", data = subset(filtered_df, highlight), size = .3) +
  scale_color_manual(values = color_vector, guide = "none") +
  labs(x = "Position", y = "-log10(p-value)",
       caption = "Magenta: BH-adjusted (alpha=0.05) outliers shared between pcadapt and outFlank.") +
  facet_wrap(~Chromosome, scales = "free_x") +
  scale_x_continuous(labels = mb_format) +
  my_theme() +
  theme(panel.spacing = unit(.2, "lines"),
        plot.margin = unit(c(1, 1, 2, 2), "lines"),
        plot.caption = element_text(face = "italic")) +
  geom_text_repel(data = min_pval_snps, aes(label = SNP),
                  size = 3, nudge_y = 0.01, segment.color = NA, max.overlaps = Inf)

ggsave(here("output", "pcadapt", "figures", "new_aut_pcadapt_outflank.pdf"),
       width = 8, height = 5, units = "in")
cluster_table(filtered_df, "new_aut_snps_pcadap_outflank.docx")
Chromosome Window_Start SNPs
1 1.4e+08 AX-583515734, AX-583518586, AX-583516491, AX-583517344, AX-583521326, AX-583531703
1 2.2e+08 AX-583750740, AX-583750871
1 2.6e+08 AX-583882507, AX-583878485
1 2.7e+08 AX-583924978
1 3.0e+08 AX-583972796
2 2.0e+07 AX-584399488, AX-584419201
2 3.0e+07 AX-584429924, AX-584428515
2 1.2e+08 AX-585110132
2 1.8e+08 AX-585196879
2 1.9e+08 AX-582457468, AX-582465641
2 2.0e+08 AX-582535272
2 3.2e+08 AX-584907866
2 3.8e+08 AX-579548089
2 3.9e+08 AX-579556477, AX-579560686, AX-579564292, AX-579565469, AX-579584369, AX-579583040, AX-579584883
2 4.0e+08 AX-579596233, AX-579602153, AX-579618571
2 4.1e+08 AX-579630465, AX-579632196, AX-579632983, AX-579632074, AX-579632927, AX-579661979
2 4.2e+08 AX-579681911, AX-579693902, AX-579697016
2 5.0e+08 AX-579940789
2 5.7e+08 AX-580238721
3 0.0e+00 AX-580344997
3 1.0e+07 AX-580398903
3 3.0e+07 AX-580480108
3 4.0e+07 AX-580560428
3 5.0e+07 AX-580575055
3 1.3e+08 AX-580917486, AX-580925458
3 1.5e+08 AX-581030099
3 1.9e+08 AX-581302901
3 2.1e+08 AX-581428225, AX-581437212, AX-581438408, AX-581442470, AX-581462648, AX-581467653
3 2.3e+08 AX-581528715, AX-581534474, AX-581543825
3 2.9e+08 AX-581810815
3 3.2e+08 AX-581929530

6. pcadapt: NON-AUTO vs NON-AUTO-FIELD (pairwise)

pcadapt performs principal component analysis and tests for selection outliers based on correlations between genetic variation and the first K principal components.

rm(list = setdiff(ls(), c("run_pcadapt", "plot_manhattan", "compare_thresholds",
                           "venn_outflank_pcadapt", "cluster_table",
                           "chrom_colors", "mb_format", "my_theme")))
gc()
##           used  (Mb) gc trigger  (Mb) max used  (Mb)
## Ncells 2508947 134.0    4241018 226.5  4241018 226.5
## Vcells 4775760  36.5   16030974 122.4 16030974 122.4
x_mn <- read.pcadapt(here("output", "outflank", "man_new.bed"), type = "bed")
plot(pcadapt(x_mn, K = 30), option = "screeplot")

K = 3 was selected based on scree plot inspection.

out_mn <- run_pcadapt("man_new.bed", "man_new.bim",
                       K = 3, out_snp_file = "man_new_SNPs_pcadapt.txt")
pcadapt_res <- out_mn$res
pcadapt_res_df <- out_mn$df
cat("BH outliers:", sum(pcadapt_res_df$highlight), "\n")
## BH outliers: 82
plot_manhattan(pcadapt_res_df,
               title_label = "pcadapt: NON-AUTO vs NON-AUTO-FIELD",
               save_path   = "man_new_pcadapt.pdf")

plot(pcadapt_res, option = "qqplot")

plot(pcadapt_res, option = "scores")

hist(pcadapt_res$pvalues, xlab = "p-values", main = NULL, breaks = 50, col = "pink")

compare_thresholds(pcadapt_res)
## q-values (alpha=0.05):     NA 
## BH adjusted (alpha=0.05):  NA 
## Bonferroni (alpha=0.05):   NA
par(mfrow = c(2, 2))
for (i in 1:2)
  plot(pcadapt_res$loadings[, i], pch = 19, cex = .3, ylab = paste0("Loadings PC", i))

6.1 Venn diagram outFlank vs pcadapt

Note: This section was not used in the final version of the manuscript.

OutFlank did not identify any outlier for this comparison; the Venn diagram is not produced.

7. Venn diagram across comparisons

Note: This section was not used in the final version of the manuscript.

SNPs shared between the pcadapt+outFlank intersection sets from the three comparisons involving AUTO are combined into a four-way Venn diagram.

NEW_MAN_AUT <- read.table(
  here("output", "pcadapt", "common_SNPs_pcadapt_outflank.txt"),
  stringsAsFactors = FALSE) |> drop_na()
MAN_AUT <- read.table(
  here("output", "pcadapt", "man_aut_common_SNPs_pcadapt_outflank.txt"),
  stringsAsFactors = FALSE)
NEW_AUT <- read.table(
  here("output", "pcadapt", "new_aut_common_SNPs_pcadapt_outflank.txt"),
  stringsAsFactors = FALSE)
NEW_MAN <- read.table(
  here("output", "pcadapt", "man_new_SNPs_pcadapt.txt"),
  stringsAsFactors = FALSE)

lst4 <- list(
  "AUTO+NON-AUTO+NON-AUTO-FIELD" = NEW_MAN_AUT$V1,
  "NON-AUTO vs AUTO"             = MAN_AUT$V1,
  "NON-AUTO-FIELD vs AUTO"       = NEW_AUT$V1,
  "NON-AUTO vs NON-AUTO-FIELD"   = NEW_MAN$V1
)

venn_diagram <- ggvenn(lst4) +
  theme(plot.margin = margin(60, 60, 60, 60, "points"),
        text = element_text(size = 5))
print(venn_diagram)

ggsave(here("output", "pcadapt", "figures", "4_way_venn_significant_snps.pdf"),
       venn_diagram, height = 8, width = 8, dpi = 300)

common_elements <- Reduce(intersect, list(NEW_MAN_AUT$V1, MAN_AUT$V1, NEW_AUT$V1))
write.table(common_elements,
  file = here("output", "pcadapt", "4_way_venn_common_SNPs_pcadapt_outflank.txt"),
  row.names = FALSE, col.names = FALSE, quote = FALSE)

cat("SNP counts per set:\n")
## SNP counts per set:
cat("  AUTO+NON-AUTO+NON-AUTO-FIELD:", nrow(NEW_MAN_AUT), "\n")
##   AUTO+NON-AUTO+NON-AUTO-FIELD: 33
cat("  NON-AUTO vs AUTO:            ", nrow(MAN_AUT), "\n")
##   NON-AUTO vs AUTO:             109
cat("  NON-AUTO-FIELD vs AUTO:      ", nrow(NEW_AUT), "\n")
##   NON-AUTO-FIELD vs AUTO:       64
cat("  NON-AUTO vs NON-AUTO-FIELD:  ", nrow(NEW_MAN), "\n")
##   NON-AUTO vs NON-AUTO-FIELD:   82
color_vector <- chrom_colors
names(color_vector) <- unique(pcadapt_res_df$Chromosome)
filtered_df <- pcadapt_res_df |> filter(SNP %in% common_elements)

ggplot(pcadapt_res_df, aes(x = Position, y = -log10(pvalues))) +
  geom_point(aes(color = Chromosome), data = subset(pcadapt_res_df, !highlight), size = .3) +
  geom_point(color = "magenta", data = subset(filtered_df, highlight), size = .3) +
  scale_color_manual(values = color_vector, guide = "none") +
  labs(x = "Position", y = "-log10(p-value)") +
  facet_wrap(~Chromosome, scales = "free_x") +
  scale_x_continuous(labels = mb_format) +
  my_theme() +
  theme(panel.spacing = unit(.2, "lines"),
        plot.margin = unit(c(1, 1, 2, 2), "lines"))

ggsave(here("output", "pcadapt", "figures", "manhattan_from_4_way_venn_snps.pdf"),
       height = 5, width = 8, dpi = 300)
cluster_table(filtered_df, "4_way_venn_pcadap_outflank.docx")
Chromosome Window_Start SNPs
2 1.8e+08 AX-585196879
2 3.8e+08 AX-579548089
2 3.9e+08 AX-579560686
2 4.1e+08 AX-579661979
3 2.1e+08 AX-581442470, AX-581462648, AX-581467653

Union of outlier SNPs across all three AUTO-involving comparisons:

combined_df <- bind_rows(MAN_AUT, NEW_AUT, NEW_MAN_AUT)
result <- combined_df |>
  group_by(V1) |>
  summarise(Count = n(), .groups = "drop") |>
  arrange(V1) |>
  dplyr::rename(SNP = V1)
cat("Unique SNPs (union):", n_distinct(result$SNP), "\n")
## Unique SNPs (union): 161
# Save union of outlier SNPs for downstream notebooks (S4, S5, S6)
write.table(result$SNP,
  file = here("output", "snpeff", "SNPs_158.txt"),
  row.names = FALSE, col.names = FALSE, quote = FALSE)

filtered_df <- pcadapt_res_df |> filter(SNP %in% result$SNP)
cluster_df <- filtered_df |>
  mutate(Window_Start = floor(Position / 1e7) * 1e7) |>
  group_by(Chromosome, Window_Start) |>
  summarise(SNP_count = n(), .groups = "drop") |>
  dplyr::filter(SNP_count >= 3)

ggplot(pcadapt_res_df, aes(x = Position, y = -log10(pvalues))) +
  geom_point(aes(color = Chromosome), data = subset(pcadapt_res_df, !highlight), size = .3) +
  geom_point(color = "magenta", data = subset(filtered_df, highlight), size = .3) +
  geom_vline(data = cluster_df, aes(xintercept = Window_Start + 5e6),
             color = "lightgray", linetype = "dashed", linewidth = 0.5) +
  geom_text_repel(data = cluster_df,
                  aes(x = Window_Start + 5e6, y = 10, label = SNP_count),
                  vjust = -1, hjust = -0.5) +
  scale_color_manual(values = color_vector, guide = "none") +
  labs(x = "Position", y = "-log10(p-value)",
       caption = "Magenta: union of BH-adjusted outliers across AUTO-involving comparisons. Dashed lines: 10 Mb windows with >= 3 SNPs.") +
  facet_wrap(~Chromosome, scales = "free_x") +
  scale_x_continuous(labels = mb_format) +
  my_theme() +
  theme(panel.spacing = unit(.2, "lines"),
        plot.margin = unit(c(1, 1, 2, 2), "lines"),
        plot.caption = element_text(face = "italic"),
        legend.position = "none")

ggsave(here("output", "pcadapt", "figures", "manhattan_from_4_way_venn_158_snps.pdf"),
       height = 5, width = 8, dpi = 300)
cluster_table(filtered_df, "4_way_venn_pcadap_outflank_158_SNPs.docx")
Chromosome Window_Start SNPs
1 0.0e+00 AX-583054970
1 1.0e+07 AX-583095890
1 9.0e+07 AX-583324654
1 1.0e+08 AX-583355868
1 1.2e+08 AX-583423493, AX-583426050
1 1.3e+08 AX-583467551, AX-583504302
1 1.4e+08 AX-583514148, AX-583518586, AX-583518994, AX-583516491, AX-583517344, AX-583517999, AX-583521326, AX-583531703
1 1.7e+08 AX-583631472
1 2.2e+08 AX-583750740, AX-583750871
1 2.6e+08 AX-583882507, AX-583878485
1 2.7e+08 AX-583889791, AX-583905588, AX-583924569, AX-583924978
1 3.0e+08 AX-583972796
1 3.4e+08 AX-584133664
2 2.0e+07 AX-584399488, AX-584419201, AX-584427590
2 3.0e+07 AX-584429924, AX-584428515, AX-584444558, AX-584453518
2 4.0e+07 AX-584473131
2 5.0e+07 AX-585010439, AX-584502184, AX-585023213
2 7.0e+07 AX-585042322
2 8.0e+07 AX-584546335
2 1.1e+08 AX-585101539
2 1.2e+08 AX-585110132
2 1.7e+08 AX-585182739
2 1.8e+08 AX-585190366, AX-585196879
2 1.9e+08 AX-582457468, AX-582465641
2 2.0e+08 AX-582494788, AX-582519052, AX-582532563, AX-582535272
2 2.6e+08 AX-584765066
2 3.2e+08 AX-585400879, AX-585413917, AX-584907866
2 3.5e+08 AX-584951049
2 3.6e+08 AX-579474142, AX-579474650
2 3.7e+08 AX-579509845
2 3.8e+08 AX-579548089
2 3.9e+08 AX-579556477, AX-579560686, AX-579564292, AX-579565469, AX-579580051
2 4.0e+08 AX-579596233, AX-579602153, AX-579603553, AX-579604110, AX-579606795, AX-579618571, AX-579619995, AX-579620627
2 4.1e+08 AX-579630462, AX-579632983, AX-579632074, AX-579636106, AX-579638540, AX-579640459, AX-579646083, AX-579661979, AX-579662371, AX-579667335
2 4.2e+08 AX-579693902, AX-579697016
2 4.5e+08 AX-579768613, AX-579778684
2 4.6e+08 AX-579808128
2 4.9e+08 AX-579909895
2 5.0e+08 AX-579940789
2 5.1e+08 AX-579955661, AX-579982534
2 5.2e+08 AX-580001880, AX-580012078, AX-580016916, AX-580026657, AX-580027537, AX-580042929
2 5.4e+08 AX-580099245
2 5.5e+08 AX-580165355
2 5.6e+08 AX-580195612
2 5.7e+08 AX-580238721
2 5.8e+08 AX-580318010, AX-580321888
3 0.0e+00 AX-580342672, AX-580344997
3 1.0e+07 AX-580398903
3 2.0e+07 AX-580425125
3 4.0e+07 AX-580560428
3 5.0e+07 AX-580575055
3 7.0e+07 AX-580620056
3 8.0e+07 AX-580673467
3 9.0e+07 AX-580719743
3 1.3e+08 AX-580917486, AX-580925458, AX-580926966
3 1.5e+08 AX-581030099
3 1.8e+08 AX-581275893
3 1.9e+08 AX-581298415, AX-581302704, AX-581302901
3 2.0e+08 AX-581399096
3 2.1e+08 AX-581408938, AX-581428225, AX-581442470, AX-581458697, AX-581461793, AX-581462648, AX-581467653, AX-581470736
3 2.2e+08 AX-581485809, AX-581488185, AX-581504582
3 2.3e+08 AX-581528715, AX-581543825, AX-581572119
3 2.9e+08 AX-581810815
3 3.2e+08 AX-581929530
3 3.5e+08 AX-582085824
3 4.4e+08 AX-582853985
3 4.5e+08 AX-582884221
3 4.8e+08 AX-583025646