Contents

1 Introduction

This vignette demonstrates how to use the package sSNAPPY to compute directional single sample pathway perturbation scores by incorporating pathway topologies, utilize sample permutation to test the significance of individual scores and compare average pathway activities across treatments.

The package also provides a function to visualise overlap between pathway genes contained in perturbed biological pathways as network plots.

2 To get ready

2.1 Installation

The package sSNAPPY can be installed using the package BiocManager

if (!"BiocManager" %in% rownames(installed.packages()))
  install.packages("BiocManager")
# Other packages required in this vignette
pkg <- c("tidyverse", "magrittr", "ggplot2", "cowplot", "DT")
BiocManager::install(pkg)
BiocManager::install("sSNAPPY")
install.packages("htmltools")

2.2 Load packages

library(sSNAPPY)
library(tidyverse)
library(magrittr)
library(ggplot2)
library(cowplot)
library(DT)
library(htmltools)

2.3 Load data

The example dataset used for this tutorial can be loaded with data() as shown below. It’s a subset of data retrieved from Singhal H et al. 2016, where ER-positive primary breast cancer tumour tissues collected from 12 patients were split into tissue fragments of equal sizes for different treatments.

For this tutorial, we are only looking at the RNA-seq data from samples that were treated with vehicle, R5020(progesterone) OR E2(estrogen) + R5020 for 48 hrs. They were from 5 different patients, giving rise to 15 samples in total. A more detailed description of the dataset can be assessed through the help page (?logCPM_example and ?metadata_example).

data(logCPM_example)
data(metadata_example)
# check if samples included in the logCPM matrix and metadata dataframe are identical
setequal(colnames(logCPM_example), metadata_example$sample)
## [1] TRUE
# View sample metadata
datatable(metadata_example,  filter = "top")

3 sSNAPPY workflow

3.1 Compute weighted single-sample logFCs (ssLogFCs)

It is expected that the logCPM matrix will be filtered to remove undetectable genes and normalised to correct for library sizes or other systematic artefacts, such as gene length or GC contents, prior to applying the sSNAPPY workflow. Filtration and normalisation have been performed on the example dataset.

Before single-sample logFCs (ssLogFCs) can be computed, row names of the logCPM matrix need to be converted to entrez ID. This is because all the pathway topology information retrieved will be in entrez ID. The conversion can be achieved through bioconductor packages AnnotationHub and ensembldb.

head(logCPM_example)

To compute the ssLogFCs, samples must be in matching pairs. In our example data, treated samples were matched to the corresponding control samples derived from the same patients. Therefore the factor parameter of the weight_ss_fc() functions will be set to be “patient”. The function also requires the control treatment level to be specified, which was “Vehicle” in this case.

weight_ss_fc() requires both the logCPM matrix and sample metadata as input. The column names of the logCPM matrix should be sample names, matching a column called sample in the metadata. The metadata must also contain a treatment column, and a column corresponding to the factor parameter (i.e.patient in this case).

#compute weighted single sample logFCs
weightedFC <- weight_ss_fc(logCPM_example, metadata = metadata_example,
factor = "patient", control = "Vehicle")

The weight_ss_fc() function firstly computes raw ssLogFCs for each gene by subtracting logCPM values of the control sample from the logCPM values of treated samples for each patient.

It has been demonstrated previously that in RNA-seq data, lowly expressed genes turn to have a larger variance, which is also demonstrated by the plots below. To reduce the impact of this artefact, weight_ss_fc also weights each ssLogFCs by estimating the relationship between the variance in ssLogFCs and mean logCPM, and defining the gene-wise weight to be the inverse of the predicted variance of that gene’s mean logCPM value.

perSample_FC <- lapply(levels(metadata_example$patient), function(x){
    temp <- logCPM_example[seq_len(1000),str_detect(colnames(logCPM_example), x)] 
    ratio <- temp[, str_detect(colnames(temp), "Vehicle", negate = TRUE)] - temp[, str_detect(colnames(temp), "Vehicle")] 
}) %>%
    do.call(cbind,.)
aveCPM <- logCPM_example[seq_len(1000),] %>%
    rowMeans() %>%
    enframe(name = "gene_id", 
            value = "aveCPM")
p1 <- perSample_FC %>%
    as.data.frame() %>%
    rownames_to_column("gene_id") %>%
    pivot_longer(cols = -"gene_id",
                 names_to = "name",
                 values_to = "logFC") %>%
    left_join(aveCPM) %>%
    ggplot(aes(aveCPM, logFC)) +
    geom_point() +
    labs(y = "sslogFC", 
         x = "Average logCPM") +
    theme(
        panel.background = element_blank()
    )
p2 <- data.frame(
    gene_id = rownames(perSample_FC),
    variance = perSample_FC %>%
        apply(1,var)) %>%
    left_join(aveCPM) %>%
    ggplot(aes(aveCPM, variance)) +
    geom_point() +
    geom_smooth(method = "loess") +
    labs(y = "Variance in ssLogFCs", 
         x = "Average logCPM") +
    theme(
        panel.background = element_blank()
    )
plot_grid(p1, p2)

The output of the weight_ss_fc() function is a list with two element, where one is the weighted ssLogFCs matrix and the other is a vector of gene-wise weights.

3.2 Retrieve pathway topologies in the required format

sSNAPPY adopts the pathway perturbation scoring algorithm proposed in SPIA, which makes use of gene-set topologies and gene-gene interaction to propagate pathway genes’ logFCs down the topologies to compute pathway perturbation scores, where signs of scores reflect the potential directions of changes.

Therefore, pathway topology information needs to be firstly retrieved from your chosen database and converted to weight adjacency matrices, the format required to apply the scoring algorithm.

This step is achieved through a chain of functions that are part of the grapghite, which have been nested into one simple function in this package called retrieve_topology(). Databases that are currently supported are:

library(graphite)
graphite::pathwayDatabases() %>%
  dplyr::filter(species ==  "hsapiens") %>%
  pander::pander()
species database
hsapiens kegg
hsapiens panther
hsapiens pathbank
hsapiens pharmgkb
hsapiens reactome
hsapiens smpdb
hsapiens wikipathways

The retrieved topology information will be a list where each element corresponds a pathway. It’s recommended to save the list as a file so this step only needs to be performed once for each database.

This vignette chose KEGG pathways as an example.

gsTopology <- retrieve_topology(database = "kegg")
head(names(gsTopology))
## [1] "Glycolysis / Gluconeogenesis"            
## [2] "Citrate cycle (TCA cycle)"               
## [3] "Pentose phosphate pathway"               
## [4] "Pentose and glucuronate interconversions"
## [5] "Fructose and mannose metabolism"         
## [6] "Galactose metabolism"

If only selected pathways are of interest to your research, it’s possible to only retrieve the topologies of those pathways by specifying the pathway names.

gsTopology_sub <- retrieve_topology(
  database = "kegg",
  pathwayName = c(
    "Glycolysis / Gluconeogenesis", 
    "Citrate cycle (TCA cycle)",
    "Pentose phosphate pathway"
    ))
names(gsTopology_sub)
## [1] "Glycolysis / Gluconeogenesis" "Citrate cycle (TCA cycle)"   
## [3] "Pentose phosphate pathway"

3.3 Score single sample pathway perturbation

Once the expression matrix, sample metadata and pathway topologies are all ready, gene-wise single-sample perturbation scores can be computed:

genePertScore <- raw_gene_pert(weightedFC$logFC, gsTopology)

and summed to derive pathway perturbation scores for each pathway in each treated samples.

pathwayPertScore <- pathway_pert(genePertScore)
head(pathwayPertScore)
##           sample           tA                                   gs_name
## 1 E2+R5020_N2_48  0.006717358 EGFR tyrosine kinase inhibitor resistance
## 2    R5020_N2_48  0.006718738 EGFR tyrosine kinase inhibitor resistance
## 3 E2+R5020_N3_48  0.002785146 EGFR tyrosine kinase inhibitor resistance
## 4    R5020_N3_48 -0.004030015 EGFR tyrosine kinase inhibitor resistance
## 5 E2+R5020_P4_48  0.004087122 EGFR tyrosine kinase inhibitor resistance
## 6    R5020_P4_48  0.015609886 EGFR tyrosine kinase inhibitor resistance

3.4 Generate null distributions of perturbation scores

To derive the empirical p-values for each single-sample perturbation scores and normalize the raw scores for comparing overall treatment effects, null distributions of scores for each pathway are generated through a sample-label permutation approach.

For each round of permutation, sample labels are randomly shuffled to derive the permuted ssLogFCs, which are then used to score pathway perturbation. We recommend performing a minimum of 1000 rounds of permutation, which means at least 8 samples are required.

The generate_permuted_scores() function does not require sample metadata but the number of treatments in the study design, including the control treatment, need to be specified. In this example data, the number of treatmentS was 3 (Vehicle, E2, and E2+R5020).

The output of the generate_permuted_scores() function is a list where each element is a vector of permuted perturbation scores for a pathway.

The permutation step relies on the parallel computing feature provided by BiocParallel. You can choose to customize the parallel back-end or stick with the default one returned by BiocParallel::bpparam(). Depending on the size of the data, this step can take some time to complete. If the sample size is large, we recommend users consider performing this step on an HPC.

permutedScore <- generate_permuted_scores(
  expreMatrix  = logCPM_example, 
  numOfTreat = 3, NB = 1000, 
  gsTopology = gsTopology, 
  weight = weightedFC$weight
)

3.5 Test significant perturbation on

3.5.1 single-sample level

After the empirical null distributions are generated, the median and mad of each distribution will be calculated for each pathway to convert the test single-sample perturbation scores derived from the compute_perturbation_score() function to robust z-scores: \[(Score - Median)/MAD\]

Two-sided p-values associated with each robust z-scores are also computed and will be corrected for multiple-testing using a user-defined approach. The default is fdr.

The normalise_by_permu() function requires the test perturbation scores and permuted perturbation scores as input.

normalisedScores <- normalise_by_permu(permutedScore, ssPertScore)

Since the permutation step takes a long time to run and the output is too large to be included as part of the package, the results of the normalise_by_permu step has been pre-computed and can be loaded with:

load(system.file("extdata", "normalisedScores.rda", package = "sSNAPPY"))

Pathways that were significant perturbed within individual samples are:

normalisedScores %>%
    dplyr::filter(adjPvalue < 0.05) %>%
    left_join(metadata_example) %>%
    mutate_at(vars(c("sample", "gs_name")), as.factor) %>%
    mutate_if(is.numeric, sprintf, fmt = '%#.4f') %>%
    mutate(Direction = ifelse(robustZ < 0, "Inhibited", "Activation")) %>%
    dplyr::select(
        sample, patient, Treatment = treatment,
        `Perturbation Score` = robustZ, Direction,
        `Gene-set name` = gs_name, 
        `P-value` = pvalue, 
        FDR = adjPvalue
    ) %>%
    datatable(
        filter = "top", 
        options = list(
            columnDefs = list(list(targets = "Direction", visible = FALSE))
        ), 
        caption = htmltools::tags$caption(
                  htmltools::em(
                      "Pathways that were significant perturbed within individual samples.")
              )
    ) %>% 
    formatStyle(
        'Perturbation Score', 'Direction',
        color = styleEqual(c("Inhibited", "Activation"), c("blue", "red"))
    )

3.5.1.1 Visualise overlap between gene-sets as networks

To visualise significantly perturbed biological pathways as a network, where edges between gene-sets reflect how much overlap those two gene-sets share, we can use the plot_gs_network function. The function can take normalise_by_permu()’s output, or a subset of it as its direct input.

Nodes in the network plots could be colored by the predicted direction of perturbation (i.e.. sign of robust z-score):

pl <- normalisedScores %>%
  dplyr::filter(adjPvalue < 0.05) %>%
  split(f = .$sample) %>%
  lapply(
    plot_gs_network, 
    # layout = "dh",
    gsTopology = gsTopology, 
    colorBy = "robustZ"
  ) %>%
    lapply(function(x){
        x + theme(
        panel.grid = element_blank(), 
        panel.background = element_blank()
        ) })
plot_grid(
    plotlist = pl, 
    labels = names(pl),
    label_size = 8,
    nrow = 1)
*Pathways significantly perturbed in individual samples, where gene-sets were colored by pathways' directions of changes.*

(#fig:sigGS_nt_zscore)Pathways significantly perturbed in individual samples, where gene-sets were colored by pathways’ directions of changes.

Or p-values:

pl <- normalisedScores %>%
    dplyr::filter(adjPvalue < 0.05) %>%
    split(f = .$sample) %>%
    lapply(
        plot_gs_network, 
        gsTopology = gsTopology, 
        colorBy = "pvalue", 
        color_lg_title = "P-value"
    ) %>%
    lapply(function(x){
        x + theme(
        panel.grid = element_blank(), 
        panel.background = element_blank()
        ) })
plot_grid(
    plotlist = pl, 
    labels = names(pl),
    label_size = 8,
    nrow = 1)
*Pathways significantly perturbed in individual samples, where gene-sets were colored by pathways' p-values*

(#fig:sigGS_nt_pvalue)Pathways significantly perturbed in individual samples, where gene-sets were colored by pathways’ p-values

The function allows you to customize the layout, colour, edge transparency and other aesthetics of the graph. More information can be found on the help page (?plot_gs_network). The output of the graph is a ggplot object and the theme of it can be changed just as any other ggplot figures.

3.5.2 treatment-level

In addition to testing pathway perturbations at single-sample level, normalised perturbation scores can also be used to model mean treatment effects within a group, such as within each treatment group. An advantage of this method is that it has a high level of flexibility that allows us to incorporate confounding factors as co-factors or co-variates to offset their effects.

For example, in the example data-set, samples were collected from patients with different progesterone receptor (PR) statuses. Knowing that PR status would affect how tumour tissues responded to estrogen and/or progesterone treatments, we can include it as a cofactor to offset its confounding effect.

fit <- normalisedScores %>%
    left_join(metadata_example) %>%
    split(f = .$gs_name) %>%
    #.["Estrogen signaling pathway"] %>%
    lapply(function(x)lm(robustZ ~ 0 + treatment + PR, data = x)) %>%
    lapply(summary)
treat_sig <- lapply(
  names(fit), 
  function(x){
    fit[[x]]$coefficients %>%
      as.data.frame() %>%
      .[seq_len(2),] %>%
      dplyr::select(Estimate, pvalue = `Pr(>|t|)` ) %>%
      rownames_to_column("Treatment") %>%
      mutate(
        gs_name = x, 
        FDR = p.adjust(pvalue, "fdr"), 
        Treatment = str_remove_all(Treatment, "treatment")
      ) 
  }) %>%
  bind_rows()

Results from the linear modelling revealed pathways that were on average perturbed due to each treatment:

treat_sig %>% 
    dplyr::filter(FDR < 0.05) %>%
    mutate_at(vars(c("Treatment", "gs_name")), as.factor) %>%
    mutate_if(is.numeric, sprintf, fmt = '%#.4f') %>%
    mutate(Direction = ifelse(Estimate < 0, "Inhibited", "Activation")) %>%
    dplyr::select(
        Treatment, `Perturbation Score` = Estimate, Direction,
        `Gene-set name` = gs_name, 
        `P-value` = pvalue, 
        FDR
    ) %>%
    datatable(
        filter = "top", 
        options = list(
            columnDefs = list(list(targets = "Direction", visible = FALSE))
        ), 
        caption = htmltools::tags$caption(
                  htmltools::em(
                      "Pathways that were significant perturbed within each treatment group.")
              )
    ) %>% 
    formatStyle(
        'Perturbation Score', 'Direction',
        color = styleEqual(c("Inhibited", "Activation"), c("blue", "red"))
    )

It is not surprising to see that the estrogen signalling pathway was significantly activated under both R5020 and E2+R5020 treatments.

3.5.2.1 Visualise genes’ contributions to pathway perturbation

If there’s a specific pathway that we would like to dig deeper into and explore which pathway genes played a key role in its perturbation, for example, activation of the “Ovarian steroidogenesis”, we can plot the gene-level perturbation scores of all genes contained in that pathway as a heatmap and annotate each column (ie. each sample) by the direction of pathway perturbation in that sample or any other sample metadata using the plot_gene_contribution() function.

From the heatmap below that we can see that the activation of this pathway was driven by a few dominating genes.

plot_gene_contribution(
    genePertScore = genePertScore, gsToPlot = "Ovarian steroidogenesis", metadata = metadata_example, 
    annotation_attribute = c("pathwayPertScore", "treatment"), pathwayPertScore = pathwayPertScore, 
    breaks = seq(-0.0001, 0.0001, length.out = 100),
    show_rownames = FALSE
)
*Gene-level perturbation scores of all genes in the Ovarian steroidogenesis gene-set. The activation of pathway was driven by a few dominating genes.*

Figure 1: Gene-level perturbation scores of all genes in the Ovarian steroidogenesis gene-set
The activation of pathway was driven by a few dominating genes.

By default, genes’ entrez IDs are used and plotted as row names, which may not be very informative. So the row names could be overwritten by providing a data.frame mapping entrez IDs to other identifiers through the mapRownameTo parameter.

Mapping between different gene identifiers could be achieved through the mapIDs() function from the Bioconductor package AnnotationDbi. But to reduce the compiling time of this vignette, mappings between entrez IDs and gene names have been provided as a data.frame called entrez2name.

Note that if the mapping information was provided and the mapping was successful for some genes but not the others, only genes that have been mapped successfully will be plotted.

load(system.file("extdata", "entrez2name.rda", package = "sSNAPPY"))
plot_gene_contribution(
    genePertScore = genePertScore, gsToPlot = "Ovarian steroidogenesis", metadata = metadata_example, 
    annotation_attribute = c("pathwayPertScore", "treatment"), pathwayPertScore = pathwayPertScore,
    breaks = seq(-0.0001, 0.0001, length.out = 100), 
    fontsize_row = 6,
    mapEntrezID = entrez2name
)
*Gene-level perturbation scores of all genes in the Ovarian steroidogenesis gene-set.Genes playing the most important roles in pathway activation was INS, CGA, GNAS.*

Figure 2: Gene-level perturbation scores of all genes in the Ovarian steroidogenesis gene-set.Genes playing the most important roles in pathway activation was INS, CGA, GNAS.

3.5.2.2 Visualise overlap between gene-sets as networks

Results of group-level perturbation can also be visualised using the plot_gs_network() function. We just need to change the name of the Estimate column we had in the treat_sig object to robustZ to colour the networks by the predicted directions of changes.

treat_sig %>% 
    dplyr::filter(FDR < 0.05, Treatment == "E2+R5020") %>%
    dplyr::rename(robustZ = Estimate) %>%
    plot_gs_network(
       
        gsTopology = gsTopology, 
        colorBy = "robustZ"
    ) +
    theme(
        panel.grid = element_blank(), 
        panel.background = element_blank()
    )
*Pathways significantly perturbed by the E2+R5020 combintation treatment, where colors of nodes reflect pathways' directions of changes.*

Figure 3: Pathways significantly perturbed by the E2+R5020 combintation treatment, where colors of nodes reflect pathways’ directions of changes.

By default, plot_gs_network() function does not include nodes that are not connected to any other nodes, which could be turnt off by setting the plotIsolated parameter to TURE.

3.5.2.3 Visualise community structure in the gene-set network

When a large number of pathways were perturbed, it is hard to answer the question “What key biological processes were perturbed?” solely by looking at all the pathway names. To solve that, we can use the plot_community() function to apply a community detection algorithm to the network we had above, and annotate each community by the biological process that most pathways assigned to that community belong to.

treat_sig %>% 
    dplyr::filter(FDR < 0.05, Treatment == "E2+R5020") %>%
    dplyr::rename(robustZ = Estimate) %>%
    plot_community(
        gsTopology = gsTopology, 
        colorBy = "community", 
        color_lg_title = "Community"
    ) +
    theme(panel.background = element_blank())
*Pathways significantly perturbed by the E2+R5020 combintation treatment, annotated by the main biological processes each pathways belong to.*

Figure 4: Pathways significantly perturbed by the E2+R5020 combintation treatment, annotated by the main biological processes each pathways belong to.

sSNAPPY was built in with categorisations of KEGG pathways. But to annotate pathways retrieved from other databases, customised annotation data.frame could be provided through the gsAnnotation parameter.

3.5.2.4 Visualise genes included in perturbed pathways networks

If we want to not only know if two pathways are connected but also the genes connecting those pathways, we can use the plot_gs2gene() function instead:

treat_sig %>% 
    dplyr::filter(FDR < 0.05, Treatment == "E2+R5020") %>%
    dplyr::rename(robustZ = Estimate) %>%
    plot_gs2gene(
        gsTopology = gsTopology, 
        colorGS_By = "robustZ", 
        label_Gene = FALSE,
        GeneNode_size = 1, 
        edgeAlpha = 0.1
    ) +
    theme(panel.background = element_blank())
*Pathways significantly perturbed by the E2+R5020 combintation treatment and all genes contained in those pathways.*

Figure 5: Pathways significantly perturbed by the E2+R5020 combintation treatment and all genes contained in those pathways.

However, since there are a large number of genes in each pathway, the plot above was quite messy and it was unrealistic to plot all genes’ names. So it is recommend to filter pathway genes by their perturbation scores or logFCs.

For example, we can rank genes by the absolute values of their mean single-sample logFCs and only focus on genes that were ranked in the top 200 of the list.

meanFC <- weightedFC$logFC %>%
    .[, str_detect(colnames(.), "E2", negate = TRUE)] %>%
    apply(1, mean )
top200_gene <- meanFC %>%
    abs() %>%
    sort(decreasing = TRUE, ) %>%
    .[1:200] %>%
    names()
top200_FC <- meanFC %>%
    .[names(.) %in% top200_gene]

When we provide genes’ logFCs as a named vector through the geneFC parameter, only pathway genes with logFCs provided will be plotted and gene nodes will be colored by genes’ directions of changes. The names of the logFC vector must be genes’ entrez IDs in the format of “ENTREZID:XXXX”, as pathway topology matrices retrieved through retrieve_topology() always use entrez IDs as identifiers.

However, it is not going to be informative to label genes with their entrez IDs. So the same as the plot_gene_contribution() function, we can provide a data.frame to match genes’ entrez IDs to our chosen identifiers through the mapEntrezID parameter in the plot_gs2gene() function too.

treat_sig %>% 
    dplyr::filter(FDR < 0.05, Treatment == "E2+R5020") %>%
    dplyr::rename(robustZ = Estimate) %>%
    plot_gs2gene(
        gsTopology = gsTopology, 
        colorGS_By = "robustZ", 
        mapEntrezID = entrez2name, 
        geneFC = top200_FC, 
        edgeAlpha = 0.3, 
        GsName_size = 4
    ) +
    theme(panel.background = element_blank())
*Pathways significantly perturbed by the E2+R5020 combintation treatment, and pathway genes with top 200 magnitudes of changes among all E2+R5020-treated samples. Both pathways and genes were colored by their directions of changes.*

Figure 6: Pathways significantly perturbed by the E2+R5020 combintation treatment, and pathway genes with top 200 magnitudes of changes among all E2+R5020-treated samples
Both pathways and genes were colored by their directions of changes.

We can also filters genes’ by their contributions to a pathway perturbations.To help with that, sSNAPPYT provides an option to rank genes’ perturbation scores within each sample for each pathway.

If in a given pathway, both positive and negative gene-wise perturbation scores exist, positive and negative scores are ranked separately, where the larger a positive rank, the more the gene contributed to the pathway’s activation, and the smaller a negative rank, the more the gene contributed to the pathways’ inhibition.

geneRank <- rank_gene_pert(genePertScore, gsTopology)

Depending on the biological question, it could be interesting to plot only the pathway genes with the same directions of changes as the pathway they belonged to, and ignoring those genes that were antagonizing the pathway perturbation.

4 References

5 Session Info

sessionInfo()
## R version 4.2.3 (2023-03-15)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.6 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.16-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.16-bioc/R/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB              LC_COLLATE=C              
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] graphite_1.44.0  htmltools_0.5.5  DT_0.27          cowplot_1.1.1   
##  [5] magrittr_2.0.3   lubridate_1.9.2  forcats_1.0.0    stringr_1.5.0   
##  [9] dplyr_1.1.1      purrr_1.0.1      readr_2.1.4      tidyr_1.3.0     
## [13] tibble_3.2.1     ggplot2_3.4.1    tidyverse_2.0.0  sSNAPPY_1.2.5   
## [17] BiocStyle_2.26.0
## 
## loaded via a namespace (and not attached):
##   [1] ggnewscale_0.4.8            colorspace_2.1-0           
##   [3] ellipsis_0.3.2              XVector_0.38.0             
##   [5] GenomicRanges_1.50.2        farver_2.1.1               
##   [7] graphlayouts_0.8.4          ggrepel_0.9.3              
##   [9] bit64_4.0.5                 AnnotationDbi_1.60.2       
##  [11] fansi_1.0.4                 codetools_0.2-19           
##  [13] splines_4.2.3               cachem_1.0.7               
##  [15] knitr_1.42                  polyclip_1.10-4            
##  [17] jsonlite_1.8.4              png_0.1-8                  
##  [19] pheatmap_1.0.12             graph_1.76.0               
##  [21] ggforce_0.4.1               BiocManager_1.30.20        
##  [23] compiler_4.2.3              httr_1.4.5                 
##  [25] Matrix_1.5-3                fastmap_1.1.1              
##  [27] limma_3.54.2                cli_3.6.1                  
##  [29] tweenr_2.0.2                tools_4.2.3                
##  [31] igraph_1.4.1                gtable_0.3.3               
##  [33] glue_1.6.2                  GenomeInfoDbData_1.2.9     
##  [35] reshape2_1.4.4              rappdirs_0.3.3             
##  [37] Rcpp_1.0.10                 Biobase_2.58.0             
##  [39] jquerylib_0.1.4             vctrs_0.6.1                
##  [41] Biostrings_2.66.0           nlme_3.1-162               
##  [43] crosstalk_1.2.0             ggraph_2.1.0               
##  [45] xfun_0.38                   timechange_0.2.0           
##  [47] lifecycle_1.0.3             org.Hs.eg.db_3.16.0        
##  [49] edgeR_3.40.2                zlibbioc_1.44.0            
##  [51] MASS_7.3-58.3               scales_1.2.1               
##  [53] tidygraph_1.2.3             hms_1.1.3                  
##  [55] MatrixGenerics_1.10.0       parallel_4.2.3             
##  [57] SummarizedExperiment_1.28.0 RColorBrewer_1.1-3         
##  [59] yaml_2.3.7                  memoise_2.0.1              
##  [61] gridExtra_2.3               pander_0.6.5               
##  [63] sass_0.4.5                  stringi_1.7.12             
##  [65] RSQLite_2.3.0               highr_0.10                 
##  [67] S4Vectors_0.36.2            BiocGenerics_0.44.0        
##  [69] BiocParallel_1.32.6         GenomeInfoDb_1.34.9        
##  [71] systemfonts_1.0.4           rlang_1.1.0                
##  [73] pkgconfig_2.0.3             matrixStats_0.63.0         
##  [75] bitops_1.0-7                evaluate_0.20              
##  [77] lattice_0.20-45             htmlwidgets_1.6.2          
##  [79] labeling_0.4.2              bit_4.0.5                  
##  [81] tidyselect_1.2.0            plyr_1.8.8                 
##  [83] bookdown_0.33               R6_2.5.1                   
##  [85] IRanges_2.32.0              generics_0.1.3             
##  [87] DelayedArray_0.24.0         DBI_1.1.3                  
##  [89] pillar_1.9.0                withr_2.5.0                
##  [91] mgcv_1.8-42                 KEGGREST_1.38.0            
##  [93] RCurl_1.98-1.12             crayon_1.5.2               
##  [95] utf8_1.2.3                  tzdb_0.3.0                 
##  [97] rmarkdown_2.21              viridis_0.6.2              
##  [99] locfit_1.5-9.7              grid_4.2.3                 
## [101] blob_1.2.4                  digest_0.6.31              
## [103] stats4_4.2.3                munsell_0.5.0              
## [105] viridisLite_0.4.1           bslib_0.4.2