Load data.
setwd("~/Dropbox/pb_shared_markdown/manuscript_code/contamination/")
library(vegan)
library(labdsv)
source("../data/pb_functions_for_manuscript.R")
pb <- read.delim("../data/otu_table_r.txt", head = TRUE, row.names = 1)
pb.tax <- pb[, ncol(pb)]
pb <- pb[, -ncol(pb)]
pb <- t(pb)
pb.map <- read.delim("../data/pb_map.txt", head = TRUE, row.names = 1)
all(row.names(pb.map) %in% row.names(pb))
## [1] TRUE
dim(pb.map)
## [1] 300 12
dim(pb)
## [1] 300 235688
pb.map <- pb.map[row.names(pb), ]
dim(pb)
## [1] 300 235688
dim(pb.map)
## [1] 300 12
length(pb.tax)
## [1] 235688
identical(row.names(pb), row.names(pb.map))
## [1] TRUE
Identify those not used in the study.
nonexp <- which(pb.map$location == "east" | pb.map$location == "west" | pb.map$location ==
"out")
Keep original copy
pb.original <- pb
pb.tax.original <- pb.tax
Take out plant sequences
plants <- grep("Streptophyt", pb.tax)
pb <- pb[, -plants]
pb.tax <- pb.tax[-plants] # remove from taxonomy to line up
Make a full rarefied dataset including controls and contaminants.
sort(rowSums(pb.original))
## uf04.2h.ab uf05.2h.ab uf09.2h.aa uf05.4h.aa cef08.4h uf07.2h.jm
## 14 41 723 781 1522 1825
## uf11.4h.jm up04.4h.aa up05.2h.aa extr.c3 of01.2h.aa uf09.4h.aa
## 1976 2520 2954 2980 3563 3624
## uf12.2h.ab of11.2h.jm extr.c4 of03.4h.aa cep06.4h of04.4h.jm
## 4248 4310 4314 4343 4436 4544
## extr.c7 uf06.4h.jm uf03.4h.ab cwp02.4h uf08.2h.aa of11.2h.ab
## 4805 4831 4932 5037 5191 5197
## neg.dish.4 uf08.2h.jm extr.c11 cef03.4h uf05.2h.jm of09.2h.aa
## 5260 5445 5452 5654 5697 6003
## out02.2h.ab uf08.4h.aa extr.c10 of02.2h.aa uf09.2h.ab up04.2h.aa
## 6241 6246 6253 6266 6389 6450
## uf03.2h.aa extr.c1 cef11.4h cwf01.4h cwp01.4h extr.c6
## 6487 6519 6589 6717 6804 6920
## op05.2h.jm uf03.2h.jm cef09.4h uf06.4h.aa op06.4h.jm up05.2h.ab
## 6931 7057 7065 7139 7214 7247
## uf09.2h.jm uf12.2h.jm uf03.4h.aa neg.swab.2 of07.4h.aa uf12.4h.aa
## 7420 7594 7691 7873 7885 7897
## uf01.2h.ab cwf12.4h cwf10.4h extr.c12 uf03.4h.jm neg.swab.3
## 7945 7995 8132 8474 8581 8730
## cwf02.4h neg.swab.4 up03.4h.ab of05.2h.aa neg.filt.4 cef12.4h
## 8747 8772 8820 8980 9038 9241
## cep01.4h up04.4h.ab of09.4h.aa of06.2h.ab uf10.4h.ab up02.4h.ab
## 9539 9593 9798 9836 10006 10009
## out01.2h.aa up02.2h.aa cep02.4h op03.2h.ab neg.filt.3 op06.2h.ab
## 10018 10025 10199 10245 10533 10967
## up06.4h.jm cwf07.4h uf07.2h.ab uf01.4h.jm op04.4h.aa of04.4h.aa
## 11058 11096 11201 11275 11413 11651
## of03.2h.aa uf11.2h.jm uf03.2h.ab neg.dish.3 uf10.4h.aa cwp04.4h
## 11746 11792 11877 12029 12067 12342
## extr.c5 out03.2h.jm uf02.2h.ab up05.4h.aa uf10.2h.ab up04.2h.jm
## 12369 12409 12525 12736 12855 12920
## extr.c15 op05.4h.ab up06.2h.jm of03.4h.jm extr.c13 uf02.4h.ab
## 13047 13199 13262 13266 13319 13667
## uf02.4h.jm uf06.2h.jm up06.4h.ab uf12.4h.jm of10.2h.jm up01.4h.ab
## 14003 14050 14219 14321 14336 14460
## of10.2h.ab of04.4h.ab uf09.4h.jm uf08.4h.jm cep04.4h uf04.4h.jm
## 14484 14562 14784 14875 15011 15186
## cwf09.4h of12.4h.aa op05.2h.ab of03.2h.jm out01.4h.jm out01.2h.jm
## 15187 15593 15660 15663 15995 16003
## uf11.4h.aa op04.2h.aa uf01.4h.aa out02.2h.jm of08.2h.aa uf08.2h.ab
## 16170 16334 16349 16354 16403 16425
## of02.2h.ab of04.2h.aa op06.2h.jm uf01.2h.jm of12.2h.ab of10.4h.ab
## 16627 16801 16933 16960 17098 17370
## cwf11.4h cef02.4h up03.4h.jm up01.2h.jm up06.4h.aa cwp06.4h
## 17528 17584 17988 18043 18173 18405
## neg.filt.2 of03.2h.ab out01.2h.ab uf09.4h.ab uf04.4h.ab neg.dish.2
## 18475 19093 19358 19378 19448 19501
## of03.4h.ab op03.4h.aa up03.2h.jm op05.2h.aa neg.dish.1 up01.4h.aa
## 19629 19834 19851 19923 20392 20597
## uf10.2h.jm up03.2h.ab uf07.4h.jm extr.c14 op03.4h.ab uf10.2h.aa
## 20620 20753 20805 21268 21729 21865
## op04.2h.jm uf07.2h.aa uf12.4h.ab of06.2h.jm up05.4h.ab of01.2h.ab
## 22021 22056 22110 22347 22419 22473
## of09.2h.ab of08.4h.aa uf11.2h.aa op04.4h.ab up03.2h.aa uf11.2h.ab
## 22768 22833 23313 23317 23579 23582
## extr.c9 cwf05.4h of07.2h.jm uf06.2h.ab neg.filt.1 uf04.2h.jm
## 23755 24256 24995 25019 25033 25113
## of06.2h.aa out01.4h.ab op01.4h.ab up02.2h.jm of11.4h.aa uf04.4h.aa
## 25960 26001 26460 27156 27305 27412
## of05.2h.ab out03.2h.ab of04.2h.jm of09.4h.ab of02.4h.aa of05.4h.aa
## 28213 28969 29019 29404 29841 30280
## up01.4h.jm neg.swab.1 out02.4h.jm extr.c8 of11.4h.ab op01.2h.ab
## 30325 30427 30760 31143 31996 32111
## cep05.4h uf01.2h.aa op02.2h.jm of05.4h.jm of12.4h.jm cef01.4h
## 32448 33241 34365 34840 35015 35256
## cwf04.4h cef05.4h uf07.4h.aa of06.4h.aa op03.2h.jm of05.2h.jm
## 35369 35792 36484 36688 36832 36981
## uf02.4h.aa op01.2h.aa of07.2h.ab op02.4h.jm of04.2h.ab uf10.4h.jm
## 37003 37791 38035 38617 39740 40604
## op01.4h.jm c.out.03.4h of05.4h.ab cef07.4h c.out.01.4h of02.4h.jm
## 40741 41508 41515 42017 42363 42711
## cep03.4h cef04.4h up05.2h.jm op04.4h.jm up02.2h.ab of12.2h.jm
## 42888 43308 43609 43673 43739 43743
## of12.2h.aa uf04.2h.aa op02.4h.aa up06.2h.ab of10.4h.jm up04.4h.jm
## 46188 47668 47853 47932 48189 48228
## up02.4h.aa of11.4h.jm of08.2h.jm out03.4h.aa out01.4h.aa uf02.2h.jm
## 48508 48928 49057 49166 50494 50814
## uf05.4h.jm uf12.2h.aa cwp03.4h uf06.2h.aa cwf06.4h op02.2h.aa
## 50884 51119 51157 51581 51733 51818
## of08.2h.ab out02.4h.aa uf05.2h.aa uf02.2h.aa up06.2h.aa of02.2h.jm
## 52636 52806 53026 53420 53830 54590
## up01.2h.ab out03.4h.ab of01.2h.jm of10.4h.aa of08.4h.ab cef10.4h
## 55121 56066 56205 58669 60012 60368
## up05.4h.jm of07.4h.jm uf05.4h.ab of10.2h.aa of02.4h.ab cwf08.4h
## 60551 61009 61050 61616 62069 62276
## cwp05.4h of07.2h.aa of11.2h.aa uf06.4h.ab op05.4h.jm of01.4h.jm
## 62606 63211 63559 64060 64454 65531
## op04.2h.ab of09.2h.jm of12.4h.ab op06.2h.aa uf01.4h.ab up02.4h.jm
## 68925 69222 71576 71693 73705 73726
## up03.4h.aa of06.4h.jm uf08.4h.ab c.out.02.4h up04.2h.ab uf11.4h.ab
## 73920 73954 74676 76285 78288 80533
## out03.4h.jm op06.4h.aa out02.2h.aa of07.4h.ab up01.2h.aa out02.4h.ab
## 80593 85142 85418 85503 87362 89350
## of09.4h.jm op01.4h.aa extr.c2 cef06.4h op01.2h.jm of08.4h.jm
## 89773 95472 100420 103197 104357 108328
## op03.4h.jm of01.4h.aa op03.2h.aa uf07.4h.ab op02.2h.ab out03.2h.aa
## 108847 113924 116025 117705 119490 121248
## op06.4h.ab op05.4h.aa of06.4h.ab of01.4h.ab cwf03.4h op02.4h.ab
## 131677 142861 186832 205424 235568 311354
sort(rowSums(pb))
## uf04.2h.ab uf05.2h.ab uf09.2h.aa uf05.4h.aa cef08.4h uf07.2h.jm
## 14 41 723 781 1522 1825
## uf11.4h.jm up04.4h.aa up05.2h.aa extr.c3 of01.2h.aa uf09.4h.aa
## 1976 2520 2951 2954 3563 3624
## uf12.2h.ab of11.2h.jm extr.c4 of03.4h.aa cep06.4h of04.4h.jm
## 4248 4310 4314 4332 4436 4516
## extr.c7 uf06.4h.jm uf03.4h.ab cwp02.4h of11.2h.ab uf08.2h.aa
## 4805 4831 4932 5015 5188 5191
## neg.dish.4 uf08.2h.jm extr.c11 cef03.4h uf05.2h.jm of09.2h.aa
## 5256 5379 5452 5654 5692 5981
## of02.2h.aa uf08.4h.aa out02.2h.ab extr.c10 uf09.2h.ab uf03.2h.aa
## 6217 6238 6241 6253 6350 6391
## up04.2h.aa cef11.4h extr.c1 cwf01.4h op05.2h.jm cwp01.4h
## 6450 6471 6518 6717 6751 6804
## op06.4h.jm extr.c6 uf03.2h.jm cef09.4h uf06.4h.aa up05.2h.ab
## 6906 6920 7057 7057 7087 7247
## uf09.2h.jm uf12.2h.jm uf03.4h.aa neg.swab.2 of07.4h.aa uf12.4h.aa
## 7387 7508 7634 7873 7873 7897
## uf01.2h.ab cwf12.4h cwf10.4h extr.c12 neg.swab.4 uf03.4h.jm
## 7914 7995 8131 8474 8478 8518
## neg.swab.3 cwf02.4h up03.4h.ab of05.2h.aa neg.filt.4 cef12.4h
## 8730 8747 8819 8880 9038 9215
## cep01.4h up04.4h.ab of09.4h.aa of06.2h.ab out01.2h.aa uf10.4h.ab
## 9539 9576 9798 9834 9951 10005
## up02.4h.ab up02.2h.aa cep02.4h op03.2h.ab neg.filt.3 op06.2h.ab
## 10008 10025 10186 10227 10517 10921
## up06.4h.jm cwf07.4h uf07.2h.ab uf01.4h.jm op04.4h.aa of03.2h.aa
## 11058 11095 11160 11260 11413 11538
## uf03.2h.ab of04.4h.aa uf11.2h.jm neg.dish.3 uf10.4h.aa cwp04.4h
## 11545 11564 11792 12028 12067 12342
## extr.c5 out03.2h.jm uf02.2h.ab up01.4h.ab up05.4h.aa uf10.2h.ab
## 12369 12390 12525 12640 12675 12836
## up04.2h.jm extr.c15 of03.4h.jm op05.4h.ab up06.2h.jm extr.c13
## 12920 12988 13188 13198 13259 13294
## uf02.4h.ab uf02.4h.jm uf06.2h.jm up06.4h.ab of10.2h.jm uf12.4h.jm
## 13666 13937 14050 14219 14301 14321
## of10.2h.ab of04.4h.ab uf09.4h.jm uf08.4h.jm cep04.4h uf04.4h.jm
## 14456 14538 14640 14784 15010 15116
## cwf09.4h out01.4h.jm of12.4h.aa op05.2h.ab of03.2h.jm out01.2h.jm
## 15186 15417 15568 15660 15663 16003
## uf11.4h.aa out02.2h.jm op04.2h.aa uf01.4h.aa of08.2h.aa uf08.2h.ab
## 16104 16221 16316 16349 16363 16425
## of02.2h.ab of04.2h.aa uf01.2h.jm cwf11.4h op06.2h.jm of12.2h.ab
## 16601 16743 16834 16914 16932 17060
## of10.4h.ab cef02.4h up03.2h.jm up03.4h.jm up01.2h.jm up06.4h.aa
## 17255 17583 17865 17951 18005 18090
## cwp06.4h neg.filt.2 of03.2h.ab out01.2h.ab op05.2h.aa uf09.4h.ab
## 18403 18472 19093 19245 19319 19378
## uf04.4h.ab neg.dish.2 of03.4h.ab op03.4h.aa neg.dish.1 uf10.2h.jm
## 19448 19501 19628 19834 20297 20469
## up01.4h.aa up03.2h.ab uf07.4h.jm extr.c14 op03.4h.ab uf10.2h.aa
## 20531 20752 20805 21268 21727 21865
## uf12.4h.ab op04.2h.jm uf07.2h.aa of06.2h.jm of01.2h.ab up05.4h.ab
## 21908 21964 22021 22346 22366 22376
## of09.2h.ab of08.4h.aa uf11.2h.aa op04.4h.ab up03.2h.aa uf11.2h.ab
## 22706 22832 23312 23317 23365 23582
## extr.c9 cwf05.4h of07.2h.jm uf06.2h.ab neg.filt.1 uf04.2h.jm
## 23754 24255 24993 25018 25032 25112
## of06.2h.aa out01.4h.ab op01.4h.ab of11.4h.aa up02.2h.jm uf04.4h.aa
## 25959 26000 26460 27047 27156 27232
## of05.2h.ab out03.2h.ab of04.2h.jm of09.4h.ab of02.4h.aa out02.4h.jm
## 28157 28771 29018 29385 29825 30148
## of05.4h.aa up01.4h.jm neg.swab.1 extr.c8 of11.4h.ab op01.2h.ab
## 30228 30325 30427 31143 31996 32111
## cep05.4h of05.2h.jm uf01.2h.aa op02.2h.jm of05.4h.jm of12.4h.jm
## 32419 32433 33129 34332 34700 34943
## cef01.4h cwf04.4h cef05.4h uf02.4h.aa uf07.4h.aa of06.4h.aa
## 35256 35367 35521 36051 36484 36516
## op03.2h.jm op01.2h.aa of07.2h.ab op02.4h.jm of04.2h.ab uf10.4h.jm
## 36521 37789 37853 38542 39565 40457
## op01.4h.jm c.out.03.4h of05.4h.ab cef07.4h c.out.01.4h of02.4h.jm
## 40679 41086 41274 42017 42229 42694
## cep03.4h cef04.4h up05.2h.jm op04.4h.jm of12.2h.jm up02.2h.ab
## 42885 43112 43608 43672 43735 43739
## of12.2h.aa uf04.2h.aa op02.4h.aa of10.4h.jm up06.2h.ab up04.4h.jm
## 46181 47666 47810 47855 47931 48227
## up02.4h.aa of11.4h.jm of08.2h.jm out03.4h.aa out01.4h.aa uf05.4h.jm
## 48339 48927 49056 49157 50320 50777
## uf02.2h.jm uf12.2h.aa cwp03.4h uf06.2h.aa cwf06.4h op02.2h.aa
## 50813 51118 51128 51580 51608 51663
## of08.2h.ab out02.4h.aa uf05.2h.aa uf02.2h.aa up06.2h.aa of02.2h.jm
## 52635 52806 53026 53181 53829 54514
## up01.2h.ab out03.4h.ab of01.2h.jm of10.4h.aa up05.4h.jm of08.4h.ab
## 55121 55624 56067 58667 59347 59620
## cef10.4h of07.4h.jm uf05.4h.ab of10.2h.aa of02.4h.ab of07.2h.aa
## 60368 61008 61050 61256 61834 62152
## cwf08.4h cwp05.4h of11.2h.aa of01.4h.jm uf06.4h.ab op05.4h.jm
## 62274 62603 63384 64005 64060 64454
## op06.2h.aa op04.2h.ab of09.2h.jm of12.4h.ab uf01.4h.ab up02.4h.jm
## 68753 68922 69219 71576 73548 73725
## up03.4h.aa of06.4h.jm uf08.4h.ab c.out.02.4h up04.2h.ab out03.4h.jm
## 73919 73954 74232 76112 77925 79188
## uf11.4h.ab op06.4h.aa out02.2h.aa of07.4h.ab up01.2h.aa out02.4h.ab
## 80253 80564 84876 85444 87130 88968
## of09.4h.jm op01.4h.aa extr.c2 cef06.4h op01.2h.jm of08.4h.jm
## 89578 94784 100418 103076 104283 108327
## op03.4h.jm of01.4h.aa op03.2h.aa uf07.4h.ab op02.2h.ab out03.2h.aa
## 108694 113922 115487 117704 119462 120473
## op06.4h.ab op05.4h.aa of06.4h.ab of01.4h.ab cwf03.4h op02.4h.ab
## 131675 142713 186451 205422 235565 311310
pb.tmp <- pb[-nonexp, ]
pb.tmp <- pb.tmp[-which(rowSums(pb.tmp) < 3000), ]
pb.3500 <- rrarefy(pb.tmp, 3500)
pb.map <- pb.map[row.names(pb.3500), c("person", "location", "duration", "sampleType")]
dim(pb.3500)
## [1] 234 234213
dim(pb.map)
## [1] 234 4
identical(row.names(pb.3500), row.names(pb.map))
## [1] TRUE
Make taxonomy data frame for indexing.
taxo <- makeTaxo(pb.tax, pb.3500)
## Warning: no non-missing arguments to max; returning -Inf
head(taxo)
## kingdom phylum class order
## 0 Bacteria Bacteroidetes Sphingobacteria Sphingobacteriales
## 1 Bacteria Proteobacteria Gammaproteobacteria
## 2 Bacteria Proteobacteria Gammaproteobacteria Oceanospirillales
## 3 Bacteria Actinobacteria Actinobacteria Actinomycetales
## 4 Bacteria Actinobacteria Actinobacteria Actinomycetales
## 5 Bacteria
## family genus abundance
## 0 Chitinophagaceae 0
## 1 0
## 2 Halomonadaceae 1
## 3 0
## 4 1
## 5 0
Get rid of empty OTUs to reduce computing demand.
pb.3500 <- pb.3500[, -which(colSums(pb.3500) == 0)]
taxo <- taxo[colnames(pb.3500), ]
head(taxo)
## kingdom phylum class order
## 2 Bacteria Proteobacteria Gammaproteobacteria Oceanospirillales
## 4 Bacteria Actinobacteria Actinobacteria Actinomycetales
## 10 Bacteria Proteobacteria Gammaproteobacteria Oceanospirillales
## 15 Bacteria Proteobacteria Betaproteobacteria
## 16 Bacteria Firmicutes Bacilli Bacillales
## 19 Bacteria
## family genus abundance
## 2 Halomonadaceae 1
## 4 1
## 10 Halomonadaceae 2
## 15 1
## 16 Bacillaceae Tumebacillus 1
## 19 1
dim(taxo)
## [1] 33504 7
dim(pb.3500)
## [1] 234 33504
identical(row.names(taxo), colnames(pb.3500))
## [1] TRUE
pb.3500.ra <- pb.3500/35
Quick refrerence vector of names for figures and whatnot. Some genus, families, etc are empty, so search for the finest taxonomic level and use that name.
cons <- function(x) {
l <- length(x)
while (x[l] == "") {
l = l - 1
}
name <- x[l]
}
consensus <- apply(taxo[, 1:6], 1, cons)
consensus[1:10]
## 2 4 10
## "Halomonadaceae" "Actinomycetales" "Halomonadaceae"
## 15 16 19
## "Betaproteobacteria" "Tumebacillus" "Bacteria"
## 21 34 70
## "Bacteria" "Bacillaceae" "Proteobacteria"
## 79
## "Comamonadaceae"
Some indexing and metadata stuff for figures and such.
occ <- which(pb.map$location == "occ")
unocc <- which(pb.map$location == "unocc")
control <- which(pb.map$location == "control")
pb.map$pch <- 21
table(pb.map$sampleType)
##
## extraction.control filter petri.dish
## 12 142 74
## reagent.control swab
## 2 4
pb.map$pch2 <- pb.map$pch
pb.map$pch2[pb.map$sampleType == "extraction.control"] <- 24
pb.map$pch2[pb.map$sampleType == "reagent.control"] <- 23
pb.map$pch2[pb.map$sampleType == "petri.dish" & pb.map$location == "control"] <- 22
pb.map$bg <- "gray60"
pb.map$bg[occ] <- "cornflowerblue"
pb.map$bg[unocc] <- "darkorange"
Make dissimilarity matrix for big dataset - this takes a while. Make with both Canberra and Bray-Curtis, and then make NMDS ordinations. The Bray-Curtis will be more useful for identifying the influence of the most abundant taxa, while Canberra is used for analysis. Here, Canberra has high stress, so is not used much. Also quick plot to make sure that BC captures differences between groups.
can <- vegdist(pb.3500, "canberra")
bc <- vegdist(pb.3500)
nmds.can <- nmds(can)
## initial value 40.307806
## iter 5 value 32.361958
## iter 10 value 31.603218
## final value 31.348427
## converged
nmds.bc <- nmds(bc)
## initial value 24.728368
## iter 5 value 18.718472
## final value 18.288598
## converged
The BC ordination clearly shows that controls are tightly grouped (mostly), while occupied and unoccupied are defintely statistically different. Even with contaminants. Strong gradient along x-axis >> will be used for regression below.
# pdf('~/Desktop/contamination_ordination.pdf', useDingbats=TRUE)
plot(nmds.bc$points, col = pb.map$location, type = "n")
points(nmds.bc$points[occ, ], pch = 16, col = pb.map$bg[occ])
points(nmds.bc$points[unocc, ], pch = 16, col = pb.map$bg[unocc])
points(nmds.bc$points[control, ], pch = 21, bg = pb.map$bg[control])
# dev.off()
plot(nmds.bc$points, type = "n", xlim = c(-0.5, -0.3), ylim = c(-0.2, 0.1))
points(nmds.bc$points[occ, ], pch = 16, col = pb.map$bg[occ])
points(nmds.bc$points[unocc, ], pch = 16, col = pb.map$bg[unocc])
points(nmds.bc$points[control, ], pch = 21, bg = pb.map$bg[control])
Find the most abundant taxa. This is done for 10, 100, 1000. Use 10 for plotting, and take out those contaminants, but use 1000 since some contaminants are really rare in actual samples.
control.taxa.10 <- rev(sort(colSums(pb.3500[control, ])))[1:10]/sum(pb.3500[control,
])
occ.taxa.10 <- rev(sort(colSums(pb.3500[occ, ])))[1:10]/sum(pb.3500[occ, ])
unocc.taxa.10 <- rev(sort(colSums(pb.3500[unocc, ])))[1:10]/sum(pb.3500[unocc,
])
control.taxa.10["other"] <- 1 - sum(control.taxa.10)
occ.taxa.10["other"] <- 1 - sum(occ.taxa.10)
unocc.taxa.10["other"] <- 1 - sum(unocc.taxa.10)
top.10 <- cbind(rev(control.taxa.10), rev(unocc.taxa.10), rev(occ.taxa.10))
# mids <- barplot(top.10, col=c('gray30', rep('gray70', 10)),
# border='gray30', space=1)
con.cum <- cumsum(rev(control.taxa.10))
occ.cum <- cumsum(rev(occ.taxa.10))
unocc.cum <- cumsum(rev(unocc.taxa.10))
control.taxa.100 <- rev(sort(colSums(pb.3500[control, ])))[1:100]/sum(pb.3500[control,
])
occ.taxa.100 <- rev(sort(colSums(pb.3500[occ, ])))[1:100]/sum(pb.3500[occ, ])
unocc.taxa.100 <- rev(sort(colSums(pb.3500[unocc, ])))[1:100]/sum(pb.3500[unocc,
])
control.taxa.100["other"] <- 1 - sum(control.taxa.100)
occ.taxa.100["other"] <- 1 - sum(occ.taxa.100)
unocc.taxa.100["other"] <- 1 - sum(unocc.taxa.100)
con.cum.100 <- cumsum(rev(control.taxa.100))
occ.cum.100 <- cumsum(rev(occ.taxa.100))
unocc.cum.100 <- cumsum(rev(unocc.taxa.100))
occ.taxa.1000 <- rev(sort(colSums(pb.3500[occ, ])))[1:1000]/sum(pb.3500[occ,
])
unocc.taxa.1000 <- rev(sort(colSums(pb.3500[unocc, ])))[1:1000]/sum(pb.3500[unocc,
])
occ.taxa.1000["other"] <- 1 - sum(occ.taxa.1000)
unocc.taxa.1000["other"] <- 1 - sum(unocc.taxa.1000)
occ.cum.1000 <- cumsum(rev(occ.taxa.1000))
unocc.cum.1000 <- cumsum(rev(unocc.taxa.1000))
Halomonas is the most abundant across the board. It drives the BC NMDS.
rev(sort(colSums(pb.3500)))[1:10]
## 190873 234992 64285 25506 121358 201443 94008 80526 151122 72518
## 279828 66921 60454 45640 37840 28903 24596 20070 17032 12480
rev(sort(colSums(pb.3500)))[1]/sum(pb.3500) # 34%
## 190873
## 0.3417
halo <- "190873"
taxo["190873", ]
## kingdom phylum class order
## 190873 Bacteria Proteobacteria Gammaproteobacteria Oceanospirillales
## family genus abundance
## 190873 Halomonadaceae Halomonas 279828
pb.3500[, "190873"]
## op01.4h.jm extr.c10 uf10.4h.aa neg.dish.3 op01.4h.ab op06.2h.aa
## 867 2151 1609 2134 350 867
## of09.2h.ab of06.4h.ab op02.4h.ab uf07.4h.ab up05.4h.jm uf03.4h.ab
## 440 591 89 871 1028 1494
## op06.4h.ab uf05.2h.aa uf02.2h.jm op06.4h.aa of12.2h.aa up05.2h.ab
## 564 1621 1609 488 1040 2310
## of06.4h.jm extr.c2 of01.2h.ab of09.2h.aa uf04.4h.aa op02.2h.ab
## 1056 183 1358 619 1245 113
## of01.2h.jm uf01.2h.aa of05.2h.ab uf11.4h.ab uf08.2h.aa op04.2h.ab
## 638 1497 1281 1240 2236 889
## of03.2h.ab uf02.2h.aa neg.filt.1 of12.4h.ab of07.2h.ab uf10.4h.jm
## 1652 2016 1866 611 771 1270
## up02.4h.jm of08.4h.jm up06.2h.aa of05.4h.jm op05.4h.aa of11.4h.jm
## 1327 1097 378 1111 741 359
## of08.2h.jm up04.2h.ab extr.c8 uf04.2h.aa uf06.4h.aa of08.2h.aa
## 503 1163 2406 2041 685 744
## of01.4h.aa up01.2h.aa of07.4h.jm uf01.2h.jm of02.4h.jm of04.2h.jm
## 798 1263 681 1552 545 587
## up02.4h.aa of10.4h.aa uf03.4h.jm of01.4h.jm neg.dish.2 of06.2h.jm
## 2188 638 649 1113 2284 490
## op05.4h.jm of12.2h.jm uf05.4h.ab uf07.2h.ab op04.4h.jm op03.4h.jm
## 220 558 1426 1721 799 541
## of08.2h.ab of09.4h.jm uf07.4h.aa op01.4h.aa of07.4h.ab of01.4h.ab
## 1543 923 1976 136 259 351
## uf04.4h.ab uf06.2h.aa of02.4h.ab uf03.4h.aa of10.2h.aa of05.2h.jm
## 2189 1115 302 910 542 685
## up05.2h.jm up03.4h.aa uf03.2h.ab uf04.2h.jm of03.4h.ab of11.4h.aa
## 292 1069 1224 1955 1223 923
## op04.4h.aa op05.2h.ab of10.2h.ab op02.2h.jm of04.2h.ab of09.4h.aa
## 1047 1263 933 909 1780 1642
## op01.2h.jm op02.4h.aa uf05.2h.jm op03.2h.jm uf12.2h.ab uf01.4h.aa
## 864 118 1759 1227 1587 749
## op05.2h.aa op03.2h.aa of03.4h.jm uf02.4h.aa uf06.2h.ab of02.2h.jm
## 1647 655 369 974 2303 1031
## of09.2h.jm of10.4h.jm uf07.2h.aa uf01.4h.ab op03.4h.ab of08.4h.ab
## 953 623 967 1373 530 566
## extr.c15 op03.4h.aa op01.2h.aa of11.2h.aa up04.4h.jm of09.4h.ab
## 2030 905 1232 599 1658 321
## uf09.4h.jm uf08.4h.ab uf05.4h.jm uf08.4h.jm uf06.4h.ab up06.2h.ab
## 1062 1423 1372 636 2099 1605
## up04.4h.ab uf11.4h.aa neg.filt.2 op02.4h.jm uf03.2h.aa op02.2h.aa
## 1601 1582 2066 906 1503 664
## op03.2h.ab uf09.4h.ab neg.swab.3 up02.2h.jm of07.2h.jm up01.2h.ab
## 2039 550 2221 1175 553 1295
## up05.4h.ab of03.2h.aa extr.c9 op04.4h.ab of05.4h.aa up06.4h.aa
## 1529 1303 2098 1191 335 1279
## of06.2h.aa uf09.2h.ab extr.c14 up01.2h.jm up03.2h.jm up04.2h.aa
## 2018 1002 1690 1289 2181 828
## of04.2h.aa uf03.2h.jm op04.2h.jm uf01.4h.jm of05.4h.ab of04.4h.ab
## 976 1781 308 1708 670 615
## of12.4h.jm uf10.2h.ab uf02.4h.ab of11.4h.ab uf10.2h.aa of02.2h.aa
## 606 1365 432 281 914 1389
## up06.2h.jm up04.2h.jm up03.4h.jm uf12.2h.jm uf10.2h.jm of10.2h.jm
## 1003 1975 1191 864 1085 334
## op06.2h.jm uf07.4h.jm of02.4h.aa uf09.2h.jm uf02.4h.jm uf04.4h.jm
## 1052 1273 812 1222 1888 393
## uf08.4h.aa of08.4h.aa op05.4h.ab of07.2h.aa uf08.2h.jm uf11.2h.ab
## 1031 774 744 832 495 1037
## up03.4h.ab uf10.4h.ab up01.4h.jm of06.2h.ab uf12.4h.jm neg.swab.1
## 497 2067 997 1812 1559 1698
## up03.2h.aa up02.2h.ab of06.4h.aa op06.2h.ab up01.4h.aa neg.dish.1
## 1355 870 938 1259 1227 2586
## uf12.2h.aa of03.2h.jm extr.c6 up06.4h.jm op05.2h.jm op01.2h.ab
## 1175 960 1490 1218 846 659
## op04.2h.aa up02.2h.aa of04.4h.aa neg.filt.4 of01.2h.aa of11.2h.ab
## 395 2075 2804 2317 1148 1561
## up05.4h.aa uf08.2h.ab extr.c5 extr.c11 of10.4h.ab uf02.2h.ab
## 2135 1908 505 2191 937 2051
## extr.c12 up03.2h.ab extr.c1 uf12.4h.aa uf01.2h.ab of02.2h.ab
## 2092 1180 1315 2037 726 1684
## of03.4h.aa up01.4h.ab extr.c13 neg.swab.4 of04.4h.jm of11.2h.jm
## 739 1159 1649 2430 1220 404
## neg.swab.2 uf11.2h.jm up02.4h.ab of07.4h.aa uf09.4h.aa of12.2h.ab
## 2385 431 2036 859 1466 1401
## uf06.2h.jm op06.4h.jm neg.filt.3 of05.2h.aa uf12.4h.ab extr.c7
## 1440 790 2246 1583 2223 2306
## of12.4h.aa uf11.2h.aa up06.4h.ab uf06.4h.jm neg.dish.4 extr.c4
## 1490 1355 2191 552 2310 1241
haloCol <- pb.3500[, halo]
plot(haloCol ~ nmds.bc$points[, 1], type = "n")
points(haloCol[occ] ~ nmds.bc$points[occ, 1], pch = 16, col = pb.map$bg[occ])
points(haloCol[unocc] ~ nmds.bc$points[unocc, 1], pch = 16, col = pb.map$bg[unocc])
points(haloCol[control] ~ nmds.bc$points[control, 1], pch = 21, bg = pb.map$bg[control],
)
# identify(haloCol[control] ~ nmds.bc$points[control, 1],
# labels=pb.map$sampleType[control]) the mixed group of controls are all
# extraction blanks
Big panel figure showing each top contaminant in the dataset, and its influence on NMDS1
# pdf('~/Desktop/contaminationTop10.pdf', height=10, width=10,
# useDingbats=TRUE)
layout(matrix(c(11:15, 1:10, 16:20), 5, 4), heights = c(1, 1, 1, 1, 1.5), widths = c(1,
2, 2, 1))
par(las = 1, mar = c(0, 4, 0, 0), fg = "gray40", col.axis = "gray40", col.lab = "gray40")
for (i in 1:10) {
id <- names(control.taxa.10)[i]
print(taxo[id, ])
idCol <- pb.3500.ra[, id]
print(sum(idCol)/sum(pb.3500.ra))
y.lim <- round(range(idCol), 2)
ybuff <- y.lim[2]/8
x.lim <- round(range(nmds.bc$points[, 1]), 2)
xbuff <- sum(abs(x.lim))/8
if (i %in% c(1:4)) {
par(las = 1, mar = c(0, 4, 0, 0))
}
if (i %in% c(5:8)) {
par(las = 1, mar = c(0, 0, 0, 4))
}
if (i == 5) {
par(las = 1, mar = c(4, 4, 0, 0))
}
if (i == 10) {
par(las = 1, mar = c(4, 0, 0, 4))
}
plot(idCol ~ nmds.bc$points[, 1], ylim = c(y.lim + c(-ybuff, 2 * ybuff)),
xlim = c(x.lim + c(-xbuff, xbuff)), type = "n", xaxt = "n", yaxt = "n",
xlab = "", ylab = "")
if (i %in% c(1:5)) {
yax <- 2
} else {
yax <- 4
}
axis(side = yax, at = y.lim, labels = TRUE)
if (i %in% c(5, 10)) {
axis(side = 1, at = x.lim)
mtext("NMDS 1", side = 1, line = 1, cex = 0.7)
}
if (i == 1) {
mtext("relative\nabundance", side = 2, cex = 0.7, line = 0.5)
}
if (i == 6) {
mtext("relative\nabundance", side = 4, cex = 0.7, line = 0.5)
legend("right", legend = c("extraction control", "reagent control",
"petri control", "dish control", "occupied", "unoccupied"), pch = c(24,
23, 22, 21, 16, 16), col = c(1, 1, 1, 1, "cornflowerblue", "darkorange"),
pt.bg = rgb(0, 0, 0, 0.3))
}
points(idCol[occ] ~ nmds.bc$points[occ, 1], pch = 16, col = pb.map$bg[occ])
points(idCol[unocc] ~ nmds.bc$points[unocc, 1], pch = 16, col = pb.map$bg[unocc])
points(idCol[control] ~ nmds.bc$points[control, 1], pch = pb.map$pch2[control],
col = "black", bg = rgb(0, 0, 0, 0.3), cex = 2)
mtext(paste(consensus[id], " (", id, ")", sep = ""), line = -1.5)
}
## kingdom phylum class order
## 190873 Bacteria Proteobacteria Gammaproteobacteria Oceanospirillales
## family genus abundance
## 190873 Halomonadaceae Halomonas 279828
## [1] 0.3417
## kingdom phylum class order
## 64285 Bacteria Proteobacteria Gammaproteobacteria Alteromonadales
## family genus abundance
## 64285 Shewanellaceae Shewanella 60454
## [1] 0.07381
## kingdom phylum class order
## 234992 Bacteria Proteobacteria Alphaproteobacteria Rhizobiales
## family genus abundance
## 234992 Methylobacteriaceae Methylobacterium 66921
## [1] 0.08171
## kingdom phylum class order
## 201443 Bacteria Proteobacteria Gammaproteobacteria Xanthomonadales
## family genus abundance
## 201443 Xanthomonadaceae Stenotrophomonas 28903
## [1] 0.03529
## kingdom phylum class order family genus
## 121358 Bacteria Firmicutes Bacilli Bacillales Bacillaceae Tumebacillus
## abundance
## 121358 37840
## [1] 0.0462
## kingdom phylum class order
## 96443 Bacteria Proteobacteria Betaproteobacteria Burkholderiales
## family genus abundance
## 96443 Burkholderiaceae Burkholderia 3582
## [1] 0.004374
## kingdom phylum class order
## 129992 Bacteria Proteobacteria Betaproteobacteria Burkholderiales
## family genus abundance
## 129992 Comamonadaceae 5273
## [1] 0.006438
## kingdom phylum class order
## 151122 Bacteria Proteobacteria Gammaproteobacteria Enterobacteriales
## family genus abundance
## 151122 Enterobacteriaceae 17032
## [1] 0.0208
## kingdom phylum class order
## 180093 Bacteria Bacteroidetes Sphingobacteria Sphingobacteriales
## family genus abundance
## 180093 Cytophagaceae Hymenobacter 782
## [1] 0.0009548
## kingdom phylum class order
## 94008 Bacteria Actinobacteria Actinobacteria Actinomycetales
## family genus abundance
## 94008 Corynebacteriaceae Corynebacterium 24596
## [1] 0.03003
for (i in 11:2) {
xup <- xdown <- yup <- ydown <- NULL
yup <- rep(con.cum[i], 2)
ydown <- rep(con.cum[i - 1], 2)
xup <- xdown <- c(1, 2)
if (names(con.cum)[i] %in% names(unocc.cum.1000)) {
ucu <- which(names(unocc.cum.1000) == names(con.cum)[i])
yup <- c(yup, rep(unocc.cum.1000[ucu], 2))
ydown <- c(ydown, rep(unocc.cum.1000[ucu - 1], 2))
xup <- xdown <- c(xup, 3, 4)
}
if (names(con.cum)[i] %in% names(occ.cum.1000)) {
ocu <- which(names(occ.cum.1000) == names(con.cum)[i])
yup <- c(yup, rep(occ.cum.1000[ocu], 2))
ydown <- c(ydown, rep(occ.cum.1000[ocu - 1], 2))
xup <- xdown <- c(xup, 5, 6)
}
par(mar = c(3, 3, 3, 3))
if (i %in% c(7, 2)) {
par(mar = c(8, 3, 3, 3))
}
barplot(top.10, col = c("gray30", rep("gray93", 10)), border = "gray30",
space = 1, yaxt = "n")
if (i %in% c(7, 2)) {
par(las = 2)
mtext(c("control", "unoccupied", "occupied"), side = 1, at = c(1.5,
3.5, 5.5), cex = 0.8, line = 0.2)
}
par(las = 1)
if (!i %in% c(7, 2)) {
par(las = 2)
mtext(c("c", "u", "o"), side = 1, at = c(1.5, 3.5, 5.5), cex = 0.8,
line = 0.2)
}
par(las = 1)
if (i %in% 11:7) {
axis(2, at = c(0, 0.5, 1), labels = c(0, 50, 100))
}
if (i %in% 2:6) {
axis(4, at = c(0, 0.5, 1), labels = c(0, 50, 100))
}
polygon(c(xup, rev(xdown)), c(yup, rev(ydown)), col = rgb(0, 0, 0, 0.3))
}
# dev.off()
Now take out those most abundant in pcr controls.
dim(pb.tmp)
## [1] 234 234213
pb.c <- pb.tmp[-c(control), ]
pb.nc <- pb.tmp[-c(control), -which(colnames(pb.tmp) %in% names(control.taxa.10[c(1,
2, 6, 7)]))]
pb.map.nc <- pb.map[-control, ]
sort(rowSums(pb.nc))
## of04.4h.aa uf08.2h.aa up05.2h.ab uf09.4h.aa uf12.2h.ab of01.2h.aa
## 1035 1209 1264 1298 1832 1981
## uf05.2h.jm of11.2h.ab uf10.4h.ab uf12.4h.aa up02.2h.aa of04.4h.jm
## 2132 2313 2361 2411 2414 2563
## uf03.4h.ab uf03.2h.jm up02.4h.ab uf03.2h.aa op03.2h.ab up06.4h.ab
## 2577 2690 2725 2755 2780 2896
## of02.2h.aa of03.4h.aa up05.4h.aa of11.2h.jm up04.2h.jm uf02.2h.ab
## 3276 3339 3429 3544 3688 3716
## of06.2h.ab uf06.4h.jm of05.2h.aa uf09.2h.jm uf09.2h.ab up04.4h.ab
## 3770 3813 3929 3956 3974 4049
## uf08.4h.aa up03.2h.jm uf06.2h.ab uf01.4h.jm of09.4h.aa uf08.2h.jm
## 4254 4255 4340 4379 4398 4431
## of09.2h.aa up04.2h.aa op06.4h.jm op05.2h.jm uf03.4h.aa uf07.2h.ab
## 4633 4664 4743 4774 4858 5012
## uf12.4h.ab uf06.4h.aa uf10.4h.aa uf02.4h.jm uf12.2h.jm uf04.4h.ab
## 5403 5463 5465 5507 5553 5567
## of07.4h.aa op06.2h.ab of02.2h.ab up06.4h.jm uf01.2h.ab uf12.4h.jm
## 5695 5766 5859 6039 6122 6242
## of03.2h.aa uf03.4h.jm uf08.2h.ab op05.2h.aa uf06.2h.jm of12.4h.aa
## 6471 6544 6603 6699 6853 6966
## uf03.2h.ab up01.4h.ab op04.4h.aa uf10.2h.ab uf04.2h.jm up03.4h.ab
## 7024 7045 7058 7083 7331 7417
## uf11.4h.aa up06.2h.jm of03.2h.ab uf01.2h.jm of12.2h.ab of06.2h.aa
## 7708 8239 8321 8341 8599 8667
## op05.2h.ab uf09.4h.jm up01.2h.jm up06.4h.aa op05.4h.ab op06.2h.jm
## 8832 8924 8944 9123 9728 9836
## of10.2h.ab uf11.2h.jm of03.2h.jm up05.4h.ab up03.2h.aa of04.2h.aa
## 10098 10105 10287 10455 10719 10737
## up03.4h.jm uf07.4h.aa of04.4h.ab of03.4h.ab uf02.4h.ab of03.4h.jm
## 10874 10982 11224 11229 11361 11362
## up03.2h.ab of10.4h.ab of01.2h.ab uf08.4h.jm uf07.4h.jm up01.4h.aa
## 11409 11440 11496 11595 11680 11794
## op04.4h.ab uf01.4h.aa uf04.2h.aa of10.2h.jm of08.2h.aa up02.4h.aa
## 11943 12079 12268 12565 12570 12576
## uf04.4h.jm op03.4h.aa uf10.2h.jm uf11.2h.aa op04.2h.aa uf01.2h.aa
## 12827 12839 12993 13104 14264 14296
## uf11.2h.ab uf07.2h.aa uf04.4h.aa uf10.2h.aa of05.2h.ab uf09.4h.ab
## 14360 14393 15082 15325 15617 15760
## of04.2h.ab uf02.2h.aa uf06.4h.ab of08.4h.aa up02.2h.jm up01.4h.jm
## 15851 15868 16233 16497 16667 17823
## op03.4h.ab of11.4h.aa of06.2h.jm of09.2h.ab op04.2h.jm of07.2h.jm
## 17842 18619 18804 19083 19087 19664
## op03.2h.jm of02.4h.aa up04.4h.jm uf02.2h.jm of05.4h.jm of08.2h.ab
## 19949 20004 20114 20435 20546 21763
## up06.2h.ab uf05.4h.jm op01.2h.aa op02.2h.jm op01.4h.ab of04.2h.jm
## 21847 22022 22716 22824 23185 23190
## op01.2h.ab uf05.2h.aa uf10.4h.jm of06.4h.aa of05.2h.jm uf02.4h.aa
## 23312 23526 23916 24057 24064 24119
## of09.4h.ab op02.4h.jm of05.4h.aa of07.2h.ab of12.4h.jm uf05.4h.ab
## 26054 26422 26816 27521 27893 27964
## of11.4h.ab up02.2h.ab op01.4h.jm up01.2h.ab of12.2h.aa uf12.2h.aa
## 28939 29254 29271 29317 30192 30442
## op04.4h.jm of05.4h.ab uf06.2h.aa uf08.4h.ab of02.4h.jm of12.2h.jm
## 30820 30826 31470 33236 34169 34801
## of02.2h.jm up02.4h.jm of10.4h.jm up05.4h.jm op02.2h.aa of01.4h.jm
## 36046 37192 37617 37713 38675 38721
## uf01.4h.ab up05.2h.jm of08.2h.jm of11.4h.jm up04.2h.ab of01.2h.jm
## 39095 39258 40386 41980 43128 43367
## uf11.4h.ab op06.2h.aa of07.2h.aa of10.4h.aa up01.2h.aa up06.2h.aa
## 44531 44663 44820 45242 45419 45429
## op02.4h.aa up03.4h.aa of06.4h.jm of08.4h.ab of07.4h.jm of09.2h.jm
## 45633 45687 45882 47273 47494 48523
## op04.2h.ab of11.2h.aa of10.2h.aa of02.4h.ab of12.4h.ab op05.4h.jm
## 48552 48748 49228 55685 56417 58353
## of09.4h.jm op06.4h.aa of08.4h.jm op01.2h.jm uf07.4h.ab of07.4h.ab
## 62013 64385 64649 69957 76818 78604
## of01.4h.aa op03.4h.jm op03.2h.aa op01.4h.aa op05.4h.aa op06.4h.ab
## 82810 87558 87843 90912 100535 104632
## op02.2h.ab of06.4h.ab of01.4h.ab op02.4h.ab
## 113890 144827 173563 298953
pb.nc.1000 <- rrarefy(pb.nc, 1000)
pb.nc.1000 <- pb.nc.1000[, -which(colSums(pb.nc.1000) == 0)]
pb.1000 <- rrarefy(pb.c, 1000)
pb.1000 <- pb.1000[, -which(colSums(pb.1000) == 0)]
Make same NMDS objects to compare.
can.nc <- vegdist(pb.nc.1000, "canberra")
bc.nc <- vegdist(pb.nc.1000)
can.c <- vegdist(pb.1000, "canberra")
bc.c <- vegdist(pb.1000)
nmds.can.nc <- nmds(can.nc)
## initial value 39.616151
## iter 5 value 30.761031
## iter 10 value 30.282891
## iter 10 value 30.278175
## iter 10 value 30.249207
## final value 30.249207
## converged
nmds.bc.nc <- nmds(bc.nc)
## initial value 33.884719
## iter 5 value 23.374712
## final value 21.536789
## converged
nmds.can.c <- nmds(can.c)
## initial value 42.609751
## iter 5 value 31.960807
## iter 10 value 31.226954
## final value 31.151002
## converged
nmds.bc.c <- nmds(bc.c)
## initial value 26.232957
## iter 5 value 20.403711
## final value 19.973245
## converged
Still see really strong pattern, but dissimilarities are pretty different. So this says probably redo all analyses to see how big of a difference.
levels(pb.map$location)
## [1] "control" "east" "occ" "out" "unocc" "west"
plot(nmds.can.nc$points, col = pb.map.nc$location)
plot(nmds.bc.nc$points, col = pb.map.nc$location)
plot(nmds.can.c$points, col = pb.map.nc$location)
plot(nmds.bc.c$points, col = pb.map.nc$location)
plot(can.nc, can.c)
plot(bc.nc, bc.c)
Still a very big statistical difference.
adonis(can.nc ~ pb.map.nc$location)
##
## Call:
## adonis(formula = can.nc ~ pb.map.nc$location)
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## pb.map.nc$location 1 1.1 1.053 2.23 0.011 0.001 ***
## Residuals 206 97.0 0.471 0.989
## Total 207 98.1 1.000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis(bc.nc ~ pb.map.nc$location)
##
## Call:
## adonis(formula = bc.nc ~ pb.map.nc$location)
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## pb.map.nc$location 1 5.5 5.46 18.6 0.083 0.001 ***
## Residuals 206 60.5 0.29 0.917
## Total 207 66.0 1.000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(nmds.bc$points, col = pb.map$location, type = "n")
points(nmds.bc$points[occ, ], pch = 16, col = pb.map$bg[occ])
points(nmds.bc$points[unocc, ], pch = 16, col = pb.map$bg[unocc])
points(nmds.bc$points[control, ], pch = 21, bg = pb.map$bg[control])
plot(nmds.bc$points, type = "n", xlim = c(-0.5, -0.3), ylim = c(-0.2, 0.1))
points(nmds.bc$points[occ, ], pch = 16, col = pb.map$bg[occ])
points(nmds.bc$points[unocc, ], pch = 16, col = pb.map$bg[unocc])
points(nmds.bc$points[control, ], pch = 21, bg = pb.map$bg[control])