R vegan NMDS vs mothur NMDS

Integrating mothur with other programs
1

Hi all,

I’m still not too experienced with mothur but thanks to my co-worker instructions, 454 SOP and this forum threads (which I’ve found really helpful), I’ve been able to process my 16S rRNA sequences, obtained by pyrosequencing.

Now that I’m performing the subsequent analyses, I’m stuck on the following issue: Running NMDS on the same data but with different programs (mothur “NMDS” command VS R vegan “metaMDS” function) ends up with different results.

The way I proceeded was as follows:

Mothur: I first generated a Bray-Curtis distance matrix (rarefying the data) to later calculate the NMDS. Then I plotted the obtained axes scores with R “plot” function.

mothur > dist.shared(shared=final.an.0.03.abund.shared, calc=thetayc-braycurtis, subsample=1955)
mothur > nmds(phylip=final.an.0.03.abund.braycurtis.0.03.lt.ave.dist, mindim=3, maxdim=3)

Number of dimensions: 3
Lowest stress : 0.252389
R-squared for configuration: 0.531698

Vegan: I imported the same Bray Curtis distance matrix to R through the function “import_mothur_dist” from phyloseq package. Then I applied vegan “metaMDS” (3D) using the following commands:

R> dist.bray <- import_mothur_dist("final.an.0.03.abund.braycurtis.0.03.lt.ave.dist")
R> dist.bray

R> m.bray <- metaMDS(dist.bray, k=3, trymax=100, wascores = T )
R> m.bray

Call:
metaMDS(comm = dist.bray, k = 3, trymax = 100, wascores = T) 

global Multidimensional Scaling using monoMDS

Data:     dist.bray 
Distance: user supplied 
Dimensions: 3 
Stress:     0.1089826 
Stress type 1, weak ties
Two convergent solutions found after 2 tries
Scaling: centring, PC rotation 
Species: scores missing

Not only is the stress clearly better when the vegan metaMDS is applied, but also plotting the results produced a much more interesting ordination with the metaMDS option.

Does someone have any idea of what is happening? Are these strong differences generated by the use of different NMDS algorithms? I appreciate any help!

Thanks,

Gema.

2

Interesting. Have you tried this with other datasets as well? Ours is modeled off of the version in ecodist (http://ftp.auckland.ac.nz/software/CRAN/doc/packages/ecodist.pdf).

3

Hi Patrick,

Thanks for your reply! I’ve just tried with another dataset from a previous pyrosequencing experiment, which yielded better quality reads.

The results of the two NMDS are different and vegan NMDS shows again a better result in terms of stress. Regarding the plots, the changes when applying vegan “metaMDS” are not so nice (probably because of the experimental design) but still different when comparing both programs figures.

I’ve also tried ecodist “nmds” function on the first dataset, with the default argument “stresscalc”, and it showed the same results than mothur nmds (as it was expected). When switching to “stresscalc=“kruskal”” results are quite similar to mothur’s :

dist.bray <- import_mothur_dist("final.an.0.03.abund.braycurtis.0.03.lt.ave.dist")
dist.bray

nmds.ecodist1 <- nmds(dist.bray, mindim = 3, maxdim = 3, nits = 10, iconf = 0, epsilon = 1e-12, 
     maxit = 500, trace = FALSE, stresscalc="default")

nmds.ecodist1

$stress
 [1] 0.2538556 0.2543998 0.2549385 0.2601659 0.2567985 0.2534177 0.2523264 0.2520055 0.2631301
[10] 0.2534806

$r2
 [1] 0.5325402 0.5217450 0.5130501 0.4569537 0.4830028 0.5381187 0.5384202 0.5424344 0.4069911
[10] 0.5203834

#----------------------

nmds.ecodist2 <- nmds(dist.bray, mindim = 3, maxdim = 3, nits = 10, iconf = 0, epsilon = 1e-12, 
                      maxit = 500, trace = FALSE, stresscalc="kruskal")

nmds.ecodist2

$stress
 [1] 0.2503991 0.2451236 0.2482832 0.2454704 0.2474642 0.2449086 0.2533447 0.2457429 0.2514687
[10] 0.2456907

$r2
 [1] 0.4633334 0.5259165 0.4908495 0.5280356 0.4923619 0.5318943 0.4238912 0.5298442 0.4502093
[10] 0.5307092

I’m using “mothur 1.34.1” and “vegan 2.3-4”. In this vegan version “metaMDS uses monoMDS as its NMDS engine” (http://www.inside-r.org/packages/cran/vegan/docs/initMDS), and monoMDS " implements Kruskal’s (1964) non-metric multidimensional scaling (NMDS) using monotone regression and primary (“weak”) treatment of ties." (http://www.inside-r.org/packages/cran/vegan/docs/monoMDS)

I’m a bit loss with all this algorithms stuff… Maybe I’m not correctly applying the “metaMDS” vegan function and this gives me a very unfortunate result (which looks as I would like it to be…)??

Thanks for your time and help!

Gema.

4

I’ve seen similar stress differences between ecodist (which I use to determine # of dimensions) and metaMDS, I don’t understand why there is this difference, I think I posted the question to R-sig-ecology but don’t remember getting any replies. metaMDS incorporates rotation of the solution to maximize variability along axis one, so I use that one as my final. This rotation is likely what’s causing your “more interesting ordination”.

Given your stress of .1 using metaMDS, I’d try 2 dimensions and see if the stress is still under .2. 2D is usually much more interpretable.

5

Hi kmitchell!

Thanks for replying to this topic. I’ll take your advices into consideration. I’m not sure if the differences in ordination plots are merely due to metaMDS rotation since the vegan approach clearly separates samples from different treatments whereas mothur NMDS don’t… Anyway, the cause of the differences in stress values is still unknown.