vegandevs / vegan

R package for community ecologists: popular ordination methods, ecological null models & diversity analysis
https://vegandevs.github.io/vegan/
GNU General Public License v2.0
449 stars 96 forks source link

Colnames problem and pairwise analysis #405

Closed elenu closed 3 years ago

elenu commented 3 years ago

Hi! I am trying to perform a PERMANOVA with metagenomics data (16S). But I have faced an error message when comparing goup distances, and I wonder how to deal with it to obtain pair-wise data.

I initially succeed if I try it with the filtared transposed otu table and metadata:

otus_t <- t(otu_table)
condition1 <- metadata$Tube2[metadata$Time.Point == 'T2']
adonis(otus_t[condition1, ] ~ metadata$Param1, perm=1000)

However, when I try to perform the same with distance parameters, I obtain always an error message:

BC.dist=vegdist(otus_t, distance="bray")
adonis(BC.dist ~ Param1, metadata, permutations = 1000)
## Error in `colnames<-`(`*tmp*`, value = colnames(lhs)) : 
##  attempt to set 'colnames' on an object with less than two dimensions

No matter what type of distance I try (we would be especially interested in Bray-Curtis, Jaccard, or Aitchinson distances), because I always obtain the same error message. Moreover, I have the same colnames in the distance parameter and in the metadata rows, so I do not understand the error.

Finally, I wonder whether you could advise me with the best way to obtain the significant otus. Given I am only interested in two categories of one group, I was considering to try a for-loop of adonis function to mimic the same analysis applied to the samples of the group. However, I do not know which otu I should use as factor because I do not have any special one to use as reference (I have more than 1000 otus).

Thank you very much in advance!

jarioksa commented 3 years ago

I cannot reproduce this, but I always succeed with anything I try, both with the CRAN release (2.5-7) and github version (2.6-0).

What does BC.dist look like? Does it contain numbers?

BTW, vegdist() does not know distance= argument: it uses method= (similarly as standard dist() of R).

elenu commented 3 years ago

Hi! Yo are wright. I introduced the wrong argument for distances. However, it's still not working. The distance parameter I obtain is similar to a correlation table. Here you have one section:

            Tube1        Tube2          Tube3         Tube4   
Tube2 0.4022656                                                                                                                        
Tube3 0.6173259 0.5280177                                                                                                              
Tube4 0.7254570 0.7269035 0.5915319                                                                                                    
Tube5 0.3682281 0.3499487 0.5685831 0.5409836 

I just realized the names of the columns and rows are, like in a correlation table, moved in order to compare samples. Might this be the reason of the error?

Thank you in advance,

jarioksa commented 3 years ago

I cannot reproduce this. I cannot find any problems. Need more information, such as a reproducible example.

DrMaestre commented 3 years ago

Did you solve it? I am having the same problem here. Thanks!

gavinsimpson commented 3 years ago

@elenu Is there a reason you didn't subset the data in the second example (with the error) in the same way as the first example?

gavinsimpson commented 3 years ago

@DrMaestre Can you provide a reproducible example or the at very least show all the code you are using and the output from str() applied to the objects in the adonis() call?

Better, if you or @elenu can share the data and code you are using privately to me (ucfagls at gmail dot com) I promise not to use them for anything other than solving the problem and then I'll delete them afterward

DrMaestre commented 3 years ago

Hi, thanks for your response. Here an example:

data(enterotype) ent = prune_taxa(!(taxa_names(enterotype) %in% "-1"), enterotype) ent = prune_taxa(taxa_sums(ent) > 0.0, ent) ent <- transform_sample_counts(ent, function(x) x/sum(x)) ent <- subset_samples(ent, !is.na(Enterotype)) df = data.frame(sample_data(ent)) BC = phyloseq::distance(ent, "bray") adonis(BC ~ SeqTech + Enterotype, data=df)

Error in colnames<-(*tmp*, value = colnames(lhs)) : attempt to set 'colnames' on an object with less than two dimensions

-phyloseq 1.30.0 -r 3.6.3 -vegan ‘2.5.7’

Thanks a lot!

On Mon, Apr 12, 2021 at 4:40 PM Gavin Simpson @.***> wrote:

@DrMaestre https://github.com/DrMaestre Can you provide a reproducible example or the at very least show all the code you are using and the output from str() applied to the objects in the adonis() call?

Better, if you or @elenu https://github.com/elenu can share the data and code you are using privately to me (ucfagls at gmail dot com) I promise not to use them for anything other than solving the problem and then I'll delete them afterward

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/vegandevs/vegan/issues/405#issuecomment-818261382, or unsubscribe https://github.com/notifications/unsubscribe-auth/AOQTT2O6Y3Y23DDN5EJCZPLTINSF7ANCNFSM42PCFDVA .

-- Thanks,

JP.

"The voyage of discovery is not in seeking new landscapes but in having new eyes"

Juan Pedro Maestre Research Associate at The University of Texas at Austin. +1 512 751 96 81 Pronouns: he, him, his Curiosity + Science + Diversity + Empathy + Development + Community

gavinsimpson commented 3 years ago

Thanks @DrMaestre but without the output from str() the code alone isn't that helpful. Particular, provide the output of str(df) and str(BC) please.

DrMaestre commented 3 years ago

Sorry, here it is.

str(df) 'data.frame': 271 obs. of 9 variables: $ Enterotype : Factor w/ 3 levels "1","2","3": 3 2 3 3 2 1 2 2 3 3 ... $ Sample_ID : Factor w/ 271 levels "AM.F10.T2","DA.AD.1",..: 1 2 3 4 5 6 7 8 9 10 ... $ SeqTech : Factor w/ 3 levels "Illumina","Pyro454",..: 3 3 3 3 3 3 3 3 3 3 ... $ SampleID : Factor w/ 32 levels "AM.F10.T2","DA.AD.1",..: 1 2 3 4 5 6 7 8 9 10 ... $ Project : Factor w/ 6 levels "kurokawa07","kurokawa07 ",..: 6 3 3 3 3 3 3 3 3 5 ... $ Nationality : Factor w/ 6 levels "american","danish",..: 1 2 2 2 2 6 6 6 6 3 ... $ Gender : Factor w/ 2 levels "F","M": 1 1 2 1 2 1 2 1 1 2 ... $ Age : num NA 59 54 49 59 25 49 47 38 63 ... $ ClinicalStatus: Factor w/ 5 levels "CD","elderly",..: 4 3 3 4 4 1 3 5 3 3 ...

str (BC) 'dist' Named num [1:36585] 0.709 0.646 0.488 0.571 0.69 ...

  • attr(*, "Size")= int 271
  • attr(*, "Labels")= chr [1:271] "AM.F10.T2" "DA.AD.1" "DA.AD.2" "DA.AD.3" ...
  • attr(*, "Diag")= logi FALSE
  • attr(*, "Upper")= logi FALSE
  • attr(*, "method")= chr "bray"
  • attr(*, "call")= language vegdist(x = c(0, 0, 0.0000363544800518205, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, | truncated ...

Thanks

On Mon, Apr 12, 2021 at 5:04 PM Gavin Simpson @.***> wrote:

Thanks @DrMaestre https://github.com/DrMaestre but without the output from str() the code alone isn't that helpful. Particular, provide the output of str(df) and str(BC) please.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/vegandevs/vegan/issues/405#issuecomment-818272636, or unsubscribe https://github.com/notifications/unsubscribe-auth/AOQTT2PWT7RAURC2CFYSOZLTINU57ANCNFSM42PCFDVA .

-- Thanks,

JP.

"The voyage of discovery is not in seeking new landscapes but in having new eyes"

Juan Pedro Maestre Research Associate at The University of Texas at Austin. +1 512 751 96 81 Pronouns: he, him, his Curiosity + Science + Diversity + Empathy + Development + Community

jarioksa commented 3 years ago

@DrMaestre Please specify also the version of vegan you are using: vegan versions 2.6-0 uses different adonis code than release versions (2.5-x and earlier). The current plan is to deprecate the adonis of the CRAN releases (2.5-x) and switch to a code that is called adonis2 in 2.5-x. Personally, I'm not interested in investing a lot in solving this problem in a code that is on its way to recycling bin: if you are using CRAN version 2.5-x, please check with adonis2 function that is more similar to the future code (that is the plan, but the best laid schemes o' mice an' men...).

Some forensic analysis: we do set colnames of beta coefficients in old CRAN adonis, but we do not set any colnames in new adonis2 of CRAN or new adonis of vegan 2.6-0 in github. Speculation warning: if these beta coefficients are based on one-dimensional fitted model, the results can drop to a vector and then we get the error attempt to set 'colnames' on an object with less than two dimensions. Solutions (subject to speculative analysis): use adonis2 in vegan 2.5-x (CRAN), or use adonis in vegan 2.6-0 (github).

elenu commented 3 years ago

Sorry for the late response. Due to we are on a rush and at least I've succeeded with the analysis of the richness (alpha diversity), we are not going to use the distances measurement in adonis.

In order to clarify my two previous comments, the second set of data was the visualization of the codes in the first subset: BC.dist. I also tried to do a matrix of the distance (as.matrix(BC.dist)), but it did not work neither in adonis.

Something I've realized is that adonis works with my normalized counts' table (otus_t) instead of the distance measurement. Adonis2 did not work neither with the distances, but with counts.

Interestingly, in my last try, the error message said Error in qr.fitted(qrhs, G) : 'qr' and 'y' must have the same number of rows. Here comes the code:

df_t <- t(otu_counts4)
#dist.bc <- vegan::vegdist(df_t, "bray") ## Bray-Curtis
dist.bc <- as.matrix(vegan::vegdist(df_t, "bray")) ## Bray-Curtis
adonis2(dist.bc ~ condition1, data = metadata)
str(df_t)
 num [1:92, 1:1180] 6.46 4.68 4.69 4.12 5.46 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:92] "Tube1" "Tube2" "Tube3" "Tube4" ...
  ..$ : chr [1:1180] "Zotu1" "Zotu2" "Zotu3" "Zotu4" ...
str(dist.bc)
 num [1:92, 1:92] 0.00 0.00 2.62e+14 4.15e+14 9.25e+13 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:92] "Tube1" "Tube2" "Tube3" "Tube4" ...
  ..$ : chr [1:92] "Tube1" "Tube2" "Tube3" "Tube4" ...

The package version I've used is ‘2.5.7’.

jarioksa commented 3 years ago

@elenu: there is something really strange in your data. Your dissimilarities dist.bc have some huge numerical values (such as 26200000000000 that is third in your output). That should not be possible, as Bray-Curtis is bound to range 0..1. The numerical data in your first example seven days ago were in correct range. I have no clue how to generate Bray-Curtis dissimilarities up to magnitude 1014.

DrMaestre commented 3 years ago

-phyloseq 1.30.0 -r 3.6.3 -vegan ‘2.5.7’

Thanks!!

On Tue, Apr 13, 2021 at 2:40 AM Jari Oksanen @.***> wrote:

@DrMaestre https://github.com/DrMaestre Please specify also the version of vegan you are using: vegan versions 2.6-0 uses different adonis code than release versions (2.5-x and earlier). The current plan is to deprecate the adonis of the CRAN releases (2.5-x) and switch to a code that is called adonis2 in 2.5-x. Personally, I'm not interested in investing a lot in solving this problem in a code that is on its way to recycling bin: if you are using CRAN version 2.5-x, please check with adonis2 function that is more similar to the future code (that is the plan, but the best laid schemes o' mice an' men...).

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/vegandevs/vegan/issues/405#issuecomment-818518040, or unsubscribe https://github.com/notifications/unsubscribe-auth/AOQTT2IAYCNMNQZTNKUB5UDTIPYP3ANCNFSM42PCFDVA .

-- Thanks,

JP.

"The voyage of discovery is not in seeking new landscapes but in having new eyes"

Juan Pedro Maestre Research Associate at The University of Texas at Austin. +1 512 751 96 81 Pronouns: he, him, his Curiosity + Science + Diversity + Empathy + Development + Community

elenu commented 3 years ago

@jarioksa I agree, that's in part why we decided to focus in a more trustable result (alpha diversity + adonis analysis).

If I calculate the distance with the counts' table without normalizing, it gives gives different values:

str(dist.bc)
 'dist' Named num [1:4186] 0.403 0.617 0.725 0.368 0.405 ...
 - attr(*, "Size")= int 92
 - attr(*, "Labels")= chr [1:92] "Tube1" "Tube2" "Tube3" "Tube4" ...
 - attr(*, "Diag")= logi FALSE
 - attr(*, "Upper")= logi FALSE
 - attr(*, "method")= chr "bray"
 - attr(*, "call")= language vegan::vegdist(x = df_t, method = "bray")

My count's data str() before normalizing:

str(otu_counts_t)
num [1:92, 1:1180] 21207 15017 5983 3497 15870 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:92] "Tube1" "Tube2" "Tube3" "Tube4" ...
  ..$ : chr [1:1180] "Zotu1" "Zotu2" "Zotu3" "Zotu4" ...

Distances calculation of counts without normalizing, do not allow adonis calculation neither.

Thank you!

DrMaestre commented 3 years ago

@DrMaestre Please specify also the version of vegan you are using: vegan versions 2.6-0 uses different adonis code than release versions (2.5-x and earlier). The current plan is to deprecate the adonis of the CRAN releases (2.5-x) and switch to a code that is called adonis2 in 2.5-x. Personally, I'm not interested in investing a lot in solving this problem in a code that is on its way to recycling bin: if you are using CRAN version 2.5-x, please check with adonis2 function that is more similar to the future code (that is the plan, but the best laid schemes o' mice an' men...).

Some forensic analysis: we do set colnames of beta coefficients in old CRAN adonis, but we do not set any colnames in new adonis2 of CRAN or new adonis of vegan 2.6-0 in github. Speculation warning: if these beta coefficients are based on one-dimensional fitted model, the results can drop to a vector and then we get the error attempt to set 'colnames' on an object with less than two dimensions. Solutions (subject to speculative analysis): use adonis2 in vegan 2.5-x (CRAN), or use adonis in vegan 2.6-0 (github).

Adonis 2 worked. Thanks so much!!!

jarioksa commented 3 years ago

The CRAN version of adonis has not changed in one month: this must be something with your data.

I cannot reproduce this: I don't have your data, I haven't even see the commands you used. The following is wild speculation

The first error in colnames probably comes from trying to set names of regression coefficient (betas) and these may be univariate or completely missing (in that case the model probably failed as well).

The second error is more cryptic as we do not have x1 anywhere in the code, and we only try to set the left-hand-side of the formula to NULL.

If you write traceback() after getting the error, we may at least see the line of code that gives the error.

hgoldspiel commented 2 years ago

I'm updating an analysis from a few years ago (before I knew about RProjects and version control with renv) and having the same error with adonis. My past analysis used vegan v2.5-5 and RVAideMemoire v0.9-73. I've updated those packages since, now using v2.5-7 for vegan and v0.9-81-2 for RVAideMemoire. I have tried installing older versions but I am running into errors that I assume are dependency issues with other packages in R that have since been updated. Before I go down the road of breaking my computer by manipulating all these package versions, I want to see if there is another fix.

I can get around the error with adonis by using adonis2, as suggested above, but I receive the same error when trying to get pairwise comparisons with pairwise.perm.manova from the RVAideMemoire package, which depends on adonis. I tried to manually edit the pairwise.perm.manova function by changing all mentions of adonis to adonis2 but that doesn't resolve the error.

Does anyone have other suggestions for me?

hgoldspiel commented 2 years ago

Well, I tried using the pairwise.adonis function from the pairwiseAdonis package and that seems to work and produce roughly the same results as my last analysis. Maybe I'll just stick with that.

gavinsimpson commented 2 years ago

@hgoldspiel Is there any chance you can provide the data (privately by email) so I can reproduce this?

We're planning to release 2.6.x to CRAN in the next month and we're unlikely to be releasing new versions too often after than for a while as Jari and I at least are very busy with other things. If we can reproduce the problem we can see how to fix this, otherwise I doubt it will get fixed