# How to make pcoa biplot in R using q2-deicode ordination
**URL:** https://forum.qiime2.org/t/how-to-make-pcoa-biplot-in-r-using-q2-deicode-ordination/8377
**Category:** Other Bioinformatics Tools
**Created:** [February 13, 2019, 11:38am UTC](https://forum.qiime2.org/t/how-to-make-pcoa-biplot-in-r-using-q2-deicode-ordination/8377 "2019-02-13T11:38:22Z")
**Posts on this page:** 19
**Page:** 1
### Author:  [@rahel\_park](https://forum.qiime2.org/u/rahel_park)
#### Post date: [February 13, 2019, 11:38am UTC](https://forum.qiime2.org/t/how-to-make-pcoa-biplot-in-r-using-q2-deicode-ordination/8377/1 "2019-02-13T11:38:22Z")
> [@Robust Aitchison PCA Beta Diversity with DEICODE](https://forum.qiime2.org/t/robust-aitchison-pca-beta-diversity-with-deicode/8333/1):
>
> Aitchison
Hi,
Thank you for a nice plugin and tutorial.
I was wondering if the visualization can be done in R as well? The emperor output is hard to export for modifications (the svg gets weirdly cut into pieces when imported to illustrator and the dots get merged together..). I imagine I can figure out how to make the PCoA ordination from the Aitchison distances, but not sure how to get the arrows for features on the plots. Have you done this and if yes, then could it be possible to have some guidance?
Thank you in advance.
---
### Author:  [@cmartino](https://forum.qiime2.org/u/cmartino)
#### Post date: [February 13, 2019, 6:22pm UTC](https://forum.qiime2.org/t/how-to-make-pcoa-biplot-in-r-using-q2-deicode-ordination/8377/2 "2019-02-13T18:22:40Z")
Hello @rahel_park ,
There is a way to move the ordination.qza output into R for biplot visualization. I would do this using the [qiime2R tutorial](https://forum.qiime2.org/t/tutorial-integrating-qiime2-and-r-for-data-visualization-and-analysis-using-qiime2r/4121) provided by @jbisanz. Specifically the command
```R
pco <- read_qza("unweighted_unifrac_pcoa_results.qza")
```
could be read in as
```R
RPCA <- read_qza("ordination.qza")
```
if you are following the tutorial above.
---
### Author:  [@rahel\_park](https://forum.qiime2.org/u/rahel_park)
#### Post date: [February 14, 2019, 9:15am UTC](https://forum.qiime2.org/t/how-to-make-pcoa-biplot-in-r-using-q2-deicode-ordination/8377/3 "2019-02-14T09:15:24Z")
Thank you for guiding me to the qiime2R tutorial.
It doesn't though look too straightforward to get to the PCoA plot including the feature arrows.
If anyone has already figured out a R script to get from the imported ordination qza file to a plot with the arrows, then could you please post your script here š
Thank you!
---
### Author:  [@cmartino](https://forum.qiime2.org/u/cmartino)
#### Post date: [March 5, 2019, 6:03pm UTC](https://forum.qiime2.org/t/how-to-make-pcoa-biplot-in-r-using-q2-deicode-ordination/8377/4 "2019-03-05T18:03:35Z")
Hello @rahel_park,
Looks as though @jbisanz solved this problem in the thread linked below. I hope this helps.
> [@Why PCoA and not PCA with DEICODE output?](https://forum.qiime2.org/t/why-pcoa-and-not-pca-with-deicode-output/8620):
>
> Can anyone explain to me why emperor biplots with DEICODE output are PCoA and not PCA? I'm second-guessing myself on this, and would like to get someone else's explanation to help clarify.
---
### Author:  [@jbisanz](https://forum.qiime2.org/u/jbisanz)
#### Post date: [March 5, 2019, 7:24pm UTC](https://forum.qiime2.org/t/how-to-make-pcoa-biplot-in-r-using-q2-deicode-ordination/8377/5 "2019-03-05T19:24:12Z")
Here is a minimal example but @cmartino can you please advise on the most appropriate way of selecting the taxa and for scaling their lengths? I am not 100% what I did is valid.
```auto
library(tidyverse)
library(qiime2R)
ord<-read_qza("~/Downloads/ordination.qza")
meta<-read_tsv(file = "~/Downloads/metadata.txt") %>%
rename(SampleID=`#SampleID`) %>%
filter(SampleID!="#q2:types")
tax<-read_qza("~/Downloads/taxonomy.qza")$data %>%
rename(FeatureID=Feature.ID)
#create the base plot with only the arrows
baseplot<-
ggplot() +
theme_bw() +
xlab(paste(round(100*ord$data$ProportionExplained[1],2),"%")) +
ylab(paste(round(100*ord$data$ProportionExplained[2],2),"%")) +
ggtitle("DEICODE biplot of Moving Pictures Data") +
geom_segment(data=ord$data$Species %>%
mutate(a=sqrt(PC1^2+PC2^2)) %>% # calculate the distance from the origin
top_n(8, a) %>% #keep 8 furthest away points
mutate(PC1=PC1*0.3, PC2=PC2*0.3) %>% # scale arrows linearly... is this ok?
left_join(tax),
aes(x=0, xend=PC1, y=0, yend=PC2, color=Taxon),
arrow = arrow(length = unit(0.3,"cm"))
)
# now overlay samples
baseplot +
geom_point(
data=ord$data$Vectors %>%
left_join(meta),
aes(x=PC1, y=PC2, fill=BodySite),
shape=21
)
```
---
### Author:  [@yoshiki](https://forum.qiime2.org/u/yoshiki)
#### Post date: [March 6, 2019, 12:09am UTC](https://forum.qiime2.org/t/how-to-make-pcoa-biplot-in-r-using-q2-deicode-ordination/8377/6 "2019-03-06T00:09:04Z")
> [@jbisanz](#):
>
> Here is a minimal example but @cmartino can you please advise on the most appropriate way of selecting the taxa and for scaling their lengths? I am not 100% what I did is valid.
The selection of the top arrows looks very similar, in the case of Emperor we use the magnitude based on all the dimensions which should be very close to the calculation you are making.
As for scaling, in Emperor the arrows and samples are scaled by the largest value in each matrix. The objective is to have dimensions that are normalized between the features and samples. Note that there are other strategies that might work for this problem.
From what I recall, in vegan the plot space is created with 2 pairs of x and y axes. One of the pairs represents the samples and the other pair represents the features. I believe doing this is not too hard in R but I am not familiar enough with ggplot.
---
### Author:  [@jbisanz](https://forum.qiime2.org/u/jbisanz)
#### Post date: [March 7, 2019, 6:01pm UTC](https://forum.qiime2.org/t/how-to-make-pcoa-biplot-in-r-using-q2-deicode-ordination/8377/7 "2019-03-07T18:01:31Z")
Using two independent sets of axes is a bit dicey in ggplot2, perhaps the solution is to just remove the axis labels all together.
```auto
plot + theme(axis.text=element_blank(), axis.ticks=element_blank())
```
---
### Author:  [@smreyes](https://forum.qiime2.org/u/smreyes)
#### Post date: [March 15, 2019, 9:08am UTC](https://forum.qiime2.org/t/how-to-make-pcoa-biplot-in-r-using-q2-deicode-ordination/8377/8 "2019-03-15T09:08:01Z")
Is it weird that the plot I made with this code only has 2 arrows? Instead of 8?
Also, random question. In a separate file, and using my own code, I plotted PC1 and PC2 for my samples (not including the taxonomy information). This plot doesn't show any separation of groups like it does using the code you provided. Do you possibly know why? I didn't do any transformations on the PC1, PC2 information, just plotted it with colors specific to my sample types. I see that in your code there is a transformation but it only seems to be in relation to the taxonomy data. Any thoughts would be helpful!
---
### Author:  [@jbisanz](https://forum.qiime2.org/u/jbisanz)
#### Post date: [March 15, 2019, 4:12pm UTC](https://forum.qiime2.org/t/how-to-make-pcoa-biplot-in-r-using-q2-deicode-ordination/8377/9 "2019-03-15T16:12:16Z")
Was this using your data or the moving pictures data? If your own data, have you checked your visualization against the q2 visualization? I would very carefully look at the values in data$Vectors and data$Species which define the X/Y plotting to see if they make sense versus what you are seeing in the biplot. If the Species data isn't scaled it might compress all of your samples into the middle of the plot making it hard to see the groupings.
---
### Author:  [@smreyes](https://forum.qiime2.org/u/smreyes)
#### Post date: [March 15, 2019, 5:25pm UTC](https://forum.qiime2.org/t/how-to-make-pcoa-biplot-in-r-using-q2-deicode-ordination/8377/10 "2019-03-15T17:25:01Z")
Hmmmm...this is a great point. I think you're right. How do you choose which scale to use? I see what you did in the previous example but not sure if there was a specific decision-making process you used.
---
### Author:  [@jbisanz](https://forum.qiime2.org/u/jbisanz)
#### Post date: [March 15, 2019, 8:16pm UTC](https://forum.qiime2.org/t/how-to-make-pcoa-biplot-in-r-using-q2-deicode-ordination/8377/11 "2019-03-15T20:16:30Z")
I picked arbitrarily, but I think as Yoshiki suggested above you could calculate a scaling factor for the X axis and Y axis by dividing the maximum PC1 and PC2 values for the Sample data by the Species Data.
---
### Author:  [@smreyes](https://forum.qiime2.org/u/smreyes)
#### Post date: [March 17, 2019, 10:32am UTC](https://forum.qiime2.org/t/how-to-make-pcoa-biplot-in-r-using-q2-deicode-ordination/8377/12 "2019-03-17T10:32:08Z")
@yoshiki, could you expand on scaling options?
Is there any other specific options that may be more useful than others? For example, I was thinking about sqrt(Eigvals$PC2/Eigvals$PC1) but don't know if that's appropriate. Can you advise?
---
### Author:  [@yoshiki](https://forum.qiime2.org/u/yoshiki)
#### Post date: [March 18, 2019, 4:01pm UTC](https://forum.qiime2.org/t/how-to-make-pcoa-biplot-in-r-using-q2-deicode-ordination/8377/13 "2019-03-18T16:01:23Z")
Hi @smreyes, I don't know what would make a scaling more appropriate than other. If you look at vegan's help for the [biplot.rda](http://cc.oulu.fi/~jarioksa/softhelp/vegan/html/biplot.rda.html) function, they list a few options:
> Scaling for species and site scores. Either species ( `2` ) or site ( `1` ) scores are scaled by eigenvalues, and the other set of scores is left unscaled, or with `3` both are scaled symmetrically by square root of eigenvalues. With negative scaling values in `rda` , species scores are divided by standard deviation of each species and multiplied with an equalizing constant. Unscaled raw scores stored in the result can be accessed with `scaling = 0` .
---
### Author:  [@11132](https://forum.qiime2.org/u/11132)
#### Post date: [November 6, 2020, 9:52am UTC](https://forum.qiime2.org/t/how-to-make-pcoa-biplot-in-r-using-q2-deicode-ordination/8377/14 "2020-11-06T09:52:44Z")
Hi Jbisanz,
Thanks for updating r code for DEICODE. Great help to R users.
I'm confused about the way the top 8 features were selected.
`mutate(a=sqrt(PC1^2+PC2^2)) %>% top_n(8, a)`
In the DEICODE tutorial, it wrote
[\*The important features with regard to sample clusters are not a single arrow but by the log ratio between features represented by arrows pointing in different directions. \*](https://forum.qiime2.org/t/robust-aitchison-pca-beta-diversity-with-deicode/8333)
Could you explain a little bit more on this?
Another problem I'm struggling with is the heatmap drawing as [the paper by Cameron](https://msystems.asm.org/content/4/1/e00016-19/figures-only). In Fig5A, rownames are sample names; colnames are features; values are feature loadings. But how can we match the sample name to feature names with loading?\
Best,
Yun
---
### Author:  [@jbisanz](https://forum.qiime2.org/u/jbisanz)
#### Post date: [November 12, 2020, 7:53am UTC](https://forum.qiime2.org/t/how-to-make-pcoa-biplot-in-r-using-q2-deicode-ordination/8377/15 "2020-11-12T07:53:24Z")
I think those are questions best answered by @cmartino!
---
### Author:  [@cmartino](https://forum.qiime2.org/u/cmartino)
#### Post date: [November 12, 2020, 5:15pm UTC](https://forum.qiime2.org/t/how-to-make-pcoa-biplot-in-r-using-q2-deicode-ordination/8377/16 "2020-11-12T17:15:07Z")
Hi Yun,
> [@11132](#):
>
> Iām confused about the way the top 8 features were selected.
> `mutate(a=sqrt(PC1^2+PC2^2)) %>% top_n(8, a)`
> In the DEICODE tutorial, it wrote
> [\*The important features with regard to sample clusters are not a single arrow but by the log ratio between features represented by arrows pointing in different directions. \*](https://forum.qiime2.org/t/robust-aitchison-pca-beta-diversity-with-deicode/8333)
> Could you explain a little bit more on this?
[Qurro](https://github.com/biocore/qurro) is a tool developed by @fedarko to help you do this form of analysis easily within QIIME2. [Here](https://nbviewer.jupyter.org/github/biocore/qurro/blob/master/example_notebooks/color_compositions/color_example.ipynb) is a tutorial that goes into depth on your question using a simple mock example dataset. [Here](https://nbviewer.jupyter.org/github/biocore/qurro/blob/master/example_notebooks/moving_pictures/moving_pictures.ipynb) is a more in-depth tutorial on choosing log-ratios of features with RPCA biplots & loadings on the moving pictures dataset. I am happy to answer any questions on those tutorials.
> [@11132](#):
>
> Another problem Iām struggling with is the heatmap drawing as [the paper by Cameron](https://msystems.asm.org/content/4/1/e00016-19/figures-only). In Fig5A, rownames are sample names; colnames are features; values are feature loadings. But how can we match the sample name to feature names with loading?\
We do not have any built-in functions for this. I would lean heavily to using log-ratios and Qurro instead. However, the procedure is not too complicated if you would like to try it:
1. Find the PC axis that `clusters` the samples by some metadata (in the paper one example is Healthy/Stressed sponges on PC1)
2. Extract the sample and feature loadings for that axis. (this is the PC axis loading used in the scatter and arrow for samples and features respectively).
3. Sort the CLR transformed `FeatureTable[Frequency]`(can use a pseudo count here) by those loadings from 2.
4. Plot the sorted table as a heatmap.
I hope this helps. Let me know if you have more questions.
Thanks,
Cameron
P.S. Thanks for tagging me @jbisanz.
---
### Author:  [@govindsah](https://forum.qiime2.org/u/govindsah)
#### Post date: [October 1, 2021, 11:11pm UTC](https://forum.qiime2.org/t/how-to-make-pcoa-biplot-in-r-using-q2-deicode-ordination/8377/17 "2021-10-01T23:11:34Z")
> [@jbisanz](#):
>
> `t`
Hi `jbisanz`
I am trying to make PCoA plot using `ordination.qza` following your script
When I try to import my metadata file `metadata_bacteriology.tsv` using the following script:
> ```
> meta<-read_tsv(file = "~/Govind/analysis_with_deblur/green_genes/lab_wise_analysis/bacteriology/metadata_bacteriology.tsv") %>%
> + rename(SampleID=`#SampleID`) %>%
> + filter(SampleID!="#q2:types")
>
> ```
it gives me the following error.
```
> Rows: 153 Columns: 10
```
> ```
> 0sāā Column specification āāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāā
> Delimiter: "\t"
> chr (10): sample-id, lab, location, lab-location, sampling-date, sampling-day, sequencing-lot, template-cleaned-before-PCR, total-input-reads, total-retained-reads
>
> ā¹ Use `spec()` to retrieve the full column specification for this data.
> ā¹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
> Error: Can't rename columns that don't exist.
>
> ```
>
> x Column `#SampleID` doesn't exist.
Can you please help?
Thank you,
Govind
here is my metadata file
[metadata\_bacteriology.tsv](https://forum.qiime2.org/uploads/short-url/72h98HbVtXULIJMCo4gUCmDxUts.tsv) (14.6 KB)
---
### Author:  [@jbisanz](https://forum.qiime2.org/u/jbisanz)
#### Post date: [October 5, 2021, 2:11am UTC](https://forum.qiime2.org/t/how-to-make-pcoa-biplot-in-r-using-q2-deicode-ordination/8377/18 "2021-10-05T02:11:51Z")
In your metadata file your sample id column is named `sample-id` not `#SampleID`. So you would want:
```auto
meta<-read_tsv(file = "~/Govind/analysis_with_deblur/green_genes/lab_wise_analysis/bacteriology/metadata_bacteriology.tsv") %>%
rename(SampleID=`sample-id`) %>%
filter(SampleID!="#q2:types")
```
You could also use:
```auto
meta<-read_q2metadata("~/Downloads/metadata_bacteriology.tsv")
```
---
### Author:  [@govindsah](https://forum.qiime2.org/u/govindsah)
#### Post date: [October 6, 2021, 4:21pm UTC](https://forum.qiime2.org/t/how-to-make-pcoa-biplot-in-r-using-q2-deicode-ordination/8377/19 "2021-10-06T16:21:17Z")
Thank you @jbisanz . It worked! š