Advertisement
GenomicsBootCamp

Localize varLD signals

Oct 12th, 2021
337
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. # varLD signal follow-up
  2. # see the video on the Genomics Boot Camp YouTube channel
  3.  
  4. # Clear workspace and load packages
  5. rm(list = ls())
  6. library(tidyverse)
  7.  
  8. # Set the location of the working directory
  9. setwd("XXXXXXXXXXXXXXXXXXXX")
  10.  
  11. ####
  12. # Examine the Top regions in detail - see exact SNP positions
  13. ####
  14.  
  15. # input parameter for the chromosome and threshold to be filtered
  16. chrOfInterest <- 13
  17. thresholdOfInterest <- "0.01%" # permissible values: "1%", "0.1%", "0.01%"
  18.  
  19.  
  20. #Read in the standardized file
  21. # PATH and FILENAME taken from the script above
  22. fullStandardFile <- read_tsv(str_c(PATH, FILENAME, chrOfInterest, "_standardized.out"))
  23.  
  24. # determine the cutoff based on thresholdOfInterest
  25. limitValue <- case_when(
  26.   thresholdOfInterest == "1%" ~ varLD.threshold[2],
  27.   thresholdOfInterest == "0.1%" ~ varLD.threshold[3],
  28.   thresholdOfInterest == "0.01%" ~ varLD.threshold[4])
  29.  
  30. # extract values above the limitValue
  31. regionsOfInterest <- fullStandardFile %>%
  32.   dplyr::filter(standardized_score >= limitValue)
  33. regionsOfInterest
  34.  
  35. # extract this region for both breeds and compare the LD
  36.  
  37. # LD table
  38. 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")
  39. 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")
  40.  
  41. # read LD matrix and plot as heatmap
  42. read_table("topRegionLd_table_breed1.ld") %>%
  43.   ggplot(., aes(SNP_A, SNP_B)) +
  44.   geom_tile(aes(fill = R2)) +
  45.   ggtitle("Charolais") +
  46.   theme_bw()
  47.  
  48.  
  49. # read LD matrix and plot as heatmap
  50. read_table("topRegionLd_table_breed2.ld") %>%
  51.   ggplot(., aes(SNP_A, SNP_B)) +
  52.   geom_tile(aes(fill = R2)) +
  53.   ggtitle("Holstein") +
  54.   theme_bw()
  55.  
  56.  
  57. ##########################################
  58.  
  59. # check on genotypes
  60. system("plink --bfile breed1_QC --cow --chr 13 --from-mb 47.15 --to-mb 47.35 --recode --out breed1LdRegion")
  61. system("plink --bfile breed2_QC --cow --chr 13 --from-mb 47.15 --to-mb 47.35 --recode --out breed2LdRegion")
  62.  
  63.  
  64. #LD matrix
  65. system("plink --bfile xbreed1_QC --cow --chr 13 --from-mb 47.15 --to-mb 47.35 --r2 square --out topRegionLd_breed1")
  66. system("plink --bfile xbreed2_QC --cow --chr 13 --from-mb 47.15 --to-mb 47.35 --r2 square --out topRegionLd_breed2")
  67.  
  68.  
  69.  
  70.  
  71.  
  72.  
  73.  
  74.  
  75.  
Advertisement
RAW Paste Data Copied
Advertisement