Localize varLD signals

# varLD signal follow-up
# see the video on the Genomics Boot Camp YouTube channel

# Clear workspace and load packages
rm(list = ls())
library(tidyverse)

# Set the location of the working directory
setwd("XXXXXXXXXXXXXXXXXXXX")

####
# Examine the Top regions in detail - see exact SNP positions
####

# input parameter for the chromosome and threshold to be filtered
chrOfInterest <- 13
thresholdOfInterest <- "0.01%" # permissible values: "1%", "0.1%", "0.01%"


#Read in the standardized file
# PATH and FILENAME taken from the script above
fullStandardFile <- read_tsv(str_c(PATH, FILENAME, chrOfInterest, "_standardized.out"))

# determine the cutoff based on thresholdOfInterest
limitValue <- case_when(
  thresholdOfInterest == "1%" ~ varLD.threshold[2],
  thresholdOfInterest == "0.1%" ~ varLD.threshold[3],
  thresholdOfInterest == "0.01%" ~ varLD.threshold[4])

# extract values above the limitValue
regionsOfInterest <- fullStandardFile %>%
  dplyr::filter(standardized_score >= limitValue)
regionsOfInterest

# extract this region for both breeds and compare the LD

# LD table
system("plink --bfile breed1_QC --cow --chr 13 --from-mb 47.15 --to-mb 47.35 --r2 --ld-window-r2 0 --out topRegionLd_table_breed1")
system("plink --bfile breed2_QC --cow --chr 13 --from-mb 47.15 --to-mb 47.35 --r2 --ld-window-r2 0 --out topRegionLd_table_breed2")

# read LD matrix and plot as heatmap
read_table("topRegionLd_table_breed1.ld") %>%
  ggplot(., aes(SNP_A, SNP_B)) +
  geom_tile(aes(fill = R2)) +
  ggtitle("Charolais") +
  theme_bw()


# read LD matrix and plot as heatmap
read_table("topRegionLd_table_breed2.ld") %>%
  ggplot(., aes(SNP_A, SNP_B)) +
  geom_tile(aes(fill = R2)) +
  ggtitle("Holstein") +
  theme_bw()


##########################################

# check on genotypes
system("plink --bfile breed1_QC --cow --chr 13 --from-mb 47.15 --to-mb 47.35 --recode --out breed1LdRegion")
system("plink --bfile breed2_QC --cow --chr 13 --from-mb 47.15 --to-mb 47.35 --recode --out breed2LdRegion")


#LD matrix
system("plink --bfile xbreed1_QC --cow --chr 13 --from-mb 47.15 --to-mb 47.35 --r2 square --out topRegionLd_breed1")
system("plink --bfile xbreed2_QC --cow --chr 13 --from-mb 47.15 --to-mb 47.35 --r2 square --out topRegionLd_breed2")