buckleylab / HTSSIP

An R package for analyzing high throughput sequence stable isotope probing data
GNU General Public License v2.0
3 stars 3 forks source link

Error in dimnames(x) <- dn : length of 'dimnames' [1] not equal to array extent #9

Closed blancaverag closed 4 years ago

blancaverag commented 4 years ago

I have been testing HTSSIP in a phyloseq object of ASVs, created as indicated in DADA2 pipeline:

ps <- phyloseq(otu_table(seqtab.nochim, taxa_are_rows=TRUE), 
               sample_data(map2), 
               tax_table(taxa))

### Use short names for ASV
dna <- Biostrings::DNAStringSet(taxa_names(ps))
names(dna) <- taxa_names(ps)
ps <- merge_phyloseq(ps, dna)
taxa_names(ps) <- paste0("ASV", seq(ntaxa(ps)))

When I apply HTSSIP on a susbset of it I get an error:

physeq_heavyHMH = prune_samples((sample_data(ps)$Fraction=='H')| (sample_data(ps)$Fraction=='MH') , ps)
physeq_heavyHMH = prune_samples(sample_sums(physeq_heavyHMH) > 0, physeq_heavyHMH)

physeq_heavyHMH_ABC = prune_samples((sample_data(physeq_heavyHMH)$Treatment=='A')| (sample_data(physeq_heavyHMH)$Treatment=='B')| (sample_data(physeq_heavyHMH)$Treatment=='C') , physeq_heavyHMH)
physeq_heavyHMH_ABC = prune_samples(sample_sums(physeq_heavyHMH_ABC) > 0, physeq_heavyHMH_ABC)

### Define controls
params=get_treatment_params(physeq_heavyHMH_ABC, 'Treatment')
params
params=dplyr::filter(params, Treatment!='A')
params
ex = "(Treatment=='A') | (Treatment=='${Treatment}')"
physeq_l_ABC_HMH= phyloseq_subset(physeq_heavyHMH_ABC, params, ex)

## Run HRSIP
df_l2fc_ABC_H_HMH = plyr::ldply(physeq_l_ABC_HMH,HRSIP, design = ~Treatment, padj_cutoff = 0.05, sparsity_threshold = c(0,0.1,0.15,0.2, 0.25, 0.3),parallel=FALSE, density_windows=data.frame(density_min=1.75, density_max=1.78), sparsity_apply='heavy')
Sparsity threshold: 0 
Density window: 1.75-1.78 
converting counts to integer mode
Sparsity threshold: 0.1 
Density window: 1.75-1.78 
converting counts to integer mode
Sparsity threshold: 0.15 
Density window: 1.75-1.78 
converting counts to integer mode
Sparsity threshold: 0.2 
Density window: 1.75-1.78 
converting counts to integer mode
Error in dimnames(x) <- dn : 
  length of 'dimnames' [1] not equal to array extent

Since manipulation of taxa_names() in phyloseq object may derive in this type of error (see post) I have tried if long names of ASV would work, but it does not solve the problem.

I have also tried ps <- phyloseq(otu_table(t(seqtab.nochim), taxa_are_rows=FALSE), sample_data(map2), tax_table(taxa)) without success.

This same code works with a phyloseq object build in qiime and based on OTUs. I can't think of any other difference among these two phyloseq objects and in fact I can apply phyloseq functions to the object based on ASVs.

Any ideas of what the problem could be here? Thank you

blancaverag commented 4 years ago

It seems that there were spaces somewhere (not sure where). Removing them before creating the phyloseq object let me proceed normally with htssip command.

taxa[taxa == " "] <- ""