Overview

This notebook implements the OutFLANK selection scan (Whitlock & Lotterhos 2015) to identify FST outlier loci in Aedes albopictus autogeny comparisons. The method fits a chi-squared distribution to a trimmed set of FST values to characterise the neutral FST distribution, then assigns q-values to identify outliers consistent with spatially heterogeneous selection.

Population label key: AUT = AUTO, MAN = NON-AUTO, NEW = NON-AUTO-FIELD.

1. Libraries

library(tidyverse)
library(here)
library(ggrepel)
# library(ggtext)  # not available during render
library(vcfR)
library(OutFLANK)
library(ggvenn)

2. Helper functions

source(here("notebooks", "helpers", "my_theme2.R"))

# Load a VCF, recode diploid genotypes to 0/1/2/9, and return a named list
# with the genotype matrix G (transposed for OutFLANK), SNP IDs, positions,
# and population labels derived from sample name prefixes.
vcf_to_genotype_matrix <- function(vcf_path) {
  vcf       <- read.vcfR(vcf_path)
  geno      <- vcfR::extract.gt(vcf)
  snpNames  <- vcfR::getID(vcf)
  positions <- vcfR::getPOS(vcf)
  popNames  <- sapply(strsplit(colnames(geno), "_"), "[", 1)

  G <- matrix(9L, nrow = nrow(geno), ncol = ncol(geno))
  G[geno %in% c("0/0", "0|0")] <- 0L
  G[geno %in% c("0/1", "1/0", "1|0", "0|1")] <- 1L
  G[geno %in% c("1/1", "1|1")] <- 2L

  list(G = t(G), snpNames = snpNames, positions = positions, popNames = popNames)
}

# Apply the calibrated neutral FST distribution to the full SNP set,
# identify outliers, and return a data frame with OutlierFlag and q-values.
run_outflank_scan <- function(vcf_path, out_trim) {
  dat  <- vcf_to_genotype_matrix(vcf_path)
  locusNames3 <- paste(dat$snpNames, dat$positions, sep = ".")
  my_fst4 <- MakeDiploidFSTMat(dat$G, locusNames = locusNames3, popNames = dat$popNames)

  P1 <- pOutlierFinderChiSqNoCorr(
    my_fst4,
    Fstbar    = out_trim$FSTNoCorrbar,
    dfInferred = out_trim$dfInferred,
    qthreshold = 0.05,
    Hmin       = 0.1
  )
  P1 |> separate(LocusName, into = c("SNP_id", "LocusName"), sep = "\\.")
}

# Produce a Manhattan plot with outlier loci highlighted in magenta,
# annotated with SNP IDs at the highest-FST outlier per chromosome.
plot_outflank_manhattan <- function(P1, bim_path, title = NULL) {
  snps <- read_delim(
    bim_path,
    col_names      = FALSE,
    show_col_types = FALSE,
    col_types      = "ccidcc"
  )
  colnames(snps) <- c("Scaffold", "SNP", "Cm", "Position", "Allele1", "Allele2")

  my_fst4_split <- P1
  merged_df <- merge(my_fst4_split, snps, by.x = "SNP_id", by.y = "SNP") |>
    rename(Chromosome = Scaffold)

  df_sub  <- subset(merged_df, He > 0.1)
  my_out  <- P1$SNP_id[P1$OutlierFlag == TRUE]
  df_sub$highlight <- df_sub$SNP_id %in% my_out

  color_vector <- c("#CCF6D6", "#F6E1CC", "#CCD8F6")
  names(color_vector) <- unique(df_sub$Chromosome)

  highest_FST <- df_sub[df_sub$highlight, ]
  highest_FST <- highest_FST[
    highest_FST$FST == ave(highest_FST$FST, highest_FST$Chromosome, FUN = max), ]

  k_format <- function(x) paste0(format(x / 1e6, big.mark = "", scientific = FALSE), "k")

  ggplot(df_sub, aes(x = Position, y = FST)) +
    geom_point(aes(color = Chromosome), data = subset(df_sub, !highlight), size = .3) +
    geom_point(color = "magenta", data = subset(df_sub, highlight), size = .3) +
    scale_color_manual(values = color_vector, guide = "none") +
    labs(x = "Position", y = expression(F[ST]), title = title) +
    facet_wrap(~ Chromosome, scales = "free_x") +
    scale_x_continuous(labels = k_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 = highest_FST,
      aes(label = SNP_id),
      size = 3, nudge_y = 0.01, segment.color = NA, max.overlaps = Inf
    )
}

3. OutFLANK: all three populations (MAN, NEW, AUT)

OutFLANK identifies FST outlier loci by fitting a chi-squared distribution to a trimmed set of FST values from quasi-independent SNPs, then applying the inferred neutral FST distribution to the complete SNP set.

3.1 LD pruning and neutral SNP calibration set

Extract intergenic SNPs and perform LD pruning to obtain a quasi-independent set.

plink2 \
--bfile output/quality_control/file7 \
--extract data/files/intergenic_SNPs.txt \
--indep-pairwise 5 1 0.1 \
--out output/outflank/indepSNP \
--silent;
grep 'samples\|variants\|remaining' output/outflank/indepSNP.log

## 60 samples (24 females, 26 males, 10 ambiguous; 60 founders) loaded from
## 110353 variants loaded from output/quality_control/file7.bim.
## --extract: 12112 variants remaining.
## 12112 variants remaining after main filters.
## --indep-pairwise (3 compute threads): 6667/12112 variants removed.

plink2 \
--bfile output/quality_control/file7 \
--extract output/outflank/indepSNP.prune.in \
--export vcf \
--make-bed \
--maf 0.1 \
--geno 0.2 \
--out output/outflank/intergenic \
--silent;
grep 'samples\|variants\|remaining' output/outflank/intergenic.log

## 60 samples (24 females, 26 males, 10 ambiguous; 60 founders) loaded from
## 110353 variants loaded from output/quality_control/file7.bim.
## --extract: 5445 variants remaining.
## --geno: 0 variants removed due to missing genotype data.
## 957 variants removed due to allele frequency threshold(s)
## 4488 variants remaining after main filters.

Compute FST values for the intergenic LD-pruned set and select low-FST loci (FST < 0.2) as the neutral calibration set.

dat_intergenic <- vcf_to_genotype_matrix(here("output", "outflank", "intergenic.vcf"))

## Scanning file to determine attributes.
## File attributes:
##   meta lines: 8
##   header_line: 9
##   variant count: 4488
##   column count: 69
## Meta line 8 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
##   Character matrix gt rows: 4488
##   Character matrix gt cols: 69
##   skip: 0
##   nrows: 4488
##   row_num: 0
## Processed variant 1000Processed variant 2000Processed variant 3000Processed variant 4000Processed variant: 4488
## All variants processed

my_fst_intergenic <- MakeDiploidFSTMat(
  dat_intergenic$G,
  locusNames = dat_intergenic$snpNames,
  popNames   = dat_intergenic$popNames
)

## Calculating FSTs, may take a few minutes...

neutral_snps <- my_fst_intergenic |>
  filter(FST < 0.2 & FST >= -0.1) |>
  pull(LocusName)

write.table(
  neutral_snps,
  file      = here("output", "outflank", "neutral_SNPs.txt"),
  row.names = FALSE,
  quote     = FALSE,
  col.names = FALSE,
  sep       = "\n"
)

plink2 \
--bfile output/outflank/intergenic \
--export vcf \
--extract output/outflank/neutral_SNPs.txt \
--out output/outflank/neutral \
--silent;
grep "samples\|variants" output/outflank/neutral.log

## 60 samples (24 females, 26 males, 10 ambiguous; 60 founders) loaded from
## 4488 variants loaded from output/outflank/intergenic.bim.
## --extract: 3722 variants remaining.
## 3722 variants remaining after main filters.

Calibrate the neutral FST distribution.

dat_neutral <- vcf_to_genotype_matrix(here("output", "outflank", "neutral.vcf"))

## Scanning file to determine attributes.
## File attributes:
##   meta lines: 8
##   header_line: 9
##   variant count: 3722
##   column count: 69
## Meta line 8 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
##   Character matrix gt rows: 3722
##   Character matrix gt cols: 69
##   skip: 0
##   nrows: 3722
##   row_num: 0
## Processed variant 1000Processed variant 2000Processed variant 3000Processed variant: 3722
## All variants processed

my_fst_neutral <- MakeDiploidFSTMat(
  dat_neutral$G,
  locusNames = dat_neutral$snpNames,
  popNames   = dat_neutral$popNames
)

## Calculating FSTs, may take a few minutes...

# NOTE: NumberOfSamples=60 is preserved from the original analysis.
# FLAG: OutFLANK documentation specifies this should be the number of
# populations (here: 3), not the number of individuals (60). This parameter
# may affect the fitted chi-squared degrees of freedom and outlier calls.
out_trim_3pop <- OutFLANK(
  FstDataFrame    = my_fst_neutral,
  Hmin            = 0.1,
  NumberOfSamples = 60,
  qthreshold      = 0.05
)

Diagnostic plots for the neutral FST fit.

OutFLANKResultsPlotter(
  out_trim_3pop,
  withOutliers      = TRUE,
  NoCorr            = TRUE,
  Hmin              = 0.1,
  binwidth          = 0.001,
  Zoom              = FALSE,
  RightZoomFraction = 0.05,
  titletext         = NULL
)

OutFLANKResultsPlotter(
  out_trim_3pop,
  withOutliers      = TRUE,
  NoCorr            = TRUE,
  Hmin              = 0.1,
  binwidth          = 0.001,
  Zoom              = TRUE,
  RightZoomFraction = 0.05,
  titletext         = NULL
)

hist(out_trim_3pop$results$pvaluesRightTail,
     xlab = "P-value (right tail)", main = "3-population neutral calibration set")

3.2 Scan across all loci

Apply the calibrated neutral FST distribution to the full LD-pruned SNP set.

plink2 \
--bfile output/quality_control/file7 \
--export vcf \
--extract output/quality_control/indepSNP.prune.in \
--make-bed \
--out output/outflank/autogenous \
--silent;
grep "samples\|variants" output/outflank/autogenous.log

## 60 samples (24 females, 26 males, 10 ambiguous; 60 founders) loaded from
## 110353 variants loaded from output/quality_control/file7.bim.
## --extract: 41849 variants remaining.
## 41849 variants remaining after main filters.

P1_3pop <- run_outflank_scan(
  here("output", "outflank", "autogenous.vcf"),
  out_trim_3pop
)

## Scanning file to determine attributes.
## File attributes:
##   meta lines: 8
##   header_line: 9
##   variant count: 41849
##   column count: 69
## Meta line 8 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
##   Character matrix gt rows: 41849
##   Character matrix gt cols: 69
##   skip: 0
##   nrows: 41849
##   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 26000Processed variant 27000Processed variant 28000Processed variant 29000Processed variant 30000Processed variant 31000Processed variant 32000Processed variant 33000Processed variant 34000Processed variant 35000Processed variant 36000Processed variant 37000Processed variant 38000Processed variant 39000Processed variant 40000Processed variant 41000Processed variant: 41849
## All variants processed
## Calculating FSTs, may take a few minutes...
## [1] "10000 done of 41849"
## [1] "20000 done of 41849"
## [1] "30000 done of 41849"
## [1] "40000 done of 41849"

my_out_3pop <- P1_3pop$SNP_id[P1_3pop$OutlierFlag == TRUE]
write.table(
  my_out_3pop,
  file      = here("output", "outflank", "SNPs_outFlank.txt"),
  row.names = FALSE,
  quote     = FALSE,
  col.names = FALSE,
  sep       = "\n"
)

p_3pop <- plot_outflank_manhattan(
  P1_3pop,
  here("output", "outflank", "autogenous.bim"),
  title = "OutFLANK: AUTO vs NON-AUTO vs NON-AUTO-FIELD"
)
print(p_3pop)

ggsave(
  here("output", "outflank", "figures", "outFlank_outliers.pdf"),
  width = 8, height = 5, units = "in"
)

4. OutFLANK: NON-AUTO vs AUTO (MAN vs AUT)

4.1 Neutral SNP calibration set

plink2 \
--bfile output/outflank/man_aut \
--extract output/outflank/indepSNP.prune.in \
--export vcf \
--make-bed \
--out output/outflank/man_aut2 \
--silent;
grep 'samples\|variants\|remaining' output/outflank/man_aut2.log

## 38 samples (13 females, 15 males, 10 ambiguous; 38 founders) loaded from
## 41332 variants loaded from output/outflank/man_aut.bim.
## --extract: 3705 variants remaining.
## 3705 variants remaining after main filters.

dat_ma2 <- vcf_to_genotype_matrix(here("output", "outflank", "man_aut2.vcf"))

## Scanning file to determine attributes.
## File attributes:
##   meta lines: 8
##   header_line: 9
##   variant count: 3705
##   column count: 47
## Meta line 8 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
##   Character matrix gt rows: 3705
##   Character matrix gt cols: 47
##   skip: 0
##   nrows: 3705
##   row_num: 0
## Processed variant 1000Processed variant 2000Processed variant 3000Processed variant: 3705
## All variants processed

my_fst_ma2 <- MakeDiploidFSTMat(
  dat_ma2$G, locusNames = dat_ma2$snpNames, popNames = dat_ma2$popNames
)

## Calculating FSTs, may take a few minutes...

neutral_ma <- my_fst_ma2 |> filter(FST < 0.2 & FST >= -0.1) |> pull(LocusName)
write.table(
  neutral_ma,
  file      = here("output", "outflank", "man_aut2_SNPs.txt"),
  row.names = FALSE, quote = FALSE, col.names = FALSE, sep = "\n"
)

plink2 \
--bfile output/outflank/man_aut2 \
--export vcf \
--extract output/outflank/man_aut2_SNPs.txt \
--out output/outflank/man_aut3 \
--silent;
grep "samples\|variants" output/outflank/man_aut3.log

## 38 samples (13 females, 15 males, 10 ambiguous; 38 founders) loaded from
## 3705 variants loaded from output/outflank/man_aut2.bim.
## --extract: 2767 variants remaining.
## 2767 variants remaining after main filters.

dat_ma3 <- vcf_to_genotype_matrix(here("output", "outflank", "man_aut3.vcf"))

## Scanning file to determine attributes.
## File attributes:
##   meta lines: 8
##   header_line: 9
##   variant count: 2767
##   column count: 47
## Meta line 8 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
##   Character matrix gt rows: 2767
##   Character matrix gt cols: 47
##   skip: 0
##   nrows: 2767
##   row_num: 0
## Processed variant 1000Processed variant 2000Processed variant: 2767
## All variants processed

my_fst_ma3 <- MakeDiploidFSTMat(
  dat_ma3$G, locusNames = dat_ma3$snpNames, popNames = dat_ma3$popNames
)

## Calculating FSTs, may take a few minutes...

# NOTE: NumberOfSamples=60 preserved. FLAG: should be 2 (populations) for pairwise comparisons.
out_trim_ma <- OutFLANK(
  FstDataFrame    = my_fst_ma3,
  Hmin            = 0.1,
  NumberOfSamples = 60,
  qthreshold      = 0.05
)

OutFLANKResultsPlotter(
  out_trim_ma,
  withOutliers = TRUE, NoCorr = TRUE, Hmin = 0.1,
  binwidth = 0.001, Zoom = FALSE, RightZoomFraction = 0.05, titletext = NULL
)

OutFLANKResultsPlotter(
  out_trim_ma,
  withOutliers = TRUE, NoCorr = TRUE, Hmin = 0.1,
  binwidth = 0.001, Zoom = TRUE, RightZoomFraction = 0.05, titletext = NULL
)

hist(out_trim_ma$results$pvaluesRightTail,
     xlab = "P-value (right tail)", main = "NON-AUTO vs AUTO neutral calibration set")

4.2 Scan across all loci

plink2 \
--bfile output/outflank/man_aut \
--export vcf \
--make-bed \
--out output/outflank/man_aut4 \
--silent;
grep "samples\|variants" output/outflank/man_aut4.log

## 38 samples (13 females, 15 males, 10 ambiguous; 38 founders) loaded from
## 41332 variants loaded from output/outflank/man_aut.bim.

P1_ma <- run_outflank_scan(
  here("output", "outflank", "man_aut4.vcf"),
  out_trim_ma
)

## Scanning file to determine attributes.
## File attributes:
##   meta lines: 8
##   header_line: 9
##   variant count: 41332
##   column count: 47
## Meta line 8 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
##   Character matrix gt rows: 41332
##   Character matrix gt cols: 47
##   skip: 0
##   nrows: 41332
##   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 26000Processed variant 27000Processed variant 28000Processed variant 29000Processed variant 30000Processed variant 31000Processed variant 32000Processed variant 33000Processed variant 34000Processed variant 35000Processed variant 36000Processed variant 37000Processed variant 38000Processed variant 39000Processed variant 40000Processed variant 41000Processed variant: 41332
## All variants processed
## Calculating FSTs, may take a few minutes...
## [1] "10000 done of 41332"
## [1] "20000 done of 41332"
## [1] "30000 done of 41332"
## [1] "40000 done of 41332"

my_out_ma <- P1_ma$SNP_id[P1_ma$OutlierFlag == TRUE]
write.table(
  my_out_ma,
  file      = here("output", "outflank", "man_aut_SNPs_outFlank.txt"),
  row.names = FALSE, quote = FALSE, col.names = FALSE, sep = "\n"
)

p_ma <- plot_outflank_manhattan(
  P1_ma,
  here("output", "outflank", "man_aut.bim"),
  title = "OutFLANK: NON-AUTO vs AUTO"
)
print(p_ma)

ggsave(
  here("output", "outflank", "figures", "man_aut_outFlank_outliers.pdf"),
  width = 8, height = 5, units = "in"
)

5. OutFLANK: NON-AUTO-FIELD vs AUTO (NEW vs AUT)

5.1 Neutral SNP calibration set

plink2 \
--bfile output/outflank/new_aut \
--extract output/outflank/indepSNP.prune.in \
--export vcf \
--make-bed \
--out output/outflank/new_aut2 \
--silent;
grep 'samples\|variants\|remaining' output/outflank/new_aut2.log

## 50 samples (24 females, 26 males; 50 founders) loaded from
## 41696 variants loaded from output/outflank/new_aut.bim.
## --extract: 3724 variants remaining.
## 3724 variants remaining after main filters.

dat_na2 <- vcf_to_genotype_matrix(here("output", "outflank", "new_aut2.vcf"))

## Scanning file to determine attributes.
## File attributes:
##   meta lines: 8
##   header_line: 9
##   variant count: 3724
##   column count: 59
## Meta line 8 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
##   Character matrix gt rows: 3724
##   Character matrix gt cols: 59
##   skip: 0
##   nrows: 3724
##   row_num: 0
## Processed variant 1000Processed variant 2000Processed variant 3000Processed variant: 3724
## All variants processed

my_fst_na2 <- MakeDiploidFSTMat(
  dat_na2$G, locusNames = dat_na2$snpNames, popNames = dat_na2$popNames
)

## Calculating FSTs, may take a few minutes...

neutral_na <- my_fst_na2 |> filter(FST < 0.2 & FST >= -0.1) |> pull(LocusName)
write.table(
  neutral_na,
  file      = here("output", "outflank", "new_aut2_SNPs.txt"),
  row.names = FALSE, quote = FALSE, col.names = FALSE, sep = "\n"
)

plink2 \
--bfile output/outflank/new_aut2 \
--export vcf \
--extract output/outflank/new_aut2_SNPs.txt \
--out output/outflank/new_aut3 \
--silent;
grep "samples\|variants" output/outflank/new_aut3.log

## 50 samples (24 females, 26 males; 50 founders) loaded from
## 3724 variants loaded from output/outflank/new_aut2.bim.
## --extract: 3078 variants remaining.
## 3078 variants remaining after main filters.

dat_na3 <- vcf_to_genotype_matrix(here("output", "outflank", "new_aut3.vcf"))

## Scanning file to determine attributes.
## File attributes:
##   meta lines: 8
##   header_line: 9
##   variant count: 3078
##   column count: 59
## Meta line 8 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
##   Character matrix gt rows: 3078
##   Character matrix gt cols: 59
##   skip: 0
##   nrows: 3078
##   row_num: 0
## Processed variant 1000Processed variant 2000Processed variant 3000Processed variant: 3078
## All variants processed

my_fst_na3 <- MakeDiploidFSTMat(
  dat_na3$G, locusNames = dat_na3$snpNames, popNames = dat_na3$popNames
)

## Calculating FSTs, may take a few minutes...

# NOTE: NumberOfSamples=60 preserved. FLAG: should be 2 (populations) for pairwise comparisons.
out_trim_na <- OutFLANK(
  FstDataFrame    = my_fst_na3,
  Hmin            = 0.1,
  NumberOfSamples = 60,
  qthreshold      = 0.05
)

OutFLANKResultsPlotter(
  out_trim_na,
  withOutliers = TRUE, NoCorr = TRUE, Hmin = 0.1,
  binwidth = 0.001, Zoom = FALSE, RightZoomFraction = 0.05, titletext = NULL
)

OutFLANKResultsPlotter(
  out_trim_na,
  withOutliers = TRUE, NoCorr = TRUE, Hmin = 0.1,
  binwidth = 0.001, Zoom = TRUE, RightZoomFraction = 0.05, titletext = NULL
)

hist(out_trim_na$results$pvaluesRightTail,
     xlab = "P-value (right tail)", main = "NON-AUTO-FIELD vs AUTO neutral calibration set")

5.2 Scan across all loci

plink2 \
--bfile output/outflank/new_aut \
--export vcf \
--make-bed \
--out output/outflank/new_aut4 \
--silent;
grep "samples\|variants" output/outflank/new_aut4.log

## 50 samples (24 females, 26 males; 50 founders) loaded from
## 41696 variants loaded from output/outflank/new_aut.bim.

P1_na <- run_outflank_scan(
  here("output", "outflank", "new_aut4.vcf"),
  out_trim_na
)

## Scanning file to determine attributes.
## File attributes:
##   meta lines: 8
##   header_line: 9
##   variant count: 41696
##   column count: 59
## Meta line 8 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
##   Character matrix gt rows: 41696
##   Character matrix gt cols: 59
##   skip: 0
##   nrows: 41696
##   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 26000Processed variant 27000Processed variant 28000Processed variant 29000Processed variant 30000Processed variant 31000Processed variant 32000Processed variant 33000Processed variant 34000Processed variant 35000Processed variant 36000Processed variant 37000Processed variant 38000Processed variant 39000Processed variant 40000Processed variant 41000Processed variant: 41696
## All variants processed
## Calculating FSTs, may take a few minutes...
## [1] "10000 done of 41696"
## [1] "20000 done of 41696"
## [1] "30000 done of 41696"
## [1] "40000 done of 41696"

my_out_na <- P1_na$SNP_id[P1_na$OutlierFlag == TRUE]
write.table(
  my_out_na,
  file      = here("output", "outflank", "new_aut_SNPs_outFlank.txt"),
  row.names = FALSE, quote = FALSE, col.names = FALSE, sep = "\n"
)

p_na <- plot_outflank_manhattan(
  P1_na,
  here("output", "outflank", "new_aut.bim"),
  title = "OutFLANK: NON-AUTO-FIELD vs AUTO"
)
print(p_na)

ggsave(
  here("output", "outflank", "figures", "new_aut_outFlank_outliers.pdf"),
  width = 8, height = 5, units = "in"
)

6. OutFLANK: NON-AUTO vs NON-AUTO-FIELD (MAN vs NEW)

6.1 Neutral SNP calibration set

plink2 \
--bfile output/outflank/man_new \
--extract output/outflank/indepSNP.prune.in \
--export vcf \
--make-bed \
--out output/outflank/man_new2 \
--silent;
grep 'samples\|variants\|remaining' output/outflank/man_new2.log

## 32 samples (11 females, 11 males, 10 ambiguous; 32 founders) loaded from
## 41206 variants loaded from output/outflank/man_new.bim.
## --extract: 3720 variants remaining.
## 3720 variants remaining after main filters.

dat_mn2 <- vcf_to_genotype_matrix(here("output", "outflank", "man_new2.vcf"))

## Scanning file to determine attributes.
## File attributes:
##   meta lines: 8
##   header_line: 9
##   variant count: 3720
##   column count: 41
## Meta line 8 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
##   Character matrix gt rows: 3720
##   Character matrix gt cols: 41
##   skip: 0
##   nrows: 3720
##   row_num: 0
## Processed variant 1000Processed variant 2000Processed variant 3000Processed variant: 3720
## All variants processed

my_fst_mn2 <- MakeDiploidFSTMat(
  dat_mn2$G, locusNames = dat_mn2$snpNames, popNames = dat_mn2$popNames
)

## Calculating FSTs, may take a few minutes...

neutral_mn <- my_fst_mn2 |> filter(FST < 0.2 & FST >= -0.1) |> pull(LocusName)
write.table(
  neutral_mn,
  file      = here("output", "outflank", "man_new2_SNPs.txt"),
  row.names = FALSE, quote = FALSE, col.names = FALSE, sep = "\n"
)

plink2 \
--bfile output/outflank/man_new2 \
--export vcf \
--extract output/outflank/man_new2_SNPs.txt \
--out output/outflank/man_new3 \
--silent;
grep "samples\|variants" output/outflank/man_new3.log

## 32 samples (11 females, 11 males, 10 ambiguous; 32 founders) loaded from
## 3720 variants loaded from output/outflank/man_new2.bim.
## --extract: 3452 variants remaining.
## 3452 variants remaining after main filters.

dat_mn3 <- vcf_to_genotype_matrix(here("output", "outflank", "man_new3.vcf"))

## Scanning file to determine attributes.
## File attributes:
##   meta lines: 8
##   header_line: 9
##   variant count: 3452
##   column count: 41
## Meta line 8 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
##   Character matrix gt rows: 3452
##   Character matrix gt cols: 41
##   skip: 0
##   nrows: 3452
##   row_num: 0
## Processed variant 1000Processed variant 2000Processed variant 3000Processed variant: 3452
## All variants processed

my_fst_mn3 <- MakeDiploidFSTMat(
  dat_mn3$G, locusNames = dat_mn3$snpNames, popNames = dat_mn3$popNames
)

## Calculating FSTs, may take a few minutes...

# NOTE: NumberOfSamples=60 preserved. FLAG: should be 2 (populations) for pairwise comparisons.
out_trim_mn <- OutFLANK(
  FstDataFrame    = my_fst_mn3,
  Hmin            = 0.1,
  NumberOfSamples = 60,
  qthreshold      = 0.05
)

OutFLANKResultsPlotter(
  out_trim_mn,
  withOutliers = TRUE, NoCorr = TRUE, Hmin = 0.1,
  binwidth = 0.001, Zoom = FALSE, RightZoomFraction = 0.05, titletext = NULL
)

OutFLANKResultsPlotter(
  out_trim_mn,
  withOutliers = TRUE, NoCorr = TRUE, Hmin = 0.1,
  binwidth = 0.001, Zoom = TRUE, RightZoomFraction = 0.05, titletext = NULL
)

hist(out_trim_mn$results$pvaluesRightTail,
     xlab = "P-value (right tail)", main = "NON-AUTO vs NON-AUTO-FIELD neutral calibration set")

6.2 Scan across all loci

No FST outliers are expected between the two non-autogenous populations.

plink2 \
--bfile output/outflank/man_new \
--export vcf \
--make-bed \
--out output/outflank/man_new4 \
--silent;
grep "samples\|variants" output/outflank/man_new4.log

## 32 samples (11 females, 11 males, 10 ambiguous; 32 founders) loaded from
## 41206 variants loaded from output/outflank/man_new.bim.

P1_mn <- run_outflank_scan(
  here("output", "outflank", "man_new4.vcf"),
  out_trim_mn
)

## Scanning file to determine attributes.
## File attributes:
##   meta lines: 8
##   header_line: 9
##   variant count: 41206
##   column count: 41
## Meta line 8 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
##   Character matrix gt rows: 41206
##   Character matrix gt cols: 41
##   skip: 0
##   nrows: 41206
##   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 26000Processed variant 27000Processed variant 28000Processed variant 29000Processed variant 30000Processed variant 31000Processed variant 32000Processed variant 33000Processed variant 34000Processed variant 35000Processed variant 36000Processed variant 37000Processed variant 38000Processed variant 39000Processed variant 40000Processed variant 41000Processed variant: 41206
## All variants processed
## Calculating FSTs, may take a few minutes...
## [1] "10000 done of 41206"
## [1] "20000 done of 41206"
## [1] "30000 done of 41206"
## [1] "40000 done of 41206"

7. Venn diagrams

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

MAN_AUT <- read.table(
  here("output", "outflank", "man_aut_SNPs_outFlank.txt"),
  stringsAsFactors = FALSE
) |> drop_na()

NEW_AUT <- read.table(
  here("output", "outflank", "new_aut_SNPs_outFlank.txt"),
  stringsAsFactors = FALSE
)

NEW_MAN_AUT <- read.table(
  here("output", "outflank", "SNPs_outFlank.txt"),
  stringsAsFactors = FALSE
)

list_of_clusters <- list(
  "NON-AUTO vs AUTO"       = MAN_AUT$V1,
  "NON-AUTO-FIELD vs AUTO" = NEW_AUT$V1,
  "All three populations"  = NEW_MAN_AUT$V1
)

venn_diagram <- ggvenn(list_of_clusters, fill_color = c("steelblue", "darkorange", "pink"))
print(venn_diagram)

common_SNPs <- Reduce(intersect, list_of_clusters)
write.table(
  common_SNPs,
  file      = here("output", "pcadapt", "common_SNPs_NEW_MAN_AUT.txt"),
  row.names = FALSE, col.names = FALSE, quote = FALSE
)

output_path <- here("output", "pcadapt", "figures", "significant_snps_NEW_MAN_AUT.pdf")
ggsave(output_path, venn_diagram, height = 5, width = 5, dpi = 300)

Autogenous - Selection scan with OutFLANK

Luciano V Cosme

2026-04-05