So, I've run dbRDA in qiime, but I need to run it in R, in part for graphics and in part because I have two treatments which could be interacting.
I ran dbRDA on my own but got different results in R from the same analysis in Qiime. The total inertia is the same, but the constrained and unconstrained axes inertias are quite different (in some cases an order of magnitude different).
I am still a bit new to R, so I assumed I loaded my data incorrectly or something. I have been working with the R scripts provided by Qiime (I'm using Qiime1.6). I am not totally confident I have done everything correctly, but I am still getting different results from Qiime and from R for the same analysis. I'm wondering if anyone can tell me why.
These are my Qiime results from dbRDA:
Inertia Proportion Rank
Total 0.3910 1.0000
Constrained 0.1180 0.3017 1
Unconstrained 0.2730 0.6983 8
Inertia is squared Unknown distance
Eigenvalues for constrained axes:
CAP1
0.1180
Eigenvalues for unconstrained axes:
MDS1 MDS2 MDS3 MDS4 MDS5 MDS6 MDS7 MDS8
0.148233 0.052007 0.029738 0.019791 0.010449 0.007996 0.003089 0.001697
and these are the R results:
Inertia Proportion Rank
Total 0.39096 1.00000
Constrained 0.01210 0.03095 1
Unconstrained 0.37886 0.96905 8
Inertia is squared Unknown distance
Eigenvalues for constrained axes:
CAP1
0.0121
Eigenvalues for unconstrained axes:
MDS1 MDS2 MDS3 MDS4 MDS5 MDS6 MDS7 MDS8
0.241903 0.053415 0.032399 0.022755 0.016641 0.006657 0.003601 0.001485
I loaded all the functions in loaddata.r included in the qiime support_files into R.
I wasn't sure, though, how to use the dbrda.r script from support_files... so once I loaded the files I just ran dbrda using commands towards the bottom of the dbrda.r script.
Here is what I did:
#put functions loaddata.r into workspace
distmat<-load.qiime.distance.matrix(filepath="Desktop/Niwot/Analyses/R/Bacteria/DB/weighted_unifrac_rarefaction_1727_29.txt")
map<-load.qiime.mapping.file(filepath="Desktop/Niwot/Analyses/R/Bacteria/DB/NiwotB_D_map.txt")
otus<-load.qiime.otu.table(filepath="Desktop/Niwot/Analyses/R/Bacteria/DB/rarefaction_1727_29.txt",include.lineages=FALSE)
qiime.data<-remove.nonoverlapping.samples(map=map,otus=otus,distmat=distmat)
N=as.factor(map$Fertilizer)
factors.frame<-data.frame(N)
capscale.results <- capscale(as.dist(qiime.data$distmat) ~ N, factors.frame)
My Qiime command was: compare_categories.py --method adonis -i beta/distance_matrix/weighted_unifrac_rarefaction_1727_29.txt -m NiwotB_D_map.txt -c Fertilizer -o stats -n 999
