TiDEomics Application: Macrophage response to six immune stimuli
27 July 2026
Package
TiDEomics 0.99.4
1 Introduction
TiDEomics is an R package designed to streamline analysis of time-course omics data, compatible with multiple group comparison and 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
Get the latest stable R release from CRAN. Then install TiDEomics from Bioconductor using the following code:
if (!requireNamespace("BiocManager", quietly = TRUE)) {
install.packages("BiocManager")
}
BiocManager::install("TiDEomics")The development version can be installed from GitHub with:
if (!require("remotes", quietly = TRUE)) install.packages("remotes")
remotes::install_github("hte123/TiDEomics")library(TiDEomics)
#> Warning: multiple methods tables found for 'transform'
#> Warning: replacing previous import 'BiocGenerics::transform' by
#> 'S4Vectors::transform' when loading 'SummarizedExperiment'
#> Warning: replacing previous import 'BiocGenerics::transform' by
#> 'S4Vectors::transform' when loading 'IRanges'
#> Warning: multiple methods tables found for 'transform'
#> Warning: replacing previous import 'BiocGenerics::transform' by
#> 'S4Vectors::transform' when loading 'Seqinfo'
#> Warning: replacing previous import 'BiocGenerics::transform' by
#> 'S4Vectors::transform' when loading 'GenomicRanges'
#> Warning: replacing previous import 'BiocGenerics::transform' by
#> 'S4Vectors::transform' when loading 'S4Arrays'
#> Warning: replacing previous import 'BiocGenerics::transform' by
#> 'S4Vectors::transform' when loading 'DelayedArray'
#> Warning: replacing previous import 'BiocGenerics::transform' by
#> 'S4Vectors::transform' when loading 'SparseArray'
#> Warning: replacing previous import 'BiocGenerics::transform' by
#> 'S4Vectors::transform' when loading 'XVector'
#> Warning: replacing previous import 'BiocGenerics::transform' by
#> 'S4Vectors::transform' when loading 'AnnotationDbi'
#> Warning: replacing previous import 'BiocGenerics::transform' by
#> 'S4Vectors::transform' when loading 'Biostrings'
#>
#> Registered S3 method overwritten by 'lme4':
#> method from
#> na.action.merMod car
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
#> Warning: package 'GenomicRanges' was built under R version 4.6.1
#> Loading required package: stats4
#> Loading required package: BiocGenerics
#> Warning: package 'BiocGenerics' was built under R version 4.6.1
#> Loading required package: generics
#> Warning: package 'generics' was built under R version 4.6.1
#>
#> 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 object is masked from 'package:utils':
#>
#> data
#> 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, scale, sequence, table,
#> tapply, transform, 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
#> 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)
#> Warning: package 'GEOquery' was built under R version 4.6.1
#> Setting options('download.file.method.GEOquery'='auto')
#> Setting options('GEOquery.inmemory.gpl'=FALSE)
library(stringr)
#> Warning: package 'stringr' was built under R version 4.6.1
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)
#> Warning: package 'ggplot2' was built under R version 4.6.13 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.
- Normalised to
log2(CPM + 1)
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, returnType = "ExpressionSet")
#> 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
tmp_dir <- tempdir()
getGEOSuppFiles("GSE263759", baseDir = tmp_dir, makeDirectory = FALSE)geo_data_tb_gz <- gzfile(file.path(tmp_dir, "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(geo_data_tb$Feature, fromType = "ENSEMBL", :
#> 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))
gene_keep <- which(gene_keep != 0)
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)
#> No Subject column specified. Samples treated as independent. For repeated-measures designs, set subject_col to the column identifying biological subjects.
#> 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.
# Normalise: log2(CPM + 1)
count_mat <- as.matrix(assays(data_obj)[[1]])
lib_sizes <- colSums(count_mat)
cpm <- t(t(count_mat) / lib_sizes * 1e6)
assays(data_obj)[[1]] <- log2(cpm + 1)
data_obj#> class: SummarizedExperiment
#> dim: 20150 76
#> metadata(0):
#> assays(1): orig
#> 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 Batchlog2(CPM + 1) normalised counts are stored in the first assay slot of data_obj.
#> RNA_Candida_0h_R1_1 RNA_Candida_0h_R2_1 RNA_Candida_2h_R1_1
#> Gnai3 5.1883751 5.058662 5.3676589
#> Cdc45 4.4337057 4.374351 4.2383557
#> H19 0.0000000 0.000000 0.0000000
#> Scml2 0.6836615 1.100543 0.3603305
#> Narf 3.7793051 3.823756 3.6911141
#> Cav2 6.7987039 6.735611 6.0744728
#> RNA_Candida_2h_R2_1 RNA_Candida_4h_R1_1 RNA_Candida_4h_R2_1
#> Gnai3 5.312290 5.8648503 5.5918556
#> Cdc45 4.137515 4.4371716 3.6071559
#> H19 0.000000 0.0000000 0.0000000
#> Scml2 1.252988 0.7866156 0.3903628
#> Narf 2.782244 4.5530972 4.8406819
#> Cav2 6.734439 5.7036025 5.9805873
#> RNA_Candida_6h_R1_1 RNA_Candida_6h_R2_1 RNA_Candida_8h_R1_2
#> Gnai3 5.7446711 5.296213 6.109646
#> Cdc45 4.2047217 3.010992 4.157170
#> H19 0.0000000 0.000000 0.000000
#> Scml2 0.9275469 1.237515 0.000000
#> Narf 4.3816199 4.112287 4.314972
#> Cav2 5.9052841 6.156645 6.028413
#> RNA_Candida_8h_R2_2 RNA_Candida_24h_R1_2 RNA_Candida_24h_R2_2
#> Gnai3 6.277461 5.3190392 5.9805382
#> Cdc45 4.655239 4.4743764 4.3381671
#> H19 0.000000 0.0000000 0.4258121
#> Scml2 0.000000 0.9142163 1.4418245
#> Narf 4.107898 4.7201120 4.2627644
#> Cav2 5.906972 6.2232698 6.4171570
#> RNA_IFNbeta_0h_R1_1 RNA_IFNbeta_0h_R2_1 RNA_IFNbeta_2h_R1_1
#> Gnai3 5.1883751 5.058662 5.425859
#> Cdc45 4.4337057 4.374351 4.361180
#> H19 0.0000000 0.000000 0.000000
#> Scml2 0.6836615 1.100543 1.190345
#> Narf 3.7793051 3.823756 2.689033
#> Cav2 6.7987039 6.735611 5.988709
#> RNA_IFNbeta_2h_R2_1 RNA_IFNbeta_4h_R1_1 RNA_IFNbeta_4h_R2_1
#> Gnai3 4.039613 5.2809666 4.6223267
#> Cdc45 3.440789 3.2742427 3.3073549
#> H19 0.000000 0.0000000 0.3861257
#> Scml2 1.571315 0.9010984 0.3861257
#> Narf 2.567882 1.5965235 1.3416418
#> Cav2 6.383446 4.9973095 5.3646132
#> RNA_IFNbeta_6h_R1_1 RNA_IFNbeta_6h_R2_1 RNA_IFNbeta_8h_R1_2
#> Gnai3 5.1936994 4.7701166 5.0597588
#> Cdc45 2.7496585 2.6867449 3.0601634
#> H19 0.0000000 0.3808889 0.0000000
#> Scml2 0.5829186 0.3808889 0.8604703
#> Narf 2.9946286 3.5701206 4.7011418
#> Cav2 4.8395637 4.4293339 3.8939166
#> RNA_IFNbeta_8h_R2_2 RNA_IFNbeta_24h_R1_2 RNA_IFNbeta_24h_R2_2
#> Gnai3 5.990905 5.162419 5.703476
#> Cdc45 2.854530 1.393298 2.655409
#> H19 0.000000 0.000000 0.000000
#> Scml2 0.000000 0.000000 0.000000
#> Narf 4.503278 5.702687 5.222827
#> Cav2 3.598178 6.604452 6.429325
#> RNA_IFNgamma_0h_R1_1 RNA_IFNgamma_0h_R2_1 RNA_IFNgamma_2h_R1_1
#> Gnai3 5.1883751 5.058662 5.4450686
#> Cdc45 4.4337057 4.374351 4.3545851
#> H19 0.0000000 0.000000 0.0000000
#> Scml2 0.6836615 1.100543 0.6853203
#> Narf 3.7793051 3.823756 4.3110545
#> Cav2 6.7987039 6.735611 6.0324148
#> RNA_IFNgamma_2h_R2_1 RNA_IFNgamma_4h_R1_1 RNA_IFNgamma_4h_R2_1
#> Gnai3 5.1608970 5.101257 4.7739087
#> Cdc45 4.2679761 3.689565 3.2754164
#> H19 0.0000000 0.000000 0.0000000
#> Scml2 0.7676835 0.000000 0.3896748
#> Narf 3.6544136 4.223695 4.4651041
#> Cav2 6.3288781 5.691258 6.0265855
#> RNA_IFNgamma_6h_R1_1 RNA_IFNgamma_6h_R2_1 RNA_IFNgamma_8h_R1_2
#> Gnai3 5.6779861 5.069496 4.8881911
#> Cdc45 3.9990410 3.656769 3.3293047
#> H19 0.0000000 0.000000 0.0000000
#> Scml2 0.4314518 0.000000 0.9076049
#> Narf 4.1763952 4.892799 3.9633639
#> Cav2 4.9693497 5.096950 4.9439949
#> RNA_IFNgamma_8h_R2_2 RNA_IFNgamma_24h_R1_2 RNA_IFNgamma_24h_R2_2
#> Gnai3 5.646327 5.467303 5.2277158
#> Cdc45 4.330867 4.148800 3.8709273
#> H19 0.000000 0.000000 0.0000000
#> Scml2 0.000000 0.763353 0.7964273
#> Narf 4.523119 4.564775 4.5974737
#> Cav2 4.937470 5.501021 5.8420939
#> RNA_LCMV_0h_R1_1 RNA_LCMV_0h_R2_1 RNA_LCMV_2h_R1_1 RNA_LCMV_2h_R2_1
#> Gnai3 5.1883751 5.058662 4.5895819 5.503052
#> Cdc45 4.4337057 4.374351 3.7931232 4.666351
#> H19 0.0000000 0.000000 0.0000000 0.000000
#> Scml2 0.6836615 1.100543 0.4629014 1.180887
#> Narf 3.7793051 3.823756 3.5814740 3.900990
#> Cav2 6.7987039 6.735611 6.2675453 7.007493
#> RNA_LCMV_4h_R1_1 RNA_LCMV_4h_R2_1 RNA_LCMV_6h_R1_1 RNA_LCMV_6h_R2_1
#> Gnai3 6.1388422 3.871637 5.938218 5.489807
#> Cdc45 4.7684290 2.507929 4.695637 4.269036
#> H19 0.0000000 0.000000 0.000000 0.000000
#> Scml2 0.3487827 1.188118 1.568938 1.698182
#> Narf 4.2308614 3.433259 4.561871 4.244837
#> Cav2 6.1444307 6.264177 6.391219 6.452689
#> RNA_LCMV_8h_R1_2 RNA_LCMV_8h_R2_2 RNA_LCMV_24h_R1_2 RNA_LCMV_24h_R2_2
#> Gnai3 5.921116 6.177433 5.6845152 6.1079876
#> Cdc45 4.547883 4.335614 3.1826626 2.6952377
#> H19 0.000000 0.000000 0.0000000 0.0000000
#> Scml2 0.000000 1.424056 0.6389806 0.4028894
#> Narf 4.132158 4.156549 3.6349631 3.2777076
#> Cav2 6.155291 5.999069 4.9190776 5.6241882
#> RNA_Listeria_0h_R1_1 RNA_Listeria_0h_R2_1 RNA_Listeria_2h_R1_1
#> Gnai3 5.1883751 5.058662 5.5544422
#> Cdc45 4.4337057 4.374351 4.5640345
#> H19 0.0000000 0.000000 0.0000000
#> Scml2 0.6836615 1.100543 0.7540288
#> Narf 3.7793051 3.823756 2.0318297
#> Cav2 6.7987039 6.735611 6.1149154
#> RNA_Listeria_2h_R2_1 RNA_Listeria_4h_R1_1 RNA_Listeria_4h_R2_1
#> Gnai3 5.6680569 5.2728133 4.463893
#> Cdc45 4.0910506 3.2939375 2.020868
#> H19 0.0000000 0.0000000 0.000000
#> Scml2 0.8948453 0.6847633 0.000000
#> Narf 2.5533587 2.1179975 1.013945
#> Cav2 6.7022317 5.3064200 5.843186
#> RNA_Listeria_6h_R1_1 RNA_Listeria_6h_R2_1 RNA_Listeria_8h_R1_2
#> Gnai3 4.967871 4.694971 5.427805
#> Cdc45 2.501007 2.407076 1.348123
#> H19 0.000000 0.000000 0.000000
#> Scml2 0.000000 0.000000 0.000000
#> Narf 2.288065 2.407076 3.125822
#> Cav2 4.776695 4.945817 4.520451
#> RNA_Listeria_8h_R2_2 RNA_Listeria_24h_R1_2 RNA_Listeria_24h_R2_2
#> Gnai3 6.097996 5.9144973 6.494620
#> Cdc45 2.587278 0.6719019 1.290451
#> H19 0.000000 0.0000000 0.000000
#> Scml2 0.000000 0.0000000 0.000000
#> Narf 3.714641 3.3070756 3.994803
#> Cav2 4.516281 5.6160907 6.147555
#> RNA_LPS_0h_R1_1 RNA_LPS_0h_R2_1 RNA_LPS_2h_R1_1 RNA_LPS_2h_R2_1
#> Gnai3 5.1883751 5.058662 5.1356768 5.1367760
#> Cdc45 4.4337057 4.374351 4.5129426 3.7885130
#> H19 0.0000000 0.000000 0.4352212 0.0000000
#> Scml2 0.6836615 1.100543 1.0400889 0.8489156
#> Narf 3.7793051 3.823756 2.1766998 2.6647052
#> Cav2 6.7987039 6.735611 6.2072770 6.6019578
#> RNA_LPS_4h_R1_1 RNA_LPS_4h_R2_1 RNA_LPS_6h_R1_1 RNA_LPS_6h_R2_1
#> Gnai3 5.131361 5.2246467 4.324832 4.8608941
#> Cdc45 3.227792 3.1101336 1.301376 1.2481841
#> H19 0.000000 0.0000000 0.000000 0.0000000
#> Scml2 0.000000 0.8184109 0.000000 0.6324255
#> Narf 2.061498 1.9257574 3.235983 3.0750672
#> Cav2 5.119434 5.8196242 4.735044 4.6557773
#> RNA_LPS_8h_R1_2 RNA_LPS_8h_R2_2 RNA_LPS_24h_R2_2 RNA_untreated_0h_R1_1
#> Gnai3 5.0724141 6.069472 5.4663065 5.1883751
#> Cdc45 1.7011189 1.534551 0.7876175 4.4337057
#> H19 0.0000000 0.000000 0.0000000 0.0000000
#> Scml2 0.6442094 0.000000 0.0000000 0.6836615
#> Narf 4.4665029 4.496857 4.9265673 3.7793051
#> Cav2 4.2488568 3.994466 6.2715237 6.7987039
#> RNA_untreated_0h_R2_1 RNA_untreated_8h_R2_2 RNA_untreated_24h_R1_2
#> Gnai3 5.058662 5.962002 5.9745566
#> Cdc45 4.374351 4.677723 4.7932730
#> H19 0.000000 0.000000 0.0000000
#> Scml2 1.100543 1.311244 0.3679547
#> Narf 3.823756 3.841368 4.5465979
#> Cav2 6.735611 6.048506 6.0204765
#> RNA_untreated_24h_R2_2
#> Gnai3 5.6846013
#> Cdc45 4.4093442
#> H19 0.0000000
#> Scml2 0.8189527
#> Narf 4.2446539
#> Cav2 6.7184035#> 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.35
#> Picking joint bandwidth of 0.343
#> Picking joint bandwidth of 0.348
#> Picking joint bandwidth of 0.349
#> Picking joint bandwidth of 0.34
#> Picking joint bandwidth of 0.341
#> Picking joint bandwidth of 0.351
5 Data processing
Normalise to time point 0
data_obj <- normalise_to_start(data_obj)
#> Normalising to group baseline at each feature's first non-NA time point.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): orig norm
#> rownames(20150): Gnai3 Cdc45 ... Gm17315 Ube2srt
#> 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() accepts either a SummarizedExperiment object directly (with assay parameter) or a pre-computed table from calc_mean_sd().
# Using pre-computed table
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
# Using SummarizedExperiment directly
plot_trend(data_obj, assay = 1,
features = c("Isg15", "Cd68", "Adgre1", "Tnf"),
title = "Example genes (from SE object)", ylab = "Log2 (count + 1)"
)
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
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
6.2.1.2 Per group
plot_umap_by_group(data_obj, nrow = 2, seed = 123)
#> Using n_neighbors = 5
#> Warning: Using size for a discrete variable is not advised.
#> Using n_neighbors = 5
#> Warning: Using size for a discrete variable is not advised.
#> Using n_neighbors = 5
#> Warning: Using size for a discrete variable is not advised.
#> Using n_neighbors = 5
#> Warning: Using size for a discrete variable is not advised.
#> Using n_neighbors = 5
#> Warning: Using size for a discrete variable is not advised.
#> Using n_neighbors = 5
#> Warning: Using size for a discrete variable is not advised.
#> Using n_neighbors = 4
#> Warning: Using size for a discrete variable is not advised.
6.2.2 Replicates merged by mean
6.2.2.1 All groups
6.2.2.2 Per group
plot_umap_by_group(data_obj_merged, nrow = 2)
#> Using n_neighbors = 5
#> Warning: Using size for a discrete variable is not advised.
#> Using n_neighbors = 5
#> Warning: Using size for a discrete variable is not advised.
#> Using n_neighbors = 5
#> Warning: Using size for a discrete variable is not advised.
#> Using n_neighbors = 5
#> Warning: Using size for a discrete variable is not advised.
#> Using n_neighbors = 5
#> Warning: Using size for a discrete variable is not advised.
#> Using n_neighbors = 5
#> Warning: Using size for a discrete variable is not advised.
#> Using n_neighbors = 2
#> Warning: failed creating initial embedding; using random embedding instead
#> Warning: Using size for a discrete variable is not advised.
6.3 PCA
using original data not normalised to time 0
6.3.1 Original data with replicates
6.3.1.1 All groups
pca_res <- plot_pca(data_obj,
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.
PC <- pca_res$pcaplot_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: Using size for a discrete variable is not advised.
#> Using size for a discrete variable is not advised.
#> Using size for a discrete variable is not advised.
#> Using size for a discrete variable is not advised.
#> Using size for a discrete variable is not advised.
#> Using size for a discrete variable is not advised.
#> Using size for a discrete variable is not advised.
6.3.2 Replicates merged by mean
6.3.2.1 All groups
pca_res_merged <- plot_pca(data_obj_merged,
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.
PC_merged <- pca_res_merged$pcaplot_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(data_obj_merged_list, nrow = 2, circle = FALSE)
#> Warning: Using size for a discrete variable is not advised.
#> Using size for a discrete variable is not advised.
#> Using size for a discrete variable is not advised.
#> Using size for a discrete variable is not advised.
#> Using size for a discrete variable is not advised.
#> Using size for a discrete variable is not advised.
#> Using size for a discrete variable is not advised.
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)) {
print(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
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
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"
)#> $dotplot_BP
#>
#> $dotplot_MF
#>
#> $dotplot_CC
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+, PC1_top-, PC2_top+, PC2_top-, PC3_top+
#> The following groups are missing in the custom palette: PC1_top+, PC1_top-, PC2_top-, PC3_top+#> $dotplot_BP
#>
#> $dotplot_MF
#>
#> $dotplot_CC
#>
#> $cnetplot_BP
#>
#> $cnetplot_MF
#>
#> $cnetplot_CC
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+, 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+, PC1_top-, PC2_top-, PC3_top+#> $dotplot_BP
#>
#> $dotplot_MF
#>
#> $dotplot_CC
#>
#> $emapplot_BP
#>
#> $emapplot_MF
#>
#> $emapplot_CC
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)) {
print(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
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
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)"
)#> $dotplot_BP
#>
#> $dotplot_MF
#>
#> $dotplot_CC
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+, PC3_top-
#> The following groups are missing in the custom palette: PC1_top+, PC1_top-, PC2_top+, PC2_top-, PC3_top+, PC3_top-
#> The following groups are missing in the custom palette: PC1_top+, PC1_top-, PC2_top-, PC3_top+#> $dotplot_BP
#>
#> $dotplot_MF
#>
#> $dotplot_CC
#>
#> $cnetplot_BP
#>
#> $cnetplot_MF
#>
#> $cnetplot_CC
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+, 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+, PC1_top-, PC2_top+, PC2_top-, PC3_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+, PC1_top-, PC2_top-, PC3_top+#> $dotplot_BP
#>
#> $dotplot_MF
#>
#> $dotplot_CC
#>
#> $emapplot_BP
#>
#> $emapplot_MF
#>
#> $emapplot_CC
8 Pairwise differential expression (DE)
8.1 Between time points within each group
Default DE criteria: p.adj < 0.05, |log2FC| > 1.
trend = TRUE is used as recommended for RNA-seq count-derived data to model the mean-variance relationship.
DE_between_time_out <- DE_between_time(data_obj, assay = 1, filter = 1, trend = TRUE)
#> 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(DE_between_time_out,
fontsize = 8, value = FALSE, nrow = 1, heatmap_width = 3
)
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. trend = TRUE for RNA-seq.
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", trend = TRUE)
#> Non-NA replicate number filter not specified. Using minimum number of replicates across all groups and time points: 0
#> Comparing Candida vs untreated: time points only in untreated: none; only in Candida: 2, 4, 6
#> 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 IFNbeta vs untreated: time points only in untreated: none; only in IFNbeta: 2, 4, 6
#> 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 IFNgamma vs untreated: time points only in untreated: none; only in IFNgamma: 2, 4, 6
#> 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 LCMV vs untreated: time points only in untreated: none; only in LCMV: 2, 4, 6
#> 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 Listeria vs untreated: time points only in untreated: none; only in Listeria: 2, 4, 6
#> 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 LPS vs untreated: time points only in untreated: none; only in LPS: 2, 4, 6
#> 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 zeroplot_volcano(DE_between_group_out,
group1 = "untreated", group2 = "IFNgamma", time = 24,
logFC_thres = 3, adjP_thres = 0.05, label = TRUE)
The group parameter of plot_DE_between_group restricts the plotting to
comparisons versus a specific group.
#> $untreated
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] 119
#>
#> $IFNbeta
#> [1] 546
#>
#> $IFNgamma
#> [1] 389
#>
#> $LCMV
#> [1] 470
#>
#> $Listeria
#> [1] 638
#>
#> $LPS
#> [1] 401
#>
#> $untreated
#> [1] 0
non_random_down_list |> lapply(length)#> $Candida
#> [1] 85
#>
#> $IFNbeta
#> [1] 1863
#>
#> $IFNgamma
#> [1] 535
#>
#> $LCMV
#> [1] 556
#>
#> $Listeria
#> [1] 2791
#>
#> $LPS
#> [1] 1447
#>
#> $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
print(plot_trend(table_mean_sd_orig,
groups = i,
features = random, title = "Overall down DEG"
))
}
9.4 GO enrichment of overall DEGs
Simplified enrichment can be customised with parameters: simplify_cutoff, simplify_by, simplify_select_fun, simplify_measure.
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
#> Processing gene list: IFNbeta
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: IFNgamma
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: LCMV
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: Listeria
#> 'select()' returned 1:many mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: LPS
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Gene list untreated is empty. Skipping.
#> Simplifying GO terms for category: BP
#> Warning: semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> 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
#> Processing gene list: IFNbeta
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: IFNgamma
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: LCMV
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: Listeria
#> 'select()' returned 1:many mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: LPS
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Gene list untreated is empty. Skipping.
#> Simplifying GO terms for category: MF
#> Warning: semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> 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
#> Processing gene list: IFNbeta
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: IFNgamma
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: LCMV
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: Listeria
#> 'select()' returned 1:many mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: LPS
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Gene list untreated is empty. Skipping.
#> Simplifying GO terms for category: CC
#> Warning: semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> 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
#> Processing gene list: IFNbeta
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: IFNgamma
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: LCMV
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: Listeria
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: LPS
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Gene list untreated is empty. Skipping.
#> Simplifying GO terms for category: BP
#> Warning: semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> 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
#> Processing gene list: IFNbeta
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: IFNgamma
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: LCMV
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: Listeria
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: LPS
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Gene list untreated is empty. Skipping.
#> Simplifying GO terms for category: MF
#> Warning: semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> 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
#> Processing gene list: IFNbeta
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: IFNgamma
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: LCMV
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: Listeria
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: LPS
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Gene list untreated is empty. Skipping.
#> Simplifying GO terms for category: CC
#> Warning: semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> 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"
)#> $dotplot_BP
#>
#> $dotplot_MF
#>
#> $dotplot_CC
plot_GO(go_list_up_merged$simplified,
plot_dotplot = TRUE,
showCategory_dotplot = 5,
label = "overall up DEGs (simplified)"
)#> $dotplot_BP
#>
#> $dotplot_MF
#>
#> $dotplot_CC
plot_GO(go_list_down_merged$all,
plot_dotplot = TRUE,
showCategory_dotplot = 5,
label = "overall down DEGs"
)#> $dotplot_BP
#>
#> $dotplot_MF
#>
#> $dotplot_CC
plot_GO(go_list_down_merged$simplified,
plot_dotplot = TRUE,
showCategory_dotplot = 5,
label = "overall down DEGs (simplified)"
)#> $dotplot_BP
#>
#> $dotplot_MF
#>
#> $dotplot_CC
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.#> $dotplot_BP
#>
#> $dotplot_MF
#>
#> $dotplot_CC
#>
#> $emapplot_BP
#>
#> $emapplot_MF
#>
#> $emapplot_CC
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.#> $dotplot_BP
#>
#> $dotplot_MF
#>
#> $dotplot_CC
#>
#> $emapplot_BP
#>
#> $emapplot_MF
#>
#> $emapplot_CC
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",
OrgDb = org.Mm.eg.db)$features
}
#> 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
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)
#> Group Candida: no missing values.
#> Group IFNbeta: no missing values.
#> Group IFNgamma: no missing values.
#> Group LCMV: no missing values.
#> Group Listeria: no missing values.
#> Group LPS: no missing values.
#> Group untreated: no missing values.
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.
run_Trendy() accepts optional group and minExp parameters for finer control over which groups to analyse and expression filtering.
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(file.path("output", "TiDEomics_application_1_trendy_results_k1_2026-07-01.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 (CPM + 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 179)
#> Warning in (function (..., deparse.level = 1) : number of columns of result is
#> not a multiple of vector length (arg 40)
#> Warning in (function (..., deparse.level = 1) : number of columns of result is
#> not a multiple of vector length (arg 153)
#> Warning in (function (..., deparse.level = 1) : number of columns of result is
#> not a multiple of vector length (arg 33)
#> Warning in (function (..., deparse.level = 1) : number of columns of result is
#> not a multiple of vector length (arg 105)
#> Warning in (function (..., deparse.level = 1) : number of columns of result is
#> not a multiple of vector length (arg 70)
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 179)
#> Warning in (function (..., deparse.level = 1) : number of columns of result is
#> not a multiple of vector length (arg 40)
#> Warning in (function (..., deparse.level = 1) : number of columns of result is
#> not a multiple of vector length (arg 153)
#> Warning in (function (..., deparse.level = 1) : number of columns of result is
#> not a multiple of vector length (arg 33)
#> Warning in (function (..., deparse.level = 1) : number of columns of result is
#> not a multiple of vector length (arg 105)
#> Warning in (function (..., deparse.level = 1) : number of columns of result is
#> not a multiple of vector length (arg 70)
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]])) {
print(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
)$features
#> 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,
interaction = TRUE,
assay = 1, core = 16
)
#> Max 2 replicate(s): interaction variance estimate may be unstable.
#> LMM: exp ~ (1|Group) + (1|Time) + (1|Group:Time) | Output: Group, Time, Group:Time, Residual
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
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
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
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
plot_trend(table_mean_sd_orig,
features = var_decomp_1 |> arrange(desc(Group)) |> head(30) |> pull(Feature),
title = "Top 30 genes with high group variance"
)
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
plot_trend(table_mean_sd_orig,
features = var_decomp_1 |> arrange(desc(Time)) |> head(30) |> pull(Feature),
title = "Top 30 genes with high time variance"
)
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
plot_trend(table_mean_sd_orig,
features = var_decomp_1 |> arrange(desc(`Group:Time`)) |> head(30) |> pull(Feature),
title = "Top 30 genes with high Group:Time variance"
)
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
11.2.2 GO enrichment of genes ranked by group/time-contributed variance
enrichGO_rank() returns a named list indexed by GO category (e.g. $BP, $MF, $CC). Access results with $BP before passing to enrichplot::gseaplot2().
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.49% 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, weight = weight, :
#> There were 16 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 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 for BP.
enrichplot::gseaplot2(gse_group[["BP"]], geneSetID = 1:3, base_size = 8)
gse_time <- 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.51% 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, weight = weight, :
#> There were 21 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, weight = weight, :
#> 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 for BP.
enrichplot::gseaplot2(gse_time[["BP"]], geneSetID = 1:3, base_size = 8)
gse_grouptime <- enrichGO_rank(var_decomp_1,
gene_rank_by = "Group: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.47% 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, weight = weight, :
#> There were 20 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 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 for BP.
enrichplot::gseaplot2(gse_grouptime[["BP"]], 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,
interaction = TRUE,
assay = 2, core = 16
)
#> Max 2 replicate(s): interaction variance estimate may be unstable.
#> LMM: exp ~ (1|Group) + (1|Time) + (1|Group:Time) | Output: Group, Time, Group:Time, Residual
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.3.1 Residual variance quantified data quality
plot_trend(table_mean_sd_orig,
features = var_decomp_2 |> 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
plot_trend(table_mean_sd_orig,
features = var_decomp_2 |> 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
plot_trend(table_mean_sd_orig,
features = var_decomp_2 |> arrange(desc(Group)) |> head(30) |> pull(Feature),
title = "Top 30 genes with high group variance"
)
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
plot_trend(table_mean_sd_orig,
features = var_decomp_2 |> arrange(desc(Time)) |> head(30) |> pull(Feature),
title = "Top 30 genes with high time variance"
)
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
plot_trend(table_mean_sd_orig,
features = var_decomp_2 |> arrange(desc(`Group:Time`)) |> head(30) |> pull(Feature),
title = "Top 30 genes with high Group:Time variance"
)
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
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.
#> 0.8112 16.5310 34.9845 41.6242 64.6950 100.0000
var_res_q3_1 <- quantile(var_decomp_1$Residual, 0.75)
var_decomp_2$Residual |> summary()#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.8112 16.5310 34.9845 41.6242 64.6950 100.0000
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 identification 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] 7831include 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] 6763
data_obj_filtered <- data_obj[filter_wgcna, ]
dim(data_obj_filtered)#> [1] 6763 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 6763 of 6763
#> Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
#> 1 1 7.48e-01 2.82000 0.871 3920.0 4200.00 4690
#> 2 2 2.50e-01 0.66600 0.419 2580.0 2830.00 3700
#> 3 3 4.39e-05 0.00665 -0.202 1830.0 2000.00 3060
#> 4 4 1.19e-01 -0.31700 -0.133 1360.0 1450.00 2600
#> 5 5 3.02e-01 -0.51700 0.152 1050.0 1080.00 2240
#> 6 6 4.28e-01 -0.65200 0.386 821.0 813.00 1950
#> 7 7 5.07e-01 -0.76300 0.529 656.0 618.00 1720
#> 8 8 5.60e-01 -0.85000 0.627 532.0 476.00 1520
#> 9 9 5.89e-01 -0.93600 0.687 437.0 369.00 1350
#> 10 10 6.25e-01 -1.00000 0.741 362.0 288.00 1210
#> 11 11 6.54e-01 -1.06000 0.781 303.0 227.00 1090
#> 12 12 6.72e-01 -1.12000 0.812 256.0 180.00 983
#> 13 13 6.82e-01 -1.17000 0.828 217.0 145.00 890
#> 14 14 7.03e-01 -1.20000 0.855 186.0 118.00 809
#> 15 15 7.06e-01 -1.26000 0.864 160.0 96.70 738
#> 16 16 7.21e-01 -1.29000 0.878 138.0 79.70 674
#> 17 17 7.32e-01 -1.32000 0.893 120.0 66.10 618
#> 18 18 7.28e-01 -1.36000 0.892 105.0 55.30 567
#> 19 19 7.24e-01 -1.40000 0.893 91.5 46.50 522
#> 20 20 7.32e-01 -1.42000 0.900 80.5 39.20 482
#> 21 21 7.36e-01 -1.45000 0.903 71.0 33.20 445
#> 22 22 7.49e-01 -1.46000 0.913 62.9 28.30 412
#> 23 23 7.60e-01 -1.47000 0.924 55.8 23.90 382
#> 24 24 7.59e-01 -1.51000 0.923 49.7 20.50 355
#> 25 25 7.67e-01 -1.51000 0.930 44.4 17.50 330
#> 26 26 7.64e-01 -1.53000 0.926 39.8 15.10 307
#> 27 27 7.74e-01 -1.54000 0.933 35.7 13.10 287
#> 28 28 7.67e-01 -1.57000 0.928 32.1 11.30 268
#> 29 29 7.75e-01 -1.57000 0.933 29.0 9.78 250
#> 30 30 7.83e-01 -1.58000 0.937 26.2 8.46 234
wgcna_input$plot
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 593 genes from module 1 because their KME is too low.
#> ..removing 214 genes from module 2 because their KME is too low.
#> ..removing 292 genes from module 3 because their KME is too low.
#> ..removing 242 genes from module 4 because their KME is too low.
#> ..removing 213 genes from module 5 because their KME is too low.
#> ..removing 128 genes from module 6 because their KME is too low.
#> ..removing 78 genes from module 7 because their KME is too low.
#> ..removing 40 genes from module 8 because their KME is too low.
#> ..removing 4 genes from module 9 because their KME is too low.
#> ..removing 3 genes from module 12 because their KME is too low.
#> ..removing 20 genes from module 14 because their KME is too low.
#> ..removing 29 genes from module 15 because their KME is too low.
#> ..removing 48 genes from module 16 because their KME is too low.
#> ..removing 15 genes from module 17 because their KME is too low.
#> ..removing 16 genes from module 19 because their KME is too low.
#> ..removing 4 genes from module 20 because their KME is too low.
#> ..removing 1 genes from module 22 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...
#> 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
WGCNA_module() converts the WGCNA network output to a standardised two-column data.frame (Feature, Module), accepted by all downstream plotting and enrichment functions. It optionally excludes the grey (unassigned) module.
gene_module <- WGCNA_module(net)
gene_module |>
dplyr::group_by(Module) |>
dplyr::summarise(n = n())12.4 Module quality
summarise_module_metrics() computes per-module statistics from the WGCNA network:
module size, proportion of features assigned, and module membership (kME) statistics.
MeanKME is the mean module membership (kME = cor(feature, module eigengene));
MeanKME2 = mean(kME^2) is the average proportion of per-feature variance explained
by the module eigengene.
12.5 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)))
print(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
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
#> Group not specified. Plotting all groups: Candida, IFNbeta, IFNgamma, LCMV, Listeria, LPS, untreated
12.6 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, down_stable, down_up
#> Module 2: NA, up_down, NA, NA, NA, NA
#> Module 3: NA, down_up, NA, NA, NA, NA
#> Module 4: NA, NA, NA, down, down, down
#> Module 5: NA, NA, NA, up, up_stable, up_down
#> Module 6: NA, NA, up_stable, NA, NA, NA
#> Module 7: NA, NA, NA, up, NA, NA
#> Module 8: NA, NA, NA, NA, up, NA12.7 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)
}enrichGO_list() accepts both named lists and data.frames directly from WGCNA_module().
go_list_wgcna <- enrichGO_list(
# gene_list = module_wgcna,
gene_list = gene_module |> filter(Module != 0),
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: 1
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: 2
#> 'select()' returned 1:many mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: 3
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: 4
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: 5
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: 6
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: 7
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> 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 for category: BP
#> Warning: semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> Performing GO enrichment for category: MF
#> Processing gene list: 1
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: 2
#> 'select()' returned 1:many mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: 3
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: 4
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: 5
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: 6
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: 7
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> 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 for category: MF
#> Warning: semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> Performing GO enrichment for category: CC
#> Processing gene list: 1
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: 2
#> 'select()' returned 1:many mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: 3
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: 4
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: 5
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: 6
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> Processing gene list: 7
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned 1:many mapping between keys and columns
#> 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 for category: CC
#> Warning: semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> semData is not provided. It will be calculated automatically.
#> Merging GO enrichment results across gene lists for each category.#> $dotplot_BP
#>
#> $dotplot_MF
#>
#> $dotplot_CC
plot_GO(go_list_wgcna$simplified,
plot_dotplot = TRUE,
showCategory_dotplot = 5,
label = "WGCNA modules (simplified)"
)#> $dotplot_BP
#>
#> $dotplot_MF
#>
#> $dotplot_CC
plot_GO(go_list_wgcna$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, 8
#> The following groups are missing in the custom palette: 1, 2, 3, 4, 5, 7, 8
#> The following groups are missing in the custom palette: 1, 3, 4, 5, 7, 8#> $dotplot_BP
#>
#> $dotplot_MF
#>
#> $dotplot_CC
#>
#> $cnetplot_BP
#>
#> $cnetplot_MF
#>
#> $cnetplot_CC
plot_GO(go_list_wgcna$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, 8
#>
#> 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, 7, 8
#>
#> Coordinate system already present.
#> ℹ Adding new coordinate system, which will replace the existing one.
#> The following groups are missing in the custom palette: 1, 3, 4, 5, 7, 8#> $dotplot_BP
#>
#> $dotplot_MF
#>
#> $dotplot_CC
#>
#> $emapplot_BP
#>
#> $emapplot_MF
#>
#> $emapplot_CC
12.8 Visualise module profiles (heatmap, mean profile, GO enrichment)
12.8.1 Vertical alignment
unscaled (same as the input)
modules_v_unscaled <- plot_modules_v(gene_module |> filter(Module != 0),
data_obj_merged, assay = 2, scale = FALSE,
ylabel = "Log2 (CPM+1)\nnormalised to Time 0",
height_ratio = 3
)
modules_v_unscaled
#> Warning: Removed 12393 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, assay = 2, scale = TRUE,
ylabel = "Z-score of log2 (CPM+1)\nnormalised to Time 0",
height_ratio = 3
)
modules_v_scaled
#> Warning: Removed 12393 rows containing non-finite outside the scale range
#> (`stat_summary()`).
12.8.2 Horizontal alignment
plot_modules_h() uses enrich_rank_by and enrich_top_n to control which
enrichment terms are displayed per module.
scaled
plot_modules_h(gene_module |> filter(Module != 0),
data_obj_merged, assay = 2, scale = TRUE,
ylabel = "Z-score",
enrich_list = go_list_wgcna$all,
enrich_category = "BP",
fontsize = 7,
heatmap_width = 8,
heatmap_height = 10,
save = "output", suffix = "scaled_BP_all",
width = 20, height = 15
)
#> Warning: Removed 8604 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 1131 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 900 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 597 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 579 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 207 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 201 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 174 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 8604 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 1131 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 900 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 597 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 579 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 207 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 201 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 174 rows containing non-finite outside the scale range
#> (`stat_summary()`).
# knitr::include_graphics(file.path("output",
# paste0("WGCNA_h_scaled_BP_all_", Sys.Date(), ".png")))The plot may have misalignment between annotations and the heatmap body in
interactive sessions. This is documented in ComplexHeatmap.
The saved figures can be used for correct display withknitr::include_graphics.
plot_modules_h(gene_module |> filter(Module != 0),
data_obj_merged, assay = 2, scale = TRUE,
ylabel = "Z-score",
enrich_list = go_list_wgcna$simplified,
enrich_category = "BP",
fontsize = 7,
heatmap_width = 8,
heatmap_height = 12,
save = "output", suffix = "scaled_BP_simplified",
width = 20, height = 15
)
#> Warning: Removed 8604 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 1131 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 900 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 597 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 579 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 207 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 201 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 174 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 8604 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 1131 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 900 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 597 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 579 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 207 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 201 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 174 rows containing non-finite outside the scale range
#> (`stat_summary()`).
# knitr::include_graphics(file.path("output",
# paste0("WGCNA_h_scaled_BP_simplified_", Sys.Date(), ".png")))enrich_category accepts a vector for multi-column enrichment display. For example,
the below chunk shows BP and CC side by side in a single heatmap:
plot_modules_h(gene_module |> filter(Module != 0),
data_obj_merged, assay = 2, scale = TRUE,
ylabel = "Z-score",
enrich_list = go_list_wgcna$all,
enrich_category = c("BP", "CC"),
enrich_top_n = 3, enrich_rank_by = "p.adjust",
fontsize = 7,
heatmap_width = 8,
heatmap_height = 12,
save = "output", suffix = "scaled_BP_CC",
width = 25, height = 15
)
#> Warning: Removed 8604 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 1131 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 900 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 597 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 579 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 207 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 201 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 174 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 8604 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 1131 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 900 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 597 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 579 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 207 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 201 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 174 rows containing non-finite outside the scale range
#> (`stat_summary()`).
12.9 Hub genes
12.9.1 WGCNA built-in: one hub per module
WGCNA::chooseTopHubInEachModule() identifies the single most connected gene in each module by evaluating intramodular connectivity against all other genes.
hub_genes_wgcna <- WGCNA::chooseTopHubInEachModule(net$input_data, net$colors,
omitColors = 0,
power = picked_power,
type = "signed"
)
hub_genes_wgcna#> 1 2 3 4 5 6 7 8
#> "Smc3" "Samd9l" "Nolc1" "Tecr" "Itgb8" "Sectm1a" "Clec4e" "Saa1"
# Show gene descriptions
clusterProfiler::bitr(hub_genes_wgcna, fromType = "SYMBOL",
toType = c("GENENAME"), OrgDb = org.Mm.eg.db)
#> 'select()' returned 1:1 mapping between keys and columns12.9.2 Top N hubs ranked by kME
extract_hubs() extracts the top N hub features per module ranked by
kME. Returns a flat character vector of feature names.
Top1 can be same or different from WGCNA::chooseTopHubInEachModule()
#> [1] "Itgb8" "Clec4e" "Sectm1a" "Nolc1" "Smc3" "Samd9l" "Cd38"
#> [8] "Tecr"
hubs_top3 <- extract_hubs(net, top_n = 3, exclude_grey = TRUE)
head(hubs_top3, 20)#> [1] "Itgb8" "Bcl2a1a" "Ptges" "Clec4e" "Mllt6" "Gpr84" "Sectm1a"
#> [8] "Gm4841" "Klrk1" "Nolc1" "Grwd1" "Rock2" "Smc3" "Rasgrp3"
#> [15] "Pank1" "Samd9l" "Nmi" "Ifit2" "Cd38" "Saa1"
clusterProfiler::bitr(hubs_top3, fromType = "SYMBOL",
toType = c("GENENAME"), OrgDb = org.Mm.eg.db)
#> 'select()' returned 1:1 mapping between keys and columns12.9.3 Label hubs on the module heatmap
plot_modules_h() accepts mark_features to label specific genes on the left
or right side of the heatmap. Passing the extracted hub genes highlights
them within each module.
plot_modules_h(gene_module |> filter(Module != 0),
data_obj_merged, assay = 2, scale = TRUE,
ylabel = "Z-score",
enrich_list = go_list_wgcna$simplified,
enrich_category = c("BP"),
mark_features = hub_genes_wgcna,
enrich_top_n = 3, enrich_rank_by = "p.adjust",
fontsize = 7,
heatmap_width = 8,
heatmap_height = 12,
save = "output", suffix = "hub_BP",
width = 25, height = 15
)
#> Warning: Removed 8604 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 1131 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 900 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 597 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 579 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 207 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 201 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 174 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 8604 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 1131 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 900 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 597 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 579 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 207 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 201 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 174 rows containing non-finite outside the scale range
#> (`stat_summary()`).
plot_modules_h(gene_module |> filter(Module != 0),
data_obj_merged, assay = 2, scale = TRUE,
ylabel = "Z-score",
enrich_list = go_list_wgcna$simplified,
enrich_category = c("BP", "Hub features"),
mark_features = hubs_top3,
enrich_top_n = 3, enrich_rank_by = "p.adjust",
fontsize = 7,
heatmap_width = 8,
heatmap_height = 12,
save = "output", suffix = "hub_BP_hubfeat",
width = 25, height = 15
)
#> Warning: Removed 8604 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 1131 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 900 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 597 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 579 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 207 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 201 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 174 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 8604 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 1131 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 900 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 597 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 579 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 207 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 201 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 174 rows containing non-finite outside the scale range
#> (`stat_summary()`).
13 Additional enrichment analyses
13.1 enrichR: drug signatures, pathways, TF targets
enrichR_list() wraps enrichR::enrichr() to query databases such as DSigDB (drug signatures), ChEA (TF targets), KEGG, and more. Results are compatible with plot_modules_h().
# Drug and pathway enrichment
enrichr_out <- enrichR_list(
gene_list = gene_module |> filter(Module != 0),
databases = c("KEGG_2019_Mouse", "WikiPathways_2024_Mouse"),
universe = background_wgcna
)
#> Welcome to enrichR
#> Checking connections ...
#> Enrichr ... Connection is Live!
#> FlyEnrichr ... Connection is Live!
#> WormEnrichr ... Connection is Live!
#> YeastEnrichr ... Connection is Live!
#> FishEnrichr ... Connection is Live!
#> OxEnrichr ... Connection is Live!
#> Connection changed to https://maayanlab.cloud/Enrichr/
#> Connection is Live!#> Uploading data to Speedrichr...
#> - Your gene set... Done.
#> - Your background... Done.
#> Getting enrichment results...
#> - KEGG_2019_Mouse... Done.
#> - WikiPathways_2024_Mouse... Done.
#> Parsing results... Done.
#> Uploading data to Speedrichr...
#> - Your gene set... Done.
#> - Your background... Done.
#> Getting enrichment results...
#> - KEGG_2019_Mouse... Done.
#> - WikiPathways_2024_Mouse... Done.
#> Parsing results... Done.
#> Uploading data to Speedrichr...
#> - Your gene set... Done.
#> - Your background... Done.
#> Getting enrichment results...
#> - KEGG_2019_Mouse... Done.
#> - WikiPathways_2024_Mouse... Done.
#> Parsing results... Done.
#> Uploading data to Speedrichr...
#> - Your gene set... Done.
#> - Your background... Done.
#> Getting enrichment results...
#> - KEGG_2019_Mouse... Done.
#> - WikiPathways_2024_Mouse... Done.
#> Parsing results... Done.
#> Uploading data to Speedrichr...
#> - Your gene set... Done.
#> - Your background... Done.
#> Getting enrichment results...
#> - KEGG_2019_Mouse... Done.
#> - WikiPathways_2024_Mouse... Done.
#> Parsing results... Done.
#> Uploading data to Speedrichr...
#> - Your gene set... Done.
#> - Your background... Done.
#> Getting enrichment results...
#> - KEGG_2019_Mouse... Done.
#> - WikiPathways_2024_Mouse... Done.
#> Parsing results... Done.
#> Uploading data to Speedrichr...
#> - Your gene set... Done.
#> - Your background... Done.
#> Getting enrichment results...
#> - KEGG_2019_Mouse... Done.
#> - WikiPathways_2024_Mouse... Done.
#> Parsing results... Done.
#> Uploading data to Speedrichr...
#> - Your gene set... Done.
#> - Your background... Done.
#> Getting enrichment results...
#> - KEGG_2019_Mouse... Done.
#> - WikiPathways_2024_Mouse... Done.
#> Parsing results... Done.
plot_modules_h(gene_module |> filter(Module != 0),
data_obj_merged, assay = 2, scale = TRUE,
ylabel = "Z-score",
enrich_list = enrichr_out,
enrich_category = "KEGG_2019_Mouse",
fontsize = 7,
heatmap_width = 8,
heatmap_height = 12,
save = "output", suffix = "enrichR_KEGG",
width = 20, height = 15
)
#> Warning: Removed 8604 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 1131 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 900 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 597 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 579 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 207 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 201 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 174 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 8604 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 1131 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 900 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 597 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 579 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 207 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 201 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 174 rows containing non-finite outside the scale range
#> (`stat_summary()`).
13.2 MSigDB enrichment
enrich_msigdb() performs hypergeometric enrichment against the Molecular
Signatures Database (MSigDB) via the msigdbr package. Set db_species to
the MSigDB species code ("MM" for mouse, "HS" for human) matching the gene
set database, and species to the query gene species. Ortholog mapping is
supported when the two differ.
# Hallmark gene sets (category "MH")
msigdb_h_res <- enrich_msigdb(
gene_list = gene_module |> filter(Module != 0),
category = "MH",
species = "Mus musculus", db_species = "MM",
universe = background_wgcna
)
#> Processing gene list: 1
#> Processing gene list: 2
#> Processing gene list: 3
#> Processing gene list: 4
#> Processing gene list: 5
#> Processing gene list: 6
#> Processing gene list: 7
#> Processing gene list: 8
plot_modules_h(gene_module |> filter(Module != 0),
data_obj_merged, assay = 2, scale = TRUE,
ylabel = "Z-score",
enrich_list = msigdb_h_res,
enrich_category = "MH",
fontsize = 7,
heatmap_width = 8,
heatmap_height = 12,
save = "output", suffix = "MSigDB_H",
width = 20, height = 15
)
#> Warning: Removed 8604 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 1131 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 900 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 597 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 579 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 207 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 201 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 174 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 8604 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 1131 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 900 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 597 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 579 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 207 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 201 rows containing non-finite outside the scale range
#> (`stat_summary()`).
#> Warning: Removed 174 rows containing non-finite outside the scale range
#> (`stat_summary()`).
14 Alternative: one-step data preparation
prepare_tide() wraps the steps above (create_input(), normalise_to_start(),
split_groups(), merge_replicates(), calc_feature_property(), variance
decomposition, residual filtering, merge_groups()) into a single
function call. It returns a list with original and filtered SEs, merged SEs,
variance decomposition results, feature properties, and empty slots for DE,
enrichment, and WGCNA results that downstream functions can fill.
tide_out <- prepare_tide(
data = geo_data_tb_new,
sample_ann = geo_sample_info_new,
filter_ratio = 0.5,
keep = "below_quantile",
residual_threshold = 0.75,
n_cores = 16,
interaction = TRUE
)
# tide_out$se # original SE (with normalised assay)
# tide_out$se_filtered # filtered SE (by group-specific features and residual variance)
# tide_out$merged_se # merged (mean of replicates)
# tide_out$merged_se_filtered # filtered + merged
# tide_out$variance # variance decomposition
# tide_out$feature_property # feature property table
# tide_out$filter_summary # per-step filtering statistics15 Session info
#> R version 4.6.0 (2026-04-24 ucrt)
#> Platform: x86_64-w64-mingw32/x64
#> Running under: Windows 11 x64 (build 26200)
#>
#> 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.81.28
#> [5] org.Mm.eg.db_3.23.0 AnnotationDbi_1.75.0
#> [7] SummarizedExperiment_1.43.0 Biobase_2.73.1
#> [9] GenomicRanges_1.65.1 Seqinfo_1.3.0
#> [11] IRanges_2.47.1 S4Vectors_0.51.2
#> [13] BiocGenerics_0.59.10 generics_0.1.4
#> [15] MatrixGenerics_1.25.0 matrixStats_1.5.0
#> [17] TiDEomics_0.99.4 BiocStyle_2.41.0
#>
#> loaded via a namespace (and not attached):
#> [1] segmented_2.2-1 fs_2.1.0
#> [3] bitops_1.0-9 enrichplot_1.33.0
#> [5] httr_1.4.8 RColorBrewer_1.1-3
#> [7] doParallel_1.0.17 ggsci_5.1.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.3 gridExtra_2.3.1
#> [17] preprocessCore_1.75.0 WGCNA_1.74
#> [19] cli_3.6.6 enrichR_3.4
#> [21] scatterpie_0.2.6 labeling_0.4.3
#> [23] sass_0.4.10 S7_0.2.2
#> [25] readr_2.2.0 askpass_1.2.1
#> [27] ggridges_0.5.7 pbapply_1.7-4
#> [29] systemfonts_1.3.2 yulab.utils_0.2.4
#> [31] gson_0.2.0 foreign_0.8-91
#> [33] R.utils_2.13.0 DOSE_4.7.2
#> [35] rentrez_1.2.4 WriteXLS_6.8.0
#> [37] limma_3.69.2 rstudioapi_0.19.0
#> [39] impute_1.87.0 RSQLite_3.53.3
#> [41] gridGraphics_0.5-1 shape_1.4.6.1
#> [43] crosstalk_1.2.2 gtools_3.9.5
#> [45] car_3.1-5 GO.db_3.23.1
#> [47] Matrix_1.7-5 abind_1.4-8
#> [49] R.methodsS3_1.8.2 PCAtools_2.25.0
#> [51] lifecycle_1.0.5 yaml_2.3.12
#> [53] carData_3.0-6 gplots_3.3.0
#> [55] qvalue_2.45.0 SparseArray_1.13.2
#> [57] grid_4.6.0 blob_1.3.0
#> [59] promises_1.5.0 dqrng_0.4.1
#> [61] crayon_1.5.3 ggtangle_0.1.2
#> [63] lattice_0.22-9 msigdbr_26.1.0
#> [65] beachmat_2.29.0 cowplot_1.2.0
#> [67] KEGGREST_1.53.1 magick_2.9.1
#> [69] pillar_1.11.1 knitr_1.51
#> [71] ComplexHeatmap_2.29.0 rjson_0.2.23
#> [73] boot_1.3-32 codetools_0.2-20
#> [75] glue_1.8.1 ggiraph_0.9.6
#> [77] ggfun_0.2.1 fontLiberation_0.1.0
#> [79] data.table_1.18.4 vctrs_0.7.3
#> [81] png_0.1-9 treeio_1.37.0
#> [83] Rdpack_2.6.6 gtable_0.3.6
#> [85] assertthat_0.2.1 cachem_1.1.0
#> [87] xfun_0.60 rbibutils_2.4.1
#> [89] S4Arrays_1.13.0 mime_0.13
#> [91] reformulas_0.4.4 survival_3.8-9
#> [93] aisdk_1.4.12 iterators_1.0.14
#> [95] statmod_1.5.2 nlme_3.1-170
#> [97] ggtree_4.3.0 bit64_4.8.2
#> [99] fontquiver_0.2.1 bslib_0.11.0
#> [101] irlba_2.3.7 KernSmooth_2.23-26
#> [103] otel_0.2.0 rpart_4.1.27
#> [105] colorspace_2.1-3 DBI_1.3.0
#> [107] Hmisc_5.2-6 nnet_7.3-20
#> [109] tidyselect_1.2.1 processx_3.9.0
#> [111] curl_7.1.0 bit_4.6.0
#> [113] compiler_4.6.0 httr2_1.3.0
#> [115] htmlTable_2.5.0 xml2_1.6.0
#> [117] plotly_4.12.0 randtests_1.0.2
#> [119] fontBitstreamVera_0.1.1 DelayedArray_0.39.3
#> [121] bookdown_0.47 checkmate_2.3.4
#> [123] scales_1.4.0 caTools_1.18.4
#> [125] callr_3.8.0 rappdirs_0.3.4
#> [127] digest_0.6.39 minqa_1.2.8
#> [129] rmarkdown_2.31 XVector_0.53.0
#> [131] htmltools_0.5.9 pkgconfig_2.0.3
#> [133] base64enc_0.1-6 umap_0.2.10.0
#> [135] lme4_2.0-6 sparseMatrixStats_1.25.0
#> [137] fastmap_1.2.0 rlang_1.3.0
#> [139] GlobalOptions_0.1.4 htmlwidgets_1.6.4
#> [141] shiny_1.14.0 DelayedMatrixStats_1.35.0
#> [143] ggh4x_0.3.1 farver_2.1.2
#> [145] jquerylib_0.1.4 jsonlite_2.0.0
#> [147] BiocParallel_1.47.0 R.oo_1.27.1
#> [149] GOSemSim_2.39.2 BiocSingular_1.29.0
#> [151] magrittr_2.0.5 Formula_1.2-5
#> [153] ggplotify_0.1.3 patchwork_1.3.2
#> [155] Rcpp_1.1.2 babelgene_22.9
#> [157] reticulate_1.46.0 ape_5.8-1
#> [159] ggnewscale_0.5.2 gdtools_0.5.1
#> [161] stringi_1.8.7 MASS_7.3-66
#> [163] plyr_1.8.9 shinyFiles_0.9.3
#> [165] parallel_4.6.0 ggrepel_0.9.8
#> [167] Biostrings_2.81.5 splines_4.6.0
#> [169] hms_1.1.4 circlize_0.4.18
#> [171] igraph_2.3.3 ggpubr_1.0.0
#> [173] fastcluster_1.3.0 ggsignif_0.6.4
#> [175] enrichit_0.2.0 reshape2_1.4.5
#> [177] ScaledMatrix_1.21.0 XML_3.99-0.23
#> [179] evaluate_1.0.5 BiocManager_1.30.27
#> [181] tzdb_0.5.0 nloptr_2.2.1
#> [183] foreach_1.5.2 tweenr_2.0.3
#> [185] httpuv_1.6.17 openssl_2.4.2
#> [187] tidyr_1.3.2 purrr_1.2.2
#> [189] polyclip_1.10-7 clue_0.3-68
#> [191] Trendy_1.35.0 ggforce_0.5.0
#> [193] rsvd_1.0.5 broom_1.0.13
#> [195] xtable_1.8-8 RSpectra_0.16-2
#> [197] tidytree_0.4.8 tidydr_0.0.6
#> [199] rstatix_1.0.0 later_1.4.8
#> [201] viridisLite_0.4.3 tibble_3.3.1
#> [203] clusterProfiler_4.21.0 aplot_0.3.1
#> [205] memoise_2.0.1 cluster_2.1.8.2References
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.