Chapter 5 Pre-processing

Performing the pre-processing steps on the phyloseq-class object. Data Processing is a critical procedure in data analysis. In this chapter, we would use the following steps to transform, normalize, and trim microbiota data.

Outline of this Chapter:

5.1 Loading Packages

library(XMAS2)
library(dplyr)
library(tibble)
library(phyloseq)
library(ggplot2)
library(ggpubr)
library(conflicted)
conflicted::conflict_prefer("normalize", "XMAS2")

5.2 Importing Data

data("dada2_ps")
dada2_ps_remove_BRS <- get_GroupPhyloseq(
                         ps = dada2_ps,
                         group = "Group",
                         group_names = "QC",
                         discard = TRUE)
dada2_ps_remove_BRS
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 896 taxa and 23 samples ]
## sample_data() Sample Data:       [ 23 samples by 1 sample variables ]
## tax_table()   Taxonomy Table:    [ 896 taxa by 7 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 896 tips and 893 internal nodes ]
## refseq()      DNAStringSet:      [ 896 reference sequences ]
metagenomic sequencing (metaphlan2/3)
data("metaphlan2_ps")
metaphlan2_ps_remove_BRS <- get_GroupPhyloseq(
                               ps = metaphlan2_ps,
                               group = "Group",
                               group_names = "QC",
                               discard = TRUE)
metaphlan2_ps_remove_BRS
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 326 taxa and 22 samples ]
## sample_data() Sample Data:       [ 22 samples by 2 sample variables ]
## tax_table()   Taxonomy Table:    [ 326 taxa by 7 taxonomic ranks ]

5.3 Transformation

Transforming the taxa abundances in otu_table by sample, which means the counts of each sample will be transformed individually. The options include:

  • “identity”, return the original data without any transformation;

  • “log10”, the transformation is log10(object), and if the data contains zeros the transformation is log10(1 + object);

  • “log10p”, the transformation is log10(1 + object).

dada2_ps_genus <- summarize_taxa(ps = dada2_ps_remove_BRS, 
                                 taxa_level = "Genus")

dada2_ps_transform <- transform_abundances(object = dada2_ps_genus, 
                                           transform = "log10")
head(dada2_ps_transform@otu_table@.Data, 3)
##                        S6030 S6032 S6033 S6035    S6036    S6037 S6040 S6043    S6045 S6046 S6048   S6049 S6050   S6054 S6055 S6058 S6059
## g__Acetanaerobacterium     0     0     0     0 0.000000 0.000000     0     0 0.000000     0     0 0.00000     0 0.30103     0     0     0
## g__Acidaminococcus         0     0     0     0 2.481443 3.066326     0     0 2.326336     0     0 0.00000     0 0.00000     0     0     0
## g__Acinetobacter           0     0     0     0 0.000000 0.000000     0     0 0.000000     0     0 0.69897     0 0.00000     0     0     0
##                          S6060 S6061    S6063 S6065 S6066 S6068
## g__Acetanaerobacterium 0.00000     0 0.000000     0     0     0
## g__Acidaminococcus     0.60206     0 1.812913     0     0     0
## g__Acinetobacter       0.00000     0 0.000000     0     0     0
transforming metagenomic sequencing
metaphlan2_ps_species <- summarize_taxa(ps = metaphlan2_ps_remove_BRS, 
                                        taxa_level = "Species")

metaphlan2_ps_transform_mgs <- transform_abundances(object = metaphlan2_ps_species, 
                                                    transform = "log10")
head(metaphlan2_ps_transform_mgs@otu_table@.Data, 3)
##                                      s1 s2 s3 s4 s5 s6        s7        s8        s9 s10 s11 s12 s13 s14 s15       s16 s17 s18 s19 s20 s21
## s__Abiotrophia_defectiva      -4.610834  0  0  0  0  0  0.000000  0.000000  0.000000   0   0   0   0   0   0  0.000000   0   0   0   0   0
## s__Acidaminococcus_fermentans  0.000000  0  0  0  0  0  0.000000 -2.957621 -2.180555   0   0   0   0   0   0 -2.326583   0   0   0   0   0
## s__Acidaminococcus_intestini  -5.376751  0  0  0  0  0 -3.222573 -2.997748  0.000000   0   0   0   0   0   0  0.000000   0   0   0   0   0
##                                     s22
## s__Abiotrophia_defectiva       0.000000
## s__Acidaminococcus_fermentans -2.616418
## s__Acidaminococcus_intestini  -2.831945

5.4 Normalization

Normalizing the OTU_table in phyloseq-class object. It is critical to normalize the feature table to eliminate any bias due to differences in the sampling sequencing depth. This function implements 7 widely-used normalization methods for microbial compositional data.

The corresponding XMAS2 argument is normMethod with the available options: “rarefy”, “TSS”, “TMM”, “RLE”, “CSS”, “CLR” and “CPM”.

Method Approach Comments
rarefy (Rarefying) Subsampling from the data to obtain samples with equal library size (also called rarefaction level) . Aims to avoid compositional effects.
TSS (Total Sum Scaling) Traditional approach for building fractions. Counts are divided by the total sum of counts in the corresponding sample. Strongly influenced by highly abundant taxa.
TMM (Trimmed Mean of M-values) First, a sample is chosen as reference. The scaling factor is then derived using a weighted trimmed mean over the differences of the log-transformed gene-count fold-change between the sample and the reference. Aims to avoid sequencing depth effects.
RLE (Relative Log Expression ) RLE uses a pseudo-reference calculated using the geometric mean of the gene-specific abundances over all samples. The scaling factors are then calculated as the median of the gene counts ratios between the samples and the reference. Similar to TMM, this normalization method is based on the hypothesis that the most genes are not DE.
CSS (Cumulative Sum Scaling) Within each sample, the counts are summed up to a predefined quantile. The counts are then divided by this sum. Aims at avoiding the influence of highly abundant taxa.
CLR (Centered log-ratio transformation) For a composition \[x = (x_{1} , ..., x_{p})\] the clr transformation is defined as \[clr = log(\frac{x_{1}}{g(x)}, ..., \frac{x_{p}}{g(x)})\], where \[g(x)=(\prod_{p}^{k=1} x_{k})^{\frac{1}{p}}\] is the geometric mean. Aims to avoid compositional effects.
CPM (Counts Per Million) pre-sample normalization of the sum of the values to 1e+06. Aims to avoid sequencing depth effects.
  • rarefy: random subsampling counts to the smallest library size in the data set. Caution: the default library size is 25000 according to our own results(Rarefaction Curves).
XMAS2::normalize(object = dada2_ps_genus, 
                 method = "rarefy")
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 198 taxa and 23 samples ]
## sample_data() Sample Data:       [ 23 samples by 1 sample variables ]
## tax_table()   Taxonomy Table:    [ 198 taxa by 6 taxonomic ranks ]
# norm_rarefy(dada2_ps_genus, size = 50000)
  • TSS: total sum scaling, also referred to as “relative abundance”, the abundances were normalized by dividing the corresponding sample library size
XMAS2::normalize(object = dada2_ps_genus, 
                 method = "TSS")
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 198 taxa and 23 samples ]
## sample_data() Sample Data:       [ 23 samples by 1 sample variables ]
## tax_table()   Taxonomy Table:    [ 198 taxa by 6 taxonomic ranks ]
  • TMM: trimmed mean of m-values. First, a sample is chosen as reference. The scaling factor is then derived using a weighted trimmed mean over the differences of the log-transformed gene-count fold-change between the sample and the reference.
XMAS2::normalize(object = dada2_ps_genus, 
                 method = "TMM")
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 198 taxa and 23 samples ]
## sample_data() Sample Data:       [ 23 samples by 1 sample variables ]
## tax_table()   Taxonomy Table:    [ 198 taxa by 6 taxonomic ranks ]
  • RLE: relative log expression, RLE uses a pseudo-reference calculated using the geometric mean of the gene-specific abundances over all samples. The scaling factors are then calculated as the median of the gene counts ratios between the samples and the reference.
XMAS2::normalize(object = dada2_ps_genus, 
                 method = "RLE")
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 198 taxa and 23 samples ]
## sample_data() Sample Data:       [ 23 samples by 1 sample variables ]
## tax_table()   Taxonomy Table:    [ 198 taxa by 6 taxonomic ranks ]
  • CSS: cumulative sum scaling, calculates scaling factors as the cumulative sum of gene abundances up to a data-derived threshold. While standard relative abundance (fraction/percentage) normalization re-scales all samples to the same total sum (100%), CSS keeps a variation in total counts between samples. CSS re-scales the samples based on a subset (quartile) of lower abundant taxa (relatively constant and independent), thereby excluding the impact of (study dominating) high abundant taxa.
XMAS2::normalize(object = dada2_ps_genus, 
                 method = "CSS")
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 198 taxa and 23 samples ]
## sample_data() Sample Data:       [ 23 samples by 1 sample variables ]
## tax_table()   Taxonomy Table:    [ 198 taxa by 6 taxonomic ranks ]
  • CLR: centered log-ratio normalization.
XMAS2::normalize(object = dada2_ps_genus, 
                 method = "CLR")
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 198 taxa and 23 samples ]
## sample_data() Sample Data:       [ 23 samples by 1 sample variables ]
## tax_table()   Taxonomy Table:    [ 198 taxa by 6 taxonomic ranks ]
  • CPM: pre-sample normalization of the sum of the values to 1e+06.
XMAS2::normalize(object = dada2_ps_genus, 
                 method = "CPM")
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 198 taxa and 23 samples ]
## sample_data() Sample Data:       [ 23 samples by 1 sample variables ]
## tax_table()   Taxonomy Table:    [ 198 taxa by 6 taxonomic ranks ]

5.5 Filtering

Whether to filter the low relative abundance or unclassified taxa by the threshold.

ps_genus_rb <- summarize_taxa(ps = dada2_ps_remove_BRS, 
                              taxa_level = "Genus", 
                              absolute = FALSE)
ps_genus_rb
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 198 taxa and 23 samples ]
## sample_data() Sample Data:       [ 23 samples by 1 sample variables ]
## tax_table()   Taxonomy Table:    [ 198 taxa by 6 taxonomic ranks ]
ps_genus_rb_filter <- run_filter(ps = ps_genus_rb, 
                                 cutoff = 1e-04, 
                                 unclass = TRUE)
ps_genus_rb_filter
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 110 taxa and 23 samples ]
## sample_data() Sample Data:       [ 23 samples by 1 sample variables ]
## tax_table()   Taxonomy Table:    [ 110 taxa by 6 taxonomic ranks ]

39 genus’s relative abundance or attributes were below 1e-04 or unclassified and they were removed by the cutoff.

5.6 Trimming

The previous function (run_filter) only focuses on the low relative abundance and unclassified taxa. Microbial data always have so many zeros. Trimming samples or taxa in otu_table by occurrences or prevalence before downstream analysis is also crucial.

  • trimming by TaxaID
run_trim(object = dada2_ps_genus, 
         cutoff = 0.4, 
         trim = "feature")
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 63 taxa and 23 samples ]
## sample_data() Sample Data:       [ 23 samples by 1 sample variables ]
## tax_table()   Taxonomy Table:    [ 63 taxa by 6 taxonomic ranks ]

Dropping the taxa whose prevalence or occurrences are less than 0.4.

  • trimming by SampleID
run_trim(object = dada2_ps_genus, 
         cutoff = 0.4, 
         trim = "sample")
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 198 taxa and 4 samples ]
## sample_data() Sample Data:       [ 4 samples by 1 sample variables ]
## tax_table()   Taxonomy Table:    [ 198 taxa by 6 taxonomic ranks ]

Dropping the samples whose prevalence or occurrences are less than 0.4.

  • trimming by TaxaID & SampleID
run_trim(object = dada2_ps_genus, 
         cutoff = 0.4, 
         trim = "both")
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 63 taxa and 4 samples ]
## sample_data() Sample Data:       [ 4 samples by 1 sample variables ]
## tax_table()   Taxonomy Table:    [ 63 taxa by 6 taxonomic ranks ]

Dropping the taxa and samples whose prevalence are less than 0.4.

filtering metagenomic sequencing
metaphlan2_ps_trim_mgs <- run_trim(object = metaphlan2_ps_remove_BRS, 
                                   cutoff = 0.4, 
                                   trim = "feature")
metaphlan2_ps_trim_mgs
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 102 taxa and 22 samples ]
## sample_data() Sample Data:       [ 22 samples by 2 sample variables ]
## tax_table()   Taxonomy Table:    [ 102 taxa by 7 taxonomic ranks ]

5.7 Imputation

The missing values in otu table maybe affect the statistical results, imputing the NAs or Zero values should taken into account.

  • limit of detection
min(dada2_ps_genus@otu_table)
## [1] 0
LOD_imputed_ps <- run_impute(object = dada2_ps_genus, 
                             impute = "LOD", 
                             LOD = 2)
min(LOD_imputed_ps@otu_table)
## [1] 2

5.8 Extracting specific taxa phyloseq-class object

The taxonomic level are Kingdom, Phylum, Class, Order, Family, Genus, Species and choosing the specific taxa to regenerate the phyloseq-class object.

  • amplicon sequencing: Phylum
dada2_ps_phylum <- summarize_taxa(ps = dada2_ps_remove_BRS, 
                                  taxa_level = "Phylum")
dada2_ps_phylum
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 10 taxa and 23 samples ]
## sample_data() Sample Data:       [ 23 samples by 1 sample variables ]
## tax_table()   Taxonomy Table:    [ 10 taxa by 2 taxonomic ranks ]
  • amplicon sequencing: Order
dada2_ps_order <- summarize_taxa(ps = dada2_ps_remove_BRS, 
                                 taxa_level = "Order")
dada2_ps_order
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 24 taxa and 23 samples ]
## sample_data() Sample Data:       [ 23 samples by 1 sample variables ]
## tax_table()   Taxonomy Table:    [ 24 taxa by 4 taxonomic ranks ]
  • amplicon sequencing: Family
dada2_ps_family <- summarize_taxa(ps = dada2_ps_remove_BRS, 
                                  taxa_level = "Family")
dada2_ps_family
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 54 taxa and 23 samples ]
## sample_data() Sample Data:       [ 23 samples by 1 sample variables ]
## tax_table()   Taxonomy Table:    [ 54 taxa by 5 taxonomic ranks ]
  • amplicon sequencing: Genus
dada2_ps_genus <- summarize_taxa(ps = dada2_ps_remove_BRS, 
                                 taxa_level = "Genus")
dada2_ps_genus
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 198 taxa and 23 samples ]
## sample_data() Sample Data:       [ 23 samples by 1 sample variables ]
## tax_table()   Taxonomy Table:    [ 198 taxa by 6 taxonomic ranks ]
extracting metagenomic sequencing
metaphlan2_ps_genus <- summarize_taxa(ps = metaphlan2_ps_remove_BRS, 
                                      taxa_level = "Genus")
metaphlan2_ps_genus
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 109 taxa and 22 samples ]
## sample_data() Sample Data:       [ 22 samples by 2 sample variables ]
## tax_table()   Taxonomy Table:    [ 109 taxa by 6 taxonomic ranks ]
metaphlan2_ps_species <- summarize_taxa(ps = metaphlan2_ps_remove_BRS, 
                                        taxa_level = "Species")
metaphlan2_ps_species
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 325 taxa and 22 samples ]
## sample_data() Sample Data:       [ 22 samples by 2 sample variables ]
## tax_table()   Taxonomy Table:    [ 325 taxa by 7 taxonomic ranks ]

5.9 Aggregating low relative abundance or unclassified taxa into others

  1. Taxa with relative abundance less than 0.0001 will be summarized into Others_LowAbundance;

  2. Unclassified taxa will be summarized into Others_Unclassified.

  • amplicon sequencing
# relative abundance
dada2_ps_genus_rb <- summarize_taxa(ps = dada2_ps_remove_BRS, 
                                    taxa_level = "Genus", 
                                    absolute = FALSE)
dada2_ps_genus_LRA <- summarize_LowAbundance_taxa(ps = ps_genus_rb, 
                                                  cutoff = 1e-04, 
                                                  unclass = TRUE)
dada2_ps_genus_LRA 
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 162 taxa and 23 samples ]
## sample_data() Sample Data:       [ 23 samples by 1 sample variables ]
## tax_table()   Taxonomy Table:    [ 162 taxa by 6 taxonomic ranks ]
tail(phyloseq::taxa_names(dada2_ps_genus_LRA))
## [1] "g__Tyzzerella_3" "g__Tyzzerella_4" "g__UBA1819"      "g__UC5_1_2E3"    "g__Veillonella"  "g__Weissella"
# absolute abundance
dada2_ps_genus_counts <- summarize_LowAbundance_taxa(ps = dada2_ps_genus, 
                                                     cutoff = 10)
dada2_ps_genus_counts
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 173 taxa and 23 samples ]
## sample_data() Sample Data:       [ 23 samples by 1 sample variables ]
## tax_table()   Taxonomy Table:    [ 173 taxa by 6 taxonomic ranks ]
aggregating metagenomic sequencing
metaphlan2_ps_genus_LRA <- summarize_LowAbundance_taxa(ps = metaphlan2_ps_genus, 
                                                       cutoff = 1e-04, 
                                                       unclass = TRUE)
metaphlan2_ps_genus_LRA
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 78 taxa and 22 samples ]
## sample_data() Sample Data:       [ 22 samples by 2 sample variables ]
## tax_table()   Taxonomy Table:    [ 78 taxa by 6 taxonomic ranks ]
tail(phyloseq::taxa_names(metaphlan2_ps_genus_LRA))
## [1] "g__Subdoligranulum" "g__Sutterella"      "g__T4likevirus"     "g__Turicibacter"    "g__Veillonella"     "g__Weissella"

5.10 Transform the abundance of taxa whose relative abundance under Limit Of Detection (LOD) into Zeros (only in metaphlan2/3)

metaphlan2_ps_Species_LOD <- aggregate_LOD_taxa(ps = metaphlan2_ps, 
                                                taxa_level = "Species", 
                                                cutoff = 1e-04)
metaphlan2_ps_Species_LOD
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 192 taxa and 23 samples ]
## sample_data() Sample Data:       [ 23 samples by 2 sample variables ]
## tax_table()   Taxonomy Table:    [ 192 taxa by 7 taxonomic ranks ]

5.11 Systematic Information

devtools::session_info()
## ─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
##  setting  value
##  version  R version 4.1.3 (2022-03-10)
##  os       macOS Monterey 12.2.1
##  system   x86_64, darwin17.0
##  ui       RStudio
##  language (EN)
##  collate  en_US.UTF-8
##  ctype    en_US.UTF-8
##  tz       Asia/Shanghai
##  date     2023-11-30
##  rstudio  2023.09.0+463 Desert Sunflower (desktop)
##  pandoc   3.1.1 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/ (via rmarkdown)
## 
## ─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
##  package              * version   date (UTC) lib source
##  abind                  1.4-5     2016-07-21 [2] CRAN (R 4.1.0)
##  ade4                   1.7-22    2023-02-06 [2] CRAN (R 4.1.2)
##  annotate               1.72.0    2021-10-26 [2] Bioconductor
##  AnnotationDbi          1.60.2    2023-03-10 [2] Bioconductor
##  ape                    5.7-1     2023-03-13 [2] CRAN (R 4.1.2)
##  backports              1.4.1     2021-12-13 [2] CRAN (R 4.1.0)
##  base64enc              0.1-3     2015-07-28 [2] CRAN (R 4.1.0)
##  Biobase                2.54.0    2021-10-26 [2] Bioconductor
##  BiocGenerics           0.40.0    2021-10-26 [2] Bioconductor
##  BiocParallel           1.28.3    2021-12-09 [2] Bioconductor
##  biomformat             1.22.0    2021-10-26 [2] Bioconductor
##  Biostrings             2.62.0    2021-10-26 [2] Bioconductor
##  bit                    4.0.5     2022-11-15 [2] CRAN (R 4.1.2)
##  bit64                  4.0.5     2020-08-30 [2] CRAN (R 4.1.0)
##  bitops                 1.0-7     2021-04-24 [2] CRAN (R 4.1.0)
##  blob                   1.2.4     2023-03-17 [2] CRAN (R 4.1.2)
##  bookdown               0.34      2023-05-09 [2] CRAN (R 4.1.2)
##  broom                  1.0.5     2023-06-09 [2] CRAN (R 4.1.3)
##  bslib                  0.6.0     2023-11-21 [1] CRAN (R 4.1.3)
##  cachem                 1.0.8     2023-05-01 [2] CRAN (R 4.1.2)
##  callr                  3.7.3     2022-11-02 [2] CRAN (R 4.1.2)
##  car                    3.1-2     2023-03-30 [2] CRAN (R 4.1.2)
##  carData                3.0-5     2022-01-06 [2] CRAN (R 4.1.2)
##  caTools                1.18.2    2021-03-28 [2] CRAN (R 4.1.0)
##  checkmate              2.2.0     2023-04-27 [2] CRAN (R 4.1.2)
##  cli                    3.6.1     2023-03-23 [2] CRAN (R 4.1.2)
##  cluster                2.1.4     2022-08-22 [2] CRAN (R 4.1.2)
##  coda                   0.19-4    2020-09-30 [2] CRAN (R 4.1.0)
##  codetools              0.2-19    2023-02-01 [2] CRAN (R 4.1.2)
##  coin                   1.4-2     2021-10-08 [2] CRAN (R 4.1.0)
##  colorspace             2.1-0     2023-01-23 [2] CRAN (R 4.1.2)
##  conflicted           * 1.2.0     2023-02-01 [2] CRAN (R 4.1.2)
##  cowplot                1.1.1     2020-12-30 [2] CRAN (R 4.1.0)
##  crayon                 1.5.2     2022-09-29 [2] CRAN (R 4.1.2)
##  data.table             1.14.8    2023-02-17 [2] CRAN (R 4.1.2)
##  DBI                    1.1.3     2022-06-18 [2] CRAN (R 4.1.2)
##  DelayedArray           0.20.0    2021-10-26 [2] Bioconductor
##  DESeq2                 1.34.0    2021-10-26 [2] Bioconductor
##  devtools               2.4.5     2022-10-11 [2] CRAN (R 4.1.2)
##  digest                 0.6.33    2023-07-07 [1] CRAN (R 4.1.3)
##  dplyr                * 1.1.2     2023-04-20 [2] CRAN (R 4.1.2)
##  DT                     0.28      2023-05-18 [2] CRAN (R 4.1.3)
##  edgeR                  3.36.0    2021-10-26 [2] Bioconductor
##  ellipsis               0.3.2     2021-04-29 [2] CRAN (R 4.1.0)
##  emmeans                1.8.7     2023-06-23 [1] CRAN (R 4.1.3)
##  estimability           1.4.1     2022-08-05 [2] CRAN (R 4.1.2)
##  evaluate               0.21      2023-05-05 [2] CRAN (R 4.1.2)
##  FactoMineR             2.8       2023-03-27 [2] CRAN (R 4.1.2)
##  fansi                  1.0.4     2023-01-22 [2] CRAN (R 4.1.2)
##  farver                 2.1.1     2022-07-06 [2] CRAN (R 4.1.2)
##  fastmap                1.1.1     2023-02-24 [2] CRAN (R 4.1.2)
##  flashClust             1.01-2    2012-08-21 [2] CRAN (R 4.1.0)
##  foreach                1.5.2     2022-02-02 [2] CRAN (R 4.1.2)
##  foreign                0.8-84    2022-12-06 [2] CRAN (R 4.1.2)
##  Formula                1.2-5     2023-02-24 [2] CRAN (R 4.1.2)
##  fs                     1.6.2     2023-04-25 [2] CRAN (R 4.1.2)
##  genefilter             1.76.0    2021-10-26 [2] Bioconductor
##  geneplotter            1.72.0    2021-10-26 [2] Bioconductor
##  generics               0.1.3     2022-07-05 [2] CRAN (R 4.1.2)
##  GenomeInfoDb           1.30.1    2022-01-30 [2] Bioconductor
##  GenomeInfoDbData       1.2.7     2022-03-09 [2] Bioconductor
##  GenomicRanges          1.46.1    2021-11-18 [2] Bioconductor
##  ggplot2              * 3.4.2     2023-04-03 [2] CRAN (R 4.1.2)
##  ggpubr               * 0.6.0     2023-02-10 [2] CRAN (R 4.1.2)
##  ggrepel                0.9.3     2023-02-03 [2] CRAN (R 4.1.2)
##  ggsignif               0.6.4     2022-10-13 [2] CRAN (R 4.1.2)
##  glmnet                 4.1-7     2023-03-23 [2] CRAN (R 4.1.2)
##  glue                   1.6.2     2022-02-24 [2] CRAN (R 4.1.2)
##  gplots                 3.1.3     2022-04-25 [2] CRAN (R 4.1.2)
##  gridExtra              2.3       2017-09-09 [2] CRAN (R 4.1.0)
##  gtable                 0.3.3     2023-03-21 [2] CRAN (R 4.1.2)
##  gtools                 3.9.4     2022-11-27 [2] CRAN (R 4.1.2)
##  highr                  0.10      2022-12-22 [2] CRAN (R 4.1.2)
##  Hmisc                  5.1-0     2023-05-08 [2] CRAN (R 4.1.2)
##  htmlTable              2.4.1     2022-07-07 [2] CRAN (R 4.1.2)
##  htmltools              0.5.7     2023-11-03 [1] CRAN (R 4.1.3)
##  htmlwidgets            1.6.2     2023-03-17 [2] CRAN (R 4.1.2)
##  httpuv                 1.6.11    2023-05-11 [2] CRAN (R 4.1.3)
##  httr                   1.4.6     2023-05-08 [2] CRAN (R 4.1.2)
##  igraph                 1.5.0     2023-06-16 [1] CRAN (R 4.1.3)
##  IRanges                2.28.0    2021-10-26 [2] Bioconductor
##  iterators              1.0.14    2022-02-05 [2] CRAN (R 4.1.2)
##  jquerylib              0.1.4     2021-04-26 [2] CRAN (R 4.1.0)
##  jsonlite               1.8.7     2023-06-29 [2] CRAN (R 4.1.3)
##  kableExtra             1.3.4     2021-02-20 [2] CRAN (R 4.1.2)
##  KEGGREST               1.34.0    2021-10-26 [2] Bioconductor
##  KernSmooth             2.23-22   2023-07-10 [2] CRAN (R 4.1.3)
##  knitr                  1.43      2023-05-25 [2] CRAN (R 4.1.3)
##  labeling               0.4.2     2020-10-20 [2] CRAN (R 4.1.0)
##  later                  1.3.1     2023-05-02 [2] CRAN (R 4.1.2)
##  lattice                0.21-8    2023-04-05 [2] CRAN (R 4.1.2)
##  leaps                  3.1       2020-01-16 [2] CRAN (R 4.1.0)
##  libcoin                1.0-9     2021-09-27 [2] CRAN (R 4.1.0)
##  lifecycle              1.0.3     2022-10-07 [2] CRAN (R 4.1.2)
##  limma                  3.50.3    2022-04-07 [2] Bioconductor
##  locfit                 1.5-9.8   2023-06-11 [2] CRAN (R 4.1.3)
##  magrittr               2.0.3     2022-03-30 [2] CRAN (R 4.1.2)
##  MASS                   7.3-60    2023-05-04 [2] CRAN (R 4.1.2)
##  Matrix                 1.6-0     2023-07-08 [2] CRAN (R 4.1.3)
##  MatrixGenerics         1.6.0     2021-10-26 [2] Bioconductor
##  matrixStats            1.0.0     2023-06-02 [2] CRAN (R 4.1.3)
##  memoise                2.0.1     2021-11-26 [2] CRAN (R 4.1.0)
##  metagenomeSeq          1.36.0    2021-10-26 [2] Bioconductor
##  mgcv                   1.8-42    2023-03-02 [2] CRAN (R 4.1.2)
##  mime                   0.12      2021-09-28 [2] CRAN (R 4.1.0)
##  miniUI                 0.1.1.1   2018-05-18 [2] CRAN (R 4.1.0)
##  modeltools             0.2-23    2020-03-05 [2] CRAN (R 4.1.0)
##  multcomp               1.4-25    2023-06-20 [2] CRAN (R 4.1.3)
##  multcompView           0.1-9     2023-04-09 [2] CRAN (R 4.1.2)
##  multtest               2.50.0    2021-10-26 [2] Bioconductor
##  munsell                0.5.0     2018-06-12 [2] CRAN (R 4.1.0)
##  mvtnorm                1.2-2     2023-06-08 [2] CRAN (R 4.1.3)
##  nlme                   3.1-162   2023-01-31 [2] CRAN (R 4.1.2)
##  nnet                   7.3-19    2023-05-03 [2] CRAN (R 4.1.2)
##  permute                0.9-7     2022-01-27 [2] CRAN (R 4.1.2)
##  phyloseq             * 1.38.0    2021-10-26 [2] Bioconductor
##  pillar                 1.9.0     2023-03-22 [2] CRAN (R 4.1.2)
##  pkgbuild               1.4.2     2023-06-26 [2] CRAN (R 4.1.3)
##  pkgconfig              2.0.3     2019-09-22 [2] CRAN (R 4.1.0)
##  pkgload                1.3.2.1   2023-07-08 [2] CRAN (R 4.1.3)
##  plyr                   1.8.8     2022-11-11 [2] CRAN (R 4.1.2)
##  png                    0.1-8     2022-11-29 [2] CRAN (R 4.1.2)
##  prettyunits            1.1.1     2020-01-24 [2] CRAN (R 4.1.0)
##  processx               3.8.2     2023-06-30 [2] CRAN (R 4.1.3)
##  profvis                0.3.8     2023-05-02 [2] CRAN (R 4.1.2)
##  promises               1.2.0.1   2021-02-11 [2] CRAN (R 4.1.0)
##  ps                     1.7.5     2023-04-18 [2] CRAN (R 4.1.2)
##  purrr                  1.0.1     2023-01-10 [2] CRAN (R 4.1.2)
##  R6                     2.5.1     2021-08-19 [2] CRAN (R 4.1.0)
##  RColorBrewer           1.1-3     2022-04-03 [2] CRAN (R 4.1.2)
##  Rcpp                   1.0.11    2023-07-06 [1] CRAN (R 4.1.3)
##  RCurl                  1.98-1.12 2023-03-27 [2] CRAN (R 4.1.2)
##  remotes                2.4.2     2021-11-30 [2] CRAN (R 4.1.0)
##  reshape2               1.4.4     2020-04-09 [2] CRAN (R 4.1.0)
##  rhdf5                  2.38.1    2022-03-10 [2] Bioconductor
##  rhdf5filters           1.6.0     2021-10-26 [2] Bioconductor
##  Rhdf5lib               1.16.0    2021-10-26 [2] Bioconductor
##  rlang                  1.1.1     2023-04-28 [1] CRAN (R 4.1.2)
##  rmarkdown              2.23      2023-07-01 [2] CRAN (R 4.1.3)
##  rpart                  4.1.19    2022-10-21 [2] CRAN (R 4.1.2)
##  RSQLite                2.3.1     2023-04-03 [2] CRAN (R 4.1.2)
##  rstatix                0.7.2     2023-02-01 [2] CRAN (R 4.1.2)
##  rstudioapi             0.15.0    2023-07-07 [2] CRAN (R 4.1.3)
##  rvest                  1.0.3     2022-08-19 [2] CRAN (R 4.1.2)
##  S4Vectors              0.32.4    2022-03-29 [2] Bioconductor
##  sandwich               3.0-2     2022-06-15 [2] CRAN (R 4.1.2)
##  sass                   0.4.6     2023-05-03 [2] CRAN (R 4.1.2)
##  scales                 1.2.1     2022-08-20 [2] CRAN (R 4.1.2)
##  scatterplot3d          0.3-44    2023-05-05 [2] CRAN (R 4.1.2)
##  sessioninfo            1.2.2     2021-12-06 [2] CRAN (R 4.1.0)
##  shape                  1.4.6     2021-05-19 [2] CRAN (R 4.1.0)
##  shiny                  1.7.4.1   2023-07-06 [2] CRAN (R 4.1.3)
##  stringi                1.7.12    2023-01-11 [2] CRAN (R 4.1.2)
##  stringr                1.5.0     2022-12-02 [2] CRAN (R 4.1.2)
##  SummarizedExperiment   1.24.0    2021-10-26 [2] Bioconductor
##  survival               3.5-5     2023-03-12 [2] CRAN (R 4.1.2)
##  svglite                2.1.1     2023-01-10 [2] CRAN (R 4.1.2)
##  systemfonts            1.0.4     2022-02-11 [2] CRAN (R 4.1.2)
##  TH.data                1.1-2     2023-04-17 [2] CRAN (R 4.1.2)
##  tibble               * 3.2.1     2023-03-20 [2] CRAN (R 4.1.2)
##  tidyr                  1.3.0     2023-01-24 [2] CRAN (R 4.1.2)
##  tidyselect             1.2.0     2022-10-10 [2] CRAN (R 4.1.2)
##  urlchecker             1.0.1     2021-11-30 [2] CRAN (R 4.1.0)
##  usethis                2.2.2     2023-07-06 [2] CRAN (R 4.1.3)
##  utf8                   1.2.3     2023-01-31 [2] CRAN (R 4.1.2)
##  vctrs                  0.6.3     2023-06-14 [1] CRAN (R 4.1.3)
##  vegan                  2.6-4     2022-10-11 [2] CRAN (R 4.1.2)
##  viridisLite            0.4.2     2023-05-02 [2] CRAN (R 4.1.2)
##  webshot                0.5.5     2023-06-26 [2] CRAN (R 4.1.3)
##  withr                  2.5.0     2022-03-03 [2] CRAN (R 4.1.2)
##  Wrench                 1.12.0    2021-10-26 [2] Bioconductor
##  xfun                   0.40      2023-08-09 [1] CRAN (R 4.1.3)
##  XMAS2                * 2.2.0     2023-11-30 [1] local
##  XML                    3.99-0.14 2023-03-19 [2] CRAN (R 4.1.2)
##  xml2                   1.3.5     2023-07-06 [2] CRAN (R 4.1.3)
##  xtable                 1.8-4     2019-04-21 [2] CRAN (R 4.1.0)
##  XVector                0.34.0    2021-10-26 [2] Bioconductor
##  yaml                   2.3.7     2023-01-23 [2] CRAN (R 4.1.2)
##  zlibbioc               1.40.0    2021-10-26 [2] Bioconductor
##  zoo                    1.8-12    2023-04-13 [2] CRAN (R 4.1.2)
## 
##  [1] /Users/zouhua/Library/R/x86_64/4.1/library
##  [2] /Library/Frameworks/R.framework/Versions/4.1/Resources/library
## 
## ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────