Diazotrophic community associations with Racomitrium lanuginosum genetic groups
Dennis Escolástico
2025-02-24
Overview
The following analyses correspond to project “Hide-and-seek in the tundra: Accounting for cryptic species uncovers molecular variation in moss clonal populations”. The diazotrophic communities of the moss Racomitrium lanuginosum were analyzed using metabarcoding of the nifH gene.
The phyloseq object and most of the files used for analyses correspond to those of the article: Escolástico-Ortiz, D.A., Blasi, C., Bellenger, JP., Derome, N. & Villarreal-A, J.C. 2023. Differentially abundant bacteria drive the N2-fixation of a widespread moss in the forest-tundra transition zone. Symbiosis 90, 193–211.
0. Load libraries and packages
library(readr);packageVersion("readr")
library(ggplot2);packageVersion("ggplot2")
library(ggpubr)
library(tidyverse)
library(RColorBrewer)
library(viridisLite)
theme_set(theme_bw())
# install.packages("FSA")
library(FSA); packageVersion("FSA")
# install.packages("multcompView")
library(multcompView); packageVersion("multcompView")
library(vegan);packageVersion("vegan")
# Installing Phyloseq
#if (!require("BiocManager", quietly = TRUE))
# install.packages("BiocManager")
#BiocManager::install("phyloseq")
#BiocManager::install("S4Vectors")
# install.packages("speedyseq")
#install.packages("devtools")
library(phyloseq); packageVersion("phyloseq")
# install.packages("remotes")
#remotes::install_github("mikemc/speedyseq")
#install.packages("picante")
library(picante); packageVersion("picante")
library(tidyr)
library(MASS)1. Data processing
We load the data to create the phyloseq objects: OTU table, Taxonomy table, Sample data and R object previously produced.
# write.table(ps_nifH.norm@tax_table,"Tax_table_moss_nifH_genotype-microbiome.tsv",sep= "\t", row.names = TRUE)
# Updated taxonomy table
# Load the data
Tax_table_moss_nifH<- read_tsv("Tax_table_moss_nifH_genotype-microbiome.tsv")## Rows: 105 Columns: 8
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (8): #OTU ID, Kingdom, Phylum, Class, Order, Family, Genus, Species
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.# Define row names
Tax_table_moss_nifH<-Tax_table_moss_nifH %>%
tibble::column_to_rownames("#OTU ID")
#Transform to matrix, otherwise the row names will not match.
Tax_table_moss_nifH<- as.matrix(Tax_table_moss_nifH)
Tax_table_moss_nifH<- tax_table(Tax_table_moss_nifH)
# Insert table
p_nifH_2@tax_table <-Tax_table_moss_nifH
p_nifH_2## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 105 taxa and 75 samples ]
## sample_data() Sample Data: [ 75 samples by 8 sample variables ]
## tax_table() Taxonomy Table: [ 105 taxa by 7 taxonomic ranks ]
## refseq() DNAStringSet: [ 105 reference sequences ]# Update the information adding the genetic group of each moss sample.
Gen_MLG_nifH <- read.table("nifH_genotype_MLG.txt",head=T)
nrow(Gen_MLG_nifH) # Check the number of rows in Gen_MLG## [1] 75sample_names(p_nifH_2) # Check the sample names in the phyloseq object## [1] "JC2788" "JC2789" "JC2790" "JC2791" "JC2792" "JC2793" "JC2794" "JC2795"
## [9] "JC2796" "JC2797" "JC2798" "JC2799" "JC2800" "JC2801" "JC2802" "JC2803"
## [17] "JC2804" "JC2805" "JC2806" "JC2807" "JC2926" "JC2927" "JC2928" "JC2929"
## [25] "JC2930" "JC2931" "JC2932" "JC2933" "JC2934" "JC2935" "JC2936" "JC2937"
## [33] "JC2938" "JC2939" "JC2940" "JC2941" "JC2942" "JC2943" "JC2944" "JC2945"
## [41] "JC2946" "JC2947" "JC2948" "JC2949" "JC2950" "JC2951" "JC2952" "JC2953"
## [49] "JC2954" "JC2955" "JC2956" "JC2957" "JC2958" "JC2959" "JC2960" "JC2961"
## [57] "JC2962" "JC2963" "JC3084" "JC3085" "JC3086" "JC3087" "JC3088" "JC3089"
## [65] "JC3090" "JC3091" "JC3210" "JC3211" "JC3212" "JC3213" "JC3214" "JC3215"
## [73] "JC3216" "JC3217" "JC3218"length(sample_names(p_nifH_2)) # Check the number of samples## [1] 75Gen_MLG_nifH <- Gen_MLG_nifH[match(p_nifH_2@sam_data$Voucher, Gen_MLG_nifH$Confir), ]
sample_data(p_nifH_2)$Confirm <- Gen_MLG_nifH[,1]
sample_data(p_nifH_2)$Genotype <- Gen_MLG_nifH[,2]
sample_data(p_nifH_2)$MLG <- as.factor(Gen_MLG_nifH[,3])
sample_data(p_nifH_2)$Lat <- Gen_MLG_nifH[,4]
sample_data(p_nifH_2)$Lon <- Gen_MLG_nifH[,5]
sample_data(p_nifH_2)$GBS <- Gen_MLG_nifH[,6]
sample_data(p_nifH_2)$ID <- Gen_MLG_nifH[,7]
sample_data(p_nifH_2)$Gen_Habitat <- Gen_MLG_nifH[,8]
# Subset data and clean
nifH_gbs<- subset_samples(p_nifH_2, GBS=="TRUE")
table(tax_table(nifH_gbs)[, "Phylum"], exclude = NULL)##
## Cyanobacteria Proteobacteria Verrucomicrobia
## 77 27 1nifH_gbs <-prune_taxa(taxa_sums(nifH_gbs) >=1,nifH_gbs)
nifH_gbs## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 105 taxa and 53 samples ]
## sample_data() Sample Data: [ 53 samples by 16 sample variables ]
## tax_table() Taxonomy Table: [ 105 taxa by 7 taxonomic ranks ]
## refseq() DNAStringSet: [ 105 reference sequences ]2. Alpha diversity
2.1 Rarefaction
Now, we rarefied the data using the minimum number of sequences for all samples to equally compare alpha diversity among groups.
# Rarefaction curves
otu_tab_nifH <- otu_table(nifH_gbs)
class(otu_tab_nifH) <- "matrix" # as.matrix() will do nothing## Warning in class(otu_tab_nifH) <- "matrix": Setting class(x) to "matrix" sets
## attribute to NULL; result will no longer be an S4 object## you get a warning here, but this is what we need to have
otu_tab_nifH <- t(otu_tab_nifH) # transpose observations to rows
rare <- rarecurve(otu_tab_nifH, step=50, lwd=1, ylab="OTU", label=F)
sample_sums(nifH_gbs)## JC2790 JC2791 JC2792 JC2793 JC2794 JC2795 JC2796 JC2797 JC2798 JC2800 JC2801
## 1436 13125 8923 3205 6701 2979 2689 5834 1369 13524 2960
## JC2802 JC2803 JC2804 JC2805 JC2806 JC2807 JC2926 JC2927 JC2929 JC2930 JC2933
## 10475 4796 6204 5303 4867 10906 8699 13575 6817 7691 9133
## JC2937 JC2938 JC2939 JC2940 JC2943 JC2944 JC2945 JC2947 JC2948 JC2949 JC2950
## 5225 3336 6289 7942 2625 8255 2671 9799 8927 8798 3867
## JC2951 JC2952 JC2953 JC2955 JC2956 JC2958 JC2959 JC2960 JC2961 JC2962 JC2963
## 8518 14466 6931 3805 2811 11300 2459 6638 8363 3016 8480
## JC3085 JC3086 JC3087 JC3088 JC3089 JC3090 JC3210 JC3212 JC3214
## 2186 3455 2772 7727 3138 11249 5509 5046 18143# Rarefied data without replacement
set.seed(1)
ps.rarefied_nifH_gbs <- rarefy_even_depth(nifH_gbs, rngseed=1, sample.size=0.9*min(sample_sums(nifH_gbs)), replace=F)## `set.seed(1)` was used to initialize repeatable random subsampling.## Please record this for your records so others can reproduce.## Try `set.seed(1); .Random.seed` for the full vector## ...## 3OTUs were removed because they are no longer
## present in any sample after random subsampling## ...# Check the number of sequences per sample
sample_sums(ps.rarefied_nifH_gbs)## JC2790 JC2791 JC2792 JC2793 JC2794 JC2795 JC2796 JC2797 JC2798 JC2800 JC2801
## 1232 1232 1232 1232 1232 1232 1232 1232 1232 1232 1232
## JC2802 JC2803 JC2804 JC2805 JC2806 JC2807 JC2926 JC2927 JC2929 JC2930 JC2933
## 1232 1232 1232 1232 1232 1232 1232 1232 1232 1232 1232
## JC2937 JC2938 JC2939 JC2940 JC2943 JC2944 JC2945 JC2947 JC2948 JC2949 JC2950
## 1232 1232 1232 1232 1232 1232 1232 1232 1232 1232 1232
## JC2951 JC2952 JC2953 JC2955 JC2956 JC2958 JC2959 JC2960 JC2961 JC2962 JC2963
## 1232 1232 1232 1232 1232 1232 1232 1232 1232 1232 1232
## JC3085 JC3086 JC3087 JC3088 JC3089 JC3090 JC3210 JC3212 JC3214
## 1232 1232 1232 1232 1232 1232 1232 1232 12322.2 Alpha diversity estimates
# Diversity indexes. Disbale phangorn or will not recognize the function.
Shannon.nifH <- diversity(t(otu_table(ps.rarefied_nifH_gbs)), index = "shannon")
Simpson.nifH <- diversity(t(otu_table(ps.rarefied_nifH_gbs)), index = "simpson")
Pielou.nifH <- Shannon.nifH/log(specnumber(t(otu_table(ps.rarefied_nifH_gbs))))
variables_nifH <- sample_data(ps.rarefied_nifH_gbs)
# Merge into one dataframe
Alpha_nifH <- data.frame(variables_nifH,Shannon.nifH,Simpson.nifH,Pielou.nifH)
names(Alpha_nifH)[names(Alpha_nifH) == "Shannon.nifH"] <- "Shannon"
names(Alpha_nifH)[names(Alpha_nifH) == "Simpson.nifH"] <- "Simpson"
names(Alpha_nifH)[names(Alpha_nifH) == "Pielou.nifH"] <- "Pielou"
Alpha_nifH$Plot <- as.factor(Alpha_nifH$Plot)
# Faith phylogenetic diversity
library(phangorn)##
## Attaching package: 'phangorn'## The following objects are masked from 'package:vegan':
##
## diversity, treedistclass(rooted_tree_nifH) # Check if the tree is rooted## [1] "phylo"getRoot(rooted_tree_nifH)## [1] 106Faith_nifH<- pd(t(otu_table(ps.rarefied_nifH_gbs)),rooted_tree_nifH,include.root=F)
Alpha_nifH <-data.frame(Alpha_nifH,Faith_nifH)
names(Alpha_nifH)[names(Alpha_nifH) == "PD"] <- "Faith"
names(Alpha_nifH)[names(Alpha_nifH) == "SR"] <- "Observed"
# Check the number of samples per category
table(Alpha_nifH$Habitat,Alpha_nifH$Genotype)##
## AB CD
## Forest tundra 12 19
## Shrub tundra 7 15# Manipulate data
data_long_nifH_gen <- Alpha_nifH %>%
pivot_longer(cols = c(Shannon,Pielou,Faith),
names_to = "Alpha_Index",
values_to = "Alpha_Value")
# Order levels
data_long_nifH_gen$Alpha_Index<-factor(data_long_nifH_gen$Alpha_Index, levels = c("Shannon", "Pielou","Faith"))
# Create plot
Plot_alpha_GBS_nifH <- ggplot(data_long_nifH_gen , aes(x = Habitat, y = Alpha_Value, fill = Genotype)) +
geom_boxplot() +
facet_wrap(~ Alpha_Index, scales = "free_y") +
labs(title = "nifH gene",
x = NULL,
y = "Alpha diversity value",
fill = "Genetic group") +
scale_fill_manual(values=c("gray","#4d4d4d")) +
theme(legend.position = "bottom", # Placing legend at the bottom
plot.title = element_text(face= "bold"), # Centering plot title
panel.grid.major = element_blank(), # Remove major grid lines
panel.grid.minor = element_blank(),
panel.background = element_blank(),
panel.border = element_blank(),
axis.line = element_line(colour = "black"),
axis.line.x.top = element_line(colour = "black"),
strip.text = element_text(size=10, face= "bold"),
strip.background=element_rect(colour="white", fill = NA),
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
Plot_alpha_GBS_nifH
# ggsave("Figure_3_Alpha_diversity_nifH_Gen_habitat_rarefied.tiff", units="mm", width=160, height=150, dpi=300, compression = "lzw")2.3 Effect of variables on alpha diversity (Linear Mixed Models)
We can evaluate the effect of Habitat and Genetic group, an their interaction on alpha diversity while accounting for potential random effects.
If the normality of residuals is not respected, we could try LMM or GLMM.
# 1. Choosing the GLMM family
# Since alpha diversity indices (e.g., Shannon, Simpson) are continuous variables, you can use:
# - Gaussian (normal) distribution with identity link (if residuals are normal).
# - Gamma or inverse Gaussian if data is right-skewed.
# - Poisson or negative binomial for count-based indices (e.g., Observed richness).
# 2. Running the LMMs (Linear Mixed Models)
library(lme4)
library(lmerTest)
# Models
lmm_Shannon <- lmer(Shannon ~ Habitat * Genotype + (1 | Plot), data = Alpha_nifH)
lmm_Pielou <- lmer(Pielou ~ Habitat * Genotype + (1 | Plot), data = Alpha_nifH)
lmm_Faith <- lmer(Faith ~ Habitat * Genotype + (1 | Plot), data = Alpha_nifH)
# Checking the model summary
summary(lmm_Shannon)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Shannon ~ Habitat * Genotype + (1 | Plot)
## Data: Alpha_nifH
##
## REML criterion at convergence: 83
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.00291 -0.32075 0.04702 0.51054 1.22155
##
## Random effects:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 0.1938 0.4402
## Residual 0.1445 0.3802
## Number of obs: 53, groups: Plot, 24
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.61984 0.17338 36.94734 9.343 2.87e-11
## HabitatShrub tundra 0.09399 0.27978 40.17915 0.336 0.739
## GenotypeCD -0.02449 0.17276 38.76063 -0.142 0.888
## HabitatShrub tundra:GenotypeCD 0.05401 0.27142 36.62610 0.199 0.843
##
## (Intercept) ***
## HabitatShrub tundra
## GenotypeCD
## HabitatShrub tundra:GenotypeCD
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) HbttSt GntyCD
## HbttShrbtnd -0.620
## GenotypeCD -0.605 0.375
## HbttStn:GCD 0.385 -0.646 -0.636summary(lmm_Pielou)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Pielou ~ Habitat * Genotype + (1 | Plot)
## Data: Alpha_nifH
##
## REML criterion at convergence: -38.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0831 -0.2391 0.1615 0.5987 1.2618
##
## Random effects:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 0.01525 0.1235
## Residual 0.01245 0.1116
## Number of obs: 53, groups: Plot, 24
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.520851 0.049720 36.441398 10.476 1.53e-12
## HabitatShrub tundra 0.041104 0.080363 39.859547 0.511 0.612
## GenotypeCD -0.002028 0.050370 38.670618 -0.040 0.968
## HabitatShrub tundra:GenotypeCD 0.029159 0.079220 36.403665 0.368 0.715
##
## (Intercept) ***
## HabitatShrub tundra
## GenotypeCD
## HabitatShrub tundra:GenotypeCD
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) HbttSt GntyCD
## HbttShrbtnd -0.619
## GenotypeCD -0.615 0.380
## HbttStn:GCD 0.391 -0.656 -0.636summary(lmm_Faith )## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Faith ~ Habitat * Genotype + (1 | Plot)
## Data: Alpha_nifH
##
## REML criterion at convergence: 70.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.76859 -0.44687 -0.07352 0.49883 1.46323
##
## Random effects:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 0.1722 0.4150
## Residual 0.1047 0.3235
## Number of obs: 53, groups: Plot, 24
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.17311 0.15590 35.67728 13.939 5.25e-16
## HabitatShrub tundra 0.18615 0.25065 38.90587 0.743 0.462
## GenotypeCD -0.09138 0.14914 36.95745 -0.613 0.544
## HabitatShrub tundra:GenotypeCD -0.11346 0.23376 34.87360 -0.485 0.630
##
## (Intercept) ***
## HabitatShrub tundra
## GenotypeCD
## HabitatShrub tundra:GenotypeCD
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) HbttSt GntyCD
## HbttShrbtnd -0.622
## GenotypeCD -0.581 0.361
## HbttStn:GCD 0.370 -0.621 -0.638# Load the DHARMa package
# install.packages("DHARMa") # If not installed
library(DHARMa)## Warning: package 'DHARMa' was built under R version 4.4.2# p > 0.05 indicates normality
# p > 0.05 means variance is equal
# Shannon
res_lmm_Shannon <- simulateResiduals(lmm_Shannon)
plot(res_lmm_Shannon)
# Pielou
res_lmm_Pielou <- simulateResiduals(lmm_Pielou)
plot(res_lmm_Pielou)
# Faith
res_lmm_Faith <- simulateResiduals(lmm_Faith)
plot(res_lmm_Faith)
# Key assumption checks
# The residuals should follow a uniform distribution.
# If deviations exist, consider non-Gaussian families.
# testUniformity(res_lmm_Shannon) # Example
# Dispersion Test (Overdispersion)
# Detects if variance is higher than expected.
# Overdispersion suggests a negative binomial or quasi-Poisson model may be better.
# testDispersion(res_lmm_Shannon) # Example
# Check Overall Significan with ANOVA
anova(lmm_Shannon)anova(lmm_Pielou)anova(lmm_Faith)3. Microbial profiling
For diazotrophic profiles, we use the non-rarefied OTU table. In order to plot the bacterial profile, we need to change the format of the data.
# Transform to relative abundance
nifH_gbs.ra<- transform_sample_counts(nifH_gbs, function(x){x / sum(x)})
# taxa_sums(nifH_gbs.ra)
# Format the data
phylo_nifH.ra_table <- psmelt(nifH_gbs.ra)Now, we identify the most abundant phyla. We will group all phyla with an abundance lower than 0.01% as “Other”.
# Calculate total relative abundance for each Phylum. This will result in more than 1.00 per phylum because it sums the relative abundance per sample not per Phylum.
phylum_abundance_2 <- aggregate(Abundance ~ Phylum, data = phylo_nifH.ra_table, FUN = sum)
# Insert a column of the relative abundance per Phylum
phylum_abundance_2 <- mutate(phylum_abundance_2,percent = Abundance / sum(Abundance))
# Check the total abundance and percentage
sum(phylum_abundance_2$Abundance)## [1] 53sum(phylum_abundance_2$percent)## [1] 1# Identify Phyla with less than 0.01% abundance
low_abundance_phyla_2 <- phylum_abundance_2$Phylum[phylum_abundance_2$percent < 0.01]
# Group Phyla with less than 0.01% abundance as "Other (<0.01%)"
phylo_nifH.ra_table$Grouped_Phylum <- ifelse(phylo_nifH.ra_table$Phylum %in% low_abundance_phyla_2, "Other", phylo_nifH.ra_table$Phylum)
phylo_nifH.ra_table$Grouped_Phylum <-as.factor(phylo_nifH.ra_table$Grouped_Phylum)
# Calculate total abundance for each grouped Phylum
grouped_phylum_abundance_2 <- aggregate(Abundance ~ Grouped_Phylum, data = phylo_nifH.ra_table, FUN = sum)
# Reorder levels of Grouped_Phylum based on abundance
phylo_nifH.ra_table$Grouped_Phylum <- factor(phylo_nifH.ra_table$Grouped_Phylum,
levels = grouped_phylum_abundance_2$Grouped_Phylum[order(desc(grouped_phylum_abundance_2$Abundance))])
# Reorder Phyla< 0.01
phylo_nifH.ra_table$Grouped_Phylum <- relevel(phylo_nifH.ra_table$Grouped_Phylum,
"Other")We plot the most abundant phyla
# Plot relative abundance per Phylum
BarPlot_nifH_GBS_phylum <- phylo_nifH.ra_table %>%
ggplot(aes(x = Sample, y = Abundance, fill = Grouped_Phylum)) +
geom_bar(stat = "identity", position="stack") +
labs(x = "Samples",
y = "Relative Abundance",title= "nifH gene",fill = "Phylum") +
facet_grid(~ Habitat + Genotype, scales = "free", space='free') +
theme_bw()+
scale_fill_manual(values=Contrasting_palette_1)+
theme(
axis.text.y = element_text(size = 8,),
axis.text.x = element_blank(),
legend.text = element_text(size = 8),
strip.text = element_text(size = 8),
panel.background = element_blank(),
panel.border = element_blank(),
axis.ticks.x = element_blank(),
panel.grid.major = element_blank(),panel.grid.minor = element_blank(),
strip.background=element_rect(colour="white")
)
BarPlot_nifH_GBS_phylum
# ggsave("Figure_Microbial_profile_Genotype_nifH_Phylum.tiff", units="mm", width=180, height=120, dpi=150, compression = "lzw")Now, we can find the most abundant genera with a similar analysis.
# Calculate total relative abundance for each Genus. This will result in more than 1.00 per genus because it sums the relative abundance per sample not per Genus.
genus_abundance_4 <- aggregate(Abundance ~ Genus, data = phylo_nifH.ra_table, FUN = sum)
# Insert a column of the total relative abundance per Genus
genus_abundance_4 <- mutate(genus_abundance_4, percent = Abundance / sum(Abundance))
# Total OTU counts
sum(genus_abundance_4$Abundance)## [1] 53# Total relative abundance of all genera. Should be 1.
sum(genus_abundance_4$percent)## [1] 1# Identify the top 20 most abundant genera
top_20_genera_Moss_nifH <- genus_abundance_4 %>%
arrange(desc(Abundance)) %>%
head(20) %>%
pull(Genus)
# Reorder levels of Grouped_Genus based on abundance
phylo_nifH.ra_table$Genus <- factor(phylo_nifH.ra_table$Genus,
levels = genus_abundance_4$Genus[order(desc(genus_abundance_4$Abundance))])
# Subgroup to obtain relative abundance per group
ps_rare_moss_nifH.ra_table_Samples <-phylo_nifH.ra_table %>%
group_by(Habitat, Genotype) %>%
mutate(Abundance = Abundance / sum(Abundance))We plot the most abundant bacterial genera.
# Produce the plot
BarPlot_nifH_GBS_genus <- phylo_nifH.ra_table %>%
ggplot(aes(x = Sample, y = Abundance, fill = Genus)) +
geom_bar(stat = "identity", position="stack") +
labs(x = "Samples",
y = "Relative Abundance",title = "nifH gene", fill = "Genus") +
facet_grid(~ Habitat + Genotype, scales = "free", space='free') +
theme_bw()+
scale_fill_manual(values=Contrasting_palette_2)+
theme(
axis.text.y = element_text(size = 8,),
axis.text.x = element_blank(),
legend.key.size = unit(0.5, 'cm'),
legend.text = element_text(size = 8),
legend.justification.right = c(0,0.5),
strip.text = element_text(size = 8),
panel.background = element_blank(),
panel.border = element_blank(),
axis.ticks.x = element_blank(),
panel.grid.major = element_blank(),panel.grid.minor = element_blank(),
strip.background=element_rect(colour="gray"),
) + guides(fill=guide_legend(ncol =1))
BarPlot_nifH_GBS_genus
# ggsave("Figure_1-1_Relative_abundance_moss_Genotype_Habitat_Genus_nifH.tiff", units="mm", width=230, height=200, dpi=300, compression = "lzw")We could present the same bar plot by pooling all the samples per habitat and genetic group
# Plot the mean per habitat and genetic group
# Subgroup to obtain relative abundance per group
phylo_nifH.ra_table_Samples_Habitat_Genetic <-phylo_nifH.ra_table %>%
group_by(Habitat,Genotype) %>%
mutate(Abundance = Abundance / sum(Abundance))
# Produce the plot
BarPlot_nifH_GBS_genus_Genotype <- phylo_nifH.ra_table_Samples_Habitat_Genetic %>%
ggplot(aes(x = Genotype, y = Abundance, fill = Genus)) +
geom_bar(stat = "identity", position="stack") +
labs(x = "Genetic group",
y = "Relative Abundance",fill = "Genus") +
facet_grid(~ Habitat, scales = "free", space='free') +
theme_bw()+
scale_fill_manual(values=Contrasting_palette_2)+
theme(
axis.text.y = element_text(size = 8,),
legend.text = element_text(size = 8),
legend.justification.right = c(0,0.5),
strip.text = element_text(size = 8),
legend.key.size = unit(0.3, 'cm'),
panel.background = element_blank(),
panel.border = element_blank(),
axis.ticks.x = element_blank(),
panel.grid.major = element_blank(),panel.grid.minor = element_blank(),
strip.background=element_rect(colour=NA),
#strip.placement = "outside",
) + guides(fill=guide_legend(ncol =1))
BarPlot_nifH_GBS_genus_Genotype
# To put the strips in the bottom use switch = "x"
# ggsave("Figure_1-2_Relative_abundance_moss_pooled_Genotype_Habitat_Genus_nifH.tiff", units="mm", width=230, height=200, dpi=300, compression = "lzw")4. Beta diversity
We analyse the beta diversity of the bacterial community.We use the non-rarefied data but we normalize the ASV counts using Geometric Mean Pairwise Ratios (GMPR).
4.1 Ordination (NMDS)
We evaluate the variation of the bacterial communities using a non-metric multidimensional scaling ordination.
# NMDS using Bray-Curtis. Genotype clustering
# Ordination using normalized data
# Normalization (GMPR)
# Extract and transform the non rarefaction ASV table
ASV_table_nifH_gbs <- as(otu_table(nifH_gbs), "matrix")
# Calculating the GMPR size factors
# install.packages("GUniFrac")
library(GUniFrac)## Registered S3 method overwritten by 'rmutil':
## method from
## print.response httrASV_table_nifH_gbs_GMPR_factors <- GMPR(ASV_table_nifH_gbs)
# Normalize using the size factors. Divide the sample counts by the size factors.
ASV_table_nifH_gbs_normalized_GMPR <- sweep(ASV_table_nifH_gbs, 2, ASV_table_nifH_gbs_GMPR_factors[colnames(ASV_table_nifH_gbs)], FUN = "/")
# Calculate Bray-Distance matrix
bray_dist_nifH <- vegdist(t(ASV_table_nifH_gbs_normalized_GMPR), method = "bray", binary = FALSE)
# Check quality of the data
library(pheatmap)
pheatmap(as.matrix(bray_dist_nifH))
# Plot relative abundance of the top taxa
# plot_bar(nifH_gbs, x = "Sample", fill = "Phylum") +
# theme(axis.text.x = element_text(angle = 90, hjust = 1))
# Transform to ASV_Table
ASV_table_nifH_gbs_normalized_GMPR<- otu_table(ASV_table_nifH_gbs_normalized_GMPR,taxa_are_rows = TRUE)
# Phyloseq normalized (GMPR)
ps_nifH_gbs_norm_GMPR <-phyloseq(ASV_table_nifH_gbs_normalized_GMPR,nifH_gbs@tax_table,nifH_gbs@sam_data)
ps_nifH_gbs_norm_GMPR ## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 105 taxa and 53 samples ]
## sample_data() Sample Data: [ 53 samples by 16 sample variables ]
## tax_table() Taxonomy Table: [ 105 taxa by 7 taxonomic ranks ]ps_nifH_gbs_norm_GMPR <- subset_taxa(ps_nifH_gbs_norm_GMPR , !is.na(Phylum))
ps_nifH_gbs_norm_GMPR ## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 105 taxa and 53 samples ]
## sample_data() Sample Data: [ 53 samples by 16 sample variables ]
## tax_table() Taxonomy Table: [ 105 taxa by 7 taxonomic ranks ]Using the Bray-Curtis matrix of the normalized ASV counts, we perform a NMDS ordination.
# Ordinate the data
set.seed(123456)
ord.nmds.bray_nifH <- ordinate(ps_nifH_gbs_norm_GMPR, method="NMDS", distance = bray_dist_nifH)## Run 0 stress 0.2097578
## Run 1 stress 0.2097578
## ... New best solution
## ... Procrustes: rmse 7.264923e-06 max resid 3.24014e-05
## ... Similar to previous best
## Run 2 stress 0.2097578
## ... Procrustes: rmse 5.863672e-05 max resid 0.0003888018
## ... Similar to previous best
## Run 3 stress 0.2146192
## Run 4 stress 0.2097578
## ... Procrustes: rmse 2.343818e-05 max resid 0.0001405265
## ... Similar to previous best
## Run 5 stress 0.2097578
## ... Procrustes: rmse 2.75302e-05 max resid 0.0001813603
## ... Similar to previous best
## Run 6 stress 0.2097578
## ... Procrustes: rmse 1.897336e-05 max resid 0.0001054131
## ... Similar to previous best
## Run 7 stress 0.209758
## ... Procrustes: rmse 4.369654e-05 max resid 0.0001823163
## ... Similar to previous best
## Run 8 stress 0.2146192
## Run 9 stress 0.2146192
## Run 10 stress 0.2097578
## ... Procrustes: rmse 5.944548e-05 max resid 0.0003938184
## ... Similar to previous best
## Run 11 stress 0.2097958
## ... Procrustes: rmse 0.003441101 max resid 0.01951287
## Run 12 stress 0.2146192
## Run 13 stress 0.2266802
## Run 14 stress 0.2097958
## ... Procrustes: rmse 0.003402881 max resid 0.01940152
## Run 15 stress 0.2097578
## ... Procrustes: rmse 1.614679e-05 max resid 0.0001067605
## ... Similar to previous best
## Run 16 stress 0.2097578
## ... Procrustes: rmse 7.26842e-05 max resid 0.0004857443
## ... Similar to previous best
## Run 17 stress 0.2705758
## Run 18 stress 0.2097578
## ... Procrustes: rmse 1.042647e-05 max resid 4.557877e-05
## ... Similar to previous best
## Run 19 stress 0.2146192
## Run 20 stress 0.2146192
## *** Best solution repeated 10 timesord.nmds.bray_nifH$stress## [1] 0.2097578stressplot(ord.nmds.bray_nifH)
NMDS.bray_nifH_gen <- plot_ordination(ps_nifH_gbs_norm_GMPR, ord.nmds.bray_nifH, color="Genotype", shape = "Habitat")
NMDS_bray_nifH_gen <- NMDS.bray_nifH_gen +
stat_ellipse(type = "norm", linetype = 2) +
theme( panel.grid.major = element_blank(),panel.grid.minor = element_blank()) +
geom_point(size=2) +
scale_color_manual(values=c("gray","#363636"))+
labs(caption="2D, stress = 0.2401", color = "Genetic group") +
guides(color = guide_legend(order = 1), shape = guide_legend(order = 2))
NMDS_bray_nifH_gen 
# ggsave("Figure_2_NMDS_moss_Genotype_Habitat_nifH.tiff", units="mm", width=150, height=150, dpi=300, compression = "lzw")4.2 PERMANOVA
We evaluate how the variation is affected by the habitat and genetic group (genotype).
# PERMANOVA
permanova_nifH_gen <- adonis2(t(ps_nifH_gbs_norm_GMPR@otu_table) ~ ps_nifH_gbs_norm_GMPR@sam_data$Habitat * ps_nifH_gbs_norm_GMPR@sam_data$Genotype, permutations = perm, method = "bray",strata = ps_nifH_gbs_norm_GMPR@sam_data$Plot)
permanova_nifH_gen# Check variance heterogenity or dispersion among groups
set.seed(123)
BC_nifH_beta_disper_1 <- betadisper(bray_dist_nifH,ps_nifH_gbs_norm_GMPR@sam_data$Gen_Habitat)
permutest(BC_nifH_beta_disper_1)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 3 0.02258 0.0075282 0.5839 999 0.619
## Residuals 49 0.63172 0.0128923# The test indicate that the there is not enough proofs
# to confirm that the homogeneity of multivariate dispersion between groups differ5. Differentially abundant bacteria
We can identify groups that are more abundant between habitats or genetic groups. In this case, we use ANCOMBC-2 to identify differential abundant diazotrophic genera with bias corrections.
We use the non-rarefied count because ANCOMBC-2 normalize the data before run the analysis.
The primary output of the ANCOM-BC2 methodology identifies taxa with differential abundance based on the chosen covariate. The results include: 1) log fold changes, 2) standard errors, 3) test statistics, 4) p-values, 5) adjusted p-values, 6) indicators denoting whether the taxon is differentially abundant (TRUE) or not (FALSE), and 7) indicators denoting whether the taxon passed the sensitivity analysis (TRUE) or not (FALSE) ().
#if (!require("BiocManager", quietly = TRUE))
# install.packages("BiocManager")
#BiocManager::install("ANCOMBC")
library(ANCOMBC)
# Check the number of samples per category
table(nifH_gbs@sam_data$Habitat,nifH_gbs@sam_data$Genotype)##
## AB CD
## Forest tundra 12 19
## Shrub tundra 7 15# Check the number of sequences per sample
sum_seq <- rowSums(t(nifH_gbs@otu_table))
sort(sum_seq)## JC2798 JC2790 JC3085 JC2959 JC2943 JC2945 JC2796 JC3087 JC2956 JC2801 JC2795
## 1369 1436 2186 2459 2625 2671 2689 2772 2811 2960 2979
## JC2962 JC3089 JC2793 JC2938 JC3086 JC2955 JC2950 JC2803 JC2806 JC3212 JC2937
## 3016 3138 3205 3336 3455 3805 3867 4796 4867 5046 5225
## JC2805 JC3210 JC2797 JC2804 JC2939 JC2960 JC2794 JC2929 JC2953 JC2930 JC3088
## 5303 5509 5834 6204 6289 6638 6701 6817 6931 7691 7727
## JC2940 JC2944 JC2961 JC2963 JC2951 JC2926 JC2949 JC2792 JC2948 JC2933 JC2947
## 7942 8255 8363 8480 8518 8699 8798 8923 8927 9133 9799
## JC2802 JC2807 JC3090 JC2958 JC2791 JC2800 JC2927 JC2952 JC3214
## 10475 10906 11249 11300 13125 13524 13575 14466 18143plot(sort(sum_seq), ylim=c(0,25000), main=c("Number of counts per each sample"), xlab=c("Sample position"))
# Run ANCOMBC
library(ANCOMBC)
set.seed(123)
out_nifH = ancombc2(
data = nifH_gbs,
tax_level = "Genus",
fix_formula = "Habitat + Genotype",
group = "Genotype",
p_adj_method = "fdr",
global = TRUE,
verbose = TRUE
)## Warning: The group variable has < 3 categories
## The multi-group comparisons (global/pairwise/dunnet/trend) will be deactivated## Obtaining initial estimates ...## Estimating sample-specific biases ...## Loading required package: foreach##
## Attaching package: 'foreach'## The following objects are masked from 'package:purrr':
##
## accumulate, when## Loading required package: rngtools## Sensitivity analysis for pseudo-count addition to 0s: ...## ANCOM-BC2 primary results ...# prv_cut = 0.05, In this case, we already filtered ASVs by prevalence (previous step in Escolastico-Ortiz et al. 2023. Symbiosis.).
# Library size cut less than 1000 reads.
# If we applied the struc_zero we will find taxa that occur differentially based on presence/absence.
# Create a new data frame with the results
res_nifH_gbs <- out_nifH$res
# By default, ANCOM-BC2 does not add pseudo-counts to zero counts, which can
# result in NAs in the bias-corrected log abundances. Users have the option to
# either leave the NAs as they are or replace them with zeros.
# This replacement is equivalent to adding pseudo-counts of ones to the zero counts.
# The p-values from ANCOM-BC2, both with and without sensitivity analysis, exhibit a nearly uniform distribution.
# Notably, enabling sensitivity analysis significantly reduces the number of significant p-values, thereby decreasing the risk of false positives.Now, we create a function that allow to plot the results.
extract_differential_abundance <- function(df_Group, res_data, taxon_col, lfc_col, q_col, grouping_var,
passed_col = NULL) {
# Create another data frame with the relevant grouping results
df_group_res <- res_data %>%
dplyr::select({{taxon_col}}, contains({{grouping_var}}), {{passed_col}})
# Add association based on significance, using custom group names
df_group_res <- df_group_res %>%
dplyr::mutate(
association = dplyr::case_when(
!!rlang::sym(q_col) < 0.05 & !!rlang::sym(lfc_col) > 0 ~ "Enriched",
!!rlang::sym(q_col) < 0.05 & !!rlang::sym(lfc_col) < 0 ~ "Depleted",
TRUE ~ "Not significant"
),
font = ifelse(!is.null(passed_col) & !!rlang::sym(passed_col), "bold", "plain"),
color = ifelse(!is.null(passed_col) & !!rlang::sym(passed_col), "blue", "black"),
)
# Extract differentially abundant groups
results_sig <- df_group_res[df_group_res[[q_col]] < 0.05,]
table(results_sig$association)
# Order the data
results_sig <- results_sig[order(results_sig$association, results_sig[[lfc_col]]),]
desired_order <- results_sig[[taxon_col]]
# Convert 'Genera' (taxon_col) to a factor with the specified order
results_sig[[taxon_col]] <- factor(results_sig[[taxon_col]], levels = desired_order)
# Return the results as a list
return(list(df_group_res = df_group_res, results_sig = results_sig))
}
# Example usage:
# extract_differential_abundance(df_Genotype, res_16S_gbs, "taxon", "lfc_GenotypeCD", "q_GenotypeCD", "Genotype","AB","CD", passed_col= "passed_ss_GenotypeCD")We use the results of our function to plot the differentially abundant taxa.
# Run Script to produce the data
Results_ANCOMBC2_nifH <- extract_differential_abundance(df_Genotype, res_nifH_gbs, "taxon", "lfc_GenotypeCD", "q_GenotypeCD", "Genotype", passed_col= "passed_ss_GenotypeCD")
# Store only the significant taxa data
Result_total_ANCOMBC2_nifH <- Results_ANCOMBC2_nifH$df_group_res
Result_sig_ANCOMBC2_nifH <- Results_ANCOMBC2_nifH$results_sig
# Plot the significant differential abundant taxa
library(ggrepel) # For better label positioning
# Extract the labels for the significant taxa
sig_labels_nifH <- ifelse(Result_total_ANCOMBC2_nifH$diff_GenotypeCD == TRUE,Result_total_ANCOMBC2_nifH$taxon, NA)
# Create plot for significant abundant taxa
Plot_significant_taxa_nifH <- ggplot(Result_total_ANCOMBC2_nifH, aes(x = lfc_GenotypeCD, y = -log10(q_GenotypeCD))) +
geom_point(aes(color = diff_GenotypeCD)) +
scale_color_manual(values = c("grey", "#4d9221")) +
theme( panel.grid.major = element_blank(),panel.grid.minor = element_blank()) +
labs(x = "Log2 Fold Change", y = "-Log10 Adjusted p-value",color= "Differentially abundant", title= "nifH gene") +
geom_vline(xintercept = 0, color = "darkgrey", linewidth = 0.5) +
geom_vline(xintercept = -1, linetype = "dashed", color = "darkgrey", linewidth = 0.5) + # LFC = -1
geom_vline(xintercept = 1, linetype = "dashed", color = "darkgrey", linewidth = 0.5) + # LFC = 1
geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "darkgrey", linewidth = 0.5) + # P-value = 0.05 threshold
geom_text_repel(aes(label = sig_labels_nifH), na.rm = TRUE,size = 2.5) # Only labels significant genes
Plot_significant_taxa_nifH
# ggsave("Figure_Significant_DA_taxa_Genotype_nifH.tiff", units="mm", width=200, height=150, dpi=150, compression = "lzw")
# THERE IS ONLY ONE DIFFERENTIALLY ABUNDANT GENUS: Pinisolibacter in the Genetic group CD
# The genus did not passed the sensitivity analyses: passed_ss_GenotypeCD = FALSE
Result_total_ANCOMBC2_nifH# Plot log2 fold changes
# Plot_lfc <- ggplot(Results_ANCOMBC2_nifH, aes(x= lfc_GenotypeCD, y = taxon, color = association)) +
# geom_point() +
# geom_segment(aes(x = 0, xend = lfc_GenotypeCD, y = taxon, yend = taxon)) +
# geom_vline(xintercept = 0, size = 0.3) +
# xlab("log2FoldChange") + ylab(NULL) +
# theme_classic() +
# scale_color_manual(values = c("AB" = "gray", "CD" = "#363636")) + # Custom colors
# theme(
# legend.title = element_blank(),
# axis.text.y = element_text(face = Result_sig_ANCOMBC2$font,color = Result_sig_ANCOMBC2$color))
#
# Plot_lfc
# ggsave("Figure_3_ANCOMBC2__Genotype_nifH.tiff", units="mm", width=230, height=120, dpi=150, compression = "lzw")THERE IS ONLY ONE DIFFERENTIALLY ABUNDANT GENUS: Pinisolibacter in the Genetic group CD. The genus did not passed the sensitivity analysis.
Supplementary analyses
We can merge the plot results for both 16S and nifH into combined plots
Microbial_profiles <- ggarrange(BarPlot_bacteria_gbs_genus + rremove ("xlab") , BarPlot_nifH_GBS_genus + rremove ("xlab"),nrow=2, ncol=1, align = "v", labels = c("a","c"))
Microbial_NMDS <- ggarrange(NMDS.BC_Bacteria_gbs, NMDS_bray_nifH_gen , nrow=2, ncol=1, align = "v", common.legend = TRUE, legend = "right", labels = c("b","d"))
Genotype_microbiome_covariation <- ggarrange(Microbial_profiles, Microbial_NMDS , nrow=1, ncol=2, align = "hv", widths = c(2, 1))
# ggsave("Figure_X_Genotype_microbiome_covariation_figure.tiff", units="mm", width=430, height=250, dpi=300, compression = "lzw")
Microbial_profiles_2 <- ggarrange(BarPlot_bacteria_gbs_genus_Habitat_genotype + rremove ("xlab") + rremove ("x.text"), BarPlot_nifH_GBS_genus_Genotype,nrow=2, ncol=1, align = "v",labels = c("a","c"))
Genotype_microbiome_covariation_2 <- Genotype_microbiome_covariation <- ggarrange(Microbial_profiles_2, Microbial_NMDS , nrow=1, ncol=2, align = "h")
# ggsave("Figure_X2_Genotype_microbiome_covariation_pooled_figure.tiff", units="mm", width=250, height=250, dpi=300, compression = "lzw")
# ggsave("Figure_X3_Genotype_microbiome_covariation_pooled_figure.tiff", units="mm", width=250, height=150, dpi=300, compression = "lzw")Session information
sessionInfo()## R version 4.4.1 (2024-06-14 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 22631)
##
## Matrix products: default
##
##
## Random number generation:
## RNG: Mersenne-Twister
## Normal: Inversion
## Sample: Rounding
##
## locale:
## [1] LC_COLLATE=English_Canada.utf8 LC_CTYPE=English_Canada.utf8
## [3] LC_MONETARY=English_Canada.utf8 LC_NUMERIC=C
## [5] LC_TIME=English_Canada.utf8
##
## time zone: America/Toronto
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggrepel_0.9.5 doRNG_1.8.6 rngtools_1.5.2
## [4] foreach_1.5.2 ANCOMBC_2.2.2 pheatmap_1.0.12
## [7] GUniFrac_1.8 DHARMa_0.4.7 lmerTest_3.1-3
## [10] lme4_1.1-35.5 Matrix_1.7-0 phangorn_2.12.1
## [13] MASS_7.3-61 picante_1.8.2 nlme_3.1-165
## [16] ape_5.8 phyloseq_1.48.0 vegan_2.6-6.1
## [19] lattice_0.22-6 permute_0.9-7 multcompView_0.1-10
## [22] FSA_0.9.5 viridisLite_0.4.2 RColorBrewer_1.1-3
## [25] lubridate_1.9.3 forcats_1.0.0 stringr_1.5.1
## [28] dplyr_1.1.4 purrr_1.0.2 tidyr_1.3.1
## [31] tibble_3.2.1 tidyverse_2.0.0 ggpubr_0.6.0
## [34] ggplot2_3.5.1 readr_2.1.5
##
## loaded via a namespace (and not attached):
## [1] fs_1.6.5 matrixStats_1.3.0
## [3] DirichletMultinomial_1.46.0 httr_1.4.7
## [5] doParallel_1.0.17 numDeriv_2016.8-1.1
## [7] tools_4.4.1 backports_1.5.0
## [9] utf8_1.2.4 R6_2.5.1
## [11] lazyeval_0.2.2 mgcv_1.9-1
## [13] rhdf5filters_1.16.0 withr_3.0.2
## [15] gridExtra_2.3 cli_3.6.3
## [17] Biobase_2.64.0 sandwich_3.1-1
## [19] labeling_0.4.3 sass_0.4.9
## [21] mvtnorm_1.3-1 proxy_0.4-27
## [23] yulab.utils_0.1.7 foreign_0.8-87
## [25] scater_1.32.1 decontam_1.24.0
## [27] readxl_1.4.3 rstudioapi_0.17.1
## [29] generics_0.1.3 gtools_3.9.5
## [31] vroom_1.6.5 car_3.1-3
## [33] inline_0.3.19 biomformat_1.32.0
## [35] ggbeeswarm_0.7.2 fansi_1.0.6
## [37] DescTools_0.99.57 S4Vectors_0.42.1
## [39] DECIPHER_3.0.0 abind_1.4-8
## [41] lifecycle_1.0.4 multcomp_1.4-26
## [43] yaml_2.3.10 carData_3.0-5
## [45] SummarizedExperiment_1.34.0 rhdf5_2.48.0
## [47] SparseArray_1.4.8 grid_4.4.1
## [49] crayon_1.5.3 beachmat_2.20.0
## [51] pillar_1.9.0 knitr_1.48
## [53] GenomicRanges_1.56.1 statip_0.2.3
## [55] boot_1.3-31 gld_2.6.6
## [57] codetools_0.2-20 fastmatch_1.1-4
## [59] glue_1.7.0 data.table_1.15.4
## [61] MultiAssayExperiment_1.30.3 vctrs_0.6.5
## [63] treeio_1.28.0 Rdpack_2.6.1
## [65] cellranger_1.1.0 gtable_0.3.6
## [67] cachem_1.1.0 xfun_0.46
## [69] rbibutils_2.3 S4Arrays_1.4.1
## [71] modeest_2.4.0 survival_3.7-0
## [73] timeDate_4041.110 SingleCellExperiment_1.26.0
## [75] iterators_1.0.14 bluster_1.14.0
## [77] statmod_1.5.0 gmp_0.7-5
## [79] TH.data_1.1-2 gap_1.6
## [81] bit64_4.0.5 GenomeInfoDb_1.40.1
## [83] fBasics_4041.97 bslib_0.8.0
## [85] irlba_2.3.5.1 vipor_0.4.7
## [87] rpart_4.1.23 DBI_1.2.3
## [89] colorspace_2.1-1 BiocGenerics_0.50.0
## [91] Hmisc_5.2-0 nnet_7.3-19
## [93] ade4_1.7-22 Exact_3.3
## [95] tidyselect_1.2.1 timeSeries_4041.111
## [97] bit_4.0.5 compiler_4.4.1
## [99] htmlTable_2.4.3 BiocNeighbors_1.22.0
## [101] expm_1.0-0 DelayedArray_0.30.1
## [103] checkmate_2.3.2 scales_1.3.0
## [105] quadprog_1.5-8 spatial_7.3-17
## [107] digest_0.6.36 minqa_1.2.8
## [109] rmarkdown_2.29 XVector_0.44.0
## [111] htmltools_0.5.8.1 pkgconfig_2.0.3
## [113] base64enc_0.1-3 sparseMatrixStats_1.16.0
## [115] MatrixGenerics_1.16.0 highr_0.11
## [117] stabledist_0.7-2 fastmap_1.2.0
## [119] rlang_1.1.4 htmlwidgets_1.6.4
## [121] UCSC.utils_1.0.0 DelayedMatrixStats_1.26.0
## [123] farver_2.1.2 jquerylib_0.1.4
## [125] zoo_1.8-12 jsonlite_1.8.8
## [127] energy_1.7-12 BiocParallel_1.38.0
## [129] BiocSingular_1.20.0 magrittr_2.0.3
## [131] Formula_1.2-5 scuttle_1.14.0
## [133] GenomeInfoDbData_1.2.12 Rhdf5lib_1.26.0
## [135] munsell_0.5.1 Rcpp_1.0.13
## [137] viridis_0.6.5 CVXR_1.0-14
## [139] stringi_1.8.4 rootSolve_1.8.2.4
## [141] stable_1.1.6 zlibbioc_1.50.0
## [143] plyr_1.8.9 parallel_4.4.1
## [145] lmom_3.2 Biostrings_2.72.1
## [147] splines_4.4.1 multtest_2.60.0
## [149] hms_1.1.3 igraph_2.0.3
## [151] ggsignif_0.6.4 reshape2_1.4.4
## [153] stats4_4.4.1 ScaledMatrix_1.12.0
## [155] rmutil_1.1.10 evaluate_1.0.1
## [157] nloptr_2.1.1 tzdb_0.4.0
## [159] clue_0.3-65 rsvd_1.0.5
## [161] broom_1.0.7 gap.datasets_0.0.6
## [163] Rmpfr_0.9-5 e1071_1.7-16
## [165] tidytree_0.4.6 rstatix_0.7.2
## [167] class_7.3-22 gsl_2.1-8
## [169] beeswarm_0.4.0 IRanges_2.38.1
## [171] cluster_2.1.6 TreeSummarizedExperiment_2.12.0
## [173] timechange_0.3.0 mia_1.12.0