---
title: "TiDEomics Tutorial"
author:
  - name: Tianen He
    affiliation:
    - University of Oxford
    email: tianen.he@ndm.ox.ac.uk
output:
  BiocStyle::html_document:
    self_contained: true
    toc: true
    toc_float:
      collapsed: true
      smooth_scroll: false
    toc_depth: 3
    number_sections: true
    code_folding: show
    highlight: tango
    df_print: paged
date: "`r BiocStyle::doc_date()`"
package: "`r BiocStyle::pkg_ver('TiDEomics')`"
vignette: >
  %\VignetteIndexEntry{TiDEomics Tutorial}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
bibliography: references.bib
nocite: '@*'
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>",
    crop = NULL,
    # Related to https://stat.ethz.ch/pipermail/bioc-devel/
    # 2020-April/016656.html
    dpi = 96,
    fig.width = 6,
    fig.height = 4,
    fig.keep = "all",
    max.height = "300px"
)

local({
    hook_output <- knitr::knit_hooks$get("output")
    knitr::knit_hooks$set(output = function(x, options) {
        if (!is.null(options$max.height)) {
            options$attr.output <- c(
                options$attr.output,
                sprintf('style="max-height: %s;"', options$max.height)
            )
        }
        hook_output(x, options)
    })
})
```

## Introduction

[**TiDEomics**](https://github.com/hte123/TiDEomics) (Time-course Differential
Expression analysis of omics data) is an R package for
**multi-group time-course omics analysis** built on `SummarizedExperiment`,
especially for comparison of more than two conditions across time.

`TiDEomics` provides:

- An integrated workflow from data processing to functional interpretation, on
    different omics data types in feature x sample matrix format, with
    handling of missing values.

- An approach to separate **time-dominant**, **group-dominant**, and
    **group-specific temporal** effects through the combination of pairwise
    differential expression, variance decomposition, and co-expression module
    analysis (WGCNA).

- A filtering strategy based on **residual variance** after mixed-model
    variance decomposition, prioritising features with structured
    differential expression over unexplained variation. In comparison,
    traditional filtering by total variance or coefficient of variation is less
    appropriate when biological variation across groups and time course is
    expected.

Existing packages on time-course omics analysis are available, typically
designed for transcriptomics data. Some packages are limited to two-group
comparisons (e.g. `splineTimeR`) or expect raw sequencing counts (e.g.
`ImpulseDE2`, `TCseq`). `maSigPro` supports multi-group time-course DE analysis
via polynomial regression, complementary to `TiDEomics` with fewer
assumptions about temporal trends. `TiDEomics` also integrates QC, filtering
and downstream functional enrichment within the same workflow.

TiDEomics incorporates `Trendy` [@Trendy] which operates on single time
series to estimate breakpoints per group per feature,
and derived variance decomposition from `PALMO` [@PALMO] and introduced
residual variance as a filtering criterion.

This tutorial demonstrates a basic workflow, explains key parameters, and gives
practical tips.

Key features covered:

- *Preparing input*: Create a `SummarizedExperiment` object with the data
    matrix and sample annotation, with checks for correct formatting.

- *Quality control*: Visualise value distributions and missingness patterns.

- *Normalisation, grouping and merging replicates*: Normalise each feature to
    the starting time point when focused on changes from baseline. Transform
    the data for analyses that require single groups or single time series.

- *Sample relationships visualisation (correlation matrix, PCA, UMAP)*: Explore
    global structure and major sources of variance, check for batch effects,
    and identify features driving PCs.

- *Pairwise differential expression (by time or group)*: Identify DE features
    between pairs of time points within each group, and between pairs of groups
    at each time point.

- *Classification by feature property (randomness and overall fold change)*:
    Classify features into different categories (e.g., non-random vs. random,
    differentially expressed vs. stable) based on properties including trend
    significance and maximum fold change.

- *Segmented regression with Trendy*: Estimate breakpoints for each feature in
    each group and summarise dynamic patterns.

- *Variance decomposition modified from PALMO*: Quantify relative contributions
    from `Time`, `Group`, and `Residual`, distinguish time or group-dependent
    DE features and "noisy" features.

- *Module identification with WGCNA*: Identify co-expression modules with
    different temporal and group-specific patterns.

- *Functional enrichment (gene ontology, drugs, etc.)*: Interpret biological
    meaning of identified DE features and modules.

## Installation

```{r install, eval = FALSE}
if (!requireNamespace("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager")
}

BiocManager::install("TiDEomics")
```

```{r load-packages}
library(TiDEomics)
library(SummarizedExperiment)
library(org.Mm.eg.db)
```

## Example input

TiDEomics expects:

- `data`: log2-transformed data.frame
    - rows: features (e.g. genes, proteins)
    - columns: First column should be feature identifiers.
        Other columns are numeric measurements of features for each sample.
    - Missing values (NA) are allowed.
- `sample_ann`: data.frame with columns:
    - `Sample`: sample names, match column names of `data`
    - `Group`: experimental group
    - `Time`: numeric time point
    - `Replicate`: replicate ID (optional, auto-generated if absent)
    - `Batch`: batch ID (optional)
    - `Subject`: biological subject ID for repeated-measures designs
        (optional, set via `subject_col`)

Example: subset of
[GSE263759](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE263759)
data set published in @Traxler2025.

- Ensembl IDs were mapped to symbols, genes with all zero counts were excluded.
- Use time 0 untreated samples for other groups' time 0.
- Sample 500 random genes.
- Normalised to `log2(CPM + 1)`

```{r load-example}
data("tutorial_sample_info", package = "TiDEomics")
data("tutorial_data", package = "TiDEomics")
```

### Preprocessing for different data types

TiDEomics expects log-transformed, normalised data as input. The core workflow
is compatible with any quantitative omics data in a feature x sample matrix.
Enrichment analysis (`enrichGO_list`, `enrichGO_rank`, `enrichR_list`,
`enrich_msigdb`) requires gene identifiers and is not applicable to non-gene
features (e.g., metabolites, lipids).

Batch correction should be applied before input if batch effects are present.

For the most common scenarios in transcriptomics and proteomics:

- RNA-seq data should be normalised for sequencing depth and composition bias
    before TiDEomics, e.g. with TMM. The tutorial uses `log2(CPM + 1)` for
    simplicity. For pairwise differential expression, `trend = TRUE` can be used
    to account for the mean-variance relationship. See below section on
    [pairwise differential expression](#pairwise-de) for details.

- MS-based proteomics data processing software may include normalisation when
    producing protein x sample matrix. Depending on software and settings,
    manual normalisation and log transformation may be required, e.g. median
    normalisation. Missing values are accepted in most TiDEomics functions,
    and a basic imputation function `impute_groups()` is provided.

## Workflow

We present the workflow as below sections, explaining the usage and parameters.

### Data preparation

```{r create-input}
data_obj <- create_input(
    data = tutorial_data,
    sample_ann = tutorial_sample_info)

data_obj
```

The `SummarizedExperiment` stores data in named **assays**: `"orig"`
(assay 1) for the original input and `"norm"` (assay 2) for time-0-normalised
values (added by `normalise_to_start()`). Downstream functions accept either the
name or the index.

```{r view-data}
assays(data_obj)[["orig"]] |> as.data.frame() |> utils::head()
```

```{r coldata}
colData(data_obj) # sample annotation
```

### Optional: custom color palette

By default, TiDEomics generates a default color palette (`scales::pal_hue()`)
based on the number of groups. A custom palette for each group can be set with
`set_custom_palette()`, which will be used in all subsequent plotting functions
where applicable.

```{r palette}
custom_palette <- c(
    "untreated" = "#1b9e77", "IFNbeta" = "#d95f02",
    "IFNgamma" = "#7570b3", "LPS" = "#e7298a"
)
set_custom_palette(custom_palette)
```

### Quality control

```{r abundance}
plot_distribution(data_obj, facet_by = "Group")
```

If distributions differ significantly, consider global normalisation (e.g.,
quantile, median) *before* using TiDEomics.

For data with missing values (e.g., proteomics), `plot_ID()` and
`plot_missing()` can be used to check missingness patterns.

### Normalisation to starting time point

`normalise_to_start()` subtracts the baseline value from each feature at
starting time point (referred to as time 0) or the first available time point
if the feature is missing at time 0.

Two modes are available for defining the baseline:

- `by_subject = FALSE` (default): group-level baseline, mean of all samples at
    starting time point in the group. This is used when no subject-level
    information is available or when between-subject baseline differences are
    of interest.

- `by_subject = TRUE`: subject-level baseline, mean of all samples at starting
    time point for each subject. This is used when subject-level information
    is available and the focus is on subject-specific changes from baseline.

```{r norm0}
data_obj <- normalise_to_start(data_obj)
```

Both original and time-0 normalised data are stored in the
`SummarizedExperiment` object and downstream functions allow the user to specify
which to use for each analysis (`assay = 1 or 'orig'` for original,
`assay = 2 or 'norm'` for time-0 normalised).

Use original data when absolute abundance
differences across groups are of interest (e.g., constitutively different
baseline levels). Use time-0 normalised data when relative changes from
baseline are of interest.

### Splitting groups and merging replicates

Use `split_groups()` to obtain group-specific `SummarizedExperiment` objects
and `merge_replicates()` to average replicates for analyses / visualisation
that require single time series for features, including
`calc_feature_property()`, `run_Trendy()`, `plot_modules_v()` and
`plot_modules_h()`.

```{r split}
data_obj_list <- split_groups(data_obj)
```

```{r merge-rep}
data_obj_merged_list <- merge_replicates(data_obj_list)
```

```{r merge-groups}
data_obj_merged <- merge_groups(data_obj_merged_list)
```

### Sample relationships (correlation, PCA, UMAP)

Correlation matrix can be plotted with `plot_cor_matrix()` to check sample
relationships and potential batch effects.

```{r cormat}
plot_cor_matrix(data_obj,
    method = "spearman",
    label_rep = TRUE, label_batch = TRUE,
    cellwidth = 2, cellheight = 2
)
```

High within-replicate correlation and clustering by `Group` or `Time` increase
confidence in data quality and downstream analysis. If batch effects are
present, consider batch correction before using TiDEomics.

There are multiple functions for PCA and UMAP:

- Use `plot_pca()` to visualise sample relationships in 2D. Output loadings
    can be used to identify features driving the principal components.

- Use `plot_pca_3D()` to visualise sample relationships in 3D.

- Use `PCAtools::eigencorplot()` to correlate PCs with Group and Time.

- Use `plot_umap()` for UMAP and set `seed` for reproducibility.

Features with missing values (NA) are automatically excluded from PCA and
UMAP. Consider imputation (e.g., group-wise minimum imputation with
`split_groups()`, `impute_groups()` and `merge_groups()`) to include features
with missing values in PCA and UMAP.

```{r pca}
PC <- plot_pca(data_obj,
    # pc1 = 1, pc2 = 2, # default to plot PC1 and PC2
    plot_screeplot = TRUE,
    plot_loadings = FALSE,
    plot_morepc = TRUE,
    circle = FALSE,
    morepc = 1:5
)
PC$p_list
```

```{r pca-3D}
plot_pca_3D(PC$pca, pcs = 1:3)
```

Note: the 3D plot may not display properly in some html, but should work in an
interactive R session.

```{r pca-eigencor}
PCAtools::eigencorplot(PC$pca,
    metavars = c("Group", "Time"),
    components = paste0("PC", 1:5),
    col = colorRampPalette(c("#3C5488FF", "white", "#E64B35FF"))(100),
    colCorval = "black"
)
```

```{r umap}
umap <- plot_umap(data_obj, seed = 1234)
umap$p_list
```

PCA and UMAP can also be performed separately on each group to explore
within-group structure and dynamics.

By default, when plotting PCA by group, circles and arrows are added to show
the trajectory of samples along time course. Set `circle = FALSE` or
`arrow = FALSE` to remove circles and arrows.

```{r pca-umap-group, fig.width = 8, fig.height = 2.5}
plot_pca_by_group(data_obj, circle = TRUE, arrow = TRUE, legend_pos = "top")
# also accepts a list: plot_pca_by_group(data_obj_list)

plot_umap_by_group(data_obj, seed = 1234, legend_pos = "top")
# also accepts a list: plot_umap_by_group(data_obj_list)
```

### Pairwise differential expression {#pairwise-de}

The two DE functions for pairwise comparison between time points and groups are
built based on the limma package [@limma].

`DE_between_time()`:

- Compare pairs of time points within each group.

- Output: Nested lists of DE statistics tables, including all features or
    only significant features.

- Use `plot_DE_between_time()` to summarise number of DE features per pair of
    time points per group.

- When a Subject column is provided via `subject_col` in `create_input()`,
    `DE_between_time()` automatically uses a paired design via
    `limma::duplicateCorrelation(block = Subject)`.

`DE_between_group()`:

- Compare pairs of groups at each time point.

- Output: Nested lists of DE statistics tables, including all features or
    only significant features.

- Use `plot_DE_between_group()` to summarise number of DE features per pair of
    groups with line plots.

Common parameters for both functions:

- `assay`: index of assay to use (1 = original, 2 = time0-normalised if
    available)

- `filter`: minimal number of non-NA replicates for the feature to be included
    in DE testing (e.g., `filter = 1` means requiring features to have at least
    1 non-NA value at both time points in the specific group, or in both groups
    at the specific time point).

- `trend`: passed to `limma::eBayes()`. Set to `TRUE` for RNA-seq
    count-derived data to model the mean-variance trend.
    `FALSE` (default) is usually appropriate for other types of data,
    where the trend is typically weak.
    The trend can be visualised with `limma::plotSA()` with elements of
    `$fit_list` of both functions' output.

- Default thresholds for DE are `p.adj < 0.05` and `|log2FC| > 1`.

- `plot_DE_between_time()` and `plot_DE_between_group()` allow re-filtering
    DE by setting new thresholds.

Please note that the example data only has two replicates per time point.
More replicates are recommended for robust pairwise DE analysis.

DE_between_time():

```{r de-between-time, fig.width = 8, fig.height = 2, warning = FALSE}
DE_between_time_out <- DE_between_time(data_obj, assay = 1,
    filter = 1, trend = FALSE)
```

```{r de-between-time-all}
DE_between_time_out$all_list$IFNbeta$`t2-t0` |> utils::head()
```

```{r de-between-time-de}
DE_between_time_out$de_list$IFNbeta$`t2-t0` |> utils::head()
```

```{r de-between-time-plot, fig.width = 7, fig.height = 2}
# Filtered with thresholds in `DE_between_time()`
plot_DE_between_time(DE_between_time_out,
    fontsize = 8, value = FALSE, nrow = 1, heatmap_width = 3
)

# Re-filtering with new thresholds
plot_DE_between_time(DE_between_time_out,
    fontsize = 8, value = FALSE, nrow = 1, heatmap_width = 3,
    adjP_thres = 0.01, logFC_thres = 1
)
```

DE_between_group():

```{r de-between-group}
DE_between_group_out <- DE_between_group(data_obj, assay = 2,
    filter = 1, trend = TRUE)
```

```{r de-between-group-plot, fig.width = 3, fig.height = 2, out.width = "70%"}
# Filtered with thresholds in `DE_between_group()`
plot_DE_between_group(DE_between_group_out)

# Re-filtering with new thresholds
plot_DE_between_group(DE_between_group_out, adjP_thres = 0.01, logFC_thres = 1)
```

```{r de-between-group-all}
DE_between_group_out$all_list$`IFNgamma-untreated`$`24` |> utils::head()
```

```{r de-between-group-de}
DE_between_group_out$de_list$`IFNgamma-untreated` |> utils::head()
```

Output of `DE_between_group()`and `DE_between_time()` can be plotted with
`plot_volcano()` to visualise DE features of selected group(s) and time
point(s).

```{r volcano}
plot_volcano(DE_between_group_out,
    group1 = "untreated", group2 = "IFNgamma", time = 24,
    logFC_thres = 0.5, adjP_thres = 0.05, label = TRUE)
```

### Feature properties and classification

`calc_feature_property()` computes per-feature properties:

- `P_trend`: calculated with `randtests::bartels.rank.test()`, testing for
    non-randomness of the time profile (i.e., whether the feature shows a
    significant overall trend across time points).

- `Max_FC`: maximum fold change across time points, calculated as the
    difference between the feature's maximum and minimum abundance across time
    points. This captures the overall magnitude of change across the time
    course, regardless of the specific time points at which changes occur.

- `Max_FC_time`: the difference between time points with maximum and minimum
    abundance, providing information about the direction (up or down) and
    duration of changing.

- `Exp_threshold`: user-set value threshold for a feature to be considered
    expressed in a sample. This can be used to filter features based on
    expression level, e.g. setting `threshold = 0` means that only 
    values > 0 (log2(CPM + 1) > 0, i.e. CPM > 0 in the tutorial dataset) are 
    considered expressed. The default is NULL, which means all non-NA
    values are considered expressed.

- `T_total`, `T_exp` and `Exp_ratio`: number of total time points, number of
    expressed time points with values > `Exp_threshold` (or non-NA values if
    `threshold = NULL`), and proportion of expressed time points (`T_exp` /
    `T_total`).

- `Rho_time`: Spearman correlation with time (positive = up, negative = down).
    Requires at least 3 expressed time points.

- `Peak_ratio`: proportion of time points that are local maxima (0 = monotonic,
    higher = more oscillatory).

- `Log2_CV`: `log2(1 + SD / |mean|)`, coefficient of variation for relative
    temporal variability.

- `AUC`: area under the time-0-normalised curve from assay 2. Positive =
    net increase from baseline, negative = net decrease. `NA` if assay 2 is not
    available (run `normalise_to_start()` first).

The property values can be used to classify features into different categories
(e.g., non-random vs. random, differentially expressed vs. stable) for
downstream analysis and interpretation. For example, features with
`P_trend < 0.05` and `Max_FC >= 1` can be candidates for non-randomly
overall-changing DE features.

Note:

- When replicates are present, the input should be the mean of replicates at
    each time point, output of `merge_replicates()`.

- The function is set to require at least 3 values > threshold to
    calculate `P_trend`.

```{r property}
data_obj_merged_list <- calc_feature_property(data_obj_merged_list,
    threshold = 0)
property_tb <- summarise_feature_property(data_obj_merged_list)

utils::head(property_tb)
```

Features expressed in only one group or a subset of groups can be extracted
with `group_specific_features()` based on the output table of
`summarise_feature_property()`.

```{r group-specific}
group_specific_features(property_tb, groups = c("untreated"),
    genename = FALSE, GO = FALSE
)
```

### Segmented regression analysis

`run_Trendy` fits segmented linear models to the time course for each feature
to estimate breakpoints with Trendy package [@Trendy], see
[Trendy](https://bioconductor.org/packages/release/bioc/html/Trendy.html)
for details.

**Notes**:

- The function requires complete input (no NA). When data contains missing
    values (e.g. proteomics), use `impute_groups()` to impute (default: group
    minimum; pass `fun = function(x) min(x) / 2` for half-minimum,
    `fun = median` for median imputation, etc.), or restrict to features with
    complete data.

- If `feature` is not specified, the function will run on all features
    filtered by expression / missing rate, which requires running
    `calc_feature_property()` before `impute_groups()`.

- Number of time points needed = [# segments] x [minimum number of samples in
    a segment]. For example, when # segments = 2 (maxK = 1) and
    minNumInSeg = 2, at least 4 time points are needed.

```{r impute}
data_obj_merged_imp_list <- impute_groups(data_obj_merged_list)
```

A subset of features is used for demonstration as `run_Trendy()` can be
time-consuming.

```{r trendy-run}
set.seed(1234)
random_features <- sample(rownames(data_obj_merged_imp_list[[1]]), 50)

example_res_list <- run_Trendy(data_obj_merged_imp_list,
    feature = random_features,
    minExp = 0.5,
    maxK = 1,
    minNumInSeg = 2, meanCut = 0, NCores = 2
)
```

Use `plot_segments()`to plot the fitted segments and breakpoints for
selected features.

```{r trendy-plot, fig.height = 2.5, fig.width = 8}
plot_segments(data_obj_merged_imp_list,
    example_res_list,
    feature = c("Slc25a51", "Aunip"), # example features
    ylab = "Log2(CPM + 1)"
)
```

Plot the distribution of breakpoints across time points with
`plot_breakpoints()` in each group, and summarise the temporal patterns
(combination of trends of each segment) of features with `summarise_Trendy()`
and `extract_segment_trends()`.

```{r trendy-summary, fig.height = 2, fig.width = 4}
plot_breakpoints(example_res_list)

trendy_summary <- summarise_Trendy(example_res_list)
trendy_summary |> utils::head()
```

```{r trendy-list}
trendy_list <- extract_segment_trends(trendy_summary)
trendy_list$IFNbeta
```

### Variance decomposition

Variance of each feature is decomposed by linear mixed models (LMM) into
contributions from `Group`, `Subject` (when present), `Time`, `Residual`, and
optionally `Group:Time` (or `Subject:Time`) interaction, which captures
group-specific temporal patterns. This helps to identify time-dependent
features, group-dependent features, and "noisy" features with high residual
variance.

**Interpreting the components:**

- High `Time`: time-dependent features, dynamic along time course consistently
    across groups. Good candidates for `run_Trendy()` or other time-focused
    analysis.

- High `Group`: group-dependent features, stable along time but differentially
    expressed between groups. May be baseline biological markers.

- High `Subject` (when present): features with baseline differences between
    individuals, biologically meaningful in patient studies.

- High `Residual`: unexplained variation, consider excluding for
    downstream analysis.

The function `decomp_variance()` is inspired by `PALMO::lmeVariance()`[@PALMO].
Group and Time are always included (auto-skipped if only 1 group or time point
present). `Subject` is included when present in colData. Set
`interaction = TRUE` to add `Group:Time` (or `Subject:Time` with Subject) as a
variance component capturing group-specific temporal patterns. Sufficient
replicates per combination are required for stable estimates.

```{r variance-decomposition, fig.height = 3}
# filter genes for variance decomposition:
# at least 50% values > 0 in at least 2 groups
decomp_filter_genes <- group_specific_features(property_tb,
    filter_ratio = 0.5,
    group_pct = 2 / 4,
    GO = FALSE, genename = FALSE
)$features

var_decomp <- decomp_variance(data_obj,
    features = decomp_filter_genes,
    assay = 1, core = 2
)

plot_variance(var_decomp, rank = "Time", top_n = 20)
plot_variance(var_decomp, rank = "Group", top_n = 20)
```

### WGCNA (co-expression modules)

Detailed tutorials of weighted gene correlation network analysis
([WGCNA](https://cran.r-project.org/web/packages/WGCNA/index.html))
[@WGCNA_1; @WGCNA_2] can be found
[online](https://github.com/edo98811/WGCNA_official_documentation/).
The following is a brief demonstration of how to prepare input and run WGCNA
with TiDEomics functions.

**Filtering strategies**:

- Features with too many missing values will be automatically removed with
    `WGCNA::goodGenes()` included in the `prepare_WGCNA()` function, and can
    also be pre-filtered with the calculated feature property (output of
    `calc_feature_property()` and `summarise_feature_property()`).

- Residual variance can be used to exclude noisy features.

- It is not recommended to filter by differential expression before WGCNA, see
[WGCNA FAQ](https://edo98811.github.io/WGCNA_official_documentation/faq.html).

```{r wgcna-prepare}
# Example filtering by residual variance < Q3
var_res_q3 <- stats::quantile(var_decomp$Residual, 0.75, na.rm = TRUE)
filter_wgcna <- var_decomp |> dplyr::filter(Residual < var_res_q3) |>
    dplyr::pull(Feature)

data_obj_wgcna <- data_obj[filter_wgcna, ]
```

**Prepare data format and choose power**:

`prepare_WGCNA()` prepares the input for WGCNA and helps to choose the
soft-thresholding power. Check the scale-free topology fit indices to confirm
that the chosen power is appropriate.

Reasonable powers are less than 15 for unsigned or signed hybrid networks, and
less than 30 for signed networks, to reach scale-free topology fit index > 0.8,
and mean connectivity (mean.k) not too high (in the hundreds or above). See
[WGCNA FAQ](https://edo98811.github.io/WGCNA_official_documentation/faq.html)
for details.

```{r choose-power}
wgcna_input <- prepare_WGCNA(data_obj_wgcna, assay = 2,
    powers = seq(1, 20),
    networkType = "signed", RsquaredCut = 0.8
)
wgcna_input$plot
picked_power <- wgcna_input$powerEstimate
picked_power
```

**Run WGCNA**:

Use `run_WGCNA()` with the output of `prepare_WGCNA()`, and visualise the
resulting modules with `plot_WGCNA()`.

`run_WGCNA()` is a wrapper of `WGCNA::blockwiseModules()`, which runs the WGCNA
analysis and returns the module assignment for each feature. Parameters for
`WGCNA::blockwiseModules()` can be set with `run_WGCNA()`, e.g., `corType`
for correlation method.

```{r run-wgcna}
net <- run_WGCNA(wgcna_input,
    power = picked_power,
    # corType = "pearson", # other option is "bicor"
    numericLabels = TRUE
)

plot_WGCNA(net, fontsize = 8)
```

**Extract modules**: show module sizes

```{r extract-wgcna-modules}
gene_module <- WGCNA_module(net, exclude_grey = TRUE)

gene_module |>
    dplyr::group_by(Module) |>
    dplyr::summarise(n = dplyr::n())
```

Module metrics (size, mean kME, etc.) can be summarised with
`summarise_module_metrics()`:

```{r module-metrics}
summarise_module_metrics(net)
```

**Plot module profiles**:

```{r plot-wgcna-modules, fig.width = 6, fig.height = 6}
plot_modules_v(gene_module,
    data_obj_merged, assay = 2, scale = TRUE,
    ylabel = "Z-score of log2(CPM + 1)",
    height_ratio = 2
)
```

### Functional enrichment

TiDEomics wrap `clusterProfiler::enrichGO()` and
`clusterProfiler::gseGO()`[@clusterProfiler_1; @clusterProfiler_2;
@clusterProfiler_3; @clusterProfiler_4]
to run GO enrichment efficiently on gene sets and ranked gene list.

For ranked gene list (e.g., ranked by Time or Group contribution from variance
decomposition), use `enrichGO_rank()` and plot with `enrichplot::gseaplot2()`.

```{r go_rank, fig.width=5, fig.height=3}
gse_group <- enrichGO_rank(var_decomp,
    gene_rank_by = "Group",
    OrgDb = org.Mm.eg.db,
    keyType = "SYMBOL", category = "BP",
    go_rank_by = "p.adjust")

enrichplot::gseaplot2(gse_group[["BP"]], geneSetID = 1:3, base_size = 8)
```

For multiple defined gene sets (e.g. WGCNA modules, DE genes by feature
classification), use `enrichGO_list()` and plot with `plot_GO()`.

Optionally, use `simplify = TRUE` with `enrichGO_list()` to simplify the GO
results by removing redundant terms (default: FALSE).

```{r go-example, fig.width = 6}
background_wgcna <- colnames(wgcna_input$data)

go_list <- enrichGO_list(
    gene_list = gene_module, OrgDb = org.Mm.eg.db,
    universe = background_wgcna,
    pvalueCutoff = 0.9, # get more results for demonstration
    qvalueCutoff = 0.9,
    simplify = FALSE,
    keyType = "SYMBOL"
)

plot_GO(go_list$all, plot_dotplot = TRUE,
    plot_emapplot = FALSE,
    plot_cnetplot = FALSE,
    showCategory_dotplot = 3)
```

Enrichment can also be performed against the Molecular Signatures Database
([MSigDB](https://www.gsea-msigdb.org/gsea/msigdb)) via `enrich_msigdb()`,
which queries the gene sets with `msigdbr::msigdbr()`
and runs a hypergeometric test across
user-selected collections (Hallmark, GO, etc.) or specific gene sets.

```{r enrich-msigdb}
# Mouse Hallmark gene sets
hallmark_msigdb <- enrich_msigdb(gene_module, universe = background_wgcna,
    minGSSize = 5, category = "MH", species = "Mus musculus", db_species = "MM")

# Human C2 (chemical and genetic perturbations)
perturb_msigdb <- enrich_msigdb(gene_module, universe = background_wgcna,
    category = "C2", subcategory = "CGP", species = "Mus musculus")

# Specific gene sets (category not required)
ifn_msigdb <- enrich_msigdb(gene_module, universe = background_wgcna,
    gene_sets = c("HALLMARK_INTERFERON_ALPHA_RESPONSE",
        "HALLMARK_INTERFERON_GAMMA_RESPONSE"),
    minGSSize = 5, species = "Mus musculus", db_species = "MM")
```

Enrichment of other gene sets (e.g., drug signatures, TF targets, pathways) can
be performed with `enrichR_list()`, which wraps `enrichR::enrichr()`.
Available databases can be checked with `enrichR::listEnrichrDbs()`.

`enrichGO_list()`, `enrich_msigdb()` and `enrichR_list()` all
return a named list of data.frames (one per category for `enrichGO_list`,
one for `enrich_msigdb()`, one per database for `enrichR_list()`) with
`Cluster` and `Description` columns, compatible with `plot_modules_h()`.

### Integrated visualisation of WGCNA modules and functional enrichment

Enriched terms can be added to the module profile plot with
`plot_modules_h()`, which is similar to `plot_modules_v()` above but
horizontally aligned with annotation of the modules with their top enriched
GO terms. The plotting function is adapted from ClusterGVis
package [@ClusterGVis].

```{r plot-wgcna-modules-go, fig.width = 8, fig.height = 5}
plot_modules_h(gene_module,
    data_obj_merged, assay = 2, scale = TRUE,
    ylabel = "Z-score of log2(CPM + 1)",
    enrich_list = go_list$all,
    enrich_category = "BP",
    heatmap_width = 6,
    heatmap_height = 8
)
```

Multiple enrichment categories (e.g. GO BP, CC, MF, MSigDB collections,
enrichR databases) can be displayed side-by-side
by passing a vector to `enrich_category`. Features of interest can
be marked on the heatmap or shown as an annotation column via
`mark_features` and `enrich_category = "Hub features"`. The
`enrich_threshold` parameter controls per-category p-value colour-coding
of enrichment terms; set to `NA` for "Hub features", which use auto-assigned
colours instead.

Hub features can be extracted by WGCNA module membership (the correlation of
the feature with the corresponding module eigengene) with `extract_hubs()`

```{r plot-wgcna-modules-multi, fig.width = 15, fig.height = 5}
hub_features <- extract_hubs(net, top_n = 3)

plot_modules_h(gene_module,
    data_obj_merged, assay = 2, scale = TRUE,
    ylabel = "Z-score of log2(CPM + 1)",
    enrich_list = c(go_list$all, hallmark_msigdb),
    enrich_category = c("BP", "CC", "MH"),
    enrich_threshold = c(0.05, 0.05, 0.05),
    mark_features = hub_features,
    heatmap_width = 6,
    heatmap_height = 8
)
```

### Universal: plot features of interest

Selected features can be plotted with `plot_trend()`, which shows the mean and
standard deviation of replicates at each time point.

```{r plot-feature}
plot_trend(data_obj, assay = 1,
    features = c("Abtb1", "Dram1", "Ifi27", "Nufip1"),
    title = "Example features")

# Or pre-calculate mean and sd with calc_mean_sd()
table_mean_sd_orig <- calc_mean_sd(data_obj)$orig
plot_trend(table_mean_sd_orig,
    features = c("Abtb1", "Dram1", "Ifi27", "Nufip1"),
    title = "Example features")
```

Visualisation of features with high & low residual variance justify the
filtering strategy for WGCNA input.

```{r plot-residual}
plot_trend(data_obj, assay = 1,
    features = var_decomp |>
        dplyr::arrange(Residual) |> utils::head(12) |> dplyr::pull(Feature),
    title = "Features with lowest residual variance"
)

plot_trend(data_obj, assay = 1,
    features = var_decomp |>
        dplyr::arrange(dplyr::desc(Residual)) |>
        utils::head(12) |> dplyr::pull(Feature),
    title = "Features with highest residual variance"
)
```

### One-step preprocessing with `prepare_tide()`

`prepare_tide()` is a wrapper function that runs multiple preprocessing steps
in a single call, including input creation, normalisation, merging replicates,
variance decomposition, and feature filtering. It is convenient for standard
analyses, while running each step individually allows for more control over
parameters.

```{r prepare-tide}
tide <- prepare_tide(
    data = tutorial_data,
    sample_ann = tutorial_sample_info,
    filter_ratio = 0.5,
    min_groups = 2,
    keep = "below_quantile",
    residual_threshold = 0.75
)
tide$filter_summary
```

Variance decomposition and filtering are run on both `"orig"` and `"norm"`
assays. The returned list `tide` contains:

| Element | Description | Use for | Step-by-step equivalent in the tutorial |
|---------|-------------|---------|-------------------------|
| `tide$se` | Unmerged SE, all features | Correlation matrix, PCA, UMAP, Pairwise DE, plot_trend | `data_obj` |
| `tide$se_filtered$orig/norm` | Unmerged SE, filtered features | WGCNA | `data_obj_wgcna` |
| `tide$merged_list` | Per-group merged SEs, all features | calc_feature_property, run_Trendy | `data_obj_merged_list` |
| `tide$merged_list_filtered$orig/norm` | Per-group merged SEs, filtered features | - | - |
| `tide$merged_se` | Single merged SE, all features | plot_modules_v/h | `data_obj_merged` |
| `tide$merged_se_filtered$orig/norm` | Single merged SE, filtered features | plot_modules_v/h | - |
| `tide$variance$orig/norm` | Variance decomposition | plot_variance, enrichGO_rank | `var_decomp` |
| `tide$feature_property` | Feature property summary | group_specific_features | `property_tb` |
| `tide$filter_summary` | Per-stage filtering stats | QC | - |
| `tide$group_specific_filter` | Features removed by cross-group filter | Review group-specific features | - |
| `tide$DE` | DE results placeholder | Centralised storage | `DE_between_time/group_out`|
| `tide$enrichment` | Enrichment results placeholder | Centralised storage | `go_list`, `gse_group`, `hallmark_msigdb` etc. |
| `tide$WGCNA` | WGCNA results placeholder | Centralised storage | `net` |

Pairwise DE, WGCNA and enrichment can be run on the returned elements, and the
results can be appended to the list for centralised storage.

```{r tide-de}
tide$DE <- list(
    between_group = DE_between_group_out,
    between_time  = DE_between_time_out
)
tide$WGCNA <- net

tide$enrichment$GO_modules <- go_list
tide$enrichment$MSigDB <- hallmark_msigdb
tide$enrichment$GO_group <- gse_group
```

## Interoperability within Bioconductor ecosystem

TiDEomics uses SummarizedExperiment as the central data container. It keeps
the feature x sample matrices, sample annotation, and feature metadata
synchronized during normalisation, splitting, and filtering. The use of
multiple assays within the same object also allows original ("orig") and
baseline-normalized ("norm") data to remain linked throughout the workflow.

As SummarizedExperiment is a standard Bioconductor container, TiDEomics
data can be passed to other Bioconductor tools for further custom analysis
and visualisation. For example, [iSEE](https://bioconductor.org/packages/iSEE)
interactive SummarizedExperiment browser can be used to explore the data
in an interactive Shiny app.

```{r isee, eval = FALSE}
iSEE::iSEE(data_obj)
```

`prepare_tide()` output can be converted to a
[`DeeDeeExperiment`](https://bioconductor.org/packages/DeeDeeExperiment)
object for formatted summary. The nested DE and enrichment
results can be flattened with `flatten_DE()` or `flatten_enrich()` before
passing to the `DeeDeeExperiment` constructor.

```{r deedee}
# Flatten nested results
de_flat <- flatten_DE(tide$DE)
enrich_flat <- flatten_enrich(tide$enrichment)

# Build DeeDeeExperiment
dde <- DeeDeeExperiment::DeeDeeExperiment(
    sce = tide$se,
    de_results = de_flat,
    enrich_results = enrich_flat
)
```

## Common usage scenarios

TiDEomics supports two experimental designs, auto-detected from the sample
annotation:

- **Cell culture / independent samples**: `create_input(data, sample_ann)`
    -- each sample is independent. LMM uses `(1|Group) + (1|Time)`. DE uses
    unpaired designs.
- **Patient / repeated-measures**:
    `create_input(data, sample_ann, subject_col = "PatientID")` -- same
    subjects tracked across time points. LMM adds `(1|Subject)`.
    `DE_between_time()` uses paired design via `limma::duplicateCorrelation()`.
    `normalise_to_start(by_subject = TRUE)` subtracts each subject's own
    baseline. `merge_replicates()` and `plot_trend()` compute statistics across
    subjects.

Users can ask different biological questions and use different functions in
TiDEomics to answer them. For example:

- Which features change over time in each sample group? -> `DE_between_time()`,
    `calc_feature_property()`, `run_Trendy()`.

- Which features differ between groups at the same time point? ->
    `DE_between_group()` (For starting time point, use `assay = 1`; for later
    time points, use `assay = 1` for original input data, or use `assay = 2` for
    change-from-baseline).

- What is the temporal pattern for individual features? -> `run_Trendy()` +
    `summarise_Trendy()` + `extract_segment_trends()`.

- Which features are noisy vs biologically driven? -> `decomp_variance()` +
    `plot_variance()`, noisiness can be represented by `Residual` variance.

- What co-expression modules are present and what are their temporal profiles?
    -> `prepare_WGCNA()` + `run_WGCNA()` + `plot_modules_v()` or
    `plot_modules_h()`.

- What biological processes are associated with the genes of interest (e.g.
    co-expression modules, differentially expressed genes in each sample
    group)? -> `enrichGO_list()` + `plot_GO()`.

- What biological processes are associated with the genes ranked by their
    property (e.g. time or group contributed variance)? -> `enrichGO_rank()`.

- How to normalise the data? -> `normalise_to_start()` to focus on changes
    from baseline; apply global normalisation (quantile, median) and batch
    correction when appropriate, before `create_input()`.

- How to handle missing values? -> `impute_groups()` for `run_Trendy`.
    `plot_pca()` and `plot_umap()` automatically exclude features with missing
    values. `prepare_WGCNA()` filter out features with too many missing values.

## Session information

```{r session-info}
sessionInfo()
```

## References