-
Notifications
You must be signed in to change notification settings - Fork 0
/
genetic_diversity_visualization_10_29_20.R
104 lines (70 loc) · 3.1 KB
/
genetic_diversity_visualization_10_29_20.R
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
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
###############################################
############### Libraries #####################
###############################################
library(adegenet)
library(poppr)
library(ggplot2)
library(poppr)
library(hierfstat)
library(car)
library(PMCMR)
library(Demerelate)
library(diveRsity)
###############################################
############### Loading Files #################
###############################################
setwd("G:\\My Drive\\Masters")
master_gendiv <- read.csv("master_gendiv.csv")
##set as factor
master_gendiv$Accession <- as.factor(master_gendiv$Accession)
##create new rows
ne_matrix <- matrix(ncol = 2, nrow = length(master_gendiv[,1]))
##create new accession names
ne_matrix <- cbind(paste0(master_gendiv$Species,"_",master_gendiv$Accession), as.numeric(master_gendiv$Ne))
##na matrix
na_matrix <- cbind(paste0(master_gendiv$Species,master_gendiv$Accession), master_gendiv$Na)
##he matrix
he_matrix <- cbind(paste0(master_gendiv$Species,master_gendiv$Accession), master_gendiv$He)
##
ho_matrix <- cbind(paste0(master_gendiv$Species,master_gendiv$Accession), master_gendiv$Ho)
#F matrix
f_matrix <- cbind(paste0(master_gendiv$Species,master_gendiv$Accession), master_gendiv$F)
##accession list
accession_list <- unique(ne_matrix[,1])
##make ne matrix results
ne_df <- matrix(nrow = length(accession_list), ncol = 3)
##colnames(gendiv_df) <- c("Ne","Na","He","Ho","F")
##write loop to calculate means
for(a in 1:length(accession_list)){
ne_df[a,3] <- mean(as.numeric(ne_matrix[ne_matrix[,1] == paste0(accession_list[[a]]),][,2]))
}
# gendiv_df[a,2] <- mean(as.numeric(na_matrix[na_matrix[,1] == paste0(accession_list[[a]]),][,2]))
# gendiv_df[a,3] <- mean(as.numeric(he_matrix[he_matrix[,1] == paste0(accession_list[[a]]),][,2]))
# gendiv_df[a,4] <- mean(as.numeric(ho_matrix[ho_matrix[,1] == paste0(accession_list[[a]]),][,2]))
# gendiv_df[a,5] <- mean(as.numeric(f_matrix[f_matrix[,1] == paste0(accession_list[[a]]),][,2]))
#}
##reformat loop
##write to data frame
##write loop to calculate means
for(a in 1:length(accession_list)){
ne_df[a,1] <- mean(as.numeric(ne_matrix[ne_matrix[,1] == paste0(accession_list[[a]]),][,2]))
ne_df[a,2] <- (ne_df[a,1]) + std.error(as.numeric(ne_matrix[ne_matrix[,1] == paste0(accession_list[[a]]),][,2]))
ne_df[a,3] <- (ne_df[a,1]) - std.error(as.numeric(ne_matrix[ne_matrix[,1] == paste0(accession_list[[a]]),][,2]))
}
##create df
ne_matrix <- cbind(paste0(master_gendiv$Species,master_gendiv$Accession), as.numeric(master_gendiv$Ne))
ne_df <- matrix(nrow = length(accession_list), ncol = 3)
##
ne_df <- data.frame(ne_df)
##now separate
ne_df[,1] <- gsub("\\_.*","",ne_df[,1])
ne_df[,2] <- gsub("^.*\\_","", ne_df[,1])
colnames(ne_df) <- c("Species", "Accession","NE")
ne_df[,2] <- as.factor(ne_df[,2])
plot(ne_df[,3])
###Plot
ggplot(data = ne_df, aes(x=Species, y=NE)) +
geom_bar(aes(fill = Accession), stat = "identity", position = "dodge") + xlab("Species") +
ylab("Effective Number of Alleles") + ggtitle("Effective Number of Alleles") +
scale_fill_manual(values = c("gray33", "gray48", "gray77")) +
geom_errorbar()