SIMPER with Phyloseq Object

[pauGuas]

Hi, I am trying to run SIMPER on a phyloseq object. The code has been running for 24 hours with no results yet. I only have 259 samples, so I'm concerned something isn't working. Here is the code I am using:

GP <- readRDS("ps_norm_seasons.rds")
# Transpose GP's OTU table, and use Basin as the the group structure
s <- vegan::simper(t(otu_table(GP)), sample_data(GP)$Basin)

Any help would be greatly appreciated!

[gavinsimpson]

How many OTUs? The method involves at least two loops over the columns (species) of comm.

[pauGuas]

@gavinsimpson I believe it's about 9000. I'm trying again on a better computer....

[jarioksa]

Which version of vegan do you use? Function simper was completely redesigned in vegan 2.6-0 (of github), and no longer has two nested loops over species within a loop of permutations – which in your case will be 999 × 9000 × 9000 ≈ 8.1 × 10¹⁰. However, with your data the greatest pressure may come from the physical memory of your computer. I think you may need at least 16GB memory, and if you have less, R will not work at its full capacity. However, even in optimal case, analysis of 9000 species will take time.

Off-topic for the technical issue: I am not sure that simper will be very useful in your case of 9000 species on 259 sampling units. The results will probably be hard to interpret with that number of species. Please consider also this text in the simper help (?simper):

The results of ‘simper’ can be very difficult to interpret and
they are often misunderstood even in publications. The method
gives the contribution of each species to overall dissimilarities,
but these are caused by variation in species abundances, and only
partly by differences among groups. Even if you make groups that
are copies of each other, the method will single out species with
high contribution, but these are not contributions to non-existing
between-group differences but to random noise variation in species
abundances. The most abundant species usually have highest
variances, and they have high contributions even when they do not
differ among groups. Permutation tests study the differences among
groups, and they can be used to find out the species for which the
differences among groups is an important component of their
contribution to dissimilarities. Analysis without ‘group’ argument
will find species contributions to the average overall
dissimilarity among sampling units. These non-grouped
contributions can be compared to grouped contributions to see how
much added value the grouping has for each species.
[cited from vegan 2.6-0 help.]

[pauGuas]

@jarioksa Thank you for your input. I was finally able to get SIMPER to run, but I have the warning message "you have empty rows; your data may be meaningless." I have tried replacing blanks with NA then deleting rows with NA to remove empty rows, but there are none in my OTU file or phyloseq objects. Any ideas how I can fix this? Unfortunately the only person in my lab who has run this is unavailable this week.

[jarioksa]

You could try:

otudat <- t(otu_table(GP))
otudat <- otudat[rowSums(otudat) > 0,]
## in case you also have empty otu
otudat <- otudat[, colSums(otudat) > 0]
simper(otudat, sample_data(GP)$Basin)

[pauGuas]

@jarioksa this worked! A question, I can only run this analysis on a few samples at a time on my computer. Does SIMPER need to be run on the entire dataset to get accurate results?

I also found out that I need to use a phyloseq object that has not been normalized, in case anyone else is looking into this in the future.