Question: How to plot UniFrac PCoA with 95% confidence Elipses in R
2
c.v.oflynn90 wrote:

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 incompatable with ordiellipse(). Does anyone know of an alternative method of doing this, or what im 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)