Load in data from pb_analysis.Rmd
load("../data/pb_analysis_to_figures.RData")
load("../data/taxon_analysis.RData")
source("../data/pb_functions_for_manuscript.R")
Load packages
library(vegan)
library(labdsv)
library(ape)
library(xtable)
Unpack list of groups.
for (i in 1:length(names(groups))) {
assign(names(groups)[i], unlist(groups[[i]]))
}
Info for manuscript:
What percentage were Methylobacterium sp. sequences?
sum(pb.1000[, "234992"])/sum(pb.1000) * 100
## [1] 13.34
How many samples were dominated by this OTU?
meth <- 0
for (i in 1:nrow(pb.1000)) {
if (names(rev(sort(pb.1000[i, ])))[1] == "234992") {
meth <- meth + 1
}
}
meth/nrow(pb.1000)
## [1] 0.327
rm(meth, i)
Most dominant in occupied samples
names(rev(sort(colSums(pb.1000[pb.map$location == "occ", ]))))[1:3]
## [1] "25506" "234992" "94008"
sum(pb.1000[pb.map$location == "occ", "25506"])/sum(pb.1000[pb.map$location ==
"occ", ])
## [1] 0.1482
sum(pb.1000[pb.map$location == "unocc", "25506"])/sum(pb.1000[pb.map$location ==
"unocc", ])
## [1] 0.03351
Create color shortcuts
s1.col <- "#91BFDB"
s2.col <- "#FC8D59"
s3.col <- "#FFFFBF"
Put plotting options in mapping table.
# assign colors
pb.map$bg[pb.map$person == "s1"] <- "#91BFDB" # blue
pb.map$bg[pb.map$person == "s2"] <- "#FC8D59" # orange
pb.map$bg[pb.map$person == "s3"] <- "#FFFFBF" # yellow
pb.map$bg2 <- pb.map$bg
pb.map$bg2[pb.map$location == "unocc"] <- "gray50"
pb.map$bor <- "gray30"
All trials together, regardless of person or collection method.
all.occ.v.unocc <- adonis(pb.can ~ pb.map$location)
dim(pb.map)
[1] 211 9
print(xtable(all.occ.v.unocc$aov.tab), type = "html")
Df | SumsOfSqs | MeanSqs | F.Model | R2 | Pr(> F) | |
---|---|---|---|---|---|---|
pb.map$location | 1 | 1.04 | 1.04 | 2.20 | 0.01 | 0.0010 |
Residuals | 209 | 98.53 | 0.47 | 0.99 | ||
Total | 210 | 99.57 | 1.00 |
All 4hr together regardless of person or method.
all.4.occ.v.unocc <- adonis(as.dist(as.matrix(pb.can)[pb.map$duration == 4,
pb.map$duration == 4]) ~ pb.map$location[pb.map$duration == 4])
print(xtable(all.4.occ.v.unocc$aov.tab), type = "html")
Df | SumsOfSqs | MeanSqs | F.Model | R2 | Pr(> F) | |
---|---|---|---|---|---|---|
pb.map$location[pb.map$duration == 4] | 1 | 0.88 | 0.88 | 1.88 | 0.02 | 0.0010 |
Residuals | 104 | 48.71 | 0.47 | 0.98 | ||
Total | 105 | 49.59 | 1.00 |
All 2hr together regardless of person or method.
all.2.occ.v.unocc <- adonis(as.dist(as.matrix(pb.can)[pb.map$duration == 2,
pb.map$duration == 2]) ~ pb.map$location[pb.map$duration == 2])
print(xtable(all.2.occ.v.unocc$aov.tab), type = "html")
Df | SumsOfSqs | MeanSqs | F.Model | R2 | Pr(> F) | |
---|---|---|---|---|---|---|
pb.map$location[pb.map$duration == 2] | 1 | 0.67 | 0.67 | 1.41 | 0.01 | 0.0010 |
Residuals | 103 | 48.79 | 0.47 | 0.99 | ||
Total | 104 | 49.46 | 1.00 |
All petris together regardless of person.
all.p.occ.v.unocc <- adonis(as.dist(as.matrix(pb.can)[pb.map$sampleType == "petri.dish",
pb.map$sampleType == "petri.dish"]) ~ pb.map$location[pb.map$sampleType ==
"petri.dish"])
print(xtable(all.p.occ.v.unocc$aov.tab), type = "html")
Df | SumsOfSqs | MeanSqs | F.Model | R2 | Pr(> F) | |
---|---|---|---|---|---|---|
pb.map$location[pb.map$sampleType == “petri.dish”] | 1 | 0.63 | 0.63 | 1.32 | 0.02 | 0.0010 |
Residuals | 69 | 32.75 | 0.47 | 0.98 | ||
Total | 70 | 33.38 | 1.00 |
All filters together regardless of person.
all.f.occ.v.unocc <- adonis(as.dist(as.matrix(pb.can)[pb.map$sampleType != "petri.dish",
pb.map$sampleType != "petri.dish"]) ~ pb.map$location[pb.map$sampleType !=
"petri.dish"])
print(xtable(all.f.occ.v.unocc$aov.tab), type = "html")
Df | SumsOfSqs | MeanSqs | F.Model | R2 | Pr(> F) | |
---|---|---|---|---|---|---|
pb.map$location[pb.map$sampleType != “petri.dish”] | 1 | 0.90 | 0.90 | 1.92 | 0.01 | 0.0010 |
Residuals | 138 | 64.74 | 0.47 | 0.99 | ||
Total | 139 | 65.64 | 1.00 |
These are all for 4hr filters.
occ.v.unocc <- adonis(occ4f.can ~ occ4f.map$location)
person.v.unocc <- adonis(occ4f.can ~ occ4f.map$occ.person2)
person.v.eachother <- adonis(as.dist(as.matrix(occ4f.can)[occ4f.map$location ==
"occ", occ4f.map$location == "occ"]) ~ subset(occ4f.map, occ4f.map$location ==
"occ")$person)
only.unocc <- adonis(as.dist(as.matrix(occ4f.can)[occ4f.map$location == "unocc",
occ4f.map$location == "unocc"]) ~ subset(occ4f.map, occ4f.map$location ==
"unocc")$person)
These are all for 4hr settling dishes.
occ.v.unocc.p <- adonis(occ4p.can ~ occ4p.map$location)
person.v.unocc.p <- adonis(occ4p.can ~ occ4p.map$occ.person2)
person.v.eachother.p <- adonis(as.dist(as.matrix(occ4p.can)[occ4p.map$location ==
"occ", occ4p.map$location == "occ"]) ~ subset(occ4p.map, occ4p.map$location ==
"occ")$person)
only.unocc.p <- adonis(as.dist(as.matrix(occ4p.can)[occ4p.map$location == "unocc",
occ4p.map$location == "unocc"]) ~ subset(occ4p.map, occ4p.map$location ==
"unocc")$person)
print(xtable(occ.v.unocc$aov.tab), type = "html")
Df | SumsOfSqs | MeanSqs | F.Model | R2 | Pr(> F) | |
---|---|---|---|---|---|---|
occ4f.map$location | 1 | 0.79 | 0.79 | 1.69 | 0.02 | 0.0010 |
Residuals | 69 | 32.13 | 0.47 | 0.98 | ||
Total | 70 | 32.92 | 1.00 |
Samples are significantly different when the room is occupied.
print(xtable(person.v.unocc$aov.tab), type = "html")
Df | SumsOfSqs | MeanSqs | F.Model | R2 | Pr(> F) | |
---|---|---|---|---|---|---|
occ4f.map$occ.person2 | 3 | 2.06 | 0.69 | 1.49 | 0.06 | 0.0010 |
Residuals | 67 | 30.86 | 0.46 | 0.94 | ||
Total | 70 | 32.92 | 1.00 |
And there is a difference when you consider who was occupying the room.
print(xtable(person.v.eachother$aov.tab), type = "html")
Df | SumsOfSqs | MeanSqs | F.Model | R2 | Pr(> F) | |
---|---|---|---|---|---|---|
subset(occ4f.map, occ4f.map$location == “occ”)$person | 2 | 1.28 | 0.64 | 1.43 | 0.08 | 0.0010 |
Residuals | 33 | 14.78 | 0.45 | 0.92 | ||
Total | 35 | 16.06 | 1.00 |
And also a big difference when you compare the 3 people against one another.
How about if we only look at the unoccupied side of the room? There should be a temporal diff, but we would expect the variation explained to be less than that when only looking at person.v.eachother.
print(xtable(only.unocc$aov.tab), type = "html")
Df | SumsOfSqs | MeanSqs | F.Model | R2 | Pr(> F) | |
---|---|---|---|---|---|---|
subset(occ4f.map, occ4f.map$location == “unocc”)$person | 2 | 0.96 | 0.48 | 1.02 | 0.06 | 0.1420 |
Residuals | 32 | 15.11 | 0.47 | 0.94 | ||
Total | 34 | 16.07 | 1.00 |
While there is a difference in the right direction, this is a bit dissappointing. There is a pretty darn significant difference between three unoccupied trails, and only somewhat better than the occupied difference.
A dataframe with all pertinent adonis results. REARRANGE FOR FINAL TABLE
adonis.df <- data.frame(rbind(all.occ.v.unocc$aov.tab[1, 4:6], all.4.occ.v.unocc$aov.tab[1,
4:6], all.2.occ.v.unocc$aov.tab[1, 4:6], all.f.occ.v.unocc$aov.tab[1, 4:6],
all.p.occ.v.unocc$aov.tab[1, 4:6], occ.v.unocc$aov.tab[1, 4:6], person.v.unocc$aov.tab[1,
4:6], person.v.eachother$aov.tab[1, 4:6], only.unocc$aov.tab[1, 4:6],
occ.v.unocc.p$aov.tab[1, 4:6], person.v.unocc.p$aov.tab[1, 4:6], person.v.eachother.p$aov.tab[1,
4:6], only.unocc.p$aov.tab[1, 4:6]))
row.names(adonis.df) <- c("AllSamplesOccVsUnocc", "All4hrOccVsUnocc", "All2hrOccVsUnocc",
"AllFiltersOccVsUnocc", "AllPetrisOccVsUnocc", "OccVsUnocc4hrFilters", "UnoccVs3people4hrFilters",
"OccupantsVsEachOtherFilters", "OnlyUnoccTestFilters", "OccVsUnocc4hrDishes",
"UnoccVs3people4hrDishes", "OccupantsVsEachOtherDishes", "OnlyUnoccTestDishes")
adonis.df$n <- c(nrow(pb.map), nrow(pb.map[pb.map$duration == "4", ]), nrow(pb.map[pb.map$duration ==
"2", ]), nrow(pb.map[pb.map$sampleType == "filter", ]), nrow(pb.map[pb.map$sampleType ==
"petri.dish", ]), nrow(occ4f.map), nrow(occ4f.map), nrow(occ4f.map[occ4f.map$location ==
"occ", ]), nrow(occ4f.map[occ4f.map$location == "unocc", ]), nrow(occ4p.map),
nrow(occ4p.map), nrow(occ4p.map[occ4p.map$location == "occ", ]), nrow(occ4p.map[occ4p.map$location ==
"unocc", ]))
adonis.df <- adonis.df[, c(4, 1, 2, 3)]
names(adonis.df) <- c("n", "F-statistic", "R2", "P-value")
print(xtable(adonis.df, digits = c(0, 0, 2, 3, 3)), type = "html")
n | F-statistic | R2 | P-value | |
---|---|---|---|---|
AllSamplesOccVsUnocc | 211 | 2.20 | 0.010 | 0.001 |
All4hrOccVsUnocc | 106 | 1.88 | 0.018 | 0.001 |
All2hrOccVsUnocc | 105 | 1.41 | 0.013 | 0.001 |
AllFiltersOccVsUnocc | 140 | 1.92 | 0.014 | 0.001 |
AllPetrisOccVsUnocc | 71 | 1.32 | 0.019 | 0.001 |
OccVsUnocc4hrFilters | 71 | 1.69 | 0.024 | 0.001 |
UnoccVs3people4hrFilters | 71 | 1.49 | 0.063 | 0.001 |
OccupantsVsEachOtherFilters | 36 | 1.43 | 0.080 | 0.001 |
OnlyUnoccTestFilters | 35 | 1.02 | 0.060 | 0.142 |
OccVsUnocc4hrDishes | 35 | 1.24 | 0.036 | 0.001 |
UnoccVs3people4hrDishes | 35 | 1.15 | 0.100 | 0.001 |
OccupantsVsEachOtherDishes | 18 | 1.13 | 0.131 | 0.001 |
OnlyUnoccTestDishes | 17 | 1.00 | 0.125 | 0.436 |
# write.table(adonis.df,
# '~/Dropbox/pb_shared_markdown/tables/adonisTable.txt', sep='\t',
# quote=FALSE)
s14f <- c(s14of, s14uf)
s14p <- c(s14op, s14up)
s24f <- c(s24of, s24uf)
s24p <- c(s24op, s24up)
s34f <- c(s34of, s34uf)
s34p <- c(s34op, s34up)
s12f <- c(s12of, s12uf)
s12p <- c(s12op, s12up)
s22f <- c(s22of, s22uf)
s22p <- c(s22op, s22up)
s32f <- c(s32of, s32uf)
s32p <- c(s32op, s32up)
adonis.groups <- c("s14f", "s24f", "s34f", "s12f", "s22f", "s32f", "s14p", "s24p",
"s34p", "s12p", "s22p", "s32p")
trial.adonis.table <- data.frame(n = 0, `F-statistic` = 0, R2 = 0, `P-value` = 0)
for (i in 1:length(adonis.groups)) {
these <- get(adonis.groups[i])
model <- adonis(vegdist(pb.1000[these, ], "canberra") ~ pb.map$location[these])
model.name <- paste(adonis.groups[i], ".occ.v.unocc", sep = "")
assign(model.name, model)
tab <- get(model.name)$aov.tab
trial.adonis.table[i, 1] <- length(these)
trial.adonis.table[i, 2:4] <- tab[1, 4:6]
row.names(trial.adonis.table)[i] <- adonis.groups[i]
}
trial.adonis.table$person <- c(rep(c("s1", "s2", "s3"), 4))
trial.adonis.table$duration <- c(rep(rep(c(4, 2), each = 3), 2))
trial.adonis.table$sampleType <- c(rep(c("air filter", "settling dish"), each = 6))
trial.adonis.table <- trial.adonis.table[, c(5, 6, 7, 1, 2, 3, 4)]
print(xtable(trial.adonis.table, digits = c(0, 0, 0, 0, 0, 2, 3, 3)), type = "html")
person | duration | sampleType | n | F.statistic | R2 | P.value | |
---|---|---|---|---|---|---|---|
s14f | s1 | 4 | air filter | 24 | 1.48 | 0.063 | 0.001 |
s24f | s2 | 4 | air filter | 23 | 1.10 | 0.050 | 0.003 |
s34f | s3 | 4 | air filter | 24 | 1.53 | 0.065 | 0.001 |
s12f | s1 | 2 | air filter | 24 | 1.40 | 0.060 | 0.001 |
s22f | s2 | 2 | air filter | 23 | 1.06 | 0.048 | 0.029 |
s32f | s3 | 2 | air filter | 22 | 1.05 | 0.050 | 0.019 |
s14p | s1 | 4 | settling dish | 12 | 1.14 | 0.102 | 0.004 |
s24p | s2 | 4 | settling dish | 11 | 1.06 | 0.105 | 0.045 |
s34p | s3 | 4 | settling dish | 12 | 1.15 | 0.103 | 0.003 |
s12p | s1 | 2 | settling dish | 12 | 1.08 | 0.097 | 0.006 |
s22p | s2 | 2 | settling dish | 12 | 1.04 | 0.094 | 0.104 |
s32p | s3 | 2 | settling dish | 12 | 1.02 | 0.093 | 0.152 |
# write.table(trial.adonis.table,
# '~/Dropbox/pb_shared_markdown/tables/adonisTrialTable.txt', sep='\t',
# quote=FALSE)
Was Staph. epideridis more common in occupied samples? - Not all of them.
for (i in 1:length(adonis.groups)) {
vect <- pb.1000[get(adonis.groups[i]), "25506"]
loc <- pb.map$location[get(adonis.groups[i])]
print(t.test(vect ~ loc))
}
t.test(pb.1000[, "25506"] ~ pb.map$location)
##
## Welch Two Sample t-test
##
## data: pb.1000[, "25506"] by pb.map$location
## t = 8.466, df = 138.8, p-value = 3.224e-14
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 87.94 141.53
## sample estimates:
## mean in group occ mean in group unocc
## 148.25 33.51
148.25/33.51
## [1] 4.424
titles <- c("4hr filters", "2hr filters", "4hr dishes", "2hr dishes")
for (i in 1:4) {
print(xtable(betad.aov.list[[i]]), type = "html", comment = FALSE)
boxplot(betad.table.list[[i]]$dists ~ betad.table.list[[i]]$person, col = "gray",
notch = TRUE, main = titles[i])
}
Df | Sum Sq | Mean Sq | F value | Pr(> F) | |
---|---|---|---|---|---|
Groups | 2 | 0.00 | 0.00 | 14.38 | 0.0000 |
Residuals | 33 | 0.00 | 0.00 |
Df | Sum Sq | Mean Sq | F value | Pr(> F) | |
---|---|---|---|---|---|
Groups | 2 | 0.00 | 0.00 | 25.24 | 0.0000 |
Residuals | 33 | 0.00 | 0.00 |
Df | Sum Sq | Mean Sq | F value | Pr(> F) | |
---|---|---|---|---|---|
Groups | 2 | 0.00 | 0.00 | 8.46 | 0.0035 |
Residuals | 15 | 0.00 | 0.00 |
Df | Sum Sq | Mean Sq | F value | Pr(> F) | |
---|---|---|---|---|---|
Groups | 2 | 0.00 | 0.00 | 2.54 | 0.1122 |
Residuals | 15 | 0.00 | 0.00 |
Differences between people?
titles <- c("4hr filters", "2hr filters", "4hr dishes", "2hr dishes")
set <- c("occ4f.map", "occ2f.map", "occ4p.map", "occ2p.map")
for (i in 1:4) {
map <- get(set[i])
map <- map[map$location == "occ", ]
hrj <- pb.hrj[row.names(map), ]
print(xtable(summary(lm(hrj$H1 ~ map$person))), type = "html", comment = FALSE)
boxplot(hrj$H1 ~ map$person, col = "gray", notch = TRUE, main = titles[i])
}
Estimate | Std. Error | t value | Pr(> |t|) | |
---|---|---|---|---|
(Intercept) | 3.6197 | 0.1639 | 22.08 | 0.0000 |
map$persons2 | -0.3490 | 0.2319 | -1.51 | 0.1417 |
map$persons3 | -0.1614 | 0.2319 | -0.70 | 0.4912 |
Estimate | Std. Error | t value | Pr(> |t|) | |
---|---|---|---|---|
(Intercept) | 3.6358 | 0.1324 | 27.46 | 0.0000 |
map$persons2 | -0.3675 | 0.1872 | -1.96 | 0.0582 |
map$persons3 | 0.1825 | 0.1872 | 0.97 | 0.3369 |
Estimate | Std. Error | t value | Pr(> |t|) | |
---|---|---|---|---|
(Intercept) | 3.0293 | 0.1848 | 16.39 | 0.0000 |
map$persons2 | -0.5483 | 0.2614 | -2.10 | 0.0533 |
map$persons3 | -0.2512 | 0.2614 | -0.96 | 0.3517 |
Estimate | Std. Error | t value | Pr(> |t|) | |
---|---|---|---|---|
(Intercept) | 3.4810 | 0.3374 | 10.32 | 0.0000 |
map$persons2 | -0.7081 | 0.4772 | -1.48 | 0.1585 |
map$persons3 | -0.4549 | 0.4772 | -0.95 | 0.3556 |
Different people 4 hr filters.
titles <- c("S1", "S2", "S3")
s123 <- c("s1", "s2", "s3")
for (i in 1:3) {
map <- occ4f.map[occ4f.map$person == s123[i], ]
hrj <- pb.hrj[row.names(map), ]
print(xtable(summary(lm(hrj$H1 ~ map$location))), type = "html", comment = FALSE)
boxplot(hrj$H1 ~ map$location, col = "gray", notch = TRUE, main = titles[i])
}
Estimate | Std. Error | t value | Pr(> |t|) | |
---|---|---|---|---|
(Intercept) | 3.6197 | 0.1707 | 21.20 | 0.0000 |
map$locationunocc | -0.1820 | 0.2414 | -0.75 | 0.4590 |
Estimate | Std. Error | t value | Pr(> |t|) | |
---|---|---|---|---|
(Intercept) | 3.2707 | 0.2379 | 13.75 | 0.0000 |
map$locationunocc | 0.5056 | 0.3440 | 1.47 | 0.1565 |
Estimate | Std. Error | t value | Pr(> |t|) | |
---|---|---|---|---|
(Intercept) | 3.4584 | 0.1703 | 20.31 | 0.0000 |
map$locationunocc | 0.4448 | 0.2408 | 1.85 | 0.0783 |
# pdf('../../figures/ordination_cluster.pdf', width=10, height=6.5)
layout(matrix(c(1:3, 7, 7, 7, 4:6, 8, 8, 9), 3, 4), widths = c(1, 1.4, 1, 1.4),
heights = c(1.1, 1, 1))
# 4 hour filters
ord <- nmds(vegdist(pb.1000[c(s14of, s14uf), ], "canberra"))
## initial value 28.336845
## iter 5 value 20.023322
## final value 19.658494
## converged
map <- pb.map[c(s14of, s14uf), ]
cols <- s1.col
par(mar = c(1, 1, 3, 1))
plot(ord$points, xaxt = "n", yaxt = "n", xlab = "", ylab = "", type = "n", bty = "n")
rect(par()$usr[1], par()$usr[3], par()$usr[2], par()$usr[4], col = "gray90",
border = "gray60")
points(ord$points, cex = ifelse(map$location == "occ", 2, 1.5), pch = ifelse(map$location ==
"occ", 21, 24), bg = map$bg2, col = "gray30")
mtext("Air filters", line = 0.2)
mtext("(a)", side = 3, adj = 0, font = 2, col = "gray30", line = 1, cex = 1.2)
mtext(c("NMDS 1", "NMDS 2"), side = c(1, 2), adj = 1, cex = 0.7, line = c(0.15,
0))
ord <- nmds(vegdist(pb.1000[c(s24of, s24uf), ], "canberra"))
## initial value 34.606164
## iter 5 value 26.582735
## iter 10 value 25.642792
## iter 15 value 25.355035
## final value 25.161056
## converged
map <- pb.map[c(s24of, s24uf), ]
cols <- s2.col
par(mar = c(1, 1, 1, 1))
plot(ord$points[, 1], ord$points[, 2], xaxt = "n", yaxt = "n", xlab = "", ylab = "",
type = "n", bty = "n")
rect(par()$usr[1], par()$usr[3], par()$usr[2], par()$usr[4], col = "gray90",
border = "gray60")
points(ord$points, cex = ifelse(map$location == "occ", 2, 1.5), pch = ifelse(map$location ==
"occ", 21, 24), bg = map$bg2, col = "gray30")
ord <- nmds(vegdist(pb.1000[c(s34of, s34uf), ], "canberra"))
## initial value 30.187730
## iter 5 value 18.966172
## iter 10 value 18.070539
## iter 10 value 18.060363
## final value 17.918395
## converged
map <- pb.map[c(s34of, s34uf), ]
cols <- s3.col
par(mar = c(1, 1, 1, 1))
plot(ord$points[, 1], ord$points[, 2], xaxt = "n", yaxt = "n", xlab = "", ylab = "",
type = "n", bty = "n")
rect(par()$usr[1], par()$usr[3], par()$usr[2], par()$usr[4], col = "gray90",
border = "gray60")
points(ord$points[, 1], ord$points[, 2], cex = ifelse(map$location == "occ",
2, 1.5), pch = ifelse(map$location == "occ", 21, 24), bg = map$bg2, col = "gray30")
# 4hr petris
ord <- nmds(vegdist(pb.1000[c(s14op, s14up), ], "canberra"))
## initial value 29.092834
## iter 5 value 16.029246
## iter 10 value 15.695975
## final value 15.653145
## converged
map <- pb.map[c(s14op, s14up), ]
cols <- s1.col
par(mar = c(1, 1, 3, 1))
plot(ord$points, xaxt = "n", yaxt = "n", xlab = "", ylab = "", type = "n", bty = "n")
rect(par()$usr[1], par()$usr[3], par()$usr[2], par()$usr[4], col = "gray90",
border = "gray60")
points(ord$points, cex = ifelse(map$location == "occ", 2, 1.5), pch = ifelse(map$location ==
"occ", 21, 24), bg = map$bg2, col = "gray30")
mtext("Settling dishes", line = 0.2)
mtext("(c)", side = 3, adj = 0, font = 2, col = "gray30", line = 1, cex = 1.2)
mtext(c("NMDS 1", "NMDS 2"), side = c(1, 2), adj = 1, cex = 0.7, line = c(0.15,
0))
ord <- nmds(vegdist(pb.1000[c(s24op, s24up), ], "canberra"))
## initial value 29.344912
## iter 5 value 15.410212
## iter 10 value 14.416981
## final value 14.269691
## converged
map <- pb.map[c(s24op, s24up), ]
cols <- s2.col
par(mar = c(1, 1, 1, 1))
plot(ord$points, xaxt = "n", yaxt = "n", xlab = "", ylab = "", type = "n", bty = "n")
rect(par()$usr[1], par()$usr[3], par()$usr[2], par()$usr[4], col = "gray90",
border = "gray60")
points(ord$points, cex = ifelse(map$location == "occ", 2, 1.5), pch = ifelse(map$location ==
"occ", 21, 24), bg = map$bg2, col = "gray30")
ord <- nmds(vegdist(pb.1000[c(s34op, s34up), ], "canberra"))
## initial value 22.684332
## iter 5 value 14.698852
## iter 10 value 14.017256
## iter 15 value 13.302352
## final value 12.974551
## converged
map <- pb.map[c(s34op, s34up), ]
cols <- s3.col
par(mar = c(1, 1, 1, 1))
plot(ord$points, xaxt = "n", yaxt = "n", xlab = "", ylab = "", type = "n", bty = "n")
rect(par()$usr[1], par()$usr[3], par()$usr[2], par()$usr[4], col = "gray90",
border = "gray60")
points(ord$points, cex = ifelse(map$location == "occ", 2, 1.5), pch = ifelse(map$location ==
"occ", 21, 24), bg = map$bg2, col = "gray30")
# 4hr filters cluster
comm <- occ4f.table
dis <- occ4f.can
map <- pb.map[row.names(occ4f.table), ]
make.file <- FALSE
off = -2.5
edge = c(0, 3)
line = TRUE
ave <- hclust(dis, "average")
phy <- ladderize(as.phylo(ave))
map$pch <- ifelse(map$location == "occ", 21, 24)
map$bg2 <- map$bg
map$bg2[map$pch == 24] <- "gray50"
map$bor <- "gray30"
par(mar = c(0, 0, 1, 0))
plot(phy, use.edge.length = FALSE, edge.color = "gray30", cex = 0.5, no.margin = FALSE,
show.tip.label = FALSE, direction = "leftwards")
if (line) {
segments(rep(edge[1], nrow(comm)), 1:nrow(comm), rep(edge[2], nrow(comm)),
1:nrow(comm), col = "gray30", lwd = 1)
}
tiplabels(pch = map$pch, bg = map$bg2, adj = off, cex = ifelse(map$pch == 21,
1.8, 1.2), col = "gray30")
legend(38, 15, legend = c("occupied", "unoccupied", "Subject 1", "Subject 2",
"Subject 3"), bty = "n", pch = c(21, 24, 22, 22, 22), col = c(rep("gray40",
5)), text.col = "gray30", pt.bg = c("gray50", "gray50", s1.col, s2.col,
s3.col), cex = 1.5, pt.cex = 2.5, y.intersp = c(1, 1, 1.4, 1.27, 1.2))
segments(44, 8, 63, 8, col = "gray70")
mtext("(b)", side = 3, at = 3, font = 2, col = "gray30", line = -1, cex = 1.2)
######## 4hr dishes cluster
comm <- occ4p.table
dis <- occ4p.can
map <- pb.map[row.names(occ4p.table), ]
make.file <- FALSE
off = 0
edge = c(0, 3)
line = FALSE
ave <- hclust(dis, "average")
phy <- ladderize(as.phylo(ave))
map$pch <- ifelse(map$location == "occ", 21, 24)
map$bg2 <- map$bg
map$bg2[map$pch == 24] <- "gray50"
map$bor <- "gray30"
par(mar = c(5, 0, 2, 0))
plot(phy, use.edge.length = FALSE, edge.color = "gray30", cex = 0.5, no.margin = FALSE,
show.tip.label = FALSE, direction = "leftwards")
if (line) {
segments(rep(edge[1], nrow(comm)), 1:nrow(comm), rep(edge[2], nrow(comm)),
1:nrow(comm), col = "gray30", lwd = 1)
}
tiplabels(pch = map$pch, bg = map$bg2, adj = off, cex = ifelse(map$pch == 21,
1.8, 1.2), col = "gray30")
par(xpd = TRUE)
mtext("(d)", side = 3, at = 3, font = 2, col = "gray30", line = 0.1, cex = 1.2)
mtext("(e)", side = 1, at = 3, font = 2, col = "gray30", cex = 1.2, line = 3)
# set up for barplot in the corner
pb.mat <- as.matrix(pb.can)
pb.matl <- pb.mat
pb.matl[upper.tri(pb.matl)] <- 0
s1 <- c(s14of, s14op)
s2 <- c(s24of, s24op)
s3 <- c(s34of, s34op)
s1f <- s14of
s2f <- s24of
s3f <- s34of
s1p <- s14op
s2p <- s24op
s3p <- s34op
self <- c(pb.matl[s1, s1], pb.matl[s2, s2], pb.matl[s3, s3])
other <- c(pb.mat[s1, c(s2, s3)], pb.mat[s2, c(s1, s3)], pb.mat[s3, c(s1, s2)])
selff <- c(pb.matl[s1f, s1f], pb.matl[s2f, s2f], pb.matl[s3f, s3f])
otherf <- c(pb.mat[s1f, c(s2f, s3f)], pb.mat[s2f, c(s1f, s3f)], pb.mat[s3f,
c(s1f, s2f)])
selfp <- c(pb.matl[s1p, s1p], pb.matl[s2p, s2p], pb.matl[s3p, s3p])
otherp <- c(pb.mat[s1p, c(s2p, s3p)], pb.mat[s2p, c(s1p, s3p)], pb.mat[s3p,
c(s1p, s2p)])
self <- self[-c(which(self == 0))]
selff <- selff[-c(which(selff == 0))]
selfp <- selfp[-c(which(selfp == 0))]
self <- 1 - self
selff <- 1 - selff
selfp <- 1 - selfp
other <- 1 - other
otherf <- 1 - otherf
otherp <- 1 - otherp
length(self)
## [1] 459
length(other)
## [1] 1944
length(selff)
## [1] 198
length(otherf)
## [1] 864
length(selfp)
## [1] 45
length(otherp)
## [1] 216
means <- data.frame(self = c(mean(self), mean(selff), mean(selfp)), other = c(mean(other),
mean(otherf), mean(otherp)), row.names = c("all", "filter", "petri"))
ses <- data.frame(self = c(se(self), se(selff), se(selfp)), other = c(se(other),
se(otherf), se(otherp)), row.names = c("all", "filter", "petri"))
# barplot in the corner
par(xpd = FALSE, las = 1, mar = c(4, 4, 0, 0.5), fg = "gray30", col.axis = "gray30",
col.lab = "gray30")
mids <- barplot(t(means), beside = TRUE, col = "gray", border = "white", ylim = c(0.03,
0.06), xpd = FALSE, axes = FALSE, xaxt = "n")
abline(h = seq(0.03, 0.06, 0.005), col = "white")
arrows(c(t(mids)), unlist(means - ses), c(t(mids)), unlist(means + ses), length = 0.03,
code = 3, angle = 90, col = "gray40", lwd = 2)
par(las = 1)
axis(2, at = c(0.03, 0.04, 0.05, 0.06), cex.axis = 1, tck = 0.02, hadj = 0.5,
col = "gray30")
par(las = 0, lheight = 0.8, xpd = TRUE)
segments(c(1, 4, 7), rep(0.03, 3), c(3, 6, 9), rep(0.03, 3), col = "gray40")
mtext("Canberra Similarity", side = 2, line = 2.2, cex = 0.8)
mtext(" self other", side = 1, at = mean(mids[c(1, 2)]), cex = 0.7, line = -0)
mtext(" self other", side = 1, at = mean(mids[c(3, 4)]), cex = 0.7, line = -0)
mtext(" self other", side = 1, at = mean(mids[c(5, 6)]), cex = 0.7, line = -0)
mtext("All\nSamples", side = 1, at = mean(mids[c(1, 2)]), cex = 0.8, line = 2.2,
font = 2)
mtext("Air\nFilters", side = 1, at = mean(mids[c(3, 4)]), cex = 0.8, line = 2.2,
font = 2)
mtext("Settling\nDishes", side = 1, at = mean(mids[c(5, 6)]), cex = 0.8, line = 2.2,
font = 2)
# dev.off()
This does the same thing as the corner barplot above, but does it for both sets.
# pdf('../../figures/similarity_barplots.pdf', width=7, height=4)
pb.mat <- as.matrix(pb.can)
pb.matl <- pb.mat
pb.matl[upper.tri(pb.matl)] <- 0
s1 <- c(s14of, s14op)
s2 <- c(s24of, s24op)
s3 <- c(s34of, s34op)
s1f <- s14of
s2f <- s24of
s3f <- s34of
s1p <- s14op
s2p <- s24op
s3p <- s34op
self <- c(pb.matl[s1, s1], pb.matl[s2, s2], pb.matl[s3, s3])
other <- c(pb.mat[s1, c(s2, s3)], pb.mat[s2, c(s1, s3)], pb.mat[s3, c(s1, s2)])
selff <- c(pb.matl[s1f, s1f], pb.matl[s2f, s2f], pb.matl[s3f, s3f])
otherf <- c(pb.mat[s1f, c(s2f, s3f)], pb.mat[s2f, c(s1f, s3f)], pb.mat[s3f,
c(s1f, s2f)])
selfp <- c(pb.matl[s1p, s1p], pb.matl[s2p, s2p], pb.matl[s3p, s3p])
otherp <- c(pb.mat[s1p, c(s2p, s3p)], pb.mat[s2p, c(s1p, s3p)], pb.mat[s3p,
c(s1p, s2p)])
self <- self[-c(which(self == 0))]
selff <- selff[-c(which(selff == 0))]
selfp <- selfp[-c(which(selfp == 0))]
self <- 1 - self
selff <- 1 - selff
selfp <- 1 - selfp
other <- 1 - other
otherf <- 1 - otherf
otherp <- 1 - otherp
length(self)
## [1] 459
length(other)
## [1] 1944
length(selff)
## [1] 198
length(otherf)
## [1] 864
length(selfp)
## [1] 45
length(otherp)
## [1] 216
means <- data.frame(self = c(mean(self), mean(selff), mean(selfp)), other = c(mean(other),
mean(otherf), mean(otherp)), row.names = c("all", "filter", "petri"))
ses <- data.frame(self = c(se(self), se(selff), se(selfp)), other = c(se(other),
se(otherf), se(otherp)), row.names = c("all", "filter", "petri"))
# barplot
par(mfrow = c(1, 2))
par(xpd = FALSE, las = 1, mar = c(3, 4, 3, 0.5), fg = "gray30", col.axis = "gray30",
col.lab = "gray30")
mids <- barplot(t(means), beside = TRUE, col = "gray", border = "white", ylim = c(0.01,
0.06), xpd = FALSE, axes = FALSE, xaxt = "n")
abline(h = seq(0.01, 0.06, 0.005), col = "white")
arrows(c(t(mids)), unlist(means - ses), c(t(mids)), unlist(means + ses), length = 0.03,
code = 3, angle = 90, col = "gray40", lwd = 2)
par(las = 1)
axis(2, at = c(0.01, 0.02, 0.03, 0.04, 0.05, 0.06), cex.axis = 1, tck = 0.02,
hadj = 0.5, col = "gray30")
par(las = 0, lheight = 0.8, xpd = TRUE)
segments(c(1, 4, 7), rep(0.01, 3), c(3, 6, 9), rep(0.01, 3), col = "gray40")
mtext("Canberra Similarity", side = 2, line = 2.2, cex = 0.8)
mtext(" self other", side = 1, at = mean(mids[c(1, 2)]), cex = 0.7, line = 0)
mtext(" self other", side = 1, at = mean(mids[c(3, 4)]), cex = 0.7, line = 0)
mtext(" self other", side = 1, at = mean(mids[c(5, 6)]), cex = 0.7, line = 0)
mtext("All\nSamples", side = 1, at = mean(mids[c(1, 2)]), cex = 0.8, line = 1.7,
font = 2)
mtext("Air\nFilters", side = 1, at = mean(mids[c(3, 4)]), cex = 0.8, line = 1.7,
font = 2)
mtext("Settling\nDishes", side = 1, at = mean(mids[c(5, 6)]), cex = 0.8, line = 1.7,
font = 2)
mtext("Occupied", font = 2, line = 1)
pb.mat <- as.matrix(pb.can)
pb.matl <- pb.mat
pb.matl[upper.tri(pb.matl)] <- 0
s1 <- c(s14uf, s14up)
s2 <- c(s24uf, s24up)
s3 <- c(s34uf, s34up)
s1f <- s14uf
s2f <- s24uf
s3f <- s34uf
s1p <- s14up
s2p <- s24up
s3p <- s34up
self <- c(pb.matl[s1, s1], pb.matl[s2, s2], pb.matl[s3, s3])
other <- c(pb.mat[s1, c(s2, s3)], pb.mat[s2, c(s1, s3)], pb.mat[s3, c(s1, s2)])
selff <- c(pb.matl[s1f, s1f], pb.matl[s2f, s2f], pb.matl[s3f, s3f])
otherf <- c(pb.mat[s1f, c(s2f, s3f)], pb.mat[s2f, c(s1f, s3f)], pb.mat[s3f,
c(s1f, s2f)])
selfp <- c(pb.matl[s1p, s1p], pb.matl[s2p, s2p], pb.matl[s3p, s3p])
otherp <- c(pb.mat[s1p, c(s2p, s3p)], pb.mat[s2p, c(s1p, s3p)], pb.mat[s3p,
c(s1p, s2p)])
self <- self[-c(which(self == 0))]
selff <- selff[-c(which(selff == 0))]
selfp <- selfp[-c(which(selfp == 0))]
self <- 1 - self
selff <- 1 - selff
selfp <- 1 - selfp
other <- 1 - other
otherf <- 1 - otherf
otherp <- 1 - otherp
length(self)
## [1] 426
length(other)
## [1] 1800
length(selff)
## [1] 187
length(otherf)
## [1] 816
length(selfp)
## [1] 40
length(otherp)
## [1] 192
means <- data.frame(self = c(mean(self), mean(selff), mean(selfp)), other = c(mean(other),
mean(otherf), mean(otherp)), row.names = c("all", "filter", "petri"))
ses <- data.frame(self = c(se(self), se(selff), se(selfp)), other = c(se(other),
se(otherf), se(otherp)), row.names = c("all", "filter", "petri"))
# barplot
par(xpd = FALSE, las = 1, mar = c(3, 4, 3, 0.5), fg = "gray30", col.axis = "gray30",
col.lab = "gray30")
mids <- barplot(t(means), beside = TRUE, col = "gray", border = "white", ylim = c(0.01,
0.06), xpd = FALSE, axes = FALSE, xaxt = "n")
abline(h = seq(0.01, 0.06, 0.005), col = "white")
arrows(c(t(mids)), unlist(means - ses), c(t(mids)), unlist(means + ses), length = 0.03,
code = 3, angle = 90, col = "gray40", lwd = 2)
par(las = 1)
axis(2, at = c(0.01, 0.02, 0.03, 0.04, 0.05, 0.06), cex.axis = 1, tck = 0.02,
hadj = 0.5, col = "gray30")
par(las = 0, lheight = 0.8, xpd = TRUE)
segments(c(1, 4, 7), rep(0.01, 3), c(3, 6, 9), rep(0.01, 3), col = "gray40")
# mtext('Canberra Similarity', side=2, line=2.2, cex=.8)
mtext(" self other", side = 1, at = mean(mids[c(1, 2)]), cex = 0.7, line = 0)
mtext(" self other", side = 1, at = mean(mids[c(3, 4)]), cex = 0.7, line = 0)
mtext(" self other", side = 1, at = mean(mids[c(5, 6)]), cex = 0.7, line = 0)
mtext("All\nSamples", side = 1, at = mean(mids[c(1, 2)]), cex = 0.8, line = 1.7,
font = 2)
mtext("Air\nFilters", side = 1, at = mean(mids[c(3, 4)]), cex = 0.8, line = 1.7,
font = 2)
mtext("Settling\nDishes", side = 1, at = mean(mids[c(5, 6)]), cex = 0.8, line = 1.7,
font = 2)
mtext("Unoccupied", font = 2, line = 1)
# dev.off()
taxon_analysis.RData
.pb.order.abs <- pb.order[, c(2, 5:10, 1, 3, 4)]
all.means <- aggregate(pb.order.abs, by = list(pb.order$gr.num), FUN = "mean")
all.ses <- aggregate(pb.order.abs, by = list(pb.order$gr.num), FUN = se)
# index of OTUs for table and rep_set blasting
names(all.means)
## [1] "Group.1" "25506" "94008" "80526" "72518" "129988" "199152"
## [8] "41465" "234992" "121358" "201443"
all.means$Group.1
## [1] "s14of" "s14op" "s14uf" "s14up" "s24of" "s24op" "s24uf" "s24up"
## [9] "s34of" "s34op" "s34uf" "s34up"
occ.f <- rev(c(1, 5, 9))
unocc.f <- rev(c(3, 7, 11))
occ.p <- rev(c(2, 6, 10))
unocc.p <- rev(c(4, 8, 12))
# make individual datasets of means and sterrors for barplots
mean.f.names <- NULL
se.f.names <- NULL
mean.p.names <- NULL
se.p.names <- NULL
for (i in 2:ncol(all.means)) {
id <- names(all.means)[i]
mean.df.id <- paste("mean.f.", id, sep = "")
mean.df <- data.frame(unocc = all.means[unocc.f, i], occ = all.means[occ.f,
i])
row.names(mean.df) <- c("s3", "s2", "s1")
assign(mean.df.id, mean.df)
mean.f.names <- c(mean.f.names, mean.df.id)
se.df.id <- paste("se.f.", id, sep = "")
se.df <- data.frame(unocc = all.ses[unocc.f, i], occ = all.ses[occ.f, i])
row.names(se.df) <- c("s3", "s2", "s1")
assign(se.df.id, se.df)
se.f.names <- c(se.f.names, se.df.id)
mean.df.id <- paste("mean.p.", id, sep = "")
mean.df <- data.frame(unocc = all.means[unocc.p, i], occ = all.means[occ.p,
i])
row.names(mean.df) <- c("s3", "s2", "s1")
assign(mean.df.id, mean.df)
mean.p.names <- c(mean.p.names, mean.df.id)
se.df.id <- paste("se.p.", id, sep = "")
se.df <- data.frame(unocc = all.ses[unocc.p, i], occ = all.ses[occ.p, i])
row.names(se.df) <- c("s3", "s2", "s1")
assign(se.df.id, se.df)
se.p.names <- c(se.p.names, se.df.id)
}
text for left margin
martext <- c("Staphylococcus epidermidis", "Corynebacterium tuberculostearicum",
"Corynebacterium mucifaciens", "Corynebacterium amycolatum", "Corynebacterium jeikeium",
"Lactobacillus crispatus", "Dolosigranulum pigrum", "Methylobacterium jeotgali",
"Tumebacillus flagellatus", "Stenotrophomonas maltophilia")
martext2 <- c("NR_074995.1 (100% similar)", "NR_028975.1 (100% similar)", "NR_026396.1 (100% similar)",
"NR_026215.1 (100% similar)", "NR_074706.1 (100% similar)", "NR_074986.1 (100% similar)",
"NR_026098.1 (99% similar)", "NR_043878.1 (100% similar)", "NR_109600.1 (92% similar)",
"NR_074875.1 (100% similar)")
# order datasets top to bottom
mean.order <- mean.f.names
se.order <- se.f.names
len <- length(mean.order)
lim <- c(0, 0.5)
# pdf('../../figures/otu_boxplots.pdf', height=8, width=8.5)
# first plot with filters
layout(matrix(c(1:(2 * len)), len, 2), heights = c(1.25, rep(1, len - 2), 1.5),
widths = c(1.8, 1))
for (i in 1:len) {
means <- get(mean.order[i])/1000
ses <- get(se.order[i])/1000
segoff <- 0
occcex <- 0.8
if (i == 1) {
par(mar = c(1, 20, 2.5, 1), xpd = TRUE)
} else if (i == len) {
par(mar = c(4, 20, 1, 1), xpd = TRUE)
} else {
par(mar = c(1, 20, 1, 1), xpd = TRUE)
}
mids <- barplot(as.matrix(means), beside = TRUE, las = 1, border = "gray80",
xlab = "", axisnames = FALSE, horiz = TRUE, col = c(s3.col, s2.col,
s1.col), xlim = lim, xaxt = "n")
abline(v = seq(0.05, 1, 0.05), col = "white")
arrows(unlist(c(means - ses)), unlist(c(mids)), unlist(c(means + ses)),
unlist(c(mids)), code = 3, angle = 90, length = 0.01, col = "gray50")
segments(rep(segoff, 2), c(1, 5), rep(segoff, 2), c(4, 8), col = "gray60",
lend = "square", lwd = 1)
mtext(c("Un", "Occ"), side = 2, at = mids[2, ], col = "gray30", las = 1,
line = 0.2, cex = occcex)
mtext(c(martext[i], martext2[i]), side = 2, line = 3.5, las = 1, cex = 0.7,
at = c(mids[1, 2], mids[3, 1]), font = c(4, 1), col = "gray30")
# mtext(c(gsub('mean.f.', '', mean.order[i]), consensus[gsub('mean.f.',
# '', mean.order[i])]), side=2, line=3.5, las=1, cex=.7, at=c(mids[1,2],
# mids[3,1]), font=c(4,1), col='gray30')
segments(-0.07, mids[1, 1], -0.07, mids[3, 2], col = "gray")
if (i == 1) {
mtext("(a) Air Filters", side = 3, adj = 0, line = 1.2, col = "gray30",
font = 2)
}
}
## Warning: zero-length arrow is of indeterminate angle and so skipped
## Warning: zero-length arrow is of indeterminate angle and so skipped
## Warning: zero-length arrow is of indeterminate angle and so skipped
par(col.axis = "gray40")
axis(1, at = c(seq(0, 0.5, 0.1)), col.ticks = "gray40", col = "gray40", tck = -0.08,
labels = c("0", "0.10", "0.20", "0.30", "0.40", "0.50"))
axis(1, at = c(seq(0.05, 0.45, 0.1)), col.ticks = "gray40", col = "transparent",
tck = -0.03, labels = FALSE)
mtext("Relative Abundance", side = 1, line = 2.5, cex = 0.8, col = "gray30")
############## petris
mean.order <- mean.p.names
se.order <- se.p.names
for (i in 1:len) {
means <- get(mean.order[i])/1000
ses <- get(se.order[i])/1000
segoff <- 0
occcex <- 0.8
if (i == 1) {
par(mar = c(1, 1, 2.5, 1), xpd = TRUE)
} else if (i == len) {
par(mar = c(4, 1, 1, 1), xpd = TRUE)
} else {
par(mar = c(1, 1, 1, 1), xpd = TRUE)
}
mids <- barplot(as.matrix(means), beside = TRUE, las = 1, border = "gray80",
xlab = "", axisnames = FALSE, horiz = TRUE, col = c(s3.col, s2.col,
s1.col), xlim = lim, xaxt = "n")
abline(v = seq(0.05, 1, 0.05), col = "white")
arrows(unlist(c(means - ses)), unlist(c(mids)), unlist(c(means + ses)),
unlist(c(mids)), code = 3, angle = 90, length = 0.01, col = "gray50")
segments(rep(segoff, 2), c(1, 5), rep(segoff, 2), c(4, 8), col = "gray60",
lend = "square", lwd = 1)
mtext(c("Un", "Occ"), side = 2, at = mids[2, ], col = "gray30", las = 1,
line = 0.2, cex = occcex)
if (i == 1) {
mtext("(b) Settling Dishes", side = 3, adj = 0, line = 1.2, col = "gray30",
font = 2)
}
}
## Warning: zero-length arrow is of indeterminate angle and so skipped
## Warning: zero-length arrow is of indeterminate angle and so skipped
## Warning: zero-length arrow is of indeterminate angle and so skipped
## Warning: zero-length arrow is of indeterminate angle and so skipped
## Warning: zero-length arrow is of indeterminate angle and so skipped
par(col.axis = "gray40")
axis(1, at = c(seq(0, 0.5, 0.1)), col.ticks = "gray40", col = "gray40", tck = -0.08,
labels = c("0", "0.10", "0.20", "0.30", "0.40", "0.50"))
axis(1, at = c(seq(0.05, 0.45, 0.1)), col.ticks = "gray40", col = "transparent",
tck = -0.03, labels = FALSE)
mtext("Relative Abundance", side = 1, line = 2.5, cex = 0.8, col = "gray30")
par(fg = "gray30")
legend("right", legend = c("Subject 1", "Subject 2", "Subject 3"), bty = "n",
pch = c(22, 22, 22), col = c(rep("gray70", 3)), pt.bg = c(s1.col, s2.col,
s3.col), cex = 1, pt.cex = 2.3, y.intersp = 0.8)
# dev.off()