Hello -
I am trying to make a biplot with my data as suggested in the Costello analysis example. However, the output of the corr.axes is a bit confusing to me. Here is an example with the taxon definitions included from the sample.final.tx.1.consensus.taxonomy file.
OTU axis1 axis2 length OTU Size Taxonomy
2 -1 -0.028571 1.000408 2 10446 Sulfurovum
13 -0.941124 -0.27323 0.979984 13 93 Thiomicrospira
So, the length value means what exactly? I understand that the values under axis 1 and 2 represent the correlation coefficients. So when I made the biplot I made arrows that start at the origin(0,0) and go to these points (e.g. -1,-0.028). So, that way the viewer can see how correlated the taxa are with each axis in the NMDS plot. However, what does the length have to do with it? Should I make that the length of the arrow?
Also, how to choose which taxa to plot? Is the length an important indicator or should I just be considering the correlation coefficient?
In this case, I chose to only plot taxa with an overall length over 0.85 and a phylotype size over 90. But this may not be a good criterion since I don’t know what “length” means.
Also, why is does that label say OTU when this is a phylotype analysis? Is this just a legacy thing or is it true that when you create phylotypes they can still be referred to as OTUs (just based on a different definition - taxonomic rather than genetic distance)?
Thanks for the clarification.
Kristina
PS. Here is the handy R script I wrote to make the NMDS plot, add the arrows and the labels. It’s not sophisticated but it works.
#R script to create an nmds or pcoa plot in R using mothur output
#Also will add arrows and labels, as in a biplot, to show which OTUs are causing shifts in microbial community
#Make sure to add a column header “group” above your point labels in .thetayc.1.lt.nmds.axes file
Kristina Fontanez, July 2011
library(calibrate) library(shape) #-------------------------------------------------------------------- #Sends the output to a pdf file #-------------------------------------------------------------------- pdf("Thetayc genus NMDS.pdf")
#--------------------------------------------------------------------
#Read in the nmds or pcoa axes file
#--------------------------------------------------------------------
A <-read.table(“fries.final.tx.thetayc.1.lt.nmds.axes”, header=T)
#--------------------------------------------------------------------
Make the plot, add axis labels and set the x,y limits (change as needed)
#--------------------------------------------------------------------
plot(A$axis1, A$axis2, col=1, pch=20, type=“p”, xlab="axis 1 ", ylab="axis 2 ", main=“thetayc genus nmds”, xlim=c(-1.5,1.5), ylim=c(-1.5,1.5))
#--------------------------------------------------------------------
#Add labels to the points
#--------------------------------------------------------------------
textxy(A$axis1, A$axis2, A$group,cx=1)
#--------------------------------------------------------------------
#Read in the correlation datafile (Result of corr.axes function) and set up variables
#--------------------------------------------------------------------
points<-read.table(“thetayc.1.nmds.corr.txt”, header=T)
x1<-points$axis1
y1<-points$axis2
labels<-points$Taxonomy
#--------------------------------------------------------------------
#Add arrows and labels showing which OTUs are causing shifts in community
#--------------------------------------------------------------------
Arrows(0,0,x1,y1,code=2,arr.length=0.3,arr.width=0.3,arr.type=“triangle”,arr.adj=1)
textxy(x1,y1,labels,cx=0.8)
dev.off()