Error: In has.phylogeny(), biom must be a BIOM-class objec

[nick-youngblut]

Reprex

library(rbiom)

# Using a custom matrix and tree
mtx <- matrix(sample.int(12*20), ncol=20)
dimnames(mtx) <- list(LETTERS[1:12], letters[1:20])
tree <- ape::as.phylo(hclust(dist(mtx)))

dm <- unifrac(mtx, tree=tree)
dm

Error

Error: In has.phylogeny(), biom must be a BIOM-class object.
Traceback:

1. unifrac(mtx, tree = tree)
2. rbiom::beta.div(biom, "unifrac", weighted, tree)
3. validate_metrics(biom, method, mode = "bdiv")
4. do.call(paste0(mode, "_metrics"), list(biom = biom))
5. bdiv_metrics(biom = structure(c(123L, 63L, 209L, 179L, 159L, 
 . 232L, 99L, 205L, 74L, 104L, 208L, 161L, 3L, 80L, 139L, 214L, 
 . 203L, 169L, 12L, 96L, 225L, 46L, 39L, 76L, 220L, 207L, 10L, 212L, 
 . 109L, 106L, 165L, 81L, 83L, 158L, 184L, 86L, 93L, 235L, 134L, 
 . 150L, 71L, 239L, 62L, 28L, 198L, 167L, 23L, 113L, 58L, 5L, 128L, 
 . 15L, 25L, 147L, 57L, 172L, 133L, 42L, 35L, 199L, 97L, 77L, 143L, 
 . 230L, 121L, 127L, 69L, 51L, 20L, 38L, 122L, 101L, 114L, 141L, 
 . 36L, 27L, 157L, 49L, 138L, 6L, 217L, 84L, 70L, 4L, 200L, 52L, 
 . 89L, 181L, 65L, 18L, 130L, 8L, 228L, 17L, 146L, 61L, 131L, 31L, 
 . 118L, 152L, 140L, 176L, 40L, 45L, 221L, 189L, 19L, 87L, 237L, 
 . 66L, 82L, 129L, 68L, 91L, 22L, 185L, 142L, 44L, 180L, 136L, 233L, 
 . 219L, 168L, 126L, 125L, 119L, 191L, 210L, 116L, 144L, 154L, 148L, 
 . 196L, 33L, 21L, 224L, 234L, 223L, 202L, 153L, 226L, 190L, 100L, 
 . 43L, 14L, 1L, 156L, 215L, 137L, 111L, 90L, 177L, 41L, 92L, 238L, 
 . 112L, 193L, 187L, 240L, 135L, 47L, 2L, 201L, 115L, 213L, 37L, 
 . 7L, 26L, 54L, 149L, 151L, 183L, 182L, 103L, 72L, 195L, 155L, 
 . 236L, 48L, 164L, 59L, 9L, 188L, 218L, 55L, 75L, 163L, 30L, 231L, 
 . 60L, 171L, 78L, 132L, 211L, 194L, 162L, 98L, 56L, 173L, 192L, 
 . 16L, 216L, 110L, 53L, 102L, 175L, 50L, 11L, 206L, 124L, 222L, 
 . 13L, 34L, 120L, 229L, 79L, 117L, 178L, 186L, 73L, 166L, 108L, 
 . 227L, 197L, 67L, 94L, 29L, 174L, 24L, 170L, 145L, 85L, 160L, 
 . 204L, 64L, 88L, 105L, 95L, 107L, 32L), .Dim = c(12L, 20L), .Dimnames = list(
 .     c("A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", 
 .     "L"), c("a", "b", "c", "d", "e", "f", "g", "h", "i", "j", 
 .     "k", "l", "m", "n", "o", "p", "q", "r", "s", "t"))))
6. metrics(biom, "bdiv")
7. has.phylogeny(biom)

sessionInfo

R version 4.1.1 (2021-08-10)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: Ubuntu 18.04.5 LTS

Matrix products: default
BLAS/LAPACK: /ebio/abt3_projects/Georg_animal_feces/envs/phyloseq2a/lib/libopenblasp-r0.3.17.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] rbiom_1.0.2.9038  LeyLabRMisc_0.2.0 vegan_2.5-7       lattice_0.20-44  
 [5] permute_0.9-5     phyloseq_1.36.0   ape_5.5           tidytable_0.6.5  
 [9] data.table_1.14.0 ggplot2_3.3.5     tidyr_1.1.3       dplyr_1.0.7      

loaded via a namespace (and not attached):
 [1] Biobase_2.52.0         jsonlite_1.7.2         splines_4.1.1         
 [4] foreach_1.5.1          RcppParallel_5.1.4     stats4_4.1.1          
 [7] GenomeInfoDbData_1.2.6 slam_0.1-48            ggrepel_0.9.1         
[10] pillar_1.6.2           glue_1.4.2             uuid_0.1-4            
[13] digest_0.6.27          XVector_0.32.0         stringfish_0.15.2     
[16] colorspace_2.0-2       htmltools_0.5.2        Matrix_1.3-4          
[19] plyr_1.8.6             pkgconfig_2.0.3        zlibbioc_1.38.0       
[22] purrr_0.3.4            scales_1.1.1           RApiSerialize_0.1.0   
[25] tibble_3.1.4           mgcv_1.8-36            generics_0.1.0        
[28] IRanges_2.26.0         ellipsis_0.3.2         withr_2.4.2           
[31] repr_1.1.3             BiocGenerics_0.38.0    survival_3.2-13       
[34] magrittr_2.0.1         crayon_1.4.1           evaluate_0.14         
[37] fansi_0.4.2            doParallel_1.0.16      nlme_3.1-153          
[40] MASS_7.3-54            tools_4.1.1            lifecycle_1.0.0       
[43] stringr_1.4.0          Rhdf5lib_1.14.0        S4Vectors_0.30.0      
[46] munsell_0.5.0          cluster_2.1.2          Biostrings_2.60.0     
[49] ade4_1.7-17            compiler_4.1.1         GenomeInfoDb_1.28.0   
[52] qs_0.25.1              rlang_0.4.11           rhdf5_2.36.0          
[55] grid_4.1.1             RCurl_1.98-1.5         pbdZMQ_0.3-5          
[58] iterators_1.0.13       IRkernel_1.2           rhdf5filters_1.4.0    
[61] biomformat_1.20.0      igraph_1.2.6           bitops_1.0-7          
[64] base64enc_0.1-3        gtable_0.3.0           codetools_0.2-18      
[67] multtest_2.48.0        reshape2_1.4.4         R6_2.5.1              
[70] fastmap_1.1.0          utf8_1.2.2             stringi_1.7.4         
[73] parallel_4.1.1         IRdisplay_1.0          Rcpp_1.0.7            
[76] vctrs_0.3.8            tidyselect_1.1.1 

[nick-youngblut]

This does not produce an error:

biom <- select(hmp50, 1:10)
dm <- unifrac(biom)

[nick-youngblut]

The error seems to be due to:

metrics <- function (biom, mode) {

  mode <- tolower(mode)[[1]]
  if (!is(biom, 'BIOM') && mode %in% c('rank', 'taxon', 'meta'))
    stop("Please provide a BIOM object when using mode='", mode, "'")

  if        (mode == 'ord')   { c("PCoA", "tSNE", "NMDS") 
  } else if (mode == 'adiv')  { c("OTUs", "Shannon", "Chao1", "Simpson", "InvSimpson") 
  } else if (mode == 'dist')  { c("correlation", "euclidean", "maximum", "manhattan", "canberra", "binary", "minkowski") 
  } else if (mode == 'clust') { c("average", "ward", "mcquitty", "single", "median", "complete", "centroid") 
  } else if (mode == 'rank')  { taxa.ranks(biom)
  } else if (mode == 'taxon') { unique(c(as.character(taxonomy(biom)), taxa.names(biom)))
  } else if (mode == 'meta')  { colnames(metadata(biom))
  } else if (mode == 'bdiv')  { 
    if (has.phylogeny(biom)) { c("Manhattan", "Euclidean", "Bray-Curtis", "Jaccard", "UniFrac")
    } else                   { c("Manhattan", "Euclidean", "Bray-Curtis", "Jaccard") }
  } else                      { NULL }
}

...where you are assuming the biom object is actually a biom-class object. Maybe !is(biom, 'BIOM') && mode %in% c('rank', 'taxon', 'meta') isn't actually checking for biom class due to && mode %in% c('rank', 'taxon', 'meta')?

[dansmith01]

Thanks for the heads up and investigation on this, @nick-youngblut. Fixed!

[nick-youngblut]

Thanks for the quick fix!

FYI: I'm planning on creating a conda-forge recipe for rbiom, unless you are planning on creating one in the near future.

[dansmith01]

That's a great idea! I didn't realize until just now that conda handles r packages too. I'm pushing CRAN's version of rbiom over there now and will update this thread when it's in place.

The next couple weeks should see the current development branch submitted to CRAN, and I'll update conda-forge then as well.

[nick-youngblut]

Great! Yeah, conda (especially with mamba) can be a one-stop-shop for installing all bioinformatics software than one generally needs, including R packages and other non-python software. I highly recommend it.

I'm looking forward to seeing rbiom on conda-forge!

[dansmith01]

@nick-youngblut Now on conda-forge! https://anaconda.org/conda-forge/r-rbiom

[nick-youngblut]

Awesome!! 🎉 Thanks for adding rbiom to conda-forge!