How to make pcoa biplot in R using q2-deicode ordination

Other Bioinformatics Tools
1

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.

6 Likes
2

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 provided by @jbisanz. Specifically the command

pco <- read_qza("unweighted_unifrac_pcoa_results.qza")

could be read in as

RPCA <- read_qza("ordination.qza")

if you are following the tutorial above.

4 Likes
3

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 :wink:
Thank you!

1 Like
4

Hello @rahel_park,

Looks as though @jbisanz solved this problem in the thread linked below. I hope this helps.

3 Likes
5

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.

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
      )

testplot2
testplot21599×806 90.6 KB

12 Likes
6

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.

6 Likes
7

Using two independent sets of axes is a bit dicey in ggplot2, perhaps the solution is to just remove the axis labels all together.

plot + theme(axis.text=element_blank(), axis.ticks=element_blank())
2 Likes
8

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!

9

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.

2 Likes
10

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.

11

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.

1 Like
12

@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?

13

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 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 .

1 Like
14

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. *
Could you explain a little bit more on this?

Another problem I’m struggling with is the heatmap drawing as the paper by Cameron. 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

15

I think those are questions best answered by @cmartino!

1 Like
16

Hi Yun,

Qurro is a tool developed by @fedarko to help you do this form of analysis easily within QIIME2. Here is a tutorial that goes into depth on your question using a simple mock example dataset. Here 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.

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.

3 Likes
17

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 (14.6 KB)

18

In your metadata file your sample id column is named sample-id not #SampleID. So you would want:

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:

meta<-read_q2metadata("~/Downloads/metadata_bacteriology.tsv")
1 Like
19

Thank you @jbisanz . It worked! :grinning: