Immune Profiling of COVID-19 Severity: A Cross-Sectional Analysis#

Dataset: Stephenson et al., Nature Medicine 2021 (E-MTAB-10026)

Background#

Stephenson et al. profiled peripheral blood mononuclear cells (PBMCs) from COVID-19 patients across severity groups and timepoints, providing insights into immune dynamics during disease progression.

Important Methodological Notes#

This is an OBSERVATIONAL study, not a randomized trial:

  • Disease severity (Mild vs Severe) is an outcome, not a treatment assignment

  • Patients are not randomly assigned to severity groups

  • Differences between groups reflect disease biology, not causal treatment effects

  • We use sctrial’s infrastructure for structured comparisons, but interpret results as descriptive associations

Time axis considerations:

  • Days From Onset (DFO): Biological time since symptom onset - captures disease progression

  • Collection_Day: Calendar time since study enrollment - captures sampling schedule

  • These are fundamentally different; we use DFO for biological interpretability

Analysis strategy: Given limited longitudinal pairing in this dataset, we focus on cross-sectional comparisons between severity groups at each timepoint, which is the most statistically appropriate approach.

1. Setup and Configuration#

[2]:

# Imports - consolidated in one cell
import warnings
warnings.filterwarnings('ignore', category=FutureWarning)
# Note: We do NOT suppress UserWarning — sctrial issues important
# statistical caveats (e.g. low-cluster reliability) as UserWarnings.

import numpy as np
import pandas as pd
import scipy.sparse as sp
import matplotlib.pyplot as plt
import seaborn as sns
import scanpy as sc
import statsmodels.formula.api as smf
from statsmodels.stats.multitest import multipletests
from scipy.stats import mannwhitneyu  # Import here for use across notebook
from pathlib import Path
import urllib.request

import sctrial as st
import itertools

# Configuration constants
MIN_GENES_FOR_SCORE = 5
MIN_PARTICIPANTS_FOR_COMPARISON = 5
MIN_CELLS_PER_CELLTYPE = 100
SEED = 42
FDR_ALPHA = 0.25  # Exploratory threshold; use 0.05 for confirmatory analyses

pd.options.mode.chained_assignment = None
print(f"sctrial version: {st.__version__ if hasattr(st, '__version__') else 'dev'}")


def _fmt_fdr(v):
    """Format FDR/p-value: scientific notation for very small values."""
    return f"{v:.2e}" if v < 0.001 else f"{v:.3f}"

sctrial version: 0.3.3

2. Data Loading and Processing#

[3]:

# Dataset loaders and helpers
# All available at top-level: st.load_stephenson_data, st.count_paired, st.categorize_celltype
from sctrial.datasets import load_stephenson_data, count_paired, categorize_celltype

[4]:

# Load data
adata = load_stephenson_data(force_reprocess=False)

print("\n=== Dataset Summary ===")
print(f"Cells: {adata.n_obs:,}")
print(f"Genes: {adata.n_vars:,}")
print(f"Severity groups: {adata.obs['severity'].unique().tolist()}")
print(f"DFO bins: {sorted(adata.obs['dfo_bin'].unique().tolist())}")
print(f"Cell types: {adata.obs['celltype'].nunique()}")
print(f"Participants: {adata.obs['participant_id'].nunique()}")


=== Dataset Summary ===
Cells: 205,202
Genes: 24,634
Severity groups: ['Severe', 'Mild']
DFO bins: ['DFO_0-7', 'DFO_15+', 'DFO_8-14']
Cell types: 50
Participants: 34

3. Sample Size Assessment#

Before analysis, we assess sample sizes to determine appropriate statistical methods.

[5]:

# Sample sizes by severity and timepoint
sample_sizes = (
    adata.obs
    .groupby(["severity", "dfo_bin"], observed=True)["participant_id"]
    .nunique()
    .unstack(fill_value=0)
)
print("Participants per severity × DFO bin:")
display(sample_sizes)

# Check for longitudinal pairing

dfo_visits = ["DFO_0-7", "DFO_8-14"]
n_paired_dfo = count_paired(adata.obs, "dfo_bin", dfo_visits)

print("")
print(f"Longitudinal pairing (DFO 0-7 → 8-14): {n_paired_dfo} participants")

if n_paired_dfo < MIN_PARTICIPANTS_FOR_COMPARISON:
    print("")
    print(f"⚠️  WARNING: Only {n_paired_dfo} paired participants available.")
    print("   Difference-in-Differences analysis requires ≥5 paired participants.")
    print("   This notebook will focus on CROSS-SECTIONAL comparisons instead.")
    CAN_DO_DID = False
else:
    CAN_DO_DID = True

# Visualize sample sizes
fig, axes = plt.subplots(1, 3, figsize=(14, 4))

# Cells by severity
adata.obs["severity"].value_counts().plot(
    kind="bar", ax=axes[0], color=["steelblue", "coral"]
)
axes[0].set_title("Cells by Severity")
axes[0].set_ylabel("Number of cells")

# Cells by DFO
adata.obs["dfo_bin"].value_counts().sort_index().plot(
    kind="bar", ax=axes[1], color="teal"
)
axes[1].set_title("Cells by Days From Onset")

# Participants per group
sample_sizes.T.plot(kind="bar", ax=axes[2])
axes[2].set_title("Participants per Group")
axes[2].set_ylabel("Number of participants")
axes[2].legend(title="Severity")

plt.tight_layout()
plt.show()

Participants per severity × DFO bin:
dfo_bin DFO_0-7 DFO_8-14 DFO_15+
severity
Mild 8 11 4
Severe 2 5 4 

Longitudinal pairing (DFO 0-7 → 8-14): 0 participants

⚠️  WARNING: Only 0 paired participants available.
   Difference-in-Differences analysis requires ≥5 paired participants.
   This notebook will focus on CROSS-SECTIONAL comparisons instead.
../_images/tutorials_example_covid19_stephenson_8_3.png

4. Study Design Configuration#

We configure the comparison structure. Note: We use sctrial’s TrialDesign for convenience, but this is an observational comparison, not a randomized trial.

[6]:

# Add log1p-CPM normalization
if "log1p_cpm" not in adata.layers:
    adata = st.add_log1p_cpm_layer(adata, counts_layer="counts", out_layer="log1p_cpm")
    print("Added log1p_cpm layer")

# Define comparison design
# Note: We call Severe the "treated" group for sctrial compatibility,
# but this is purely for software purposes - there is no treatment!
design = st.TrialDesign(
    participant_col="participant_id",
    visit_col="dfo_bin",
    arm_col="severity",
    arm_treated="Severe",   # Reference group for comparison
    arm_control="Mild",
    celltype_col="celltype",
)

# Available visits for analysis
available_visits = sorted(adata.obs["dfo_bin"].unique().tolist())
print(f"Available DFO bins: {available_visits}")

# Covariates - match actual column names in the dataset
covariate_mapping = {
    "Age_interval": "age",
    "Sex": "sex",
    "Site": "site",
    "Smoker": "smoker",
}
covariates = [c for c in covariate_mapping.keys() if c in adata.obs.columns]
print(f"Available covariates: {covariates}")

# Check covariate distributions
if covariates:
    print("\nCovariate summary by severity:")
    for cov in covariates[:3]:  # Show first 3
        print(f"\n{cov}:")
        print(adata.obs.groupby("severity")[cov].value_counts().unstack(fill_value=0))

design

Available DFO bins: ['DFO_0-7', 'DFO_15+', 'DFO_8-14']
Available covariates: ['Age_interval', 'Sex', 'Site', 'Smoker']

Covariate summary by severity:

Age_interval:
Age_interval  (20, 29]  (30, 39]  (40, 49]  (50, 59]  (60, 69]  (70, 79]  \
severity
Mild              7446     16252     22525     21651     35854     17798
Severe               0      5589     15433     32828      4164      2233

Age_interval  (80, 89]
severity
Mild             14410
Severe            9019

Sex:
Sex       Female   Male
severity
Mild       84262  51674
Severe     49556  19710

Site:
Site      Cambridge    Ncl  Sanger
severity
Mild          26277  87676   21983
Severe         6397  44275   18594

[6]:

TrialDesign(participant_col='participant_id', visit_col='dfo_bin', arm_col='severity', arm_treated='Severe', arm_control='Mild', celltype_col='celltype', crossover_col=None, baseline_visit=None, followup_visit=None)

4.1 Baseline Covariate Balance (Treated vs Control)#

Here we assess baseline covariate balance between severity groups at the earliest visit. Standardized mean differences (SMD) below 0.1 indicate good balance.

[7]:

# Baseline covariate balance (select visit with most balanced counts)
covariates = [c for c in ["Age_interval", "Sex", "Site"] if c in adata.obs.columns]
if covariates:
    if "dfo_bin" in adata.obs.columns:
        # Pick visit with max min(n_treated, n_control) at participant-level
        visit_scores = []
        for v in sorted(adata.obs["dfo_bin"].dropna().unique()):
            sub = adata.obs[adata.obs["dfo_bin"] == v]
            if design.participant_col in sub.columns:
                counts = sub.groupby(design.arm_col)[design.participant_col].nunique()
            else:
                counts = sub.groupby(design.arm_col).size()
            n_t = int(counts.get(design.arm_treated, 0))
            n_c = int(counts.get(design.arm_control, 0))
            visit_scores.append((v, n_t, n_c, min(n_t, n_c)))
        # choose visit with largest min(n_t, n_c); tie-breaker by total participants
        visit_scores.sort(key=lambda x: (x[3], x[1]+x[2]), reverse=True)
        visit0, n_t, n_c, _ = visit_scores[0]
        print(f"Selected baseline visit: {visit0} (treated={n_t}, control={n_c})")

        try:
            balance = st.check_covariate_balance(adata, design, covariates, visit=visit0)
            if not balance.empty:
                # Numeric covariates: one figure
                num_bal = balance[balance["level"].isna()].copy()
                if not num_bal.empty:
                    plt.figure(figsize=(6, max(2, 0.4 * len(num_bal))))
                    sns.barplot(data=num_bal, y="covariate", x="smd", color="steelblue")
                    plt.axvline(0.1, color="red", linestyle="--", linewidth=1)
                    plt.axvline(-0.1, color="red", linestyle="--", linewidth=1)
                    plt.title(f"Numeric covariate balance at {visit0}")
                    plt.xlabel("Standardized Mean Difference (SMD)")
                    plt.ylabel("Covariate")
                    plt.tight_layout()
                    plt.show()

                # Categorical covariates: one panel per covariate
                cat_bal = balance[balance["level"].notna()].copy()
                for cov in cat_bal["covariate"].unique():
                    sub = cat_bal[cat_bal["covariate"] == cov].copy()
                    sub = sub.reindex(sub["smd"].abs().sort_values(ascending=False).index)
                    plt.figure(figsize=(6, max(2, 0.35 * len(sub))))
                    sns.barplot(data=sub, y="level", x="smd", color="orchid")
                    plt.axvline(0.1, color="red", linestyle="--", linewidth=1)
                    plt.axvline(-0.1, color="red", linestyle="--", linewidth=1)
                    plt.title(f"Categorical balance: {cov} at {visit0}")
                    plt.xlabel("Standardized Mean Difference (SMD)")
                    plt.ylabel("Level")
                    plt.tight_layout()
                    plt.show()
        except Exception as e:
            print(f"Covariate balance check failed: {e}")
    else:
        print("dfo_bin not found; skipping covariate balance check.")
else:
    print("No baseline covariates available for balance checking.")

Selected baseline visit: DFO_8-14 (treated=5, control=11)
../_images/tutorials_example_covid19_stephenson_12_1.png
../_images/tutorials_example_covid19_stephenson_12_2.png
../_images/tutorials_example_covid19_stephenson_12_3.png 
[8]:

# Covariates for adjusted analyses (baseline balance indicated imbalance)
covariates_adj = [c for c in ["Age_interval", "Sex", "Site"] if c in adata.obs.columns]
print("Adjusted covariates:", covariates_adj)

Adjusted covariates: ['Age_interval', 'Sex', 'Site']

5. Cell Type Quality Control#

Filter to well-represented cell types for robust analysis.

[9]:

# Cell type counts
ct_counts = adata.obs["celltype"].value_counts()
print(f"Total cell types: {len(ct_counts)}")

# Filter to well-represented types
valid_celltypes = ct_counts[ct_counts >= MIN_CELLS_PER_CELLTYPE].index.tolist()
print(f"Cell types with ≥{MIN_CELLS_PER_CELLTYPE} cells: {len(valid_celltypes)}")

# Show top cell types
print("\nTop 15 cell types:")
display(ct_counts.head(15).to_frame("n_cells"))

# Visualize
fig, ax = plt.subplots(figsize=(12, 5))
ct_counts.head(20).plot(kind="bar", ax=ax, color="steelblue")
ax.axhline(MIN_CELLS_PER_CELLTYPE, color="red", linestyle="--", label=f"Min threshold ({MIN_CELLS_PER_CELLTYPE})")
ax.set_title("Cell Type Abundance")
ax.set_ylabel("Number of cells")
ax.legend()
plt.xticks(rotation=45, ha="right")
plt.tight_layout()
plt.show()

Total cell types: 50
Cell types with ≥100 cells: 41

Top 15 cell types:
n_cells
celltype
NK_16hi 29417
CD14_mono 20207
CD4.Naive 19076
CD4.CM 18804
B_naive 17063
CD83_CD14_mono 16326
CD8.TE 14366
CD8.Naive 10469
CD4.Tfh 5196
gdT 4938
CD8.EM 4864
Platelets 4824
CD4.IL22 4817
CD16_mono 4790
MAIT 3475
../_images/tutorials_example_covid19_stephenson_15_2.png

6. UMAP Visualization#

[10]:

# Compute UMAP if missing
if "X_umap" not in adata.obsm:
    print("Computing UMAP...")
    sc.pp.pca(adata)
    sc.pp.neighbors(adata)
    sc.tl.umap(adata)

# Global UMAP colored by cell type
fig, axes = plt.subplots(1, 2, figsize=(16, 6))

sc.pl.umap(adata, color="celltype", ax=axes[0], show=False,
           legend_loc="right margin", legend_fontsize=7, title="Cell Types")
sc.pl.umap(adata, color="severity", ax=axes[1], show=False,
           palette=["steelblue", "coral"], title="Severity")

plt.tight_layout()
plt.show()

# UMAP by severity and DFO
fig, axes = plt.subplots(2, 3, figsize=(15, 10))
for i, severity in enumerate(["Mild", "Severe"]):
    for j, dfo in enumerate(["DFO_0-7", "DFO_8-14", "DFO_15+"]):
        ax = axes[i, j]
        mask = (adata.obs["severity"] == severity) & (adata.obs["dfo_bin"] == dfo)
        if mask.sum() > 0:
            sub = adata[mask].copy()
            sc.pl.umap(sub, color="celltype", ax=ax, show=False,
                      legend_loc="none", title=f"{severity} - {dfo} (n={mask.sum():,})")
        else:
            ax.set_title(f"{severity} - {dfo} (no cells)")
            ax.axis("off")

plt.tight_layout()
plt.show()
../_images/tutorials_example_covid19_stephenson_17_0.png
../_images/tutorials_example_covid19_stephenson_17_1.png

7. COVID-19 Immune Signatures#

We define the 19 canonical gene signatures (matching the shared signature library used in the manuscript figures) using st.score_gene_sets with the z-mean method. Key signatures for COVID-19 profiling include:

  • Type I / Type II Interferon: Key antiviral pathways, often elevated in severe disease

  • Inflammatory Response: Cytokines and S100 alarmins associated with cytokine storm

  • Cytotoxic T Cell Activity: NK/CD8 T cell effector function, important for viral clearance

  • T Cell Exhaustion: Dysfunction markers elevated in severe/prolonged infection

  • Monocyte Macrophage: Monocyte/macrophage markers, drivers of immunopathology

  • Humoral Plasma Cell: B cell and plasma cell responses

  • (plus 13 additional canonical signatures)

[11]:

available_genes = set(adata.var_names)

# Canonical 19-signature dictionary (matches manuscript figures _shared.py)
gene_signatures = {
    "Cytotoxic T Cell Activity": [
        "GZMB", "GZMA", "GZMH", "GZMK", "GZMM", "PRF1", "GNLY", "NKG7",
        "KLRK1", "KLRD1", "KLRG1", "FASLG", "IFNG", "CST7", "CCL5",
        "CX3CR1", "FGFBP2",
    ],
    "T Cell Exhaustion": [
        "TOX", "TOX2", "PDCD1", "HAVCR2", "LAG3", "TIGIT", "ENTPD1",
        "CXCL13", "LAYN", "CD38", "PRDM1", "BATF", "NR4A2", "NR4A3",
        "CD160", "CD244",
    ],
    "Checkpoint Receptors": [
        "PDCD1", "CTLA4", "LAG3", "HAVCR2", "TIGIT", "BTLA", "CD160",
        "CD244", "CD96",
    ],
    "Type I Interferon": [
        "ISG15", "IFI27", "IFI44", "IFI44L", "IFI35", "MX1", "OAS1",
        "OASL", "IFIT2", "IFIT3", "IFITM3", "RSAD2", "LY6E", "USP18",
        "IRF7", "STAT2", "ISG20", "BST2", "HERC6", "EPSTI1", "SAMD9L",
        "PARP9", "RTP4", "LGALS3BP", "CMPK2",
    ],
    "Type II Interferon": [
        "STAT1", "IRF1", "GBP1", "GBP2", "GBP4", "GBP5", "CXCL9",
        "CXCL10", "CXCL11", "IDO1", "SOCS1", "SOCS3", "JAK2", "UBE2L6",
        "WARS", "IL18BP", "VCAM1", "SERPING1", "FGL2", "IFI30",
    ],
    "Memory T Cell": [
        "IL7R", "TCF7", "LEF1", "CCR7", "SELL", "CD27", "CD28", "ID3",
        "REL", "FOXP1", "BACH2", "MYB", "FOXO1", "KLF2", "S1PR1",
    ],
    "T Cell Activation": [
        "CD69", "CD44", "IL2RA", "ICOS", "TNFRSF4", "TNFRSF9", "CD40LG",
        "EGR1", "EGR2", "EGR3", "NR4A1", "MYC", "IL2RB", "IRF4",
        "TNFRSF18",
    ],
    "Inflammatory Response": [
        "IL1B", "IL6", "TNF", "CXCL8", "CCL2", "CCL3", "CCL4", "NFKB1",
        "NLRP3", "CASP1", "CXCL1", "CXCL2", "CXCL3", "CCL20", "PTGS2",
        "IL1A", "TNFAIP3", "NFKBIA", "IL18", "S100A8", "S100A9",
    ],
    "Antigen Presentation": [
        "HLA-A", "HLA-B", "HLA-C", "B2M", "TAP1", "TAP2", "TAPBP",
        "NLRC5", "PSMB8", "PSMB9", "PSMB10", "CIITA", "CD74", "HLA-DRA",
        "HLA-DRB1", "HLA-DPA1", "HLA-DPB1", "HLA-DQA1", "HLA-DMA",
        "HLA-DMB", "CTSS",
    ],
    "Cell Proliferation": [
        "MKI67", "TOP2A", "PCNA", "CDK1", "CCNB1", "CCNA2", "MCM2",
        "MCM7", "TYMS", "BIRC5", "UBE2C", "CENPF", "CENPE", "CCNB2",
        "RRM2", "HMGB2", "TPX2", "NUSAP1", "CKS2", "AURKA", "CKS1B",
        "STMN1",
    ],
    "Regulatory T Cell": [
        "FOXP3", "IKZF2", "IKZF4", "IL2RA", "CTLA4", "TNFRSF18",
        "TNFRSF4", "TNFRSF9", "TIGIT", "ENTPD1", "LRRC32", "RTKN2",
        "CCR8", "LAYN", "IL1R2",
    ],
    "NK Cell Activity": [
        "FCGR3A", "NCR1", "NCR3", "NCR2", "KLRF1", "KLRC1", "KLRC2",
        "KLRD1", "KLRB1", "IL2RB", "TBX21", "GZMM", "SH2D1B", "GNLY",
        "NKG7",
    ],
    "Apoptosis": [
        "BAX", "BAK1", "BID", "BCL2L11", "BBC3", "PMAIP1", "CASP3",
        "CASP7", "CASP8", "CASP9", "FAS", "TNFRSF10B", "DIABLO", "APAF1",
        "CDKN1A", "GADD45A", "GADD45B", "DDIT3", "FDXR", "BTG2",
    ],
    "Oxidative Stress Response": [
        "NFE2L2", "HMOX1", "NQO1", "GCLC", "GCLM", "GSR", "SOD1", "SOD2",
        "CAT", "GPX1", "TXN", "TXNRD1", "PRDX1", "PRDX2", "SRXN1",
        "G6PD", "GPX4", "GLRX", "SLC7A11",
    ],
    "Humoral Plasma Cell": [
        "MS4A1", "CD79A", "CD79B", "CD19", "CD37", "BANK1", "TCL1A",
        "FCRL5", "POU2AF1", "MZB1", "JCHAIN", "XBP1", "PRDM1", "DERL3",
        "SDC1", "TNFRSF17", "SEC11C", "HERPUD1", "FKBP11", "SLAMF7",
        "IGHM", "IGHG1", "IGKC", "IGLC2", "FDCSP",
    ],
    "Monocyte Macrophage": [
        "CD14", "CD68", "LYZ", "CSF1R", "AIF1", "TYROBP", "FCER1G",
        "CTSS", "FCGR3A", "ITGAM", "MNDA", "CD163", "MRC1", "MSR1",
        "MARCO", "APOE", "C1QA", "C1QB", "C1QC", "TREM2", "SELENOP",
        "LGALS3", "FCN1", "VCAN", "S100A9",
    ],
    "Tissue-Resident Memory": [
        "ITGAE", "ITGA1", "CXCR6", "ZNF683", "CD101", "RGS1", "RGS2",
        "DUSP6", "RBPJ", "CD69", "ID2", "PRDM1", "NR4A1", "DUSP1",
        "EGR1",
    ],
    "T Follicular Helper": [
        "CXCL13", "CXCR5", "PDCD1", "ICOS", "BCL6", "TOX", "TOX2",
        "IL21", "BTLA", "CD200", "MAF", "ASCL2", "SH2D1A", "CD40LG",
        "TNFRSF4", "BATF",
    ],
    "Hypoxia": [
        "VEGFA", "SLC2A1", "SLC2A3", "HK1", "HK2", "PGK1", "LDHA",
        "ALDOA", "ENO1", "PDK1", "BNIP3L", "NDRG1", "ADM", "ANGPTL4",
        "DDIT4", "PPP1R15A", "ERO1A", "P4HA1", "CA12", "IGFBP3", "PLIN2",
    ],
}

# Filter to available genes and report coverage
print("Gene signature coverage (19 canonical signatures):")
print("-" * 50)
filtered_signatures = {}
for name, genes in gene_signatures.items():
    found = [g for g in genes if g in available_genes]
    pct = len(found) / len(genes) * 100
    print(f"{name}: {len(found)}/{len(genes)} genes ({pct:.0f}%)")
    if len(found) >= MIN_GENES_FOR_SCORE:
        filtered_signatures[name] = found
    else:
        print(f"  Skipping (need >={MIN_GENES_FOR_SCORE} genes)")

# Score gene sets using z-mean method (accounts for different expression scales)
if filtered_signatures:
    adata = st.score_gene_sets(
        adata,
        filtered_signatures,
        layer="log1p_cpm",
        method="zmean",
        prefix="sig_"
    )
    print(f"\n Scored {len(filtered_signatures)} signatures")

# Get signature columns
signature_cols = [c for c in adata.obs.columns if c.startswith("sig_")]
print(f"\nSignature scores available: {signature_cols}")

Gene signature coverage (19 canonical signatures):
--------------------------------------------------
Cytotoxic T Cell Activity: 17/17 genes (100%)
T Cell Exhaustion: 16/16 genes (100%)
Checkpoint Receptors: 9/9 genes (100%)
Type I Interferon: 25/25 genes (100%)
Type II Interferon: 19/20 genes (95%)
Memory T Cell: 14/15 genes (93%)
T Cell Activation: 15/15 genes (100%)
Inflammatory Response: 21/21 genes (100%)
Antigen Presentation: 21/21 genes (100%)
Cell Proliferation: 22/22 genes (100%)
Regulatory T Cell: 15/15 genes (100%)
NK Cell Activity: 15/15 genes (100%)
Apoptosis: 20/20 genes (100%)
Oxidative Stress Response: 18/19 genes (95%)
Humoral Plasma Cell: 24/25 genes (96%)
Monocyte Macrophage: 25/25 genes (100%)
Tissue-Resident Memory: 15/15 genes (100%)
T Follicular Helper: 16/16 genes (100%)
Hypoxia: 21/21 genes (100%)

 Scored 19 signatures

Signature scores available: ['sig_Cytotoxic T Cell Activity', 'sig_T Cell Exhaustion', 'sig_Checkpoint Receptors', 'sig_Type I Interferon', 'sig_Type II Interferon', 'sig_Memory T Cell', 'sig_T Cell Activation', 'sig_Inflammatory Response', 'sig_Antigen Presentation', 'sig_Cell Proliferation', 'sig_Regulatory T Cell', 'sig_NK Cell Activity', 'sig_Apoptosis', 'sig_Oxidative Stress Response', 'sig_Humoral Plasma Cell', 'sig_Monocyte Macrophage', 'sig_Tissue-Resident Memory', 'sig_T Follicular Helper', 'sig_Hypoxia']

8. Cross-Sectional Comparisons: Severity Differences at Each Timepoint#

This is the primary analysis. We compare immune signatures between Mild and Severe patients at each DFO timepoint, using:

  • Participant-level aggregation (to avoid pseudoreplication)

  • Cell-type adjustment (to account for compositional differences)

  • FDR correction for multiple testing

[12]:

# Run comparisons at each timepoint using sctrial's built-in function
# This properly aggregates to participant level and handles the statistics correctly
print("="*60)
print("CROSS-SECTIONAL ANALYSIS: Severe vs Mild at each DFO bin")
print("="*60)
print("\nPositive beta = higher in Severe; Negative beta = higher in Mild")
print("Using participant-level aggregation to avoid pseudoreplication\n")

# Use a lower threshold (3) to include more timepoints for exploratory analysis
MIN_CROSS_SECTIONAL = 3

all_results = []
for visit in available_visits:
    print(f"\n{visit}:")

    # Check sample size first
    ad_visit = adata[adata.obs["dfo_bin"] == visit]
    n_per_group = ad_visit.obs.groupby("severity")["participant_id"].nunique()
    n_min = n_per_group.min()

    if n_min < MIN_CROSS_SECTIONAL:
        print(f"  Skipped: Insufficient participants (Mild={n_per_group.get('Mild', 0)}, Severe={n_per_group.get('Severe', 0)})")
        continue

    if n_min < MIN_PARTICIPANTS_FOR_COMPARISON:
        print(f"  Note: Small sample (Mild={n_per_group.get('Mild', 0)}, Severe={n_per_group.get('Severe', 0)}) — interpret with caution")

    # Use sctrial's built-in between_arm_comparison
    # This properly aggregates to participant level
    res = st.between_arm_comparison(
        adata,
        visit=visit,
        features=signature_cols,
        design=design,
        aggregate="participant_visit",
        standardize=True,
        method="ols",
    )

    if not res.empty:
        # Rename columns for consistency
        res = res.rename(columns={"beta_arm": "beta", "p_arm": "p_value", "FDR_arm": "fdr"})
        res["timepoint"] = visit
        all_results.append(res)

        # Display results
        display_cols = ["feature", "beta", "p_value", "fdr", "n_units"]
        display(res[display_cols].round(4))

        # Highlight significant
        sig = res[res["fdr"] < FDR_ALPHA]
        if not sig.empty:
            print(f"  Significant (FDR<{FDR_ALPHA}): {sig['feature'].tolist()}")

# Combine results
if all_results:
    combined_results = pd.concat(all_results, ignore_index=True)
else:
    combined_results = pd.DataFrame()

============================================================
CROSS-SECTIONAL ANALYSIS: Severe vs Mild at each DFO bin
============================================================

Positive beta = higher in Severe; Negative beta = higher in Mild
Using participant-level aggregation to avoid pseudoreplication


DFO_0-7:
  Skipped: Insufficient participants (Mild=8, Severe=2)

DFO_15+:
  Note: Small sample (Mild=4, Severe=4) — interpret with caution
feature beta p_value fdr n_units
0 sig_Cytotoxic T Cell Activity 0.7073 0.3557 0.9375 8
1 sig_T Cell Exhaustion 0.4132 0.5991 0.9375 8
2 sig_Checkpoint Receptors 0.0826 0.9173 0.9375 8
3 sig_Type I Interferon 0.4334 0.5810 0.9375 8
4 sig_Type II Interferon 0.6399 0.4069 0.9375 8
5 sig_Memory T Cell -0.9252 0.2128 0.9375 8
6 sig_T Cell Activation -0.0940 0.9059 0.9375 8
7 sig_Inflammatory Response 0.7047 0.3577 0.9375 8
8 sig_Antigen Presentation -0.5003 0.5220 0.9375 8
9 sig_Cell Proliferation 0.3280 0.6780 0.9375 8
10 sig_Regulatory T Cell 0.4100 0.6020 0.9375 8
11 sig_NK Cell Activity 0.6846 0.3727 0.9375 8
12 sig_Apoptosis -0.2463 0.7560 0.9375 8
13 sig_Oxidative Stress Response -0.0624 0.9375 0.9375 8
14 sig_Humoral Plasma Cell -0.3339 0.6724 0.9375 8
15 sig_Monocyte Macrophage 0.9365 0.2064 0.9375 8
16 sig_Tissue-Resident Memory -0.1212 0.8789 0.9375 8
17 sig_T Follicular Helper -0.1505 0.8498 0.9375 8
18 sig_Hypoxia -0.2393 0.7628 0.9375 8 

DFO_8-14:
feature beta p_value fdr n_units
0 sig_Cytotoxic T Cell Activity 0.1083 0.8488 0.9947 16
1 sig_T Cell Exhaustion 0.2592 0.6471 0.9947 16
2 sig_Checkpoint Receptors 0.2125 0.7077 0.9947 16
3 sig_Type I Interferon -0.0038 0.9947 0.9947 16
4 sig_Type II Interferon -0.3955 0.4825 0.9947 16
5 sig_Memory T Cell -0.3468 0.5389 0.9947 16
6 sig_T Cell Activation -0.1168 0.8370 0.9947 16
7 sig_Inflammatory Response -0.1321 0.8160 0.9947 16
8 sig_Antigen Presentation -0.1378 0.8082 0.9947 16
9 sig_Cell Proliferation 1.0609 0.0446 0.4237 16
10 sig_Regulatory T Cell 0.7941 0.1464 0.9271 16
11 sig_NK Cell Activity 0.0613 0.9141 0.9947 16
12 sig_Apoptosis -0.2132 0.7069 0.9947 16
13 sig_Oxidative Stress Response 0.6504 0.2405 0.9947 16
14 sig_Humoral Plasma Cell 1.2357 0.0158 0.3001 16
15 sig_Monocyte Macrophage -0.0170 0.9761 0.9947 16
16 sig_Tissue-Resident Memory -0.3109 0.5822 0.9947 16
17 sig_T Follicular Helper 0.3914 0.4871 0.9947 16
18 sig_Hypoxia 0.3308 0.5580 0.9947 16

8.1 Adjusted Cross-Sectional Comparisons#

[13]:

# Adjusted cross-sectional analysis (includes covariates)
print('\nAdjusted cross-sectional analysis (covariates):', covariates_adj)
cross_sectional_adj = []
features_adj = signature_cols
visits_adj = tuple(sorted(
    adata.obs[design.visit_col].dropna().unique(),
    key=lambda c: int(c.split('_')[1].split('-')[0].rstrip('+')),
)[:2])
if features_adj:
    for v in visits_adj:
        print(f'\nAdjusted results at visit: {v}')
        res = st.between_arm_comparison(
            adata,
            visit=v,
            features=features_adj,
            design=design,
            aggregate='participant_visit',
            standardize=True,
            method='ols',
            covariates=covariates_adj if covariates_adj else None,
        )
        if res is not None and not res.empty:
            res['visit'] = v
            cross_sectional_adj.append(res)
            display(res[['feature', 'beta_arm', 'p_arm', 'FDR_arm', 'n_units']].round(4))
if cross_sectional_adj:
    all_cross_adj = __import__('pandas').concat(cross_sectional_adj, ignore_index=True)
else:
    all_cross_adj = __import__('pandas').DataFrame()


Adjusted cross-sectional analysis (covariates): ['Age_interval', 'Sex', 'Site']

Adjusted results at visit: DFO_0-7
feature beta_arm p_arm FDR_arm n_units
0 sig_Cytotoxic T Cell Activity 1.6032 0.2774 0.5583 10
1 sig_T Cell Exhaustion 0.0658 0.9558 0.9942 10
2 sig_Checkpoint Receptors -0.0629 0.9665 0.9942 10
3 sig_Type I Interferon -0.0310 0.9529 0.9942 10
4 sig_Type II Interferon -1.0476 0.0435 0.1999 10
5 sig_Memory T Cell -0.7990 0.0668 0.2115 10
6 sig_T Cell Activation -0.8079 0.0923 0.2505 10
7 sig_Inflammatory Response -0.8251 0.4945 0.7829 10
8 sig_Antigen Presentation -1.0444 0.3007 0.5583 10
9 sig_Cell Proliferation 0.0112 0.9942 0.9942 10
10 sig_Regulatory T Cell -0.2294 0.6239 0.8214 10
11 sig_NK Cell Activity 0.7910 0.5856 0.8214 10
12 sig_Apoptosis -1.3513 0.0059 0.1122 10
13 sig_Oxidative Stress Response -1.4389 0.0365 0.1999 10
14 sig_Humoral Plasma Cell -0.3907 0.6485 0.8214 10
15 sig_Monocyte Macrophage -0.8579 0.2538 0.5583 10
16 sig_Tissue-Resident Memory -1.8448 0.0526 0.1999 10
17 sig_T Follicular Helper -0.7806 0.3233 0.5583 10
18 sig_Hypoxia -1.5210 0.0207 0.1969 10 

Adjusted results at visit: DFO_8-14
feature beta_arm p_arm FDR_arm n_units
0 sig_Cytotoxic T Cell Activity -0.2901 0.7362 0.9142 16
1 sig_T Cell Exhaustion 0.6987 0.2205 0.4656 16
2 sig_Checkpoint Receptors -0.0466 0.9075 0.9579 16
3 sig_Type I Interferon -0.2702 0.5857 0.8560 16
4 sig_Type II Interferon -0.3901 0.4237 0.7319 16
5 sig_Memory T Cell -0.0927 0.7926 0.9142 16
6 sig_T Cell Activation 0.3347 0.3558 0.6761 16
7 sig_Inflammatory Response 0.7218 0.1863 0.4425 16
8 sig_Antigen Presentation -0.3277 0.6905 0.9142 16
9 sig_Cell Proliferation 1.4313 0.0131 0.1480 16
10 sig_Regulatory T Cell 1.0010 0.0296 0.1480 16
11 sig_NK Cell Activity 0.1767 0.8179 0.9142 16
12 sig_Apoptosis 0.6596 0.0450 0.1480 16
13 sig_Oxidative Stress Response 0.8233 0.0337 0.1480 16
14 sig_Humoral Plasma Cell 1.3473 0.0313 0.1480 16
15 sig_Monocyte Macrophage 0.0157 0.9663 0.9663 16
16 sig_Tissue-Resident Memory 0.3744 0.4975 0.7876 16
17 sig_T Follicular Helper 0.2950 0.0688 0.1868 16
18 sig_Hypoxia 0.9004 0.0467 0.1480 16

9. Visualization of Severity Differences#

[14]:

# Visualize signature distributions by severity and timepoint
fig, axes = plt.subplots(2, 3, figsize=(15, 10))

for i, sig in enumerate(signature_cols[:6]):
    ax = axes.flat[i]

    # Aggregate to participant level for cleaner visualization
    df_plot = (
        adata.obs
        .groupby(["participant_id", "severity", "dfo_bin"], observed=True)[sig]
        .mean()
        .reset_index()
    )

    sns.boxplot(
        data=df_plot, x="dfo_bin", y=sig, hue="severity",
        palette={"Mild": "steelblue", "Severe": "coral"},
        ax=ax
    )
    ax.set_title(sig.replace("sig_", ""))
    ax.set_xlabel("Days From Onset")
    ax.set_ylabel("Score (z-mean)")
    if i > 0:
        ax.get_legend().remove()

plt.tight_layout()
plt.show()

# Heatmap of effect sizes across ALL visits (computed directly from participant means)
# Uses a lower threshold (≥2 per group) for visualization purposes
heatmap_rows = []
for visit in available_visits:
    ad_v = adata[adata.obs["dfo_bin"] == visit]
    df_agg = (
        ad_v.obs
        .groupby(["participant_id", "severity"], observed=True)[signature_cols]
        .mean()
        .reset_index()
    )
    n_mild = (df_agg["severity"] == "Mild").sum()
    n_severe = (df_agg["severity"] == "Severe").sum()
    if min(n_mild, n_severe) < 2:
        continue
    for sig in signature_cols:
        mild_vals = df_agg.loc[df_agg["severity"] == "Mild", sig].values
        severe_vals = df_agg.loc[df_agg["severity"] == "Severe", sig].values
        diff = severe_vals.mean() - mild_vals.mean()
        heatmap_rows.append({"feature": sig, "timepoint": visit, "diff": diff})

if heatmap_rows:
    df_heatmap = pd.DataFrame(heatmap_rows)
    pivot = df_heatmap.pivot(index="feature", columns="timepoint", values="diff")
    # Sort columns chronologically (alphabetical puts DFO_15-27 before DFO_8-14)
    pivot = pivot[sorted(pivot.columns, key=lambda c: int(c.split("_")[1].split("-")[0].rstrip("+")))]

    plt.figure(figsize=(8, 6))
    sns.heatmap(
        pivot, cmap="RdBu_r", center=0, annot=True, fmt=".2f",
        cbar_kws={"label": "Mean difference (Severe - Mild)"}
    )
    plt.title("Severity Differences Across Timepoints\n(participant-level mean difference)")
    plt.tight_layout()
    plt.show()
elif not combined_results.empty:
    pivot = combined_results.pivot(index="feature", columns="timepoint", values="beta")
    pivot = pivot[sorted(pivot.columns, key=lambda c: int(c.split("_")[1].split("-")[0].rstrip("+")))]
    plt.figure(figsize=(8, 6))
    sns.heatmap(
        pivot, cmap="RdBu_r", center=0, annot=True, fmt=".2f",
        cbar_kws={"label": "Effect size (Severe - Mild)"}
    )
    plt.title("Severity Effect Sizes Across Timepoints")
    plt.tight_layout()
    plt.show()
../_images/tutorials_example_covid19_stephenson_25_0.png
../_images/tutorials_example_covid19_stephenson_25_1.png

10. Cell-Type Specific Analysis#

Examine severity differences within specific immune cell populations.

[15]:

# Focus on major immune populations
# Map fine cell types to major lineages for robust analysis
adata.obs["lineage"] = adata.obs["celltype"].apply(categorize_celltype)
print("Lineage distribution:")
print(adata.obs["lineage"].value_counts())

focus_lineages = ["CD4_T", "CD8_T", "Monocytes", "NK"]

# Auto-select the visit with the best participant balance across arms
visit_scores = []
for v in sorted(adata.obs["dfo_bin"].dropna().unique()):
    sub = adata.obs[adata.obs["dfo_bin"] == v]
    counts = sub.groupby("severity")["participant_id"].nunique()
    n_mild = int(counts.get("Mild", 0))
    n_severe = int(counts.get("Severe", 0))
    visit_scores.append((v, n_mild, n_severe, min(n_mild, n_severe)))
visit_scores.sort(key=lambda x: (x[3], x[1] + x[2]), reverse=True)
focus_visit = visit_scores[0][0]

print(f"\n=== Cell-Type Specific Analysis at {focus_visit} ===")
print(f"    (auto-selected: best participant balance across arms)")

# Use sctrial between_arm_comparison per lineage
lineage_results = []
for lineage in focus_lineages:
    ad_sub = adata[(adata.obs["lineage"] == lineage) & (adata.obs["dfo_bin"] == focus_visit)].copy()
    if ad_sub.n_obs < 100:
        print(f"{lineage}: Insufficient cells ({ad_sub.n_obs})")
        continue
    n_per_group = ad_sub.obs.groupby("severity")["participant_id"].nunique()
    if n_per_group.min() < 3:
        print(f"{lineage}: Insufficient participants")
        continue
    print(f"{lineage} (n={ad_sub.n_obs:,} cells, {n_per_group.sum()} participants):")

    res = st.between_arm_comparison(
        ad_sub,
        visit=focus_visit,
        features=signature_cols[:4],
        design=design,
        method="wilcoxon",
    )
    res["lineage"] = lineage
    lineage_results.append(res)

if lineage_results:
    df_lineage = pd.concat(lineage_results, ignore_index=True)
    # Add FDR correction across all tests
    from statsmodels.stats.multitest import multipletests
    df_lineage["FDR"] = multipletests(df_lineage["p_arm"], method="fdr_bh")[1]
    print("\nResults (participant-level Mann-Whitney U via sctrial):")
    display(df_lineage[["lineage", "feature", "beta_arm", "p_arm", "FDR"]].round(4))

    # Pivot for heatmap
    pivot = df_lineage.pivot(index="feature", columns="lineage", values="beta_arm")
    plt.figure(figsize=(8, 5))
    sns.heatmap(pivot, cmap="RdBu_r", center=0, annot=True, fmt=".2f")
    plt.title(f"Severity Effect by Cell Type ({focus_visit})")
    plt.tight_layout()
    plt.show()

Lineage distribution:
lineage
CD4_T        49451
Monocytes    42563
NK           35643
CD8_T        33613
B_cells      28582
Other        11316
DCs           4034
Name: count, dtype: int64

=== Cell-Type Specific Analysis at DFO_8-14 ===
    (auto-selected: best participant balance across arms)
CD4_T (n=20,684 cells, 16 participants):
CD8_T (n=15,382 cells, 16 participants):
Monocytes (n=19,550 cells, 16 participants):
NK (n=14,538 cells, 16 participants):

Results (participant-level Mann-Whitney U via sctrial):
lineage feature beta_arm p_arm FDR
0 CD4_T sig_Cytotoxic T Cell Activity -0.0332 0.7427 1.0
1 CD4_T sig_T Cell Exhaustion 0.0153 0.4409 1.0
2 CD4_T sig_Checkpoint Receptors 0.0417 0.1804 1.0
3 CD4_T sig_Type I Interferon 0.0393 0.9130 1.0
4 CD8_T sig_Cytotoxic T Cell Activity 0.2510 0.3773 1.0
5 CD8_T sig_T Cell Exhaustion 0.0273 0.7427 1.0
6 CD8_T sig_Checkpoint Receptors 0.0516 0.4409 1.0
7 CD8_T sig_Type I Interferon 0.0324 0.8269 1.0
8 Monocytes sig_Cytotoxic T Cell Activity -0.0265 0.1804 1.0
9 Monocytes sig_T Cell Exhaustion -0.0340 0.6612 1.0
10 Monocytes sig_Checkpoint Receptors -0.0148 0.6612 1.0
11 Monocytes sig_Type I Interferon 0.0377 0.6612 1.0
12 NK sig_Cytotoxic T Cell Activity -0.0238 0.5833 1.0
13 NK sig_T Cell Exhaustion -0.0271 0.9130 1.0
14 NK sig_Checkpoint Receptors -0.0231 1.0000 1.0
15 NK sig_Type I Interferon 0.0001 1.0000 1.0
../_images/tutorials_example_covid19_stephenson_27_2.png

11. Cell Type Abundance Differences#

Compare cell type composition between severity groups.

[16]:

# Calculate cell type proportions per participant
# Note: Each row in 'props' after grouping is one participant-visit-lineage combination
# When we test, each participant contributes ONE proportion value per lineage
props = (
    adata.obs
    .groupby(["participant_id", "severity", "dfo_bin", "lineage"], observed=True)
    .size()
    .reset_index(name="n")
)
totals = props.groupby(["participant_id", "dfo_bin"])["n"].transform("sum")
props["proportion"] = props["n"] / totals

# Visualize
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# By severity
props_by_sev = props.groupby(["severity", "lineage"])["proportion"].mean().unstack()
props_by_sev.plot(kind="bar", ax=axes[0], width=0.8)
axes[0].set_title("Cell Type Proportions by Severity")
axes[0].set_ylabel("Mean Proportion")
axes[0].legend(title="Lineage", bbox_to_anchor=(1.02, 1))

# Monocyte proportion over time
mono_props = props[props["lineage"] == "Monocytes"].copy()
sns.boxplot(data=mono_props, x="dfo_bin", y="proportion", hue="severity",
            palette={"Mild": "steelblue", "Severe": "coral"}, ax=axes[1])
axes[1].set_title("Monocyte Proportion by Severity and Time")
axes[1].set_ylabel("Proportion")

plt.tight_layout()
plt.show()

# Statistical test for abundance differences
# CORRECT: Each observation in props_early is already at participant level
# so Mann-Whitney U is appropriate here (no pseudoreplication)
print("\nAbundance differences (Severe vs Mild) at DFO 0-7:")
print("Note: Each data point = one participant's cell type proportion")
print("-" * 50)

props_early = props[props["dfo_bin"] == "DFO_0-7"]
abundance_results = []
for lineage in ["Monocytes", "CD8_T", "NK", "B_cells"]:
    sub = props_early[props_early["lineage"] == lineage]
    mild = sub.loc[sub["severity"] == "Mild", "proportion"]
    severe = sub.loc[sub["severity"] == "Severe", "proportion"]

    if len(mild) >= 3 and len(severe) >= 3:
        stat, pval = mannwhitneyu(mild, severe, alternative="two-sided")
        diff = severe.mean() - mild.mean()
        abundance_results.append({
            "lineage": lineage,
            "diff": diff,
            "p_value": pval,
            "n_mild": len(mild),
            "n_severe": len(severe),
        })
        print(f"{lineage}: diff={diff:+.3f}, p={pval:.4f} (n_mild={len(mild)}, n_severe={len(severe)})")

if abundance_results:
    df_abund = pd.DataFrame(abundance_results)
    df_abund["fdr"] = multipletests(df_abund["p_value"], method="fdr_bh")[1]
    print(f"\nFDR-corrected results:")
    display(df_abund.round(4))
../_images/tutorials_example_covid19_stephenson_29_0.png 

Abundance differences (Severe vs Mild) at DFO 0-7:
Note: Each data point = one participant's cell type proportion
--------------------------------------------------

11.1 Adjusted Abundance DiD#

[17]:

# Adjusted abundance DiD (covariates)
print('\nAdjusted abundance DiD (covariates):', covariates_adj)
# Build numeric covariates for abundance models when possible
covariates_ab = []
if 'Age_interval' in adata.obs.columns:
    # Convert intervals like '(50, 59]' to midpoint
    def _age_mid(x):
        try:
            s = str(x).strip().strip('()[]')
            lo, hi = s.split(',')
            return (float(lo) + float(hi)) / 2
        except Exception:
            return float('nan')
    adata.obs['age_mid'] = adata.obs['Age_interval'].map(_age_mid)
    if adata.obs['age_mid'].notna().any():
        covariates_ab.append('age_mid')
if 'Sex' in adata.obs.columns:
    covariates_ab.append('Sex')
if 'Site' in adata.obs.columns:
    # high-cardinality sites can destabilize models; include only if <=3 levels
    if adata.obs['Site'].nunique() <= 3:
        covariates_ab.append('Site')
print('Adjusted abundance covariates:', covariates_ab)

if design.celltype_col:
    # Diagnostic: participant counts per celltype at selected visits
    visits_ab = sorted(
        adata.obs[design.visit_col].dropna().unique(),
        key=lambda c: int(c.split('_')[1].split('-')[0].rstrip('+')),
    )[:2]
    diag = (
        adata.obs[adata.obs[design.visit_col].isin(visits_ab)]
        .groupby([design.celltype_col, design.arm_col])[design.participant_col]
        .nunique()
        .unstack(fill_value=0)
    )
    print('Participant counts per cell type (treated/control):')
    display(diag)

    def _run_ab(covs, label):
        print(f'Attempt: {label} covariates -> {covs}')
        return st.abundance_did(
            adata,
            design,
            visits=tuple(visits_ab),
            covariates=covs if covs else None,
            min_units=3,
            use_bootstrap=True,
            n_boot=199,
            seed=SEED,
        )

    ab_adj = _run_ab(covariates_ab, 'full')
    if ab_adj is None or ab_adj.empty:
        # retry with numeric-only covariates
        num_covs = [c for c in covariates_ab if c == 'age_mid']
        ab_adj = _run_ab(num_covs, 'numeric-only')
    if ab_adj is None or ab_adj.empty:
        # final fallback: unadjusted
        ab_adj = _run_ab([], 'unadjusted fallback')

    if ab_adj is not None and not ab_adj.empty:
        display(ab_adj.round(4))
        try:
            plt.figure(figsize=(10, max(3, 0.4 * len(ab_adj))))
            sns.barplot(data=ab_adj, y='celltype', x='beta_DiD', color='teal')
            plt.axvline(0, color='black', linewidth=0.8)
            plt.title('Adjusted Abundance DiD (covariates)')
            plt.xlabel('DiD effect (beta)')
            plt.ylabel('Cell type')
            plt.tight_layout()
            plt.show()
        except Exception as e:
            print(f'Adjusted abundance plot failed: {e}')
    else:
        print('No adjusted abundance DiD results after retries.')


    # If DiD is not feasible (no paired participants), run cross-sectional abundance at baseline
    if ab_adj is None or ab_adj.empty:
        # Check paired participants overall
        visits_use = tuple(visits_ab)
        obs_use = adata.obs[adata.obs[design.visit_col].isin(visits_use)].copy()
        totals_all = (
            obs_use.groupby([design.participant_col, design.visit_col, design.arm_col], observed=True)
            .size()
            .reset_index(name='n_cells')
        )
        wide_tot = totals_all.pivot_table(
            index=design.participant_col,
            columns=design.visit_col,
            values='n_cells',
            aggfunc='mean',
            observed=True,
        )
        paired_units = wide_tot[wide_tot[visits_use[0]].notna() & wide_tot[visits_use[1]].notna()].index
        if len(paired_units) == 0:
            print('No paired participants for abundance DiD. Running baseline cross-sectional abundance instead.')
            baseline_visit = visits_use[0]
            obs_base = adata.obs[adata.obs[design.visit_col] == baseline_visit].copy()
            # proportions per participant × celltype
            counts_base = (
                obs_base.groupby([design.participant_col, design.arm_col, design.celltype_col], observed=True)
                .size()
                .reset_index(name='n_cells')
            )
            totals_base = (
                counts_base.groupby([design.participant_col, design.arm_col], observed=True)['n_cells']
                .sum()
                .reset_index(name='total_cells')
            )
            counts_base = counts_base.merge(totals_base, on=[design.participant_col, design.arm_col], how='left')
            counts_base['prop'] = counts_base['n_cells'] / counts_base['total_cells'].clip(lower=1)

            # attach covariates (participant-level)
            cov_df = obs_base[[design.participant_col] + covariates_ab].drop_duplicates()
            counts_base = counts_base.merge(cov_df, on=design.participant_col, how='left')
            counts_base['arm_bin'] = design.arm_bin(counts_base)
            counts_base = counts_base.dropna(subset=['arm_bin'])

            rows = []
            for ct in sorted(counts_base[design.celltype_col].unique()):
                tmp = counts_base[counts_base[design.celltype_col] == ct].copy()
                if tmp[design.participant_col].nunique() < 3:
                    continue
                if tmp['prop'].nunique() < 2:
                    continue
                formula = 'prop ~ arm_bin'
                if covariates_ab:
                    formula += ' + ' + ' + '.join(covariates_ab)
                import statsmodels.formula.api as smf
                fit = smf.ols(formula, data=tmp).fit()
                if 'arm_bin' not in fit.params or fit.params['arm_bin'] is None:
                    continue
                rows.append({
                    'celltype': ct,
                    'n_participants': int(tmp[design.participant_col].nunique()),
                    'beta_arm': float(fit.params['arm_bin']),
                    'p_arm': float(fit.pvalues['arm_bin']),
                })
            if rows:
                import pandas as pd
                from statsmodels.stats.multitest import multipletests
                res_xs = pd.DataFrame(rows)
                mask = res_xs['p_arm'].notna()
                res_xs['FDR_arm'] = pd.NA
                if mask.any():
                    res_xs.loc[mask, 'FDR_arm'] = multipletests(res_xs.loc[mask,'p_arm'], method='fdr_bh')[1]
                display(res_xs.sort_values('p_arm').head(20))
                try:
                    plt.figure(figsize=(10, max(3, 0.4 * len(res_xs))))
                    sns.barplot(data=res_xs, y='celltype', x='beta_arm', color='teal')
                    plt.axvline(0, color='black', linewidth=0.8)
                    plt.title(f'Baseline abundance differences at {baseline_visit} (covariate-adjusted)')
                    plt.xlabel('Arm effect (beta)')
                    plt.ylabel('Cell type')
                    plt.tight_layout()
                    plt.show()
                except Exception as e:
                    print(f'Baseline abundance plot failed: {e}')


Adjusted abundance DiD (covariates): ['Age_interval', 'Sex', 'Site']
Adjusted abundance covariates: ['age_mid', 'Sex', 'Site']
Participant counts per cell type (treated/control):
severity Mild Severe
celltype
ASDC 13 4
B_exhausted 19 7
B_immature 19 7
B_naive 19 7
B_non-switched_memory 19 7
B_switched_memory 19 7
C1_CD16_mono 11 5
CD4.CM 19 7
CD4.EM 19 7
CD4.IL22 19 7
CD4.Naive 19 7
CD4.Prolif 16 7
CD4.Tfh 19 7
CD4.Th1 16 4
CD4.Th2 5 3
CD4.Th17 2 2
CD8.EM 19 7
CD8.Naive 19 7
CD8.Prolif 18 6
CD8.TE 19 7
CD14_mono 18 7
CD16_mono 16 7
CD83_CD14_mono 19 7
DC1 15 4
DC2 18 7
DC3 19 7
DC_prolif 6 2
HSC_CD38neg 11 7
HSC_CD38pos 18 5
HSC_MK 4 4
HSC_erythroid 15 7
HSC_myeloid 6 3
HSC_prolif 7 2
ILC1_3 15 7
ILC2 7 2
MAIT 19 7
Mono_prolif 9 5
NKT 19 7
NK_16hi 19 7
NK_56hi 19 7
NK_prolif 19 7
Plasma_cell_IgA 18 7
Plasma_cell_IgG 19 7
Plasma_cell_IgM 17 7
Plasmablast 19 7
Platelets 19 7
RBC 13 7
Treg 11 3
gdT 19 7
pDC 18 7 
Attempt: full covariates -> ['age_mid', 'Sex', 'Site']
Attempt: numeric-only covariates -> ['age_mid']
Attempt: unadjusted fallback covariates -> []
No adjusted abundance DiD results after retries.
No paired participants for abundance DiD. Running baseline cross-sectional abundance instead.
celltype n_participants beta_arm p_arm FDR_arm
21 DC1 7 0.001090 0.032445 0.781814
11 CD4.IL22 10 -0.007421 0.044069 0.781814
36 NK_16hi 10 0.075020 0.055844 0.781814
23 DC3 10 -0.001242 0.143460 0.875708
35 NKT 10 0.002001 0.147275 0.875708
19 CD8.TE 10 0.130125 0.165944 0.875708
2 B_immature 10 -0.005815 0.167006 0.875708
1 B_exhausted 10 0.001560 0.173708 0.875708
28 HSC_erythroid 9 -0.000836 0.218194 0.875708
6 C1_CD16_mono 8 -0.010052 0.234704 0.875708
46 gdT 10 -0.033003 0.295915 0.875708
47 pDC 10 0.015888 0.312960 0.875708
8 CD16_mono 10 -0.034212 0.324744 0.875708
12 CD4.Naive 10 -0.034840 0.337689 0.875708
42 Plasmablast 10 0.003829 0.337901 0.875708
39 Plasma_cell_IgA 10 0.003787 0.340136 0.875708
3 B_naive 10 -0.025556 0.354453 0.875708
16 CD8.EM 10 -0.011124 0.413749 0.909467
40 Plasma_cell_IgG 10 0.003477 0.427501 0.909467
34 Mono_prolif 6 0.001219 0.434371 0.909467
../_images/tutorials_example_covid19_stephenson_31_4.png

12. Individual Gene Analysis#

Examine specific genes of interest for COVID-19.

[18]:

# Key COVID-19 genes
covid_genes = [
    "ISG15", "IFIT1", "MX1",       # IFN response
    "S100A8", "S100A9", "S100A12",  # Alarmins
    "GZMB", "PRF1", "NKG7",         # Cytotoxicity
    "IL1B", "TNF", "CXCL8",         # Cytokines
]
covid_genes = [g for g in covid_genes if g in adata.var_names]
print(f"Analyzing {len(covid_genes)} COVID-relevant genes")

# Extract expression for visualization
gene_expr = pd.DataFrame(
    adata[:, covid_genes].layers["log1p_cpm"].toarray() if sp.issparse(adata[:, covid_genes].layers["log1p_cpm"])
    else adata[:, covid_genes].layers["log1p_cpm"],
    columns=covid_genes,
    index=adata.obs_names
)
gene_expr = gene_expr.join(adata.obs[["severity", "dfo_bin", "participant_id", "lineage"]])

# Aggregate to participant level
gene_means = gene_expr.groupby(["participant_id", "severity", "dfo_bin"])[covid_genes].mean().reset_index()

# Visualize key genes
fig, axes = plt.subplots(2, 3, figsize=(14, 8))
for i, gene in enumerate(covid_genes[:6]):
    ax = axes.flat[i]
    sns.boxplot(
        data=gene_means, x="dfo_bin", y=gene, hue="severity",
        palette={"Mild": "steelblue", "Severe": "coral"}, ax=ax
    )
    ax.set_title(gene)
    ax.set_xlabel("")
    if i > 0:
        ax.get_legend().remove()

plt.tight_layout()
plt.show()

# Dotplot of genes by cell type
if len(covid_genes) > 0:
    adata_dot = adata[adata.obs["lineage"].isin(["Monocytes", "CD8_T", "NK", "CD4_T"])].copy()
    sc.pl.dotplot(
        adata_dot,
        var_names=covid_genes[:8],
        groupby="lineage",
        standard_scale="var",
        title="COVID-19 Genes by Cell Type"
    )

Analyzing 12 COVID-relevant genes
../_images/tutorials_example_covid19_stephenson_33_1.png
../_images/tutorials_example_covid19_stephenson_33_2.png

13. Advanced Statistical Analyses#

Statistical tools for cross-sectional observational data:

  • Effect sizes: Cohen’s d for severity group comparisons

  • Power analysis: Sample size planning for future studies

  • Effective sample size: Understanding clustering effects

13.1 Effect Sizes for Cross-Sectional Comparisons#

Effect sizes quantify the magnitude of differences between Mild and Severe cases.

[19]:

print("=" * 60)
print("EFFECT SIZE ANALYSIS (Cross-Sectional)")
print("=" * 60)

if signature_cols:
    # Calculate effect sizes for severity comparison at each timepoint
    effect_results = []

    for visit in available_visits:
        ad_visit = adata[adata.obs["dfo_bin"] == visit]

        # Aggregate to participant level
        df_agg = (
            ad_visit.obs
            .groupby(["participant_id", "severity"], observed=True)[signature_cols]
            .mean()
            .reset_index()
        )

        for sig in signature_cols:
            mild_vals = df_agg.loc[df_agg["severity"] == "Mild", sig].dropna().values
            severe_vals = df_agg.loc[df_agg["severity"] == "Severe", sig].dropna().values

            if len(mild_vals) >= 3 and len(severe_vals) >= 3:
                # Cohen's d and Hedge's g
                d = st.cohens_d(severe_vals, mild_vals)  # Positive = higher in Severe
                g = st.hedges_g(severe_vals, mild_vals)

                # Bootstrap CI (returns: effect_size, ci_lower, ci_upper)
                try:
                    _, ci_low, ci_high = st.bootstrap_effect_size_ci(
                        severe_vals, mild_vals,
                        n_boot=999,
                        alpha=0.05,
                        seed=SEED
                    )
                except Exception:
                    ci_low, ci_high = np.nan, np.nan

                effect_results.append({
                    "timepoint": visit,
                    "feature": sig,
                    "cohens_d": d,
                    "hedges_g": g,
                    "ci_lower": ci_low,
                    "ci_upper": ci_high,
                    "n_mild": len(mild_vals),
                    "n_severe": len(severe_vals),
                })

    if effect_results:
        df_effect = pd.DataFrame(effect_results)

        print("\nEffect sizes (Severe vs Mild):")
        print("  Positive = higher in Severe, Negative = higher in Mild")
        print("")
        display(df_effect.round(3))

        # Heatmap of effect sizes
        pivot = df_effect.pivot(index="feature", columns="timepoint", values="hedges_g")
        # Sort columns chronologically
        pivot = pivot[sorted(pivot.columns, key=lambda c: int(c.split("_")[1].split("-")[0].rstrip("+")))]

        plt.figure(figsize=(8, 6))
        sns.heatmap(
            pivot, cmap="RdBu_r", center=0, annot=True, fmt=".2f",
            cbar_kws={"label": "Hedge's g"}
        )
        plt.title("Effect Sizes: Severe vs Mild by Timepoint")
        plt.tight_layout()
        plt.show()

        # Forest plot for one timepoint
        visit_data = df_effect[df_effect["timepoint"] == "DFO_8-14"]
        if not visit_data.empty:
            fig, ax = plt.subplots(figsize=(8, 5))
            y_pos = np.arange(len(visit_data))

            ax.barh(y_pos, visit_data["hedges_g"], xerr=[
                visit_data["hedges_g"] - visit_data["ci_lower"],
                visit_data["ci_upper"] - visit_data["hedges_g"]
            ], capsize=5, color=["coral" if g > 0 else "steelblue" for g in visit_data["hedges_g"]])
            ax.axvline(0, color="black", linewidth=0.5)
            ax.set_yticks(y_pos)
            ax.set_yticklabels(visit_data["feature"])
            ax.set_xlabel("Hedge's g (95% CI)")
            ax.set_title("Effect Sizes at DFO 8-14 (Severe vs Mild)")

            # Reference lines
            for thresh in [-0.8, -0.5, -0.2, 0.2, 0.5, 0.8]:
                ax.axvline(thresh, color="gray", linestyle=":", alpha=0.5)

            plt.tight_layout()
            plt.show()
else:
    print("No signature columns for effect size analysis.")

============================================================
EFFECT SIZE ANALYSIS (Cross-Sectional)
============================================================

Effect sizes (Severe vs Mild):
  Positive = higher in Severe, Negative = higher in Mild
timepoint feature cohens_d hedges_g ci_lower ci_upper n_mild n_severe
0 DFO_15+ sig_Cytotoxic T Cell Activity 0.707 0.614 -1.688 6.606 4 4
1 DFO_15+ sig_T Cell Exhaustion 0.392 0.341 -1.053 2.652 4 4
2 DFO_15+ sig_Checkpoint Receptors 0.077 0.066 -2.077 1.728 4 4
3 DFO_15+ sig_Type I Interferon 0.412 0.358 -0.849 4.844 4 4
4 DFO_15+ sig_Type II Interferon 0.631 0.548 -0.715 3.076 4 4
5 DFO_15+ sig_Memory T Cell -0.986 -0.856 -5.441 0.473 4 4
6 DFO_15+ sig_T Cell Activation -0.087 -0.076 -3.461 4.153 4 4
7 DFO_15+ sig_Inflammatory Response 0.704 0.612 -0.583 2.820 4 4
8 DFO_15+ sig_Antigen Presentation -0.481 -0.418 -1.982 1.007 4 4
9 DFO_15+ sig_Cell Proliferation 0.308 0.268 -1.369 1.878 4 4
10 DFO_15+ sig_Regulatory T Cell 0.389 0.338 -0.912 3.330 4 4
11 DFO_15+ sig_NK Cell Activity 0.681 0.592 -1.521 8.783 4 4
12 DFO_15+ sig_Apoptosis -0.230 -0.200 -2.273 1.247 4 4
13 DFO_15+ sig_Oxidative Stress Response -0.058 -0.050 -1.523 2.062 4 4
14 DFO_15+ sig_Humoral Plasma Cell -0.314 -0.273 -2.832 0.973 4 4
15 DFO_15+ sig_Monocyte Macrophage 1.002 0.870 -0.233 10.230 4 4
16 DFO_15+ sig_Tissue-Resident Memory -0.112 -0.098 -2.338 3.204 4 4
17 DFO_15+ sig_T Follicular Helper -0.140 -0.121 -2.726 2.598 4 4
18 DFO_15+ sig_Hypoxia -0.223 -0.194 -2.013 2.287 4 4
19 DFO_8-14 sig_Cytotoxic T Cell Activity 0.105 0.099 -0.823 1.337 11 5
20 DFO_8-14 sig_T Cell Exhaustion 0.252 0.239 -0.971 1.476 11 5
21 DFO_8-14 sig_Checkpoint Receptors 0.206 0.195 -0.918 1.419 11 5
22 DFO_8-14 sig_Type I Interferon -0.004 -0.003 -0.683 1.337 11 5
23 DFO_8-14 sig_Type II Interferon -0.389 -0.368 -1.071 0.896 11 5
24 DFO_8-14 sig_Memory T Cell -0.340 -0.321 -1.391 0.540 11 5
25 DFO_8-14 sig_T Cell Activation -0.113 -0.107 -1.372 0.977 11 5
26 DFO_8-14 sig_Inflammatory Response -0.128 -0.121 -1.615 1.277 11 5
27 DFO_8-14 sig_Antigen Presentation -0.133 -0.126 -0.919 0.777 11 5
28 DFO_8-14 sig_Cell Proliferation 1.190 1.125 0.032 2.878 11 5
29 DFO_8-14 sig_Regulatory T Cell 0.829 0.784 -0.104 1.926 11 5
30 DFO_8-14 sig_NK Cell Activity 0.059 0.056 -1.002 1.352 11 5
31 DFO_8-14 sig_Apoptosis -0.207 -0.196 -1.658 0.851 11 5
32 DFO_8-14 sig_Oxidative Stress Response 0.661 0.625 -0.110 2.867 11 5
33 DFO_8-14 sig_Humoral Plasma Cell 1.481 1.400 0.465 3.630 11 5
34 DFO_8-14 sig_Monocyte Macrophage -0.016 -0.016 -0.810 0.974 11 5
35 DFO_8-14 sig_Tissue-Resident Memory -0.304 -0.287 -1.634 0.864 11 5
36 DFO_8-14 sig_T Follicular Helper 0.385 0.364 -0.683 1.598 11 5
37 DFO_8-14 sig_Hypoxia 0.324 0.306 -0.764 1.513 11 5
../_images/tutorials_example_covid19_stephenson_36_2.png
../_images/tutorials_example_covid19_stephenson_36_3.png

13.2 Power Analysis#

Power analysis helps understand what effect sizes this study could detect and plan future longitudinal studies.

Caveat: This dataset has no longitudinally paired participants for the primary Mild vs Severe comparison. The power curves below use the average group size as a rough guide for what a future paired longitudinal study of this scale could detect — they do not represent the power of the current cross-sectional design. Cross-sectional comparisons generally have lower power than paired longitudinal designs.

[20]:

print("=" * 60)
print("POWER ANALYSIS")
print("=" * 60)

# Current sample sizes by group
n_per_group = (
    adata.obs
    .groupby(["severity", "dfo_bin"], observed=True)["participant_id"]
    .nunique()
)
print("\nParticipants per group:")
display(n_per_group.unstack())

# Get typical sample sizes
n_mild = adata.obs[adata.obs["severity"] == "Mild"]["participant_id"].nunique()
n_severe = adata.obs[adata.obs["severity"] == "Severe"]["participant_id"].nunique()
n_avg = (n_mild + n_severe) / 2

print(f"\nOverall: {n_mild} Mild, {n_severe} Severe participants")
print(f"Average per group: {n_avg:.0f}")

# Power with current sample
print(f"\nHypothetical longitudinal DiD power (n≈{n_avg:.0f} per arm, IF paired):")
print("  Note: This dataset has no paired participants; values are for study planning only.")
for effect_size in [0.5, 0.8, 1.0, 1.5]:
    power = st.power_did(n_per_group=int(n_avg), effect_size=effect_size)
    print(f"  Effect size d={effect_size}: {power:.1%} power")

# Sample size needed for longitudinal DiD
print("\nSample size for longitudinal DiD (80% power):")
print("  (Future study planning)")
for effect_size in [0.5, 0.8, 1.0]:
    n_needed = st.sample_size_did(effect_size=effect_size, power=0.80)
    print(f"  Effect size d={effect_size}: {n_needed} per arm ({2*n_needed} total)")

# Power curve visualization
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# Power curves
n_range = np.arange(5, 51)
for effect_size, color in [(0.5, "blue"), (0.8, "green"), (1.0, "orange")]:
    powers = [st.power_did(n_per_group=n, effect_size=effect_size) for n in n_range]
    axes[0].plot(n_range, powers, label=f"d={effect_size}", color=color, linewidth=2)

axes[0].axhline(0.8, color="black", linestyle="--", alpha=0.5, label="80% power")
axes[0].axvline(n_avg, color="gray", linestyle=":", label=f"Current n≈{n_avg:.0f}")
axes[0].set_xlabel("Sample size per arm")
axes[0].set_ylabel("Power")
axes[0].set_title("Power Curves for Future Longitudinal Study")
axes[0].legend(loc="lower right")
axes[0].set_ylim(0, 1)
axes[0].grid(True, alpha=0.3)

# Effective sample size illustration
total_cells = adata.n_obs
cells_per_pt = adata.obs.groupby("participant_id").size().mean()
print(f"\nAverage cells per participant: {cells_per_pt:.0f}")

iccs = np.linspace(0.01, 0.30, 50)
design_effects = [st.design_effect(cells_per_pt, icc) for icc in iccs]
n_participants = adata.obs["participant_id"].nunique()
effective_ns = [st.effective_sample_size(n_participants, cells_per_pt, icc) for icc in iccs]

axes[1].plot(iccs, design_effects, linewidth=2, color="steelblue")
axes[1].set_xlabel("Intraclass Correlation (ICC)")
axes[1].set_ylabel("Design Effect")
axes[1].set_title(f"Design Effect vs ICC\n(avg {cells_per_pt:.0f} cells per participant)")
axes[1].grid(True, alpha=0.3)

# Add effective n on secondary axis
ax2 = axes[1].twinx()
ax2.plot(iccs, effective_ns, linewidth=2, color="coral", linestyle="--")
ax2.set_ylabel("Effective Sample Size", color="coral")
ax2.tick_params(axis='y', labelcolor='coral')

plt.tight_layout()
plt.show()

print("\n" + "=" * 60)
print("KEY INSIGHTS")
print("=" * 60)
print(f"""
1. Cross-sectional comparison:
   - Current sample provides good power for medium-large effects
   - Participant-level aggregation is critical for valid inference

2. Longitudinal limitation:
   - Only {n_paired_dfo} participants have paired observations
   - Future studies need ≥10 paired per arm for DiD

3. Cell-level clustering:
   - With ICC=0.10 and {cells_per_pt:.0f} cells/participant:
     Design effect = {st.design_effect(cells_per_pt, 0.10):.0f}
   - This is why participant-level analysis is essential
""")

============================================================
POWER ANALYSIS
============================================================

Participants per group:
dfo_bin DFO_0-7 DFO_8-14 DFO_15+
severity
Mild 8 11 4
Severe 2 5 4 

Overall: 23 Mild, 11 Severe participants
Average per group: 17

Hypothetical longitudinal DiD power (n≈17 per arm, IF paired):
  Note: This dataset has no paired participants; values are for study planning only.
  Effect size d=0.5: 30.8% power
  Effect size d=0.8: 64.5% power
  Effect size d=1.0: 83.0% power
  Effect size d=1.5: 99.2% power

Sample size for longitudinal DiD (80% power):
  (Future study planning)
  Effect size d=0.5: 63 per arm (126 total)
  Effect size d=0.8: 25 per arm (50 total)
  Effect size d=1.0: 16 per arm (32 total)

Average cells per participant: 6035
../_images/tutorials_example_covid19_stephenson_38_3.png 

============================================================
KEY INSIGHTS
============================================================

1. Cross-sectional comparison:
   - Current sample provides good power for medium-large effects
   - Participant-level aggregation is critical for valid inference

2. Longitudinal limitation:
   - Only 0 participants have paired observations
   - Future studies need ≥10 paired per arm for DiD

3. Cell-level clustering:
   - With ICC=0.10 and 6035 cells/participant:
     Design effect = 604
   - This is why participant-level analysis is essential

Single-Gene Cross-Sectional Check (Participant-Level)#

Here, we use a focused single-gene comparison to validate cross-sectional differences at the participant level. This avoids pseudoreplication and provides an interpretable effect for one gene in a specific lineage.

[21]:

print("\n=== Single-Gene Cross-Sectional Check (Participant-Level) ===")
gene = "ISG15"
celltypes = sorted(adata.obs["lineage"].dropna().unique())

# Compute participant-level means per celltype and group
gene_name = st.resolve_feature(adata, gene)
# Use log1p-CPM (library-size normalized) to avoid confounding by sequencing depth
expr = st.extract_gene_vector(adata, gene_name, layer="log1p_cpm")
df_expr = adata.obs[["participant_id", "severity", "lineage"]].copy()
df_expr["expr"] = expr

df_part = (
    df_expr.groupby(["participant_id", "severity", "lineage"], observed=True)["expr"]
    .mean()
    .reset_index()
)

# Run per-celltype tests
rows = []
for ct in celltypes:
    result, _ = st.compare_gene_in_celltype(
        adata,
        gene=gene,
        celltypes=ct,
        group_col="severity",
        group1="Severe",
        group2="Mild",
        participant_col="participant_id",
        celltype_col="lineage",
        min_cells_per_patient=10,
        min_patients_per_group=3,
    )
    rows.append(result)

df_res = pd.DataFrame(rows)
df_res["celltype"] = df_res["celltypes"].apply(lambda x: x[0] if isinstance(x, list) and len(x) == 1 else str(x))
display(df_res)

# Plot participant-level distributions by celltype
g = sns.catplot(
    data=df_part,
    x="lineage",
    y="expr",
    hue="severity",
    kind="box",
    height=4,
    aspect=2.2,
    palette={"Mild": "steelblue", "Severe": "coral"},
)
g.set_xticklabels(rotation=45, ha="right")
g.set_axis_labels("Cell type", f"{gene} (log1p CPM)")
g.fig.suptitle(f"{gene} by Cell Type and Severity (Participant-Level)", y=1.05)
plt.tight_layout()
plt.show()

# Barplot of effect sizes with -log10(p) overlay
df_plot = df_res.dropna(subset=["p_value"]).copy()
if not df_plot.empty:
    df_plot["neglog10p"] = -np.log10(df_plot["p_value"].clip(lower=1e-12))
    df_plot = df_plot.sort_values("delta")
    fig, ax1 = plt.subplots(figsize=(8, 4))
    ax1.bar(df_plot["celltype"], df_plot["delta"], color="gray")
    ax1.axhline(0, color="black", linewidth=0.8)
    ax1.set_ylabel("Mean(Severe) - Mean(Mild)")
    ax1.set_xlabel("Cell type")
    ax1.set_title(f"{gene} Cross-Sectional Effect by Cell Type")
    ax1.tick_params(axis="x", rotation=45)
    ax2 = ax1.twinx()
    ax2.plot(df_plot["celltype"], df_plot["neglog10p"], color="darkred", marker="o")
    ax2.set_ylabel("-log10(p)")
    plt.tight_layout()
    plt.show()
else:
    print("No valid p-values available for plotting.")


=== Single-Gene Cross-Sectional Check (Participant-Level) ===
gene celltypes group1 group2 n_group1 n_group2 mean_group1 mean_group2 delta p_value celltype
0 ISG15 [B_cells] Severe Mild 11 23 0.416558 0.465953 -0.049395 0.632249 B_cells
1 ISG15 [CD4_T] Severe Mild 11 23 0.349447 0.416083 -0.066635 0.912062 CD4_T
2 ISG15 [CD8_T] Severe Mild 11 23 0.406572 0.488814 -0.082242 0.941309 CD8_T
3 ISG15 [DCs] Severe Mild 10 21 0.699574 0.906295 -0.206721 0.688090 DCs
4 ISG15 [Monocytes] Severe Mild 11 23 0.780016 0.974678 -0.194662 0.768375 Monocytes
5 ISG15 [NK] Severe Mild 11 23 0.480815 0.571585 -0.090770 0.712779 NK
6 ISG15 [Other] Severe Mild 11 23 0.221410 0.340090 -0.118680 0.357408 Other
../_images/tutorials_example_covid19_stephenson_40_2.png
../_images/tutorials_example_covid19_stephenson_40_3.png

Module-Score Pseudobulk DiD by Lineage#

Here, we collapse each participant×visit×lineage into a pseudobulk mean of module scores, then compute DiD per lineage. This highlights which lineages drive the severity-associated changes.

Note: This dataset has no longitudinally paired participants, so the code falls back to a repeated-cross-sectional (unpaired) OLS DiD. This is a weaker design than a true paired DiD — treat any discoveries as exploratory hypothesis-generating only, not confirmatory.

[22]:

if signature_cols:
    # Choose visit pair with best paired coverage across arms
    visit_levels = sorted(adata.obs["dfo_bin"].dropna().unique())
    best_pair = None
    best_score = -1
    min_per_arm = 2
    for v1, v2 in itertools.combinations(visit_levels, 2):
        wide = (
            adata.obs
            .groupby(["participant_id", "dfo_bin"], observed=True)
            .size()
            .unstack(fill_value=0)
        )
        paired_ids = wide[(wide.get(v1, 0) > 0) & (wide.get(v2, 0) > 0)].index
        if len(paired_ids) == 0:
            continue
        arm_counts = (
            adata.obs[adata.obs["participant_id"].isin(paired_ids)]
            .groupby("severity")["participant_id"]
            .nunique()
        )
        if len(arm_counts) < 2:
            continue
        min_arm = min(arm_counts.get("Severe", 0), arm_counts.get("Mild", 0))
        if min_arm < min_per_arm:
            continue
        if min_arm > best_score:
            best_score = min_arm
            best_pair = (v1, v2)

    allow_unpaired = False
    if best_pair is None:
        print("No visit pair has >=2 paired participants per arm. DiD may be underpowered.")
        visits_pair = tuple(visit_levels[:2])
        allow_unpaired = True
        # Fallback: repeated-cross-sectional OLS DiD (weaker than paired)
    else:
        visits_pair = best_pair
    print("Using visits:", visits_pair)

    # Filter lineages with enough paired participants for the chosen pair
    paired_counts = (
        adata.obs[adata.obs["dfo_bin"].isin(visits_pair)]
        .groupby(["participant_id", "lineage"])["dfo_bin"]
        .nunique()
        .reset_index(name="n_visits")
    )
    paired_lineages = paired_counts[paired_counts["n_visits"] >= 2]["lineage"].value_counts()
    keep_lineages = paired_lineages[paired_lineages >= 2].index.tolist()
    if not keep_lineages:
        keep_lineages = sorted(adata.obs["lineage"].dropna().unique())

    print("Lineages with paired participants:", keep_lineages)

    adata_use = adata[adata.obs["lineage"].isin(keep_lineages)].copy()
    pb_mod = st.module_score_pseudobulk(
        adata_use,
        module_cols=signature_cols,
        design=design,
        visits=visits_pair,
        pool_col="lineage",
        min_cells_per_group=1,
    )
    display(pb_mod.head())

    def _run_mod_did(pb_df, allow_unpaired_flag):
        return st.module_score_did_by_pool(
            pb_df,
            design=design,
            visits=visits_pair,
            min_paired=2,
            n_perm=300,
            seed=SEED,
            fdr_within="module",
            allow_unpaired=allow_unpaired_flag,
        )

    res_mod = _run_mod_did(pb_mod, allow_unpaired)
    if res_mod.empty and not allow_unpaired:
        print("No valid paired lineage-level DiD results. Retrying with unpaired OLS DiD.")
        res_mod = _run_mod_did(pb_mod, True)

    if res_mod.empty:
        print("No valid lineage-level DiD results. Falling back to pooled lineage.")
        pb_all = st.module_score_pseudobulk(
            adata_use,
            module_cols=signature_cols,
            design=design,
            visits=visits_pair,
            pool_map={"All_Lineages": keep_lineages},
            celltype_col="lineage",
            min_cells_per_group=1,
        )
        res_mod = st.module_score_did_by_pool(
            pb_all,
            design=design,
            visits=visits_pair,
            min_paired=2,
            n_perm=300,
            seed=SEED,
            fdr_within=None,
            allow_unpaired=True,
        )

    if not res_mod.empty:
        display(res_mod.sort_values("p_DiD").head(20))
        pivot = res_mod.pivot(index="module", columns="pool", values="beta_DiD")
        plt.figure(figsize=(8, 5))
        sns.heatmap(pivot, cmap="RdBu_r", center=0)
        plt.title("Module-Score DiD by Lineage")
        plt.tight_layout()
        plt.show()
    else:
        print("No valid module-score DiD results after fallback.")
else:
    print("No module scores available for pseudobulk DiD.")

No visit pair has >=2 paired participants per arm. DiD may be underpowered.
Using visits: ('DFO_0-7', 'DFO_15+')
Lineages with paired participants: ['B_cells', 'CD4_T', 'CD8_T', 'DCs', 'Monocytes', 'NK', 'Other']
participant_id dfo_bin severity pool module module_score n_cells
0 AP1 DFO_15+ Severe B_cells sig_Antigen Presentation 0.255384 209
1 AP1 DFO_15+ Severe B_cells sig_Apoptosis 0.055282 209
2 AP1 DFO_15+ Severe B_cells sig_Cell Proliferation 0.126665 209
3 AP1 DFO_15+ Severe B_cells sig_Checkpoint Receptors -0.074214 209
4 AP1 DFO_15+ Severe B_cells sig_Cytotoxic T Cell Activity -0.429571 209 
/var/folders/xl/pr40r62n3wx54s8_6cgjzg2c0000gt/T/ipykernel_99779/2438914800.py:67: UserWarning: FDR correction is applied within each 'module' group (column FDR_DiD). Per-group FDR does not control the overall false discovery rate across all tests. Consult FDR_DiD_global (when fdr_global=True) for a globally corrected q-value.
  return st.module_score_did_by_pool(
pool module mean_delta_treated mean_delta_control beta_DiD p_DiD p_treated p_control n_units FDR_DiD FDR_DiD_global
9 B_cells sig_Monocyte Macrophage 0.012421 -0.035665 0.048086 0.006500 1.361194e-01 0.020210 18 0.045500 0.718733
37 CD4_T sig_Type II Interferon 0.057539 -0.116316 0.173855 0.019131 4.203365e-01 0.000354 18 0.050437 0.718733
94 Monocytes sig_Type II Interferon -0.121601 -0.441053 0.319452 0.020166 2.247482e-01 0.000007 18 0.050437 0.718733
56 CD8_T sig_Type II Interferon 0.080186 -0.099801 0.179987 0.021616 2.476298e-01 0.021417 18 0.050437 0.718733
28 CD4_T sig_Monocyte Macrophage 0.004282 -0.041806 0.046088 0.039128 7.902376e-01 0.009327 18 0.136947 0.875250
50 CD8_T sig_Regulatory T Cell 0.060579 -0.028357 0.088936 0.046339 9.530340e-02 0.317496 18 0.215342 0.875250
31 CD4_T sig_Regulatory T Cell 0.071323 -0.042156 0.113480 0.061526 1.924073e-01 0.191786 18 0.215342 0.875250
18 B_cells sig_Type II Interferon -0.004134 -0.097863 0.093729 0.073625 9.330677e-01 0.000089 18 0.095046 0.875250
113 NK sig_Type II Interferon 0.051381 -0.066763 0.118144 0.080356 3.456472e-01 0.124200 18 0.095046 0.875250
75 DCs sig_Type II Interferon 0.099679 -0.180316 0.279996 0.081468 4.828776e-01 0.043298 18 0.095046 0.875250
104 NK sig_Monocyte Macrophage 0.026409 -0.010286 0.036695 0.088830 2.208636e-02 0.569794 18 0.165750 0.875250
25 CD4_T sig_Hypoxia 0.159209 0.022439 0.136769 0.093014 5.103041e-02 0.444655 18 0.475952 0.875250
107 NK sig_Regulatory T Cell 0.118349 0.039809 0.078540 0.094175 2.974803e-05 0.289469 18 0.219741 0.875250
51 CD8_T sig_T Cell Activation 0.191334 0.058516 0.132819 0.102992 9.464380e-03 0.172901 18 0.368526 0.875250
43 CD8_T sig_Humoral Plasma Cell 0.087774 -0.011407 0.099180 0.104823 1.240236e-01 0.701375 18 0.449847 0.875250
108 NK sig_T Cell Activation 0.199677 0.097204 0.102473 0.105293 2.076444e-09 0.068507 18 0.368526 0.875250
123 Other sig_Monocyte Macrophage 0.032530 -0.025789 0.058319 0.130064 1.300227e-01 0.418430 18 0.165750 0.933005
44 CD8_T sig_Hypoxia 0.167578 0.058348 0.109231 0.135986 2.161701e-02 0.033359 18 0.475952 0.933005
66 DCs sig_Monocyte Macrophage 0.495102 0.218690 0.276412 0.140132 6.280042e-04 0.083729 18 0.165750 0.933005
47 CD8_T sig_Monocyte Macrophage -0.007196 -0.037414 0.030218 0.142071 5.721417e-01 0.021685 18 0.165750 0.933005
../_images/tutorials_example_covid19_stephenson_42_4.png