Skip to contents
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.133
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")

Circle of ASVs

circle_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

library("grid")
venn_pq(data_fungi, fact = "Height")

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.

Library requirement for Debian Linux OS

sudo apt-get install libgsl-dev libmpfr-dev

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.9.1 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