vegandevs / vegan

R package for community ecologists: popular ordination methods, ecological null models & diversity analysis
https://vegandevs.github.io/vegan/
GNU General Public License v2.0
449 stars 96 forks source link

adonis() on Bray-Curtis distance for any ordination? #462

Closed marwa38 closed 2 years ago

marwa38 commented 2 years ago

Hello .. running adonis test on bray-curtis distance is enough to know if there is statistical difference using PCoA or Unifrac ordination? so in the below example the p value is 0.002 in both cases (PCoA & Unweighted Unifrac), right?

bray.pcoa <- ordinate(ps.3.rare, method = "PCoA", distance = "bray")
# Plot using plot_ordination function of PhyloSeq:
plot_ordination(ps.3.rare, bray.pcoa, color = "treatment", axes = c(1,2)) +
  geom_point(size = 2) +
  labs(title = "PCoA of Bray Curtis Distances", color = "treatment") +
  stat_ellipse()

image

meta.ps.3.rare <- as(sample_data(ps.3.rare), "data.frame")
adonis(distance(ps.3.rare, method="bray") ~ treatment, data = meta.ps.3.rare)
# Call:
#   adonis(formula = distance(ps.3.rare, method = "bray") ~ treatment,      data = meta.ps.3.rare) 
# 
# Permutation: free
# Number of permutations: 999
# 
# Terms added sequentially (first to last)
# 
# Df SumsOfSqs  MeanSqs F.Model   R2 Pr(>F)   
# treatment  1   0.04418 0.044179  1.8818 0.07  0.002 **
#   Residuals 25   0.58694 0.023478         0.93          
# Total     26   0.63112                  1.00          
# ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

unwt.unifrac <- plot_ordination(ps.3.rare, 
                                ordu.unwt.uni, color="treatment")
unwt.unifrac <- unwt.unifrac + ggtitle("Unweighted UniFrac") + geom_point(size = 2)
unwt.unifrac <- unwt.unifrac + theme_classic() + scale_color_brewer("treatment", palette = "Set2")
unwt.unifrac <- unwt.unifrac + stat_ellipse()  
print(unwt.unifrac)

image

your advice and/or comment is much appreciated. many thanks Marwa

jarioksa commented 2 years ago

I think you have to consult Unifrac documentation: I have no idea what is "unifrac ordination". I have several books on multivariate statistics, but I cannot find "unifrac ordination" in any of them. NB., the graphs are wrongly made: ordination needs equal aspect ratio: 0.1 distance should be the same length vertically and horizontally.

marwa38 commented 2 years ago

Many thanks for the info @jarioksa could you please advise on what is wrong with the code (I am using vegan) that produced the wrong ordination?

bray.pcoa <- ordinate(ps.3.rare, method = "PCoA", distance = "bray")
plot_ordination(ps.3.rare, bray.pcoa, color = "treatment", axes = c(1,2)) +
  geom_point(size = 2)
gavinsimpson commented 2 years ago

The problem is with ordinate() which is not a vegan function. I have repeatedly engaged with the lead developer of phyloseq about this but they don't change the code as of yet. It needs to have a + coord_equal() applied to the ggplot object to get the correct aspect ratio of the plot.

gavinsimpson commented 2 years ago

Apologies I meant this is an issue in plot_ordinate() which also isn't a vegan function.

marwa38 commented 2 years ago

Thanks very much, Gavin @gavinsimpson I used coord_equal() with plot_ordination and worked fine. I was wondering if you think rooting the tree will have any effect on the ordination and distance? Which vegan function I could use for plotting the ordination that doesn't have this problem? cheers

gavinsimpson commented 2 years ago

It depends - someone tweeted about this when I was asking and said they didn't recall noticing any qualitative differences between ordinations using weight unifrac or Bray Curtis. Of course there could be differences induced if your tree has issues. Which is why I don't understand the motivation for unifrac and weighted unifrac distance - it simply doesn't make any sense to me to move away from simply computing a distance/dissimilarity in the data as provided. If you're using the tree to get around problems in calling OTUs then you should fix the OTU calling part, not bolt on some fudge from the phylogeny. But then again I don't do too much with HTPS data.

Vegan has plot() methods plus lots of other helper functions for plotting ordinations; see the documentation for the package. But you can just use plot_ordinate() if you prefer the ggplot2 output. I have a ggvegan package on github to help interface vegan with ggplot2 but it's bare bones at the moment and only users that know what they are doing should really use it.