PCAdapt.R 1.91 KB
Newer Older
eboulanger's avatar
eboulanger committed
1
2
3
4
5
#----- PCAdapt -----#
# vignette: https://cran.r-project.org/web/packages/pcadapt/vignettes/pcadapt.html

#!/usr/bin/env Rscript
args = commandArgs(trailingOnly=TRUE)
eboulanger's avatar
eboulanger committed
6
#args <- c("../01-SNPfilters/01-Diplodus/dip_all_filtered.vcf", "dip")
eboulanger's avatar
eboulanger committed
7
8
9
10
11

# load libraries
library(pcadapt)
library(vcfR)

eboulanger's avatar
eboulanger committed
12
# define arguments
eboulanger's avatar
eboulanger committed
13
path_to_file <- args[1]
eboulanger's avatar
eboulanger committed
14
15
16
speciesCode <- args[2]

# import vcf file to pcadapt format
eboulanger's avatar
eboulanger committed
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
45
46
47
48
49
50
51
52
53
54
dat <- read.pcadapt(path_to_file, type = "vcf")

# choose number K of principal components
x <- pcadapt(input = dat, K = 20)
# set plot output file name
pdf(file=paste0("RPlots_", args[2],".pdf"))
# scree plot of principal components
plot(x, option = "screeplot", K = 10)
# score plot of principal components
plot(x, option = "scores")
plot(x, option = "scores", i = 3, j = 4) # check PC axis 3 and 4

# compute outliers detection with K = 4
x <- pcadapt(dat, K = 4)
# Manhattan plot
#jpeg(filename = paste0(args[2],"_manhattan_pcadpt",".jpeg"))
plot(x, option = "manhattan")
#dev.off()
# QQ plot
plot(x , option = "qqplot", threshold = 0.05)
# histogram of p-values
#jpeg(filename = paste0(args[2],"_histP_pcadpt",".jpeg"))
hist(x$pvalues, xlab = "p-values", main = NULL, breaks = 50, col = "orange")
#dev.off()
# histogram of test statistic
plot(x, option = "stat.distribution")

# choose cutoff for outlier detection
# Bonferroni correction
padj <- p.adjust(x$pvalues, method="bonferroni")
alpha <- 0.05
outliers <- which(padj < alpha)
print(paste(length(outliers), "outliers detected"),quote=0)

# extract corresponding SNP names and positions
vcf <- read.vcfR(path_to_file, verbose=F)
loci <- as.data.frame(vcf@fix[,1:2])
outlier_loci <- loci[outliers,]
eboulanger's avatar
eboulanger committed
55
outlier_positions <- droplevels(outlier_loci$POS)
eboulanger's avatar
eboulanger committed
56
57

# output positions table
eboulanger's avatar
eboulanger committed
58
59
write.table(outlier_positions, file = paste0("outl_pos_pcadpt_",speciesCode,".txt"), sep = "\t", quote = F, row.names = F, col.names = F)
print(paste0("outlier positions table exported to outl_pos_pcadpt_",speciesCode), quote = 0)