Error in saved_plots[[label]][[plot_number]] when running MaAsLin2
versions :
R v4.2.3
Maaslin2 v1.14.1
Hi, I do have an error message that I don’t understand when I am using the following code. The problem is that I don’t have anything called “saved_plots[[label]][[plot_number]]” :
lm.s ← MaAsLin2.plus(ps=phyloseq(otu_table(ra.ps.s)/100,
sample_data(ra.ps.s)),
output=‘temp_directory’,
metadata=c(‘Case_status’,‘collection_method’,‘seqs_scaled’),
min_prevalence=0.05,
normalization=‘NONE’,
max_significance=0.05,
standardize=FALSE,
plot_heatmap=FALSE,
plot_scatter=FALSE,
save_models = FALSE,
reference = NULL)
…
2023-09-11 05:36:43 INFO::Creating boxplot for categorical data, Case_status vs k__Bacteria.p__Firmicutes.c__Bacilli.o__Lactobacillales.f__Streptococcaceae.g__Streptococcus.s__Streptococcus_mitis
2023-09-11 05:36:43 INFO::Creating boxplot for categorical data, Case_status vs k__Bacteria.p__Firmicutes.c__Clostridia.o__Clostridiales.f__Lachnospiraceae.g__Lachnoclostridium.s__Clostridium_hylemonae
2023-09-11 05:36:43 INFO::Creating boxplot for categorical data, Case_status vs k__Bacteria.p__Firmicutes.c__Bacilli.o__Lactobacillales.f__Streptococcaceae.g__Streptococcus.s__Streptococcus_vestibularis
Error in saved_plots[[label]][[plot_number]] :
attempt to select less than one element in integerOneIndex
For information, here is the code for the function that I’ve called MaAsLin2.plus :
MaAsLin2.plus ← function(ps, output, metadata,
min_abundance=0,
min_prevalence=0.1,
min_variance=0,
normalization=‘TSS’,
transform=‘LOG’,
analysis_method=‘LM’,
max_significance=0.25,
random_effects=NULL,
fixed_effects=NULL,
correction=‘BH’,
standardize=TRUE,
cores=1,
plot_heatmap=FALSE,
plot_scatter=FALSE,
heatmap_first_n=50,
reference=NULL,
max_pngs=10,
save_scatter = FALSE,
save_models = FALSE){
library(phyloseq)
library(Maaslin2)
ci ← function(coef, se){
lower.ci ← coef - 1.96se
upper.ci ← coef + 1.96se
return(c(lower.ci=lower.ci,upper.ci=upper.ci))
}
make sure samples are rows and features are columns
if (taxa_are_rows(ps)){
ps ← phyloseq(t(otu_table(ps)), sample_data(ps))
}
make sure only LOG is chosen for MaAsLin transformation
if (!(transform == “LOG”)){
stop(‘function only supports LOG transformation at this time’)
}
run MaAsLin2
input_data ← data.frame(otu_table(ps))
input_metadata ← data.frame(sample_data(ps)[,colnames(sample_data(ps)) %in% metadata])
fits ← Maaslin2(input_data, input_metadata, output, min_abundance, min_prevalence,
min_variance,normalization, transform, analysis_method, max_significance,
random_effects, fixed_effects, correction, standardize, cores,
plot_heatmap, plot_scatter, max_pngs, save_scatter, heatmap_first_n, save_models, reference)
put back original feature names
for (feat in seq_along(fits$results$feature)){
fits$results$feature[feat] ← taxa_names(ps)[make.names(taxa_names(ps)) ==
fits$results$feature[feat]]
}
filter for samples with data for variables included in model
for (var in seq_along(unique(fits$results$metadata))){
ps ← phyloseq(otu_table(ps),
sample_data(ps)[!is.na(sample_data(ps)[,metadata[var]]),])
}
get summary of results
res ← data.frame()
for (var in seq_along(unique(fits$results$metadata))){
# get variable name
var.name ← metadata[var]
# get N in each group
if (length(table(sample_data(ps)[,var.name])) == 2){
group.1.index ← sample_data(ps)[,var.name] ==
names(table(sample_data(ps)[,var.name]))[2]
group.1.index[is.na(group.1.index)] ← FALSE
group.2.index ← sample_data(ps)[,var.name] ==
names(table(sample_data(ps)[,var.name]))[1]
group.2.index[is.na(group.2.index)] ← FALSE
n1 ← colSums(otu_table(ps)[group.1.index,] > 0)
n2 ← colSums(otu_table(ps)[group.2.index,] > 0)
mean1 ← colMeans(otu_table(ps)[group.1.index,])
mean2 ← colMeans(otu_table(ps)[group.2.index,])
}else{
n1 ← rep(sum(table(sample_data(ps)[,var.name])), ntaxa(ps))
names(n1) ← taxa_names(ps)
n2 ← rep(NA, ntaxa(ps))
names(n2) ← taxa_names(ps)
mean1 ← colMeans(otu_table(ps))
mean2 ← rep(NA, ntaxa(ps))
names(mean2) ← taxa_names(ps)
}
# calculate fold change and confidence interval of fold change
if(length(table(sample_data(ps)[,var.name])) == 2){
FC ← 2^(fits$results$coef[fits$results$metadata == var.name])
FC.lower ← c()
FC.upper ← c()
for (coef in seq_along(fits$results$coef[fits$results$metadata == var.name])){
FC.lower ← c(FC.lower, 2^(ci(fits$results$coef[fits$results$metadata ==
var.name][coef],
fits$results$stderr[fits$results$metadata ==
var.name][coef])[‘lower.ci’]))
FC.upper ← c(FC.upper, 2^(ci(fits$results$coef[fits$results$metadata ==
var.name][coef],
fits$results$stderr[fits$results$metadata ==
var.name][coef])[‘upper.ci’]))
}
}else{
FC ← NA
FC.lower ← NA
FC.upper ← NA
}
# summarize results for variable
rvar ← data.frame(Variable=var.name,
Feature=fits$results$feature[fits$results$metadata == var.name],
N1=n1[fits$results$feature[fits$results$metadata == var.name]],
N2=n2[fits$results$feature[fits$results$metadata == var.name]],
Mean1=mean1[fits$results$feature[fits$results$metadata == var.name]],
Mean2=mean2[fits$results$feature[fits$results$metadata == var.name]],
Beta=fits$results$coef[fits$results$metadata == var.name],
SE=fits$results$stderr[fits$results$metadata == var.name],
P=fits$results$pval[fits$results$metadata == var.name],
FDR=p.adjust(fits$results$pval[fits$results$metadata == var.name],
method=correction),
FC=FC, FC_lower=FC.lower, FC_upper=FC.upper,
check.names=FALSE)
res ← rbind(res, rvar[order(rvar$P),])
# add untested features if they exist
if (nrow(rvar) != ntaxa(ps)){
res ← rbind(res,
data.frame(Variable=var.name,
Feature=taxa_names(ps)[!(taxa_names(ps) %in%
fits$results$feature[fits$results$metadata == var.name])],
N1=n1[taxa_names(ps)[!(taxa_names(ps) %in%
fits$results$feature[fits$results$metadata == var.name])]],
N2=n2[taxa_names(ps)[!(taxa_names(ps) %in%
fits$results$feature[fits$results$metadata == var.name])]],
Mean1=mean1[taxa_names(ps)[!(taxa_names(ps) %in%
fits$results$feature[fits$results$metadata == var.name])]],
Mean2=mean2[taxa_names(ps)[!(taxa_names(ps) %in%
fits$results$feature[fits$results$metadata == var.name])]],
Beta=NA, SE=NA, P=NA, FDR=NA, FC=NA, FC_lower=NA, FC_upper=NA,
check.names=FALSE)
)
}
}
return(list(result.summary=res, Maaslin2.output=fits))
}
The only thing I have found in the documentation about this specific in the viz.R file (located at https://github.com/biobakery/Maaslin2/tree/master/R/viz.R) but I don’t see how it is related to my problem… I am very confused and I don’t understand the package enought to solve this by myself. Maybe more experienced people could figure what is the problem here? Thanks a lot by advance…
Julien