forked from sydneyg/OTUvESV
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathNMDS_zotusvotus_16SandITS2_T3.R
289 lines (235 loc) · 13.8 KB
/
NMDS_zotusvotus_16SandITS2_T3.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
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
#Feb 19, 2018
#make beta diversity panels for fungi and bacteria separate OTU v ESV
#try it in base R and with ggplot
#Reset R's Brain
rm(list=ls())
#load libraries
library(plyr )
library(vegan)
library(ggplot2)
library(ggpubr)
library(plyr )
#set working directory
setwd("~/Dropbox/StatsandProgramming/16SElevationGradient/")
################################################################################################
######Fungi: avg.dist bray-curtis dissimilarity matrices - median and square root transformed for T1 T2 T3 trasplant
################################################################################################
#bray curtis
fungi_T3_bray_otu <- read.csv("data/DissMatricesforPRIMER/braycurtis_rarefied_dm_sqrt_ITS2_T3_otu.csv", row.names=1 )
fungi_T3_bray_zotu <- read.csv("data/DissMatricesforPRIMER/braycurtis_rarefied_dm_sqrt_ITS2_T3_zotu.csv", row.names=1 )
#check that row and column names all match
row.names(fungi_T3_bray_otu) == row.names(fungi_T3_bray_zotu)
row.names(fungi_T3_bray_otu) ==colnames(fungi_T3_bray_otu)
colnames(fungi_T3_bray_otu) ==colnames(fungi_T3_bray_zotu)
################################################################################################
######Bacteria: avg.dist bray-curtis dissimilarity matrices - median and square root transformed for T1 T2 T3 trasplant
################################################################################################
#bray-curtis
bac_T3_bray_otu <- read.csv("data/DissMatricesforPRIMER/braycurtis_rarefied_dm_sqrt_16S_T3_otu.csv", row.names=1 )
bac_T3_bray_zotu <- read.csv("data/DissMatricesforPRIMER/braycurtis_rarefied_dm_sqrt_16S_T3_zotu.csv", row.names=1 )
#check that row and column names all match
row.names(bac_T3_bray_otu) == row.names(bac_T3_bray_zotu)
row.names(bac_T3_bray_otu) ==colnames(bac_T3_bray_otu)
colnames(bac_T3_bray_otu) ==colnames(bac_T3_bray_zotu)
######################################################################################################################
###### Make example NMDS plot and stress plots - Bac and Fungi, OTU and ZOTU
######################################################################################################################
#make the nmds
plotnmds_fungi_T3_otu <- metaMDS(fungi_T3_bray_otu, k=2, trymax=100)
plotnmds_fungi_T3_otu #stress = 0.15
stressplot(plotnmds_fungi_T3_otu)
plotnmds_fungi_T3_zotu <- metaMDS(fungi_T3_bray_zotu, k=2, trymax=100)
plotnmds_fungi_T3_zotu #stress= 0.14
stressplot(plotnmds_fungi_T3_zotu)
plotnmds_bac_T3_otu <- metaMDS(bac_T3_bray_otu, k=3, trymax=100)
plotnmds_bac_T3_otu #stress = 0.067
stressplot(plotnmds_bac_T3_otu)
plotnmds_bac_T3_zotu <- metaMDS(bac_T3_bray_zotu, k=3, trymax=100)
plotnmds_bac_T3_zotu #stress = 0.078
stressplot(plotnmds_bac_T3_zotu)
######################################################################################################################
###### read in predictor files and match thtem
######################################################################################################################
#read in predictor files
map_T3_bac <- read.csv("data/16S_usearch10/16S_T3_map.csv", row.names=1)
##Discard all samples in map that don't appear otu table
map_T3_bac<- map_T3_bac [map_T3_bac$samplename2 %in% rownames(bac_T3_bray_otu), ]
dim(map_T3_bac)
dim(bac_T3_bray_otu)
##Match row names in tables
map_T3_bac<- map_T3_bac[match(rownames(bac_T3_bray_otu),map_T3_bac$samplename2),]
map_T3_bac$samplename2 == rownames(bac_T3_bray_otu)
rownames(bac_T3_bray_otu) ==colnames(bac_T3_bray_otu)
#read in predictor files
map_T3_fungi <- read.csv("data/fungi_ITS2/map_T3_transplant.csv", row.names=1)
##Discard all samples in map that don't appear otu table
map_T3_fungi<- map_T3_fungi [map_T3_fungi$samplename2 %in% rownames(fungi_T3_bray_otu), ]
dim(map_T3_fungi)
dim(fungi_T3_bray_otu)
##Match row names in tables
map_T3_fungi<- map_T3_fungi[match(rownames(fungi_T3_bray_otu),map_T3_fungi$samplename2 ),]
map_T3_fungi$samplename2 == rownames(fungi_T3_bray_otu)
rownames(fungi_T3_bray_otu) ==colnames(fungi_T3_bray_otu)
######################################################################################################################
###### colored by site
######################################################################################################################
names(map_T3_fungi)
map_T3_fungi$colorbysite <- as.character(map_T3_fungi$colorbysite )
map_T3_fungi$colorbyinoc <- as.character(map_T3_fungi$colorbyinoc )
names(map_T3_bac)
map_T3_bac$colorbysite <- as.character(map_T3_bac$colorbysite )
map_T3_bac$colorbyinoc <- as.character(map_T3_bac$colorbyinoc )
plot(scores(plotnmds_fungi_T3_otu), col=map_T3_fungi$colorbyinoc, pch=map_T3_fungi$shapesbysite, cex=1)
plot(scores(plotnmds_fungi_T3_zotu), col=map_T3_fungi$colorbyinoc, pch=map_T3_fungi$shapesbysite, cex=1)
plot(scores(plotnmds_bac_T3_otu), col=map_T3_bac$colorbysite, pch=map_T3_bac$ShapesbyInoculum, cex=1)
plot(scores(plotnmds_bac_T3_zotu), col=map_T3_bac$colorbysite, pch=map_T3_bac$ShapesbyInoculum, cex=1)
######################################################################################################################
###### Figure for Publication A) colored by site B) colored by Inoc with bray curtis sqrt OTU
######################################################################################################################
plot.new()
#change it to make it simpler, no shapes
#make panels C and D
pdf("Figures/otuvzotu_dissimilarity/NMDS_ITS2_T3_OTUvESV.pdf",height=6, width=10)
par(mfrow=c(1,2), mai=c(1,0.8,1,0.2))
#OTU: plot the scores, color by inoculum, pch 16 ffor all, same y and x limits
plot(scores(plotnmds_fungi_T3_otu), col=map_T3_fungi$colorbyinoc, pch=16, cex=1, ylim=c(-0.4,0.6), xlim=c(-0.5,0.6))
#add in stress
mtext("Stress = 0.15", line=-2, adj=0.9, cex=1, side=3)
#add in OTU
mtext("Fungal OTU", line=1, adj=0.5, cex=2, side=3)
#add in panel letter and make it bold
mtext(expression(bold("C")), line=1, adj=0.05, cex=2, side=3)
#ESV: plot the scores, color by inoculum, pch 16 ffor all, same y and x limits
plot(scores(plotnmds_fungi_T3_zotu), col=map_T3_fungi$colorbyinoc, pch=16, cex=1,ylim=c(-0.4,0.6), xlim=c(-0.5,0.6))
#add in stress
mtext("Stress = 0.14", line=-2, adj=0.9, cex=1, side=3)
#add in OTU
mtext("Fungal ESV", line=1, adj=0.5, cex=2, side=3)
#add in panel letter and make it bold
mtext(expression(bold("D")), line=1, adj=0.05, cex=2, side=3)
dev.off()
pdf("Figures/otuvzotu_dissimilarity/NMDS_16S_T3_OTUvESV.pdf",height=6, width=10)
par(mfrow=c(1,2), mai=c(1,0.8,1,0.2))
#OTU: plot the scores, color by site, pch 16 ffor all, same y and x limits
plot(scores(plotnmds_bac_T3_otu), col=map_T3_bac$colorbysite, pch=16, cex=1, ylim=c(-0.4,0.4), xlim=c(-0.55,0.5))
#add in stress
mtext("Stress = 0.067", line=-2, adj=0.9, cex=1, side=3)
#add in OTU
mtext("Bacterial OTU", line=1, adj=0.5, cex=2, side=3)
#add in panel letter and make it bold
mtext(expression(bold("C")), line=1, adj=0.05, cex=2, side=3)
#ESV: plot the scores for the ZOTU, color by site, pch 16 for all, same x and y limits
plot(scores(plotnmds_bac_T3_zotu), col=map_T3_bac$colorbysite, pch=16, cex=1,ylim=c(-0.4,0.4),xlim=c(-0.55,0.5))
#add in stress
mtext("Stress = 0.078", line=-2, adj=0.9, cex=1, side=3)
#add in zotu
mtext("Bacterial ESV", line=1, adj=0.5, cex=2, side=3)
#add in panel letter and make it bold
mtext(expression(bold("D")), line=1, adj=0.05, cex=2, side=3)
dev.off()
par(mfrow=c(1,1),par(mai=c(1.02,0.82,0.82,0.42)))
######################################################################################################################
###### Make same figure with ggplot so I can arrange it with the richness figure
######################################################################################################################
#make data frame for fungi
fungi_T3_otu <- as.data.frame(scores(plotnmds_fungi_T3_otu))
fungi_T3_zotu <- as.data.frame(scores(plotnmds_fungi_T3_zotu))
#bind them togeter and get string for site and inoculum and put in correct order
fungi_T3_otuvzotu <- rbind(fungi_T3_otu,fungi_T3_zotu)
#make a column for Amplicon
fungi_T3_otuvzotu$Amplicon <- c(rep("OTU",nrow(fungi_T3_otu)), rep("ESV",nrow(fungi_T3_zotu)))
#order it to be OTU first
fungi_T3_otuvzotu$Amplicon <- factor(fungi_T3_otuvzotu$Amplicon,levels=c("OTU","ESV"))
#add timepoint
fungi_T3_otuvzotu$Timepoint <- rep("T3",nrow(fungi_T3_otuvzotu))
#add kingdom
fungi_T3_otuvzotu$Kingdom <- rep("Fungi",nrow(fungi_T3_otuvzotu))
#make function to get data frame with site and inoculum in correct order
dataframefunction <- function(df) {
df$Site <- str_sub(row.names(df ), 2, 2)
df$Inoculum <- str_sub(row.names(df), 3, 3)
df$Site <- factor(df$Site,levels=c(1,4,2,3,5))
df$Inoculum <- factor(df$Inoculum,levels=c("D","W","G","P","S"))
return(df)
}
fungi_T3_otuvzotu <- dataframefunction(fungi_T3_otuvzotu)
fungi_T3_otuvzotu
##############################
#make data frame for bacteria
##############################
bac_T3_otu <- as.data.frame(scores(plotnmds_bac_T3_otu))
bac_T3_zotu <- as.data.frame(scores(plotnmds_bac_T3_zotu))
#bind them togeter and get string for site and inoculum and put in correct order
bac_T3_otuvzotu <- rbind(bac_T3_otu,bac_T3_zotu)
#make a column for Amplicon
bac_T3_otuvzotu$Amplicon <- c(rep("OTU",nrow(bac_T3_otu)), rep("ESV",nrow(bac_T3_zotu)))
#order it to be OTU first
bac_T3_otuvzotu$Amplicon <- factor(bac_T3_otuvzotu$Amplicon,levels=c("OTU","ESV"))
#add timepoint
bac_T3_otuvzotu$Timepoint <- rep("T3",nrow(bac_T3_otuvzotu))
#add kingdom
bac_T3_otuvzotu$Kingdom <- rep("Bacteria",nrow(bac_T3_otuvzotu))
bac_T3_otuvzotu <- dataframefunction(bac_T3_otuvzotu)
bac_T3_otuvzotu
#save files
write.csv(bac_T3_otuvzotu, "Figures/otuvzotu_dissimilarity/bacteria_NMDS_scores_otuzotu.csv")
write.csv(fungi_T3_otuvzotu , "Figures/otuvzotu_dissimilarity/fungi_NMDS_scores_otuzotu.csv")
######################################################################################################################
###### Make NMDS figure function in ggplot - colored by site and shapes by inoculum
######################################################################################################################
#make vectors for legend
sitenames <- c("Desert","Scrubland","Grassland","Pine-Oak","Subalpine")
sitecolors <- c("red","orange","green","blue","purple")
shapes <- c(17,3,16,15,18)
#includes an x axis
NMDS_sitecolor_function <- function(df, figurename) {
#make ggplot object
ggplot(df, aes(x=NMDS1, y=NMDS2, col=Inoculum, group=Inoculum)) +
geom_point(size=2) +
theme_bw() +
labs(title = figurename)+
scale_color_manual(values=sitecolors, labels=sitenames) + #add in manual colors and change the legend names
#scale_shape_manual(values=shapes, labels=sitenames) + #add in manual shapes and change the legend name
theme(strip.text.x = element_text(size = 14, colour = "black"), #make T1, T2, T3 labels bigger
axis.text=element_text(size=10,angle=70, hjust=1), #change size of x and y axis tick labels
axis.title=element_text(size=12),#change size of x and y axis labels
legend.spacing = unit(0,"cm"), #change spacing between legends
legend.text=element_text(size=10), #change size of legend text
panel.grid.major = element_blank(), panel.grid.minor = element_blank(), #remove major and minor grid panels
legend.position="none") #remove legend
}
fungalplot <- ggplot(fungi_T3_otuvzotu, aes(x=NMDS1, y=NMDS2, col=Inoculum, group=Inoculum)) + geom_point(size=2) +
theme_bw() +
labs(title = "Fungal NMDS")+
scale_color_manual(values=sitecolors, labels=sitenames) + #add in manual colors and change the legend names
#scale_shape_manual(values=shapes, labels=sitenames) + #add in manual shapes and change the legend name
theme(strip.text.x = element_text(size = 14, colour = "black"), #make T1, T2, T3 labels bigger
axis.text=element_text(size=12), #change size of x and y axis tick labels
axis.title=element_text(size=14),#change size of x and y axis labels
legend.spacing = unit(0,"cm"), #change spacing between legends
legend.text=element_text(size=10), #change size of legend text
panel.grid.major = element_blank(), panel.grid.minor = element_blank(), #remove major and minor grid panels
legend.position="none") #remove legend
fungalplot + facet_wrap(~Amplicon) +stat_ellipse() #include ellipse around the group
#make panel names function
panelnamesfungi <- c('OTU'="C. Fungal OTU",'ESV'="D. Fungal ESV")
fungalplot + facet_wrap(~Amplicon,labeller=as_labeller(panelnamesfungi )) +stat_ellipse() #include ellipse around the group
#make bacteria plot
bacterialplot <- ggplot(bac_T3_otuvzotu, aes(x=NMDS1, y=NMDS2, col=Site, group=Site)) + geom_point(size=2) +
theme_bw() +
#labs(title = c("Bacterial NMDS")+
scale_color_manual(values=sitecolors, labels=sitenames) + #add in manual colors and change the legend names
#scale_shape_manual(values=shapes, labels=sitenames) + #add in manual shapes and change the legend name
theme(strip.text.x = element_text(size = 14, colour = "black"), #make T1, T2, T3 labels bigger
axis.text=element_text(size=12), #change size of x and y axis tick labels
axis.title=element_text(size=14),#change size of x and y axis labels
legend.spacing = unit(0,"cm"), #change spacing between legends
legend.text=element_text(size=10), #change size of legend text
panel.grid.major = element_blank(), panel.grid.minor = element_blank(), #remove major and minor grid panels
legend.position="none") #remove legend
bacterialplot + facet_wrap(~Amplicon) +stat_ellipse() #include ellipse around the group
#make panel names function
panelnamesbac <- c('OTU'="C. Bacterial OTU",'ESV'="D. Bacterial ESV")
#make NMDS for colors site and shapes inoculum
#bacteria all time points faceted
bacterialplot + facet_wrap(~Amplicon,labeller=as_labeller(panelnamesbac)) +stat_ellipse() #include ellipse around the group