Compute Fst with PLINK - grouped data --within

# Fst analysis with pooled samples from many breeds
# Analysis done using PLINK and the --within option

# Clear workspace
rm(list = ls())
# Set working directory
setwd("d:/analysis/2020_GenomicsBootCamp_Demo/")
library(tidyverse)
library(qqman)

# Read in the data - ADAPTmap
goats <- read_tsv("ADAPTmap_genotypeTOP_20160222_full.fam", col_names = F)

# select goats from target countries and apply country label
prodLabels <- goats %>%
  mutate(productionType = case_when(
    substr(X1,1,3) == "SAA" ~ "Dairy",          # Saanen goat
    substr(X1,1,3) == "TOG" ~ "Dairy",          # Toggenburg goat
    substr(X1,1,3) == "RAN" ~ "Meat",           # Rangeland goat
    substr(X1,1,3) == "TED" ~ "Meat"            # Teddi goat
  )) %>%
  drop_na() %>%
  select(X1, X2, productionType) %>%
  write_delim("productionComparison.txt")

# Keep only these animals and do the quality control
system(str_c("plink --bfile ADAPTmap_genotypeTOP_20160222_full --chr-set 29 --autosome ",
             "--keep productionComparison.txt --nonfounders ",
             "--geno 0.1 --mind 0.1 ",
             "--make-bed --out afterQC"))

# do the Fst analysis between the "Dairy" and "Meat" breeds using --within
system(str_c("plink --bfile afterQC --chr-set 29 --autosome ",
             "--fst --within productionComparison.txt ",
             "--out productionTypeFst"))


# load and visualize Fst results

data <- read.delim("productionTypeFst.fst") %>%
  drop_na()

manhattan(data, chr = "CHR", bp = "POS", p = "FST", logp = F)