Autogenous - quality control.

1. Getting started: R libraries and software for QC
1.1 R libraries and software
1.2 Software
2. The data
2.1 Use Plink2 to convert to bed format
2.2 Use R to update the .fam file
2.3 Checking the number of samples and localities
3. Quality control steps
3.1 Missingness
3.2 Minor allele frequency
3.3 HWE test
3.4 LD pruning
3.6 Heterozygosity
3.7 Relatedness
3.8 Quick PCA with Plink using the LD pruned data
4. Create chromosomal scale
5. Plot SNP density
6. Create data sets for downstream analysis

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. Getting started: R libraries and software for QC

1.1 R libraries and software

library(tidyverse)

## Warning in system("timedatectl", intern = TRUE): running command 'timedatectl'
## had status 1 and error message 'Function not implemented'

library(here)
library(scales)
library(RColorBrewer)
library(ggrepel)
library(adegenet)

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

1.2 Software

All analyses were performed inside a reproducible container (ghcr.io/cosmelab/albopictus-autogeny:latest) with PLINK 1.9, PLINK 2, bcftools, samtools, bedtools, snpEff, and R 4.2. See the Dockerfile and pixi.toml for the full dependency list.

Quality control follows the general guidelines of Marees et al. (2018).

2. The data

The data was prepared with family and individual IDs assigned, and reference alleles set to match the AalbF3 genome assembly.

bcftools query -l data/genotype_calls/autogenous.vcf | wc -l

## 61

2.1 Use Plink2 to convert to bed format

plink2 \
--allow-extra-chr \
--vcf data/genotype_calls/autogenous.vcf \
--const-fid \
--make-bed \
--exclude data/files/albopictus_SNPs_fail_segregation.txt \
--fa data/genome/albo.fasta.gz \
--ref-from-fa 'force' `# sets REF alleles when it can be done unambiguously, we use force to change the alleles` \
--out output/quality_control/file1 \
--silent;
grep "variants" output/quality_control/file1.log

## --vcf: 131048 variants scanned.
## 131048 variants loaded from output/quality_control/file1-temporary.pvar.zst.
## --exclude: 129381 variants remaining.
## 129381 variants remaining after main filters.
## --ref-from-fa force: 0 variants changed, 129381 validated.

2.2 Use R to update the .fam file

Sample metadata was imported to update the .fam file with population assignments.

samples <-
  read.delim(
    file   = here(
      "data",
      "files",
      "meta_data.txt"
    ),
    header = TRUE
  )
head(samples)

##   Sample.Filename Family_ID Individual_ID Father_ID Mother_ID Sex
## 1 306_MAN_USA.CEL       MAN           306         0         0   0
## 2 307_MAN_USA.CEL       MAN           307         0         0   0
## 3 308_MAN_USA.CEL       MAN           308         0         0   0
## 4 309_MAN_USA.CEL       MAN           309         0         0   0
## 5 310_MAN_USA.CEL       MAN           310         0         0   0
## 6 311_MAN_USA.CEL       MAN           311         0         0   0
##   Affection.Status
## 1               -9
## 2               -9
## 3               -9
## 4               -9
## 5               -9
## 6               -9

fam1 <-
  read.delim(
    file   = here(
      "output", "quality_control", "file1.fam"
    ),
    header = FALSE,

  )
head(fam1)

##   V1              V2 V3 V4 V5 V6
## 1  0 306_MAN_USA.CEL  0  0  0 -9
## 2  0 307_MAN_USA.CEL  0  0  0 -9
## 3  0 308_MAN_USA.CEL  0  0  0 -9
## 4  0 309_MAN_USA.CEL  0  0  0 -9
## 5  0 310_MAN_USA.CEL  0  0  0 -9
## 6  0 311_MAN_USA.CEL  0  0  0 -9

# Extract the number part from the columns
fam1_temp <- fam1 |>
  mutate(num_id = as.numeric(str_extract(V2, "^\\d+")))

samples_temp <- samples |>
  mutate(num_id = as.numeric(str_extract(Sample.Filename, "^\\d+")))

# Perform the left join using the num_id columns and keep the order of fam1
df <- fam1_temp |>
  dplyr::left_join(samples_temp, by = "num_id") |>
  dplyr::select(-num_id) |>
  dplyr::select(8:13)

head(df)

##   Family_ID Individual_ID Father_ID Mother_ID Sex Affection.Status
## 1       MAN           306         0         0   0               -9
## 2       MAN           307         0         0   0               -9
## 3       MAN           308         0         0   0               -9
## 4       MAN           309         0         0   0               -9
## 5       MAN           310         0         0   0               -9
## 6       MAN           311         0         0   0               -9

nrow(df)

## [1] 61

The original .fam is overwritten with updated metadata.

# Update .fam with sample metadata
write.table(
  df,
  file      = here(
    "output", "quality_control", "file1.fam"
  ),
  sep       = "\t",
  row.names = FALSE,
  col.names = FALSE,
  quote     = FALSE
)

2.3 Checking the number of samples and localities

awk '{print $1}' output/quality_control/file1.fam | sort | uniq -c | awk '{print $2, $1}'

## AUT 29
## MAN 10
## NEW 22

3. Quality control steps

3.1 Missingness

plink2 \
--allow-extra-chr \
--bfile output/quality_control/file1 \
--geno 0.2                `# we set genotype missiningness to 20% with this option` \
--make-bed \
--out output/quality_control/file2 \
--silent \
--missing;                 # --missing produces sample-based and variant-based missing data reports. If run with --within/--family, the variant-based report is stratified by cluster.
grep "variants" output/quality_control/file2.log

## 129381 variants loaded from output/quality_control/file1.bim.
## --geno: 2271 variants removed due to missing genotype data.
## 127110 variants remaining after main filters.

indmiss <-
  read.delim(
    file   = here(
      "output", "quality_control", "file2.smiss"
    ),
    header = TRUE
  )

snpmiss <-
  read.delim(
    file   = here(
      "output", "quality_control", "file2.vmiss"
    ),
    header = TRUE
  )

ggplot(
  indmiss,
  aes(
    x = F_MISS
  )
) +
  geom_histogram(
    color            = "black",
    fill             = "#B6FAD7",
    bins             = 6
  ) +
  geom_text(
    stat             = "bin",
    aes(
    label = after_stat(count)
  ),
    vjust            = -0.5,
    color            = "purple",
    size             = 3,
    bins             = 6
  ) +
  geom_vline(
    aes(
    xintercept = mean(F_MISS)),
    color            = "orange",
    linetype         = "dotted",
    linewidth        = .5
  ) +
  annotate(
    "text",
    x = mean(indmiss$F_MISS),
    y = 30,
    label = paste0(
      "Mean \n",
      scales::percent(mean(indmiss$F_MISS),
        accuracy = 0.01
      )
    ),
    size = 3,
    color = "orange",
    hjust = -.1
  ) +
  labs(
    x                = "Individual Missingness (%)",
    y                = "Frequency (n)"
  ) +
  my_theme() +
  scale_x_continuous(
    labels           = scales::percent,
    n.breaks         = 6
  )

ggsave(
  here(
    "output", "quality_control", "figures" , "individual_missingness.pdf"
  ),
  width              = 7,
  height             = 5,
  units              = "in"
)

ggplot(
  snpmiss,
  aes(
    x = F_MISS
  )
) +
  geom_histogram(
    color = "black",
    fill = "#B6FAD7",
    bins = 6
  ) +
  stat_bin(
    geom = "text",
    aes(
      label = format(
        after_stat(count),
        big.mark = ",",
        scientific = FALSE
      )
    ),
    vjust = -0.5,
    color = "purple",
    size = 2,
    bins = 6
  ) +
  geom_vline(
    aes(
      xintercept = mean(F_MISS)
    ),
    color = "orange",
    linetype = "dotted",
    linewidth = 0.5
  ) +
  annotate(
    "text",
    x = mean(snpmiss$F_MISS),
    y = 16000,
    label = paste0(
      "Mean \n",
      scales::percent(mean(snpmiss$F_MISS),
        accuracy = 0.01
      )
    ),
    size = 3,
    color = "orange",
    vjust = -.2
  ) +
  labs(
    x = "Variant Missingness (%)",
    y = "Frequency (n)"
  ) +
  scale_x_continuous(
    labels = scales::percent,
    n.breaks = 6
  ) +
  scale_y_continuous(
    labels = scales::label_comma(),
    n.breaks = 5
  ) +
  my_theme()

ggsave(
  here(
    "output", "quality_control", "figures", "SNPs_missingness.pdf"
  ),
  width  = 7,
  height = 5,
  units  = "in"
)

Individuals with >20% missing genotypes were removed.

plink2 \
--allow-extra-chr \
--bfile output/quality_control/file2 \
--mind 0.2               `# here we set the individual missingness threshold of 20%`\
--make-bed \
--out output/quality_control/file3 \
--silent;
grep "samples\|variants" output/quality_control/file3.log

## 61 samples (25 females, 26 males, 10 ambiguous; 61 founders) loaded from
## 127110 variants loaded from output/quality_control/file2.bim.
## 0 samples removed due to missing genotype data (--mind).
## 61 samples (25 females, 26 males, 10 ambiguous; 61 founders) remaining after

No individuals were removed due to missingness.

3.2 Minor allele frequency

Allele frequency distribution was estimated with PLINK2.

plink2 \
--allow-extra-chr \
--bfile output/quality_control/file3 \
--freq \
--out output/quality_control/MAF_check \
--silent

maf_freq <-
  read.delim(
    here(
      "output", "quality_control", "MAF_check.afreq"
    ),
    header = TRUE
  )

ggplot(
  maf_freq,
  aes(ALT_FREQS)
) +
  geom_histogram(
    colour = "black",
    fill = "#C4F3F5",
    bins = 40
  ) +
  labs(
    x = "Minor Allele Frequency (MAF)",
    y = "Frequency (n)",
    caption = "Red line at MAF 10% threshold."
  ) +
  annotate(
    "text",
    x = .1,
    y = 8000,
    label = "15,890 SNPs",
    size = 3,
    color = "blue",
    vjust = -.2
  ) +
  geom_vline(xintercept = 0.1, color = "red") +
  my_theme() +
  theme(plot.caption = element_text(face = "italic")) +
  scale_y_continuous(label = scales::number_format(big.mark = ",")) +
  scale_x_continuous(breaks = c(0, 0.1, 0.2, 0.4, 0.6, 0.8, 1))

ggsave(
  here(
    "output", "quality_control", "figures", "MAF.pdf"
  ),
  width  = 5,
  height = 4,
  units  = "in"
)

Now we apply the MAF filter.

plink2 \
--allow-extra-chr \
--bfile output/quality_control/file3 \
--maf 0.01 \
--make-bed \
--out output/quality_control/file4 \
--silent;
grep "variants" output/quality_control/file4.log

## 127110 variants loaded from output/quality_control/file3.bim.
## 15890 variants removed due to allele frequency threshold(s)
## 111220 variants remaining after main filters.

HWE was tested per population following Marees et al. (2018).

mkdir -p output/quality_control/hardy;
cp output/quality_control/file4.* output/quality_control/hardy/

3.3 HWE test

HWE was tested per population with a p-value threshold of 1e-6. SNP missingness was re-applied within each population to stabilize HWE estimates.

for fam in $(awk '{print $1}' output/quality_control/hardy/file4.fam | sort | uniq);
do
echo $fam | \
plink2 \
--allow-extra-chr \
--silent \
--keep-allele-order \
--bfile output/quality_control/hardy/file4 \
--keep-fam /dev/stdin \
--make-bed \
--out output/quality_control/hardy/$fam \
--hwe 0.000001 \
--geno 0.1 \
--write-snplist; \
done

SNPs passing HWE in all populations were retained.

cat output/quality_control/hardy/*.snplist | awk '!a[$0]++' > output/quality_control/passed_hwe.txt;
wc -l output/quality_control/passed_hwe.txt

## 110353 output/quality_control/passed_hwe.txt

How many variants we had before

cat output/quality_control/file4.bim | awk '{print $2}'| awk '!a[$0]++' | wc -l

## 111220

Variants not passing HWE test

111220 - 110353

## [1] 867

Remove the SNPs that failed

plink2 \
--allow-extra-chr \
--bfile output/quality_control/file4 \
--extract output/quality_control/passed_hwe.txt \
--make-bed \
--out output/quality_control/file4a \
--silent;
grep "variants" output/quality_control/file4a.log

## 111220 variants loaded from output/quality_control/file4.bim.
## --extract: 110353 variants remaining.
## 110353 variants remaining after main filters.

3.4 LD pruning

plink2 \
--allow-extra-chr \
--bfile output/quality_control/file4a \
--indep-pairwise 5 1 0.1 \
--out output/quality_control/indepSNP \
--silent;
grep 'pairwise\|variants\|samples' output/quality_control/indepSNP.log

##   --indep-pairwise 5 1 0.1
## 61 samples (25 females, 26 males, 10 ambiguous; 61 founders) loaded from
## 110353 variants loaded from output/quality_control/file4a.bim.
## --indep-pairwise (11 compute threads): 68504/110353 variants removed.

3.6 Heterozygosity

plink2 \
--allow-extra-chr \
--bfile output/quality_control/file4a \
--extract output/quality_control/indepSNP.prune.in \
--het \
--out output/quality_control/R_check \
--silent;
grep 'variants' output/quality_control/R_check.log

## 110353 variants loaded from output/quality_control/file4a.bim.
## --extract: 41849 variants remaining.
## 41849 variants remaining after main filters.

het <- read.delim(
  here(
    "output", "quality_control", "R_check.het"
  ),
  head = TRUE
)
colnames(het)

## [1] "X.FID"  "IID"    "O.HOM." "E.HOM." "OBS_CT" "F"

Heterozygosity rate was calculated per individual. Individuals deviating more than 4 SD from the mean were flagged.

het$HET_RATE <- (het$"OBS_CT" - het$"O.HOM") / het$"OBS_CT"

het_fail <-
  subset(
    het, (het$HET_RATE < mean(
      het$HET_RATE
    ) - 4 * sd(
      het$HET_RATE
    )) |
      (het$HET_RATE > mean(
        het$HET_RATE
      ) + 4 * sd(
        het$HET_RATE
      ))
  )

het_fail$HET_DST <-
  (het_fail$HET_RATE - mean(
    het$HET_RATE
  )) / sd(
    het$HET_RATE
  )

write.table(
  het_fail,
  here(
    "output", "quality_control", "fail-het-qc.txt"
  ),
  row.names = FALSE
)

ggplot(
  het,
  aes(
    HET_RATE
  )
) +
  geom_histogram(
    colour           = "black",
    fill             = "#CDFAF8",
    bins             = 40
  ) +
  labs(
    x                = "Heterozygosity Rate",
    y                = "Number of Individuals"
  ) +
  geom_vline(
    aes(
      xintercept     = mean(
        HET_RATE
      )
    ),
    col              = "#F2C46F",
    linewidth        = 1.5
  ) +
  geom_vline(
    aes(
      xintercept     = mean(
        HET_RATE
      ) + 4 * sd(
        HET_RATE
      )
    ),
    col              = "#BFB9B9",
    linewidth        = 1
  ) +
  geom_vline(
    aes(
      xintercept     = mean(
        HET_RATE
      ) - 4 * sd(
        HET_RATE
      )
    ),
    col              = "#BFB9B9",
    linewidth        = 1
  ) +
  my_theme() +
  scale_y_continuous(
    labels           = comma
  )

ggsave(
  here(
    "output", "quality_control", "figures", "Heterozygosity.pdf"
  ),
  width  = 5,
  height = 4,
  units  = "in"
)

The orange line indicates the mean; gray lines indicate +/- 4 SD. Individuals identified as outliers were listed for removal.

sed 's/"// g' output/quality_control/fail-het-qc.txt | awk '{print$1, $2}'> output/quality_control/het_fail_ind.txt;
echo 'How many mosquitoes we need to remove from our data set:';
cat output/quality_control/het_fail_ind.txt | tail -n +2 | wc -l;
echo 'Which mosquitoes we have to remove:';
tail -n +2 output/quality_control/het_fail_ind.txt

## How many mosquitoes we need to remove from our data set:
## 1
## Which mosquitoes we have to remove:
## AUT 399

Two individuals were identified as heterozygosity outliers and removed from the dataset.

Heterozygosity outliers were removed:

plink2 \
--allow-extra-chr \
--bfile output/quality_control/file4a \
--remove output/quality_control/het_fail_ind.txt \
--make-bed \
--out output/quality_control/file5 \
--silent;
grep 'variants\|samples' output/quality_control/file5.log

## 61 samples (25 females, 26 males, 10 ambiguous; 61 founders) loaded from
## 110353 variants loaded from output/quality_control/file4a.bim.
## --remove: 60 samples remaining.
## 60 samples (24 females, 26 males, 10 ambiguous; 60 founders) remaining after

3.7 Relatedness

Cryptic relatedness was assessed with KING kinship coefficients (Manichaikul et al. 2010). KING coefficients are scaled such that duplicate samples have kinship 0.5, first-degree relations ~0.25, and second-degree relations ~0.125. A cutoff of 0.354 (geometric mean of 0.5 and 0.25) was used to screen for duplicates.

plink2 \
--allow-extra-chr \
--bfile output/quality_control/file5 \
--extract output/quality_control/indepSNP.prune.in \
--make-king-table rel-check \
--king-table-filter 0.354 \
--out output/quality_control/file5a \
--silent;
grep 'variants\|samples' output/quality_control/file5a.log

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

Individuals flagged by the KING relatedness filter:

head output/quality_control/file5a.kin0

## #FID1    IID1    FID2    IID2    NSNP    HETHET  IBS0    KINSHIP

No related individuals were found.

plink2 \
--allow-extra-chr \
--bfile output/quality_control/file5 \
--extract output/quality_control/indepSNP.prune.in \
--make-king triangle bin \
--out output/quality_control/file6 \
--silent;
grep 'variants\|samples' output/quality_control/file6.log

## 60 samples (24 females, 26 males, 10 ambiguous; 60 founders) loaded from
## 110353 variants loaded from output/quality_control/file5.bim.
## --extract: 41849 variants remaining.
## 41849 variants remaining after main filters.
## --make-king pass 1/1: Scanning for rare variants... done.
## 1245 variants handled by initial scan (40604 remaining).
## --make-king: 41849 variants processed.

PLINK2 removed related individuals while maximizing sample retention.

plink2 \
--allow-extra-chr \
--bfile output/quality_control/file5 \
--king-cutoff output/quality_control/file6 0.354 \
--make-bed \
--out output/quality_control/file7 \
--silent;
grep 'samples\|variants\|remaining' output/quality_control/file7.log

## 60 samples (24 females, 26 males, 10 ambiguous; 60 founders) loaded from
## 110353 variants loaded from output/quality_control/file5.bim.
## output/quality_control/file7.king.cutoff.out.id , and 60 remaining sample IDs

head -n 10 output/quality_control/file7.king.cutoff.out.id

## #FID IID

No individuals were removed due to relatedness.

3.8 Quick PCA with Plink using the LD pruned data

plink2 \
--allow-extra-chr \
--bfile output/quality_control/file7 \
--pca allele-wts \
--freq \
--out output/quality_control/pca_pops \
--silent;
grep 'samples\|variants' output/quality_control/pca_pops.log

## 60 samples (24 females, 26 males, 10 ambiguous; 60 founders) loaded from
## 110353 variants loaded from output/quality_control/file7.bim.

pca <- read.delim(
  here(
    "output", "quality_control", "pca_pops.eigenvec"
  ),
  head = TRUE
)

head(pca)

##   X.FID IID        PC1        PC2        PC3         PC4        PC5        PC6
## 1   MAN 306 -0.0913239 -0.0883668 -0.1575640  0.03060270  0.0736423 -0.0966122
## 2   MAN 307 -0.1165830 -0.1604110 -0.1928870  0.00426111 -0.0392307 -0.1805130
## 3   MAN 308 -0.1092400 -0.1071860 -0.1721780 -0.00431477  0.0581054 -0.4130700
## 4   MAN 309 -0.1105640 -0.1247490 -0.1573170 -0.00275033  0.0039986 -0.0441429
## 5   MAN 310 -0.0691896 -0.0441948 -0.0806318  0.19202600 -0.0199620  0.0510735
## 6   MAN 311 -0.0899829 -0.1201360 -0.1003930  0.01009320 -0.0576055  0.0219653
##            PC7       PC8        PC9       PC10
## 1  0.019931800 -0.220985  0.1620490  0.0959124
## 2 -0.000233315  0.233686 -0.1521200 -0.0316382
## 3  0.019968600  0.168698 -0.1586580  0.0386871
## 4  0.010466500 -0.225898  0.1217150  0.0941331
## 5 -0.000256229  0.140714  0.0297651 -0.2804500
## 6 -0.023091400 -0.230319 -0.0835833 -0.2497100

# Shapes
N = 100
M = 1000
good.shapes = c(1:25, 33:127)

# Colors
palette1 <- brewer.pal(12, "Paired")
palette2 <- brewer.pal(11, "Set3")
palette23 <- c(palette1, palette2)

# Recode population labels to manuscript names
pca$X.FID <- recode(pca$X.FID,
  "AUT" = "AUTO", "MAN" = "NON-AUTO", "NEW" = "NON-AUTO-FIELD")

# Compute the count
family_count <- pca |>
  group_by(X.FID) |>
  summarize(count = n())

# Merge the count back to the main data
df4 <- pca |>
  left_join(family_count, by = "X.FID")

# Create a custom label for the legend
df4$family_label <-
  paste(df4$X.FID, " (", df4$count, ")", sep = "")

# Define the color and shape manually
palette23 <- c(palette1, palette2)
colors <- setNames(palette23, unique(df4$family_label))
shapes <-
  setNames(good.shapes[c(1:25, 58:67)], unique(df4$family_label))


# Compute the center of ellipses for each continent
ellipse_centers <- df4 |>
  group_by(X.FID) |>
  summarise(PC1_center = mean(PC1), PC2_center = mean(PC2))

# Create the plot
ggplot(df4, aes(PC1, PC2)) +
  geom_point(aes(shape = family_label, color = family_label)) +
  stat_ellipse(
    aes(fill = X.FID, group = X.FID),
    geom = "polygon",
    alpha = 0.2,
    level = 0.8,
    segments = 40,
    color = "transparent",
    show.legend = FALSE
  ) +
  geom_text_repel(
    data = ellipse_centers,
    aes(x = PC1_center, y = PC2_center, label = X.FID),
    color = c("black", "red", "blue")
  ) +
  xlab("PC1 (7.94% Variance)") +
  ylab("PC2 (4.09% Variance)") +
  guides(
  color = guide_legend(title = "Family", order = 1, ncol = 3),
  shape = guide_legend(title = "Family", order = 1, ncol = 3),
  fill = guide_legend(title = "Family", order = 1, ncol = 3)
  ) +
  scale_fill_manual(values = c("gray", "red", "lightblue"), labels = df4$family_label) +
  scale_color_manual(values = c("black", "red", "blue")) +
  scale_shape_manual(values = shapes) +
  theme_classic() +
  theme(
    plot.caption = element_text(face = "italic"),
    legend.position = "top",
    legend.direction = "horizontal",
    legend.justification = "center",
    legend.box.just = "center",
    legend.box.background = element_blank(),
    plot.margin = margin(5.5, 30, 5.5, 5.5, "points"),
    legend.margin = margin(10, 10, 10, 10)
)

ggsave(
  here(
    "output", "quality_control", "figures", "PCA_plink.pdf"
  ),
  width  = 7,
  height = 5,
  units  = "in"
)

objects <- ls()
save(list = objects, file = here(
    "output", "quality_control", "figures", "pca.RData"
  ))

4. Create chromosomal scale

Import the .bim file with the SNPs to create a new chromosomal scale.

snps <-
  read_delim(
    here(
      "output", "quality_control", "file7.bim"
    ),
    col_names      = FALSE,
    show_col_types = FALSE,
    col_types      = "ccidcc"
  )

colnames(snps) <- c("Scaffold", "SNP", "Cm", "Position", "Allele1", "Allele2")

head(snps)

## # A tibble: 6 x 6
##   Scaffold SNP             Cm Position Allele1 Allele2
##   <chr>    <chr>        <int>    <dbl> <chr>   <chr>  
## 1 1.1      AX-581444870     0    97856 C       T      
## 2 1.1      AX-583035083     0   305518 A       G      
## 3 1.1      AX-583035102     0   308124 A       G      
## 4 1.1      AX-583033342     0   315059 C       G      
## 5 1.1      AX-583035163     0   315386 A       G      
## 6 1.1      AX-583033356     0   315674 C       T

SNP data were separated by chromosome.

chr_snps_list <- list()
for (chr in c("1", "2", "3")) {
  chr_snps_list[[chr]] <- snps |>
    filter(str_detect(Scaffold, paste0("^", chr, "\\."))) |>
    as_tibble()
}

The reference genome was indexed with samtools.

# index the genome
samtools faidx data/genome/albo.fasta.gz

Scaffold sizes were extracted from the .fai index file.

awk 'BEGIN{FS=" "; OFS="\t"} {print $1, $2}' data/genome/albo.fasta.gz.fai > data/genome/scaffold_sizes.txt;
head data/genome/scaffold_sizes.txt

## 1.1  351198
## 1.10 11939576
## 1.100    3389100
## 1.101    470438
## 1.102    2525157
## 1.103    150026
## 1.104    51905
## 1.105    662742
## 1.106    7592381
## 1.107    204538

Import the file with sizes of each scaffold.

sizes <-
  read_delim(
    here(
      "data", "genome", "scaffold_sizes.txt"
    ),
    col_names      = FALSE,
    show_col_types = FALSE,
    col_types      = "cd"
  )

colnames(sizes) <- c("Scaffold", "Size")

sizes <-
  sizes |>
  mutate(
    Chromosome = case_when(
      startsWith(Scaffold, "1") ~ "1",
      startsWith(Scaffold, "2") ~ "2",
      startsWith(Scaffold, "3") ~ "3"
    )
  ) |>
  arrange(Scaffold)

head(sizes)

## # A tibble: 6 x 3
##   Scaffold     Size Chromosome
##   <chr>       <dbl> <chr>     
## 1 1.1        351198 1         
## 2 1.10     11939576 1         
## 3 1.100     3389100 1         
## 4 1.101      470438 1         
## 5 1.102     2525157 1         
## 6 1.103      150026 1

Create new scale. Get the scaffolds for each chromosome.

chr_scaffolds_list <- list()
for (chr in c("1", "2", "3")) {
  chr_scaffolds_list[[chr]] <- sizes |>
    filter(str_detect(Scaffold, paste0("^", chr))) |>
    as_tibble()
}

Create a scale for each chromosome.

for (chr in c("1", "2", "3")) {
  chr_scaffolds_list[[chr]]$overall_size_before_bp <- 0
  for (i in 2:nrow(chr_scaffolds_list[[chr]])) {
    chr_scaffolds_list[[chr]]$overall_size_before_bp[i] <-
      chr_scaffolds_list[[chr]]$overall_size_before_bp[i - 1] +
      chr_scaffolds_list[[chr]]$Size[i - 1]
  }
}

Merge the data frames scaffolds and SNPs.

chr_scale_list <- list()
for (chr in c("1", "2", "3")) {
  chr_scale_list[[chr]] <- chr_snps_list[[chr]] |>
    left_join(chr_scaffolds_list[[chr]], by = "Scaffold") |>
    na.omit() |>
    mutate(
      midPos_fullseq = as.numeric(Position) + as.numeric(overall_size_before_bp)
    )
}

Merge all chromosome scales.

chroms <- bind_rows(chr_scale_list) |>
  dplyr::select(
    Chromosome, SNP, Cm, midPos_fullseq, Allele1, Allele2
  )
head(chroms)

## # A tibble: 6 x 6
##   Chromosome SNP             Cm midPos_fullseq Allele1 Allele2
##   <chr>      <chr>        <int>          <dbl> <chr>   <chr>  
## 1 1          AX-581444870     0          97856 C       T      
## 2 1          AX-583035083     0         305518 A       G      
## 3 1          AX-583035102     0         308124 A       G      
## 4 1          AX-583033342     0         315059 C       G      
## 5 1          AX-583035163     0         315386 A       G      
## 6 1          AX-583033356     0         315674 C       T

write.table(
  chroms,
  file      = here(
    "output", "quality_control", "file7B.bim"
  ),
  sep       = "\t",
  row.names = FALSE,
  col.names = FALSE,
  quote     = FALSE
)

Rename the .bim files

mv output/quality_control/file7.bim output/quality_control/file7_backup.bim;
mv output/quality_control/file7B.bim output/quality_control/file7.bim

Validation check with PLINK2 to confirm variant ordering:

plink2 \
--bfile output/quality_control/file7 \
--make-bed \
--out output/quality_control/test01;
rm output/quality_control/test01.*

## PLINK v2.0.0-a.6.9LM 64-bit Intel (29 Jan 2025)    cog-genomics.org/plink/2.0/
## (C) 2005-2025 Shaun Purcell, Christopher Chang   GNU General Public License v3
## Logging to output/quality_control/test01.log.
## Options in effect:
##   --bfile output/quality_control/file7
##   --make-bed
##   --out output/quality_control/test01
## 
## Start time: Sun Apr  5 23:08:42 2026
## 32011 MiB RAM detected, ~26062 available; reserving 16005 MiB for main
## workspace.
## Using up to 12 threads (change this with --threads).
## 60 samples (24 females, 26 males, 10 ambiguous; 60 founders) loaded from
## output/quality_control/file7.fam.
## 110353 variants loaded from output/quality_control/file7.bim.
## 1 binary phenotype loaded (28 cases, 22 controls).
## Writing output/quality_control/test01.fam ... done.
## Writing output/quality_control/test01.bim ... done.
## Writing output/quality_control/test01.bed ... 0%59%done.
## End time: Sun Apr  5 23:08:42 2026

No warnings from PLINK2.

5. Plot SNP density

After quality control with 110,353 SNPs

snp_den_qc <- import_bim(
  here(
    "output", "quality_control", "file7.bim"
  )
)

Make plot of the SNP density

snp_den_qc |>
  rename(
    Chromosome = 1
  ) |>
  mutate(
    Position           = as.numeric(
      Position
    )
  ) |>
  ggplot(
    aes(
      x                = Position
    )
  ) +
  geom_histogram(
    aes(
      y                = after_stat(
        count
      )
    ),
    binwidth           = 1e6
  ) +
  facet_wrap(
    vars(
      Chromosome
    ),
    scales             = "free_x"
  ) +
  labs(
    title              = "SNP Density after QC",
    x                  = expression(
      "Position in the genome (Mb)"
    ),
    y                  = expression(
      "Number of SNPs"
    )
  ) +
  scale_x_continuous(
    labels             = function(x) {
      format(
        x / 1e6,
        big.mark       = ",",
        scientific     = FALSE
      )
    }
  ) +
  geom_density(
  aes(
    y = 1e6 * after_stat(count)
  ),
  color = "red",
  linewidth = .75,
  alpha = .4,
  fill = "pink"
  ) +
  theme_minimal(
    base_size          = 12
  ) +
  theme(
    panel.grid.major   = element_line(
      linetype         = "dashed",
      linewidth        = 0.2,
      color            = "#fabbe2"
    ),
    panel.grid.minor   = element_line(
      linetype         = "dashed",
      linewidth        = 0.2,
      color            = "#fabbe2"
    ),
    panel.spacing      = unit(0.5, "lines"),
    strip.text         = element_text(
      face             = "bold", hjust = .5
    ),
    strip.background.x = element_rect(
      color            = "gray"
    ),
    plot.margin        = margin(10, 10, 10, 10),
    axis.text          = element_text(size = 12),
    axis.title         = element_text(size = 14)
  )

ggsave(
  here(
    "output", "quality_control","figures", "snp_density_after_qc.pdf"
  ),
  width  = 10,
  height = 6,
  units  = "in"
)

SNPs per chromosome

snps_per_chrm <-
  snp_den_qc |>
  count(
    Scaffold) |>
  rename(
    Chromosome = 1,
    "SNPs (N) " = 2
  )

knitr::kable(snps_per_chrm, caption = "SNPs per chromosome after quality control")

SNPs per chromosome after quality control
Chromosome	SNPs (N)
1	25105
2	45902
3	39346

Mean SNP density per 1 Mb window by chromosome:

albo_den <-
  snp_den_qc |>
  dplyr::select(
    Scaffold, Position
  ) |>
  group_by(
    Scaffold,
    windows               = cut_width(
      Position,
      width               = 1e6,
      boundary            = 0
    )
  ) |>
  summarise(
    n                     = n(),
    .groups               = "keep"
  ) |>
  group_by(
    Scaffold
  ) |>
  summarise(
    mean                  = mean(n),
    n                     = n(),
    .groups               = "keep"
  ) |>
  rename(
    Chromosome            = 1,
    "SNPs per 1Mb window" = 2,
    "Number of windows"   = 3
  )
knitr::kable(albo_den, caption = "SNPs per 1Mb window")

SNPs per 1Mb window
Chromosome	SNPs per 1Mb window	Number of windows
1	67.85135	370
2	79.00516	581
3	80.46217	489

after_qc <-
  snps_per_chrm |>
  left_join(
    albo_den,
    by = "Chromosome"
  )
knitr::kable(after_qc, caption = "QC summary per chromosome")

QC summary per chromosome
Chromosome	SNPs (N)	SNPs per 1Mb window	Number of windows
1	25105	67.85135	370
2	45902	79.00516	581
3	39346	80.46217	489

6. Create data sets for downstream analysis

Population subsets were created for pairwise comparisons (see label key in Section 1).

echo "MAN
AUT" > output/quality_control/man_aut.txt;
echo "NEW
AUT" > output/quality_control/new_aut.txt;
echo "MAN
NEW" > output/quality_control/man_new.txt
echo "AUT" > output/quality_control/aut.txt

VCF and BED files were created for each pairwise comparison using LD-pruned data.

MAN vs AUT

plink2 \
--bfile output/quality_control/file7 \
--keep-fam output/quality_control/man_aut.txt \
--export vcf \
--extract output/quality_control/indepSNP.prune.in \
--make-bed \
--geno 0.2 \
--out output/outflank/man_aut \
--silent;
grep "samples\|variants" output/outflank/man_aut.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.
## --keep-fam: 38 samples remaining.
## 38 samples (13 females, 15 males, 10 ambiguous; 38 founders) remaining after
## --geno: 517 variants removed due to missing genotype data.
## 41332 variants remaining after main filters.

NEW vs AUT

plink2 \
--bfile output/quality_control/file7 \
--keep-fam output/quality_control/new_aut.txt \
--export vcf \
--extract output/quality_control/indepSNP.prune.in \
--make-bed \
--geno 0.2 \
--out output/outflank/new_aut \
--silent;
grep "samples\|variants" output/outflank/new_aut.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.
## --keep-fam: 50 samples remaining.
## 50 samples (24 females, 26 males; 50 founders) remaining after main filters.
## --geno: 153 variants removed due to missing genotype data.
## 41696 variants remaining after main filters.

MAN vs NEW

plink2 \
--bfile output/quality_control/file7 \
--keep-fam output/quality_control/man_new.txt \
--export vcf \
--extract output/quality_control/indepSNP.prune.in \
--make-bed \
--geno 0.2 \
--out output/outflank/man_new \
--silent;
grep "samples\|variants" output/outflank/man_new.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.
## --keep-fam: 32 samples remaining.
## 32 samples (11 females, 11 males, 10 ambiguous; 32 founders) remaining after
## --geno: 643 variants removed due to missing genotype data.
## 41206 variants remaining after main filters.

AUT

plink2 \
--bfile output/quality_control/file7 \
--keep-fam output/quality_control/aut.txt \
--export vcf \
--extract output/quality_control/indepSNP.prune.in \
--make-bed \
--geno 0.2 \
--out output/outflank/aut \
--silent;
grep "samples\|variants" output/outflank/aut.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.
## --keep-fam: 28 samples remaining.
## 28 samples (13 females, 15 males; 28 founders) remaining after main filters.
## --geno: 844 variants removed due to missing genotype data.
## 41005 variants remaining after main filters.