library(MiscMetabar)
data(data_fungi)Alpha diversity
Hill number
Numerous metrics of diversity exist. Hill numbers 1 is a kind of general framework for alpha diversity index.
renyi_res <- vegan::renyi(data_fungi@otu_table)
head(renyi_res)
#> 0 0.25 0.5 1 2
#> A10-005-B_S188_MERGED.fastq.gz 4.204693 3.615323 3.044244 2.0754183 1.1561862
#> A10-005-H_S189_MERGED.fastq.gz 4.248495 3.150337 2.063712 0.9938545 0.6464533
#> A10-005-M_S190_MERGED.fastq.gz 3.988984 3.607773 3.308361 2.9304156 2.5781603
#> A12-007_S191_MERGED.fastq.gz 5.036953 4.423283 3.751946 2.6309045 1.8468497
#> A12-007-B_S2_MERGED.fastq.gz 3.850148 3.349461 2.874269 2.1881145 1.6533262
#> A15-004_S3_MERGED.fastq.gz 4.025352 3.890614 3.747555 3.4533352 2.9618789
#> 4 8 16 32
#> A10-005-B_S188_MERGED.fastq.gz 0.8007816 0.6864775 0.6407124 0.6200442
#> A10-005-H_S189_MERGED.fastq.gz 0.5146170 0.4469617 0.4172228 0.4037640
#> A10-005-M_S190_MERGED.fastq.gz 2.3022639 2.1015018 1.9744485 1.9110311
#> A12-007_S191_MERGED.fastq.gz 1.5920328 1.4947479 1.4421368 1.4109287
#> A12-007-B_S2_MERGED.fastq.gz 1.3979355 1.2774402 1.2106906 1.1733049
#> A15-004_S3_MERGED.fastq.gz 2.4997916 2.2378349 2.1045646 2.0377518
#> 64 Inf
#> A10-005-B_S188_MERGED.fastq.gz 0.6102022 0.6006678
#> A10-005-H_S189_MERGED.fastq.gz 0.3973550 0.3911464
#> A10-005-M_S190_MERGED.fastq.gz 1.8806976 1.8513117
#> A12-007_S191_MERGED.fastq.gz 1.3916365 1.3700776
#> A12-007-B_S2_MERGED.fastq.gz 1.1547034 1.1366612
#> A15-004_S3_MERGED.fastq.gz 2.0054156 1.9740810Test for difference in diversity (hill number)
One way to keep into account for difference in the number of sequences per samples is to use a Tukey test on a linear model with the square roots of the number of sequence as the first explanatory variable of the linear model 2.
p <- MiscMetabar::hill_pq(data_fungi, variable = "Height")
p$plot_Hill_0
Hill number 1
p$plot_tuckey
Result of the Tuckey post-hoc test
See also the tutorial of the microbiome package for an alternative using the non-parametric Kolmogorov-Smirnov test for two-group comparisons when there are no relevant covariates.
Alpha diversity using package MicrobiotaProcess
library("MicrobiotaProcess")
library("ggh4x")
clean_pq(subset_samples_pq(data_fungi, !is.na(data_fungi@sam_data$Height))) %>%
as.MPSE() %>%
mp_cal_alpha() %>%
mp_plot_alpha(.group = "Height")
#> Warning in wilcox.test.default(c(4, 4, 5, 5, 4, 4, 3, 6, 3, 4, 3, 2, 3, :
#> cannot compute exact p-value with ties
#> Warning in wilcox.test.default(c(3, 4, 3, 5, 4, 6, 3, 5, 4, 4, 3, 3, 3, :
#> cannot compute exact p-value with ties
#> Warning in wilcox.test.default(c(3, 4, 3, 5, 4, 6, 3, 5, 4, 4, 3, 3, 3, :
#> cannot compute exact p-value with ties
#> Warning in wilcox.test.default(c(1.32966134885476, 1.242453324894,
#> 1.56071040904141, : cannot compute exact p-value with ties
#> Warning in wilcox.test.default(c(0.867563228481461, 1.32966134885476,
#> 0.867563228481461, : cannot compute exact p-value with ties
#> Warning in wilcox.test.default(c(0.867563228481461, 1.32966134885476,
#> 0.867563228481461, : cannot compute exact p-value with ties
Session information
sessionInfo()
#> R version 4.3.3 (2024-02-29)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Debian GNU/Linux 11 (bullseye)
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.13.so; LAPACK version 3.9.0
#>
#> locale:
#> [1] LC_CTYPE=fr_FR.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=fr_FR.UTF-8 LC_COLLATE=fr_FR.UTF-8
#> [5] LC_MONETARY=fr_FR.UTF-8 LC_MESSAGES=fr_FR.UTF-8
#> [7] LC_PAPER=fr_FR.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: Europe/Paris
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] ggh4x_0.2.8 MicrobiotaProcess_1.14.1 MiscMetabar_0.8.00
#> [4] purrr_1.0.2 dplyr_1.1.4 dada2_1.30.0
#> [7] Rcpp_1.0.12 ggplot2_3.5.0 phyloseq_1.46.0
#>
#> loaded via a namespace (and not attached):
#> [1] libcoin_1.0-10 RColorBrewer_1.1-3
#> [3] rstudioapi_0.16.0 jsonlite_1.8.8
#> [5] magrittr_2.0.3 TH.data_1.1-2
#> [7] modeltools_0.2-23 farver_2.1.1
#> [9] rmarkdown_2.26 fs_1.6.3
#> [11] zlibbioc_1.48.2 ragg_1.3.0
#> [13] vctrs_0.6.5 multtest_2.58.0
#> [15] memoise_2.0.1 Rsamtools_2.18.0
#> [17] RCurl_1.98-1.14 ggtree_3.10.1
#> [19] htmltools_0.5.8.1 S4Arrays_1.2.1
#> [21] Rhdf5lib_1.24.2 gridGraphics_0.5-1
#> [23] SparseArray_1.2.4 rhdf5_2.46.1
#> [25] sass_0.4.9 bslib_0.7.0
#> [27] desc_1.4.3 plyr_1.8.9
#> [29] sandwich_3.1-0 zoo_1.8-12
#> [31] cachem_1.0.8 GenomicAlignments_1.38.2
#> [33] igraph_2.0.3 lifecycle_1.0.4
#> [35] iterators_1.0.14 pkgconfig_2.0.3
#> [37] Matrix_1.6-5 R6_2.5.1
#> [39] fastmap_1.1.1 GenomeInfoDbData_1.2.11
#> [41] MatrixGenerics_1.14.0 aplot_0.2.2
#> [43] digest_0.6.35 ggnewscale_0.4.10
#> [45] colorspace_2.1-0 ShortRead_1.60.0
#> [47] patchwork_1.2.0 S4Vectors_0.40.2
#> [49] textshaping_0.3.7 GenomicRanges_1.54.1
#> [51] hwriter_1.3.2.1 vegan_2.6-4
#> [53] labeling_0.4.3 fansi_1.0.6
#> [55] abind_1.4-5 mgcv_1.9-1
#> [57] compiler_4.3.3 withr_3.0.0
#> [59] BiocParallel_1.36.0 highr_0.10
#> [61] ggsignif_0.6.4 MASS_7.3-60.0.1
#> [63] DelayedArray_0.28.0 biomformat_1.30.0
#> [65] permute_0.9-7 tools_4.3.3
#> [67] ape_5.8 glue_1.7.0
#> [69] nlme_3.1-164 rhdf5filters_1.14.1
#> [71] grid_4.3.3 cluster_2.1.6
#> [73] reshape2_1.4.4 ade4_1.7-22
#> [75] generics_0.1.3 gtable_0.3.4
#> [77] tidyr_1.3.1 data.table_1.15.4
#> [79] coin_1.4-3 utf8_1.2.4
#> [81] XVector_0.42.0 BiocGenerics_0.48.1
#> [83] ggrepel_0.9.5 foreach_1.5.2
#> [85] pillar_1.9.0 stringr_1.5.1
#> [87] yulab.utils_0.1.4 splines_4.3.3
#> [89] gghalves_0.1.4 treeio_1.26.0
#> [91] lattice_0.22-6 survival_3.5-8
#> [93] deldir_2.0-4 tidyselect_1.2.1
#> [95] Biostrings_2.70.3 knitr_1.46
#> [97] gridExtra_2.3 IRanges_2.36.0
#> [99] SummarizedExperiment_1.32.0 ggtreeExtra_1.12.0
#> [101] stats4_4.3.3 xfun_0.43
#> [103] Biobase_2.62.0 matrixStats_1.3.0
#> [105] stringi_1.8.3 lazyeval_0.2.2
#> [107] ggfun_0.1.4 yaml_2.3.8
#> [109] evaluate_0.23 codetools_0.2-19
#> [111] interp_1.1-6 tibble_3.2.1
#> [113] ggplotify_0.1.2 cli_3.6.2
#> [115] RcppParallel_5.1.7 systemfonts_1.0.6
#> [117] munsell_0.5.1 jquerylib_0.1.4
#> [119] GenomeInfoDb_1.38.8 png_0.1-8
#> [121] parallel_4.3.3 pkgdown_2.0.7
#> [123] latticeExtra_0.6-30 jpeg_0.1-10
#> [125] bitops_1.0-7 ggstar_1.0.4
#> [127] mvtnorm_1.2-4 tidytree_0.4.6
#> [129] scales_1.3.0 crayon_1.5.2
#> [131] rlang_1.1.3 multcomp_1.4-25