Performing the robust Aitchison transformation prior to interpretative multivariate analysis (db-RDA)

[magibc]

Dear @cameronmartino Cameron,

The RPCA tool inside Gemelli performs these two steps:

1. Compositional preprocessing using the centered log ratio transform on only the non-zero values of the data (no pseudo count)
2. Dimensionality reduction through Robust PCA on only the non-zero values of the data ( matrix completion).

I can also perform Step 1 using the decostand() function from the R package vegan with the method rclr.

However, for the second step of creating the PCA using only non-zero values (as Gemelli creates a PCA), I would like to know if any R package could perform this second step to graph a db-RDA. I was not capable to found one R package.

The function cmdscale allows me to use the robust Aitchison method, but I'm unsure if it performs Step 2 as Gemelli does. That's why I'm asking you.

Thank you again,

Magi

[cameronmartino]

Hi @magibc,

There is an RPCA function (step 2) in the works for R (will be in vegan) but it is not ready yet and won't match Gemelli, hopefully it will be soon. To my knowledge there is no working version in R anywhere else yet. For now I advise relying on the python version and I can report back here when the R version of step 2 (RPCA) is done.

Thanks!

Cameron

[magibc]

Thank you very much,

Ok, for the moment I will use your robust.aitchison method included in vegan package and I will find Robust PCA in another R Package. Yes, keep me informed when it is implemented in Gemelli. Thanks again.

Magí.

[cameronmartino]

Hi @magibc,

Please use the Python implementation in Gemelli. The robust.aitchison method in vegan does not currently match the output in this implementation. I am working on getting the vegan implementation fixed but it is not ready yet and I would not suggest using it.

[magibc]

Ok Cameron. Yes I will try to defend the reviewer comments received because I have tried the function fill.OptSpace() of filling R package but I get an error that I am not capable to solve it. Thanks again.

[cameronmartino]

Hi @magibc,

The R OptSpace function also does not replicate this package exactly. There have been alterations and bug fixes over the years in this implementation that are not reflected in the OptSpace R package. I am working on making sure that is reflected in the vegan implementation but that is not done yet. For now, stick to this Python package.

[cameronmartino]

Hi @magibc,

RPCA is now in vegan, see this PR: vegandevs/vegan#667

It will be released in the next vegan version update, but can be installed through the github now. Here is an example of the use of the Python and R versions:

wget https://github.com/biocore/gemelli/raw/refs/heads/master/ipynb/tutorials/qiime2-moving-pictures-tutorial/table.biom
# export to tsv for R
biom convert -i table.biom -o table.tsv --to-tsv
sed -i '' 1d table.tsv 

library(devtools)
install_github("antagomir/vegan", ref = "nafix")
library(vegan)
# load table
table <- read.table("table.tsv", header = T, comment = '', check = F, sep = '\t', row.names = 1)
table <- t(as.matrix(table))
rclr_table <- vegan::decostand(table, method="rclr")
# calculate the distance
distance = vegan::vegdist(rclr_table, method = "euclidean")
write.table(as.matrix(distance), file = "R-rpca-distances.tsv", row.names=rownames(table), col.names=rownames(table), sep="\t")
# save to test RPCA
svd.res <- svd(rclr_table)
write.table(svd.res$u, file = "R-rpca-sample-loadings.tsv", row.names=rownames(table), sep="\t")
write.table(svd.res$d, file = "R-rpca-eig-vals.tsv", row.names=FALSE, sep="\t")
write.table(svd.res$v, file = "R-rpca-feature-loadings.tsv", row.names=colnames(table), sep="\t")

import numpy as np
import pandas as pd
from biom import load_table
from gemelli.rpca import rpca
from gemelli import __version__
print('Version: %s' % __version__)

# load table
bt = load_table('table.biom')
# run RPCA (this runs RCLR internally)
py_ordination, py_distance = rpca(bt)
# extract RPCA results from skbio type
py_rpca_U = py_ordination.samples.values
py_rpca_V = py_ordination.features.values
py_rpca_S = py_ordination.eigvals.values
py_dist = py_distance.to_data_frame().values
# check the OptSpace results
R_rpca_U = pd.read_csv('R-rpca-sample-loadings.tsv', sep='\t', index_col=None).values[:, :3]
R_rpca_S = pd.read_csv('R-rpca-eig-vals.tsv', sep='\t', index_col=None).values[:3]
R_rpca_V = pd.read_csv('R-rpca-feature-loadings.tsv', sep='\t', index_col=None).values[:, :3]
R_rpca_dist = pd.read_csv('R-rpca-distances.tsv', sep='\t', index_col=None).values
# check
print(' Sample loadings match: %s \n' % str(abs((np.abs(R_rpca_U) - np.abs(py_rpca_U)).sum()) < 1e-6),
        'Feature loadings match: %s \n' % str(abs((np.abs(R_rpca_V) - np.abs(py_rpca_V)).sum()) < 1e-6),
        'Singular values match: %s \n' % str(abs((np.abs(R_rpca_S) - np.abs(py_rpca_S)).sum()) < 1e-6),
        'Distances match: %s \n' % str(abs((R_rpca_dist - R_rpca_dist).sum()) < 1e-6))

Version: 0.0.12
 Sample loadings match: True 
 Feature loadings match: True 
 Singular values match: True 
 Distances match: True