-
Notifications
You must be signed in to change notification settings - Fork 0
/
RelAbun_BPlot.R
148 lines (121 loc) · 5.61 KB
/
RelAbun_BPlot.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
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
library(tidyverse)
library(reshape2)
library(hrbrthemes)
library(viridis)
Tax.sum <- function(OTU.Table, Tax.Table, Tax.lvl ){
z <- NULL
y <- NULL
for (i in 1:length(unique(Tax.Table[colnames(OTU.Table),Tax.lvl]))) {
if (length(OTU.Table[,which(Tax.Table[colnames(OTU.Table),Tax.lvl]==unique(Tax.Table[colnames(OTU.Table),Tax.lvl])[i])])!=length(rownames(OTU.Table))) {
z <- which(Tax.Table[colnames(OTU.Table),Tax.lvl]==unique(Tax.Table[colnames(OTU.Table),Tax.lvl])[i])
y <- cbind(y, apply(OTU.Table[,which(Tax.Table[colnames(OTU.Table),Tax.lvl]==unique(Tax.Table[colnames(OTU.Table),Tax.lvl])[i])], 1, function(x) sum(x)))
} else {
y <- cbind(y, OTU.Table[,which(Tax.Table[colnames(OTU.Table),Tax.lvl]==unique(Tax.Table[colnames(OTU.Table),Tax.lvl])[i])])
}
}
colnames(y) <- unique(Tax.Table[colnames(OTU.Table),Tax.lvl])
invisible((y))
}
rltv.Otu.Table <- function(x){
x.Data.rltv <- NULL
for (i in 1:dim(x)[1]) {
x.Data.rltv <- rbind(x.Data.rltv, x[i,]/apply(x, 1, function(x) sum(x))[i])
}
rownames(x.Data.rltv) <- rownames(x)
invisible(x.Data.rltv)
}
ShaASVs <- read.csv2("https://raw.githubusercontent.com/cmlglvz/datasets/master/Data/eAnalisis/ShaASVs.csv", header = TRUE, sep = ";", dec = ".", row.names = 1, skip = 0, fill = TRUE)
ShaTXs <- read.csv2("https://raw.githubusercontent.com/cmlglvz/datasets/master/Data/eAnalisis/ShaTXs.csv", header = TRUE, sep = ";", dec = ".", row.names = 1, skip = 0, fill = TRUE)
colnames(ShaASVs) <- ShaTXs$OTU
ShaRAs <- as.data.frame(t(ShaASVs))
ShaRAs <- mutate(ShaRAs,
Cha = rowSums(ShaRAs[1:12]),
Fla = rowSums(ShaRAs[13:24]),
Hu = rowSums(ShaRAs[25:36]),
Pc = rowSums(ShaRAs[37:41]),
Total = rowSums(ShaRAs[1:41])
)
sRAs <- ShaRAs[, c(42:45)]
sRAs <- rltv.Otu.Table(sRAs)
sRAs <- apply(sRAs, 2, function(x) x * 100) %>% as.data.frame()
sRAs <- as.data.frame(t(sRAs))
ChlTXs <- filter(ShaTXs, Division == "Chlorophyta")
Chl_I <- filter(ChlTXs, Class == "Mamiellophyceae")
Chl_II <- filter(ChlTXs, Class == "Chloropicophyceae" | Class == "Trebouxiophyceae" | Class == "Pyramimonadophyceae" | Class == "Chlorodendrophyceae")
fChl_I <- filter(Chl_I, Genus == "Micromonas")
sChl_I <- filter(Chl_I, Genus == "Ostreococcus" | Genus == "Bathycoccus" | Genus == "Crustomastigaceae_X")
fIRAs <- select(sRAs, fChl_I$OTU)
fIRAs <- as.data.frame(t(fIRAs))
afChlI <- mutate(fIRAs,
OTU = rownames(fIRAs),
.before = "Cha")
bfChlI <- gather(afChlI, key = "Site", value = "RAbund", -1)
fChlIgg <- ggplot(bfChlI, aes(x = OTU, y = RAbund, fill = Site)) +
geom_boxplot() +
scale_fill_viridis(discrete = TRUE, alpha = 1) +
geom_jitter(color = "grey", size = 0.4, alpha = 0.9) +
theme_ipsum() +
theme(
legend.position = "right",
plot.title = element_text(size = 11)
) +
ggtitle("Chlorophyta relative abundance") +
xlab("")
vln <- bfChlI %>%
ggplot(aes(fill = Site, y = RAbund, x = OTU)) +
geom_violin(position = "dodge", alpha = 0.8, outlier.colour = "transparent") +
scale_fill_viridis(discrete = TRUE, name = "") +
theme_ipsum() +
xlab("") +
ylab("Relative Abundance (%)")
vln
OchTXs <- filter(ShaTXs, Division == "Ochrophyta")
OchASVs <- select(ShaASVs, all_of(OchTXs$Seq))
colnames(OchASVs) <- OchTXs$OTU
OchASVs <- mutate(OchASVs, Site = as.factor(c(rep("Cha", 12), rep("Fla", 12), rep("Hu", 12), rep("Pc", 5))), .before = "ASV_18")
HapTXs <- filter(ShaTXs, Division == "Haptophyta")
HapASVs <- select(ShaASVs, all_of(HapTXs$Seq))
colnames(HapASVs) <- HapTXs$OTU
HapASVs <- mutate(HapASVs, Site = as.factor(c(rep("Cha", 12), rep("Fla", 12), rep("Hu", 12), rep("Pc", 5))), .before = "ASV_38")
CKTXs <- filter(ShaTXs, Division == "Cryptophyta" | Division == "Katablepharidophyta" | Division == "Cryptophyta:nucl")
CKASVs <- select(ShaASVs, all_of(CKTXs$Seq))
colnames(CKASVs) <- CKTXs$OTU
CKASVs <- mutate(CKASVs, Site = as.factor(c(rep("Cha", 12), rep("Fla", 12), rep("Hu", 12), rep("Pc", 5))), .before = "ASV_31")
bChl <- gather(aChl, key = "OTU", value = "Rel.Abun", -c(1, 2))
ggChl <- ggplot(bChl, aes(x = variable, y = value, fill = Site)) +
geom_boxplot() +
scale_fill_viridis(discrete = TRUE, alpha = 0.6) +
geom_jitter(color = "black", size = 0.4, alpha = 0.9) +
theme_ipsum() +
theme(
legend.position = "none",
plot.title = element_text(size = 11)
) +
ggtitle("Chlorophyta relative abundance") +
xlab("")
ggChl
gg <- ggplot(bfChlI, aes(x = Site, y = RAbund, fill = OTU)) +
geom_boxplot(colour = "black", position = position_dodge(0.5)) +
geom_vline(xintercept = c(1.5,2.5,3.5), colour = "grey85", size = 1.2) +
theme(legend.title = element_text(size = 12, face = "bold"),
legend.text = element_text(size = 10, face = "bold"), legend.position = "right",
axis.text.x = element_text(face = "bold",colour = "black", size = 12),
axis.text.y = element_text(face = "bold", size = 12, colour = "black"),
axis.title.y = element_text(face = "bold", size = 14, colour = "black"),
panel.background = element_blank(),
panel.border = element_rect(fill = NA, colour = "black"),
legend.key=element_blank()) +
labs(x= "", y = "Relative Abundance (%)", fill = "OTU") +
scale_fill_viridis(discrete = TRUE, alpha = 0.6, option = "D")
gg
vln <- ggplot(bChl, aes(x = variable, y = value, fill = Site)) +
geom_violin() +
scale_fill_viridis(discrete = TRUE, alpha = 1, option = "A") +
theme_ipsum() +
theme(
legend.position = "right",
plot.title = element_text(size = 11)
) +
ggtitle("Violin chart") +
xlab("")
vln