KasperSkytte / ampvis2

Tools for visualising microbial community amplicon data
https://kasperskytte.github.io/ampvis2/
GNU General Public License v3.0
66 stars 23 forks source link

amp_rarecurve() bug with 0 reads observations #117

Closed brendbech closed 2 years ago

brendbech commented 2 years ago

Hello

I've run into a bug in your amp_rarecurve() function . With the following code i get the following error:

amp_rarecurve(data = amp, #A large phyloseq object converted to ampvis
              stepsize = 1000)

+ Error in seq.default(1, tot[i], by = stepsize) : 
+ wrong sign in 'by' argument 

The error occurred only for certain phyloseq objects, so i was pretty sure it had something to do with the data.

I did a little digging and found out that the error occurs when tot[i] has no reads, meaning it is 0, and therefore the sequence tries to go from 1 to 0 by a stepsize of 1000. Obviously this does not work. I fixed it by using trace(edit = TRUE), editing line 17 (It might not be the true line number). This fixed the bug for me, but i'm not sure if it will work for everybody.

Here is the fix:

out <- lapply(seq_len(nr), function(i) {
    n <- seq(1, if (tot[i] == 0) 1 else tot[i],  #ADDED THIS IF STATEMENT. The rest is as before.
    by = stepsize)
    if (n[length(n)] != tot[i]) {
      n <- c(n, tot[i])
    }
    drop(vegan::rarefy(abund[i, ], n))
  })
KasperSkytte commented 2 years ago

Hi there

Thanks for the proposed fix, I've implemented something similar. But empty samples shouldn't really be present in the first place, the problem is upstream in the data processing. But I've added a check for empty samples and a warning in amp_load, and made amp_rarecurve handle empty samples, though nothing will be visible in the plot anyways.