otu_table not working with vegan rarecurve

[ireneadm]

Hi,

I have recently installed the 1.42.0 version of phyloseq and I am trying to get rarefaction curves for my datasets. However, till I am using this version when I try for example:

data("GlobalPatterns") rare<-rarecurve(otu_table(GlobalPatterns), step=100, lwd=2, ylab="OTU", label=F)

I get this error: Error in as(x, "matrix")[i, j, drop = FALSE] : (subscript) logical subscript too long

and before I could make it work without problems. seems there is somenthing wrong with otu_table? has somenthing change in the function?

[ycl6]

Hi @ire1990 I didn't have the error message you have, though, one thing you should change in your code is to transpose the OTU table of GlobalPatterns.

From what I can tell from vegan's demo dataset such as BCI and dune, the rarecurve() expects the samples in rows and taxa in columns. Therefore you need to do t(otu_table(GlobalPatterns)) to make the orientation of the matrix right. You probably don't have to do this to another dataset previously if taxa_are_rows(your_ps_obj) returns FALSE.

library(phyloseq)
library(vegan)
#> Loading required package: permute
#> Loading required package: lattice
#> This is vegan 2.6-2

packageVersion("phyloseq")
#> [1] '1.38.0'
packageVersion("vegan")
#> [1] '2.6.2'

data(GlobalPatterns)

taxa_are_rows(GlobalPatterns)
#> [1] TRUE

mat <- t(otu_table(GlobalPatterns))
raremax <- min(rowSums(mat))

system.time(rarecurve(mat, step = 100, sample = raremax, col = "blue", label = FALSE))

#>    user  system elapsed
#> 507.948   0.337 508.340

^{Created on 2022-11-22 by the reprex package (v2.0.1)}

[ireneadm]

Hi @ycl6

Since I opened the issue I also re run the code transposing the otu table to have samples in rows but I get the same error.. I see you have a previous version of phyloseq and with that version I did not have problems running the rarecurve function with otu_table. I recently downloaded the 1.42 version and since then I can't run it anymore..

[ycl6]

Hi @ire1990

There isn't any change, code-wise, to the phyloseq package for a long time. The last 4 commits updates phyloseq from 1.38.0 to 1.41.1. Version 1.42.0 is probably just a version bump in Bioconductor. https://github.com/joey711/phyloseq/commits/master

I don't see any significant code change to rarecurve() in vegan too. Could it be R version? I was testing this in R-4.1.3.

[ireneadm]

Hi @ycl6

Yes could be the R version as I am using the R.4.2.2 that I recently updated.. however it seems really weird that if nothing has changed in the packages the code wont work anymore...

[ycl6]

Hi @ire1990

I saw vegan's developer has replied to you on stackoverflow. Indeed as Jari points out, the otu_table() function returns a otu_table class object, which rarecurve() wasn't able to convert it to a matrix class object using as.matrix() internally.

Before, all the operations in rarecurve() will run alright using a otu_table class object as input, but an update in Oct (https://github.com/vegandevs/vegan/commit/177aee1c18bce59e9b10ead1ab5e4a8adbda1ca0), a piece of code was added to count taxa in the input. In thise case, the operation can't be performed on the otu_table class object (it expects a matrix structure), resulting in the error message you saw.

    ## should be observed counts
    minobs <- min(x[x > 0])
    if (minobs > 1)
        warning(gettextf("most observed count data have counts 1, but smallest count is %d", minobs))

You can use the as() function to coerce the otu_table class object to a matrix class. This prevent the warning message from appearing when changing the class explicitly using class(x) <-.

library(phyloseq)
library(vegan)
#> Loading required package: permute
#> Loading required package: lattice
#> This is vegan 2.6-4

R.version.string
#> [1] "R version 4.2.2 (2022-10-31)"
packageVersion("phyloseq")
#> [1] '1.42.0'
packageVersion("vegan")
#> [1] '2.6.4'

data(GlobalPatterns)

taxa_are_rows(GlobalPatterns)
#> [1] TRUE
mat <- t(otu_table(GlobalPatterns))
class(mat) <- "matrix"
#> Warning message:
#> In class(mat) <- "matrix" :
#>  Setting class(x) to "matrix" sets attribute to NULL; result will no longer be an S4 object
class(mat)
#> [1] "matrix" "array"

mat <- as(t(otu_table(GlobalPatterns)), "matrix")
class(mat)
#> [1] "matrix" "array"

raremax <- min(rowSums(mat))

system.time(rarecurve(mat, step = 100, sample = raremax, col = "blue", label = FALSE))

#>    user  system elapsed
#> 457.546   0.105 457.694

^{Created on 2022-11-23 with reprex v2.0.2}