unable to run power_analysis

[melissaszy]

Hello,

I am trying to run power_analysis on my phyloseq object:

power_analysis(physeq)

I received these messages:
INFO [2023-11-20 14:53:09] Using median depth from reference data (6472).
INFO [2023-11-20 14:53:09] Estimating corncob model parameters for 941 taxa...
Error in as(x, "matrix")[i, j, drop = FALSE] : subscript out of bounds

May I ask what is wrong?

thanks so much!

[cdiener]

Hmm, looks like an error in corncob, but hard to diagnose without the data at hand or a more detailed traceback. I would start by checking that your OTU table does not include NAs and contains counts (integers).

[melissaszy]

Hello @cdiener,
Thank you for your suggestion.
I extracted the otu table from my phyloseq object and checked for NAs and to ensure that the data is in integers. Both conditions are fulfilled. How should I check if there is an error with corncob? Sorry, I am unable to upload my data here.

I am using mbtools version 0.49.0 and corncob version 0.4.1

x<-as.data.frame(physeq@otu_table)

sum(x$LF006e)
[1] 6472
sum(is.na(x))
[1] 0

Thank you for your time.

[cdiener]

Hmm okay, maybe something with estimating the parameters. What is the output after running traceback() right after the error occurs?

[melissaszy]

Here is what I get:

power_analysis(physeq)
INFO [2023-11-30 21:34:46] Using median depth from reference data (6472).
INFO [2023-11-30 21:34:46] Estimating corncob model parameters for 941 taxa...
Error in as(x, "matrix")[i, j, drop = FALSE] : subscript out of bounds
traceback()
9: otu_table(clean)[, taxon]
8: otu_table(clean)[, taxon]
7: FUN(X[[i]], ...)
6: apfun(taxa_names(clean), function(taxon) {
r <- tryCatch(corncob::bbdml(formula = reformulate("1", taxon),
phi.formula = ~1, data = clean), error = function(e) list(mu.resp = NA,
phi.resp = NA))
cn <- otu_table(clean)[, taxon]
data.table(mu = r$mu.resp[1], phi = r$phi.resp[1], mean_reads = mean(cn),
min_reads = min(cn), max_reads = max(cn), prevalence = mean(cn >
0), taxon = taxon)
})
5: apfun(taxa_names(clean), function(taxon) {
r <- tryCatch(corncob::bbdml(formula = reformulate("1", taxon),
phi.formula = ~1, data = clean), error = function(e) list(mu.resp = NA,
phi.resp = NA))
cn <- otu_table(clean)[, taxon]
data.table(mu = r$mu.resp[1], phi = r$phi.resp[1], mean_reads = mean(cn),
min_reads = min(cn), max_reads = max(cn), prevalence = mean(cn >
0), taxon = taxon)
})
4: rbindlist(.)
3: apfun(taxa_names(clean), function(taxon) {
r <- tryCatch(corncob::bbdml(formula = reformulate("1", taxon),
phi.formula = ~1, data = clean), error = function(e) list(mu.resp = NA,
phi.resp = NA))
cn <- otu_table(clean)[, taxon]
data.table(mu = r$mu.resp[1], phi = r$phi.resp[1], mean_reads = mean(cn),
min_reads = min(cn), max_reads = max(cn), prevalence = mean(cn >
0), taxon = taxon)
}) %>% rbindlist()
2: get_corncob_pars(ps, config$threads)
1: power_analysis(physeq)

Thank you so much for your help!

[cdiener]

Thanks that helps. I think I have a suspicion. Does taxa_are_rows(physeq) return TRUE in your case? If that is the case, power_ananlysis(orient_taxa(physeq, as="columns")) might fix this, and I can also add a fix to avoid this in the future.

[melissaszy]

Yes, taxa_are_rows(physeq) does return TRUE. Are you referring to orient_taxa from the package speedyseq?

I downloaded speedyseq package and got the following:

library(speedyseq)
Attaching package: ‘speedyseq’
The following objects are masked from ‘package:phyloseq’:
filter_taxa, plot_bar, plot_heatmap, plot_tree, psmelt, tax_glom, tip_glom, transform_sample_counts

power_analysis(orient_taxa(physeq, as="columns"))
INFO [2023-12-01 17:52:45] Using median depth from reference data (6472).
INFO [2023-12-01 17:52:45] Estimating corncob model parameters for 941 taxa...
Error in validObject(.Object) : invalid class “otu_table” object:
OTU abundance data must have non-zero dimensions.

I looked at my OTU table and it is non-zero:

physeq@otu_table
OTU Table: [ 941 taxa and 53 samples ]:
LF004e LF005e LF006e LF008e LF009e LF011e LF060e LF061e LF062e LF063e LF113e LF114e LF115e LF116e
1 4ab6c3b1b8eb737… 26 0 0 0 0 0 0 0 0 0 0 0 99 0
2 8f741a0a153cdae… 0 0 0 0 0 0 0 0 0 0 0 0 12 0

I ran traceback() right after getting the validObject(.Object) error and got the following:

traceback()
19: stop(msg, ": ", errors, domain = NA)
18: validObject(.Object)
17: .nextMethod(.Object = .Object, ... = ...)
16: callNextMethod()
15: .nextMethod(.Object = .Object, ... = ...)
14: callNextMethod()
13: initialize(value, ...)
12: initialize(value, ...)
11: new("otu_table", object, taxa_are_rows = taxa_are_rows)
10: .local(object, taxa_are_rows)
9: otu_table(newx, taxa_are_rows(x))
8: otu_table(newx, taxa_are_rows(x))
7: x[, taxa, drop = FALSE]
6: x[, taxa, drop = FALSE]
5: prune_taxa(taxa, otu_table(x))
4: prune_taxa(taxa, otu_table(x))
3: prune_taxa(pars$taxon, ps)
2: prune_taxa(pars$taxon, ps)
1: power_analysis(orient_taxa(physeq, as = "columns"))

Do you have any suggestions on how to proceed?

Thank you once again.

[cdiener]

Uh, bummer. Looks like none of the taxa have usable parameters. Will be hard to diagnose without the data set, but you could start by checking if mbtools:::get_corncob_pars(phyloseq) gives you sensible output (most mu and phi are not NA and phi is not much larger than mu).