yiluheihei / microbiomeMarker

R package for microbiome biomarker discovery
https://yiluheihei.github.io/microbiomeMarker
GNU General Public License v3.0
171 stars 40 forks source link

plot_cladogram error #57

Open SonWende opened 2 years ago

SonWende commented 2 years ago

Hello,

Thank you for your support and creating and maintaining this package.

When I try to make cladograms for more my data after running lefse with the command plot_cladogram(mm.lefse.root.noT, color=c("red","green", "yellow", "blue")) i get the following error message:

Error in if (n <= length(x)) { : Missing Value, where TRUE/FALSE is required
3.
get_unique_id(sum(ind))
2.
get_short_label_id(clade_label, clade_label_level)
1.
plot_cladogram(mm.lefse.root.noT, color = c("red", "green", "yellow", "blue")) 

When running plot_cladogram() with the supplied kostic_crc testdata is works fine.

I attached an rds file with the mm object used. example.zip

Here is the Session Info:


R version 4.1.2 (2021-11-01)

Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 18363)

Matrix products: default

locale:
[1] LC_COLLATE=German_Germany.1252  LC_CTYPE=German_Germany.1252    LC_MONETARY=German_Germany.1252 LC_NUMERIC=C                   
[5] LC_TIME=German_Germany.1252    

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

other attached packages:
 [1] ade4_1.7-18            microbiome_1.16.0      reshape2_1.4.4         vegan_2.5-7            lattice_0.20-45        permute_0.9-7         
 [7] forcats_0.5.1          stringr_1.4.0          purrr_0.3.4            readr_2.1.2            tidyr_1.2.0            tibble_3.1.6          
[13] ggplot2_3.3.5          tidyverse_1.3.1        dplyr_1.0.8            phyloseq_1.38.0        ANCOMBC_1.4.0          microbiomeMarker_1.1.2
[19] qiime2R_0.99.6         BiocManager_1.30.16   

loaded via a namespace (and not attached):
  [1] utf8_1.2.2                  tidyselect_1.1.2            RSQLite_2.2.10              AnnotationDbi_1.56.2        htmlwidgets_1.5.4          
  [6] grid_4.1.2                  BiocParallel_1.28.3         Rtsne_0.15                  devtools_2.4.3              munsell_0.5.0              
 [11] codetools_0.2-18            DT_0.21                     withr_2.5.0                 colorspace_2.0-3            Biobase_2.54.0             
 [16] knitr_1.37                  rstudioapi_0.13             stats4_4.1.2                MatrixGenerics_1.6.0        Rdpack_2.1.4               
 [21] GenomeInfoDbData_1.2.7      plotROC_2.2.1               bit64_4.0.5                 rhdf5_2.38.0                rprojroot_2.0.2            
 [26] vctrs_0.3.8                 treeio_1.18.1               generics_0.1.2              xfun_0.30                   R6_2.5.1                   
 [31] doParallel_1.0.17           GenomeInfoDb_1.30.1         clue_0.3-60                 locfit_1.5-9.5              bitops_1.0-7               
 [36] rhdf5filters_1.6.0          cachem_1.0.6                gridGraphics_0.5-1          DelayedArray_0.20.0         assertthat_0.2.1           
 [41] scales_1.1.1                nnet_7.3-16                 gtable_0.3.0                processx_3.5.2              rlang_1.0.2                
 [46] genefilter_1.76.0           GlobalOptions_0.1.2         splines_4.1.2               lazyeval_0.2.2              broom_0.7.12               
 [51] checkmate_2.0.0             modelr_0.1.8                yaml_2.3.5                  backports_1.4.1             Hmisc_4.6-0                
 [56] tools_4.1.2                 usethis_2.1.5               zCompositions_1.4.0         ggplotify_0.1.0             ellipsis_0.3.2             
 [61] gplots_3.1.1                biomformat_1.22.0           RColorBrewer_1.1-2          BiocGenerics_0.40.0         sessioninfo_1.2.2          
 [66] Rcpp_1.0.8                  plyr_1.8.6                  base64enc_0.1-3             zlibbioc_1.40.0             RCurl_1.98-1.6             
 [71] ps_1.6.0                    prettyunits_1.1.1           rpart_4.1-15                Wrench_1.12.0               GetoptLong_1.0.5           
 [76] S4Vectors_0.32.3            haven_2.4.3                 SummarizedExperiment_1.24.0 cluster_2.1.2               fs_1.5.2                   
 [81] magrittr_2.0.2              data.table_1.14.2           circlize_0.4.14             reprex_2.0.1                truncnorm_1.0-8            
 [86] matrixStats_0.61.0          pkgload_1.2.4               hms_1.1.1                   patchwork_1.1.1             xtable_1.8-4               
 [91] XML_3.99-0.9                jpeg_0.1-9                  readxl_1.3.1                IRanges_2.28.0              gridExtra_2.3              
 [96] shape_1.4.6                 testthat_3.1.2              compiler_4.1.2              KernSmooth_2.23-20          crayon_1.5.0               
[101] htmltools_0.5.2             tzdb_0.2.0                  ggfun_0.0.5                 mgcv_1.8-38                 Formula_1.2-4              
[106] geneplotter_1.72.0          aplot_0.1.2                 lubridate_1.8.0             DBI_1.1.2                   dbplyr_2.1.1               
[111] ComplexHeatmap_2.10.0       MASS_7.3-54                 Matrix_1.3-4                brio_1.1.3                  cli_3.2.0                  
[116] rbibutils_2.2.7             parallel_4.1.2              igraph_1.2.11               GenomicRanges_1.46.1        pkgconfig_2.0.3            
[121] foreign_0.8-81              xml2_1.3.3                  foreach_1.5.2               ggtree_3.2.1                annotate_1.72.0            
[126] multtest_2.50.0             XVector_0.34.0              rvest_1.0.2                 NADA_1.6-1.1                yulab.utils_0.0.4          
[131] callr_3.7.0                 digest_0.6.29               Biostrings_2.62.0           cellranger_1.1.0            tidytree_0.3.9             
[136] htmlTable_2.4.0             curl_4.3.2                  gtools_3.9.2                rjson_0.2.21                nloptr_2.0.0               
[141] lifecycle_1.0.1             nlme_3.1-153                jsonlite_1.8.0              Rhdf5lib_1.16.0             desc_1.4.1                 
[146] limma_3.50.1                fansi_1.0.2                 pillar_1.7.0                KEGGREST_1.34.0             fastmap_1.1.0              
[151] httr_1.4.2                  pkgbuild_1.3.1              survival_3.2-13             glue_1.6.2                  remotes_2.4.2              
[156] png_0.1-7                   iterators_1.0.14            glmnet_4.1-3                bit_4.0.4                   stringi_1.7.6              
[161] metagenomeSeq_1.36.0        blob_1.2.2                  DESeq2_1.34.0               latticeExtra_0.6-29         caTools_1.18.2             
[166] memoise_2.0.1      ```

Maybe you can help 
Thanks a lot in advance
SonWende commented 2 years ago

Okay, its working for me now

Just to give you a little update, Here' s what i managed to find out:

I noticed that the tax_table() from the marker results did not contain the full taxonomy string when i ran lefse with a tax_rank like this

run_lefse(ps,
                taxa_rank = "Family",
                wilcoxon_cutoff = 0.01,
                group = "Treatment", 
                multigrp_strat = TRUE, 
                kw_cutoff = 0.01,
                norm = "CPM"))

For all other plots except plot_cladogram() this did not pose an issue,but i guess without the full tax string there can be no cladogram

When i ran a tax_glom() before lefse, the marker tax_table() had the full path, but did not actually work on the desired Family level. tax_glom() aggregates taxa on a certain levels, but the higher taxonomic levels still remain as columns of "NA" in the taxonomy table, which seems to cause issues.

tax_glom(ps, "Family") %>% run_lefse(
                wilcoxon_cutoff = 0.01,
                group = "Trichoderma", 
                multigrp_strat = TRUE, 
                kw_cutoff = 0.01,
                norm = "CPM"))

So the removal of the leftover columns after running tax_glom fixed this issue for me.

my_tax_glom <- function(ps, level = "Family"){
  cols <-  colnames(tax_table(ps))
  remove_cols <- cols[seq(from = which(level == cols)+1, to = length(cols))]
  prep <- tax_glom(ps, level) 
  tax_table(prep) <- tax_table(prep) %>% data.frame %>% select(-remove_cols) %>% as.matrix 
  return(prep)
}

mm.lefse<- my_tax_glom(ps, "Family) %>% run_lefse(...)

All the best wishes, Sonja