-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathtest_distance_methods_on_MDS.R
executable file
·87 lines (85 loc) · 3.37 KB
/
test_distance_methods_on_MDS.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
# Test the impact of different distance methods
# NOT: requires a phyloseq object!
# In plot_ordination, we coloured based on an env variable. replace as you see fit.
dists <- unlist(distanceMethodList)
dists <- dists[-(1:3)]
dists = dists[-which(dists=="ANY")]
dists = dists[-which(dists=="mountford")]
plist <- vector("list", length=length(dists))
names(plist) = dists
for( i in dists ){
# Calculate distance matrix
iDist <- distance(physeq_agg_rel_n0, method=i)
# Calculate ordination
iMDS <- ordinate(physeq_agg_rel_n0, "MDS", distance=iDist)
## Make plot
# Don't carry over previous plot (if error, p will be blank)
p <- NULL
# Create plot, store as temp variable, p
p <- plot_ordination(physeq_agg_rel_n0, iMDS, color="Depth.category")
# Add title to each plot
p <- p + ggtitle(paste("MDS using distance method ", i, sep=""))
# Save the graphic to file.
plist[[i]] = p
}
df = ldply(plist, function(x) x$data)
names(df)[1] <- "distance"
p = ggplot(df, aes(Axis.1, Axis.2, color=Depth.category)) +
geom_point(size=0.5, alpha=0.5) +
facet_wrap(~distance, scales="free") +
theme_bw() +
theme(
axis.text = element_blank(),
axis.title = element_blank(),
axis.ticks = element_blank()
) +
scale_color_manual(values = c("red2", "forestgreen", "blue3", "yellow2", "grey"))
ggsave("/Users/christopherhempel/Desktop/distance-summary-plot-rsde.png", p, width = 13, height = 9, units = "in")
# Check distance effect on clr transformed data
# Note:
# If any values are zero, the clr transform routine first adds a small pseudocount of min(relative abundance)/2 to all values. To avoid this, you can replace any zeros in advance by setting zero_replace to a number > 0.
# Better do manually
physeq_agg_rel_n0_clr <- microbiome::transform(physeq_agg_rel_n0, "clr")
# Test the impact of different distance methods
dists <- unlist(distanceMethodList)
dists <- dists[-(1:3)]
dists = dists[-which(dists=="ANY")]
dists = dists[-which(dists=="mountford")]
dists = dists[-which(dists=="bray")]
dists = dists[-which(dists=="kulczynski")]
dists = dists[-which(dists=="jaccard")]
dists = dists[-which(dists=="morisita")]
dists = dists[-which(dists=="horn")]
dists = dists[-which(dists=="binomial")]
dists = dists[-which(dists=="chao")]
dists = dists[-which(dists=="binary")]
plist <- vector("list", length=length(dists))
names(plist) = dists
for( i in dists ){
# Calculate distance matrix
iDist <- distance(physeq_agg_rel_n0_clr, method=i)
# Calculate ordination
iMDS <- ordinate(physeq_agg_rel_n0_clr, "MDS", distance=iDist)
## Make plot
# Don't carry over previous plot (if error, p will be blank)
p <- NULL
# Create plot, store as temp variable, p
p <- plot_ordination(physeq_agg_rel_n0_clr, iMDS, color="Depth.category")
# Add title to each plot
p <- p + ggtitle(paste("MDS using distance method ", i, sep=""))
# Save the graphic to file.
plist[[i]] = p
}
df = ldply(plist, function(x) x$data)
names(df)[1] <- "distance"
p = ggplot(df, aes(Axis.1, Axis.2, color=Depth.category)) +
geom_point(size=0.5, alpha=0.5) +
facet_wrap(~distance, scales="free") +
theme_bw() +
theme(
axis.text = element_blank(),
axis.title = element_blank(),
axis.ticks = element_blank()
) +
scale_color_manual(values = c("red2", "forestgreen", "blue3", "yellow2", "grey"))
ggsave("/Users/christopherhempel/Desktop/distance-summary-plot-rsde_clr.png", p, width = 13, height = 9, units = "in")