Question: How to plot UniFrac PCoA with 95% confidence Elipses in R
2
gravatar for c.v.oflynn
4.8 years ago by
c.v.oflynn90
United Kingdom
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) 

 

ADD COMMENTlink modified 4.6 years ago by Biostar ♦♦ 20 • written 4.8 years ago by c.v.oflynn90
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1721 users visited in the last hour