---
title: "scProkaR core workflow for microbial single-cell RNA-seq"
author:
- name: Ning Liu
  affiliation: University of Adelaide
- name: Jinxin Zhao
  affiliation: Monash University
date: "`r Sys.Date()`"
package: scProkaR
output:
  BiocStyle::html_document:
    toc: true
    toc_float: true
    toc_depth: 2
abstract: |
  This vignette demonstrates scProkaR using a bacterial single-cell RNA-seq
  time course with two polymyxin B concentrations and an untreated control.
  It covers data import, quality control, integration, clustering,
  visualisation and pseudobulk differential expression with edgeR. Cell-level
  data are stored in a SingleCellExperiment. The integration and
  differential-expression sections explain the limitations of the example
  design.
vignette: >
  %\VignetteIndexEntry{1. scProkaR core workflow}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>",
    crop = NULL,
    fig.width = 6,
    fig.height = 4.5,
    dpi = 96,
    fig.retina = 1
)
```

# Introduction

`scProkaR` provides tools for analysing bacterial single-cell RNA sequencing
(scRNA-seq) data. The workflow covers quality control, integration,
clustering and differential expression in a polymyxin B time course.

Bacterial data require quality-control criteria suited to the organism and
assay. Relevant measurements include library size, detected genes, and the
fractions of counts assigned to ribosomal RNA and ribosomal-protein genes.
Gene identifiers may include both locus tags and gene symbols.

The main functions used in this workflow are:

* `CreateBacObject()` converts count matrices, 10x Genomics output,
    `SingleCellExperiment` objects and Seurat objects into a
    `SingleCellExperiment`.
* `RunBacQC()` computes per-cell quality-control (QC) metrics, and
    `FilterBacCells()` applies user-specified thresholds.
* `IntegrateBacData()` computes integration embeddings;
    `RegisterIntegrationEmbedding()` and
    `RegisterExistingIntegrationEmbeddings()`
    register externally computed embeddings.
* `BenchmarkIntegration()` compares embeddings using scIB-inspired metrics.
* `aggregate_pseudobulk()`, `normalize_pseudobulk()`,
    `run_edger_pairwise_de()` and `run_edger_spline_de()` support pseudobulk
    differential expression with `r Biocpkg("edgeR")`.
* Time Aware Trajectory Analysis (TATA) is described in the companion vignette.

The cell-level workflow uses `r Biocpkg("SingleCellExperiment")` as its main
data container. Counts are stored in `counts`, log-normalised expression in
`logcounts`, cell annotations in `colData()`, and embeddings in `reducedDims()`.
Processing records are stored in `metadata(sce)$scProkaR`.

These objects can be used with `r Biocpkg("scater")`, `r Biocpkg("scran")`
and other Bioconductor tools that accept a `SingleCellExperiment`.
scProkaR uses `r Biocpkg("scrapper")` for feature selection when available,
`r Biocpkg("batchelor")` for MNN correction, and `r Biocpkg("edgeR")` for
pseudobulk testing. Benchmarking, modelling and plotting functions return
other result types, as described in their help pages.

# Installation

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

```{r load-packages}
library(scProkaR)
library(SingleCellExperiment)
library(ggplot2)
```

# The example data

The package includes `sce1`, a bacterial scRNA-seq dataset with untreated
cells and cells exposed to two polymyxin B (PMB) concentrations.
The available preparation notes are in `inst/scripts/make-sce1.R`; source
accession and protocol details still need to be completed there.

```{r load-data}
data("sce1", package = "scProkaR")
sce1
```

The design is a dose-by-time grid:

```{r describe-design}
table(sce1$treatments, sce1$timepoints)
```

The untreated `control` was sampled at 0, 1, 4 and 7 hours.
The `PMB0.5` and `PMB2` treatment groups were sampled at 1, 4 and 7 hours;
only the control has a time-zero sample. Each of the ten samples contains
500 cells, giving 5000 cells in total. `colData()` records the sample,
treatment, timepoint and precomputed `clusters` labels. The labels provide
a reference for the integration examples; their biological interpretation
depends on the source annotation.

```{r describe-coldata}
head(as.data.frame(colData(sce1)), 4)
```

The 3722 feature identifiers include both gene symbols and locus tags:

```{r describe-features}
head(rownames(sce1), 8)
```

Inspect the per-cell library sizes before choosing thresholds for filtering:

```{r describe-depth}
summary(Matrix::colSums(assay(sce1, "counts")))
```

## Subsetting for this vignette

To reduce build time, this vignette samples 250 cells from each of the ten
samples, retaining 2500 cells in total. Equal sample sizes are preserved,
but subsampling can change feature selection, neighbourhoods, clusters and
statistical estimates. Assess the full dataset for the final analysis.

```{r subset-cells}
set.seed(101)
keep_cells <- unlist(lapply(
    split(seq_len(ncol(sce1)), sce1$sample),
    function(idx) sample(idx, 250)
))
sce_small <- sce1[, sort(keep_cells)]
dim(sce_small)
```

# Creating a scProkaR object

The example below supplies a count matrix and a per-cell metadata
`data.frame` to `CreateBacObject()`. The `sample_col`, `batch_col`,
`condition_col` and `time_col` arguments identify the corresponding metadata
columns. Here, treatment is also used as the batch variable; the integration
section explains the resulting confounding.

Setting `run_unintegrated = TRUE` runs log-normalisation, feature selection,
principal component analysis (PCA) and graph-based clustering. When
`r CRANpkg("uwot")` is installed, it also stores a UMAP embedding as
`umap.unintegrated`. This provides an uncorrected reference for the
integration examples.

```{r create-object}
counts_matrix <- assay(sce_small, "counts")
cell_meta <- as.data.frame(colData(sce_small))

set.seed(101)
sce <- CreateBacObject(
    counts_matrix,
    cell_metadata = cell_meta,
    sample_col = "sample",
    batch_col = "treatments",
    condition_col = "treatments",
    time_col = "timepoints",
    organism = "bacteria",
    run_unintegrated = TRUE,
    unintegrated_dims = seq_len(20),
    unintegrated_k = 20,
    unintegrated_resolution = 1
)
sce
```

The result is an ordinary `SingleCellExperiment`, with the analysis history
recorded in its metadata:

```{r inspect-object}
reducedDimNames(sce)
names(S4Vectors::metadata(sce)$scProkaR)
```

## Other entry points

Two further input types are accepted by the same function.

**10x Genomics directories.** Pass the path to a CellRanger output folder
(the one containing `matrix.mtx.gz`, `barcodes.tsv.gz` and `features.tsv.gz`,
whether or not it is nested inside `filtered_feature_bc_matrix/`) and
scProkaR reads it directly, using `r Biocpkg("DropletUtils")` when that
package is available and falling back to `r CRANpkg("Matrix")` otherwise:

```{r tenx-entry-point}
tenx_call <- quote(CreateBacObject(
    "path/to/filtered_feature_bc_matrix",
    sample_id_value = "PMB2_7",
    sample_col = "sample"
))
tenx_call
```

**Seurat objects.** `CreateBacObject()` imports the requested assay and
layer, cell metadata and existing reductions. The following example runs
when the optional `r CRANpkg("SeuratObject")` package is installed:

```{r seurat-entry-point}
if (requireNamespace("SeuratObject", quietly = TRUE)) {
    demo_cells <- seq_len(200)
    seu <- SeuratObject::CreateSeuratObject(
        counts = counts_matrix[, demo_cells],
        meta.data = cell_meta[demo_cells, , drop = FALSE]
    )
    sce_from_seurat <- CreateBacObject(
        seu,
        seurat_assay = "RNA",
        seurat_layer = "counts",
        sample_col = "sample",
        batch_col = "treatments"
    )
    print(dim(sce_from_seurat))
    print(head(colnames(colData(sce_from_seurat))))
} else {
    message("SeuratObject is not installed; skipping the Seurat entry point.")
}
```

A `Layer 'data' is empty` warning may appear because this demonstration
object contains raw counts without a normalised data layer. Normalisation
is performed later when required by the workflow.

Each supported input is converted to a `SingleCellExperiment`.
Use `MergeBacObjects()` to combine objects prepared separately for each sample.
It aligns features, makes duplicate cell identifiers unique, retains the
original identifiers in `colData()$original_cell_id`, and carries forward
compatible reductions.

# Quality control

`RunBacQC()` writes per-cell QC metrics to `colData()`, including
`total_counts`, `detected_features` and two ribosomal fractions:

* **`rrna_fraction`** is the fraction of counts assigned to ribosomal RNA
    genes. Its interpretation depends on the annotation and rRNA depletion
    protocol; a high value alone does not establish cell damage.
* **`ribo_fraction`** is the fraction assigned to ribosomal-protein genes
    such as `rpl*` and `rps*`. Inspect its distribution across conditions before
    deciding whether it is relevant to filtering.

Genes are classified by regular expression on the row names by default
(`rrna_pattern = "^(rrs|rrl|rrf)"`), which matches the standard bacterial
rRNA nomenclature. If your annotation carries an explicit gene-class column in
`rowData()`, pass its name as `gene_class_col` and it will be used instead.

```{r run-qc}
sce <- RunBacQC(sce)
head(as.data.frame(colData(sce))[, c(
    "total_counts",
    "detected_features",
    "rrna_fraction",
    "ribo_fraction"
)], 4)
```

No feature name in `sce1` matches the default rRNA pattern, so the computed
`rrna_fraction` is zero for every cell. This naming check cannot establish
whether rRNA was absent from the original libraries or removed during
processing.

```{r check-rrna}
sum(grepl("^(rrs|rrl|rrf)", rownames(sce), ignore.case = TRUE))
summary(sce$rrna_fraction)
```

The following summary compares median ribosomal-protein fractions across
the three treatment groups:

```{r summarise-ribo}
round(tapply(sce$ribo_fraction, sce$treatments, median), 4)
```

These medians pool cells across timepoints and do not isolate a dose effect.
The example therefore retains `ribo_fraction` for inspection and does not
use it to filter cells.

`FilterBacCells()` applies thresholds to the QC metrics. The rRNA cutoff
removes no cells here because all computed rRNA fractions are zero.
Filtering in this example therefore depends on counts and detected features.

```{r filter-cells}
sce <- FilterBacCells(
    sce,
    min_counts = 50,
    min_features = 30,
    max_rrna_fraction = 0.2
)
dim(sce)
table(sce$treatments)
```

Choose thresholds from the distributions and assay characteristics of your
own data. The values above illustrate filtering for this dataset and are
not universal cutoffs.

# Integration across treatment groups

This example uses `treatments` as the batch variable, making batch and
condition identical in the analysis. Any removed differences may therefore
include the drug response. The examples demonstrate the integration functions;
they cannot distinguish technical batch effects from treatment effects.
Use recorded technical batches, with conditions represented across batches
where possible, when evaluating correction in an experimental study.

## MNN correction

`IntegrateBacData(method = "mnn")` delegates to the mutual-nearest-neighbours
correction in `r Biocpkg("batchelor")`. That package is suggested rather than
required, so the call is guarded:

```{r integrate-mnn}
if (requireNamespace("batchelor", quietly = TRUE)) {
    set.seed(101)
    sce <- IntegrateBacData(
        sce,
        batch_col = "treatments",
        method = "mnn",
        dims = seq_len(20)
    )
    print(reducedDimNames(sce))
} else {
    message("batchelor is not installed; skipping MNN integration.")
}
```

## Harmony correction

`method = "harmony"` calls `r CRANpkg("harmony")`. The example runs when
this optional package is installed:

```{r integrate-harmony}
if (requireNamespace("harmony", quietly = TRUE)) {
    set.seed(101)
    sce <- IntegrateBacData(
        sce,
        batch_col = "treatments",
        method = "harmony",
        dims = seq_len(20)
    )
    print(reducedDimNames(sce))
} else {
    message("harmony is not installed; skipping Harmony integration.")
}
```

Each backend writes its embedding to `reducedDims()` under
`integrated_<method>` and records the call in `metadata(sce)$scProkaR`.

## Registering your own embedding

`RegisterIntegrationEmbedding()` registers a cell-by-dimension matrix from
an external analysis for use in downstream integration comparisons.
The matrix must correspond to the cells in the object.

Here, the registered matrix contains PCA coordinates centred within each
treatment group. This removes differences in group means and illustrates
registration without an additional integration package. As above, centring
by treatment can also remove biological differences.

```{r register-embedding}
pca_emb <- reducedDim(sce, "PCA")[, seq_len(20), drop = FALSE]
batch_vec <- as.character(sce$treatments)
batch_centres <- t(vapply(
    split(seq_len(nrow(pca_emb)), batch_vec),
    function(idx) colMeans(pca_emb[idx, , drop = FALSE]),
    numeric(ncol(pca_emb))
))
centred_emb <- pca_emb - batch_centres[batch_vec, , drop = FALSE]

sce <- RegisterIntegrationEmbedding(
    sce,
    embedding = centred_emb,
    method_name = "centred"
)
reducedDimNames(sce)
```

For an object that already contains external latent spaces, use
`RegisterExistingIntegrationEmbeddings(sce, reductions = "X_scVI")`, or pass
several names such as `c("seurat_cca", "X_scVI")`. The default copies them to
`integrated_cca` and `integrated_scvi`; `copy = FALSE` registers the original
reductions directly. Omitting `reductions` detects likely integration names.
This registers existing coordinates; it does not run the external methods.

Both registration functions protect existing embeddings from replacement
unless `overwrite = TRUE` is supplied. Registering the same values again is
allowed. `BenchmarkIntegration()` includes registered methods and other
`integrated_*` reductions by default, and accepts names without regard to case
(for example, `methods = "SCVI"`).

The remaining examples use MNN when its embedding is available and the
centred embedding otherwise. This choice keeps the examples runnable across
installations; it does not establish which integration method is best.

```{r choose-reduction}
integrated_names <- grep(
    "^integrated_",
    reducedDimNames(sce),
    value = TRUE
)
main_reduction <- if ("integrated_mnn" %in% integrated_names) {
    "integrated_mnn"
} else {
    "integrated_centred"
}
main_reduction
```

# UMAP and clustering on the integrated space

`RunIntegratedUMAP()` runs a UMAP on any stored embedding and saves it as
`umap_<reduction>`. It requires `r CRANpkg("uwot")`, so we guard it and fall
back to plotting the first two dimensions of the integrated embedding itself
when `uwot` is unavailable.

```{r run-umap}
if (requireNamespace("uwot", quietly = TRUE)) {
    set.seed(101)
    sce <- RunIntegratedUMAP(
        sce,
        reduction = main_reduction,
        n_neighbors = 20
    )
    vis_reduction <- paste0("umap_", main_reduction)
} else {
    message("uwot is not installed; plotting the embedding directly.")
    vis_reduction <- main_reduction
}
vis_reduction
```

The uncorrected counterpart for side-by-side comparison is chosen the same
way:

```{r choose-unintegrated}
unintegrated_reduction <- if (
    "umap.unintegrated" %in% reducedDimNames(sce)
) {
    "umap.unintegrated"
} else {
    "PCA"
}
unintegrated_reduction
```

`RunIntegratedClustering()` builds and partitions a shared-nearest-neighbour
graph from the selected embedding. `k` controls neighbourhood size and
`resolution` controls partition granularity. Compare the resulting clusters
with the uncorrected analysis before using them downstream.

```{r run-clustering}
set.seed(101)
sce <- RunIntegratedClustering(
    sce,
    reduction = main_reduction,
    cluster_col = "integrated_clusters",
    k = 20,
    algorithm = "louvain",
    resolution = 1
)
table(sce$integrated_clusters)
```

Cross-tabulation compares the inferred clusters with the supplied reference
labels. A reference group spread across several rows has been split;
several reference groups in one row have been merged.

```{r compare-clusters}
table(sce$integrated_clusters, sce$clusters)
```

# Visualisation

`PlotReduction()` draws any stored embedding, coloured by any `colData()`
column, and returns a plain `ggplot` object you can modify further.

```{r plot-batch}
PlotReduction(
    sce,
    reduction = vis_reduction,
    colour_by = "treatments",
    point_size = 0.5
)
```

```{r plot-labels}
PlotReduction(
    sce,
    reduction = vis_reduction,
    colour_by = "clusters",
    point_size = 0.5
)
```

Time can be overlaid the same way. Coercing the timepoint to a factor with an
explicit level order keeps the legend in chronological rather than
lexicographic order:

```{r plot-time}
sce$timepoint_factor <- factor(
    sce$timepoints,
    levels = c("0", "1", "4", "7")
)
PlotReduction(
    sce,
    reduction = vis_reduction,
    colour_by = "timepoint_factor",
    point_size = 0.5
)
```

`facet_by` displays each treatment group using the same embedding
coordinates, allowing occupied regions to be compared across groups:

```{r plot-facet}
PlotReduction(
    sce,
    reduction = vis_reduction,
    colour_by = "integrated_clusters",
    facet_by = "treatments",
    point_size = 0.5
)
```

`PlotIntegrationOverview()` assembles the standard before/after panel set in
one call and returns them as a named list.

```{r integration-overview}
overview <- PlotIntegrationOverview(
    sce,
    unintegrated_reduction = unintegrated_reduction,
    integrated_reduction = vis_reduction,
    batch_col = "treatments",
    label_col = "clusters",
    split_by = "treatments",
    point_size = 0.5
)
names(overview)
```

```{r overview-before-batch}
overview$unintegrated_by_batch
```

```{r overview-after-batch}
overview$integrated_by_batch
```

```{r overview-after-label}
overview$integrated_by_label
```

Compare `unintegrated_by_batch` and `integrated_by_batch` for changes in
group separation. Use `integrated_by_label` to inspect the reference labels
and `integrated_split` to compare occupied regions across treatments.
Because treatment is the batch variable here, increased mixing may reflect
loss of the drug response.

# Benchmarking the integration

`BenchmarkIntegration()` complements the plots with metrics of batch removal
and biological conservation, then ranks the candidate embeddings.

The metrics are inspired by scIB (Luecken et al., 2022), but the
implementations differ from the reference scIB software. Compare methods
within the same run; the reported values should not be presented as
reference scIB scores or compared directly across studies.

```{r run-benchmark}
bench_methods <- c("PCA", integrated_names)
set.seed(101)
bench <- BenchmarkIntegration(
    sce,
    batch_col = "treatments",
    label_col = "clusters",
    methods = bench_methods
)
bench$ranking
```

The two families are:

* **Batch removal** includes `pcr_batch`, one minus the mean fraction of
    embedding-coordinate variance explained by batch; `ilisi`, local batch
    mixing; and `kbet_like`, neighbourhood agreement with overall batch
    composition. This category also includes `graph_connectivity`, which
    measures connectivity within reference labels in the nearest-neighbour
    graph.
* **Biological conservation** includes `clisi`, local label purity; `ari`
    and `nmi`, agreement between k-means clusters and reference labels; and
    `silhouette`, separation of reference labels in the embedding.

The `score` column is oriented so higher values indicate better performance
under each metric. `ranking` averages the available scores for each method;
the `value` column retains each metric's returned value.

```{r benchmark-scores}
head(bench$scores, 8)
```

```{r benchmark-overall}
bench$plots$overall
```

The trade-off plot displays batch removal on the x-axis and biological
conservation on the y-axis. Higher values indicate better scores in each
category. With treatment used as batch, these scores alone cannot establish
whether biological variation has been preserved.

```{r benchmark-tradeoff}
bench$plots$tradeoff
```

```{r benchmark-key-metrics}
bench$plots$key_metrics
```

The remaining panels, `bench$plots$batch_removal`,
`bench$plots$bio_conservation` and `bench$plots$heatmap`, show the scores
by category and individual metric.

Inspect individual metrics alongside the experimental design and plots
before selecting an embedding. The overall ranking averages different goals
and cannot resolve the batch-treatment confounding in this example.

# Pseudobulk differential expression

Pseudobulk analysis sums cell counts within independent biological replicates
and tests the resulting sample-level counts with `r Biocpkg("edgeR")`.
The experimental replicate remains the unit of inference.

## Defining replicates

The packaged metadata identify one sample per treatment and timepoint,
without independent biological replicate labels. To demonstrate the code,
the cells from each sample are split into three artificial groups.
These groups do not provide biological replication, so the resulting
p-values and false discovery rates (FDRs) cannot support biological conclusions.
For experimental analyses, aggregate cells from independent samples such as
separately grown cultures.

```{r assign-replicates}
set.seed(101)
replicate_id <- integer(ncol(sce))
for (idx in split(seq_len(ncol(sce)), sce$sample)) {
    replicate_id[idx] <- sample(rep(seq_len(3), length.out = length(idx)))
}
sce$replicate <- factor(replicate_id)
sce$pb_id <- paste(sce$sample, sce$replicate, sep = "_r")
table(sce$sample, sce$replicate)
```

## Aggregating and normalising

`aggregate_pseudobulk()` sums counts for each unique combination of
`sample_cols` and stores the corresponding annotations in the result's
`colData()`. It returns a `SingleCellExperiment`; `min_cells` excludes
aggregates with fewer than the specified number of cells.

```{r aggregate-pseudobulk}
pb <- aggregate_pseudobulk(
    sce,
    sample_cols = c(
        "pb_id", "sample", "replicate", "treatments", "timepoints"
    ),
    aggregation = "sum",
    min_cells = 10
)
dim(pb)
head(as.data.frame(colData(pb))[, c(
    "pb_id", "treatments", "timepoints", "ncells"
)], 4)
```

`normalize_pseudobulk()` computes edgeR normalisation factors (TMM by
default), stores them and the library sizes in `colData()`, and adds a
`logcounts` assay of log-CPM values for visualisation. The statistical tests
below re-normalise internally from raw counts, so this step is for inspection
and plotting rather than a prerequisite.

```{r normalize-pseudobulk}
pb <- normalize_pseudobulk(pb, method = "TMM")
summary(pb$norm.factors)
```

PCA of the pseudobulk profiles helps identify sample differences and
outliers. Here, the artificial groups share their source library and should
not be interpreted as independent replicates.

```{r pseudobulk-pca}
pb_pca <- stats::prcomp(
    t(as.matrix(assay(pb, "logcounts"))),
    scale. = TRUE
)
reducedDim(pb, "pb_pca") <- pb_pca$x[, seq_len(4), drop = FALSE]
pb$timepoint_factor <- factor(pb$timepoints, levels = c("0", "1", "4", "7"))
PlotReduction(pb, reduction = "pb_pca", colour_by = "treatments",
    point_size = 3)
```

```{r pseudobulk-pca-time}
PlotReduction(pb, reduction = "pb_pca", colour_by = "timepoint_factor",
    point_size = 3)
```

## Pairwise contrasts

`run_edger_pairwise_de()` fits an edgeR quasi-likelihood negative-binomial
model and extracts the requested contrasts. Here, each timepoint is compared
with hour 0 in an additive model that includes treatment as a covariate.

Hour 0 is represented only by the control group. The time contrasts
therefore rely on the additive model assumption and do not estimate separate
treatment-specific time responses. The artificial grouping above further
limits these results to demonstration.

```{r pairwise-de}
pairwise_time <- run_edger_pairwise_de(
    pb,
    group_col = "timepoints",
    covariates = "treatments",
    contrast_type = "reference",
    reference_level = "0",
    normalization = "TMM",
    fdr_cutoff = 0.05
)
names(pairwise_time$tables)
```

`$tables` contains one `data.frame` per contrast; `$combined_table` combines
them with a `contrast` column. edgeR's `filterByExpr()` removes genes with
insufficient expression for the design, so the tables can contain fewer
than the 3722 input features.

```{r pairwise-summary}
table(
    pairwise_time$combined_table$contrast,
    pairwise_time$combined_table$significant
)
```

```{r pairwise-top-genes}
tab_4_vs_0 <- pairwise_time$tables[["4_vs_0"]]
head(tab_4_vs_0[order(tab_4_vs_0$FDR), c("gene", "logFC", "FDR")], 10)
```

The volcano plot displays log-fold change against statistical significance.
The guide lines mark the fold-change and FDR thresholds; results here remain
illustrative because the groups are not biological replicates.

```{r pairwise-volcano}
plot_pairwise_de_volcano(
    pairwise_time,
    contrast = "4_vs_0",
    fdr_cutoff = 0.05,
    lfc_cutoff = 1,
    top_n_labels = 8
)
```

The MA plot displays log-fold change against average abundance, helping
identify large estimated changes among low-expression genes:

```{r pairwise-ma}
plot_pairwise_de_ma(
    pairwise_time,
    contrast = "4_vs_0",
    fdr_cutoff = 0.05
)
```

To compare composite groups, such as PMB2 versus control at 7 hours,
combine the grouping columns and supply the contrast directly:

```{r pairwise-composite}
pairwise_dose <- run_edger_pairwise_de(
    pb,
    group_col = "treatments",
    combine_group_cols = c("treatments", "timepoints"),
    contrast_type = "manual",
    manual_pairs = data.frame(
        group1 = c("PMB2__7", "PMB0.5__7"),
        group2 = c("control__7", "control__7"),
        stringsAsFactors = FALSE
    ),
    normalization = "TMM"
)
names(pairwise_dose$tables)
```

```{r pairwise-composite-volcano}
plot_pairwise_de_volcano(
    pairwise_dose,
    contrast = "PMB2__7_vs_control__7",
    fdr_cutoff = 0.05,
    lfc_cutoff = 1,
    top_n_labels = 8
)
```

## Time-course testing with splines

`run_edger_spline_de()` models time with a natural cubic spline basis and
jointly tests its coefficients for evidence of temporal variation.
The model can represent non-monotonic responses within the chosen spline
basis. With `condition_col` specified, it fits each treatment group separately.

Choose spline degrees of freedom (`df`) in relation to the number of
observed timepoints and the model design. The treated groups have three
timepoints, so this example uses `df = 2`. Check model rank and residual
degrees of freedom before interpreting tests.

```{r spline-de}
spline_de <- run_edger_spline_de(
    pb,
    time_col = "timepoints",
    condition_col = "treatments",
    df = 2,
    normalization = "TMM",
    return_curve_for = "significant",
    fdr_cutoff = 0.05
)
names(spline_de$by_condition)
```

```{r spline-summary}
table(
    spline_de$combined_table$treatments,
    spline_de$combined_table$significant
)
```

The counts above summarise genes passing the selected FDR cutoff in
each fit. They illustrate the output format; artificial replicate groups
prevent interpretation of these counts as evidence of a treatment effect.

```{r spline-top-genes}
spline_top <- spline_de$combined_table[
    order(spline_de$combined_table$FDR),
    c("gene", "treatments", "FDR")
]
head(spline_top, 10)
```

`plot_time_series_deg_curves()` displays the expression curves for selected
genes, allowing the temporal patterns to be inspected:

```{r spline-curves}
curve_genes <- head(unique(spline_top$gene), 6)
plot_time_series_deg_curves(spline_de, genes = curve_genes, ncol = 3)
```

If you prefer a single call, `run_scproka_de()` wraps aggregation,
normalisation and testing into one function.

# Where to go next

The companion vignette, `vignette("tata-workflow", package = "scProkaR")`,
covers Time Aware Trajectory Analysis for perturbation time courses.

The cell-level `SingleCellExperiment` can also be used with
`r Biocpkg("scater")` for plotting, `r Biocpkg("scran")` for marker analysis,
and other tools that accept this container.

# Further reading

* Amezquita RA, Lun ATL, Becht E, *et al.* (2020). Orchestrating single-cell
    analysis with Bioconductor. *Nature Methods* 17, 137-145.
    <https://bioconductor.org/books/release/OSCA/>
* Robinson MD, McCarthy DJ, Smyth GK (2010). edgeR: a Bioconductor package for
    differential expression analysis of digital gene expression data.
    *Bioinformatics* 26, 139-140.
    <https://bioconductor.org/packages/edgeR>
* Chen Y, Lun ATL, Smyth GK (2016). From reads to genes to pathways:
    differential expression analysis of RNA-Seq experiments using Rsubread and
    the edgeR quasi-likelihood pipeline. *F1000Research* 5, 1438.
* Haghverdi L, Lun ATL, Morgan MD, Marioni JC (2018). Batch effects in
    single-cell RNA-sequencing data are corrected by matching mutual nearest
    neighbors. *Nature Biotechnology* 36, 421-427.
    <https://bioconductor.org/packages/batchelor>
* Korsunsky I, Millard N, Fan J, *et al.* (2019). Fast, sensitive and accurate
    integration of single-cell data with Harmony. *Nature Methods* 16,
    1289-1296. <https://CRAN.R-project.org/package=harmony>
* Luecken MD, Buttner M, Chaichoompu K, *et al.* (2022). Benchmarking
    atlas-level data integration in single-cell genomics. *Nature Methods* 19,
    41-50. The metrics in `BenchmarkIntegration()` are inspired by this work but
    are independent reimplementations.
* McInnes L, Healy J, Melville J (2018). UMAP: Uniform Manifold Approximation
    and Projection for dimension reduction. *arXiv:1802.03426*.
    <https://CRAN.R-project.org/package=uwot>
* Lun ATL, McCarthy DJ, Marioni JC (2016). A step-by-step workflow for
    low-level analysis of single-cell RNA-seq data with Bioconductor.
    *F1000Research* 5, 2122. <https://bioconductor.org/packages/scran>
* Imdahl F, Vafadarnejad E, Homberger C, Saliba AE, Vogel J (2020).
    Single-cell RNA-sequencing reports growth-condition-specific global
    transcriptomes of individual bacteria. *Nature Microbiology* 5, 1202-1206.

# Session information

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