Bayescan_input.R 885 Bytes
Newer Older
eboulanger's avatar
eboulanger committed
1
2
3
4
5
#----- BayeScan -----#
## produce Bayescan input file

#!/usr/bin/env Rscript
args = commandArgs(trailingOnly=TRUE)
eboulanger's avatar
eboulanger committed
6
7
#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
8
9
10
11
12
13
14
15
16
17
18

start <- Sys.time()

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

# define arguments
inputVcf <- args[1]
speciesCode <- args[2]
eboulanger's avatar
eboulanger committed
19
outputBayeScan <- paste0(speciesCode,"_bayesc_input_3pop")
eboulanger's avatar
eboulanger committed
20
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
grp <- find.clusters(genind, max.n.clust = 40)
eboulanger's avatar
eboulanger committed
29
30
400
2
eboulanger's avatar
eboulanger committed
31
32
33
34
35
36
37
38
39
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