Differential Abundance Analysis
Yixuan Yang
Set Up
rm(list = ls())
set.seed(2024)
library(tidyverse)
library(phyloseq)
library(here)
library(gt)
library(GGally)
library(ggcorrplot)
theme_set(theme_bw())
# DA tools
library(ALDEx2)
library(ANCOMBC)
library(corncob)
library(MicrobiomeStat) # for Linda
ps <- readRDS(here('data','following_study','ps.rds'))
ps.rare <- readRDS(here('data','following_study','ps_rarefied.rds'))
ps.present <- ps %>% transform_sample_counts(sign)
alpha <- 0.05
max.cores <- parallel::detectCores()Taxa Filtering
Before the DA analysis, we want to filter the data by prevlence and minimum read depth. Here we are showing a tracing plot of the remained read and taxa for different prevlence threshold with minimum read depth of 50.
# test for filtering threshold
test.filter <- function(ps, prevlence, min_read = 0) {
sub.ps <- ps %>%
filter_taxa(function(x) sum(x)>min_read & sum(x>0)/length(x) >= prevlence, prune = TRUE)
retained.read <- otu_table(sub.ps) %>% rowSums() %>% mean()
retained.taxa <- otu_table(sub.ps) %>% ncol()
orig.read <- otu_table(ps) %>% rowSums() %>% mean()
orig.taxa <- otu_table(ps) %>% ncol()
return(data.frame(taxon_size = retained.taxa/orig.taxa, average.read.depth = retained.read/orig.read))
}
trace <- map_dfr(seq(0,0.9, 0.1), function(i) test.filter(ps=ps, prevlence=i))
ggplot(trace, aes(x = average.read.depth, y = taxon_size)) +
geom_point() +
geom_line() +
scale_x_reverse(labels = scales::percent) +
scale_y_continuous(labels = scales::percent) +
xlab('percentage of reads remained') +
ylab('percentage of ASV remained') 
By setting the minimum prevalence to 0.1 and minimum read depth of 100, we kept 13% ASVs with 96% reads remained.
head(trace)## taxon_size average.read.depth
## 1 1.00000000 1.0000000
## 2 0.13442167 0.9576704
## 3 0.07502605 0.9280344
## 4 0.05071205 0.8980429
## 5 0.03577631 0.8594170
## 6 0.02952414 0.8426915keep <- ps %>%
filter_taxa(function(x) sum(x>0)/length(x)>0.1, prune = TRUE) %>%
tax_table() %>% rownames()Basic information of non-rarefied data
ps <- ps %>% prune_taxa(keep, .)
ps## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 387 taxa and 98 samples ]
## sample_data() Sample Data: [ 98 samples by 28 sample variables ]
## tax_table() Taxonomy Table: [ 387 taxa by 7 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 387 tips and 386 internal nodes ]
## refseq() DNAStringSet: [ 387 reference sequences ]Basic information of rarefied data
ps.rare <- ps.rare %>% prune_taxa(keep, .)
ps.rare## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 387 taxa and 98 samples ]
## sample_data() Sample Data: [ 98 samples by 28 sample variables ]
## tax_table() Taxonomy Table: [ 387 taxa by 7 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 387 tips and 386 internal nodes ]
## refseq() DNAStringSet: [ 387 reference sequences ]Differential Abundance Analysis
ALDEx2
pairwise T-test
clr.data <- aldex.clr(t(otu_table(ps)), as.character(sample_data(ps)$Epileptic.or.Control), mc.samples=128, denom="all")
aldex.t.test.p <- aldex.ttest(clr.data, paired.test = TRUE)
aldex.t.test.effect <- aldex.effect(clr.data, paired.test = TRUE, useMC = TRUE, verbose = FALSE)
if (!identical(rownames(aldex.t.test.p), rownames(aldex.t.test.effect))) stop('taxa names are not matched')
aldex.t.test <- data.frame(aldex.t.test.p, aldex.t.test.effect) %>%
rownames_to_column('taxon') %>%
mutate(LFC = rab.win.Epileptic - rab.win.Control) %>%
mutate(p.value = we.ep) %>%
dplyr::select(taxon, LFC, p.value)
aldex.t.test %>% arrange(str_rank(taxon, numeric = TRUE)) %>% head()## taxon LFC p.value
## 1 ASV1 -0.05952564 0.686929409
## 2 ASV2 0.28095678 0.007961295
## 3 ASV3 0.17524459 0.848390423
## 4 ASV4 -0.04808613 0.538993838
## 5 ASV5 -1.14280561 0.276980164
## 6 ASV6 0.65206133 0.032268942GLM test
model <- model.matrix(~ Epileptic.or.Control + Household, data.frame(sample_data(ps)))
clr.data <- aldex.clr(t(otu_table(ps)), model, mc.samples=128, denom="all")
aldex.glm.test.p <- aldex.glm(clr.data, model)## |------------(25%)----------(50%)----------(75%)----------|aldex.glm.test.effect <- aldex.glm.effect(clr.data, verbose = FALSE, useMC = TRUE)$Epileptic.or.ControlEpileptic
if (!identical(rownames(aldex.glm.test.p), rownames(aldex.glm.test.effect))) stop('taxa names are not matched')
aldex.glm.test <- data.frame(p.value = aldex.glm.test.p$`Epileptic.or.ControlEpileptic:pval`, aldex.glm.test.effect) %>%
rownames_to_column('taxon') %>%
mutate(LFC = rab.win.1 - rab.win.0) %>% # 0 is control, 1 is epileptic, set is `data-cleaning.Rmd`
dplyr::select(taxon, LFC, p.value)
aldex.glm.test %>% arrange(str_rank(taxon, numeric = TRUE)) %>% head()## taxon LFC p.value
## 1 ASV1 -0.03414324 0.681977030
## 2 ASV2 0.24792837 0.007521401
## 3 ASV3 0.18789559 0.859989678
## 4 ASV4 -0.03754511 0.542065670
## 5 ASV5 -1.14198778 0.241079468
## 6 ASV6 0.66082144 0.030846682ANCOMBC
ancombc.test <- ancombc(data = ps,
formula = 'Epileptic.or.Control + Household',
taxa_are_rows = FALSE,
prv_cut = 0,
conserve = TRUE,
n_cl = max.cores)
ancombc.test <- data.frame(taxon = ancombc.test$res$lfc$taxon,
# convert logfold change to log2fold change
LFC = log2(exp(ancombc.test$res$lfc$Epileptic.or.ControlEpileptic)),
p.value = ancombc.test$res$p_val$Epileptic.or.ControlEpileptic)
ancombc.test %>% arrange(str_rank(taxon, numeric = TRUE)) %>% head## taxon LFC p.value
## 1 ASV1 -0.02483695 0.8989784944
## 2 ASV2 0.76488405 0.0005426785
## 3 ASV3 -0.14919417 0.5991345418
## 4 ASV4 -0.35090932 0.1478109441
## 5 ASV5 -0.80988944 0.0177297871
## 6 ASV6 0.44661012 0.0156575090corncob
corncob.test <- differentialTest(formula= ~ Epileptic.or.Control + Household,
phi.formula = ~ 1,
formula_null = ~ Household,
phi.formula_null = ~ 1,
test = 'LRT',
fdr_cutoff = 1, # don't filter by FDR
data = ps,
taxa_are_rows = FALSE)
# convert corcob result to log2fold difference
corncob.extract.result <- function(corncob.da, i) {
regression.result <- corncob.da$all_models[[i]]$coefficients
# p value
p <- corncob.da$p[i]
# relative abundance
RA.ctl <- regression.result['mu.(Intercept)','Estimate'] %>% invlogit()
RA.epi <- regression.result[c('mu.(Intercept)','mu.Epileptic.or.ControlEpileptic'),'Estimate'] %>% sum() %>% invlogit()
LFC <- RA.epi - RA.ctl # NOTE: this is not log2fold change, but the difference of relative abundance, since corncob does not provide lfc.
summary <- data.frame(LFC = LFC, p.value = p) %>% rownames_to_column('taxon')
return(summary)
}
corncob.test <- map_dfr(1:length(taxa_names(ps)), function(i) corncob.extract.result(corncob.test, i))
corncob.test %>%
arrange(str_rank(taxon, numeric = TRUE)) %>%
mutate(`relative abundance difference`= LFC) %>% dplyr::select(-LFC) %>% head## taxon p.value relative abundance difference
## 1 ASV1 0.7569841965 0.006657252
## 2 ASV2 0.0002559817 0.019875019
## 3 ASV3 0.6072699142 0.005305560
## 4 ASV4 0.4849536685 0.002897453
## 5 ASV5 0.0101062035 -0.008041232
## 6 ASV6 0.0087398620 0.010164788LinDA
linda.test <- linda(t(data.frame(otu_table(ps))),
data.frame(sample_data(ps)),
formula = '~ Epileptic.or.Control + Household',
feature.dat.type = 'count',
n.cores = max.cores,
verbose = FALSE)
linda.test <- linda.test$output$Epileptic.or.ControlEpileptic %>%
dplyr::select(log2FoldChange, pvalue) %>%
rownames_to_column('taxon')
colnames(linda.test)[2:3] <- c('LFC', 'p.value')
linda.test %>% arrange(str_rank(taxon, numeric = TRUE)) %>% head()## taxon LFC p.value
## 1 ASV1 0.008832133 0.97362763
## 2 ASV2 0.828567789 0.01110301
## 3 ASV3 -0.113867747 0.77997430
## 4 ASV4 -0.306886536 0.36471738
## 5 ASV5 -0.727514458 0.12941569
## 6 ASV6 0.505157615 0.04603249pairwise T-test on rarefied data
paired.t.test <- function(data, pseudo.count = 0.5) {
splited.data <- data %>%
arrange(Household, Epileptic.or.Control) %>%
split(.$Epileptic.or.Control)
t.test.result <- t.test(splited.data[['Epileptic']]$Abundance,
splited.data[['Control']]$Abundance,
paired = TRUE, alternative = 'two.sided')
data.frame(taxon = unique(data$OTU),
LFC = mean(log2(splited.data[['Epileptic']]$Abundance + pseudo.count) - log2(splited.data[['Control']]$Abundance + pseudo.count)),
p.value = t.test.result$p.value)
}
pairwise.t.test <- psmelt(ps.rare) %>%
group_by(OTU) %>%
group_split() %>%
map_dfr(paired.t.test)
pairwise.t.test %>% arrange(str_rank(taxon, numeric = TRUE)) %>% head()## taxon LFC p.value
## 1 ASV1 0.08532570 0.88527495
## 2 ASV2 0.88582424 0.01093835
## 3 ASV3 -0.03860217 0.93203722
## 4 ASV4 -0.27082143 0.76325632
## 5 ASV5 -0.66961126 0.03099468
## 6 ASV6 0.53771013 0.17916855Comparative Analysis
Here we compare the results of different methods. We first show the top 10 ASVs with the smallest p-value for each method.
test.results <- list("ALDEx2 t-test" = aldex.t.test,
"ALDEx2 GLM" = aldex.glm.test,
"ANCOM-BC" = ancombc.test,
"corncob" = corncob.test,
"LinDA" = linda.test,
"paired t-test" = pairwise.t.test) %>%
bind_rows(.id = 'method') %>%
group_by(method) %>% arrange(str_rank(taxon, numeric = TRUE)) %>% ungroup()
tax <- tax_table(ps) %>%
as.data.frame() %>%
dplyr::select(Phylum, Genus) %>%
rownames_to_column('taxon')
p.value.tb <- test.results %>%
dplyr::select(taxon, method, p.value) %>%
pivot_wider(names_from = method, values_from = p.value) %>%
rowwise() %>%
mutate(SignificantCount = sum(c_across(-taxon) < alpha))
top.p <- p.value.tb %>%
left_join(tax,., by = 'taxon') %>%
group_by(SignificantCount) %>%
arrange(str_rank(taxon, numeric = TRUE), .by_group = TRUE) %>%
arrange(desc(SignificantCount)) %>%
ungroup() %>%
dplyr::select(-SignificantCount)
tb <- top.p %>%
head(10) %>%
gt() %>%
fmt_number(decimals = 4) %>%
tab_style_body(
columns = where(is.numeric),
style = cell_text(weight = "bold"),
fn = function(x) x <= alpha
)
tb| taxon | Phylum | Genus | ALDEx2 t-test | ALDEx2 GLM | ANCOM-BC | corncob | LinDA | paired t-test |
|---|---|---|---|---|---|---|---|---|
| ASV2 | Actinobacteriota | Collinsella | 0.0080 | 0.0075 | 0.0005 | 0.0003 | 0.0111 | 0.0109 |
| ASV6 | Firmicutes | [Ruminococcus] gnavus group | 0.0323 | 0.0308 | 0.0157 | 0.0087 | 0.0460 | 0.1792 |
| ASV42 | Firmicutes | [Ruminococcus] torques group | 0.0892 | 0.0830 | 0.0001 | 0.0403 | 0.0158 | 0.0398 |
| ASV100 | Firmicutes | Lachnoclostridium | 0.1795 | 0.1688 | 0.0001 | 0.0001 | 0.0151 | 0.0111 |
| ASV420 | Firmicutes | Clostridium sensu stricto 1 | 0.1596 | 0.1502 | 0.0010 | 0.0001 | 0.0183 | 0.0288 |
| ASV446 | NA | NA | 0.1204 | 0.1121 | 0.0008 | 0.0002 | 0.0124 | 0.0340 |
| ASV645 | Firmicutes | NA | 0.3842 | 0.4280 | 0.0004 | 0.0002 | 0.0158 | 0.0292 |
| ASV702 | Firmicutes | Holdemania | 0.5177 | 0.4372 | 0.0000 | 0.0000 | 0.0086 | 0.0209 |
| ASV5 | Fusobacteriota | Fusobacterium | 0.2770 | 0.2411 | 0.0177 | 0.0101 | 0.1294 | 0.0310 |
| ASV12 | Firmicutes | Blautia | 0.1156 | 0.1222 | 0.0003 | 0.0355 | 0.0172 | 0.1086 |
There are 2 ASVs that are not classified to genus level by the naive bayes classifier. Here we search it with the BLAST nt/nr database
refseq(ps)[c('ASV446', 'ASV645')] %>% dada2:::pfasta(ids = c('ASV446', 'ASV645'))## >ASV446
## TCCTGTTTGCTCCCCACGCTTTCGAGCCTCAACGTCAGTCATCGTCCAGAAAGCCGCCTTCGCCACTGGTGTTCCTCCTAATATCTACGCATTTCACCGCTACACTAGGAATTCCGCTTTCCTCTCCGACACTCTAGCCTGACAGTTCCAAATGCAGTCCCGGGGTTGAGCCCCGGGCTTTCACATCTGGCTTGCCATGCCGTCTACGCTCCCTTTACACCCAGTAAATCCGGATAACGCTTGCCCCCTACGT
## >ASV645
## TACGTAGGGGGCAAGCGTTATCCGGATTTACTGGGTGTAAAGGGAGCGTAGACGGCGCAGCAAGTCTGAAGTGAAATCCCCGGGCTCAACCCGGGAACTGCTTTGGAAACTGTTGGGCTAGAGTGCTGGAGAGGCAAGCGGAATTCCTAGTGTAGCGGTGAAATGCGTAGATATTAGGAGGAACACCAGTGGCGAAGGCGGCTTGCTGGACAGTAACTGACGTTGAGGCTCGAAAGCGTGGGGAGCAAACAGGFor ASV446, the top hit is Firmicutes, Mediterraneibacter (Phylum, Genus)
For ASV645, the top hit is Firmicutes, Lachnospiraceae (Phylum, Family)
Does the effect size of each method agreed in direction? The plot
below shows that including corncob all methods have
consistent direction of effect size.
test.results$taxon <- factor(test.results$taxon,
levels = top.p$taxon)
test.results %>%
mutate(direction = if_else(LFC > 0, 'positive', 'negative')) %>%
filter(taxon %in% head(top.p$taxon, 10)) %>%
ggplot() +
geom_tile(aes(x = method, y = taxon, fill = direction), color = 'black') +
theme(axis.text.x = element_text(angle = 60, hjust = 1))
Effect size of top 10 ASVs for each method.
effect.plot <- test.results %>%
filter(method != 'corncob') %>%
filter(taxon %in% head(top.p$taxon, 10)) %>%
ggplot() +
geom_bar(aes(x = taxon, y = LFC, fill = method),
position="dodge", stat = 'identity', width = 0.7) +
theme(axis.text.x = element_text(angle = 60, hjust = 1))
effect.plot
correlation of p-value for each method.
top.p %>%
dplyr::select(-c(Phylum, Genus)) %>%
column_to_rownames('taxon') %>%
ggpairs(progress = FALSE)
cor.plot <- top.p %>%
dplyr::select(-c(Phylum, Genus)) %>%
column_to_rownames('taxon') %>%
cor() %>%
ggcorrplot(lab = TRUE,
lab_size = 3,
type = 'lower',
show.diag = TRUE,
show.legend = FALSE,
tl.cex = 10)
cor.plot
Volcano plots
test.results %>%
filter(method != 'corncob') %>%
ggplot() +
geom_point(aes(x = LFC, y = -log10(p.value))) +
facet_wrap(~method, scales = 'free')
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)
## ALDEx2 * 1.38.0 2024-10-29 [2] Bioconductor 3.20 (R 4.4.1)
## ANCOMBC * 2.8.0 2024-10-29 [2] Bioconductor 3.20 (R 4.4.1)
## 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)
## base64enc 0.1-3 2015-07-28 [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)
## bit 4.5.0.1 2024-12-03 [2] CRAN (R 4.4.1)
## bit64 4.5.2 2024-09-22 [2] CRAN (R 4.4.1)
## bitops 1.0-9 2024-10-03 [2] CRAN (R 4.4.1)
## boot 1.3-31 2024-08-28 [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)
## cellranger 1.1.0 2016-07-27 [2] CRAN (R 4.4.0)
## checkmate 2.3.2 2024-07-29 [2] CRAN (R 4.4.0)
## class 7.3-22 2023-05-03 [2] CRAN (R 4.4.1)
## cli 3.6.4 2025-02-13 [2] CRAN (R 4.4.1)
## clue 0.3-66 2024-11-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)
## corncob * 0.4.1 2024-01-10 [2] CRAN (R 4.4.0)
## crayon 1.5.3 2024-06-20 [2] CRAN (R 4.4.0)
## CVXR 1.0-15 2024-11-07 [2] CRAN (R 4.4.1)
## 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)
## DescTools 0.99.58 2024-11-08 [2] CRAN (R 4.4.1)
## detectseparation 0.3 2022-08-26 [2] CRAN (R 4.4.0)
## digest 0.6.37 2024-08-19 [2] CRAN (R 4.4.1)
## directlabels 2024.1.21 2024-01-24 [2] CRAN (R 4.4.0)
## doParallel 1.0.17 2022-02-07 [2] CRAN (R 4.4.0)
## doRNG 1.8.6 2023-01-16 [2] CRAN (R 4.4.0)
## dplyr * 1.1.4 2023-11-17 [2] CRAN (R 4.4.0)
## e1071 1.7-16 2024-09-16 [2] CRAN (R 4.4.1)
## energy 1.7-12 2024-08-24 [2] CRAN (R 4.4.1)
## evaluate 1.0.3 2025-01-10 [1] CRAN (R 4.4.1)
## Exact 3.3 2024-07-21 [2] CRAN (R 4.4.0)
## expm 1.0-0 2024-08-19 [2] 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)
## fBasics 4041.97 2024-08-19 [2] CRAN (R 4.4.1)
## 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)
## foreign 0.8-87 2024-06-26 [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)
## GGally * 2.2.1 2024-02-14 [2] CRAN (R 4.4.0)
## ggcorrplot * 0.1.4.1 2023-09-05 [2] CRAN (R 4.4.0)
## ggplot2 * 3.5.1 2024-04-23 [2] CRAN (R 4.4.0)
## ggrepel 0.9.6 2024-09-07 [2] CRAN (R 4.4.1)
## ggstats 0.7.0 2024-09-22 [2] CRAN (R 4.4.1)
## gld 2.6.6 2022-10-23 [2] CRAN (R 4.4.0)
## glue 1.8.0 2024-09-30 [2] CRAN (R 4.4.1)
## gmp 0.7-5 2024-08-23 [2] CRAN (R 4.4.1)
## gridExtra 2.3 2017-09-09 [2] CRAN (R 4.4.0)
## gsl 2.1-8 2023-01-24 [2] CRAN (R 4.4.0)
## 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)
## gtools 3.9.5 2023-11-20 [2] CRAN (R 4.4.0)
## haven 2.5.4 2023-11-30 [2] CRAN (R 4.4.0)
## here * 1.0.1 2020-12-13 [2] CRAN (R 4.4.0)
## Hmisc 5.2-1 2024-12-02 [2] CRAN (R 4.4.1)
## hms 1.1.3 2023-03-21 [2] CRAN (R 4.4.0)
## htmlTable 2.4.3 2024-07-21 [2] CRAN (R 4.4.0)
## htmltools 0.5.8.1 2024-04-04 [2] CRAN (R 4.4.0)
## htmlwidgets 1.6.4 2023-12-06 [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)
## lme4 1.1-35.5 2024-07-03 [2] CRAN (R 4.4.0)
## lmerTest 3.1-3 2020-10-23 [2] CRAN (R 4.4.0)
## lmom 3.2 2024-09-30 [2] CRAN (R 4.4.1)
## lpSolveAPI 5.5.2.0-17.12 2024-07-19 [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)
## microbiome 1.28.0 2024-10-29 [2] Bioconductor 3.20 (R 4.4.1)
## MicrobiomeStat * 1.2 2024-04-01 [2] CRAN (R 4.4.0)
## minqa 1.2.8 2024-08-17 [2] CRAN (R 4.4.0)
## modeest 2.4.0 2019-11-18 [2] CRAN (R 4.4.0)
## multcomp 1.4-26 2024-07-18 [2] CRAN (R 4.4.0)
## 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)
## mvtnorm 1.3-2 2024-11-04 [2] CRAN (R 4.4.1)
## NADA * 1.6-1.1 2020-03-22 [2] CRAN (R 4.4.0)
## nlme 3.1-166 2024-08-14 [2] CRAN (R 4.4.0)
## nloptr 2.1.1 2024-06-25 [2] CRAN (R 4.4.0)
## nnet 7.3-19 2023-05-03 [2] CRAN (R 4.4.1)
## numDeriv 2016.8-1.1 2019-06-06 [2] CRAN (R 4.4.0)
## patchwork * 1.3.0 2024-09-16 [2] CRAN (R 4.4.1)
## 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)
## proxy 0.4-27 2022-06-09 [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)
## quadprog 1.5-8 2019-11-20 [2] CRAN (R 4.4.0)
## 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)
## rbibutils 2.3 2024-10-04 [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)
## RcppZiggurat 0.1.6 2020-10-20 [2] CRAN (R 4.4.0)
## Rdpack 2.6.2 2024-11-15 [2] CRAN (R 4.4.1)
## readr * 2.1.5 2024-01-10 [2] CRAN (R 4.4.0)
## readxl 1.4.3 2023-07-06 [2] CRAN (R 4.4.0)
## registry 0.5-1 2019-03-05 [2] CRAN (R 4.4.0)
## reshape2 1.4.4 2020-04-09 [2] CRAN (R 4.4.0)
## Rfast 2.1.0 2023-11-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)
## Rmpfr 1.0-0 2024-11-18 [2] CRAN (R 4.4.1)
## rmutil 1.1.10 2022-10-27 [2] CRAN (R 4.4.0)
## rngtools 1.5.2 2021-09-20 [2] CRAN (R 4.4.0)
## ROI 1.0-1 2023-04-20 [2] CRAN (R 4.4.0)
## ROI.plugin.lpsolve 1.0-2 2023-07-07 [2] CRAN (R 4.4.0)
## rootSolve 1.8.2.4 2023-09-21 [2] CRAN (R 4.4.0)
## rpart 4.1.23 2023-12-05 [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)
## rstudioapi 0.17.1 2024-10-22 [2] CRAN (R 4.4.1)
## Rtsne 0.17 2023-12-07 [2] CRAN (R 4.4.0)
## 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)
## sandwich 3.1-1 2024-09-15 [2] CRAN (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)
## slam 0.1-55 2024-11-13 [2] CRAN (R 4.4.1)
## SparseArray 1.6.0 2024-10-29 [2] Bioconductor 3.20 (R 4.4.1)
## spatial 7.3-17 2023-07-20 [2] CRAN (R 4.4.1)
## stable 1.1.6 2022-03-02 [2] CRAN (R 4.4.0)
## stabledist 0.7-2 2024-08-17 [2] CRAN (R 4.4.0)
## statip 0.2.3 2019-11-17 [2] CRAN (R 4.4.0)
## statmod 1.5.0 2023-01-06 [2] CRAN (R 4.4.0)
## 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)
## TH.data 1.1-2 2023-04-17 [2] CRAN (R 4.4.0)
## 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)
## timeDate 4041.110 2024-09-22 [2] CRAN (R 4.4.1)
## timeSeries 4041.111 2024-09-22 [2] CRAN (R 4.4.1)
## truncnorm * 1.0-9 2023-03-20 [2] CRAN (R 4.4.0)
## trust 0.1-8 2020-01-10 [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)
## 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)
## zCompositions * 1.5.0-4 2024-06-19 [2] CRAN (R 4.4.0)
## zlibbioc 1.52.0 2024-10-29 [2] Bioconductor 3.20 (R 4.4.1)
## zoo 1.8-12 2023-04-13 [2] CRAN (R 4.4.0)
##
## [1] /Users/yixuanyang/Library/R/arm64/4.4/library
## [2] /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library
##
## ──────────────────────────────────────────────────────────────────────────────