Bioinfo_scripts/run_ihs.R at main · CarolinaPB/Bioinfo_scripts · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
library(vcfR)
library(rehh)
library(data.table)
library(R.utils)
library(optparse)

option_list = list(
  make_option(c("-r", "--rds"), type="character", default=NULL,
              help="first rds to be used", metavar="character"),
  make_option(c("-p", "--pop"), type="character", default=NULL,
              help="name of population 1", metavar="character"),
);

opt_parser = OptionParser(option_list=option_list);
opt = parse_args(opt_parser);


print(paste0("Running iHS on ", opt$pop))


wgscan <- readRDS(file = opt$rds)


#iHS ANALYISIS
###########################################################################
# Chromosome list for your species
# Excluding chr29-not available for galGal6 ref genome
chrs <- c(seq(1,28), seq(30,33), "W", "Z")
hap_file = opt$vcf

# calculate genome-wide iHS values
print("calculating ihs")
palette("default")
wgscan.ihs <- ihh2ihs(wgscan)

print("saving ihs")
saveRDS(wgscan.ihs, file= paste0("iHS_", opt$pop, ".rds"))

pdf(paste0("iHS_", opt$pop,".pdf"))
manhattanplot(wgscan.ihs,
              pval = TRUE,
              threshold = 4,
              main = paste0("iHS population (", opt$pop, ")"))
dev.off()