Advertisement
GenomicsBootCamp

Compute Fst with PLINK - grouped data --within

Jun 22nd, 2021
476
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. # Fst analysis with pooled samples from many breeds
  2. # Analysis done using PLINK and the --within option
  3.  
  4. # Clear workspace
  5. rm(list = ls())
  6. # Set working directory
  7. setwd("d:/analysis/2020_GenomicsBootCamp_Demo/")
  8. library(tidyverse)
  9. library(qqman)
  10.  
  11. # Read in the data - ADAPTmap
  12. goats <- read_tsv("ADAPTmap_genotypeTOP_20160222_full.fam", col_names = F)
  13.  
  14. # select goats from target countries and apply country label
  15. prodLabels <- goats %>%
  16.   mutate(productionType = case_when(
  17.     substr(X1,1,3) == "SAA" ~ "Dairy",          # Saanen goat
  18.     substr(X1,1,3) == "TOG" ~ "Dairy",          # Toggenburg goat  
  19.     substr(X1,1,3) == "RAN" ~ "Meat",           # Rangeland goat
  20.     substr(X1,1,3) == "TED" ~ "Meat"            # Teddi goat
  21.   )) %>%
  22.   drop_na() %>%
  23.   select(X1, X2, productionType) %>%
  24.   write_delim("productionComparison.txt")
  25.  
  26. # Keep only these animals and do the quality control
  27. system(str_c("plink --bfile ADAPTmap_genotypeTOP_20160222_full --chr-set 29 --autosome ",
  28.              "--keep productionComparison.txt --nonfounders ",
  29.              "--geno 0.1 --mind 0.1 ",
  30.              "--make-bed --out afterQC"))
  31.  
  32. # do the Fst analysis between the "Dairy" and "Meat" breeds using --within
  33. system(str_c("plink --bfile afterQC --chr-set 29 --autosome ",
  34.              "--fst --within productionComparison.txt ",
  35.              "--out productionTypeFst"))
  36.  
  37.  
  38. # load and visualize Fst results
  39.  
  40. data <- read.delim("productionTypeFst.fst") %>%
  41.   drop_na()
  42.  
  43. manhattan(data, chr = "CHR", bp = "POS", p = "FST", logp = F)
  44.  
  45.  
Advertisement
RAW Paste Data Copied
Advertisement