Bayescan_input.R 1.07 KB
Newer Older
eboulanger's avatar
eboulanger committed
1
2
#----- BayeScan -----#
## produce Bayescan input file
eboulanger's avatar
eboulanger committed
3
# setwd("~/Documents/project_SEACONNECT/seaConnect--dataPrep/02-Bayescan/")
eboulanger's avatar
eboulanger committed
4
5
6

#!/usr/bin/env Rscript
args = commandArgs(trailingOnly=TRUE)
eboulanger's avatar
eboulanger committed
7
8
#args <- c("../01-SNPfilters/01-Diplodus/dip_all_filtered.vcf","dip")
#args <- c("../01-SNPfilters/02-Mullus/mul_all_filtered.vcf","mul")
eboulanger's avatar
eboulanger committed
9
10
11
12
13
14
15
16
17
18
19

start <- Sys.time()

# load libraries
library(vcfR)
library(adegenet)
library(radiator)

# define arguments
inputVcf <- args[1]
speciesCode <- args[2]
eboulanger's avatar
eboulanger committed
20
outputBayeScan <- paste0(speciesCode,"_bayesc_input_2pop")
eboulanger's avatar
eboulanger committed
21
22
23
24
25
26
27
28
  
# import vcf file to genind
vcfFile <- read.vcfR(inputVcf, verbose = F)
!is.null(vcfFile)
genind <- vcfR2genind(vcfFile)
!is.null(genind)

# use DAPC groups as pop
eboulanger's avatar
eboulanger committed
29
grp <- find.clusters(genind, max.n.clust = 40, n.pca = round(nrow(genind$tab)/3))
eboulanger's avatar
eboulanger committed
30
2
eboulanger's avatar
eboulanger committed
31
32
grp$size 
pop(genind) <- grp$grp
eboulanger's avatar
eboulanger committed
33
pop(genind) <- paste0("pop", pop(genind))
eboulanger's avatar
eboulanger committed
34
35

# convert to bayescan format
eboulanger's avatar
eboulanger committed
36
37
tidy_genind <- tidy_genomic_data(genind)
write_bayescan(tidy_genind, pop.select = levels(pop(genind)), filename = outputBayeScan)
eboulanger's avatar
eboulanger committed
38
39
40
41

end <- Sys.time()
print(paste("run ended at", end))
end - start