TiDEomics Application: Macrophage response to six immune stimuli
29 April 2026
Package
TiDEomics 0.99.0
1 Introduction
TiDEomics is an R package designed to streamline analysis of time-course omics data (e.g. proteomics, transcriptomics), compatible with both single group analysis and multiple group comparison, and accommodating data with missing values.
Here we present a full analysis workflow using TiDEomics on a time-course transcriptomics dataset of murine macrophage exposed to six immune stimuli, with samples collected at 0, 2, 4, 6, 8, and 24 hours post-stimulation. The data set is from GSE263759 published in Traxler et al. (2025).
2 Package dependencies
if (!require("remotes", quietly = TRUE)) install.packages("remotes")
remotes::install_github("hte123/TiDEomics")library(TiDEomics)
library(magrittr)
library(SummarizedExperiment)
#> Loading required package: MatrixGenerics
#> Loading required package: matrixStats
#>
#> Attaching package: 'MatrixGenerics'
#> The following objects are masked from 'package:matrixStats':
#>
#> colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
#> colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
#> colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
#> colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
#> colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
#> colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
#> colWeightedMeans, colWeightedMedians, colWeightedSds,
#> colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
#> rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
#> rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
#> rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
#> rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
#> rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
#> rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
#> rowWeightedSds, rowWeightedVars
#> Loading required package: GenomicRanges
#> Loading required package: stats4
#> Loading required package: BiocGenerics
#> Loading required package: generics
#>
#> Attaching package: 'generics'
#> The following objects are masked from 'package:base':
#>
#> as.difftime, as.factor, as.ordered, intersect, is.element, setdiff,
#> setequal, union
#>
#> Attaching package: 'BiocGenerics'
#> The following objects are masked from 'package:stats':
#>
#> IQR, mad, sd, var, xtabs
#> The following objects are masked from 'package:base':
#>
#> anyDuplicated, aperm, append, as.data.frame, basename, cbind,
#> colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
#> get, grep, grepl, is.unsorted, lapply, Map, mapply, match, mget,
#> order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
#> rbind, Reduce, rownames, sapply, saveRDS, table, tapply, unique,
#> unsplit, which.max, which.min
#> Loading required package: S4Vectors
#>
#> Attaching package: 'S4Vectors'
#> The following object is masked from 'package:utils':
#>
#> findMatches
#> The following objects are masked from 'package:base':
#>
#> expand.grid, I, unname
#> Loading required package: IRanges
#>
#> Attaching package: 'IRanges'
#> The following object is masked from 'package:grDevices':
#>
#> windows
#> Loading required package: Seqinfo
#>
#> Attaching package: 'GenomicRanges'
#> The following object is masked from 'package:magrittr':
#>
#> subtract
#> Loading required package: Biobase
#> Welcome to Bioconductor
#>
#> Vignettes contain introductory material; view with
#> 'browseVignettes()'. To cite Bioconductor, see
#> 'citation("Biobase")', and for packages 'citation("pkgname")'.
#>
#> Attaching package: 'Biobase'
#> The following object is masked from 'package:MatrixGenerics':
#>
#> rowMedians
#> The following objects are masked from 'package:matrixStats':
#>
#> anyMissing, rowMedians
library(org.Mm.eg.db)
#> Loading required package: AnnotationDbi
#>
library(GEOquery)
#> Setting options('download.file.method.GEOquery'='auto')
#> Setting options('GEOquery.inmemory.gpl'=FALSE)
library(stringr)
library(dplyr)
#>
#> Attaching package: 'dplyr'
#> The following object is masked from 'package:AnnotationDbi':
#>
#> select
#> The following object is masked from 'package:Biobase':
#>
#> combine
#> The following objects are masked from 'package:GenomicRanges':
#>
#> intersect, setdiff, union
#> The following object is masked from 'package:Seqinfo':
#>
#> intersect
#> The following objects are masked from 'package:IRanges':
#>
#> collapse, desc, intersect, setdiff, slice, union
#> The following objects are masked from 'package:S4Vectors':
#>
#> first, intersect, rename, setdiff, setequal, union
#> The following objects are masked from 'package:BiocGenerics':
#>
#> combine, intersect, setdiff, setequal, union
#> The following object is masked from 'package:generics':
#>
#> explain
#> The following object is masked from 'package:matrixStats':
#>
#> count
#> The following objects are masked from 'package:stats':
#>
#> filter, lag
#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, union
library(ggplot2)3 Input data
GSE263759 dataset is downloaded and processed:
- Ensembl IDs were mapped to symbols, genes with zero counts in all samples were excluded.
- Use time 0 untreated samples for other groups’ time 0.
A subset of the data is stored in data/tutorial_sample_info.rda and data/tutorial_data.rda for use in the Tutorial, as recorded in data-raw/tutorial_input.R.
# download sample information from GEO
geo_data <- getGEO("GSE263759", GSEMatrix = TRUE)
#> Found 1 file(s)
#> GSE263759_series_matrix.txt.gz
geo_sample_info <- pData(geo_data[[1]]) %>%
dplyr::select(
"title", "stimulus:ch1", "time point:ch1",
"bio-replicate:ch1", "batch:ch1"
) %>%
dplyr::rename(
"Group" = "stimulus:ch1",
"Time" = "time point:ch1",
"Replicate" = "bio-replicate:ch1",
"Batch" = "batch:ch1",
"Sample" = "title"
) %>%
mutate(
Time = Time %>% str_replace("h", "") %>% as.numeric(),
Replicate = Replicate %>% str_replace("R", "") %>% as.numeric(),
Batch = Batch %>% as.numeric()
)
row.names(geo_sample_info) <- NULL
geo_sample_info %>%
dplyr::group_by(Group) %>%
dplyr::summarise(n = n())
# download expression matrix from GEO
getGEOSuppFiles("GSE263759", baseDir = "..//..//data-raw", makeDirectory = FALSE)#
#> Using locally cached version of supplementary file(s) GSE263759 found here:
#> ..//..//data-raw/GSE263759_RNA_raw_counts.csv.gzgeo_data_tb_gz <- gzfile("..//..//data-raw//GSE263759_RNA_raw_counts.csv.gz", 'rt')
geo_data_tb <- read.csv(geo_data_tb_gz) %>%
dplyr::rename("Feature" = "gene")
# map ensembl IDs to gene symbols
geo_data_tb_symbol <- geo_data_tb$Feature %>%
clusterProfiler::bitr(fromType = "ENSEMBL", toType = "SYMBOL", OrgDb = org.Mm.eg.db) %>%
merge(geo_data_tb, by.x = "ENSEMBL", by.y = "Feature") %>%
distinct(SYMBOL, .keep_all = TRUE) %>%
dplyr::select(-ENSEMBL) %>%
dplyr::rename("Feature" = "SYMBOL")
#>
#> 'select()' returned 1:many mapping between keys and columns
#> Warning in clusterProfiler::bitr(., fromType = "ENSEMBL", toType = "SYMBOL", :
#> 41.32% of input gene IDs are fail to map...
# remove genes with all zero counts
gene_keep <- geo_data_tb_symbol[, -1] %>%
apply(1, function(x) sum(x)) %>%
`!=`(0) %>%
which()
geo_data_tb_filtered <- geo_data_tb_symbol[gene_keep, ]
# apply time 0 of untreated to all groups
# duplicate the time 0 untreated samples to all other samples' time 0
geo_sample_info_new <- geo_sample_info
geo_data_tb_new <- geo_data_tb_filtered
for (i in unique(geo_sample_info$Group)) {
if (i != "untreated") {
group_t0 <- geo_data_tb_filtered %>%
dplyr::select(starts_with("RNA_untreated_0h"))
group_t0_sample <- geo_sample_info %>%
filter(Group == "untreated" & Time == 0)
colnames(group_t0) <- colnames(group_t0) %>%
str_replace("untreated", i)
geo_data_tb_new <- cbind(geo_data_tb_new, group_t0)
group_t0_sample <- group_t0_sample %>%
mutate(
Sample = Sample %>%
str_replace("RNA_untreated", paste0("RNA_", i)),
Group = i
)
geo_sample_info_new <- rbind(geo_sample_info_new, group_t0_sample)
}
}data_obj <- create_input(
data = geo_data_tb_new,
sample_ann = geo_sample_info_new)
#> Converting 'Group' column to factor. Default order is alphabetical.
#> Converting 'Replicate' column to factor. Default order is numerical.
#> Converting 'Batch' column to factor. Default order is numerical.
# Log-transform the data when not already done
assays(data_obj)[[1]] <- log2(assays(data_obj)[[1]] + 1) # + 1 to avoid log(0)
data_obj#> class: SummarizedExperiment
#> dim: 20150 76
#> metadata(0):
#> assays(1): ''
#> rownames(20150): Gnai3 Cdc45 ... Gm17315 Ube2srt
#> rowData names(0):
#> colnames(76): RNA_Candida_0h_R1_1 RNA_Candida_0h_R2_1 ...
#> RNA_untreated_24h_R1_2 RNA_untreated_24h_R2_2
#> colData names(5): Sample Group Time Replicate BatchThe log2 (raw count + 1) is used as the input, stored in the first assay slot of data_obj.
#> DataFrame with 76 rows and 5 columns
#> Sample Group Time Replicate
#> <character> <factor> <numeric> <factor>
#> RNA_Candida_0h_R1_1 RNA_Candida_0h_R1_1 Candida 0 1
#> RNA_Candida_0h_R2_1 RNA_Candida_0h_R2_1 Candida 0 2
#> RNA_Candida_2h_R1_1 RNA_Candida_2h_R1_1 Candida 2 1
#> RNA_Candida_2h_R2_1 RNA_Candida_2h_R2_1 Candida 2 2
#> RNA_Candida_4h_R1_1 RNA_Candida_4h_R1_1 Candida 4 1
#> ... ... ... ... ...
#> RNA_untreated_0h_R1_1 RNA_untreated_0h_R1_1 untreated 0 1
#> RNA_untreated_0h_R2_1 RNA_untreated_0h_R2_1 untreated 0 2
#> RNA_untreated_8h_R2_2 RNA_untreated_8h_R2_2 untreated 8 2
#> RNA_untreated_24h_R1_2 RNA_untreated_24h_R1_2 untreated 24 1
#> RNA_untreated_24h_R2_2 RNA_untreated_24h_R2_2 untreated 24 2
#> Batch
#> <factor>
#> RNA_Candida_0h_R1_1 1
#> RNA_Candida_0h_R2_1 1
#> RNA_Candida_2h_R1_1 1
#> RNA_Candida_2h_R2_1 1
#> RNA_Candida_4h_R1_1 1
#> ... ...
#> RNA_untreated_0h_R1_1 1
#> RNA_untreated_0h_R2_1 1
#> RNA_untreated_8h_R2_2 2
#> RNA_untreated_24h_R1_2 2
#> RNA_untreated_24h_R2_2 23.1 Set custom color palette
custom_palette <- c(
"untreated" = "#7570b3",
"IFNbeta" = "#e6ab02",
"IFNgamma" = "#a6761d",
"LPS" = "#66a61e",
"Candida" = "#e7298a",
"LCMV" = "#d95f02",
"Listeria" = "#1b9e77"
)
set_custom_palette(custom_palette)
#> Custom palette has been set successfully for groups untreated, IFNbeta, IFNgamma, LPS, Candida, LCMV, Listeria.4 Quality control
Abundance distribution
plot_distribution(data_obj, facet_by = "Group")
#> Picking joint bandwidth of 0.433
#> Picking joint bandwidth of 0.427
#> Picking joint bandwidth of 0.428
#> Picking joint bandwidth of 0.43
#> Picking joint bandwidth of 0.415
#> Picking joint bandwidth of 0.426
#> Picking joint bandwidth of 0.434
5 Data processing
Normalise to time point 0
Split object into one per group
Merge replicates by mean
Merge groups into one object with merged replicates
#> class: SummarizedExperiment
#> dim: 20150 39
#> metadata(0):
#> assays(2): '' ''
#> rownames(20150): 0610005C13Rik 0610009E02Rik ... Zzef1 Zzz3
#> rowData names(0):
#> colnames(39): Candida_0 Candida_2 ... untreated_8 untreated_24
#> colData names(5): Sample Group Time Replicate Batch5.1 Calculate mean and SD for plotting
5.2 Visualise example genes expression
plot_trend(table_mean_sd_orig,
features = c("Isg15", "Cd68", "Adgre1", "Tnf"),
title = "Example genes", ylab = "Log2 (count + 1)"
)
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Warning: Removed 12 rows containing missing values or values outside the scale range
#> (`geom_point()`).
6 Sample relationships
6.1 Correlation between samples
plot_cor_matrix(data_obj,
method = "spearman",
label_rep = TRUE, label_batch = TRUE,
cellwidth = 1.5, cellheight = 1.5
)
6.2 UMAP
using original data not normalised to time 0
6.2.1 Original data with replicates
6.2.1.1 All groups
#> Warning: Using size for a discrete variable is not advised.
6.2.2 Replicates merged by mean
6.2.2.1 All groups
#> Warning: Using size for a discrete variable is not advised.
6.2.2.2 Per group
plot_umap_by_group(data_obj_merged, nrow = 2)
#> Using n_neighbors = 5
#> Using n_neighbors = 5
#> Using n_neighbors = 5
#> Using n_neighbors = 5
#> Using n_neighbors = 5
#> Using n_neighbors = 5
#> Using n_neighbors = 2
#> Warning: failed creating initial embedding; using random embedding instead
6.3 PCA
using original data not normalised to time 0
6.3.1 Original data with replicates
6.3.1.1 All groups
PC <- plot_pca(data_obj,
plot = TRUE,
plot_screeplot = TRUE,
plot_loadings = FALSE,
plot_morepc = TRUE,
circle = FALSE,
morepc = 1:5
)
#> Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
#> ℹ Please use `linewidth` instead.
#> ℹ The deprecated feature was likely used in the PCAtools package.
#> Please report the issue to the authors.
#> This warning is displayed once per session.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.
#> Warning: Using size for a discrete variable is not advised.
#> Scale for colour is already present.
#> Adding another scale for colour, which will replace the existing scale.
#> Scale for colour is already present.
#> Adding another scale for colour, which will replace the existing scale.
#> Scale for colour is already present.
#> Adding another scale for colour, which will replace the existing scale.
#> Scale for colour is already present.
#> Adding another scale for colour, which will replace the existing scale.
#> Scale for colour is already present.
#> Adding another scale for colour, which will replace the existing scale.
#> Scale for colour is already present.
#> Adding another scale for colour, which will replace the existing scale.
#> Scale for colour is already present.
#> Adding another scale for colour, which will replace the existing scale.
#> Scale for colour is already present.
#> Adding another scale for colour, which will replace the existing scale.
#> Scale for colour is already present.
#> Adding another scale for colour, which will replace the existing scale.
#> Scale for colour is already present.
#> Adding another scale for colour, which will replace the existing scale.
plot_pca_3D(PC, pcs = 1:3)
#> Warning: `line.width` does not currently support multiple values.
#> Warning: `line.width` does not currently support multiple values.
#> Warning: `line.width` does not currently support multiple values.
#> Warning: `line.width` does not currently support multiple values.
#> Warning: `line.width` does not currently support multiple values.
#> Warning: `line.width` does not currently support multiple values.
#> Warning: `line.width` does not currently support multiple values.Note: the 3D plot may not display properly in some html, but should work in an interactive R session.
PCAtools::eigencorplot(PC,
metavars = c("Group", "Time", "Replicate", "Batch"),
components = paste0("PC", 1:5),
col = colorRampPalette(c("#3C5488FF", "white", "#E64B35FF"))(100),
colCorval = "black"
)
#> Warning in PCAtools::eigencorplot(PC, metavars = c("Group", "Time",
#> "Replicate", : Group is not numeric - please check the source data as
#> non-numeric variables will be coerced to numeric
#> Warning in PCAtools::eigencorplot(PC, metavars = c("Group", "Time",
#> "Replicate", : Replicate is not numeric - please check the source data as
#> non-numeric variables will be coerced to numeric
#> Warning in PCAtools::eigencorplot(PC, metavars = c("Group", "Time",
#> "Replicate", : Batch is not numeric - please check the source data as
#> non-numeric variables will be coerced to numeric
6.3.1.2 Per group
plot_pca_by_group(data_obj, circle = TRUE, arrow = TRUE, nrow = 2)
#> Warning: Removed 1 row containing missing values or values outside the scale range
#> (`geom_segment()`).
#> Removed 1 row containing missing values or values outside the scale range
#> (`geom_segment()`).
#> Removed 1 row containing missing values or values outside the scale range
#> (`geom_segment()`).
#> Removed 1 row containing missing values or values outside the scale range
#> (`geom_segment()`).
#> Removed 1 row containing missing values or values outside the scale range
#> (`geom_segment()`).
#> Removed 1 row containing missing values or values outside the scale range
#> (`geom_segment()`).
#> Removed 1 row containing missing values or values outside the scale range
#> (`geom_segment()`).
#> Removed 1 row containing missing values or values outside the scale range
#> (`geom_segment()`).
6.3.2 Replicates merged by mean
6.3.2.1 All groups
PC_merged <- plot_pca(data_obj_merged,
plot = TRUE,
plot_screeplot = TRUE,
plot_loadings = FALSE,
plot_morepc = TRUE,
circle = FALSE,
morepc = 1:5
)
#> Warning: Using size for a discrete variable is not advised.
#> Scale for colour is already present.
#> Adding another scale for colour, which will replace the existing scale.
#> Scale for colour is already present.
#> Adding another scale for colour, which will replace the existing scale.
#> Scale for colour is already present.
#> Adding another scale for colour, which will replace the existing scale.
#> Scale for colour is already present.
#> Adding another scale for colour, which will replace the existing scale.
#> Scale for colour is already present.
#> Adding another scale for colour, which will replace the existing scale.
#> Scale for colour is already present.
#> Adding another scale for colour, which will replace the existing scale.
#> Scale for colour is already present.
#> Adding another scale for colour, which will replace the existing scale.
#> Scale for colour is already present.
#> Adding another scale for colour, which will replace the existing scale.
#> Scale for colour is already present.
#> Adding another scale for colour, which will replace the existing scale.
#> Scale for colour is already present.
#> Adding another scale for colour, which will replace the existing scale.
plot_pca_3D(PC_merged, pcs = 1:3)
#> Warning: `line.width` does not currently support multiple values.
#> Warning: `line.width` does not currently support multiple values.
#> Warning: `line.width` does not currently support multiple values.
#> Warning: `line.width` does not currently support multiple values.
#> Warning: `line.width` does not currently support multiple values.
#> Warning: `line.width` does not currently support multiple values.
#> Warning: `line.width` does not currently support multiple values.Note: the 3D plot may not display properly in some html, but should work in an interactive R session.
PCAtools::eigencorplot(PC_merged,
metavars = c("Group", "Time"),
components = paste0("PC", 1:5),
col = colorRampPalette(c("#3C5488FF", "white", "#E64B35FF"))(100),
colCorval = "black"
)
#> Warning in PCAtools::eigencorplot(PC_merged, metavars = c("Group", "Time"), :
#> Group is not numeric - please check the source data as non-numeric variables
#> will be coerced to numeric
6.3.2.2 Per group
plot_pca_by_group_list(data_obj_merged_list, nrow = 2, circle = FALSE)
#> Warning: Removed 1 row containing missing values or values outside the scale range
#> (`geom_segment()`).
#> Removed 1 row containing missing values or values outside the scale range
#> (`geom_segment()`).
#> Removed 1 row containing missing values or values outside the scale range
#> (`geom_segment()`).
#> Removed 1 row containing missing values or values outside the scale range
#> (`geom_segment()`).
#> Removed 1 row containing missing values or values outside the scale range
#> (`geom_segment()`).
#> Removed 1 row containing missing values or values outside the scale range
#> (`geom_segment()`).
#> Removed 1 row containing missing values or values outside the scale range
#> (`geom_segment()`).
#> Removed 1 row containing missing values or values outside the scale range
#> (`geom_segment()`).
7 Top correlated genes with PCs
top correlated (loading) genes with each principal components (PC)
7.1 Original data with replicates
pca_loadings_complete <- PC$loadings %>% as.data.frame()
loadings_list_complete <- list()
top_n <- 30
for (i in 1:3) {
loadings_list_complete[[paste0("PC", i, "_top+")]] <- pca_loadings_complete %>%
arrange(desc(get(paste0("PC", i)))) %>%
head(top_n) %>%
row.names()
loadings_list_complete[[paste0("PC", i, "_top-")]] <- pca_loadings_complete %>%
arrange(get(paste0("PC", i))) %>%
head(top_n) %>%
row.names()
}
for (i in names(loadings_list_complete)) {
plot_trend(table_mean_sd_orig,
features = loadings_list_complete[[i]],
title = paste0(i, " loadings"),
ylab = "Log2 (count + 1)"
)
}
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Warning: Removed 90 rows containing missing values or values outside the scale range
#> (`geom_point()`).
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Warning: Removed 90 rows containing missing values or values outside the scale range
#> (`geom_point()`).
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Warning: Removed 90 rows containing missing values or values outside the scale range
#> (`geom_point()`).
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Warning: Removed 90 rows containing missing values or values outside the scale range
#> (`geom_point()`).
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Warning: Removed 90 rows containing missing values or values outside the scale range
#> (`geom_point()`).
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Warning: Removed 90 rows containing missing values or values outside the scale range
#> (`geom_point()`).
7.1.1 GO enrichment
background_pca <- row.names(pca_loadings_complete)
go_list_pca <- enrichGO_list(
gene_list = loadings_list_complete,
OrgDb = org.Mm.eg.db,
universe = background_pca,
keyType = "SYMBOL",
simplify = FALSE
)
#> GO category not specified. Using all three: BP, MF, CC.
#> Performing GO enrichment for category: BP
#> Processing gene list: PC1_top+
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: PC1_top-
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: PC2_top+
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: PC2_top-
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: PC3_top+
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: PC3_top-
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Performing GO enrichment for category: MF
#> Processing gene list: PC1_top+
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: PC1_top-
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: PC2_top+
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: PC2_top-
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: PC3_top+
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: PC3_top-
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Performing GO enrichment for category: CC
#> Processing gene list: PC1_top+
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: PC1_top-
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: PC2_top+
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: PC2_top-
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: PC3_top+
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: PC3_top-
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Merging GO enrichment results across gene lists for each category.plot_GO(go_list_pca$all,
plot_dotplot = TRUE,
showCategory_dotplot = 5,
label = "top PCA loading genes"
)
#> No significant GO CC terms found.plot_GO(go_list_pca$all,
plot_cnetplot = TRUE,
showCategory_cnetplot = 5,
label = "top PCA loading genes"
)
#> The following groups are missing in the custom palette: PC1_top+, PC1_top-, PC2_top+, PC2_top-, PC3_top+
#> The following groups are missing in the custom palette: PC1_top+, PC2_top-, PC3_top+
#> No significant GO CC terms found.plot_GO(go_list_pca$all,
plot_emapplot = TRUE,
showCategory_emapplot = 5,
label = "top PCA loading genes")
#> Coordinate system already present.
#> ℹ Adding new coordinate system, which will replace the existing one.
#> The following groups are missing in the custom palette: PC1_top+, PC1_top-, PC2_top+, PC2_top-, PC3_top+
#> Coordinate system already present.
#> ℹ Adding new coordinate system, which will replace the existing one.
#> The following groups are missing in the custom palette: PC1_top+, PC2_top-, PC3_top+
#> No significant GO CC terms found.7.2 Replicates merged by mean
pca_loadings_merged <- PC_merged$loadings %>% as.data.frame()
loadings_list_merged <- list()
top_n <- 30
for (i in 1:3) {
loadings_list_merged[[paste0("PC", i, "_top+")]] <- pca_loadings_merged %>%
arrange(desc(get(paste0("PC", i)))) %>%
head(top_n) %>%
row.names()
loadings_list_merged[[paste0("PC", i, "_top-")]] <- pca_loadings_merged %>%
arrange(get(paste0("PC", i))) %>%
head(top_n) %>%
row.names()
}
for (i in names(loadings_list_merged)) {
plot_trend(table_mean_sd_orig,
features = loadings_list_merged[[i]],
title = paste0(i, " loadings (mean of triplicates)"),
ylab = "Log2 (count + 1)"
)
}
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Warning: Removed 90 rows containing missing values or values outside the scale range
#> (`geom_point()`).
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Warning: Removed 90 rows containing missing values or values outside the scale range
#> (`geom_point()`).
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Warning: Removed 90 rows containing missing values or values outside the scale range
#> (`geom_point()`).
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Warning: Removed 90 rows containing missing values or values outside the scale range
#> (`geom_point()`).
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Warning: Removed 90 rows containing missing values or values outside the scale range
#> (`geom_point()`).
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Warning: Removed 90 rows containing missing values or values outside the scale range
#> (`geom_point()`).
7.2.1 GO enrichment
background_pca_merged <- row.names(pca_loadings_merged)
go_list_pca_merged <- enrichGO_list(
gene_list = loadings_list_merged,
OrgDb = org.Mm.eg.db,
universe = background_pca_merged,
keyType = "SYMBOL",
simplify = FALSE,
)
#> GO category not specified. Using all three: BP, MF, CC.
#> Performing GO enrichment for category: BP
#> Processing gene list: PC1_top+
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: PC1_top-
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: PC2_top+
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: PC2_top-
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: PC3_top+
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: PC3_top-
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Performing GO enrichment for category: MF
#> Processing gene list: PC1_top+
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: PC1_top-
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: PC2_top+
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: PC2_top-
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: PC3_top+
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: PC3_top-
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Performing GO enrichment for category: CC
#> Processing gene list: PC1_top+
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: PC1_top-
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: PC2_top+
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: PC2_top-
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: PC3_top+
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: PC3_top-
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Merging GO enrichment results across gene lists for each category.plot_GO(go_list_pca_merged$all,
plot_dotplot = TRUE,
showCategory_dotplot = 5,
label = "top PCA loading genes (mean of replicates)"
)
#> No significant GO CC terms found.plot_GO(go_list_pca_merged$all,
plot_cnetplot = TRUE,
showCategory_cnetplot = 5,
label = "top PCA loading genes (mean of replicates)"
)
#> The following groups are missing in the custom palette: PC1_top+, PC1_top-, PC2_top+, PC2_top-, PC3_top+
#> The following groups are missing in the custom palette: PC1_top+, PC2_top+, PC3_top+
#> No significant GO CC terms found.plot_GO(go_list_pca_merged$all,
plot_emapplot = TRUE,
showCategory_emapplot = 5,
label = "top PCA loading genes (mean of replicates)")
#> Coordinate system already present.
#> ℹ Adding new coordinate system, which will replace the existing one.
#> The following groups are missing in the custom palette: PC1_top+, PC1_top-, PC2_top+, PC2_top-, PC3_top+
#> Coordinate system already present.
#> ℹ Adding new coordinate system, which will replace the existing one.
#> The following groups are missing in the custom palette: PC1_top+, PC2_top+, PC3_top+
#> No significant GO CC terms found.8 Pairwise differential expression (DE)
8.1 Between time points within each group
Default DE criteria: p.adj < 0.05, |log2FC| > 1
DE_between_time_out <- DE_between_time(data_obj, assay = 1, filter = 1)
#> Comparing group Candida time 2 to 0: keeping 20150 of 20150 features (100.0%)
#> Comparing group Candida time 4 to 0: keeping 20150 of 20150 features (100.0%)
#> Comparing group Candida time 6 to 0: keeping 20150 of 20150 features (100.0%)
#> Comparing group Candida time 8 to 0: keeping 20150 of 20150 features (100.0%)
#> Comparing group Candida time 24 to 0: keeping 20150 of 20150 features (100.0%)
#> Comparing group Candida time 4 to 2: keeping 20150 of 20150 features (100.0%)
#> Comparing group Candida time 6 to 2: keeping 20150 of 20150 features (100.0%)
#> Comparing group Candida time 8 to 2: keeping 20150 of 20150 features (100.0%)
#> Comparing group Candida time 24 to 2: keeping 20150 of 20150 features (100.0%)
#> Comparing group Candida time 6 to 4: keeping 20150 of 20150 features (100.0%)
#> Comparing group Candida time 8 to 4: keeping 20150 of 20150 features (100.0%)
#> Comparing group Candida time 24 to 4: keeping 20150 of 20150 features (100.0%)
#> Comparing group Candida time 8 to 6: keeping 20150 of 20150 features (100.0%)
#> Comparing group Candida time 24 to 6: keeping 20150 of 20150 features (100.0%)
#> Comparing group Candida time 24 to 8: keeping 20150 of 20150 features (100.0%)
#> Comparing group IFNbeta time 2 to 0: keeping 20150 of 20150 features (100.0%)
#> Comparing group IFNbeta time 4 to 0: keeping 20150 of 20150 features (100.0%)
#> Comparing group IFNbeta time 6 to 0: keeping 20150 of 20150 features (100.0%)
#> Comparing group IFNbeta time 8 to 0: keeping 20150 of 20150 features (100.0%)
#> Comparing group IFNbeta time 24 to 0: keeping 20150 of 20150 features (100.0%)
#> Comparing group IFNbeta time 4 to 2: keeping 20150 of 20150 features (100.0%)
#> Comparing group IFNbeta time 6 to 2: keeping 20150 of 20150 features (100.0%)
#> Comparing group IFNbeta time 8 to 2: keeping 20150 of 20150 features (100.0%)
#> Comparing group IFNbeta time 24 to 2: keeping 20150 of 20150 features (100.0%)
#> Comparing group IFNbeta time 6 to 4: keeping 20150 of 20150 features (100.0%)
#> Comparing group IFNbeta time 8 to 4: keeping 20150 of 20150 features (100.0%)
#> Comparing group IFNbeta time 24 to 4: keeping 20150 of 20150 features (100.0%)
#> Comparing group IFNbeta time 8 to 6: keeping 20150 of 20150 features (100.0%)
#> Comparing group IFNbeta time 24 to 6: keeping 20150 of 20150 features (100.0%)
#> Comparing group IFNbeta time 24 to 8: keeping 20150 of 20150 features (100.0%)
#> Comparing group IFNgamma time 2 to 0: keeping 20150 of 20150 features (100.0%)
#> Comparing group IFNgamma time 4 to 0: keeping 20150 of 20150 features (100.0%)
#> Comparing group IFNgamma time 6 to 0: keeping 20150 of 20150 features (100.0%)
#> Comparing group IFNgamma time 8 to 0: keeping 20150 of 20150 features (100.0%)
#> Comparing group IFNgamma time 24 to 0: keeping 20150 of 20150 features (100.0%)
#> Comparing group IFNgamma time 4 to 2: keeping 20150 of 20150 features (100.0%)
#> Comparing group IFNgamma time 6 to 2: keeping 20150 of 20150 features (100.0%)
#> Comparing group IFNgamma time 8 to 2: keeping 20150 of 20150 features (100.0%)
#> Comparing group IFNgamma time 24 to 2: keeping 20150 of 20150 features (100.0%)
#> Comparing group IFNgamma time 6 to 4: keeping 20150 of 20150 features (100.0%)
#> Comparing group IFNgamma time 8 to 4: keeping 20150 of 20150 features (100.0%)
#> Comparing group IFNgamma time 24 to 4: keeping 20150 of 20150 features (100.0%)
#> Comparing group IFNgamma time 8 to 6: keeping 20150 of 20150 features (100.0%)
#> Comparing group IFNgamma time 24 to 6: keeping 20150 of 20150 features (100.0%)
#> Comparing group IFNgamma time 24 to 8: keeping 20150 of 20150 features (100.0%)
#> Comparing group LCMV time 2 to 0: keeping 20150 of 20150 features (100.0%)
#> Comparing group LCMV time 4 to 0: keeping 20150 of 20150 features (100.0%)
#> Comparing group LCMV time 6 to 0: keeping 20150 of 20150 features (100.0%)
#> Comparing group LCMV time 8 to 0: keeping 20150 of 20150 features (100.0%)
#> Comparing group LCMV time 24 to 0: keeping 20150 of 20150 features (100.0%)
#> Comparing group LCMV time 4 to 2: keeping 20150 of 20150 features (100.0%)
#> Comparing group LCMV time 6 to 2: keeping 20150 of 20150 features (100.0%)
#> Comparing group LCMV time 8 to 2: keeping 20150 of 20150 features (100.0%)
#> Comparing group LCMV time 24 to 2: keeping 20150 of 20150 features (100.0%)
#> Comparing group LCMV time 6 to 4: keeping 20150 of 20150 features (100.0%)
#> Comparing group LCMV time 8 to 4: keeping 20150 of 20150 features (100.0%)
#> Comparing group LCMV time 24 to 4: keeping 20150 of 20150 features (100.0%)
#> Comparing group LCMV time 8 to 6: keeping 20150 of 20150 features (100.0%)
#> Comparing group LCMV time 24 to 6: keeping 20150 of 20150 features (100.0%)
#> Comparing group LCMV time 24 to 8: keeping 20150 of 20150 features (100.0%)
#> Comparing group Listeria time 2 to 0: keeping 20150 of 20150 features (100.0%)
#> Comparing group Listeria time 4 to 0: keeping 20150 of 20150 features (100.0%)
#> Comparing group Listeria time 6 to 0: keeping 20150 of 20150 features (100.0%)
#> Comparing group Listeria time 8 to 0: keeping 20150 of 20150 features (100.0%)
#> Comparing group Listeria time 24 to 0: keeping 20150 of 20150 features (100.0%)
#> Comparing group Listeria time 4 to 2: keeping 20150 of 20150 features (100.0%)
#> Comparing group Listeria time 6 to 2: keeping 20150 of 20150 features (100.0%)
#> Comparing group Listeria time 8 to 2: keeping 20150 of 20150 features (100.0%)
#> Comparing group Listeria time 24 to 2: keeping 20150 of 20150 features (100.0%)
#> Comparing group Listeria time 6 to 4: keeping 20150 of 20150 features (100.0%)
#> Comparing group Listeria time 8 to 4: keeping 20150 of 20150 features (100.0%)
#> Comparing group Listeria time 24 to 4: keeping 20150 of 20150 features (100.0%)
#> Comparing group Listeria time 8 to 6: keeping 20150 of 20150 features (100.0%)
#> Comparing group Listeria time 24 to 6: keeping 20150 of 20150 features (100.0%)
#> Comparing group Listeria time 24 to 8: keeping 20150 of 20150 features (100.0%)
#> Comparing group LPS time 2 to 0: keeping 20150 of 20150 features (100.0%)
#> Comparing group LPS time 4 to 0: keeping 20150 of 20150 features (100.0%)
#> Comparing group LPS time 6 to 0: keeping 20150 of 20150 features (100.0%)
#> Comparing group LPS time 8 to 0: keeping 20150 of 20150 features (100.0%)
#> Comparing group LPS time 24 to 0: keeping 20150 of 20150 features (100.0%)
#> Comparing group LPS time 4 to 2: keeping 20150 of 20150 features (100.0%)
#> Comparing group LPS time 6 to 2: keeping 20150 of 20150 features (100.0%)
#> Comparing group LPS time 8 to 2: keeping 20150 of 20150 features (100.0%)
#> Comparing group LPS time 24 to 2: keeping 20150 of 20150 features (100.0%)
#> Comparing group LPS time 6 to 4: keeping 20150 of 20150 features (100.0%)
#> Comparing group LPS time 8 to 4: keeping 20150 of 20150 features (100.0%)
#> Comparing group LPS time 24 to 4: keeping 20150 of 20150 features (100.0%)
#> Comparing group LPS time 8 to 6: keeping 20150 of 20150 features (100.0%)
#> Comparing group LPS time 24 to 6: keeping 20150 of 20150 features (100.0%)
#> Comparing group LPS time 24 to 8: keeping 20150 of 20150 features (100.0%)
#> Comparing group untreated time 8 to 0: keeping 20150 of 20150 features (100.0%)
#> Comparing group untreated time 24 to 0: keeping 20150 of 20150 features (100.0%)
#> Comparing group untreated time 24 to 8: keeping 20150 of 20150 features (100.0%)plot_DE_between_time(data_obj,
de_list = DE_between_time_out$de_list,
fontsize = 8, value = FALSE, nrow = 1, heatmap_width = 3
)
#> Warning: The input is a data frame-like object, convert it to a matrix.
#> Warning: The input is a data frame-like object, convert it to a matrix.
#> Warning: The input is a data frame-like object, convert it to a matrix.
#> Warning: The input is a data frame-like object, convert it to a matrix.
#> Warning: The input is a data frame-like object, convert it to a matrix.
#> Warning: The input is a data frame-like object, convert it to a matrix.
#> Warning: The input is a data frame-like object, convert it to a matrix.
#> Warning: Note: not all columns in the data frame are numeric. The data frame
#> will be converted into a character matrix.
plot_volcano(DE_between_time_out,
group = "IFNgamma", time1 = 0, time2 = 24,
logFC_thres = 3, adjP_thres = 0.05, label = TRUE)
8.2 Between groups at each time point
Default DE criteria: p.adj < 0.05, |log2FC| > 1
All groups use untreated time 0 data, so no DE features at time point 0
DE_between_group_out <- DE_between_group(data_obj, assay = 2,
group = "untreated")
#> Non-NA replicate number filter not specified. Using minimum number of replicates across all groups and time points: 0
#> Comparing group Candida to untreated at Time 0: keeping 20150 of 20150 features (100.0%)
#> Warning: Zero sample variances detected, have been offset away from zero
#> Comparing group Candida to untreated at Time 8: keeping 20150 of 20150 features (100.0%)
#> Warning: Zero sample variances detected, have been offset away from zero
#> Comparing group Candida to untreated at Time 24: keeping 20150 of 20150 features (100.0%)
#> Warning: Zero sample variances detected, have been offset away from zero
#> Comparing group IFNbeta to untreated at Time 0: keeping 20150 of 20150 features (100.0%)
#> Warning: Zero sample variances detected, have been offset away from zero
#> Comparing group IFNbeta to untreated at Time 8: keeping 20150 of 20150 features (100.0%)
#> Warning: Zero sample variances detected, have been offset away from zero
#> Comparing group IFNbeta to untreated at Time 24: keeping 20150 of 20150 features (100.0%)
#> Warning: Zero sample variances detected, have been offset away from zero
#> Comparing group IFNgamma to untreated at Time 0: keeping 20150 of 20150 features (100.0%)
#> Warning: Zero sample variances detected, have been offset away from zero
#> Comparing group IFNgamma to untreated at Time 8: keeping 20150 of 20150 features (100.0%)
#> Warning: Zero sample variances detected, have been offset away from zero
#> Comparing group IFNgamma to untreated at Time 24: keeping 20150 of 20150 features (100.0%)
#> Warning: Zero sample variances detected, have been offset away from zero
#> Comparing group LCMV to untreated at Time 0: keeping 20150 of 20150 features (100.0%)
#> Warning: Zero sample variances detected, have been offset away from zero
#> Comparing group LCMV to untreated at Time 8: keeping 20150 of 20150 features (100.0%)
#> Warning: Zero sample variances detected, have been offset away from zero
#> Comparing group LCMV to untreated at Time 24: keeping 20150 of 20150 features (100.0%)
#> Warning: Zero sample variances detected, have been offset away from zero
#> Comparing group Listeria to untreated at Time 0: keeping 20150 of 20150 features (100.0%)
#> Warning: Zero sample variances detected, have been offset away from zero
#> Comparing group Listeria to untreated at Time 8: keeping 20150 of 20150 features (100.0%)
#> Warning: Zero sample variances detected, have been offset away from zero
#> Comparing group Listeria to untreated at Time 24: keeping 20150 of 20150 features (100.0%)
#> Warning: Zero sample variances detected, have been offset away from zero
#> Comparing group LPS to untreated at Time 0: keeping 20150 of 20150 features (100.0%)
#> Warning: Zero sample variances detected, have been offset away from zero
#> Comparing group LPS to untreated at Time 8: keeping 20150 of 20150 features (100.0%)
#> Warning: Zero sample variances detected, have been offset away from zero
#> Comparing group LPS to untreated at Time 24: keeping 20150 of 20150 features (100.0%)
#> Warning: Zero sample variances detected, have been offset away from zero
plot_volcano(DE_between_group_out,
group1 = "untreated", group2 = "IFNgamma", time = 24,
logFC_thres = 3, adjP_thres = 0.05, label = TRUE)
9 Differentially expressed genes (DEG) classified by feature properties
9.1 Calculate expression ratio, randomness, maximum overall fold change
use data not normalised to time 0, replicates merged by mean
data_obj_merged_list <- calc_feature_property(data_obj_merged_list, threshold = 0)
property_random_fc <- summarise_feature_property(data_obj_merged_list)
head(property_random_fc)9.2 Non-random DEG in each group
DEG: max overall fold change >= fold change threshold (2 fold)
filter_list <- list() # genes with at least certain values > threshold (0)
non_random_list <- list()
non_random_up_list <- list()
non_random_down_list <- list()
filter_ratio <- 0.5 # minimum expressed in 50% samples
max_fc_thres <- 1 # two fold changes in log2 scale
for (i in unique(data_obj$Group)) {
filter_list[[i]] <- property_random_fc %>%
filter(Group == i) %>%
filter(Exp_ratio > filter_ratio) %>%
pull(Feature)
non_random_list[[i]] <- property_random_fc %>%
filter(Group == i) %>%
filter(P_trend < 0.05 & Exp_ratio > filter_ratio) %>%
pull(Feature)
non_random_up_list[[i]] <- property_random_fc %>%
filter(Group == i) %>%
filter(P_trend < 0.05 & Exp_ratio > filter_ratio &
Max_FC >= max_fc_thres & Max_FC_time > 0) %>%
pull(Feature)
non_random_down_list[[i]] <- property_random_fc %>%
filter(Group == i) %>%
filter(P_trend < 0.05 & Exp_ratio > filter_ratio &
Max_FC >= max_fc_thres & Max_FC_time < 0) %>%
pull(Feature)
}#> $Candida
#> [1] 168
#>
#> $IFNbeta
#> [1] 579
#>
#> $IFNgamma
#> [1] 336
#>
#> $LCMV
#> [1] 591
#>
#> $Listeria
#> [1] 372
#>
#> $LPS
#> [1] 279
#>
#> $untreated
#> [1] 0
non_random_down_list %>% lapply(length)#> $Candida
#> [1] 146
#>
#> $IFNbeta
#> [1] 2036
#>
#> $IFNgamma
#> [1] 808
#>
#> $LCMV
#> [1] 621
#>
#> $Listeria
#> [1] 3353
#>
#> $LPS
#> [1] 1704
#>
#> $untreated
#> [1] 0The untreated group with only 3 time points doesn’t have non-random DEG.
9.3 Visualise example non-random DEGs
for (i in unique(data_obj$Group)) {
random <- sample(non_random_down_list[[i]],
min(30, length(non_random_down_list[[i]])))
if (length(random) == 0) {
next
}
# mean, sd
plot_trend(table_mean_sd_orig,
groups = i,
features = random, title = "Overall down DEG"
)
}
9.4 GO enrichment of overall DEGs
go_list_up_merged <- enrichGO_list(
gene_list = non_random_up_list,
OrgDb = org.Mm.eg.db,
universe_list = filter_list,
keyType = "SYMBOL",
simplify = TRUE
)
#> GO category not specified. Using all three: BP, MF, CC.
#> Performing GO enrichment for category: BP
#> Processing gene list: Candida
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: IFNbeta
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: IFNgamma
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: LCMV
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: Listeria
#> 'select()' returned 1:many mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: LPS
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Gene list untreated is empty. Skipping.
#> Performing GO enrichment for category: MF
#> Processing gene list: Candida
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: IFNbeta
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: IFNgamma
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: LCMV
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: Listeria
#> 'select()' returned 1:many mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: LPS
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Gene list untreated is empty. Skipping.
#> Performing GO enrichment for category: CC
#> Processing gene list: Candida
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: IFNbeta
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: IFNgamma
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: LCMV
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: Listeria
#> 'select()' returned 1:many mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: LPS
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Gene list untreated is empty. Skipping.
#> Merging GO enrichment results across gene lists for each category.
go_list_down_merged <- enrichGO_list(
gene_list = non_random_down_list,
OrgDb = org.Mm.eg.db,
universe_list = filter_list,
keyType = "SYMBOL",
simplify = TRUE
)
#> GO category not specified. Using all three: BP, MF, CC.
#> Performing GO enrichment for category: BP
#> Processing gene list: Candida
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: IFNbeta
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: IFNgamma
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: LCMV
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: Listeria
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: LPS
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Gene list untreated is empty. Skipping.
#> Performing GO enrichment for category: MF
#> Processing gene list: Candida
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: IFNbeta
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: IFNgamma
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: LCMV
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: Listeria
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: LPS
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Gene list untreated is empty. Skipping.
#> Performing GO enrichment for category: CC
#> Processing gene list: Candida
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: IFNbeta
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: IFNgamma
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: LCMV
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: Listeria
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: LPS
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Gene list untreated is empty. Skipping.
#> Merging GO enrichment results across gene lists for each category.plot_GO(go_list_up_merged$all,
plot_dotplot = TRUE,
showCategory_dotplot = 5,
label = "overall up DEGs"
)
plot_GO(go_list_up_merged$simplified,
plot_dotplot = TRUE,
showCategory_dotplot = 5,
label = "overall up DEGs (simplified)"
)
plot_GO(go_list_down_merged$all,
plot_dotplot = TRUE,
showCategory_dotplot = 5,
label = "overall down DEGs"
)
plot_GO(go_list_down_merged$simplified,
plot_dotplot = TRUE,
showCategory_dotplot = 5,
label = "overall down DEGs (simplified)"
)
plot_GO(go_list_up_merged$all,
plot_emapplot = TRUE,
showCategory_emapplot = 5,
label = "overall up DEGs"
)
#> Coordinate system already present.
#> ℹ Adding new coordinate system, which will replace the existing one.
#> Coordinate system already present.
#> ℹ Adding new coordinate system, which will replace the existing one.
#> Coordinate system already present.
#> ℹ Adding new coordinate system, which will replace the existing one.
plot_GO(go_list_down_merged$all,
plot_emapplot = TRUE,
showCategory_emapplot = 5,
label = "overall down DEGs"
)
#> Coordinate system already present.
#> ℹ Adding new coordinate system, which will replace the existing one.
#> Coordinate system already present.
#> ℹ Adding new coordinate system, which will replace the existing one.
#> Coordinate system already present.
#> ℹ Adding new coordinate system, which will replace the existing one.
9.5 Group-specifically expressed genes
genes expressed in only one group with >50% values >0
filter_ratio <- 0.5
group_unique_genes_list = list()
for (i in unique(property_random_fc$Group)) {
group_unique_genes_list[[i]] = group_specific_features(property_random_fc,
groups = i,
filter_ratio = filter_ratio,
group_pct = 1,
GO = FALSE, genename = FALSE,
keytype = "SYMBOL",
org.db = "org.Hs.eg.db")
}
#> Filtering criteria: >=50% values >0 in >=1 of groups: Candida
#> Filtering criteria: >=50% values >0 in >=1 of groups: IFNbeta
#> Filtering criteria: >=50% values >0 in >=1 of groups: IFNgamma
#> Filtering criteria: >=50% values >0 in >=1 of groups: LCMV
#> Filtering criteria: >=50% values >0 in >=1 of groups: Listeria
#> Filtering criteria: >=50% values >0 in >=1 of groups: LPS
#> Filtering criteria: >=50% values >0 in >=1 of groups: untreatedplot_trend(table_mean_sd_orig,
features = group_unique_genes_list[["untreated"]],
title = paste0("Genes with >", 100 * filter_ratio, "% expression in only untreated group"))
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Warning: Removed 213 rows containing missing values or values outside the scale range
#> (`geom_point()`).
Features with enough values in only one group will be excluded from variance decomposition, because the group contributed variance is unknown.
10 Segmented regression analysis
10.1 Impute with group minimum (replicates merged by mean)
no NA values in this dataset so imputation will not change the data
data_obj_merged_imp_list <- impute_groups(data_obj_merged_list)
data_obj_merged_imp_list %>% lapply(function(x) ncol(x))#> $Candida
#> [1] 6
#>
#> $IFNbeta
#> [1] 6
#>
#> $IFNgamma
#> [1] 6
#>
#> $LCMV
#> [1] 6
#>
#> $Listeria
#> [1] 6
#>
#> $LPS
#> [1] 6
#>
#> $untreated
#> [1] 310.2 Run Trendy
time consuming, load from presaved results
assays(data_obj_merged_imp_list$IFNbeta)[[1]] %>%
rowMeans() %>%
summary()
res_list <- run_Trendy(data_obj_merged_imp_list,
maxK = 1,
minNumInSeg = 2,
meanCut = 0,
NCores = 16
)# save(res_list, file = paste0("output\\TiDEomics_application_1_trendy_results_k1_", Sys.Date(), ".rda"))
load("output\\TiDEomics_application_1_trendy_results_k1_2026-04-24.rda")When # segments = 2 (maxK = 1) and minNumInSeg = 2, at least 4 time points are needed.
Untreated group with only 3 time points were not analysed by run_Trendy()
10.3 Visualise results
Segments of example genes
plot_segments(data_obj_merged_imp_list,
res_list,
feature = c("Isg15", "Tnf"), # example features
ylab = "Log2 (raw count + 1)"
)
#> Plotting segmented regression for group: Candida
#> Plotting segmented regression for group: IFNbeta
#> Plotting segmented regression for group: IFNgamma
#> Plotting segmented regression for group: LCMV
#> Plotting segmented regression for group: Listeria
#> Plotting segmented regression for group: LPS
Breakpoint distribution and details per group
plot_breakpoints(res_list)
#> Warning in (function (..., deparse.level = 1) : number of columns of result is
#> not a multiple of vector length (arg 11)
#> Warning in (function (..., deparse.level = 1) : number of columns of result is
#> not a multiple of vector length (arg 12)
#> Warning in (function (..., deparse.level = 1) : number of columns of result is
#> not a multiple of vector length (arg 116)
#> Warning in (function (..., deparse.level = 1) : number of columns of result is
#> not a multiple of vector length (arg 34)
#> Warning in (function (..., deparse.level = 1) : number of columns of result is
#> not a multiple of vector length (arg 90)
#> Warning in (function (..., deparse.level = 1) : number of columns of result is
#> not a multiple of vector length (arg 40)
trendy_summary <- summarise_Trendy(res_list)
#> Warning in (function (..., deparse.level = 1) : number of columns of result is
#> not a multiple of vector length (arg 11)
#> Warning in (function (..., deparse.level = 1) : number of columns of result is
#> not a multiple of vector length (arg 12)
#> Warning in (function (..., deparse.level = 1) : number of columns of result is
#> not a multiple of vector length (arg 116)
#> Warning in (function (..., deparse.level = 1) : number of columns of result is
#> not a multiple of vector length (arg 34)
#> Warning in (function (..., deparse.level = 1) : number of columns of result is
#> not a multiple of vector length (arg 90)
#> Warning in (function (..., deparse.level = 1) : number of columns of result is
#> not a multiple of vector length (arg 40)
trendy_summary %>% head()Example genes with different segment patterns
# plot example genes with different trend patterns
i <- "LPS"
for (j in names(trendy_list[[i]])) {
if (!is.null(trendy_list[[i]][[j]])) {
plot_trend(table_mean_sd_orig,
groups = i,
features = trendy_list[[i]][[j]] %>%
sample(min(15, length(trendy_list[[i]][[j]]))),
title = paste0(j)
)
}
}
11 Variance decomposition analysis
11.1 Filter input
Input features are filtered to have at least 50% values > 0 in at least 2 groups, as comparison between groups is needed.
filter_ratio = 0.5
decomp_filter_genes <- group_specific_features(property_random_fc,
filter_ratio = filter_ratio,
group_pct = 2 / length(unique(data_obj$Group)),
GO = FALSE, genename = FALSE
)
#> Filtering criteria: >=50% values >0 in >=2 of groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated11.2 Original data (not time 0 normalised) with replicates
Compared to normalised data, only using original data can identify genes with stable but different levels of expression in groups.
var_decomp_1 <- decomp_variance(data_obj,
features = decomp_filter_genes,
assay = 1, core = 16
)
plot_variance(var_decomp_1, rank = "Time", top_n = 20)
#> Features not specified. Plotting top 20 features ranked by Time.
plot_variance(var_decomp_1, rank = "Group", top_n = 20)
#> Features not specified. Plotting top 20 features ranked by Group.
plot_trend(table_mean_sd_orig,
features = var_decomp_1 %>% arrange(desc(Time)) %>% head(30) %>% pull(Feature),
title = "Top 30 time dependent genes"
)
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Warning: Removed 90 rows containing missing values or values outside the scale range
#> (`geom_point()`).
plot_trend(table_mean_sd_orig,
features = var_decomp_1 %>% arrange(desc(Group)) %>% head(30) %>% pull(Feature),
title = "Top 30 group dependent genes"
)
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Warning: Removed 90 rows containing missing values or values outside the scale range
#> (`geom_point()`).
11.2.1 Residual variance quantified data quality
plot_trend(table_mean_sd_orig,
features = var_decomp_1 %>% arrange(desc(Residual)) %>% head(30) %>% pull(Feature),
title = "Top 30 genes with high residual variance"
)
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Warning: Removed 90 rows containing missing values or values outside the scale range
#> (`geom_point()`).
plot_trend(table_mean_sd_orig,
features = var_decomp_1 %>% arrange(Residual) %>% head(30) %>% pull(Feature),
title = "Top 30 genes with low residual variance"
)
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Warning: Removed 90 rows containing missing values or values outside the scale range
#> (`geom_point()`).
11.2.2 GO enrichment of genes ranked by group/time-contributed variance
gse_group <- enrichGO_rank(var_decomp_1,
gene_rank_by = "Group",
OrgDb = org.Mm.eg.db,
keyType = "SYMBOL", category = "BP",
go_rank_by = "p.adjust")
#> Warning in prepare_gsea_inputs(geneList, scoreType, exponent): There are ties
#> in the preranked stats (0.52% of the list). The order of those tied genes will
#> be arbitrary, which may produce unexpected results.
#> Warning in gsea(geneList = geneList, gene_sets = geneSets, minGSSize =
#> minGSSize, : There were 18 pathways for which P-values were not calculated
#> properly due to unbalanced gene-level statistic values. For such pathways
#> pvalue, NES and log2err are set to NA. You can try to increase nPermSimple.
#> Warning in gsea(geneList = geneList, gene_sets = geneSets, minGSSize =
#> minGSSize, : For some pathways, in reality P-values are less than 1e-10. You
#> can set the eps argument to zero for better estimation.
#> Warning in calculate_qvalue(gsea_res$pvalue): Invalid p-values detected (NA,
#> non-finite, <0, or >1). qvalue will be computed on valid p-values only.
#> Warning in enrichit::gsea_gson(geneList = geneList, exponent = exponent, : NA
#> values detected in gene set IDs. Replacing with string 'NA'.
#> Warning in enrichit::gsea_gson(geneList = geneList, exponent = exponent, :
#> Duplicate gene set IDs detected: NA... (Total 1). Unique suffixes added.
#> Removing NA ID gene sets.
enrichplot::gseaplot2(gse_group, geneSetID = 1:3, base_size = 8) 
gse_group <- enrichGO_rank(var_decomp_1,
gene_rank_by = "Time",
OrgDb = org.Mm.eg.db,
keyType = "SYMBOL", category = "BP",
go_rank_by = "p.adjust")
#> Warning in prepare_gsea_inputs(geneList, scoreType, exponent): There are ties
#> in the preranked stats (0.53% of the list). The order of those tied genes will
#> be arbitrary, which may produce unexpected results.
#> Warning in gsea(geneList = geneList, gene_sets = geneSets, minGSSize =
#> minGSSize, : There were 17 pathways for which P-values were not calculated
#> properly due to unbalanced gene-level statistic values. For such pathways
#> pvalue, NES and log2err are set to NA. You can try to increase nPermSimple.
#> Warning in gsea(geneList = geneList, gene_sets = geneSets, minGSSize =
#> minGSSize, : For some pathways, in reality P-values are less than 1e-10. You
#> can set the eps argument to zero for better estimation.
#> Warning in calculate_qvalue(gsea_res$pvalue): Invalid p-values detected (NA,
#> non-finite, <0, or >1). qvalue will be computed on valid p-values only.
#> Warning in enrichit::gsea_gson(geneList = geneList, exponent = exponent, : NA
#> values detected in gene set IDs. Replacing with string 'NA'.
#> Warning in enrichit::gsea_gson(geneList = geneList, exponent = exponent, :
#> Duplicate gene set IDs detected: NA... (Total 1). Unique suffixes added.
#> Removing NA ID gene sets.
enrichplot::gseaplot2(gse_group, geneSetID = 1:3, base_size = 8)
11.3 Time 0 normalised data with replicates
Time 0 normalised data will be used as WGCNA input, so variance decomposition is also performed on time 0 normalised data to filter for WGCNA input
var_decomp_2 <- decomp_variance(data_obj,
features = decomp_filter_genes,
assay = 2, core = 16
)
plot_variance(var_decomp_2, rank = "Time", top_n = 20)
#> Features not specified. Plotting top 20 features ranked by Time.
plot_variance(var_decomp_2, rank = "Group", top_n = 20)
#> Features not specified. Plotting top 20 features ranked by Group.
11.4 Compare variance with / without time 0 normalisation
var_decomp_1: not normalised to time 0
var_decomp_2: normalised to time 0
In this dataset, all groups have the same starting status (untreated), so after normalisation to time 0, the variance decomposition is the same as the original data.
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 17.71 44.19 57.96 60.51 76.15 100.00
var_res_q3_1 <- quantile(var_decomp_1$Residual, 0.75)
var_decomp_2$Residual %>% summary()#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 17.71 44.19 57.96 60.51 76.15 100.00
var_res_q3_2 <- quantile(var_decomp_2$Residual, 0.75)var_decomp_merged <- merge(var_decomp_1, var_decomp_2,
by = "Feature",
suffixes = c(".orig", ".norm")
)
summary(var_decomp_merged$Residual.orig / var_decomp_merged$Residual.norm)#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 1 1 1 1 1 112 Module identificiation with WGCNA
WGCNA reference: https://github.com/edo98811/WGCNA_official_documentation/
12.1 Filter input
exclude genes with high residual variance
# time 0 normalised
gene_resid_2 <- var_decomp_2 %>%
filter(Residual < var_res_q3_2) %>%
pull(Feature)
length(gene_resid_2) # only including genes with >=50% values >0 in >=2 groups#> [1] 10909exclude genes that are random in all groups
non_random_union <- Reduce(union, non_random_list)
length(non_random_union) # all non-random genes have >=50% non NA#> [1] 7628include genes with at least 50% values > 0 in only one group, which were excluded from variance decomposition analysis
group_unique_genes = Reduce(union, group_unique_genes_list)
stopifnot(group_unique_genes %>% unique() %>% length() == length(group_unique_genes))
length(group_unique_genes) #> [1] 769use time 0 normalised data with replicates
filter_wgcna <- intersect(gene_resid_2, non_random_union) %>%
union(group_unique_genes)
length(filter_wgcna)#> [1] 7292
data_obj_filtered <- data_obj[filter_wgcna, ]
dim(data_obj_filtered)#> [1] 7292 7612.2 Choose power
wgcna_input <- prepare_WGCNA(data_obj_filtered,
assay = 2,
blockSize = 20000,
powers = c(seq(1, 30, by = 1)),
networkType = "signed",
RsquaredCut = 0.8,
MeanConnectivity = 100
)#> Allowing multi-threading with up to 24 threads.
#> pickSoftThreshold: calculating connectivity for given powers...
#> ..working on genes 1 through 7292 of 7292
#> Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
#> 1 1 0.8560 3.500 0.8220 4390.0 4690.00 5240
#> 2 2 0.5040 1.020 0.4130 2910.0 3210.00 4100
#> 3 3 0.0547 0.214 -0.1110 2060.0 2290.00 3360
#> 4 4 0.0303 -0.131 -0.1970 1530.0 1680.00 2830
#> 5 5 0.1870 -0.312 -0.0443 1170.0 1250.00 2420
#> 6 6 0.3340 -0.425 0.1540 918.0 950.00 2100
#> 7 7 0.4470 -0.536 0.3360 735.0 730.00 1850
#> 8 8 0.5130 -0.618 0.4470 598.0 567.00 1630
#> 9 9 0.5610 -0.707 0.5550 493.0 444.00 1450
#> 10 10 0.5880 -0.783 0.6250 410.0 351.00 1300
#> 11 11 0.6240 -0.835 0.6900 345.0 279.00 1170
#> 12 12 0.6360 -0.901 0.7300 292.0 223.00 1060
#> 13 13 0.6510 -0.955 0.7640 249.0 180.00 956
#> 14 14 0.6700 -0.998 0.7970 213.0 145.00 869
#> 15 15 0.6860 -1.040 0.8230 184.0 118.00 792
#> 16 16 0.6930 -1.080 0.8420 159.0 95.90 724
#> 17 17 0.6930 -1.130 0.8500 139.0 79.00 663
#> 18 18 0.6970 -1.170 0.8600 121.0 65.30 609
#> 19 19 0.6950 -1.220 0.8630 107.0 54.70 560
#> 20 20 0.7110 -1.240 0.8790 93.8 45.90 517
#> 21 21 0.7100 -1.270 0.8820 82.9 38.40 477
#> 22 22 0.7160 -1.310 0.8900 73.5 32.10 442
#> 23 23 0.7280 -1.320 0.9010 65.4 27.20 410
#> 24 24 0.7200 -1.360 0.8980 58.3 23.20 380
#> 25 25 0.7300 -1.380 0.9070 52.2 19.90 354
#> 26 26 0.7300 -1.410 0.9090 46.8 17.10 329
#> 27 27 0.7410 -1.420 0.9160 42.0 14.70 307
#> 28 28 0.7380 -1.450 0.9150 37.9 12.70 287
#> 29 29 0.7470 -1.460 0.9210 34.2 11.10 268
#> 30 30 0.7550 -1.470 0.9270 30.9 9.68 251
The R2 didn’t reach 0.8 under power 30, so choose the lowest power that reached the plateau.
12.3 Run WGCNA
picked_power <- 12
minCoreKME <- 0.9
minKMEtoStay <- 0.7
reassignThreshold <- 0
net <- run_WGCNA(wgcna_input,
power = picked_power,
maxBlockSize = 20000,
numericLabels = TRUE,
minCoreKME = minCoreKME,
minKMEtoStay = minKMEtoStay,
reassignThreshold = reassignThreshold,
verbose = 3
)#> Allowing multi-threading with up to 24 threads.
#> Calculating module eigengenes block-wise from all genes
#> Flagging genes and samples with too many missing values...
#> ..step 1
#> ..Working on block 1 .
#> TOM calculation: adjacency..
#> ..will not use multithreading.
#> Fraction of slow calculations: 0.000000
#> ..connectivity..
#> ..matrix multiplication (system BLAS)..
#> ..normalization..
#> ..done.
#> ....clustering..
#> ....detecting modules..
#> ....calculating module eigengenes..
#> ....checking kME in modules..
#> ..removing 979 genes from module 1 because their KME is too low.
#> ..removing 295 genes from module 2 because their KME is too low.
#> ..removing 373 genes from module 3 because their KME is too low.
#> ..removing 121 genes from module 4 because their KME is too low.
#> ..removing 72 genes from module 5 because their KME is too low.
#> ..removing 178 genes from module 6 because their KME is too low.
#> ..removing 102 genes from module 7 because their KME is too low.
#> ..removing 132 genes from module 8 because their KME is too low.
#> ..removing 68 genes from module 9 because their KME is too low.
#> ..removing 48 genes from module 10 because their KME is too low.
#> ..removing 30 genes from module 12 because their KME is too low.
#> ..removing 14 genes from module 13 because their KME is too low.
#> ..removing 27 genes from module 14 because their KME is too low.
#> ..removing 13 genes from module 16 because their KME is too low.
#> ..removing 1 genes from module 18 because their KME is too low.
#> ..merging modules that are too close..
#> mergeCloseModules: Merging modules whose distance is less than 0.15
#> Calculating new MEs...
# save(net, file = paste0("output\\GSE263759_WGCNA_net_p", picked_power, "_", Sys.Date(), ".rda"))
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
gene_module <- data.frame(Module = as.factor(net$colors)) %>%
tibble::rownames_to_column("Feature") %>%
arrange(Module)
gene_module %>%
dplyr::group_by(Module) %>%
dplyr::summarise(n = n())12.4 Visualise example genes in each module
for (i in seq(1, n_module)) {
module_genes <- gene_module %>%
filter(Module == i) %>%
dplyr::select(Feature) %>%
pull()
random <- sample(module_genes, min(30, length(module_genes)))
plot_trend(table_mean_sd_norm,
errorbar = FALSE,
features = random,
title = paste0("Module ", i, " random genes")
)
}
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Warning: Removed 90 rows containing missing values or values outside the scale range
#> (`geom_point()`).
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Warning: Removed 90 rows containing missing values or values outside the scale range
#> (`geom_point()`).
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Warning: Removed 90 rows containing missing values or values outside the scale range
#> (`geom_point()`).
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Warning: Removed 90 rows containing missing values or values outside the scale range
#> (`geom_point()`).
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Warning: Removed 90 rows containing missing values or values outside the scale range
#> (`geom_point()`).
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Warning: Removed 90 rows containing missing values or values outside the scale range
#> (`geom_point()`).
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Warning: Removed 90 rows containing missing values or values outside the scale range
#> (`geom_point()`).
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Warning: Removed 90 rows containing missing values or values outside the scale range
#> (`geom_point()`).
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Warning: Removed 90 rows containing missing values or values outside the scale range
#> (`geom_point()`).
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Warning: Removed 90 rows containing missing values or values outside the scale range
#> (`geom_point()`).
12.5 Summarise module segment pattern
integrate with output from Trendy
trendy_genes_pattern_list <- summarise_module_pattern(
module = gene_module,
trendy_summary = trendy_summary
)
#> Most common pattern in each module (Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS):
#> Module 1: NA, NA, NA, NA, stable_stable, NA
#> Module 2: NA, down, down, NA, down, down
#> Module 3: NA, down_stable, NA, NA, stable_stable, NA
#> Module 4: NA, up_down, NA, stable_up, NA, NA
#> Module 5: NA, NA, NA, up, NA, NA
#> Module 6: NA, NA, NA, NA, up, NA
#> Module 7: down, down_stable, NA, NA, up_down, up_down
#> Module 8: NA, NA, NA, NA, stable_stable, NA
#> Module 9: NA, up_down, NA, NA, NA, NA
#> Module 10: NA, up_stable, NA, NA, NA, NA
#> Module 1 top 5 patterns:#> Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS Count
#> 1 NA, NA, NA, NA, stable_stable, NA 28
#> 2 NA, down_up, NA, NA, down_stable, NA 22
#> 3 NA, down_stable, NA, NA, stable_stable, NA 20
#> 4 NA, down_up, NA, down, down_stable, stable_stable 18
#> 5 NA, down_stable, NA, down, down_stable, down_stable 14
#> Module 2 top 5 patterns:#> Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS Count
#> 1 NA, down, down, NA, down, down 13
#> 2 NA, down, NA, NA, down, down 8
#> 3 NA, down, down, down, down, down 7
#> 4 NA, down, NA, down, down, down 5
#> 5 NA, down, stable_stable, down, down, down 4
#> Module 3 top 5 patterns:#> Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS Count
#> 1 NA, down_stable, NA, NA, stable_stable, NA 11
#> 2 NA, down_up, NA, NA, down_stable, NA 11
#> 3 NA, down_up, NA, NA, stable_stable, NA 11
#> 4 NA, NA, NA, NA, stable_stable, NA 7
#> 5 NA, down, NA, NA, down, NA 5
#> Module 4 top 5 patterns:#> Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS Count
#> 1 NA, up_down, NA, stable_up, NA, NA 5
#> 2 NA, up_stable, NA, NA, NA, NA 5
#> 3 NA, up_stable, NA, up, NA, NA 4
#> 4 NA, up_stable, up_stable, NA, NA, NA 4
#> 5 NA, up_down, NA, NA, NA, NA 3
#> Module 5 top 5 patterns:#> Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS Count
#> 1 NA, NA, NA, up, NA, NA 23
#> 2 NA, up, NA, NA, NA, NA 7
#> 3 up_stable, NA, NA, up, NA, NA 7
#> 4 up, NA, NA, NA, NA, NA 6
#> 5 NA, NA, NA, NA, stable_stable, NA 5
#> Module 6 top 5 patterns:#> Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS Count
#> 1 NA, NA, NA, NA, up, NA 5
#> 2 NA, NA, NA, NA, up_stable, up_down 4
#> 3 NA, NA, NA, stable_up, up_stable, up_stable 3
#> 4 NA, NA, NA, NA, NA, stable_stable 2
#> 5 NA, NA, NA, NA, up, up 2
#> Module 7 top 5 patterns:#> Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS Count
#> 1 down, down_stable, NA, NA, up_down, up_down 2
#> 2 NA, down, NA, down_up, up_down, up_down 2
#> 3 NA, NA, NA, up, NA, NA 2
#> 4 down, NA, down_stable, stable_stable, NA, up_stable 1
#> 5 NA, down, down, stable_up, NA, up_stable 1
#> Module 8 top 5 patterns:#> Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS Count
#> 1 NA, NA, NA, NA, stable_stable, NA 5
#> 2 NA, down_up, NA, NA, stable_stable, stable_stable 3
#> 3 NA, down_stable, NA, NA, stable_stable, NA 2
#> 4 up, NA, NA, NA, NA, NA 2
#> 5 up_stable, NA, NA, NA, NA, NA 2
#> Module 9 top 5 patterns:#> Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS Count
#> 1 NA, up_down, NA, NA, NA, NA 2
#> 2 NA, up_down, up, NA, NA, NA 2
#> 3 NA, up_stable, NA, NA, NA, NA 2
#> 4 NA, up_stable, up_stable, NA, NA, NA 2
#> 5 NA, up_up, up_stable, NA, NA, NA 2
#> Module 10 top 5 patterns:#> Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS Count
#> 1 NA, up_stable, NA, NA, NA, NA 4
#> 2 NA, up_down, NA, NA, NA, NA 3
#> 3 up, up_down, NA, up, NA, NA 3
#> 4 NA, NA, NA, up, NA, NA 2
#> 5 up_stable, NA, NA, up, NA, NA 212.6 GO enrichment of modules
background: all genes included in WGCNA
module_wgcna <- list()
background_wgcna <- gene_module$Feature
for (i in seq(0, n_module)) {
module_wgcna[[as.character(i)]] <- gene_module %>%
filter(Module == i) %>%
pull(Feature)
}go_list_wgcna_merged <- enrichGO_list(
gene_list = module_wgcna,
OrgDb = org.Mm.eg.db,
universe = background_wgcna,
keyType = "SYMBOL",
simplify = TRUE,
minGSSize = 10,
maxGSSize = 300
)
#> GO category not specified. Using all three: BP, MF, CC.
#> Performing GO enrichment for category: BP
#> Processing gene list: 0
#> 'select()' returned 1:many mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: 1
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: 2
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: 3
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: 4
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: 5
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: 6
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: 7
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: 8
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: 9
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: 10
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Performing GO enrichment for category: MF
#> Processing gene list: 0
#> 'select()' returned 1:many mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: 1
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: 2
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: 3
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: 4
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: 5
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: 6
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: 7
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: 8
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: 9
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: 10
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Performing GO enrichment for category: CC
#> Processing gene list: 0
#> 'select()' returned 1:many mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: 1
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: 2
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: 3
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: 4
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: 5
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: 6
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: 7
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: 8
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: 9
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Processing gene list: 10
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Simplifying GO terms.
#> Warning: semData is not provided. It will be calculated automatically.
#> Merging GO enrichment results across gene lists for each category.go_list_wgcna_merged_no0 <- lapply(go_list_wgcna_merged, function(y) {
lapply(y, function (x) {
x@compareClusterResult <- x@compareClusterResult %>%
filter(Cluster != "0")
return(x)
})
})plot_GO(go_list_wgcna_merged_no0$all,
plot_dotplot = TRUE,
showCategory_dotplot = 5,
label = "WGCNA modules"
)
plot_GO(go_list_wgcna_merged_no0$simplified,
plot_dotplot = TRUE,
showCategory_dotplot = 5,
label = "WGCNA modules (simplified)"
)
plot_GO(go_list_wgcna_merged_no0$all,
plot_cnetplot = TRUE,
showCategory_cnetplot = 5,
label = "WGCNA modules"
)
#> The following groups are missing in the custom palette: 1, 2, 3, 4, 5, 6, 7, 9, 10
#> The following groups are missing in the custom palette: 1, 2, 3, 6
#> The following groups are missing in the custom palette: 1, 2, 3, 4, 5, 6
plot_GO(go_list_wgcna_merged_no0$all,
plot_emapplot = TRUE,
showCategory_emapplot = 5,
label = "WGCNA modules"
)
#> Coordinate system already present.
#> ℹ Adding new coordinate system, which will replace the existing one.
#> The following groups are missing in the custom palette: 1, 2, 3, 4, 5, 6, 7, 9, 10
#> Coordinate system already present.
#> ℹ Adding new coordinate system, which will replace the existing one.
#> The following groups are missing in the custom palette: 1, 2, 3, 6
#> Coordinate system already present.
#> ℹ Adding new coordinate system, which will replace the existing one.
#> The following groups are missing in the custom palette: 1, 2, 3, 4, 5, 6
12.7 Visualise module profiles (heatmap, mean profile, GO enrichment)
12.7.1 Vertical alignment
unscaled (same as the input)
modules_v_unscaled <- plot_modules_v(gene_module %>% filter(Module != 0),
data_obj_merged, scale = FALSE,
ylabel = "Log2 (count+1)\nnormalised to Time 0",
height_ratio = 3
)
modules_v_unscaled
#> Warning: Removed 13752 rows containing non-finite outside the scale range
#> (`stat_summary()`).
scaled
modules_v_scaled <- plot_modules_v(gene_module %>% filter(Module != 0),
data_obj_merged, scale = TRUE,
ylabel = "Z-score of log2 (count+1)\nnormalised to Time 0",
height_ratio = 3
)
modules_v_scaled
#> Warning: Removed 13752 rows containing non-finite outside the scale range
#> (`stat_summary()`).
12.7.2 Horizontal alignment
scaled
plot_modules_h(gene_module %>% filter(Module != 0),
data_obj_merged, scale = TRUE,
ylabel = "Z-score",
go_list = go_list_wgcna_merged$all,
go_category = "BP",
fontsize = 7,
heatmap_width = 8,
heatmap_height = 12,
width = 15,
height = 10
)
#> Warning: Removed 8568 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 1326 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 1248 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 822 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 582 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 423 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 240 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 225 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 177 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 141 rows containing non-finite outside the scale range
#> (`stat_summary()`).
plot_modules_h(gene_module %>% filter(Module != 0),
data_obj_merged, scale = TRUE,
ylabel = "Z-score",
go_list = go_list_wgcna_merged$simplified,
go_category = "BP",
fontsize = 7,
heatmap_width = 8,
heatmap_height = 12,
width = 15,
height = 10
)
#> Warning: Removed 8568 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 1326 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 1248 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 822 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 582 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 423 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 240 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 225 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 177 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 141 rows containing non-finite outside the scale range
#> (`stat_summary()`).
plot_modules_h(gene_module %>% filter(Module != 0),
data_obj_merged, scale = TRUE,
ylabel = "Z-score",
go_list = go_list_wgcna_merged$all,
go_category = "CC",
fontsize = 7,
heatmap_width = 8,
heatmap_height = 12,
width = 15,
height = 10
)
#> Warning: Removed 8568 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 1326 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 1248 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 822 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 582 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 423 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 240 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 225 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 177 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 141 rows containing non-finite outside the scale range
#> (`stat_summary()`).
12.8 Hub genes
hub_genes <- WGCNA::chooseTopHubInEachModule(net$input_data, net$colors,
omitColors = 0,
power = picked_power,
type = "signed"
)
clusterProfiler::bitr(hub_genes, fromType = "SYMBOL", toType = c("GENENAME"),
OrgDb = org.Mm.eg.db)
#> 'select()' returned 1:1 mapping between keys and columns13 Session info
#> R version 4.6.0 (2026-04-24 ucrt)
#> Platform: x86_64-w64-mingw32/x64
#> Running under: Windows 11 x64 (build 26100)
#>
#> Matrix products: default
#> LAPACK version 3.12.1
#>
#> locale:
#> [1] LC_COLLATE=English_United Kingdom.utf8
#> [2] LC_CTYPE=English_United Kingdom.utf8
#> [3] LC_MONETARY=English_United Kingdom.utf8
#> [4] LC_NUMERIC=C
#> [5] LC_TIME=English_United Kingdom.utf8
#>
#> time zone: Europe/London
#> tzcode source: internal
#>
#> attached base packages:
#> [1] stats4 stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] ggplot2_4.0.3 dplyr_1.2.1
#> [3] stringr_1.6.0 GEOquery_2.79.0
#> [5] org.Mm.eg.db_3.23.0 AnnotationDbi_1.75.0
#> [7] SummarizedExperiment_1.43.0 Biobase_2.73.0
#> [9] GenomicRanges_1.65.0 Seqinfo_1.3.0
#> [11] IRanges_2.47.0 S4Vectors_0.51.0
#> [13] BiocGenerics_0.59.0 generics_0.1.4
#> [15] MatrixGenerics_1.25.0 matrixStats_1.5.0
#> [17] magrittr_2.0.5 TiDEomics_0.99.0
#> [19] BiocStyle_2.41.0
#>
#> loaded via a namespace (and not attached):
#> [1] segmented_2.2-1 bitops_1.0-9
#> [3] fs_2.1.0 enrichplot_1.33.0
#> [5] httr_1.4.8 RColorBrewer_1.1-3
#> [7] doParallel_1.0.17 ggsci_5.0.0
#> [9] dynamicTreeCut_1.63-1 tools_4.6.0
#> [11] backports_1.5.1 R6_2.6.1
#> [13] lazyeval_0.2.3 GetoptLong_1.1.1
#> [15] withr_3.0.2 gridExtra_2.3
#> [17] preprocessCore_1.75.0 WGCNA_1.74
#> [19] cli_3.6.6 scatterpie_0.2.6
#> [21] labeling_0.4.3 sass_0.4.10
#> [23] S7_0.2.2 readr_2.2.0
#> [25] pbapply_1.7-4 ggridges_0.5.7
#> [27] askpass_1.2.1 systemfonts_1.3.2
#> [29] yulab.utils_0.2.4 foreign_0.8-91
#> [31] gson_0.1.0 DOSE_4.7.0
#> [33] R.utils_2.13.0 rentrez_1.2.4
#> [35] limma_3.69.0 impute_1.87.0
#> [37] rstudioapi_0.18.0 RSQLite_2.4.6
#> [39] gridGraphics_0.5-1 shape_1.4.6.1
#> [41] gtools_3.9.5 crosstalk_1.2.2
#> [43] car_3.1-5 GO.db_3.23.1
#> [45] Matrix_1.7-5 abind_1.4-8
#> [47] PCAtools_2.25.0 R.methodsS3_1.8.2
#> [49] lifecycle_1.0.5 yaml_2.3.12
#> [51] carData_3.0-6 gplots_3.3.0
#> [53] qvalue_2.45.0 SparseArray_1.13.0
#> [55] grid_4.6.0 blob_1.3.0
#> [57] promises_1.5.0 dqrng_0.4.1
#> [59] crayon_1.5.3 ggtangle_0.1.2
#> [61] lattice_0.22-9 beachmat_2.29.0
#> [63] cowplot_1.2.0 KEGGREST_1.53.0
#> [65] pillar_1.11.1 knitr_1.51
#> [67] ComplexHeatmap_2.29.0 boot_1.3-32
#> [69] rjson_0.2.23 codetools_0.2-20
#> [71] glue_1.8.1 ggiraph_0.9.6
#> [73] ggfun_0.2.0 fontLiberation_0.1.0
#> [75] data.table_1.18.2.1 Rdpack_2.6.6
#> [77] vctrs_0.7.3 png_0.1-9
#> [79] treeio_1.37.0 gtable_0.3.6
#> [81] cachem_1.1.0 xfun_0.57
#> [83] rbibutils_2.4.1 S4Arrays_1.13.0
#> [85] mime_0.13 reformulas_0.4.4
#> [87] survival_3.8-6 aisdk_1.1.0
#> [89] iterators_1.0.14 statmod_1.5.1
#> [91] nlme_3.1-169 ggtree_4.3.0
#> [93] bit64_4.8.0 fontquiver_0.2.1
#> [95] bslib_0.10.0 irlba_2.3.7
#> [97] rpart_4.1.27 KernSmooth_2.23-26
#> [99] otel_0.2.0 Hmisc_5.2-5
#> [101] colorspace_2.1-2 DBI_1.3.0
#> [103] nnet_7.3-20 tidyselect_1.2.1
#> [105] processx_3.9.0 bit_4.6.0
#> [107] compiler_4.6.0 curl_7.1.0
#> [109] httr2_1.2.2 htmlTable_2.5.0
#> [111] xml2_1.5.2 randtests_1.0.2
#> [113] fontBitstreamVera_0.1.1 DelayedArray_0.39.0
#> [115] plotly_4.12.0 bookdown_0.46
#> [117] checkmate_2.3.4 caTools_1.18.3
#> [119] scales_1.4.0 callr_3.7.6
#> [121] rappdirs_0.3.4 digest_0.6.39
#> [123] minqa_1.2.8 rmarkdown_2.31
#> [125] XVector_0.53.0 base64enc_0.1-6
#> [127] htmltools_0.5.9 pkgconfig_2.0.3
#> [129] lme4_2.0-1 umap_0.2.10.0
#> [131] sparseMatrixStats_1.25.0 fastmap_1.2.0
#> [133] rlang_1.2.0 GlobalOptions_0.1.4
#> [135] htmlwidgets_1.6.4 shiny_1.13.0
#> [137] DelayedMatrixStats_1.35.0 ggh4x_0.3.1
#> [139] farver_2.1.2 jquerylib_0.1.4
#> [141] jsonlite_2.0.0 BiocParallel_1.47.0
#> [143] GOSemSim_2.39.0 R.oo_1.27.1
#> [145] BiocSingular_1.29.0 Formula_1.2-5
#> [147] ggplotify_0.1.3 patchwork_1.3.2
#> [149] Rcpp_1.1.1-1.1 ape_5.8-1
#> [151] ggnewscale_0.5.2 gdtools_0.5.0
#> [153] reticulate_1.46.0 stringi_1.8.7
#> [155] MASS_7.3-65 plyr_1.8.9
#> [157] shinyFiles_0.9.3 parallel_4.6.0
#> [159] ggrepel_0.9.8 Biostrings_2.81.0
#> [161] splines_4.6.0 hms_1.1.4
#> [163] circlize_0.4.18 fastcluster_1.3.0
#> [165] igraph_2.3.0 ggpubr_0.6.3
#> [167] ggsignif_0.6.4 enrichit_0.1.4
#> [169] reshape2_1.4.5 ScaledMatrix_1.21.0
#> [171] XML_3.99-0.23 evaluate_1.0.5
#> [173] BiocManager_1.30.27 nloptr_2.2.1
#> [175] tzdb_0.5.0 foreach_1.5.2
#> [177] tweenr_2.0.3 httpuv_1.6.17
#> [179] tidyr_1.3.2 openssl_2.4.0
#> [181] purrr_1.2.2 polyclip_1.10-7
#> [183] clue_0.3-68 Trendy_1.35.0
#> [185] ggforce_0.5.0 rsvd_1.0.5
#> [187] xtable_1.8-8 broom_1.0.12
#> [189] RSpectra_0.16-2 tidytree_0.4.7
#> [191] tidydr_0.0.6 rstatix_0.7.3
#> [193] later_1.4.8 viridisLite_0.4.3
#> [195] tibble_3.3.1 clusterProfiler_4.21.0
#> [197] aplot_0.2.9 memoise_2.0.1
#> [199] cluster_2.1.8.2
Bacher, Rhonda, Ning Leng, Li-Fang Chu, et al. 2018. “Trendy: Segmented Regression Analysis of Expression Dynamics in High-Throughput Ordered Profiling Experiments.” BMC Bioinformatics. https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-018-2405-x.
Gu, Zuguang. 2022. “Complex Heatmap Visualization.” iMeta, ahead of print. https://doi.org/10.1002/imt2.43.
Gu, Zuguang, Roland Eils, and Matthias Schlesner. 2016. “Complex Heatmaps Reveal Patterns and Correlations in Multidimensional Genomic Data.” Bioinformatics, ahead of print. https://doi.org/10.1093/bioinformatics/btw313.
Huber, W., Carey, et al. 2015. “Orchestrating High-Throughput Genomic Analysis with Bioconductor.” Nature Methods 12 (2): 115–21. http://www.nature.com/nmeth/journal/v12/n2/full/nmeth.3252.html.
Langfelder, Peter, and Steve Horvath. 2008. “WGCNA: An r Package for Weighted Correlation Network Analysis.” BMC Bioinformatics, no. 1: 559. https://link.springer.com/article/10.1186/1471-2105-9-559.
Langfelder, Peter, and Steve Horvath. 2012. “Fast R Functions for Robust Correlations and Hierarchical Clustering.” Journal of Statistical Software 46 (11): 1–17. https://www.jstatsoft.org/v46/i11/.
Ritchie, Matthew E, Belinda Phipson, Di Wu, et al. 2015. “limma Powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies.” Nucleic Acids Research 43 (7): e47. https://doi.org/10.1093/nar/gkv007.
Traxler, Peter, Stephan Reichl, Lukas Folkman, et al. 2025. “Integrated Time-Series Analysis and High-Content CRISPR Screening Delineate the Dynamics of Macrophage Immune Regulation.” Cell Systems 16 (8): 101346. https://doi.org/https://doi.org/10.1016/j.cels.2025.101346.
Vasaikar, Suhas V, Adam K Savage, Qiuyu Gong, et al. 2023. “A Comprehensive Platform for Analyzing Longitudinal Multi-Omics Data.” Nature Communications 14 (1): 1684. https://www.nature.com/articles/s41467-023-37432-w.
Wickham, Hadley. 2016. Ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. https://ggplot2.tidyverse.org.
Wu, Tianzhi, Erqiang Hu, Shuangbin Xu, et al. 2021. “clusterProfiler 4.0: A Universal Enrichment Tool for Interpreting Omics Data.” The Innovation 2 (3): 100141. https://doi.org/10.1016/j.xinn.2021.100141.
Xu, Shuangbin, Erqiang Hu, Yantong Cai, et al. 2024. “Using clusterProfiler to Characterize Multiomics Data.” Nature Protocols 19 (11): 3292–320. https://doi.org/10.1038/s41596-024-01020-z.
Yu, Guangchuang. 2024. “Thirteen Years of clusterProfiler.” The Innovation 5 (6): 100722. https://doi.org/10.1016/j.xinn.2024.100722.
Yu, Guangchuang, Li-Gen Wang, Yanyan Han, and Qing-Yu He. 2012. “clusterProfiler: An r Package for Comparing Biological Themes Among Gene Clusters.” OMICS: A Journal of Integrative Biology 16 (5): 284–87. https://doi.org/10.1089/omi.2011.0118.
Zhang, Jun, Hongyuan Li, Wenjun Tao, and Jun Zhou. 2026. “ClusterGVis: An Advanced Visualization and Clustering Tool for Gene Expression Analysis.” Genomics, Proteomics & Bioinformatics, January, qzag005. https://doi.org/10.1093/gpbjnl/qzag005.