Description
Joey-
I've been using CCA to plot functional pathways. I'd like to plot the environmental vectors in a similar way that vegan does. Is there a way to do this in phyloseq?
I can make a pretty plot in phyloseq, which has everything I want, but it lacks the biplot arrows:
c=plot_ordination(TOP5,cca.deptreat,shape="TREATMENT",label="Definition",type="biplot",color="DEPTH",title="CCA of Top5 KEGG orthologs constrained by treatment*depth")
c+scale_color_gradient(limits=c(150,500),low="red",high="blue",expression(Depth[meters]))
In vegan I use the output of your ordinate function to make a plot with the arrows included:
cca.deptreat
Call: cca(formula = OTU ~ DEPTH * TREATMENT, data = data)
Inertia Proportion Rank
Total 0.4983 1.0000
Constrained 0.3211 0.6443 3
Unconstrained 0.1773 0.3557 3
Inertia is mean squared contingency coefficient
Eigenvalues for constrained axes:
CCA1 CCA2 CCA3
0.14973 0.11848 0.05285
Eigenvalues for unconstrained axes:
CA1 CA2 CA3
0.08029 0.05853 0.03846
plot(cca.deptreat,choices=c(1,2),display=c("wa","bp"),type="points",xlim=c(-4,1.5),scaling=2)
points(cca.deptreat,disp="sites",pch=21,col="red",bg="red",cex=1.3)
text(cca.deptreat,"sites",pos=3,axis.bp=TRUE)
I'd like to got those biplot arrows onto the phyloseq plot. Can this be done?
Thanks,
Kristina
Activity
joey711 commentedon Dec 4, 2013
Yes, see the
geom_segment
documentation for ggplot2:http://docs.ggplot2.org/0.9.3.1/geom_segment.html
There is an
arrow
argument.Tt would look something like:
You would have to get the arrow location data from the relevant ordination results to define
arrowdata
. I'll look into it more when I get a chance, but if you already know the coordinates, the above is how you would add it to your ggplot2 object.kfontanez commentedon Dec 4, 2013
Joey-
I might have figured it out but the results aren't as pretty as I'd like.
To get the biplot scores out of the ordination object:
To get the appropriate scaling multiplier used in the vegan plot in my previous post:
Next, multiply your biplot scores by the scaling multiplier:
Use the scaled biplot scores to draw the DEPTH arrow.
The arrowhead is missing and the colors look odd so I must have done something wrong here. I wasn't able to use your data= flag correctly so I could use some tips about how to use the arrowdata variable as input to the geom_segment function.
When I tried it, I got this:
I must be missing something basic. I tried changing the arrowdata to numeric with as.numeric but that failed as well.
What class does it have to be to work?
Thanks,
Kristina
joey711 commentedon Dec 5, 2013
That's a basic ggplot2 thing. Most data input to ggplot commands must be a
data.frame
. Try theas.data.frame
command. The columns of thedata.frame
will need names. It will probably help to read a little more about how ggplot2 expects data and arguments to work.kfontanez commentedon Dec 6, 2013
I think I've got it working.
First, change arrowdata to data.frame:
Then, set the x and y limits of the arrow for the first variable, Depth:
Now, plot:
The line is a bit wavy, which is not ideal but it works! You can do multiple arrows at a time by just adding multiple geom_segment calls like this:
joey711 commentedon Dec 9, 2013
That's cool!
The line waviness is probably an artifact of how you are viewing/saving the graphic images. ggplot2 creates pure vector graphics, so you should always be able to avoid pixelation/waviness if you need to / want to.
kfontanez commentedon Dec 10, 2013
Great, thanks for your help, as always!
Kristina
joey711 commentedon Dec 10, 2013
I'm going to keep this open until I post some more general code for adding this to a plot. I think it is a commonly-needed graphical annotation.
joey711 commentedon Dec 13, 2013
Fully Reproducible Example for CCA
How to make labels have same angle as arrows
joey711 commentedon Dec 13, 2013
Thanks again, @kfontanez for the feedback, and for getting this example started. I think it is potentially very helpful for others wanting to use this display.
chadsmith123 commentedon Aug 7, 2015
Can either you comment on how to interpret the arrows for the interaction terms in this analysis? Thanks!
montoyah commentedon Mar 16, 2017
Hi all,
I'm trying to reproduce this example to add arrows to my CAP plot, following this tutorial but have no idea where the arguments for "display" like "sites" or "species" in "vegan::scores()" function come from. Must they be variable names in the ps object? I even replaced one of my sam_data columns to "sites" but throws an error:
Error
in scores.default(ordcapPlot, display = "sites") : Can't find scores`In my case, I have taxa down to genus, no species. The Vegan help pages I found didn't clarify my doubts.
abundancesTop20NoNA <- abundancesTop20 %>% subset_samples( !is.na(phosphate..mM..alkaline) & !is.na(phosphate..mM..acidic) & !is.na(Acetate..mM.) & !is.na(sulfate..mM.) )
bray_no_na <- phyloseq::distance(physeq = abundancesTop20NoNA, method = "bray")
ordcap = ordinate(abundancesTop20NoNA, #abundances, method = "CAP", distance = bray_no_na, formula = ~ sites * phosphate..mM..alkaline)
Create plot:
ordcapPlot <- plot_ordination(abundancesTop20NoNA, ordination = ordcap, #type = "biplot", color="sites", axes = c(1, 2)) + aes(shape = sampleReceived) + geom_point(aes(color = type),alpha = 0.4, size = 4) + geom_point(colour = "grey90", size = 1.5) + scale_color_manual(values = c("#a65628", "red")) #, "#ffae19", "#4daf4a", "#1919ff", "darkorchid3", "magenta")
Here'e were the error appears:
arrowmat <- scores(ordcapPlot, display = "sites")
Thanks in advance for any help,
Oscar.
psObject.zip
jflater commentedon Jul 26, 2018
Try
arrowmat <- vegan::scores(ordcapPlot, display = "sites")
or try changing display to bpdisplay = "bp"