Set Up

rm(list = ls())
set.seed(2024)
library(tidyverse)
library(here)
library(rstatix)
library(phyloseq)
library(gt)
library(broom)
theme_set(theme_bw())

ps.rare <- readRDS(here('data','following_study','ps_rarefied.rds'))

Collinsella

In DA analyses section, we found one Collinsella ASV that was consistently identified differentially abundant by all DA methods.

Differential Abundance

The log2 fold difference in relative abundance between IE and control dogs from the same household for the five most abundant Collinsella ASVs and total Collinsella. Absent bars indicate that the Collinsella ASV was not detected in either dog in that household.

ps.collinsella <- ps.rare %>%
    subset_taxa(Genus == 'Collinsella') 
genus.collinsella <- ps.collinsella %>% 
    tax_glom('Genus') %>% 
    otu_table() %>% data.frame() %>% 
    rename(Collinsella = ASV2)

# add genus level collinsella count to otu table
collinsella.df <- ps.collinsella %>% 
    filter_taxa(function(x) sum(x>0)/length(x)>=0.1, prune = TRUE) %>%  # filter out taxa with prevalence lower than 10%
    otu_table() %>% 
    cbind(genus.collinsella) %>%
    cbind(data.frame(sample_data(ps.collinsella))) %>% 
    pivot_longer(c(starts_with('ASV'), 'Collinsella'), names_to = 'ASV', values_to = 'Abundance')


lfc <- collinsella.df %>%
    group_by(Household, ASV) %>%
    summarise(LFC = log2(Abundance[Epileptic.or.Control == "Epileptic"] + 0.5) -
                  log2(Abundance[Epileptic.or.Control == "Control"] + 0.5),
              .groups = 'drop')

ggplot(lfc) +
    geom_bar(aes(x = as.numeric(Household), y = LFC, fill = (LFC > 0)), stat = 'identity') +
    facet_wrap(~ASV, scales = 'free_y') +
    labs(x = 'Household', y = 'log 2 fold change', fill = 'LFC > 0')

ggsave(here('figures','Figure5.png'))
## Saving 7 x 5 in image
ggsave(here('figures','Figure5.pdf'))
## Saving 7 x 5 in image

We can visually see some ASVs are more abunant in IE group. Here, we use pairwise sign test to check.

collinsella.df %>% 
    group_by(ASV) %>%
    arrange(Household) %>% # make sure data are paired
    pairwise_sign_test(Abundance ~ Epileptic.or.Control, ref.group = 'Epileptic', alternative = 'greater', paired = TRUE) %>% gt
ASV .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
ASV2 Abundance Epileptic Control 49 49 34 49 0.005 0.005 **
ASV44 Abundance Epileptic Control 49 49 15 20 0.021 0.021 *
ASV70 Abundance Epileptic Control 49 49 14 23 0.202 0.202 ns
ASV79 Abundance Epileptic Control 49 49 13 29 0.771 0.771 ns
ASV83 Abundance Epileptic Control 49 49 15 30 0.572 0.572 ns
Collinsella Abundance Epileptic Control 49 49 35 49 0.002 0.002 **

ASV2, ASV44 and the total Collinsella are significantly more abundant in IE dogs.

Chi-Square test

Association between sex and the presence of Collinsella can have a confounding effect

ftable(xtabs( ~ Sex  + sign(Abundance) + ASV, data = collinsella.df))
##                     ASV ASV2 ASV44 ASV70 ASV79 ASV83 Collinsella
## Sex sign(Abundance)                                             
## F   0                      1    45    34    25    29           1
##     1                     57    13    24    33    29          57
## M   0                      0    26    28    28    23           0
##     1                     40    14    12    12    17          40
collinsella.df %>% 
    split(.$ASV) %>% 
    map_dfr(function(df) tidy(chisq.test(xtabs( ~ Sex  + sign(Abundance), data = df))), .id = 'ASV') %>%
    gt %>% fmt_number(decimals = 3)
## Warning in chisq.test(xtabs(~Sex + sign(Abundance), data = df)): Chi-squared
## approximation may be incorrect
## Warning in chisq.test(xtabs(~Sex + sign(Abundance), data = df)): Chi-squared
## approximation may be incorrect
ASV statistic p.value parameter method
ASV2 0.000 1.000 1.000 Pearson's Chi-squared test with Yates' continuity correction
ASV44 1.301 0.254 1.000 Pearson's Chi-squared test with Yates' continuity correction
ASV70 0.875 0.350 1.000 Pearson's Chi-squared test with Yates' continuity correction
ASV79 5.856 0.016 1.000 Pearson's Chi-squared test with Yates' continuity correction
ASV83 0.276 0.599 1.000 Pearson's Chi-squared test with Yates' continuity correction
Collinsella 0.000 1.000 1.000 Pearson's Chi-squared test with Yates' continuity correction

Presence of ASV79 is associated with sex.

Abundance Change Across other factors

Age

Does the abundance of Collinsella different by age?

collinsella.df %>% 
    split(.$ASV) %>% 
    map_dfr(function(df) tidy(lm(Abundance ~ Age..months., data = df)), .id = 'ASV') %>% 
    group_by(ASV) %>% gt %>% fmt_number(decimals = 3)
term estimate std.error statistic p.value
ASV2
(Intercept) 3,893.154 591.636 6.580 0.000
Age..months. −0.771 7.078 −0.109 0.913
ASV44
(Intercept) 4.294 187.426 0.023 0.982
Age..months. 2.921 2.242 1.303 0.196
ASV70
(Intercept) −34.719 51.240 −0.678 0.500
Age..months. 1.400 0.613 2.284 0.025
ASV79
(Intercept) 19.143 60.824 0.315 0.754
Age..months. 0.748 0.728 1.028 0.307
ASV83
(Intercept) −0.162 19.133 −0.008 0.993
Age..months. 0.726 0.229 3.172 0.002
Collinsella
(Intercept) 3,927.815 586.870 6.693 0.000
Age..months. 4.794 7.021 0.683 0.496

The abundance of ASV83 was significantly higher in older dogs.

Seizure Freedom/Seizure Control

Does the abundance of Collinsella change by seizure freedom or seizure controlled dogs?

collinsella.df %>% 
    filter(Epileptic.or.Control == 'Epileptic') %>%
    split(.$ASV) %>%
    map_dfr(function(df) tidy(t.test(Abundance ~ Seizure.Freedom..Y.N., data = df)), .id = 'ASV') %>%
    gt %>% fmt_number(decimals = 3)
ASV estimate estimate1 estimate2 statistic p.value parameter conf.low conf.high method alternative
ASV2 −296.688 4,176.125 4,472.812 −0.390 0.699 31.314 −1,849.271 1,255.896 Welch Two Sample t-test two.sided
ASV44 −180.312 268.812 449.125 −0.517 0.611 22.078 −904.025 543.400 Welch Two Sample t-test two.sided
ASV70 69.656 114.031 44.375 1.118 0.271 36.866 −56.637 195.950 Welch Two Sample t-test two.sided
ASV79 −185.094 20.719 205.812 −1.385 0.186 15.093 −469.799 99.612 Welch Two Sample t-test two.sided
ASV83 −40.688 48.812 89.500 −1.210 0.240 20.213 −110.781 29.406 Welch Two Sample t-test two.sided
Collinsella −641.688 4,652.375 5,294.062 −0.832 0.413 27.717 −2,222.698 939.323 Welch Two Sample t-test two.sided
 
collinsella.df %>% 
    filter(Epileptic.or.Control == 'Epileptic') %>%
    split(.$ASV) %>%
    map_dfr(function(df) tidy(t.test(Abundance ~ Seizure.Control..Satisfactory.Unsatisfactory., data = df)), .id = 'ASV') %>%
    gt %>% fmt_number(decimals = 3)
ASV estimate estimate1 estimate2 statistic p.value parameter conf.low conf.high method alternative
ASV2 154.789 4,316.943 4,162.154 0.178 0.860 19.465 −1,658.623 1,968.201 Welch Two Sample t-test two.sided
ASV44 −149.389 288.457 437.846 −0.420 0.679 18.537 −894.329 595.551 Welch Two Sample t-test two.sided
ASV70 −92.354 65.800 158.154 −0.673 0.513 12.794 −389.207 204.500 Welch Two Sample t-test two.sided
ASV79 98.681 109.143 10.462 1.581 0.123 34.555 −28.061 225.424 Welch Two Sample t-test two.sided
ASV83 28.892 70.200 41.308 1.130 0.267 31.857 −23.203 80.987 Welch Two Sample t-test two.sided
Collinsella 71.791 4,885.714 4,813.923 0.087 0.932 20.149 −1,654.774 1,798.356 Welch Two Sample t-test two.sided

none of the Collinsella ASV were found to be significantly differentially abundant between seizure free or seizure controlled dogs.

Lactobacillus

Lactobacillus has attracted growing attention for its potential beneficial role in epilepsy.

Differential Abundance

ps.lactobacillus <- ps.rare %>%
    subset_taxa(Genus == 'Lactobacillus') 
genus.lactobacillus <- ps.lactobacillus %>% 
    tax_glom('Genus') %>% 
    otu_table() %>% data.frame() %>% 
    rename(Lactobacillus = ASV55)

# add genus level lactobacillus count to otu table
lactobacillus.df <- ps.lactobacillus %>% 
    filter_taxa(function(x) sum(x>0)/length(x)>=0.1, prune = TRUE) %>%  # filter out taxa with prevalence lower than 10%
    otu_table() %>% 
    cbind(genus.lactobacillus) %>%
    cbind(data.frame(sample_data(ps.lactobacillus))) %>% 
    pivot_longer(c(starts_with('ASV'), 'Lactobacillus'), names_to = 'ASV', values_to = 'Abundance')


lfc <- lactobacillus.df %>%
    group_by(Household, ASV) %>%
    summarise(LFC = log2(Abundance[Epileptic.or.Control == "Epileptic"] + 1) -
                  log2(Abundance[Epileptic.or.Control == "Control"] + 1),
              .groups = 'drop')

ggplot(lfc) +
    geom_bar(aes(x = as.numeric(Household), y = LFC, fill = (LFC > 0)), stat = 'identity') +
    facet_wrap(~ASV, scales = 'free_y') +
    labs(x = 'Household', y = 'log 2 fold change', fill = 'LFC > 0')

lactobacillus.df  %>% 
    group_by(ASV) %>%
    arrange(Household) %>% # make sure data are paired
    pairwise_sign_test(Abundance ~ Epileptic.or.Control, ref.group = 'Epileptic', alternative = 'greater', paired = TRUE) %>% gt
ASV .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
ASV105 Abundance Epileptic Control 49 49 4 8 0.637 0.637 ns
ASV55 Abundance Epileptic Control 49 49 9 13 0.133 0.133 ns
ASV906 Abundance Epileptic Control 49 49 7 8 0.035 0.035 *
Lactobacillus Abundance Epileptic Control 49 49 16 26 0.163 0.163 ns

ASV906 is significantly more abundant in epileptic dogs

Chi-Square Test

Does the presense of Lactobacillus associated with sex?

ftable(xtabs(sign(Abundance) ~ Sex + Epileptic.or.Control + ASV, data = lactobacillus.df))
##                          ASV ASV105 ASV55 ASV906 Lactobacillus
## Sex Epileptic.or.Control                                      
## F   Control                       5     4      4            13
##     Epileptic                     2     5      6            12
## M   Control                       2     2      0             5
##     Epileptic                     3     6      1             8
lactobacillus.df %>% 
    split(.$ASV) %>% 
    map_dfr(function(df) tidy(chisq.test(xtabs( ~ Sex  + sign(Abundance), data = df))), .id = 'ASV') %>%
    gt %>% fmt_number(decimals = 3)
## Warning in chisq.test(xtabs(~Sex + sign(Abundance), data = df)): Chi-squared
## approximation may be incorrect
## Warning in chisq.test(xtabs(~Sex + sign(Abundance), data = df)): Chi-squared
## approximation may be incorrect
ASV statistic p.value parameter method
ASV105 0.000 1.000 1.000 Pearson's Chi-squared test with Yates' continuity correction
ASV55 0.093 0.761 1.000 Pearson's Chi-squared test with Yates' continuity correction
ASV906 3.789 0.052 1.000 Pearson's Chi-squared test with Yates' continuity correction
Lactobacillus 0.719 0.396 1.000 Pearson's Chi-squared test with Yates' continuity correction

ASV906 is marginally significant

Abundance Change

Age

Does the abundance change by age?

lactobacillus.df %>% 
    split(.$ASV) %>% 
    map_dfr(function(df) tidy(lm(Abundance ~ Age..months., data = df)), .id = 'ASV') %>% 
    group_by(ASV) %>% gt %>% fmt_number(decimals = 3)
term estimate std.error statistic p.value
ASV105
(Intercept) −26.663 60.499 −0.441 0.660
Age..months. 0.905 0.724 1.250 0.214
ASV55
(Intercept) −329.842 173.680 −1.899 0.061
Age..months. 5.587 2.078 2.689 0.008
ASV906
(Intercept) 0.338 0.291 1.162 0.248
Age..months. 0.000 0.003 −0.043 0.966
Lactobacillus
(Intercept) −357.026 196.484 −1.817 0.072
Age..months. 6.574 2.351 2.797 0.006

ASV55 and the total Lactobacillus abundance is significantly more abundant in epileptic dogs.

Seizure Free / Seizure Control

Does abundance change by these groups?

lactobacillus.df %>% 
    filter(Epileptic.or.Control == 'Epileptic') %>%
    split(.$ASV) %>%
    map_dfr(function(df) tidy(t.test(Abundance ~ Seizure.Freedom..Y.N., data = df)), .id = 'ASV') %>%
    gt %>% fmt_number(decimals = 3)
ASV estimate estimate1 estimate2 statistic p.value parameter conf.low conf.high method alternative
ASV105 52.156 52.406 0.250 1.491 0.146 31.003 −19.166 123.479 Welch Two Sample t-test two.sided
ASV55 215.719 233.594 17.875 0.927 0.361 31.233 −258.556 689.994 Welch Two Sample t-test two.sided
ASV906 −1.062 0.188 1.250 −1.671 0.113 17.014 −2.404 0.279 Welch Two Sample t-test two.sided
Lactobacillus 266.438 286.375 19.938 1.083 0.287 31.212 −235.323 768.198 Welch Two Sample t-test two.sided
 
lactobacillus.df %>% 
    filter(Epileptic.or.Control == 'Epileptic') %>%
    split(.$ASV) %>%
    map_dfr(function(df) tidy(t.test(Abundance ~ Seizure.Control..Satisfactory.Unsatisfactory., data = df)), .id = 'ASV') %>%
    gt %>% fmt_number(decimals = 3)
ASV estimate estimate1 estimate2 statistic p.value parameter conf.low conf.high method alternative
ASV105 47.185 47.800 0.615 1.473 0.150 34.025 −17.921 112.291 Welch Two Sample t-test two.sided
ASV55 221.321 221.629 0.308 1.043 0.304 34.000 −209.805 652.447 Welch Two Sample t-test two.sided
ASV906 0.637 0.714 0.077 1.927 0.062 37.651 −0.032 1.307 Welch Two Sample t-test two.sided
Lactobacillus 269.044 270.429 1.385 1.198 0.239 34.001 −187.249 725.337 Welch Two Sample t-test two.sided

ASV906 was marginally less abundant in IE dogs that are seizure controlled

More Information

Prevalence

prev <- ps.rare %>%
    filter_taxa(function(x) sum(x>0)/length(x)>=0.1, prune = TRUE) %>% 
    psmelt() %>% 
    filter(Genus %in% c('Collinsella', 'Lactobacillus')) %>%
    group_by(OTU, Genus) %>%
    summarise(Frequency = sum(Abundance > 0),
              Prevalence = sum(Abundance > 0) / n()) %>% 
    arrange(Genus)
## `summarise()` has grouped output by 'OTU'. You can override using the `.groups`
## argument.
prev
## # A tibble: 8 × 4
## # Groups:   OTU [8]
##   OTU    Genus         Frequency Prevalence
##   <chr>  <chr>             <int>      <dbl>
## 1 ASV2   Collinsella          97      0.990
## 2 ASV44  Collinsella          27      0.276
## 3 ASV70  Collinsella          36      0.367
## 4 ASV79  Collinsella          45      0.459
## 5 ASV83  Collinsella          46      0.469
## 6 ASV105 Lactobacillus        12      0.122
## 7 ASV55  Lactobacillus        17      0.173
## 8 ASV906 Lactobacillus        11      0.112
ps.rare %>%
    tax_glom('Genus') %>% 
    psmelt() %>% 
    filter(Genus %in% c('Collinsella', 'Lactobacillus')) %>%
    group_by(OTU, Genus) %>%
    summarise(Frequency = sum(Abundance > 0),
              Prevalence = sum(Abundance > 0) / n()) %>% 
    arrange(Genus)
## `summarise()` has grouped output by 'OTU'. You can override using the `.groups`
## argument.
## # A tibble: 2 × 4
## # Groups:   OTU [2]
##   OTU   Genus         Frequency Prevalence
##   <chr> <chr>             <int>      <dbl>
## 1 ASV2  Collinsella          97      0.990
## 2 ASV55 Lactobacillus        38      0.388

Taxonomy

tax_table(ps.rare)[prev$OTU, c('Genus','Species')]
## Taxonomy Table:     [8 taxa by 2 taxonomic ranks]:
##        Genus           Species      
## ASV2   "Collinsella"   NA           
## ASV44  "Collinsella"   "aerofaciens"
## ASV70  "Collinsella"   "tanakaei"   
## ASV79  "Collinsella"   NA           
## ASV83  "Collinsella"   "phocaeensis"
## ASV105 "Lactobacillus" NA           
## ASV55  "Lactobacillus" NA           
## ASV906 "Lactobacillus" NA

Some of the species were not identified. We can use blast to search against the nt/nr database and get their percentage identity.

# extract sequence of these ASVs
refseq(ps.rare)[prev$OTU] %>% dada2:::pfasta(ids = prev$OTU)
## >ASV2
## TACGTAGGGGGCGAGCGTTATCCGGATTCATTGGGCGTAAAGCGCGCGTAGGCGGCCCGGCAGGCAGGGGGTCAAATGGCGGGGCTCAACCCCGTCCCGCCCCCTGAACCGCCGGGCTCGGGTCCGGTAGGGGAGGGTGGAACACCCGGTGTAGCGGTGGAATGCGCAGATATCGGGTGGAACACCGGTGGCGAAGGCGGCCCTCTGGGCCGAGACCGACGCTGAGGCGCGAAAGCTGGGGGAGCGAACAGG
## >ASV44
## TACGTAGGGGGCGAGCGTTATCCGGATTCATTGGGCGTAAAGCGCGCGTAGGCGGCCCGGCAGGCCGGGGGTCGAAGCGGGGGGCTCAACCCCCCGAAGCCCCCGGAACCTCCGCGGCTTGGGTCCGGTAGGGGAGGGTGGAACACCCGGTGTAGCGGTGGAATGCGCAGATATCGGGTGGAACACCGGTGGCGAAGGCGGCCCTCTGGGCCGAGACCGACGCTGAGGCGCGAAAGCTGGGGGAGCGAACAGG
## >ASV70
## TACGTAGGGGGCGAGCGTTATCCGGATTCATTGGGCGTAAAGCGCGCGTAGGCGGCCGCGTAGGCGGGGGGTCAAATCCCGGGGCTCAACCCCGGTCCGCCCCCCGAACCCCGCGGCTCGGGTCCGGTAGGGGAGGGTGGAATTCCCGGTGTAGCGGTGGAATGCGCAGATATCGGGAGGAACACCGGTGGCGAAGGCGGCCCTCTGGGCCGAGACCGACGCTGAGGCGCGAAAGCTGGGGGAGCGAACAGG
## >ASV79
## TACGTAGGGGGCGAGCGTTATCCGGATTCATTGGGCGTAAAGCGCGCGTAGGCGGCCGCGTAGGCGGGGGGTCAAATCCCGGGGCTCAACCCCGGTCCGCCCCCCGAACCCCGCGGCTCGGGTCCGGTAGGGGAGGGTGGAACACCCGGTGTAGCGGTGGAATGCGCAGATATCGGGTGGAACACCGGTGGCGAAGGCGGCCCTCTGGGCCGAGACCGACGCTGAGGCGCGAAAGCTGGGGGAGCGAACAGG
## >ASV83
## TACGTAGGGGGCGAGCGTTATCCGGATTCATTGGGCGTAAAGCGCGCGTAGGCGGCCCGGCAGGCAGGGGGTCAAATGGCGGGGCTCAACCCCGTCCCGCCCCCTGAACCGCCGGGCTTGGGTCCGGTAGGGGAGGGTGGAACACCCGGTGTAGCGGTGGAATGCGCAGATATCGGGTGGAACACCGGTGGCGAAGGCGGCCCTCTGGGCCGAGACCGACGCTGAGGCGCGAAAGCTGGGGGAGCGAACAGG
## >ASV105
## TACGTAGGTGGCAAGCGTTGTCCGGATTTATTGGGCGTAAAGCGAGTGCAGGCGGTTCAATAAGTCTGATGTGAAAGCCTTCGGCTCAACCGGAGAATTGCATCAGAAACTGTTGAACTTGAGTGCAGAAGAGGAGAGTGGAACTCCATGTGTAGCGGTGGAATGCGTAGATATATGGAAGAACACCAGTGGCGAAGGCGGCTCTCTGGTCTGCAACTGACGCTGAGGCTCGAAAGCATGGGTAGCGAACAGG
## >ASV55
## TACGTAGGTGGCAAGCGTTGTCCGGATTTATTGGGCGTAAAGCGAGCGCAGGCGGAAGAATAAGTCTGATGTGAAAGCCCTCGGCTTAACCGAGGAACTGCATCGGAAACTGTTTTTCTTGAGTGCAGAAGAGGAGAGTGGAACTCCATGTGTAGCGGTGGAATGCGTAGATATATGGAAGAACACCAGTGGCGAAGGCGGCTCTCTGGTCTGCAACTGACGCTGAGGCTCGAAAGCATGGGTAGCGAACAGG
## >ASV906
## TACGTAGGTGGCAAGCGTTGTCCGGATTTATTGGGCGTAAAGCGAGCGCAGGCGGAAAAATAAGTCTAATGTGAAAGCCCTCGGCTTAACCGAGGAACTGCATCGGAAACTGTTTTTCTTGAGTGCAGAAGAGGAGAGTGGAACTCCATGTGTAGCGGTGGAATGCGTAGATATATGGAAGAACACCAGTGGCGAAGGCGGCTCTCTGGTCTGCAACTGACGCTGAGGCTCGAAAGCATGGGTAGCGAACAGG

Here’s the result:

taxon species identity
ASV2 Collinsella stercoris 100%
ASV44 Collinsella aerofaciens 100%
ASV70 Collinsella tanakaei 100%
ASV79 Collinsella tanakaei 98.81%
ASV83 Collinsella phocaeensis 100%
ASV105 Lactobacillus gasseri 100%
ASV55 Lactobacillus crispatus 100%
ASV906 Lactobacillus amylovorus 100% 
sessioninfo::session_info()
## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value
##  version  R version 4.4.1 (2024-06-14)
##  os       macOS 15.5
##  system   aarch64, darwin20
##  ui       X11
##  language (EN)
##  collate  en_US.UTF-8
##  ctype    en_US.UTF-8
##  tz       America/New_York
##  date     2025-06-17
##  pandoc   3.5 @ /Users/yixuanyang/miniforge3/bin/ (via rmarkdown)
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  package              * version date (UTC) lib source
##  abind                  1.4-8   2024-09-12 [2] CRAN (R 4.4.1)
##  ade4                   1.7-22  2023-02-06 [2] CRAN (R 4.4.0)
##  ape                    5.8-1   2024-12-16 [2] CRAN (R 4.4.1)
##  backports              1.5.0   2024-05-23 [2] CRAN (R 4.4.0)
##  Biobase                2.66.0  2024-10-29 [2] Bioconductor 3.20 (R 4.4.1)
##  BiocGenerics           0.52.0  2024-10-29 [2] Bioconductor 3.20 (R 4.4.1)
##  BiocParallel           1.39.0  2024-05-04 [2] Bioconductor 3.20 (R 4.4.0)
##  biomformat             1.34.0  2024-10-29 [2] Bioconductor 3.20 (R 4.4.1)
##  Biostrings             2.74.0  2024-10-29 [2] Bioconductor 3.20 (R 4.4.1)
##  bitops                 1.0-9   2024-10-03 [2] CRAN (R 4.4.1)
##  broom                * 1.0.7   2024-09-26 [2] CRAN (R 4.4.1)
##  bslib                  0.9.0   2025-01-30 [1] CRAN (R 4.4.1)
##  cachem                 1.1.0   2024-05-16 [2] CRAN (R 4.4.0)
##  car                    3.1-3   2024-09-27 [2] CRAN (R 4.4.1)
##  carData                3.0-5   2022-01-06 [2] CRAN (R 4.4.0)
##  cli                    3.6.4   2025-02-13 [2] CRAN (R 4.4.1)
##  cluster                2.1.8   2024-12-11 [2] CRAN (R 4.4.1)
##  codetools              0.2-20  2024-03-31 [2] CRAN (R 4.4.1)
##  colorspace             2.1-1   2024-07-26 [2] CRAN (R 4.4.0)
##  crayon                 1.5.3   2024-06-20 [2] CRAN (R 4.4.0)
##  dada2                  1.33.0  2024-05-04 [2] Bioconductor 3.20 (R 4.4.0)
##  data.table             1.16.4  2024-12-06 [2] CRAN (R 4.4.1)
##  DelayedArray           0.32.0  2024-10-29 [2] Bioconductor 3.20 (R 4.4.1)
##  deldir                 2.0-4   2024-02-28 [2] CRAN (R 4.4.0)
##  digest                 0.6.37  2024-08-19 [2] CRAN (R 4.4.1)
##  dplyr                * 1.1.4   2023-11-17 [2] CRAN (R 4.4.0)
##  evaluate               1.0.3   2025-01-10 [1] CRAN (R 4.4.1)
##  farver                 2.1.2   2024-05-13 [2] CRAN (R 4.4.0)
##  fastmap                1.2.0   2024-05-15 [2] CRAN (R 4.4.0)
##  forcats              * 1.0.0   2023-01-29 [2] CRAN (R 4.4.0)
##  foreach                1.5.2   2022-02-02 [2] CRAN (R 4.4.0)
##  Formula                1.2-5   2023-02-24 [2] CRAN (R 4.4.0)
##  generics               0.1.3   2022-07-05 [2] CRAN (R 4.4.0)
##  GenomeInfoDb           1.42.0  2024-10-29 [2] Bioconductor 3.20 (R 4.4.1)
##  GenomeInfoDbData       1.2.13  2024-10-05 [2] Bioconductor
##  GenomicAlignments      1.42.0  2024-10-29 [2] Bioconductor 3.20 (R 4.4.1)
##  GenomicRanges          1.58.0  2024-10-29 [2] Bioconductor 3.20 (R 4.4.1)
##  ggplot2              * 3.5.1   2024-04-23 [2] CRAN (R 4.4.0)
##  glue                   1.8.0   2024-09-30 [2] CRAN (R 4.4.1)
##  gt                   * 0.11.1  2024-10-04 [2] CRAN (R 4.4.1)
##  gtable                 0.3.6   2024-10-25 [2] CRAN (R 4.4.1)
##  here                 * 1.0.1   2020-12-13 [2] CRAN (R 4.4.0)
##  hms                    1.1.3   2023-03-21 [2] CRAN (R 4.4.0)
##  htmltools              0.5.8.1 2024-04-04 [2] CRAN (R 4.4.0)
##  httr                   1.4.7   2023-08-15 [2] CRAN (R 4.4.0)
##  hwriter                1.3.2.1 2022-04-08 [2] CRAN (R 4.4.0)
##  igraph                 2.1.2   2024-12-07 [2] CRAN (R 4.4.1)
##  interp                 1.1-6   2024-01-26 [2] CRAN (R 4.4.0)
##  IRanges                2.40.0  2024-10-29 [2] Bioconductor 3.20 (R 4.4.1)
##  iterators              1.0.14  2022-02-05 [2] CRAN (R 4.4.0)
##  jpeg                   0.1-10  2022-11-29 [2] CRAN (R 4.4.0)
##  jquerylib              0.1.4   2021-04-26 [2] CRAN (R 4.4.0)
##  jsonlite               2.0.0   2025-03-27 [1] CRAN (R 4.4.1)
##  knitr                  1.50    2025-03-16 [1] CRAN (R 4.4.1)
##  labeling               0.4.3   2023-08-29 [2] CRAN (R 4.4.0)
##  lattice                0.22-6  2024-03-20 [2] CRAN (R 4.4.1)
##  latticeExtra           0.6-30  2022-07-04 [2] CRAN (R 4.4.0)
##  lifecycle              1.0.4   2023-11-07 [2] CRAN (R 4.4.0)
##  lubridate            * 1.9.4   2024-12-08 [2] CRAN (R 4.4.1)
##  magrittr               2.0.3   2022-03-30 [2] CRAN (R 4.4.0)
##  MASS                   7.3-61  2024-06-13 [2] CRAN (R 4.4.0)
##  Matrix                 1.7-1   2024-10-18 [2] CRAN (R 4.4.1)
##  MatrixGenerics         1.18.0  2024-10-29 [2] Bioconductor 3.20 (R 4.4.1)
##  matrixStats            1.4.1   2024-09-08 [2] CRAN (R 4.4.1)
##  mgcv                   1.9-1   2023-12-21 [2] CRAN (R 4.4.1)
##  multtest               2.62.0  2024-10-29 [2] Bioconductor 3.20 (R 4.4.1)
##  munsell                0.5.1   2024-04-01 [2] CRAN (R 4.4.0)
##  nlme                   3.1-166 2024-08-14 [2] CRAN (R 4.4.0)
##  permute                0.9-7   2022-01-27 [2] CRAN (R 4.4.0)
##  phyloseq             * 1.50.0  2024-10-29 [2] Bioconductor 3.20 (R 4.4.1)
##  pillar                 1.10.2  2025-04-05 [1] CRAN (R 4.4.1)
##  pkgconfig              2.0.3   2019-09-22 [2] CRAN (R 4.4.0)
##  plyr                   1.8.9   2023-10-02 [2] CRAN (R 4.4.0)
##  png                    0.1-8   2022-11-29 [2] CRAN (R 4.4.0)
##  purrr                * 1.0.4   2025-02-05 [2] CRAN (R 4.4.1)
##  pwalign                1.2.0   2024-10-29 [2] Bioconductor 3.20 (R 4.4.1)
##  R6                     2.6.1   2025-02-15 [1] CRAN (R 4.4.1)
##  ragg                   1.3.3   2024-09-11 [2] CRAN (R 4.4.1)
##  RColorBrewer           1.1-3   2022-04-03 [2] CRAN (R 4.4.0)
##  Rcpp                   1.0.14  2025-01-12 [1] CRAN (R 4.4.1)
##  RcppParallel           5.1.9   2024-08-19 [2] CRAN (R 4.4.1)
##  readr                * 2.1.5   2024-01-10 [2] CRAN (R 4.4.0)
##  reshape2               1.4.4   2020-04-09 [2] CRAN (R 4.4.0)
##  rhdf5                  2.49.0  2024-05-04 [2] Bioconductor 3.20 (R 4.4.0)
##  rhdf5filters           1.17.0  2024-05-04 [2] Bioconductor 3.20 (R 4.4.0)
##  Rhdf5lib               1.27.0  2024-05-04 [2] Bioconductor 3.20 (R 4.4.0)
##  rlang                  1.1.5   2025-01-17 [2] CRAN (R 4.4.1)
##  rmarkdown              2.29    2024-11-04 [2] CRAN (R 4.4.1)
##  rprojroot              2.0.4   2023-11-05 [2] CRAN (R 4.4.0)
##  Rsamtools              2.21.2  2024-09-26 [2] Bioconductor 3.20 (R 4.4.1)
##  rstatix              * 0.7.2   2023-02-01 [2] CRAN (R 4.4.0)
##  rstudioapi             0.17.1  2024-10-22 [2] CRAN (R 4.4.1)
##  S4Arrays               1.6.0   2024-10-29 [2] Bioconductor 3.20 (R 4.4.1)
##  S4Vectors              0.44.0  2024-10-29 [2] Bioconductor 3.20 (R 4.4.1)
##  sass                   0.4.9   2024-03-15 [2] CRAN (R 4.4.0)
##  scales                 1.3.0   2023-11-28 [2] CRAN (R 4.4.0)
##  sessioninfo            1.2.2   2021-12-06 [2] CRAN (R 4.4.0)
##  ShortRead              1.63.2  2024-09-26 [2] Bioconductor 3.20 (R 4.4.1)
##  SparseArray            1.6.0   2024-10-29 [2] Bioconductor 3.20 (R 4.4.1)
##  stringi                1.8.7   2025-03-27 [1] CRAN (R 4.4.1)
##  stringr              * 1.5.1   2023-11-14 [2] CRAN (R 4.4.0)
##  SummarizedExperiment   1.36.0  2024-10-29 [2] Bioconductor 3.20 (R 4.4.1)
##  survival               3.8-3   2024-12-17 [2] CRAN (R 4.4.1)
##  systemfonts            1.1.0   2024-05-15 [2] CRAN (R 4.4.0)
##  textshaping            0.4.1   2024-12-06 [2] CRAN (R 4.4.1)
##  tibble               * 3.2.1   2023-03-20 [2] CRAN (R 4.4.0)
##  tidyr                * 1.3.1   2024-01-24 [2] CRAN (R 4.4.0)
##  tidyselect             1.2.1   2024-03-11 [2] CRAN (R 4.4.0)
##  tidyverse            * 2.0.0   2023-02-22 [2] CRAN (R 4.4.0)
##  timechange             0.3.0   2024-01-18 [2] CRAN (R 4.4.0)
##  tzdb                   0.4.0   2023-05-12 [2] CRAN (R 4.4.0)
##  UCSC.utils             1.2.0   2024-10-29 [2] Bioconductor 3.20 (R 4.4.1)
##  utf8                   1.2.4   2023-10-22 [2] CRAN (R 4.4.0)
##  vctrs                  0.6.5   2023-12-01 [2] CRAN (R 4.4.0)
##  vegan                  2.6-8   2024-08-28 [2] CRAN (R 4.4.1)
##  withr                  3.0.2   2024-10-28 [2] CRAN (R 4.4.1)
##  xfun                   0.52    2025-04-02 [1] CRAN (R 4.4.1)
##  xml2                   1.3.6   2023-12-04 [2] CRAN (R 4.4.0)
##  XVector                0.46.0  2024-10-29 [2] Bioconductor 3.20 (R 4.4.1)
##  yaml                   2.3.10  2024-07-26 [2] CRAN (R 4.4.0)
##  zlibbioc               1.52.0  2024-10-29 [2] Bioconductor 3.20 (R 4.4.1)
## 
##  [1] /Users/yixuanyang/Library/R/arm64/4.4/library
##  [2] /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library
## 
## ──────────────────────────────────────────────────────────────────────────────