---
title: "Auditing donor influence in pseudobulk analysis"
author: "Fuhao Jiang"
date: "2026-10-06"
output:
  BiocStyle::html_document:
    toc: true
    toc_depth: 2
vignette: >
  %\VignetteIndexEntry{Auditing donor influence in pseudobulk analysis}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4)
```

## Introduction

scDonorAudit reports how a **fixed two-condition pseudobulk analysis** changes
when one biological donor is omitted. It uses
[edgeR](https://bioconductor.org/packages/edgeR) for the specified contrast and
links the planned omissions to design viability, execution events, and observed
gene-effect changes.

### Related work and scope

The following tools and analyses address different parts of this workflow:

| Tool or work | Primary task | Relationship to this package |
|:--|:--|:--|
| [scuttle](https://bioconductor.org/packages/scuttle) | General cell aggregation with `aggregateAcrossCells()` | `preparePseudobulk()` records registered sample/cell-type coverage, including combinations with no observed counts. |
| [muscat](https://bioconductor.org/packages/muscat) and [scran `pseudoBulkDGE()`](https://bioconductor.org/packages/scran) | Multi-sample or cell-type-wise differential expression | scDonorAudit specifies a two-condition independent or paired edgeR analysis and records changes across planned whole-donor omissions. |
| [dreamlet](https://bioconductor.org/packages/dreamlet) | Pseudobulk analysis for complex designs and weights, with assay- and gene-level errors | scDonorAudit links membership, design diagnosis, warning/error events, per-sample normalization, and per-gene validity to each planned omission. |
| [scFastDE](https://bioconductor.org/packages/devel/bioc/html/scFastDE.html) | Donor-weighted pseudobulk differential expression, including paired designs | Weighting and deletion sensitivity answer different questions. |
| [SampleQC](https://github.com/wmacnair/SampleQC) | Multi-sample single-cell quality control | Quality metrics and changes to a specified DE contrast are distinct outputs. |
| [Study-specific LODO pseudobulk analysis](https://pmc.ncbi.nlm.nih.gov/articles/PMC12991174/) | Refit after excluding each donor | scDonorAudit supplies the deletion ledger and summaries of observed change and valid coverage. |

The source repository is
[Emr-27/scDonorAudit](https://github.com/Emr-27/scDonorAudit).

## Installation

For the current source version, follow the installation instructions in the
[repository README](https://github.com/Emr-27/scDonorAudit#installation).

After acceptance and a successful Bioconductor build, use the command below
in the matching R/Bioconductor environment. The package will first be available
in devel and later in the applicable release; see the
[Bioconductor installation guide](https://bioconductor.org/install/).

```{r bioconductor-installation, eval=FALSE}
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install("scDonorAudit")
library(scDonorAudit)
```

## Aggregate single-cell counts

This small example shows how cell counts become pseudobulk counts and how
the sample registry records missing combinations. Two cells from `s1` are
summed into its `T`-cell column; registered sample `s3` has no observed cells.

```{r single-cell-aggregation}
sce_counts <- matrix(c(1L, 4L, 2L, 5L, 3L, 6L, 4L, 7L), nrow = 2L,
    dimnames = list(c("g1", "g2"), paste0("cell", 1:4)))
sce_cells <- S4Vectors::DataFrame(
    sample_id = c("s1", "s1", "s2", "s2"),
    cell_type = c("T", "T", "T", "B"))
sce_registry <- data.frame(sample_id = c("s1", "s2", "s3"),
    donor_id = c("d1", "d2", "d3"),
    condition = c("ctrl", "stim", "ctrl"))
sce_input <- SingleCellExperiment::SingleCellExperiment(
    assays = list(counts = sce_counts), colData = sce_cells)
cell_pb <- scDonorAudit::preparePseudobulk(sce_input, "sample_id", "donor_id",
    "condition", "cell_type", sample_table = sce_registry)
SummarizedExperiment::assay(cell_pb, "counts")
knitr::kable(as.data.frame(SummarizedExperiment::colData(cell_pb))[
    , c("sample_id", "cell_type", "n_cells")])
knitr::kable(S4Vectors::metadata(cell_pb)$scdonoraudit$coverage)
```

The three pseudobulk columns correspond to the observed sample/cell-type
combinations. The coverage table has six rows: registered samples crossed
with `B` and `T`. Missing combinations have `observed = FALSE` and
`n_cells = NA`, with no extra count columns.

## Construct an auditable example

The simulated counts below demonstrate registered coverage and paired-design
checks. Each column of the `SummarizedExperiment` is one observed
sample and cell type. Donor `d4` lacks a `B`-cell pseudobulk under `stim`.
The full sample registry makes that absence visible.

```{r input}
set.seed(1103)
sample_table <- data.frame(sample_id = paste0("s", seq_len(8)), donor_id = rep(paste0("d",
    seq_len(4)), each = 2), condition = rep(c("ctrl", "stim"), 4))
source_rows <- rbind(transform(sample_table, cell_type = "T"), transform(sample_table[seq_len(7),
    ], cell_type = "B"))
counts <- matrix(rnbinom(250 * nrow(source_rows), mu = 25, size = 5), nrow = 250,
    dimnames = list(paste0("g", seq_len(250)), paste0("column", seq_len(nrow(source_rows)))))
source_rows$n_cells <- 30L
se <- SummarizedExperiment::SummarizedExperiment(assays = list(counts = counts),
    colData = S4Vectors::DataFrame(source_rows))
pb <- scDonorAudit::preparePseudobulk(se, "sample_id", "donor_id", "condition",
    "cell_type", sample_table = sample_table, n_cells = "n_cells")
contrast <- c(numerator = "stim", denominator = "ctrl")
```

The named contrast fixes the sign: positive log2 fold-change means higher
expression in `stim` than `ctrl`. The count assay must contain raw,
non-negative integer counts. The package checks numeric validity; users must
confirm that the assay contains unnormalized counts. Factor coding is fixed inside
the design, so a user's global `options("contrasts")` does not change the
named effect's direction or scale.

Supported assays are in-memory numeric base matrices or numeric `Matrix`
objects. Logical/pattern matrices and `DelayedArray` or disk-backed assays
are unsupported. Sparse cell counts remain sparse during aggregation; only
the smaller cell-type pseudobulk matrices are converted for edgeR fitting.
For inputs that cannot fit in memory, aggregate raw counts upstream and use
the preaggregated `SummarizedExperiment` entry point, retaining the full sample
registry and `n_cells`. Output memory also grows with genes and planned fits.

## Review the design before fitting

For either pairing policy, the registry must contain one sample per condition
for every donor. `complete_pairs` handles incomplete cell-type eligibility;
the registry itself must remain paired.

The default paired policy is strict. An incomplete eligible pair blocks that
cell type without affecting other cell types.

```{r audit}
audit <- scDonorAudit::auditDesign(pb, "paired", contrast)
as.data.frame(audit[["spec"]])
subset(as.data.frame(audit[["sample_membership"]]), cell_type == "B" &
    reason_code != "OK")
```

`plan_status = "not_created"` means a deletion plan was never established for
the blocked cell type. It does not mean that all donor deletions succeeded or
that the result was stable.

## Fit and inspect the ledger

```{r fit}
result <- scDonorAudit::assessDonorInfluence(pb, "paired", contrast)
subset(as.data.frame(result[["fits"]]), select = c(cell_type, run_id, omitted_donor,
    execution_status, reason_code))
head(as.data.frame(scDonorAudit::summarizeInfluence(result, effect_threshold = 0.5))[,
    c("cell_type", "gene_id", "n_planned", "n_effect_valid", "max_abs_delta_observed")])
```

The summary's `n_planned` includes all eligible donor deletions, even failed
or skipped refits. Read the maximum effect change alongside valid-effect
coverage: it is the maximum among refits with comparable effects. Missing
effects remain missing, rather than becoming zero changes or
non-significant tests. Adjusted p values use a fixed baseline gene family
within each cell type; they do not control error across cell types or deletion
runs.

The analysis first fixes the eligible samples and complete-pair policy, then
plans all whole-donor deletions. Genes retained by the full-data filter remain
the comparison family in each refit. The named contrast and covariate meanings
also stay fixed. Each successful deletion refit **re-estimates** library-size
normalization and the edgeR model on its remaining samples.

Read each gene in this order: is its baseline effect valid; how many deletions
were planned; how many effects are comparable; then how large were the observed
changes? `n_sign_reversal` counts even tiny crossings of zero. With a positive
`effect_threshold`, `n_material_reversal` requires both effects to cross that
threshold in opposite directions, and `baseline_near_zero` flags a baseline
for which a raw sign change is especially easy to overread. The minimum and
maximum deletion effects are an **observed range**, not a confidence interval.
An unavailable effect may reflect design failure, backend failure, or a
gene-level status; it is not evidence of stability.

Effect changes describe sensitivity to donor omission. A large change can
reflect donor heterogeneity, limited precision, or other data properties; it
does not establish that a donor is invalid or should be excluded. The p values
come from the specified edgeR model, not a test that a donor is aberrant.

### Output dictionary

| Component | Meaning |
|:--|:--|
| `spec` | Per-cell-type plan status and number of planned donor omissions. |
| `sample_membership` / `coverage` | Registered sample-by-cell-type combinations and the actual baseline inclusion decisions. |
| `fits` / `fit_samples` | Baseline and deletion ledger. `design_reason_code` preserves the original design diagnosis; `execution_status`, `reason_code`, and `stage` describe the later execution outcome. The separate `fit_samples` table aligns sample IDs with library sizes and normalization factors for each completed run. |
| `gene_filter` / `results` | Fixed baseline gene family and gene-by-run effects, tests, validity flags, and gene status codes. |
| `issues` | Design problems and execution warnings or errors; a run can have both a backend warning and an error. |
| `config` / `provenance` | Requested design, contrast and thresholds; source preparation, package/R/edgeR versions, and the actual named contrast. |

The complete field/type/status dictionary is installed with the package:

```{r schema-location}
system.file("OUTPUT_SCHEMA.md", package = "scDonorAudit")
```

It also documents the 16 typed summary columns, which remain available in
zero-row summaries when all genes are filtered or all cell types are blocked.

Design rank deficiency can coexist with an estimable target contrast. A rank
problem is still a reason to skip fitting and is not relabeled as target
non-estimability. When a baseline cannot run, the planned deletion rows remain
in `fits`; `design_reason_code` keeps their independent audit results even if
their execution outcome becomes `BASELINE_UNAVAILABLE`.

To analyze the `B`-cell data from complete pairs, request the alternative
policy explicitly. This changes the analysis cohort, so report the excluded
pair when presenting the result.

```{r complete-pairs}
reduced <- scDonorAudit::assessDonorInfluence(pb, "paired", contrast, pair_policy = "complete_pairs")
subset(as.data.frame(reduced[["spec"]]), select = c(cell_type, plan_status,
    n_planned, n_samples))
```

## Inspect figures

The coverage plot shows included and excluded sample/cell-type combinations.
Gene plots use the explicit, stable gene IDs supplied by the user.

```{r coverage-plot}
scDonorAudit::plotInfluence(result, "coverage")
```

```{r gene-plots}
scDonorAudit::plotInfluence(result, "influence", cell_type = "T", gene_id = c("g1",
    "g2", "g3"))
scDonorAudit::plotInfluence(result, "effect", cell_type = "T", gene_id = "g1")
```

Before interpreting a gene, inspect `result[["fits"]]`,
`result[["sample_membership"]]`, and the number of valid deletion runs.
A failed fit or unestimable contrast limits the observed sensitivity analysis.
Grey influence tiles mean the effect is unavailable; consult the fit ledger
and gene status for the specific cause.

## Report missing coverage, design failure and valid effects together

This independent-donor example has eight donors, a missing `B` pseudobulk
for `d8`, and a batch category observed only in `d1`. The baseline models
are estimable. Omitting `d1` leaves the frozen batch column all zero, causing
rank deficiency; the target contrast can still be estimable. The failed
design is recorded as a skipped run. Other deletions provide valid effects.

```{r mixed-case}
set.seed(995201)
mixed_registry <- data.frame(sample_id = paste0("s", seq_len(8)), donor_id = paste0("d",
    seq_len(8)), condition = rep(c("ctrl", "stim"), each = 4), batch = c("special",
    rep("common", 7)))
mixed_rows <- rbind(transform(mixed_registry, cell_type = "T"), transform(mixed_registry[seq_len(7),
    ], cell_type = "B"))
mixed_rows$n_cells <- 30L
mixed_counts <- matrix(rnbinom(250 * nrow(mixed_rows), mu = 25, size = 5),
    nrow = 250, dimnames = list(paste0("g", seq_len(250)), NULL))
mixed_source <- SummarizedExperiment::SummarizedExperiment(assays = list(counts = mixed_counts),
    colData = S4Vectors::DataFrame(mixed_rows))
mixed_pb <- scDonorAudit::preparePseudobulk(mixed_source, "sample_id",
    "donor_id", "condition", "cell_type", sample_vars = "batch", sample_table = mixed_registry,
    n_cells = "n_cells")
mixed_result <- scDonorAudit::assessDonorInfluence(mixed_pb, "independent",
    contrast, covariates = "batch")
mixed_fits <- as.data.frame(mixed_result[["fits"]])
stopifnot(all(mixed_fits$execution_status[mixed_fits$omitted_donor %in%
    "d1"] == "skipped"), all(mixed_fits$reason_code[mixed_fits$omitted_donor %in%
    "d1"] == "DESIGN_RANK_DEFICIENT"))
knitr::kable(subset(as.data.frame(mixed_result[["sample_membership"]]),
    !in_baseline))
knitr::kable(mixed_fits[, c("cell_type", "run_id", "omitted_donor", "execution_status",
    "design_reason_code", "target_estimable")])
mixed_summary <- as.data.frame(scDonorAudit::summarizeInfluence(mixed_result))
knitr::kable(subset(mixed_summary, gene_id == "g1", select = c(cell_type,
    gene_id, n_planned, n_effect_valid, complete_effect_coverage, max_abs_delta_observed)),
    digits = 3)
scDonorAudit::plotInfluence(mixed_result, "influence", cell_type = "T",
    gene_id = c("g1", "g2", "g3"))
```

For `T`, eight omissions are planned but only seven can supply effects;
for `B`, seven are planned and six can supply effects. The missing `d8/B`
combination was never eligible for that cell type's deletion plan. Omitting
`d1` remains in each plan and appears as unavailable in the plot. The observed
maximum covers the valid deletions; the unavailable deletion has no measured
effect change. Backend failures use
`execution_status="failed"`, `reason_code="BACKEND_ERROR"` and an error event;
the same planned-denominator and missing-effect rules apply.

## Reproduce a real single-nucleus example

The package includes a compact pseudobulk derived from the public
[Crowell19 mouse cortex dataset](https://doi.org/10.6084/m9.figshare.8976473.v1)
(ExperimentHub `EH3297`). Raw counts from 1,740 Astrocyte nuclei were summed
within each of eight mice, retaining all 11,076 genes. Four mice received LPS
and four received Vehicle. The original single-nucleus object is available
through ExperimentHub and is not bundled here. The included
`extdata/README-crowell19.md` records the transformation, attribution, and
source's CC BY 4.0 license.

```{r real-data}
count_file <- system.file("extdata", "crowell19_astrocyte_counts.csv",
    package = "scDonorAudit")
sample_file <- system.file("extdata", "crowell19_astrocyte_samples.csv",
    package = "scDonorAudit")
real_counts <- as.matrix(utils::read.csv(count_file, row.names = 1, check.names = FALSE))
real_samples <- utils::read.csv(sample_file)
real_source <- SummarizedExperiment::SummarizedExperiment(assays = list(counts = real_counts),
    colData = S4Vectors::DataFrame(real_samples))
real_registry <- real_samples[, c("sample_id", "donor_id", "condition")]
real_pb <- scDonorAudit::preparePseudobulk(real_source, "sample_id", "donor_id",
    "condition", "cell_type", sample_table = real_registry, n_cells = "n_cells")
real_contrast <- c(numerator = "LPS", denominator = "Vehicle")
real_audit <- scDonorAudit::auditDesign(real_pb, "independent", real_contrast,
    min_cells = 20L)
as.data.frame(real_audit[["spec"]])
knitr::kable(as.data.frame(real_audit[["sample_membership"]])[, c("sample_id",
    "donor_id", "condition", "n_cells", "in_baseline")])
real_result <- scDonorAudit::assessDonorInfluence(real_pb, "independent",
    real_contrast, min_cells = 20L)
table(as.data.frame(real_result[["fits"]])$execution_status)
knitr::kable(as.data.frame(real_result[["fits"]])[, c("run_id", "omitted_donor",
    "execution_status", "design_reason_code", "reason_code")])
head(as.data.frame(scDonorAudit::summarizeInfluence(real_result))[, c("gene_id",
    "n_planned", "n_effect_valid", "max_abs_delta_observed")])
```

Eight mice offer limited replication, especially after one is omitted. Read
the observed changes with this sample size and the specified model in mind.

## Read the donor-deletion results

The fit ledger above has one full-data fit and eight donor-deletion fits. To
illustrate what the summary adds, consider genes with valid effects in all
eight deletions and an absolute full-data log2 fold-change of at least 1.
Order these genes by their largest observed absolute change, breaking ties by
gene ID, and show the genes at the 10th and 90th percentile ranks. This rule
selects examples with different sensitivity for interpreting the audit output.

```{r real-interpretation-selection}
real_summary <- as.data.frame(scDonorAudit::summarizeInfluence(real_result))
candidates <- real_summary[real_summary$complete_effect_coverage & is.finite(real_summary$baseline_logFC) &
    abs(real_summary$baseline_logFC) >= 1 & is.finite(real_summary$max_abs_delta_observed),
    ]
candidates <- candidates[order(candidates$max_abs_delta_observed, candidates$gene_id),
    ]
if (nrow(candidates) < 2L) {
    stop("The Crowell19 example requires at least two genes to illustrate ",
        "different levels of donor-deletion sensitivity; the stated ",
        "selection rule returned ", nrow(candidates), ".")
}
chosen <- candidates[c(ceiling(0.1 * nrow(candidates)), ceiling(0.9 * nrow(candidates))),
    ]
chosen$illustration <- c("lower change", "higher change")
knitr::kable(chosen[, c("illustration", "gene_id", "baseline_logFC", "n_planned",
    "n_effect_valid", "max_abs_delta_observed")], digits = 3)
```

For the higher-change example, inspect every deletion alongside the maximum.
A missing or failed deletion would appear in the fit ledger and reduce
`n_effect_valid`.

```{r real-interpretation-refits}
real_se <- real_result[["results"]][[chosen$cell_type[2L]]]
gene_row <- match(chosen$gene_id[2L], rownames(real_se))
refits <- as.data.frame(SummarizedExperiment::colData(real_se))
refits$logFC <- as.numeric(SummarizedExperiment::assay(real_se, "logFC")[gene_row,
    ])
refits$change_from_full <- refits$logFC - refits$logFC[1L]
largest <- which.max(abs(refits$change_from_full[-1L])) + 1L
knitr::kable(refits[, c("run_id", "omitted_donor", "logFC", "change_from_full")],
    digits = 3)
```

In this dataset, `r nrow(candidates)` genes meet the stated selection rule.
For `r chosen$gene_id[2L]`, the full-data LPS-versus-Vehicle log2 fold-change
is `r sprintf("%.3f", chosen$baseline_logFC[2L])`; omitting
`r refits$omitted_donor[largest]` changes it by
`r sprintf("%.3f", refits$change_from_full[largest])` to
`r sprintf("%.3f", refits$logFC[largest])`. Of the
`r chosen$n_planned[2L]` planned deletion effects,
`r chosen$n_effect_valid[2L]` are valid here. By comparison, the
lower-change example
`r chosen$gene_id[1L]` has a maximum absolute change of
`r sprintf("%.3f", chosen$max_abs_delta_observed[1L])`. These observations
describe sensitivity to donor omission for the specified model, these eight
mice, and the selected genes. The two genes illustrate the stated percentile
ranks; conclusions about other genes or datasets require their own analysis.

## Method and data references

- Chen Y, Chen L, Lun ATL, Baldoni PL, Smyth GK (2025). edgeR v4: powerful
  differential analysis of sequencing data with expanded functionality and
  improved support for small counts and larger datasets.
  *Nucleic Acids Research* 53(2), gkaf018.
  [doi:10.1093/nar/gkaf018](https://doi.org/10.1093/nar/gkaf018).
- Crowell HL et al. (2020). muscat detects subpopulation-specific state
  transitions from multi-sample multi-condition single-cell transcriptomics
  data. *Nature Communications* 11, 6077.
  [doi:10.1038/s41467-020-19894-4](https://doi.org/10.1038/s41467-020-19894-4).
  Dataset attribution and transformations are in `extdata/README-crowell19.md`.

Use `citation("edgeR")` for citations matching the installed backend version
and `citation("scDonorAudit")` for the package citation. The help-example data
are synthetic; their generation
is documented in `extdata/README-example.md`.

## Session information

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