Bayescan_input.R 887 Bytes
Newer Older
eboulanger's avatar
eboulanger committed
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
#----- BayeScan -----#
## produce Bayescan input file

#!/usr/bin/env Rscript
args = commandArgs(trailingOnly=TRUE)
#args <- c("../01-SNPfilters/01-Diplodus/diplodus_all_filtered.vcf","dip")
#args <- c("../01-SNPfilters/02-Mullus/mullus_all_filtered.vcf","mul")

start <- Sys.time()

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

# define arguments
inputVcf <- args[1]
speciesCode <- args[2]
outputBayeScan <- paste0(speciesCode,"_bayesc_input_2pop")
  
# 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
grp <- find.clusters(genind, max.n.clust = 40)
grp$size 
pop(genind) <- grp$grp

# convert to bayescan format
genomic_converter(genind, output="bayescan", filename = outputBayeScan)

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