How to plot UniFrac PCoA with 95% confidence Elipses in R
0
2
Entering edit mode
7.9 years ago
c.v.oflynn ▴ 100

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)
metagenomics ordiellipse vegan R phyloseq • 9.8k views
ADD COMMENT

Login before adding your answer.

Traffic: 1041 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6