library(MiscMetabar)
data(data_fungi)Permanova
data_fungi_woNA4height <- subset_samples(data_fungi, !is.na(data_fungi@sam_data$Height))
res_ado <- adonis_pq(data_fungi_woNA4height, "Tree_name+Height")
knitr::kable(res_ado)| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| Tree_name | 61 | 36.0741288 | 0.5746127 | 1.532526 | 0.001 |
| Height | 2 | 0.8514605 | 0.0135626 | 1.103256 | 0.134 |
| Residual | 67 | 25.8543072 | 0.4118246 | NA | NA |
| Total | 130 | 62.7798964 | 1.0000000 | NA | NA |
Graph Test
data_fungi_woNA4height <- subset_samples(data_fungi, !is.na(data_fungi@sam_data$Height))
graph_test_pq(data_fungi_woNA4height, "Height")
Compare two (group of) samples
Biplot
data_fungi_low_high <- subset_samples(
data_fungi,
data_fungi@sam_data$Height %in%
c("Low", "High")
)
data_fungi_low_high <- subset_taxa_pq(
data_fungi_low_high,
taxa_sums(data_fungi_low_high) > 5000
)
biplot_pq(data_fungi_low_high, fact = "Height", merge_sample_by = "Height")
#> Warning: Removed 1 row containing missing values or values outside the scale range
#> (`geom_rect()`).
Compare two (group of) samples with a table
compare_pairs_pq(data_fungi_low_high,
bifactor = "Height",
merge_sample_by = "Height",
modality = "Time"
)
#> # A tibble: 4 × 13
#> modality nb_ASV_High nb_ASV_Low nb_shared_ASV div_High div_Low nb_shared_seq
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 0 12 16 9 1.8 1.37 57639
#> 2 5 20 18 14 1.95 1.98 76006
#> 3 10 11 13 10 1.18 1.25 47042
#> 4 15 17 19 12 2 2.04 161348
#> # ℹ 6 more variables: percent_shared_seq_High <dbl>,
#> # percent_shared_seq_Low <dbl>, percent_shared_ASV_High <dbl>,
#> # percent_shared_ASV_Low <dbl>, ratio_nb_High_Low <dbl>,
#> # ratio_div_High_Low <dbl>Venn diagram
ggvenn_pq(data_fungi, fact = "Height") +
ggplot2::scale_fill_distiller(palette = "BuPu", direction = 1) +
labs(title = "Share number of ASV among Height in tree")
ggvenn_pq(data_fungi, fact = "Height", min_nb_seq = 5000) +
ggplot2::scale_fill_distiller(palette = "BuPu", direction = 1) +
labs(title = "Share number of ASV with more than 5000 seqs")
ggvenn_pq(data_fungi,
fact = "Height", taxonomic_rank = "Genus",
min_nb_seq = 100
) +
ggplot2::scale_fill_distiller(palette = "BuPu", direction = 1) +
labs(title = "Share number of Genus represented by at least one ASV with more than 100 seqs")
Upset plot
Venn diagram can quickly become complex to read when the number of modalities increase. One graphical solution is upset plot. MiscMetabar propose a solution based on the package ComplexUpset.
upset_pq(data_fungi, fact = "Height")
upset_pq(data_fungi, fact = "Time")
ComplexUpset package allow powerful configuration of you
plot as you can see in the following figure.
upset_pq(
data_fungi,
fact = "Time",
width_ratio = 0.2,
annotations = list(
"Sequences per ASV \n (log10)" = (
ggplot(mapping = aes(y = log10(Abundance)))
+
geom_jitter(aes(
color =
Abundance
), na.rm = TRUE)
+
geom_violin(alpha = 0.5, na.rm = TRUE) +
theme(legend.key.size = unit(0.2, "cm")) +
theme(axis.text = element_text(size = 12))
),
"ASV per phylum" = (
ggplot(mapping = aes(fill = Phylum))
+
geom_bar() +
ylab("ASV per phylum") +
theme(legend.key.size = unit(0.2, "cm")) +
theme(axis.text = element_text(size = 12))
)
)
)
Change in abundance across a factor
Benchdamic
There is a lot of available methods. Please refer to R package benchdamic for a list of method and a implementation of a benchmark for your data.
Using Deseq2 package
data("GlobalPatterns", package = "phyloseq")
GP <- subset_samples(
GlobalPatterns,
GlobalPatterns@sam_data$SampleType %in% c("Soil", "Skin")
)
plot_deseq2_pq(GP, c("SampleType", "Soil", "Skin"), pval = 0.001)
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=fr_FR.UTF-8
#> [9] LC_ADDRESS=fr_FR.UTF-8 LC_TELEPHONE=fr_FR.UTF-8
#> [11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=fr_FR.UTF-8
#>
#> time zone: Europe/Paris
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] grid stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] MiscMetabar_0.8.00 purrr_1.0.2 dplyr_1.1.4 dada2_1.30.0
#> [5] Rcpp_1.0.12 ggplot2_3.5.0 phyloseq_1.46.0
#>
#> loaded via a namespace (and not attached):
#> [1] bitops_1.0-7 pbapply_1.7-2
#> [3] deldir_2.0-4 gridExtra_2.3
#> [5] permute_0.9-7 rlang_1.1.3
#> [7] magrittr_2.0.3 ade4_1.7-22
#> [9] matrixStats_1.3.0 compiler_4.3.3
#> [11] mgcv_1.9-1 png_0.1-8
#> [13] systemfonts_1.0.6 vctrs_0.6.5
#> [15] reshape2_1.4.4 phyloseqGraphTest_0.1.1
#> [17] stringr_1.5.1 shape_1.4.6.1
#> [19] pkgconfig_2.0.3 crayon_1.5.2
#> [21] fastmap_1.1.1 XVector_0.42.0
#> [23] labeling_0.4.3 utf8_1.2.4
#> [25] Rsamtools_2.18.0 rmarkdown_2.26
#> [27] ragg_1.3.0 xfun_0.43
#> [29] zlibbioc_1.48.2 cachem_1.0.8
#> [31] GenomeInfoDb_1.38.8 jsonlite_1.8.8
#> [33] biomformat_1.30.0 highr_0.10
#> [35] rhdf5filters_1.14.1 DelayedArray_0.28.0
#> [37] Rhdf5lib_1.24.2 BiocParallel_1.36.0
#> [39] jpeg_0.1-10 parallel_4.3.3
#> [41] cluster_2.1.6 R6_2.5.1
#> [43] bslib_0.7.0 stringi_1.8.3
#> [45] RColorBrewer_1.1-3 ComplexUpset_1.3.3
#> [47] GenomicRanges_1.54.1 jquerylib_0.1.4
#> [49] SummarizedExperiment_1.32.0 iterators_1.0.14
#> [51] knitr_1.46 IRanges_2.36.0
#> [53] Matrix_1.6-5 splines_4.3.3
#> [55] igraph_2.0.3 tidyselect_1.2.1
#> [57] rstudioapi_0.16.0 abind_1.4-5
#> [59] yaml_2.3.8 ggVennDiagram_1.5.2
#> [61] vegan_2.6-4 codetools_0.2-19
#> [63] hwriter_1.3.2.1 lattice_0.22-6
#> [65] tibble_3.2.1 plyr_1.8.9
#> [67] Biobase_2.62.0 withr_3.0.0
#> [69] ShortRead_1.60.0 evaluate_0.23
#> [71] desc_1.4.3 survival_3.5-8
#> [73] rJava_1.0-11 RcppParallel_5.1.7
#> [75] circlize_0.4.16 Biostrings_2.70.3
#> [77] pillar_1.9.0 MatrixGenerics_1.14.0
#> [79] foreach_1.5.2 stats4_4.3.3
#> [81] generics_0.1.3 RCurl_1.98-1.14
#> [83] S4Vectors_0.40.2 munsell_0.5.1
#> [85] scales_1.3.0 ggnetwork_0.5.13
#> [87] glue_1.7.0 tools_4.3.3
#> [89] interp_1.1-6 data.table_1.15.4
#> [91] locfit_1.5-9.9 GenomicAlignments_1.38.2
#> [93] fs_1.6.3 rhdf5_2.46.1
#> [95] tidyr_1.3.1 ape_5.8
#> [97] latticeExtra_0.6-30 colorspace_2.1-0
#> [99] patchwork_1.2.0 nlme_3.1-164
#> [101] GenomeInfoDbData_1.2.11 cli_3.6.2
#> [103] textshaping_0.3.7 fansi_1.0.6
#> [105] S4Arrays_1.2.1 gtable_0.3.4
#> [107] DESeq2_1.42.1 sass_0.4.9
#> [109] digest_0.6.35 BiocGenerics_0.48.1
#> [111] SparseArray_1.2.4 farver_2.1.1
#> [113] venneuler_1.1-4 memoise_2.0.1
#> [115] htmltools_0.5.8.1 pkgdown_2.0.7
#> [117] multtest_2.58.0 lifecycle_1.0.4
#> [119] GlobalOptions_0.1.2 MASS_7.3-60.0.1