Abstract

The analysis of time-series data is non-trivial as in involves dependencies between data points. Our package helps researchers analyze RNA-seq datasets which include gene expression measurements taken over time. The methods are specifically designed for at datasets with a small number of replicates and time points – a typical case for RNA-seq time course studies. Short time courses are more difficult to analyze, as many statistical methods designed for time-series data might require a minimum number of time points, e.g. functional data analysis (FDA) and goodness of fit methods might be ineffective. Our approach is non-parametric and gives the user a flexibility to incorporate different normalization techniques, and distance metrics. TimeSeriesExperiment is a comprehensive time course analysis toolbox with methods for data visualization, clustering and differential expression analysis. Additionally, the package can perform enrichment analysis for DE genes.

Introduction

We will demonstrate the effectiveness of TimeSeriesExperiment package using the in-house dataset exploring the role of Cop1 in pro-inflammatory response. The experiments were designed to study the induction and repression kinetics of genes in bone marrow derived macrophages (BMDMs).

The dataset includes cells from 6 mice. The cells were divided into two equal groups. For one group the Cop1 gene was in vitro knock out with tamoxifen. All samples where then subject to LPS treatment to induce an inflammatory response. Bulk RNA-seq was performed at 6 time-points: one at time 0 before LPS treatment, then at time 2.5, 4, 6, 9 and 13 hours after LPS was added.

Obtain and process data

First we load the necessary packages.

.packages <- c(
  "edgeR", "viridis", "Biobase", "BiocFileCache", "SummarizedExperiment", 
  "ggplot2", "dplyr", "tidyr", "tibble", "readr", 
  "TimeSeriesExperiment")
sapply(.packages, require, character.only = TRUE)

##                edgeR              viridis              Biobase        BiocFileCache SummarizedExperiment              ggplot2                dplyr                tidyr               tibble 
##                 TRUE                 TRUE                 TRUE                 TRUE                 TRUE                 TRUE                 TRUE                 TRUE                 TRUE 
##                readr TimeSeriesExperiment 
##                 TRUE                 TRUE

theme_set(theme_bw())
theme_update(
  text = element_text(size = 15),
  strip.background = element_rect(fill = "grey70"),
  strip.text = element_text(size = 10))

The dataset we will study is available from GEO repositories under accession number: GSE114762. We can import the processed read counts saved in a supplementary file ‘GSE114762_raw_counts.csv.gz’. Similarly both the phenotype and feature data are also available for download as ‘GSE114762_sample_info.csv.gz’, and ‘GSE114762_gene_data.csv.gz’ files respectively.

# Remote location of the data
urls <- paste0(
    "https://www.ncbi.nlm.nih.gov/geo/download/",
    "?acc=GSE114762&format=file&file=",
    c(
        "GSE114762_raw_counts.csv.gz",
        "GSE114762_gene_data.csv.gz",
        "GSE114762_sample_info.csv.gz"
    )
)

We will save the files to cache using BiocFileCache utilities to avoid unnecessary multiple downloading.

library(BiocFileCache)
bfc <- BiocFileCache(ask = FALSE)

cnts <- read_csv(bfcrpath(rnames = urls[1])) %>%
  remove_rownames() %>%
  as.matrix()
rownames(cnts) <- cnts[,1]
cnts <- cnts[,-1]

gene.data <- read_csv(bfcrpath(rnames = urls[2])) %>%
  as.data.frame() %>%
  remove_rownames() 
rownames(gene.data) <- gene.data[,1]
gene.data <- gene.data[,-1]

pheno.data <- read_csv(bfcrpath(rnames = urls[3])) %>%
  as.data.frame() %>%
  remove_rownames() 
rownames(pheno.data) <- pheno.data[,1]
pheno.data <- pheno.data[,-1]

Building TimeSeriesExperiment object

We can now combine all the data, cnts, gene.data and pheno.data, into a single TimeSeriesExperiment (S4) object. The data will store everything together in a way that is easier to perform further time course data analysis. The most important fields in the object are “assays”, “colData” which will contain information on group, replicate and time associated with each sample. In TimeSeriesExperiment the time variable MUST be numeric.

cop1.te <- TimeSeriesExperiment(
  assays = list(cnts),
  rowData = gene.data,
  colData = pheno.data,
  timepoint = "time",
  replicate = "replicate",
  group = "group"
)
cop1.te

## class: TimeSeriesExperiment 
## dim: 36528 36 
## metadata(0):
## assays(1): raw
## rownames(36528): 100009600 100009609 ... 99929 99982
## rowData names(5): feature symbol size type desc
## colnames(36): SAM24331086 SAM24331087 ... SAM24331084 SAM24331085
## colData names(9): sample group ... timepoint category
## ========== 
## timepoints(36): 0 2.5 ... 9 13
## groups(36): Loxp Loxp ... WT WT
## replicates(36): Loxp_1 Loxp_1 ... WT_3 WT_3
## -----

## Warning in sprintf(fmt, length(vals), lbls): 2 arguments not used by format 'assayCollapsed: NULL
## '

## assayCollapsed: NULL

## Warning in sprintf(fmt, length(vals), lbls): 2 arguments not used by format 'colDataCollapsed: NULL
## '

## colDataCollapsed: NULL

## Warning in sprintf(fmt, length(vals), lbls): 2 arguments not used by format 'timeSeries: NULL
## '

## timeSeries: NULL

## Warning in sprintf(fmt, length(vals), lbls): 2 arguments not used by format 'dimensionReduction: NULL
## '

## dimensionReduction: NULL

## Warning in sprintf(fmt, length(vals), lbls): 2 arguments not used by format 'clusterRows: NULL
## '

## clusterRows: NULL

## Warning in sprintf(fmt, length(vals), lbls): 2 arguments not used by format 'differentialExpression: NULL
## '

## differentialExpression: NULL

Alternatively, we can build TimeSeriesExperiment object from a ExpressionSet or SummarizedExperiment objects.

To show how this is done, we first combine the three data tables into an ExpressionSet objects.

cop1.es <- ExpressionSet(
  as.matrix(cnts),
  phenoData = AnnotatedDataFrame(pheno.data),
  featureData = AnnotatedDataFrame(gene.data))

Then, we can easily convert an ExpressionSet to a TimeSeriesExperiment object using makeTimeSeriesExperimentFromExpressionSet() function, and indicating columns names with relevant information.

(cop1.te <- makeTimeSeriesExperimentFromExpressionSet(
  cop1.es, timepoint = "time", group = "group",
  replicate = "individual"))

## class: TimeSeriesExperiment 
## dim: 36528 36 
## metadata(0):
## assays(1): raw
## rownames(36528): 100009600 100009609 ... 99929 99982
## rowData names(5): feature symbol size type desc
## colnames(36): SAM24331086 SAM24331087 ... SAM24331084 SAM24331085
## colData names(9): sample group ... timepoint category
## ========== 
## timepoints(36): 0 2.5 ... 9 13
## groups(36): Loxp Loxp ... WT WT
## replicates(36): 1 1 ... 3 3
## -----

## Warning in sprintf(fmt, length(vals), lbls): 2 arguments not used by format 'assayCollapsed: NULL
## '

## assayCollapsed: NULL

## Warning in sprintf(fmt, length(vals), lbls): 2 arguments not used by format 'colDataCollapsed: NULL
## '

## colDataCollapsed: NULL

## Warning in sprintf(fmt, length(vals), lbls): 2 arguments not used by format 'timeSeries: NULL
## '

## timeSeries: NULL

## Warning in sprintf(fmt, length(vals), lbls): 2 arguments not used by format 'dimensionReduction: NULL
## '

## dimensionReduction: NULL

## Warning in sprintf(fmt, length(vals), lbls): 2 arguments not used by format 'clusterRows: NULL
## '

## clusterRows: NULL

## Warning in sprintf(fmt, length(vals), lbls): 2 arguments not used by format 'differentialExpression: NULL
## '

## differentialExpression: NULL

Repeating the same for SummarizedExperiment:

cop1.se <- SummarizedExperiment(
  assays = list(counts = cnts),
  rowData = gene.data, colData = pheno.data)

(cop1.te <- makeTimeSeriesExperimentFromSummarizedExperiment(
  cop1.se, time = "time", group = "group",
  replicate = "individual"))

## class: TimeSeriesExperiment 
## dim: 36528 36 
## metadata(0):
## assays(1): raw
## rownames(36528): 100009600 100009609 ... 99929 99982
## rowData names(5): feature symbol size type desc
## colnames(36): SAM24331086 SAM24331087 ... SAM24331084 SAM24331085
## colData names(9): sample group ... timepoint category
## ========== 
## timepoints(36): 0 2.5 ... 9 13
## groups(36): Loxp Loxp ... WT WT
## replicates(36): 1 1 ... 3 3
## -----

## Warning in sprintf(fmt, length(vals), lbls): 2 arguments not used by format 'assayCollapsed: NULL
## '

## assayCollapsed: NULL

## Warning in sprintf(fmt, length(vals), lbls): 2 arguments not used by format 'colDataCollapsed: NULL
## '

## colDataCollapsed: NULL

## Warning in sprintf(fmt, length(vals), lbls): 2 arguments not used by format 'timeSeries: NULL
## '

## timeSeries: NULL

## Warning in sprintf(fmt, length(vals), lbls): 2 arguments not used by format 'dimensionReduction: NULL
## '

## dimensionReduction: NULL

## Warning in sprintf(fmt, length(vals), lbls): 2 arguments not used by format 'clusterRows: NULL
## '

## clusterRows: NULL

## Warning in sprintf(fmt, length(vals), lbls): 2 arguments not used by format 'differentialExpression: NULL
## '

## differentialExpression: NULL

# Remove raw data
rm(cnts, gene.data, pheno.data, cop1.se, cop1.es)

# Compute CPMs 
cop1.te <- normalizeData(cop1.te)

## Normalizing data...

# Fine genes with minimum mean count of 5 at least in one of the two groups
min_mean_cpm <- 5
group_cpm_means <- data.frame(row.names = rownames(cop1.te))
norm.cnts <- assays(cop1.te)$norm
for(g in unique(groups(cop1.te))) {
  g_cnts <- norm.cnts[ , which(groups(cop1.te) == g)]
  group_cpm_means[, g] <- apply(g_cnts, 1, mean)
}
group_cpm_max <- apply(as.vector(group_cpm_means), 1, max)
genes_expressed <- rownames(cop1.te)[group_cpm_max > min_mean_cpm]

# Filter out the noisy genes
cop1.te <- filterFeatures(cop1.te, genes_expressed)
cop1.te

## class: TimeSeriesExperiment 
## dim: 9356 36 
## metadata(0):
## assays(2): raw norm
## rownames(9356): 100017 100019 ... 99929 99982
## rowData names(5): feature symbol size type desc
## colnames(36): SAM24331086 SAM24331087 ... SAM24331084 SAM24331085
## colData names(9): sample group ... timepoint category
## ========== 
## timepoints(36): 0 2.5 ... 9 13
## groups(36): Loxp Loxp ... WT WT
## replicates(36): 1 1 ... 3 3
## -----

## Warning in sprintf(fmt, length(vals), lbls): 2 arguments not used by format 'assayCollapsed: NULL
## '

## assayCollapsed: NULL

## Warning in sprintf(fmt, length(vals), lbls): 2 arguments not used by format 'colDataCollapsed: NULL
## '

## colDataCollapsed: NULL

## Warning in sprintf(fmt, length(vals), lbls): 2 arguments not used by format 'timeSeries: NULL
## '

## timeSeries: NULL

## Warning in sprintf(fmt, length(vals), lbls): 2 arguments not used by format 'dimensionReduction: NULL
## '

## dimensionReduction: NULL

## Warning in sprintf(fmt, length(vals), lbls): 2 arguments not used by format 'clusterRows: NULL
## '

## clusterRows: NULL

## Warning in sprintf(fmt, length(vals), lbls): 2 arguments not used by format 'differentialExpression: NULL
## '

## differentialExpression: NULL

# Collapsed sample data stored in cop1.te@sample.data.collapsed,
# and mean expression values in cop1.te@data.collapsed
cop1.te <- collapseReplicates(cop1.te, FUN = mean)

## Aggregating across replicates...

Time course format

The main focus on TimeSeriesExperiment is to analyze and visualize time-series data efficiently. For this reason, we convert the expression data in a form of a rectangular matrix into a “time-course format” where each row stores a single time series corresponding to a specified combination of group membership and and replicate id (here mouse id). This data wrangling step can be performed with makeTimeSeries() function, the “time-course” will be stored in a slot timeSeries. This slot contains a list containing data.frames tc and tc_collapsed (if assayCollapsed is defined).

Before converting data to “time-course” format, gene transformation should be performed. Transformation is thus allows for a fair gene-to-gene comparison. For example, when clustering genes, one uses pairwise distances to estimate dissimilarities between genes. Since, we are generally more interested in grouping genes based on similarities in their trajectories rather than their absolute expression levels, the read counts must be transformed before computing the distances.

Currently, gene transformation methods available in TimeSeriesExperiment are “scale_feat_sum” (scaling by gene sum) or “var_stab” (variance stabilization). The user can specify a variance stabilization method if “var_stab” is used. VST methods supported are: “log1p” (log plus one), “asinh” (inverse hyperbolic sine) or “deseq” (DESeq2::varianceStabilizingTransformation).

Usually simply scaling by the gene sum, that is normalizing so that the total abundance the same (and equal to 1) for all genes gives good clustering of gene trajectories.

# Before conversion, scales gene expression to an even sum, for a fair
# gene-to-gene comparison.
cop1.te <- makeTimeSeries(
    cop1.te, feature.trans.method = "scale_feat_sum")

## Converting to timeseries format...

## Warning: `as_data_frame()` was deprecated in tibble 2.0.0.
## Please use `as_tibble()` instead.
## The signature and semantics have changed, see `?as_tibble`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

# untransfomed
head(timeSeries(cop1.te, "ts"))

# transfomed
head(timeSeries(cop1.te, "ts_trans"))

# Add lags to time-course data
cop1.te <- addLags(cop1.te, lambda = c(0.5, 0.25))

## Adding lags with coefficients: 0.5 0.25...

head(timeSeries(cop1.te, "ts_trans"))

Data Visualization

In this section we show plotting utilities available in TimeSeriesExperiment. Visualizations are data exploration tools and serve as the first step in our data analysis. In the following subsections we will describe more in details how heatmaps and PCA plots can be generated.

plotHeatmap(cop1.te, num.feat = 100)

PCA

Another way to explore the dataset is through dimensionality reduction. Here we will project the data into a space of principal components. With PCA, you can visualize both samples and features in the same coordinates space with a biplot. Here we will keep these two maps separately, as the visualization can become overcrowded with points which obscure the inherent structure. Even though, we are plotting the feature and sample projections separately, they can be compared side by side to see which groups of features are more correlated with which group of samples.

cop1.te <- runPCA(cop1.te, var.stabilize.method = "asinh")

plotSamplePCA(cop1.te, col.var = "group", size = 5)

## Warning: Use of `pca.scores[[1]]` is discouraged. Use `.data[[1]]` instead.

## Warning: Use of `pca.scores[[2]]` is discouraged. Use `.data[[2]]` instead.

plotSamplePCA(cop1.te, col.var = "timepoint", size = 5)

## Warning: Use of `pca.scores[[1]]` is discouraged. Use `.data[[1]]` instead.

## Warning: Use of `pca.scores[[2]]` is discouraged. Use `.data[[2]]` instead.

plotTimeSeriesPCA(
    cop1.te,
    linecol = c("WT" = "#e31a1c", "Loxp" = "#1f78b4"),
    main = paste("Visualizing gene profiles with PCA"),
    m  = 15, n = 15, col = adjustcolor("grey77", alpha=0.7), 
    cex.main = 3, cex.axis = 2, cex.lab = 2)

Gene Clustering

To cluster the gene trajectories we will use the data stored in timeSeries slot of cop1.te. These are time-series of transformed expression values together with appended time lags which resolving differences between the temporal trends.

We use hierarchical clustering define gene groupings. Either static or dynamic branch cutting (from dynamicTreeCut package) algorithms can be used to assign clusters. Since hierarchical clustering is computationally intensive (with \(O(n^3)\) complexity for standard implementations), we apply it only to a subset of genes. Specifically, we pick n.top.feat with average (over replicates) most variable expression over time in each of selected groups (here we use both wild type and knock-out) to perform clustering. Remaining genes are, then, assigned to a cluster with the closest centroid. An additional advantage of using only a subset of most variable genes for clustering is that, the core genes which exhibit negligible changes over time (which might be the majority of genes) will not much effect on clustering results.

Here we pick 1000 genes for clustering, which is roughly 1/3 of the number of genes after filtering.

params_for_clustering <- list(
  dynamic = TRUE, 
  dynamic_cut_params = list(deepSplit = TRUE, minModuleSize = 150))

cop1.te <- clusterTimeSeries(
  cop1.te, n.top.feat = 3000, groups = "all",
  clust.params = params_for_clustering)

## Averaging timecourses over all 'groups' selected and recomputing lags with coefficients: 0.5 0.25

## Warning: `data_frame()` was deprecated in tibble 1.1.0.
## Please use `tibble()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

We can see the size of each of clusters computed

# See the count of genes in each cluster
cluster_map <- clusterMap(cop1.te)
table(cluster_map$cluster)

## 
##   C1   C2   C3   C4   C5   C6 
## 2873 2418 2201  979  478  407

We can plot the hclust dendrogram obtained from hierarchical clustering performed.

# Plot the hierarchical clustering of genes2cluster
hclust_obj <- clusterAssignment(cop1.te, "hclust")
plot(x = hclust_obj, labels = FALSE, xlab = "genes", sub = "")

Here we plot average (over replicates) gene trajectories grouped into 10 clusters found using the above described approach. The expression profiles for wild type and knock-out are plotted separately, side by side.

plotTimeSeriesClusters(cop1.te, transparency = 0.2) +
    scale_color_manual(values = c("WT" = "#1f78b4", "Loxp" = "#e31a1c")) +
    theme(strip.text = element_text(size = 10)) +
    ylim(NA, 0.55)

Timecourse plots above show genes clustered clearly into groups related to their pattern of response to LPS treatment. We see cluster C1, C2, and C4 generally inhibited. Genes in C6 spike right after LPS was applied. C3 shows moderate secondary response, and C5 exhibit late increase in expression.

Differential Expression Ranking

In the previous section we grouped the genes based on the trajectories recorded for both experimental groups. In this section we will describe how to find specific genes that exhibit different expression patterns between two experimental groups over the time-course of the study.

cop1.te <- timepointDE(cop1.te, timepoints = "all", alpha = 0.05)

## testing timepoint: 0

## testing timepoint: 2.5

## testing timepoint: 4

## testing timepoint: 6

## testing timepoint: 9

## testing timepoint: 13

tmp_de <- differentialExpression(cop1.te, "timepoint_de")
# First 6 DE genes at timepoint = 2.5:
head(tmp_de$`2.5`)

top10_genes <- sapply(tmp_de, function(tab) {
    if(nrow(tab) == 0)
      return(rep(NA, 10))
    top10 <- tab %>% 
      arrange(desc(abs(logFC))) %>%
      .[["symbol"]] %>% .[1:10]
    if(length(top10) < 10)
      top10 <- c(top10, rep(NA, 10 - length(top10)))
    return(top10)
  }) %>% 
  as.data.frame()
top10_genes

de_genes <- lapply(tmp_de, function(x) x$symbol)
de_any_tmp <- unique(unlist(de_genes))
cat("Out of all", nrow(cop1.te), ", there were",
    length(de_any_tmp), "were found differentially expressed at any timepoint.")

## Out of all 9356 , there were 997 were found differentially expressed at any timepoint.

library(UpSetR)
upset(fromList(de_genes),
      sets = rev(c("2.5", "4", "6", "9", "13")), keep.order = TRUE,
      number.angles = 30, #point.size = 3.5, line.size = 1,
      mainbar.y.label = "DE Gene Intersections",
      sets.x.label = "DE Genes Per Timepoint")

cop1.te <- trajectoryDE(cop1.te, dist_method = "euclidean")
de_trajectory_res <- differentialExpression(cop1.te, "trajectory_de")
saveRDS(de_trajectory_res, "de_trajectory_cop1.rds")

head(de_trajectory_res, 20)

# Select top most different genes across two groups
n_genes <- 20
genes_with_de_trajectory <- de_trajectory_res %>%
    filter(pval <= max(0.05, min(pval)), R2 > 0.7) %>%
    arrange(-R2)
(genes_to_plot <- genes_with_de_trajectory$feature[1:n_genes])

##  [1] "26374"     "320782"    "104156"    "27355"     "16763"     "13039"     "17394"     "66455"     "240913"    "219144"    "15370"     "213956"    "15530"     "100038947" "14118"     "17105"    
## [17] "67465"     "18600"     "237256"    "58217"

plotTimeSeries(cop1.te, features = genes_to_plot) +
    scale_color_manual(values = c("WT" =  "#1f78b4", "Loxp" = "#e31a1c"))

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

plotTimeSeries(cop1.te, features = genes_to_plot) +
    scale_color_manual(values = c("WT" =  "#1f78b4", "Loxp" = "#e31a1c"))

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Functional pathways

In this section we will describe procedures to find functional pathways/gene sets corresponding to genes exhibiting differential expression profiles. We use publicly available reference databases to find relevant pathways.

You can select any list of genes for this step. However, the method is intended for sets of genes found differentially expressed using methods provided by TimeSeriesExperiment, discussed in the previous section.

Below, we use genes with DE trajectory, selected in the previous section. There are 372 in the set:

length(genes_with_de_trajectory$feature)

## [1] 372

Multiple functional pathways might be affected by knocking-out Cop1 gene. Therefore, we expect that genes found differential expressed can be parts of distinct pathways. Selected DE genes can be grouped by their cluster membership found earlier.

Below we plot expression profiles of selected DE genes, separated according to their cluster assignment.

plotTimeSeriesClusters(
  cop1.te, transparency = 0.5,
  features = genes_with_de_trajectory$feature) +
  scale_color_manual(values = c("WT" = "#1f78b4", "Loxp" = "#e31a1c"))

We can test sets (clusters) of DE genes individually for pathway enrichment using pathwayEnrichment() which is built on top of goanna() and kegga() functions from limma package. The features argument provided to pathwayEnrichment() must be Entrez Gene IDs. The user must also specify the relevant species. Optionally, the user can modify the filtering parameters applied to enrichment results: fltr_DE, fltr_N, and fltr_P.DE. For details see the documentation ?pathwayEnrichment.

enrich_res <- pathwayEnrichment(
  object = cop1.te,
  features = genes_with_de_trajectory$feature,
  species = "Mm", ontology ="BP")
saveRDS(enrich_res, "enrich_res_cop1.rds")

We can plot enrichment results for genes with differential trajectory in each cluster separately using plotEnrichment(). Here we plot the top n_max = 15 terms for each cluster. Next to enrichment terms, in the brackets, we have included information on the exact number of DE genes in a set (DE) and the size of the reference enrichment term (N). The points are colored by this fraction. The size of the points correspond the enrichment term size (N).

Below we print we print for example a plot for cluster C4 and C5

clst <- "C4"
plotEnrichment(
    enrich = enrich_res[[clst]], n_max = 15) +
  ggtitle(paste0("Cluster , ", clst, " enrichment"))

clst <- "C5"
plotEnrichment(
    enrich = enrich_res[[clst]], n_max = 15) +
  ggtitle(paste0("Cluster , ", clst, " enrichment"))

Session Information

sessionInfo()

## R version 4.1.1 (2021-08-10)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows Server x64 (build 17763)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=C                           LC_CTYPE=English_United States.1252    LC_MONETARY=English_United States.1252 LC_NUMERIC=C                           LC_TIME=English_United States.1252    
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] UpSetR_1.4.0                TimeSeriesExperiment_1.12.0 readr_2.0.2                 tibble_3.1.5                tidyr_1.1.4                 dplyr_1.0.7                 ggplot2_3.3.5              
##  [8] SummarizedExperiment_1.24.0 GenomicRanges_1.46.0        GenomeInfoDb_1.30.0         IRanges_2.28.0              S4Vectors_0.32.0            MatrixGenerics_1.6.0        matrixStats_0.61.0         
## [15] BiocFileCache_2.2.0         dbplyr_2.1.1                Biobase_2.54.0              BiocGenerics_0.40.0         viridis_0.6.2               viridisLite_0.4.0           edgeR_3.36.0               
## [22] limma_3.50.0                rmarkdown_2.11             
## 
## loaded via a namespace (and not attached):
##   [1] backports_1.2.1        circlize_0.4.13        Hmisc_4.6-0            plyr_1.8.6             splines_4.1.1          BiocParallel_1.28.0    digest_0.6.28          foreach_1.5.1         
##   [9] htmltools_0.5.2        magick_2.7.3           fansi_0.5.0            magrittr_2.0.1         checkmate_2.0.0        memoise_2.0.0          cluster_2.1.2          doParallel_1.0.16     
##  [17] tzdb_0.1.2             ComplexHeatmap_2.10.0  Biostrings_2.62.0      annotate_1.72.0        vroom_1.5.5            jpeg_0.1-9             colorspace_2.0-2       blob_1.2.2            
##  [25] rappdirs_0.3.3         xfun_0.27              crayon_1.4.1           RCurl_1.98-1.5         jsonlite_1.7.2         genefilter_1.76.0      survival_3.2-13        iterators_1.0.13      
##  [33] glue_1.4.2             gtable_0.3.0           zlibbioc_1.40.0        XVector_0.34.0         GetoptLong_1.0.5       DelayedArray_0.20.0    shape_1.4.6            scales_1.1.1          
##  [41] DBI_1.1.1              Rcpp_1.0.7             xtable_1.8-4           htmlTable_2.3.0        clue_0.3-60            foreign_0.8-81         bit_4.0.4              proxy_0.4-26          
##  [49] Formula_1.2-4          htmlwidgets_1.5.4      httr_1.4.2             RColorBrewer_1.1-2     ellipsis_0.3.2         pkgconfig_2.0.3        XML_3.99-0.8           farver_2.1.0          
##  [57] nnet_7.3-16            sass_0.4.0             locfit_1.5-9.4         utf8_1.2.2             dynamicTreeCut_1.63-1  tidyselect_1.1.1       labeling_0.4.2         rlang_0.4.12          
##  [65] AnnotationDbi_1.56.0   munsell_0.5.0          tools_4.1.1            cachem_1.0.6           cli_3.0.1              generics_0.1.1         RSQLite_2.2.8          evaluate_0.14         
##  [73] stringr_1.4.0          fastmap_1.1.0          yaml_2.2.1             knitr_1.36             bit64_4.0.5            purrr_0.3.4            KEGGREST_1.34.0        nlme_3.1-153          
##  [81] compiler_4.1.1         rstudioapi_0.13        filelock_1.0.2         curl_4.3.2             png_0.1-7              geneplotter_1.72.0     bslib_0.3.1            stringi_1.7.5         
##  [89] highr_0.9              lattice_0.20-45        Matrix_1.3-4           vegan_2.5-7            permute_0.9-5          vctrs_0.3.8            pillar_1.6.4           lifecycle_1.0.1       
##  [97] jquerylib_0.1.4        GlobalOptions_0.1.2    data.table_1.14.2      bitops_1.0-7           R6_2.5.1               latticeExtra_0.6-29    gridExtra_2.3          codetools_0.2-18      
## [105] MASS_7.3-54            assertthat_0.2.1       DESeq2_1.34.0          rjson_0.2.20           withr_2.4.2            GenomeInfoDbData_1.2.7 mgcv_1.8-38            parallel_4.1.1        
## [113] hms_1.1.1              grid_4.1.1             rpart_4.1-15           base64enc_0.1-3