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 edgeR.
Time Aware Trajectory Analysis (TATA) is described in the companion vignette.

The cell-level workflow uses 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 scater, scran and other Bioconductor tools that accept a SingleCellExperiment. scProkaR uses scrapper for feature selection when available, batchelor for MNN correction, and edgeR for pseudobulk testing. Benchmarking, modelling and plotting functions return other result types, as described in their help pages.

Installation

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("scProkaR")

library(scProkaR)
library(SingleCellExperiment)
#> Loading required package: SummarizedExperiment
#> Loading required package: MatrixGenerics
#> Loading required package: matrixStats
#> 
#> Attaching package: 'MatrixGenerics'
#> The following objects are masked from 'package:matrixStats':
#> 
#>     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
#>     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
#>     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
#>     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
#>     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
#>     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
#>     colWeightedMeans, colWeightedMedians, colWeightedSds,
#>     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
#>     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
#>     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
#>     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
#>     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
#>     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
#>     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
#>     rowWeightedSds, rowWeightedVars
#> Loading required package: GenomicRanges
#> Loading required package: stats4
#> Loading required package: BiocGenerics
#> Loading required package: generics
#> 
#> Attaching package: 'generics'
#> The following objects are masked from 'package:base':
#> 
#>     as.difftime, as.factor, as.ordered, intersect, is.element, setdiff,
#>     setequal, union
#> 
#> Attaching package: 'BiocGenerics'
#> The following objects are masked from 'package:stats':
#> 
#>     IQR, mad, sd, var, xtabs
#> The following 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
#> 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(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.

data("sce1", package = "scProkaR")
sce1
#> class: SingleCellExperiment 
#> dim: 3722 5000 
#> metadata(0):
#> assays(1): counts
#> rownames(3722): dnaA ABUW-RS00010 ... ABUW-RS20075 ABUW-RS20460
#> rowData names(0):
#> colnames(5000): P1_1_GGAGCTGAGAGGTCCGTCTA P1_1_TATTGTCATCGGAGCGCTGC ...
#>   control4_TCGCGAAGCAGAGCTGCAAC control4_CGGATGCGAACGATTGATGA
#> colData names(4): sample clusters treatments timepoints
#> reducedDimNames(0):
#> mainExpName: NULL
#> altExpNames(0):

The design is a dose-by-time grid:

table(sce1$treatments, sce1$timepoints)
#>          
#>             0   1   4   7
#>   control 500 500 500 500
#>   PMB0.5    0 500 500 500
#>   PMB2      0 500 500 500

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.

head(as.data.frame(colData(sce1)), 4)
#>                             sample clusters treatments timepoints
#> P1_1_GGAGCTGAGAGGTCCGTCTA PMB0.5_1        F     PMB0.5          1
#> P1_1_TATTGTCATCGGAGCGCTGC PMB0.5_1        F     PMB0.5          1
#> P1_1_GTATCTTCATCGACCGTAAT PMB0.5_1        F     PMB0.5          1
#> P1_1_GCAGGCAATCGGAGCGATAT PMB0.5_1        F     PMB0.5          1

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

head(rownames(sce1), 8)
#> [1] "dnaA"         "ABUW-RS00010" "recF"         "gyrB"         "cybC"        
#> [6] "ABUW-RS00030" "ABUW-RS00035" "ABUW-RS00045"

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

summary(Matrix::colSums(assay(sce1, "counts")))
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>    28.0    74.0   127.0   226.9   244.0  2828.0

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.

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)
#> [1] 3722 2500

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 uwot is installed, it also stores a UMAP embedding as umap.unintegrated. This provides an uncorrected reference for the integration examples.

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
#> class: SingleCellExperiment 
#> dim: 3722 2500 
#> metadata(1): scProkaR
#> assays(2): counts logcounts
#> rownames(3722): dnaA ABUW-RS00010 ... ABUW-RS20075 ABUW-RS20460
#> rowData names(0):
#> colnames(2500): P1_1_TATTGTCATCGGAGCGCTGC P1_1_GTATCTTCATCGACCGTAAT ...
#>   control4_TTCGAGTAGAGCTCTGATAT control4_TCGCGAAGCAGAGCTGCAAC
#> colData names(5): sample clusters treatments timepoints
#>   unintegrated_clusters
#> reducedDimNames(2): PCA umap.unintegrated
#> mainExpName: NULL
#> altExpNames(0):

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

reducedDimNames(sce)
#> [1] "PCA"               "umap.unintegrated"
names(S4Vectors::metadata(sce)$scProkaR)
#> [1] "package"       "version"       "created"       "organism"     
#> [5] "columns"       "assays"        "integration"   "visualization"
#> [9] "unintegrated"

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 DropletUtils when that package is available and falling back to Matrix otherwise:

tenx_call <- quote(CreateBacObject(
    "path/to/filtered_feature_bc_matrix",
    sample_id_value = "PMB2_7",
    sample_col = "sample"
))
tenx_call
#> CreateBacObject("path/to/filtered_feature_bc_matrix", sample_id_value = "PMB2_7", 
#>     sample_col = "sample")

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

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.")
}
#> Warning: Layer 'data' is empty
#> [1] 3722  200
#> [1] "orig.ident"   "nCount_RNA"   "nFeature_RNA" "sample"       "clusters"    
#> [6] "treatments"

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.

sce <- RunBacQC(sce)
head(as.data.frame(colData(sce))[, c(
    "total_counts",
    "detected_features",
    "rrna_fraction",
    "ribo_fraction"
)], 4)
#>                           total_counts detected_features rrna_fraction
#> P1_1_TATTGTCATCGGAGCGCTGC          583               176             0
#> P1_1_GTATCTTCATCGACCGTAAT          377                98             0
#> P1_1_GCAGGTAAGGCGCCTGGAGC          538               168             0
#> P1_1_GGTCCAGTTGGCGAAGATGA         1544               359             0
#>                           ribo_fraction
#> P1_1_TATTGTCATCGGAGCGCTGC    0.09262436
#> P1_1_GTATCTTCATCGACCGTAAT    0.00265252
#> P1_1_GCAGGTAAGGCGCCTGGAGC    0.10966543
#> P1_1_GGTCCAGTTGGCGAAGATGA    0.09196891

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.

sum(grepl("^(rrs|rrl|rrf)", rownames(sce), ignore.case = TRUE))
#> [1] 0
summary(sce$rrna_fraction)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>       0       0       0       0       0       0

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

round(tapply(sce$ribo_fraction, sce$treatments, median), 4)
#> control  PMB0.5    PMB2 
#>  0.0615  0.0520  0.0843

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.

sce <- FilterBacCells(
    sce,
    min_counts = 50,
    min_features = 30,
    max_rrna_fraction = 0.2
)
dim(sce)
#> [1] 3722 2341
table(sce$treatments)
#> 
#> control  PMB0.5    PMB2 
#>     974     750     617

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 batchelor. That package is suggested rather than required, so the call is guarded:

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.")
}
#> [1] "PCA"               "umap.unintegrated" "integrated_mnn"

Harmony correction

method = "harmony" calls harmony. The example runs when this optional package is installed:

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.")
}
#> [1] "PCA"                "umap.unintegrated"  "integrated_mnn"    
#> [4] "integrated_harmony"

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.

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)
#> [1] "PCA"                "umap.unintegrated"  "integrated_mnn"    
#> [4] "integrated_harmony" "integrated_centred"

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.

integrated_names <- grep(
    "^integrated_",
    reducedDimNames(sce),
    value = TRUE
)
main_reduction <- if ("integrated_mnn" %in% integrated_names) {
    "integrated_mnn"
} else {
    "integrated_centred"
}
main_reduction
#> [1] "integrated_mnn"

UMAP and clustering on the integrated space

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

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
#> [1] "umap_integrated_mnn"

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

unintegrated_reduction <- if (
    "umap.unintegrated" %in% reducedDimNames(sce)
) {
    "umap.unintegrated"
} else {
    "PCA"
}
unintegrated_reduction
#> [1] "umap.unintegrated"

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.

set.seed(101)
sce <- RunIntegratedClustering(
    sce,
    reduction = main_reduction,
    cluster_col = "integrated_clusters",
    k = 20,
    algorithm = "louvain",
    resolution = 1
)
table(sce$integrated_clusters)
#> 
#>   1   2   3   4   5   6   7 
#> 433 307 395 395 334 333 144

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.

table(sce$integrated_clusters, sce$clusters)
#>    
#>       A   B   C   D   E   F
#>   1 129   1  41   0   9 253
#>   2  38   8 217   0  12  32
#>   3 235  28  37  11  76   8
#>   4 183  66 127   0   2  17
#>   5   4 266  19  14  29   2
#>   6   6   3   3 237  80   4
#>   7   1   7   1 115   8  12

Visualisation

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

PlotReduction(
    sce,
    reduction = vis_reduction,
    colour_by = "treatments",
    point_size = 0.5
)

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:

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:

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.

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)
#> [1] "unintegrated_by_batch" "integrated_by_batch"   "unintegrated_by_label"
#> [4] "integrated_by_label"   "integrated_split"

overview$unintegrated_by_batch

overview$integrated_by_batch

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.

bench_methods <- c("PCA", integrated_names)
set.seed(101)
bench <- BenchmarkIntegration(
    sce,
    batch_col = "treatments",
    label_col = "clusters",
    methods = bench_methods
)
bench$ranking
#>               method     score rank
#> 1     integrated_mnn 0.5678011    1
#> 2 integrated_harmony 0.5489021    2
#> 3                PCA 0.5221084    3
#> 4 integrated_centred 0.5182816    4

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.

head(bench$scores, 8)
#>    method reduction             metric       value     score higher_better
#> 1     PCA       PCA          pcr_batch  0.94364443 0.9436444          TRUE
#> 2     PCA       PCA              ilisi  1.16230360 0.0811518          TRUE
#> 3     PCA       PCA              clisi  1.65773927 0.6790128          TRUE
#> 4     PCA       PCA          kbet_like  0.41498002 0.4149800          TRUE
#> 5     PCA       PCA graph_connectivity  0.98154067 0.9815407          TRUE
#> 6     PCA       PCA                ari  0.21684337 0.2168434          TRUE
#> 11    PCA       PCA                nmi  0.37719573 0.3771957          TRUE
#> 12    PCA       PCA         silhouette -0.03500392 0.4824980          TRUE

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.

bench$plots$tradeoff

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 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.

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)
#>            
#>              1  2  3
#>   PMB0.5_1  84 83 83
#>   PMB0.5_4  84 83 83
#>   PMB0.5_7  84 83 83
#>   PMB2_1    84 83 83
#>   PMB2_4    39 39 39
#>   PMB2_7    84 83 83
#>   control_0 84 83 83
#>   control_1 84 83 83
#>   control_4 77 77 76
#>   control_7 82 81 81

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.

pb <- aggregate_pseudobulk(
    sce,
    sample_cols = c(
        "pb_id", "sample", "replicate", "treatments", "timepoints"
    ),
    aggregation = "sum",
    min_cells = 10
)
dim(pb)
#> [1] 3722   30
head(as.data.frame(colData(pb))[, c(
    "pb_id", "treatments", "timepoints", "ncells"
)], 4)
#>             pb_id treatments timepoints ncells
#> pb_1 control_0_r1    control          0     84
#> pb_2 control_0_r2    control          0     83
#> pb_3 control_0_r3    control          0     83
#> pb_4 control_1_r1    control          1     84

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.

pb <- normalize_pseudobulk(pb, method = "TMM")
#> calcNormFactors has been renamed to normLibSizes
summary(pb$norm.factors)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>  0.7796  0.8554  0.9546  1.0278  1.0745  1.8018

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.

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)

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.

pairwise_time <- run_edger_pairwise_de(
    pb,
    group_col = "timepoints",
    covariates = "treatments",
    contrast_type = "reference",
    reference_level = "0",
    normalization = "TMM",
    fdr_cutoff = 0.05
)
#> calcNormFactors has been renamed to normLibSizes
names(pairwise_time$tables)
#> [1] "1_vs_0" "4_vs_0" "7_vs_0"

$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.

table(
    pairwise_time$combined_table$contrast,
    pairwise_time$combined_table$significant
)
#>         
#>          FALSE TRUE
#>   1_vs_0   221  269
#>   4_vs_0   228  262
#>   7_vs_0   207  283

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)
#>                      gene     logFC          FDR
#> ABUW-RS11870 ABUW-RS11870  8.963791 1.106667e-17
#> omp33-36         omp33-36 -5.059676 8.840756e-13
#> ABUW-RS10220 ABUW-RS10220  4.550533 1.015487e-12
#> ABUW-RS03355 ABUW-RS03355 -7.426328 2.643033e-12
#> ABUW-RS17440 ABUW-RS17440 -4.146133 2.643033e-12
#> ABUW-RS19720 ABUW-RS19720  3.309762 7.440514e-11
#> ABUW-RS10635 ABUW-RS10635  6.364087 1.672197e-10
#> ABUW-RS03350 ABUW-RS03350 -6.809939 7.212523e-10
#> carO                 carO  3.584620 8.334559e-10
#> basD                 basD  4.442057 8.334559e-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.

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:

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:

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"
)
#> calcNormFactors has been renamed to normLibSizes
names(pairwise_dose$tables)
#> [1] "PMB2__7_vs_control__7"   "PMB0.5__7_vs_control__7"

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.

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
)
#> calcNormFactors has been renamed to normLibSizes
#> calcNormFactors has been renamed to normLibSizes
#> calcNormFactors has been renamed to normLibSizes
names(spline_de$by_condition)
#> [1] "control" "PMB0.5"  "PMB2"

table(
    spline_de$combined_table$treatments,
    spline_de$combined_table$significant
)
#>          
#>           FALSE TRUE
#>   control    98   91
#>   PMB0.5    284  389
#>   PMB2       11  140

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.

spline_top <- spline_de$combined_table[
    order(spline_de$combined_table$FDR),
    c("gene", "treatments", "FDR")
]
head(spline_top, 10)
#>                             gene treatments          FDR
#> PMB2.ABUW-RS17530   ABUW-RS17530       PMB2 5.572203e-37
#> PMB2.ssrS                   ssrS       PMB2 1.128559e-24
#> PMB2.ABUW-RS03195   ABUW-RS03195       PMB2 8.840879e-24
#> PMB2.ABUW-RS04445   ABUW-RS04445       PMB2 1.061374e-22
#> PMB2.ABUW-RS13935   ABUW-RS13935       PMB2 7.486060e-22
#> PMB0.5.ABUW-RS14490 ABUW-RS14490     PMB0.5 1.686219e-21
#> PMB2.secY                   secY       PMB2 4.676025e-21
#> PMB0.5.basJ                 basJ     PMB0.5 8.695325e-21
#> PMB0.5.katE                 katE     PMB0.5 9.209130e-21
#> PMB0.5.entE                 entE     PMB0.5 2.509724e-20

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

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 scater for plotting, scran for marker analysis, and other tools that accept this container.

Session information

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats4    stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#>  [1] ggplot2_4.0.3               SingleCellExperiment_1.35.2
#>  [3] SummarizedExperiment_1.43.0 Biobase_2.73.2             
#>  [5] GenomicRanges_1.65.4        Seqinfo_1.3.2              
#>  [7] IRanges_2.47.5              S4Vectors_0.51.10          
#>  [9] BiocGenerics_0.59.12        generics_0.1.4             
#> [11] MatrixGenerics_1.25.0       matrixStats_1.5.0          
#> [13] scProkaR_0.99.2             BiocStyle_2.41.0           
#> 
#> loaded via a namespace (and not attached):
#>  [1] rlang_1.3.0               magrittr_2.0.5           
#>  [3] RcppAnnoy_0.0.23          otel_0.2.0               
#>  [5] compiler_4.6.1            DelayedMatrixStats_1.35.0
#>  [7] vctrs_0.7.3               pkgconfig_2.0.3          
#>  [9] fastmap_1.2.0             XVector_0.53.0           
#> [11] scuttle_1.23.2            labeling_0.4.3           
#> [13] rmarkdown_2.32            xfun_0.61                
#> [15] bluster_1.23.1            cachem_1.1.0             
#> [17] beachmat_2.29.3           jsonlite_2.0.0           
#> [19] DelayedArray_0.39.8       BiocParallel_1.47.0      
#> [21] irlba_2.4.1               parallel_4.6.1           
#> [23] cluster_2.1.8.3           R6_2.6.1                 
#> [25] bslib_0.12.0              RColorBrewer_1.1-3       
#> [27] limma_3.99.1              parallelly_1.48.0        
#> [29] jquerylib_0.1.4           Rcpp_1.1.2               
#> [31] knitr_1.52                future.apply_1.20.2      
#> [33] splines_4.6.1             batchelor_1.29.0         
#> [35] Matrix_1.7-6              igraph_2.3.4             
#> [37] tidyselect_1.2.1          abind_1.4-8              
#> [39] yaml_2.3.12               codetools_0.2-20         
#> [41] listenv_1.1.0             lattice_0.23-1           
#> [43] tibble_3.3.1              withr_3.0.3              
#> [45] S7_0.2.2                  evaluate_1.0.5           
#> [47] future_1.76.0             pillar_1.11.1            
#> [49] BiocManager_1.30.27       sp_2.2-3                 
#> [51] sparseMatrixStats_1.25.0  scales_1.4.0             
#> [53] scrapper_1.7.9            globals_0.19.1           
#> [55] RhpcBLASctl_0.23-42       glue_1.8.1               
#> [57] maketools_1.3.2           tools_4.6.1              
#> [59] BiocNeighbors_2.7.3       sys_3.4.3                
#> [61] ScaledMatrix_1.21.0       RSpectra_0.16-2          
#> [63] locfit_1.5-9.12           RANN_2.6.3               
#> [65] buildtools_1.0.0          dotCall64_1.2            
#> [67] cowplot_1.2.0             grid_4.6.1               
#> [69] edgeR_4.99.6              BiocSingular_1.29.2      
#> [71] cli_3.6.6                 rsvd_1.0.5               
#> [73] spam_2.11-4               S4Arrays_1.13.2          
#> [75] dplyr_1.2.1               uwot_0.2.5               
#> [77] ResidualMatrix_1.23.0     gtable_0.3.6             
#> [79] sass_0.4.10               digest_0.6.39            
#> [81] progressr_1.0.0           SparseArray_1.13.4       
#> [83] SeuratObject_5.4.0        farver_2.1.2             
#> [85] htmltools_0.5.9           lifecycle_1.0.5          
#> [87] statmod_1.5.2             harmony_2.0.5

scProkaR core workflow for microbial single-cell RNA-seq