How to plot UniFrac PCoA with 95% confidence Elipses in R

Question

[c.v.oflynn]

Hi Guys,

I'm having difficulty plotting a PCoA for UniFrac distances with elipses. I'm using phyloseq to compute an ordination object and then creating elipses with ordiellipse() from vegan package. I can do almost exactly what I want for correspondence analysis (CCA), as in example below, or princomp() or other methods to create an ordination object. But for PCoA with unifrac the ordination object I create is incompatible with ordiellipse(). Does anyone know of an alternative method of doing this, or what I'm doing wrong?

Thanks

My example; borrowed code from https://github.com/joey711/phyloseq/issues/323

working example CCA

library(phyloseq)
library(vegan)

data(GlobalPatterns)

GP = GlobalPatterns
wh0 = genefilter_sample(GP, filterfun_sample(function(x) x > 5), A = 0.5 * nsamples(GP))
GP1 = prune_taxa(wh0, GP)

ord = ordinate(GP1~SampleType, "CCA", "bray")
(ordplot <- plot_ordination(GP1, ord, "samples", color="SampleType"))

# get data from the plotted ordination. Not strictly necessary but can be useful
orddata <- ordplot$data

# add a fake variable to the experiment
sample_data(GP1)$fake <- 1:nsamples(GP1)

# get map out of phyloseq object 
map.df <- data.frame(sample_data(GP1))

# use the envfit function of vegan
(ord.ef <- envfit(ord~fake,permu=1000,data = map.df))

# make a blank plot
plot(x=orddata$CCA1, y=orddata$CCA2,type="n") 

# plot the points for each type
points(ord, disp="si", cex = 0.5, col=map.df$SampleType,pch=19)

# draw ellipses
ordiellipse(ord, map.df$SampleType,label=TRUE,cex=0.7,conf=0.95,col="brown") 

# show the fitted fake variable
plot(ord.ef)

not working PCoA with unifrac

library(phyloseq)
library(vegan)

data(GlobalPatterns)

GP = GlobalPatterns
wh0 = genefilter_sample(GP, filterfun_sample(function(x) x > 5), A = 0.5 * nsamples(GP))
GP1 = prune_taxa(wh0, GP)

ord = ordinate(GP1, "PCoA", "unifrac", weighted = TRUE)
(ordplot <- plot_ordination(GP1, ord, "samples", color="SampleType"))

# get data from the plotted ordination. Not strictly necessary but can be useful
orddata <- ordplot$data

# add a fake variable to the experiment
sample_data(GP1)$fake <- 1:nsamples(GP1)

# get map out of phyloseq object 
map.df <- data.frame(sample_data(GP1))

##throws error
##Error in scores.default(ord, display = display, choices = choices, ...) : 
##  Can't find scores
# use the envfit function of vegan
(ord.ef <- envfit(ord~fake,permu=1000,data = map.df))

# make a blank plot
plot(x=orddata$Axis.1, y=orddata$Axis.2,type="n")

# plot the points for each type
points(ord$vectors[,1], ord$vectors[,2], disp="si", cex = 0.5, col=map.df$SampleType,pch=19)

##throws error
##Error in scores.default(ord, display = display, ...) : Can't find scores
# draw ellipses
ordiellipse(ord, map.df$SampleType,label=TRUE,cex=0.7,conf=0.95,col="brown") 

# show the fitted fake variable
plot(ord.ef)