-
Notifications
You must be signed in to change notification settings - Fork 0
/
PCA_Figure5.sh
72 lines (48 loc) · 2.02 KB
/
PCA_Figure5.sh
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
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
# PCA - Figure 5
ml R_packages/3.6.0
cd /home/vpeona/2021SatCorvides/Intermediate/RE2/PresAbs/Hybrids
CODE=/home/vpeona/2021SatCorvides/Code
Rscript --vanilla $CODE/pca_corvus_hybrids.R proportions.txt
cd /home/vpeona/2021SatCorvides/Intermediate/RE2/PresAbs
Rscript --vanilla $CODE/pca_all_species.R proportions.txt
## Figure 5a
#!/usr/bin/env Rscript
# Usage: Rscript --vanilla pca_all_species.R proportions.txt
args = commandArgs(trailingOnly=TRUE)
# test if there is at least one argument: if not, return an error
if (length(args) < 1) {
stop("3 arguments must be supplied", call.=FALSE)
} else if (length(args) == 1) {
propfile = args[1]
}
library(tidyr)
library(data.table)
library(ggbiplot)
data = fread(propfile, header = TRUE)
new = data[,c(1,2,4)]
new_spread = spread(new, Satellite, Proportion)
new_spread = as.data.frame(new_spread)
sat.pca <- prcomp(new_spread[,c(2:length(new_spread))], center = TRUE,scale. = TRUE)
sat.groups = c(rep("BOP", 4), rep("Corvus", 7), rep("BOP", 12))
plot = ggbiplot(sat.pca, ellipse=TRUE, labels=new_spread$Species, groups=sat.groups)
ggsave(plot = plot, filename = "PCA_all_species.pdf", units = "cm", width = 30, height = 20, dpi = 300, device= "pdf")
## Figure 5b
#!/usr/bin/env Rscript
# Usage: Rscript --vanilla pca_corvus_hybrids.R proportions.txt
args = commandArgs(trailingOnly=TRUE)
# test if there is at least one argument: if not, return an error
if (length(args) < 1) {
stop("3 arguments must be supplied", call.=FALSE)
} else if (length(args) == 1) {
propfile = args[1]
}
library(tidyr)
library(data.table)
library(ggbiplot)
data = fread(propfile, header = TRUE)
new = data[,c(1,2,4)]
new_spread = spread(new, Satellite, Proportion)
sat.pca <- prcomp(new_spread[,c(2:23)], center = TRUE,scale. = TRUE)
sat.groups = c(rep("corCon", 3), rep("corCor", 3), rep("Hybrids", 6))
plot = ggbiplot(sat.pca, ellipse=TRUE, labels=new_spread$Species, groups=sat.groups)
ggsave(plot = plot, filename = "PCA_corvus_hybrids_2.2.pdf", units = "cm", width = 30, height = 20, dpi = 300, device= "pdf")