Personal Microbial Cloud

Beta-diversity Data Analysis

Read in prepared workspace and load packages.

load("../data/pb_data.RData")
library(vegan)
## Loading required package: permute This is vegan 2.0-8
library(labdsv)
## Loading required package: mgcv This is mgcv 1.7-22. For overview type
## 'help("mgcv-package")'. Loading required package: MASS
## 
## Attaching package: 'labdsv'
## 
## The following object is masked from 'package:stats':
## 
## density
library(ape)
library(xtable)
source("../data/pb_functions_for_manuscript.R")

Unpack list of groups.

for (i in 1:length(names(groups))) {
    assign(names(groups)[i], unlist(groups[[i]]))
}
rm(i)

How do samples fall into groups? I.e., do we have balance in treatments?

table(pb.map$person, pb.map$sampleType)
##     
##      filter petri.dish
##   s1     48         24
##   s2     46         23
##   s3     46         24
table(pb.map$location)
## 
##   occ unocc 
##   108   103
table(pb.map$duration)
## 
##   2   4 
## 105 106
table(pb.map$sampleType)
## 
##     filter petri.dish 
##        140         71

The canberra distance metric is used to calculate beta diversity. Then 2d NMDS objects are created for the whole experiment and for each treatment individually.

pb.can <- vegdist(pb.1000, "canberra")

Since most abundant taxa were removed as contaminants, try Bray-Curtis.

pb.bc <- vegdist(pb.1000, "bray")
# Filters 4h diff people
occ4f <- c(s24of, s34of, s14of, s24uf, s34uf, s14uf)
occ4f.table <- pb.1000[occ4f, ]
occ4f.map <- pb.map[occ4f, ]
occ4f.can <- as.dist(as.matrix(pb.can)[occ4f, occ4f])
occ4f.nmds <- bestnmds(occ4f.can)
occ <- which(occ4f.map$location == "occ")
unocc <- which(occ4f.map$location == "unocc")

## 2h diff people
occ2f <- c(s22of, s32of, s12of, s22uf, s32uf, s12uf)
occ2f.table <- pb.1000[occ2f, ]
occ2f.map <- pb.map[occ2f, ]
occ2f.can <- as.dist(as.matrix(pb.can)[occ2f, occ2f])
occ2f.nmds <- bestnmds(occ2f.can)
occ2 <- which(occ2f.map$location == "occ")
unocc2 <- which(occ2f.map$location == "unocc")

# Petri dishes 4h diff people
occ4p <- c(s24op, s34op, s14op, s24up, s34up, s14up)
occ4p.table <- pb.1000[occ4p, ]
occ4p.map <- pb.map[occ4p, ]
occ4p.can <- as.dist(as.matrix(pb.can)[occ4p, occ4p])
occ4p.nmds <- bestnmds(occ4p.can)
occ.p <- which(occ4p.map$location == "occ")
unocc.p <- which(occ4p.map$location == "unocc")

## 2h diff people
occ2p <- c(s22op, s32op, s12op, s22up, s32up, s12up)
occ2p.table <- pb.1000[occ2p, ]
occ2p.map <- pb.map[occ2p, ]
occ2p.can <- as.dist(as.matrix(pb.can)[occ2p, occ2p])
occ2p.nmds <- bestnmds(occ2p.can)
occ2.p <- which(occ2p.map$location == "occ")
unocc2.p <- which(occ2p.map$location == "unocc")

Pretty much all of these NMDS ordinations have high stress levels:


nmds.stress <- data.frame(treatment = c("4hr filters", "2hr filters", "4hr petri dishes", 
    "2hr petri dishes"), stress = c(occ4f.nmds$stress, occ2f.nmds$stress, occ4p.nmds$stress, 
    occ2p.nmds$stress), n = c(length(occ4f), length(occ2f), length(occ4p), length(occ2p)))
print(xtable(nmds.stress), type = "html")
treatment stress n
1 4hr filters 27.91 71
2 2hr filters 30.94 69
3 4hr petri dishes 26.00 35
4 2hr petri dishes 25.73 36

But since statistics back up the visible patterns, it doesn't present much of a problem.

The functions file (pb_functions_for_mansucript.R) contains a simple function to calculate richness, Shannon-Weiner diversity, and Evenness for each sample.

pb.hrj <- Evenness(pb.1000)

occ4f.map$R <- pb.hrj[occ4f, "R"]
occ4f.map$H1 <- pb.hrj[occ4f, "H1"]

occ2f.map$R <- pb.hrj[occ2f, "R"]
occ2f.map$H1 <- pb.hrj[occ2f, "H1"]

occ4p.map$R <- pb.hrj[occ4p, "R"]
occ4p.map$H1 <- pb.hrj[occ4p, "H1"]

occ2p.map$R <- pb.hrj[occ2p, "R"]
occ2p.map$H1 <- pb.hrj[occ2p, "H1"]

One last piece to assess the spread of points with beta-dispersion tests. It appears that s1 and s3 points are tightly grouped but s2 is not - this gives us a metric to report.

# 4of
pb.onlyocc.can <- as.dist(as.matrix(pb.can)[c(s24of, s34of, s14of), c(s24of, 
    s34of, s14of)])
pb.onlyocc.map <- pb.map[c(s24of, s34of, s14of), ]
occ4f.bd <- betadisper(pb.onlyocc.can, pb.onlyocc.map$person)
occ4f.betad <- data.frame(dists = occ4f.bd$distances, person = pb.onlyocc.map$person)
occ4f.betad.aov <- anova(occ4f.bd)
# boxplot(occ4f.bd$distances~pb.onlyocc.map$person)

# 2of
pb.onlyocc2.can <- as.dist(as.matrix(pb.can)[c(s22of, s32of, s12of), c(s22of, 
    s32of, s12of)])
pb.onlyocc2.map <- pb.map[c(s22of, s32of, s12of), ]
occ2f.bd <- betadisper(pb.onlyocc2.can, pb.onlyocc2.map$person)
occ2f.betad <- data.frame(dists = occ2f.bd$distances, person = pb.onlyocc2.map$person)
occ2f.betad.aov <- anova(occ2f.bd)
# boxplot(occ2f.bd$distances~pb.onlyocc2.map$person)

# 4op
pb.onlyoccp.can <- as.dist(as.matrix(pb.can)[c(s24op, s34op, s14op), c(s24op, 
    s34op, s14op)])
pb.onlyoccp.map <- pb.map[c(s24op, s34op, s14op), ]
occ4p.bd <- betadisper(pb.onlyoccp.can, pb.onlyoccp.map$person)
occ4p.betad <- data.frame(dists = occ4p.bd$distances, person = pb.onlyoccp.map$person)
occ4p.betad.aov <- anova(occ4p.bd)
# boxplot(occ4p.bd$distances~pb.onlyoccp.map$person)

# 2of
pb.onlyocc2p.can <- as.dist(as.matrix(pb.can)[c(s22op, s32op, s12op), c(s22op, 
    s32op, s12op)])
pb.onlyocc2p.map <- pb.map[c(s22op, s32op, s12op), ]
occ2p.bd <- betadisper(pb.onlyocc2p.can, pb.onlyocc2p.map$person)
occ2p.betad <- data.frame(dists = occ2p.bd$distances, person = pb.onlyocc2p.map$person)
occ2p.betad.aov <- anova(occ2p.bd)
# boxplot(occ2p.bd$distances~pb.onlyocc2p.map$person)

betad.aov.list <- list(occ4f.betad.aov, occ2f.betad.aov, occ4p.betad.aov, occ2p.betad.aov)
betad.table.list <- list(occ4f.betad, occ2f.betad, occ4p.betad, occ2p.betad)

rm(pb.onlyocc.can, pb.onlyocc.map, occ4f.bd, pb.onlyocc2.can, pb.onlyocc2.map, 
    occ2f.bd, pb.onlyoccp.can, pb.onlyoccp.map, occ4p.bd, pb.onlyocc2p.can, 
    pb.onlyocc2p.map, occ2p.bd, occ4f.betad, occ4f.betad.aov, occ4p.betad, occ4p.betad.aov, 
    occ2f.betad, occ2f.betad.aov, occ2p.betad, occ2p.betad.aov)

After all of these data manipulations, some of these objects now go to a new script to make figures.

ls()
##  [1] "betad.aov.list"   "betad.table.list" "cons"            
##  [4] "Evenness"         "gr.name"          "gr2"             
##  [7] "groups"           "makeTaxo"         "nmds.stress"     
## [10] "occ"              "occ.p"            "occ2"            
## [13] "occ2.p"           "occ2f"            "occ2f.can"       
## [16] "occ2f.map"        "occ2f.nmds"       "occ2f.table"     
## [19] "occ2p"            "occ2p.can"        "occ2p.map"       
## [22] "occ2p.nmds"       "occ2p.table"      "occ4f"           
## [25] "occ4f.can"        "occ4f.map"        "occ4f.nmds"      
## [28] "occ4f.table"      "occ4p"            "occ4p.can"       
## [31] "occ4p.map"        "occ4p.nmds"       "occ4p.table"     
## [34] "pb.1000"          "pb.bc"            "pb.can"          
## [37] "pb.hrj"           "pb.map"           "rarefied.otus"   
## [40] "rarefied.samples" "rarefied.seqs"    "s12of"           
## [43] "s12op"            "s12uf"            "s12up"           
## [46] "s14of"            "s14op"            "s14uf"           
## [49] "s14up"            "s22of"            "s22op"           
## [52] "s22uf"            "s22up"            "s24of"           
## [55] "s24op"            "s24uf"            "s24up"           
## [58] "s32of"            "s32op"            "s32uf"           
## [61] "s32up"            "s34of"            "s34op"           
## [64] "s34uf"            "s34up"            "se"              
## [67] "taxo"             "total.otus"       "total.samples"   
## [70] "total.seqs"       "unocc"            "unocc.p"         
## [73] "unocc2"           "unocc2.p"

# clean groupings.
rm(list = ls()[which(ls() %in% names(groups))])
rm(occ, unocc, occ.p, unocc.p, occ2, unocc2, occ2.p, unocc2.p, occ2f, occ2p, 
    occ4f, occ4p)
ls()
##  [1] "betad.aov.list"   "betad.table.list" "cons"            
##  [4] "Evenness"         "gr.name"          "gr2"             
##  [7] "groups"           "makeTaxo"         "nmds.stress"     
## [10] "occ2f.can"        "occ2f.map"        "occ2f.nmds"      
## [13] "occ2f.table"      "occ2p.can"        "occ2p.map"       
## [16] "occ2p.nmds"       "occ2p.table"      "occ4f.can"       
## [19] "occ4f.map"        "occ4f.nmds"       "occ4f.table"     
## [22] "occ4p.can"        "occ4p.map"        "occ4p.nmds"      
## [25] "occ4p.table"      "pb.1000"          "pb.bc"           
## [28] "pb.can"           "pb.hrj"           "pb.map"          
## [31] "rarefied.otus"    "rarefied.samples" "rarefied.seqs"   
## [34] "se"               "taxo"             "total.otus"      
## [37] "total.samples"    "total.seqs"
# save the rest.
save.image(file = "../data/pb_analysis_to_figures.RData")
getwd()
## [1] "/Users/jfmeadow/Dropbox/pb_manuscript/manuscript_code/analysis"