Microbiome analysis in R - Part III

Note

For a free consultation, please visit our consultation scheduler!

Workshop and Data Overview

The third installment of the microbiome analysis is to explore phyloseq object for abundance visualization and diversity analysis. For this session, we will be using the full phyloseq object from Salomon, JD et al., 2023 so we have more covariants to work with. The data has two main covariates, group (CNT and CPB representing control samples and cardiopulmonary bypass samples) and time (pre and post representing pre-operation and post-operation time points).

Getting Ready

Load R packages we will be using for today’s workshop.

suppressMessages({
        library("tidyverse")
        library("phyloseq")
        library("rstatix")
        library("vegan")
        library("picante")
})
no_of_cores = 1

Copy all needed files to $WORK, create the 16S_Workshop directory if it does not exist, and set it as working directory using commands below.

system("cp -Ru /work/riethoven/jeanjack/16S_Workshop $WORK/", intern=TRUE)
setwd(file.path(Sys.getenv("WORK"), "16S_Workshop"))

Load the phyloseq object constructed from whole dataset of Salomon, JD et al., 2023.

load("./intermediate/SalomonJD_et_al.rda")

Inspect objects loaded into your environment.

ls()

[1] "create_dt"   "no_of_cores" "ps"         

ps

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 4578 taxa and 24 samples ]
sample_data() Sample Data:       [ 24 samples by 6 sample variables ]
tax_table()   Taxonomy Table:    [ 4578 taxa by 7 taxonomic ranks ]
phy_tree()    Phylogenetic Tree: [ 4578 tips and 4576 internal nodes ]
refseq()      DNAStringSet:      [ 4578 reference sequences ]

Inspect and Filter phyloseq object

Let’s have a look at the 5 components of phyloseq object.

• sample_data(ps): sample meta data

• otu_table(ps): ASV table after removing chimera

• taxa_table(ps): taxonomic assignments for those ASVs

• refseq(ps) is the result of multiple sequence alignment with MAFFT

• phy_tree(ps) is the result of phylogenetic tree reconstruction with FastTree

A few screening and filtering steps before final phyloseq object is ready for abundance visualization and diversity analysis.

• Step 1, remove Archaea, unknowns, chloroplasts;

ps.clean <- subset_taxa(ps, Kingdom == "Bacteria") %>%
        subset_taxa(!is.na(Phylum)) %>%
        subset_taxa(!Class %in% c("Chloroplast")) %>%
        subset_taxa(!Family %in% c("mitochondria"))
ps.clean

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 4541 taxa and 24 samples ]
sample_data() Sample Data:       [ 24 samples by 6 sample variables ]
tax_table()   Taxonomy Table:    [ 4541 taxa by 7 taxonomic ranks ]
phy_tree()    Phylogenetic Tree: [ 4541 tips and 4539 internal nodes ]
refseq()      DNAStringSet:      [ 4541 reference sequences ]

We are only interested in bacteria. Removed 37 taxa.

• Step 2, filter taxa (ASV) with low abundance (< 2);

ps.clean.p0 <- filter_taxa(ps.clean, function (x) {sum(x > 0) >= 2}, prune=TRUE)
ps.clean.p0

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 3564 taxa and 24 samples ]
sample_data() Sample Data:       [ 24 samples by 6 sample variables ]
tax_table()   Taxonomy Table:    [ 3564 taxa by 7 taxonomic ranks ]
phy_tree()    Phylogenetic Tree: [ 3564 tips and 3562 internal nodes ]
refseq()      DNAStringSet:      [ 3564 reference sequences ]

Abundance is the sum of the number of times an ASV is observed across all samples. Removed 977 taxa.

Complete metadata are shown below:

sample_meta = sample_data(ps.clean.p0)
sample_meta

Visualizing Abundance

Often an early step in many microbiome analysis is to visualize the abundance of organisms at specific taxonomic ranks. Stacked bar plots and faceted box plots are two ways of doing this.

Absolute abundance (stacked bar plots)

By group (CPB VS CNT)

phyloseq::plot_bar(ps.clean.p0, fill = "Phylum") +
        geom_bar(aes(color = Phylum, fill = Phylum), stat = "identity", position = "stack") +
        facet_wrap(~group, scales = "free_x")

By time (pre VS post)

phyloseq::plot_bar(ps.clean.p0, fill = "Phylum") + 
        geom_bar(aes(color = Phylum, fill = Phylum), stat = "identity", position = "stack") +
        facet_wrap(~time, scales = "free_x")

Note

Notice that post comes before pre, that is because most packages arrange levels of a factor using alphabetical order. To rearrange the levels of factor or to identify which level of a factor is reference, we can go back to getting-ready phase and do this:

sample_meta = sample_data(ps.clean.p0) %>%
        dplyr::as_tibble() %>%
        dplyr::mutate(time = factor(time)) %>% # this step make `time` a factor so we can rearrange its levels
        within(., time <- relevel(time, ref = "pre")) %>% # this step set `pre` as ref level
        as.data.frame()

rownames(sample_meta) = sample_meta$sample 

sample_data(ps.clean.p0) = data.frame(sample_meta) # update the sample_data() of our `phyloseq` object

phyloseq::plot_bar(ps.clean.p0, fill = "Phylum") + 
        geom_bar(aes(color = Phylum, fill = Phylum), stat = "identity", position = "stack") +
        facet_wrap(~time, scales = "free_x")

Absolute abundance (faceted box plot)

Another way to roughly visualize changes is the faceted box plot. Below we have a faceted box plot showing absolute abundance on phylum level facet by time and group.

# Agglomerate to phylum-level and rename
ps.clean.phylum = phyloseq::tax_glom(ps.clean.p0, taxrank = rank_names(ps.clean.p0)[2])
phyloseq::taxa_names(ps.clean.phylum) = phyloseq::tax_table(ps.clean.phylum)[,"Phylum"]
otu_table(ps.clean.phylum)[1:5,1:5]

OTU Table:          [5 taxa and 5 samples]
                     taxa are columns
       Bacteroidota Synergistota Spirochaetota Actinobacteriota Elusimicrobiota
12post         1857            0            25               23               0
12pre         18579            1          1823              461               0
14post        23601            0           490              511               0
14pre          7502            3           190              167               0
20post        24655            2          1844             3691               0

By group (CPB VS CNT)

#Melt and plot
phyloseq::psmelt(ps.clean.phylum) %>%
  ggplot(data = ., aes(x = group, y = Abundance)) +
  geom_boxplot(outlier.shape  = NA) +
  geom_jitter(aes(color = OTU), height = 0, width = .2) +
  labs(x = "", y = "Abundance\n") +
  facet_wrap(~ OTU, scales = "free", nrow = 2) +
  ggtitle("Phylum absolute abundance by time")

By time (pre VS post)

#Melt and plot
phyloseq::psmelt(ps.clean.phylum) %>%
  ggplot(data = ., aes(x = time, y = Abundance)) +
  geom_boxplot(outlier.shape  = NA) +
  geom_jitter(aes(color = OTU), height = 0, width = .2) +
  labs(x = "", y = "Abundance\n") +
  facet_wrap(~ OTU, scales = "free", nrow = 2) +
  ggtitle("Phylum absolute abundance by group")

Relative abundance

One of the challenges we face working with NGS-derived sequence data is that the total number of reads for each sample is not directly tied to the starting quantity of DNA. The total count does not carry any information on the absolute abundance of taxa. As long as the count is sufficiently large, it is just a factor that we want to account for in our analysis and is not of particular interest other than differences across samples can be a source of bias. Therefore, we usually transform absolute abundance to relative abundance as count normalization before any diversity analysis.

ps.clean.re <- transform_sample_counts(ps.clean.p0, function(x) x / sum(x))

By group (CPB VS CNT)

phyloseq::plot_bar(ps.clean.re, fill = "Phylum") + 
        geom_bar(aes(color = Phylum, fill = Phylum), stat = "identity", position = "stack") +
        facet_wrap(~group, scales = "free_x")

By time (pre VS post)

phyloseq::plot_bar(ps.clean.re, fill = "Phylum") + 
        geom_bar(aes(color = Phylum, fill = Phylum), stat = "identity", position = "stack") +
        facet_wrap(~time, scales = "free_x")

Alpha Diversity

Basically, alpha diversity is the within-sample diversity and includes how many organisms are observed (i.e. observed OTUs) and how evenly they are distributed. However, the issue of how to best estimate these quantities using data derived from next-generation sequencing (NGS) is controversial. This is due to two reasons:

• The observed richness in a sample/site is typically underestimated due to inexhaustive sampling. Thus, valid estimators of diversity require extrapolating from the available observations to provide estimates of the unobserved taxa (and also account for the sampling variability).

• Extrapolation estimators require an accurate count of the rare taxa (including singletons) in each sample, which for NGS-based metagenomics studies we typically do not have, as singletons generally cannot be differentiated from sequencing errors using the current best bioinformatics workflows.

It has been argued, however, that diversity metrics can nevertheless be compared between samples because the errors and biases are mostly systematic (i.e. occur with similar rates in all samples). A major underlying assumption here is that abundance structures are the same for the two groups being compared. This is a reasonable assumption when comparing similar environments, but it is hard to know without exhaustive sampling.

Below we will estimate and test for differences according to interested covariates (group or time) using the plug-in estimates for 7 different alpha diversity metrics.

# alpha diversity should be calculated before filtering on abundance and prevalence
tree = phyloseq::phy_tree(ps)
samp = data.frame(phyloseq::otu_table(ps))

#Generate a data.frame with adiv measures
adiv <- data.frame(
  phyloseq::estimate_richness(ps, measures = c("Observed", "Shannon", "Chao1", "Simpson", "InvSimpson", "Fisher")),
  "PD" = picante::pd(samp, tree, include.root=FALSE)[,1],
  dplyr::select(as_tibble(phyloseq::sample_data(ps)), sample, group, time, ID )) %>%
        dplyr::select(-se.chao1)

measures = colnames(adiv)[1:7]
#head(adiv)

By group (CPB VS CNT)

adiv %>%
  group_by( group ) %>%
  dplyr::summarise(median_Observed = median(Observed),
                   median_Chao1 = median(Chao1),
                   median_Shannon = median(Shannon),
                   median_Simpson = median(Simpson),
                   median_InvSimpson = median(InvSimpson),
                   median_Fisher = median(Fisher),
                   median_PD = median(PD))

# A tibble: 2 × 8
  group median_Observed median_Chao1 median_Shannon median_Simpson
  <chr>           <dbl>        <dbl>          <dbl>          <dbl>
1 CNT             1408.        1843.           5.46          0.987
2 CPB             1459         1845.           5.43          0.989
# ℹ 3 more variables: median_InvSimpson <dbl>, median_Fisher <dbl>,
#   median_PD <dbl>

#Plot adiv measures
adiv %>%
  gather(key = metric, value = value, measures) %>% #c("Observed", "Shannon", "Chao1", "Simpson", "PD")
  mutate(metric = factor(metric, levels = measures)) %>% #c("Observed", "Shannon", "Chao1", "Simpson", "PD")
  ggplot(aes(x = group, y = value)) +
  geom_boxplot(outlier.color = NA) +
  geom_jitter(aes(color = group), height = 0, width = .2) +
  labs(x = "", y = "Alpha Diversity Measures") +
  facet_wrap(~ metric, scales = "free") +
  theme(legend.position="none") 

Wilcoxon rank sum tests for various alpha diversity metrics between group

res = adiv %>%
        gather(key = metric, value = value, measures) %>% 
        mutate(metric = factor(metric, levels = measures)) %>%
        group_by(metric) %>%
        rstatix::wilcox_test( value ~ group, p.adjust.method = "none") %>%
        rstatix::add_significance() %>%
        rstatix::add_xy_position(x = "group", scales = "free_y")

res

By time (pre VS post)

adiv %>%
  group_by( time ) %>%
  dplyr::summarise(median_Observed = median(Observed),
                   median_Chao1 = median(Chao1),
                   median_Shannon = median(Shannon),
                   median_Simpson = median(Simpson),
                   median_InvSimpson = median(InvSimpson),
                   median_Fisher = median(Fisher),
                   median_PD = median(PD))

# A tibble: 2 × 8
  time  median_Observed median_Chao1 median_Shannon median_Simpson
  <chr>           <dbl>        <dbl>          <dbl>          <dbl>
1 post            1400.        1833.           5.44          0.987
2 pre             1452.        1855.           5.49          0.987
# ℹ 3 more variables: median_InvSimpson <dbl>, median_Fisher <dbl>,
#   median_PD <dbl>

#Plot adiv measures
adiv %>%
  gather(key = metric, value = value, measures) %>% #c("Observed", "Shannon", "Chao1", "Simpson", "PD")
  mutate(metric = factor(metric, levels = measures)) %>% #c("Observed", "Shannon", "Chao1", "Simpson", "PD")
  ggplot(aes(x = time, y = value)) +
  geom_boxplot(outlier.color = NA) +
  geom_jitter(aes(color = time), height = 0, width = .2) +
  labs(x = "", y = "Alpha Diversity Measures") +
  facet_wrap(~ metric, scales = "free") +
  theme(legend.position="none") 

Wilcoxon rank sum tests for various alpha diversity metrics between time

res = adiv %>%
        gather(key = metric, value = value, measures) %>% 
        mutate(metric = factor(metric, levels = measures)) %>%
        group_by(metric) %>%
        rstatix::wilcox_test( value ~ time, p.adjust.method = "none") %>%
        rstatix::add_significance() %>%
        rstatix::add_xy_position(x = "time", scales = "free_y")

res

Beta Diversity

Alpha diversity describes the diversity within a sample. Beta diversity, on the other hand, describes the diversity between samples. Beta-diversity provides a measure of similarity, or dissimilarity, of one microbial composition to another. This implies that beta diversity measures can not be calculated for a single sample, but is calculated for all pairs of samples in the dataset.

Beta-diversity is typically calculated on the OTU/ASV/species composition tables directly, but can be calculated using abundances at higher taxonomic levels. Combining beta diversity measures with dimension reduction models such as Principal Coordinate Analysis (PCoA), comprise a very powerful tool for visualizing similarities and disimilarities between samples.

Principal Coordinate Analysis (PCoA)

Step 1, from ASV table to distance matrix

Many ecologically informative measures are also commonly used and include:

• Bray-Curtis similarity
• Jaccard similarity
• Unifrac and weighted Unifrac metrics

Both, Jaccard and Bray Curtis makes a flat comparison of the samples, meaning, that any two different OTU/ASV are alike. This, however we know is not the case, as OTU/ASV belonging to the same taxonomic group (say Family) are more similar than two belonging to different groups. Utilizing this actively in the computation of similarities have some advances, as it removes uncertainty due to somewhat wrong annotation, reduces emphasis on differences between phylogenetic similar OTU/ASV and enhances differences due to differences on a diverse phylogenetic level. Unifrac can be computed focusing on presence/absence as well as abundance. These are called unweigted and weighted unifrac respectively.

# beta diversity should be calculated using relative abundance
ps.beta = ps.clean.re

ordBC <- ordinate(ps.beta, "PCoA", "bray")
#str(ordBC)
ordJC <- ordinate(ps.beta, "PCoA", "jaccard")
ordUF <- ordinate(ps.beta, "PCoA", "unifrac")
ordwUF <- ordinate(ps.beta, "PCoA", "wunifrac")

smpID <- sample_data(ps.beta)$sample

# Keep first 2 vectors (latent variables, PCs) of each distance matrix
df <- rbind(data.frame(ordBC$vectors[,1:2], sample = smpID, method = 'BC'),
            data.frame(ordJC$vectors[,1:2], sample = smpID,method = 'Jaccard'),
            data.frame(ordUF$vectors[,1:2], sample = smpID,method = 'unifrac'),
            data.frame(ordwUF$vectors[,1:2], sample = smpID,method = 'wunifrac'))

# add sample_data info
df <- merge(df, data.frame(sample_data(ps.beta)), by = 'sample')

df

Step 2, from distance matrix to PCoA plot.

Beta diversity measures has a similar structure as correlation/covariance matrices. Such matrices can be projected into a low dimensional representation which can be visualized, showing the largest common variation from complex data in a simple plot. These plots are often referred to as ordination plots.

ggplot(data = df, aes(Axis.1,Axis.2, 
                      color = group, shape=time) ) + 
  geom_point() + 
  stat_ellipse() + 
  facet_wrap(~method,scales = 'free')

PERMANOVA analysis for Beta Diversity

Now we have visualize beta diversity of our samples, most of the graphs above have considerable overlaps. How do we test it? We can use Permutational multivariate analysis of variance (PERMANOVA) to test the effect of group (CNT vs CPB), or time (pre-op vs post-op), using the vegan::adonis2 function. This analysis is ANOVA using distance matrices - for partitioning distance matrices among sources of variation (e.g., our group and time) and fitting linear models to distances matrices. There are no analytical null distribution for PERMANOVA, because this is constructed by random permutation.

D_BC <- phyloseq::distance(ps.beta, "bray")
D_JC <- phyloseq::distance(ps.beta, "jaccard")
D_UF <- phyloseq::distance(ps.beta,  "unifrac")
D_wUF <- phyloseq::distance(ps.beta, "wunifrac")

dist_list = list( "Bray Curtis" = D_BC, 'Jaccard' = D_JC, 'Unifrac' = D_UF, 'weighted Unifrac' = D_wUF)

By group (CPB VS CNT)

We are testing community similarities between CNT and CPB, i.e., comparing how dissimilar the communities are by group.

for (i in 1:length(dist_list)) {
        res_group = vegan::adonis2(dist_list[[i]] ~ phyloseq::sample_data(ps.beta)$group)
        
        cat(names(dist_list[i]))
        print(res_group)
        cat("\n", "--------------------------------------------------------------------------", "\n")
}

Bray CurtisPermutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

vegan::adonis2(formula = dist_list[[i]] ~ phyloseq::sample_data(ps.beta)$group)
                                     Df SumOfSqs      R2     F Pr(>F)  
phyloseq::sample_data(ps.beta)$group  1   0.4231 0.07528 1.791  0.014 *
Residual                             22   5.1970 0.92472               
Total                                23   5.6201 1.00000               
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 -------------------------------------------------------------------------- 
JaccardPermutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

vegan::adonis2(formula = dist_list[[i]] ~ phyloseq::sample_data(ps.beta)$group)
                                     Df SumOfSqs      R2      F Pr(>F)   
phyloseq::sample_data(ps.beta)$group  1   0.5431 0.07076 1.6752  0.004 **
Residual                             22   7.1328 0.92924                 
Total                                23   7.6760 1.00000                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 -------------------------------------------------------------------------- 
UnifracPermutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

vegan::adonis2(formula = dist_list[[i]] ~ phyloseq::sample_data(ps.beta)$group)
                                     Df SumOfSqs      R2      F Pr(>F)  
phyloseq::sample_data(ps.beta)$group  1  0.25493 0.08254 1.9791  0.029 *
Residual                             22  2.83379 0.91746                
Total                                23  3.08873 1.00000                
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 -------------------------------------------------------------------------- 
weighted UnifracPermutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

vegan::adonis2(formula = dist_list[[i]] ~ phyloseq::sample_data(ps.beta)$group)
                                     Df SumOfSqs      R2      F Pr(>F)  
phyloseq::sample_data(ps.beta)$group  1  0.03701 0.10194 2.4973   0.02 *
Residual                             22  0.32603 0.89806                
Total                                23  0.36303 1.00000                
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 -------------------------------------------------------------------------- 

By time (pre VS post)

We are testing community similarities between pre and post, i.e., comparing how dissimilar the communities are by time.

for (i in 1:length(dist_list)) {
        res_time = vegan::adonis2(dist_list[[i]] ~ phyloseq::sample_data(ps.beta)$time)
        
        cat(names(dist_list[i]))
        print(res_time)
        cat("\n", "--------------------------------------------------------------------------", "\n")
}

Bray CurtisPermutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

vegan::adonis2(formula = dist_list[[i]] ~ phyloseq::sample_data(ps.beta)$time)
                                    Df SumOfSqs      R2      F Pr(>F)
phyloseq::sample_data(ps.beta)$time  1   0.0751 0.01336 0.2979  0.998
Residual                            22   5.5450 0.98664              
Total                               23   5.6201 1.00000              

 -------------------------------------------------------------------------- 
JaccardPermutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

vegan::adonis2(formula = dist_list[[i]] ~ phyloseq::sample_data(ps.beta)$time)
                                    Df SumOfSqs      R2      F Pr(>F)
phyloseq::sample_data(ps.beta)$time  1   0.1316 0.01714 0.3837  0.997
Residual                            22   7.5444 0.98286              
Total                               23   7.6760 1.00000              

 -------------------------------------------------------------------------- 
UnifracPermutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

vegan::adonis2(formula = dist_list[[i]] ~ phyloseq::sample_data(ps.beta)$time)
                                    Df SumOfSqs      R2      F Pr(>F)
phyloseq::sample_data(ps.beta)$time  1  0.11204 0.03627 0.8281  0.639
Residual                            22  2.97668 0.96373              
Total                               23  3.08873 1.00000              

 -------------------------------------------------------------------------- 
weighted UnifracPermutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

vegan::adonis2(formula = dist_list[[i]] ~ phyloseq::sample_data(ps.beta)$time)
                                    Df SumOfSqs      R2      F Pr(>F)
phyloseq::sample_data(ps.beta)$time  1  0.00520 0.01433 0.3198  0.977
Residual                            22  0.35783 0.98567              
Total                               23  0.36303 1.00000              

 -------------------------------------------------------------------------- 

References and Further Readings

1. phyloseq website: http://joey711.github.io/phyloseq/index.html
2. phyloseq paper: https://doi.org/10.1371/journal.pone.0061217
3. Xia, Yinglin, Jun Sun, and Ding-Geng Chen. Statistical analysis of microbiome data with R. Springer, 2018.
4. Microbiome Data Analysis
5. 16S rDNA amplicon sequencing analysis using R

Reproducibility

R session information:

sessionInfo()

R version 4.2.1 (2022-06-23 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19045)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.utf8 
[2] LC_CTYPE=English_United States.utf8   
[3] LC_MONETARY=English_United States.utf8
[4] LC_NUMERIC=C                          
[5] LC_TIME=English_United States.utf8    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] picante_1.8.2   nlme_3.1-162    ape_5.7-1       vegan_2.6-4    
 [5] lattice_0.21-8  permute_0.9-7   rstatix_0.7.2   phyloseq_1.40.0
 [9] lubridate_1.9.2 forcats_1.0.0   stringr_1.5.0   dplyr_1.1.2    
[13] purrr_1.0.1     readr_2.1.4     tidyr_1.3.0     tibble_3.2.1   
[17] ggplot2_3.4.2   tidyverse_2.0.0

loaded via a namespace (and not attached):
 [1] bitops_1.0-7           GenomeInfoDb_1.32.4    tools_4.2.1           
 [4] backports_1.4.1        bslib_0.4.2            utf8_1.2.3            
 [7] R6_2.5.1               DT_0.27                DBI_1.1.3             
[10] BiocGenerics_0.42.0    mgcv_1.8-42            colorspace_2.1-0      
[13] rhdf5filters_1.8.0     ade4_1.7-22            withr_2.5.0           
[16] tidyselect_1.2.0       compiler_4.2.1         cli_3.6.1             
[19] Biobase_2.56.0         labeling_0.4.2         sass_0.4.6            
[22] scales_1.2.1           digest_0.6.31          rmarkdown_2.21        
[25] XVector_0.36.0         pkgconfig_2.0.3        htmltools_0.5.5       
[28] fastmap_1.1.1          htmlwidgets_1.6.2      rlang_1.1.1           
[31] rstudioapi_0.14        farver_2.1.1           jquerylib_0.1.4       
[34] generics_0.1.3         jsonlite_1.8.4         crosstalk_1.2.0       
[37] car_3.1-2              RCurl_1.98-1.12        magrittr_2.0.3        
[40] GenomeInfoDbData_1.2.8 biomformat_1.24.0      Matrix_1.5-4          
[43] Rcpp_1.0.10            munsell_0.5.0          S4Vectors_0.34.0      
[46] Rhdf5lib_1.18.2        fansi_1.0.4            abind_1.4-5           
[49] lifecycle_1.0.3        stringi_1.7.12         yaml_2.3.7            
[52] carData_3.0-5          MASS_7.3-60            zlibbioc_1.42.0       
[55] rhdf5_2.40.0           plyr_1.8.8             grid_4.2.1            
[58] parallel_4.2.1         crayon_1.5.2           Biostrings_2.64.1     
[61] splines_4.2.1          multtest_2.52.0        hms_1.1.3             
[64] knitr_1.42             pillar_1.9.0           igraph_1.4.2          
[67] reshape2_1.4.4         codetools_0.2-19       stats4_4.2.1          
[70] glue_1.6.2             evaluate_0.21          data.table_1.14.8     
[73] vctrs_0.6.2            tzdb_0.4.0             foreach_1.5.2         
[76] gtable_0.3.3           cachem_1.0.8           xfun_0.39             
[79] broom_1.0.4            survival_3.5-5         iterators_1.0.14      
[82] IRanges_2.30.1         cluster_2.1.4          timechange_0.2.0      
[85] ellipsis_0.3.2