Skip to content

Instantly share code, notes, and snippets.

@darmitage
Created November 1, 2011 20:49
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save darmitage/1331859 to your computer and use it in GitHub Desktop.
Save darmitage/1331859 to your computer and use it in GitHub Desktop.
Nestedness Code
############################################################################
#CODE FOR ANALYZING NESTEDNESS OF A COMMUNITY MATRIX & PLOTTING THE RESULTS#
############################################################################
require(vegan)
require(ade4)
#Here, samp_85 is simply a community matrix of my 85% OTU cyano bins. The [5:9] just chooses
#specific rows for analysis (here, a depth gradient).
samp <- samp_85[c(5:9),]
sums <- subset(colSums(samp), colSums(samp) > 0)
samp <- samp[,names(sums)]
#Here's the code to calculate the NODF value for nestedness. It then compares
#this value to simulated values after shuffling. The "quasiswap" method is the most conservative
nf <- oecosimu(samp,nestednodf, method = "quasiswap", nsimul = 999, order = FALSE, weighted = FALSE)
nf
#This gives another, less informative measure of nestedness called "temperature
#its main benefit is that it makes a nice graph of the packed matrix, which is all I use it for
out <- nestedtemp(samp, order = FALSE)
#This does the same thing but for the amplicon libraries from the other cyano mat sections
samp <- samp_[c(4,10,1,11,12),]
sums <- subset(colSums(samp), colSums(samp) > 0)
samp <- samp[,names(sums)]
nf <- oecosimu(samp,nestednodf, method = "quasiswap", nsimul = 999, order = FALSE, weighted = FALSE)
nf
out1 <- nestedtemp(samp, order = FALSE)
## Finally, we can plot the packed matrix. Keep in mind that I have changed the row and
#column names for aesthetics, but you should simply start with plot(out, kind = "incid")
layout(1:2)
plot(out, kind="incid", col=rev(c("grey35", "white")), names = FALSE, adj = 1)
axis(2, at=seq(1,0,length.out = nrow(samp)), labels = rownames(samp), tick = FALSE, cex.axis = 0.75)
axis(3, at=seq(0,1,length.out = ncol(samp)), labels = seq(1,ncol(samp)), cex.axis = 0.75 , padj = 1, tick = FALSE)
title(main="Cyano OTUs (97%)\n New Mat", cex.main = 0.9, sub = "p = 0.9", cex.sub = 1, line = 2)
plot(out1, kind="incid", col=rev(c("grey35", "white")), names = FALSE, adj = 1)
axis(2, at=seq(1,0,length.out = nrow(samp)), labels = rownames(samp), tick = FALSE, cex.axis = 0.75)
axis(3, at=seq(0,1,length.out = ncol(samp)), labels = seq(1,ncol(samp)), cex.axis = 0.75 , padj = 1, tick = FALSE)
title(main="Cyano OTUs (97%)\n Old Mat", cex.main = 0.9, sub = "p = 0.30", cex.sub = 1, line = 2)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment