This page was generated from docs/tutorial/replogle/replogle-py.ipynb. Interactive online version: .

Batch-wise causarray analysis: Replogle K562

This tutorial analyses a subset of the Replogle et al. CRISPRi screen, compares adjusted log-fold changes with a marginal Wilcoxon analysis, and checks propensity-score support.

1_prep_tutorial_data.py prepares raw counts for the most abundant perturbations and a shared non-targeting control pool. gcate_lfc_batch fits subsets of perturbations with shared controls to limit the size of each fit.

The workflow is prepare inputs → select a factor rank → fit batches → check propensity support and refit where it fails → read the effects. Support is checked before the effects are read, because an arm whose inverse-probability weights collapse can report a large discovery count that reflects the weights rather than the biology. The comparison and propensity sections use crispyx and local caches; recompute dependent caches after changing data, settings, or package versions.

import sys
sys.path.insert(0, '../../..')

import time
from pathlib import Path
import numpy as np
import pandas as pd
import scipy.sparse as sp
import matplotlib.pyplot as plt
import scanpy as sc

from causarray import prep_causarray_data, gcate_lfc_batch, align_test_mask
from causarray.gcate import plot_r

import causarray
print('causarray version:', causarray.__version__)

# Analysis-wide significance threshold, used by the diagnostics and the results.
FDR = 0.05

causarray version: 0.1.1

1. Load data

Load the subset prepared by 1_prep_tutorial_data.py. The output reports its dimensions and treatment-group sizes.

adata = sc.read_h5ad('data/replogle_subset.h5ad')
print(adata)

CTRL_LABEL = adata.uns['ctrl_label']   # 'non-targeting'
PERT_COL   = adata.uns['pert_col']     # 'gene'

vc = adata.obs[PERT_COL].value_counts()
print(f'\nCtrl cells  : {vc[CTRL_LABEL]}')
print(f'Pert cells  : {vc.drop(CTRL_LABEL).sum():,}  across {len(vc) - 1} perturbations')
print(f'Cells/pert  : median {vc.drop(CTRL_LABEL).median():.0f}, '
      f'range {vc.drop(CTRL_LABEL).min()}-{vc.drop(CTRL_LABEL).max()}')

AnnData object with n_obs × n_vars = 79865 × 8563
    obs: 'batch', 'gene', 'gene_id', 'transcript', 'gene_transcript', 'guide_id', 'percent_mito', 'UMI_count', 'z_gemgroup_UMI', 'core_scale_factor', 'core_adjusted_UMI_count', 'disease', 'cancer', 'cell_line', 'sex', 'age', 'perturbation', 'organism', 'perturbation_type', 'tissue_type', 'ncounts', 'ngenes', 'nperts', 'percent_ribo'
    var: 'chr', 'start', 'end', 'class', 'strand', 'length', 'in_matrix', 'mean', 'std', 'cv', 'fano', 'ensembl_id', 'ncounts', 'ncells'
    uns: 'ctrl_label', 'pert_col'

Ctrl cells  : 2000
Pert cells  : 77,865  across 200 perturbations
Cells/pert  : median 346, range 273-1996

2. Prepare causarray inputs

prep_causarray_data caps extreme counts, adds an intercept, and includes standardized log-library size in the propensity design X_A. The batch fit uses that design; section 5 explains why the final propensity model leaves library size out. The outcome matrix contains counts, not log-normalized expression.

Y_raw = adata.X.toarray() if sp.issparse(adata.X) else np.array(adata.X)
Y = pd.DataFrame(Y_raw, columns=adata.var_names.tolist())
del Y_raw

A = (pd.get_dummies(adata.obs[PERT_COL].astype(str), drop_first=False)
       .drop(columns=[CTRL_LABEL]))

Y, A, X, X_A = prep_causarray_data(Y, A)

print(f'Y : {Y.shape}  (cells x genes)')
print(f'A : {A.shape}  (cells x perturbations)')
print(f'X : {X.shape}  X_A : {X_A.shape}')

Y : (79865, 8563)  (cells x genes)
A : (79865, 200)  (cells x perturbations)
X : (79865, 1)  X_A : (79865, 2)

3. Number of latent factors

The JIC table from 2_estimate_r.py evaluates candidate ranks on a control-heavy subsample of the raw-count data. Select the rank with the smallest JIC and inspect nearby ranks when the curve is flat. Subsampling reduces computation; it need not select the same rank as a full-data fit.

A flat JIC curve does not establish that confounding is absent.

# JIC table on the raw counts (2_estimate_r.py)
df_r = pd.read_csv('results/replogle-r.csv')
fig = plot_r(df_r)
plt.tight_layout()
plt.show()

best_r = int(df_r.loc[df_r['JIC'].idxmin(), 'r'])
print(f'\nSelected r = {best_r}  (min JIC)')
print(df_r.to_string(index=False))

Selected r = 10  (min JIC)
 r  deviance       nu       JIC     time_s
 0 -3.931351 0.291434 -3.639917  34.124589
 5 -4.042865 0.298683 -3.744181 105.917385
10 -4.053505 0.305933 -3.747573 110.224151
15 -4.059356 0.313183 -3.746173 113.079680
20 -4.064079 0.320432 -3.743646 120.366626
25 -4.068070 0.327682 -3.740389 124.748423
30 -4.071637 0.334931 -3.736706 144.562551

4. Batch fitting with gcate_lfc_batch

batch_size controls the number of perturbations per fit; max_cells caps perturbed cells and n_ctrl sets the shared control pool. The function distributes perturbations across batches and writes results to cache_path so a run can resume.

Use a fresh cache when changing the data, rank, fitting settings, or inference version. A cached result is not validation of the current implementation.

R = best_r

t0 = time.perf_counter()
df_res = gcate_lfc_batch(
    Y, X, A, R,
    W_A=X_A,
    batch_size=15,
    max_cells=2000,
    n_ctrl=2000,
    family='nb',
    cache_path='results/replogle_results.h5',
    random_state=0,
    verbose=True,
    gcate_kwargs=dict(
        kwargs_es_1=dict(rel_tol=2e-4, max_iters=30),
        kwargs_es_2=dict(rel_tol=2e-4, max_iters=30),
    ),
)
t_total = time.perf_counter() - t0

print(f'\nTotal wall time: {t_total/60:.1f} min')
print(f'Result shape   : {df_res.shape}')
df_res.head()

[gcate_lfc_batch] Resuming: 14 batches already cached in 'results/replogle_results.h5'
'Pre-estimating dispersion on ctrl cell subsample...'

GCATE batches: 100%|██████████| 14/14 [00:00<00:00, 26973.02batch/s]

Total wall time: 0.1 min
Result shape   : (1712600, 22)

	gene_names	tau	std	log2fc	log2fc_se	stat	rej	pvalue	padj	pvalue_emp_null_adj	...	n_control	count_treated	count_control	mean_control	mean_treated	estimable	var_floored	std_raw	trt	batch
0	LINC01409	0.323101	0.281189	0.466136	0.405669	1.149054	0.0	0.250666	0.535747	0.743786	...	2000	12.0	248.0	0.114105	0.157625	True	True	0.278146	AATF	0
1	LINC01128	-0.095688	0.247474	-0.138048	0.357030	-0.386657	0.0	0.699050	0.865696	0.450997	...	2000	20.0	546.0	0.256268	0.232883	True	False	0.247474	AATF	0
2	NOC2L	0.096539	0.081480	0.139277	0.117551	1.184821	0.0	0.236223	0.518468	0.724830	...	2000	219.0	4278.0	1.964351	2.163443	True	False	0.081480	AATF	0
3	KLHL17	0.090186	0.286134	0.130111	0.412804	0.315187	0.0	0.752651	0.896647	0.794920	...	2000	16.0	300.0	0.137459	0.150432	True	True	0.271956	AATF	0
4	HES4	0.214933	0.206852	0.310083	0.298425	1.039068	0.0	0.298894	0.588576	0.803012	...	2000	30.0	499.0	0.233575	0.289582	True	True	0.198972	AATF	0

5 rows × 22 columns

5. Propensity model: diagnostics and factor selection

Read this before the effects. Where the propensity model separates an arm from the controls, the inverse-probability weights rest on a few cells, and the arm’s discoveries then reflect those cells rather than the perturbation.

Guides are assigned at random, so the propensity model only has to absorb chance imbalance between an arm and the controls. Two consequences for its covariates:

• Library size is left out, as a judgement call. It could be a confounder, since capture depth and cell quality affect both guide detection and measured expression, or a consequence of the knockdown: knocking down transcription machinery shrinks a cell’s total RNA, and SUPT5H cells have about a third of the control library size. These data cannot tell which. With library size in the model such arms separate from the controls almost completely, so their weights would rest on a few cells. Leaving it out keeps them supported, at the cost of not adjusting for depth-related confounding; the outcome model still normalizes for depth through the size factors. The LFC documentation discusses this choice.

• Latent factors are included, then checked arm by arm. Factors estimated from expression can also track a perturbation’s own effect. When one does, the model separates that arm, and the weights pile onto the perturbed cells that look least perturbed.

The check flags an arm when its treated effective sample size (ESS) falls below half the arm, its histogram overlap with the controls below 0.3, or its AUC above 0.9. For a flagged arm, the factor most imbalanced between the arm and the controls (largest absolute standardized mean difference) is dropped from that arm’s model, and the step repeats until ESS ≥ 0.5 and overlap ≥ 0.3. Other arms keep every factor.

4_refit_propensity.py applies this to all 14 batches and produces the results this notebook reads (section 4). The cells below walk through the batch containing SUPT5H. They read that batch’s latent factors from results/replogle_results_r10.nuisances.h5, a large local file written by 3_run_batch.py; if it is absent, they show the per-arm report from 4_refit_propensity.py instead.

import h5py
from causarray import (estimate_propensity_scores, summarize_propensity_scores,
                       summarize_treatment_associations, select_propensity_factors,
                       refit_propensity_scores)

# The batch fit saved each batch's latent factors and cells in a large local
# file (tens of GB, not tracked). Without it, skip the walk-through and show
# the per-arm report that 4_refit_propensity.py wrote for all batches.
NUISANCE_PATH = Path(f'results/replogle_results_r{R}.nuisances.h5')
HAS_FACTORS = NUISANCE_PATH.exists()
if HAS_FACTORS:
    # Read the batch that contains SUPT5H (its outcome predictions are not needed).
    with h5py.File(NUISANCE_PATH, 'r') as handle:
        for key in sorted(handle.keys()):
            names = [n.decode() if isinstance(n, bytes) else n for n in handle[key]['pert_names'][:]]
            if 'SUPT5H' in names:
                U_focal = handle[key]['U'][:]
                cells_focal = handle[key]['cell_idx'][:]
                focal_perts = names
                break
    A_focal = A.iloc[cells_focal][focal_perts].to_numpy(dtype=float)
    W_A_focal = np.c_[X_A[cells_focal, :1], U_focal]     # intercept + latent factors
    focal_covariates = ['intercept'] + [f'U{k + 1}' for k in range(U_focal.shape[1])]
    print(f'{key}: {len(focal_perts)} perturbations, {len(cells_focal):,} cells')

    pi_before = estimate_propensity_scores(
        A_focal, W_A_focal, K=1, C=1.0, class_weight=None, clip=None, random_state=0)
    support_before = summarize_propensity_scores(
        A_focal, pi_before, treatment_names=focal_perts, clip_bounds=None).set_index('treatment')
    display(support_before[['prevalence', 'auc', 'overlap_ratio', 'ess_treated_fraction']]
            .sort_values('ess_treated_fraction').head(8).round(3))
else:
    print(f'{NUISANCE_PATH} not found; skipping the propensity walk-through.')
    report_path = Path('results/replogle_propensity_selection.csv')
    if report_path.exists():
        report = pd.read_csv(report_path, index_col='treatment')
        display(report[report['flagged']].round(3))

batch_0012: 14 perturbations, 4,000 cells

	prevalence	auc	overlap_ratio	ess_treated_fraction
treatment
SUPT6H	0.046	0.689	0.331	0.250
SUPT5H	0.070	0.714	0.385	0.358
TIMM23B	0.078	0.595	0.624	0.830
TSR2	0.083	0.614	0.622	0.851
SSRP1	0.093	0.586	0.536	0.853
SON	0.075	0.455	0.632	0.888
SNRPG	0.068	0.512	0.574	0.914
TRRAP	0.063	0.571	0.797	0.936

if HAS_FACTORS:
    # The four arms whose weights are most concentrated, with every factor included.
    ps_focus = support_before.nsmallest(4, 'ess_treated_fraction').index.tolist()
    print('lowest treated ESS:', ps_focus)

    def overlap_panels(pi, title):
        ctrl_mask = A_focal.sum(axis=1) == 0
        fig, axes = plt.subplots(1, len(ps_focus), figsize=(3.6 * len(ps_focus), 3.4), sharey=True)
        for ax, pert in zip(np.atleast_1d(axes), ps_focus):
            j = focal_perts.index(pert)
            treated = A_focal[:, j] == 1
            hi = max(pi[ctrl_mask | treated, j].max(), 1e-3)
            bins = np.linspace(0, hi, 40)
            ax.hist(pi[ctrl_mask, j], bins=bins, density=True, histtype='step',
                    color='#4c78a8', label='control')
            ax.hist(pi[treated, j], bins=bins, density=True, histtype='step',
                    color='#e45756', label='treated')
            ax.set_title(pert, fontsize=10)
            ax.set_xlabel('propensity score')
            for side in ('top', 'right'):
                ax.spines[side].set_visible(False)
        np.atleast_1d(axes)[0].set_ylabel('density')
        np.atleast_1d(axes)[0].legend(frameon=False, fontsize=8)
        fig.suptitle(title, y=1.02)
        fig.tight_layout()
        plt.show()

    overlap_panels(pi_before, 'Before: every latent factor included')

lowest treated ESS: ['SUPT6H', 'SUPT5H', 'TIMM23B', 'TSR2']

Which factors separate the flagged arms

The standardized mean difference (SMD) of each factor between an arm and the controls shows what the propensity model is picking up. With random guides a factor should differ little between the arm and the controls. Here no single factor differs much (absolute SMD below 0.3), but together they let the propensity model pick out part of the arm, and that is enough to concentrate its weights. The table lists, for each flagged arm, its most imbalanced factors; the selection drops them in that order.

if HAS_FACTORS:
    smd = summarize_treatment_associations(
        A_focal, W_A_focal[:, 1:], treatment_names=focal_perts, covariate_names=focal_covariates[1:])
    smd['abs_smd'] = smd['standardized_mean_difference'].abs()
    flag = ((support_before['ess_treated_fraction'] < 0.5) | (support_before['overlap_ratio'] < 0.3)
            | (support_before['auc'] > 0.9))
    flagged = support_before.index[flag].tolist()
    display(smd[smd['treatment'].isin(flagged)]
            .sort_values(['treatment', 'abs_smd'], ascending=[True, False])
            .groupby('treatment').head(3)
            [['treatment', 'covariate', 'standardized_mean_difference']].round(2))

    drops, selection = select_propensity_factors(
        A_focal, W_A_focal, treatment_names=focal_perts, covariate_names=focal_covariates)
    pi_after, _ = refit_propensity_scores(
        A_focal, W_A_focal, pi_hat=pi_before.copy(), treatment_names=focal_perts,
        covariate_names=focal_covariates, drop_by_treatment=drops,
        K=1, C=1.0, class_weight=None, clip=None, random_state=0)
    display(selection.set_index('treatment')
            [['dropped', 'auc_all', 'auc_chosen', 'overlap_ratio_all', 'overlap_ratio_chosen',
              'ess_treated_fraction_all', 'ess_treated_fraction_chosen']].round(3))

	treatment	covariate	standardized_mean_difference
69	SUPT5H	U10	-0.26
68	SUPT5H	U9	0.15
64	SUPT5H	U5	-0.15
79	SUPT6H	U10	-0.26
78	SUPT6H	U9	0.19
74	SUPT6H	U5	-0.14

	dropped	auc_all	auc_chosen	overlap_ratio_all	overlap_ratio_chosen	ess_treated_fraction_all	ess_treated_fraction_chosen
treatment
SUPT5H	U10	0.714	0.639	0.385	0.550	0.358	0.63
SUPT6H	U10,U9	0.689	0.668	0.331	0.536	0.250	0.62

if HAS_FACTORS:
    overlap_panels(pi_after, 'After: imbalanced factors dropped on flagged arms only')

6. Results

6.1 Discovery summary

sig = df_res[df_res['padj'] < FDR]

print(f'Significant (padj < {FDR}): {len(sig):,} gene x pert pairs')
print(f'  Perturbations with >= 1 hit: {sig["trt"].nunique()}')
print(f'  Unique genes affected      : {sig["gene_names"].nunique():,}')

disc_per_pert = sig.groupby('trt').size().sort_values(ascending=False)
print(f'\nTop-10 perts by discovery count:')
print(disc_per_pert.head(10).to_string())

Significant (padj < 0.05): 145,340 gene x pert pairs
  Perturbations with >= 1 hit: 194
  Unique genes affected      : 8,388

Top-10 perts by discovery count:
trt
CSE1L    3078
MED12    3068
TSR2     2940
LSM6     2524
NCBP2    2396
DHX15    2302
HSPA9    2169
INTS2    2138
TRRAP    2137
MED30    1983

6.2 Volcano plot (all perturbations combined)

fig, ax = plt.subplots(figsize=(7, 5))

colors = np.where(df_res['padj'] < FDR, '#e74c3c', '#aaaaaa')
ax.scatter(
    df_res['tau'], -np.log10(df_res['pvalue'].clip(1e-300)),
    c=colors, s=1, alpha=0.4, rasterized=True,
)
ax.axhline(-np.log10(0.05 / len(df_res)), color='navy', lw=0.8, ls='--',
           label='Bonferroni')
ax.set_ylim(0, 10)
ax.set_xlabel('Estimated LFC (tau)')
ax.set_ylabel('-log10(p-value)')
ax.set_title('Volcano plot — Replogle-E-K562 (top-200 perts)')
ax.legend(fontsize=8)
plt.tight_layout()
plt.show()

6.3 Discovery count per perturbation

Each bar counts genes passing the BH-adjusted threshold for one perturbation. Differences can reflect biological response, sample size, expression support, or model fit; the counts alone do not distinguish these explanations.

top_n = 30
top_disc = disc_per_pert.head(top_n)

fig, ax = plt.subplots(figsize=(10, 4))
top_disc.plot(kind='bar', ax=ax, color='steelblue', edgecolor='none')
ax.set_ylabel(f'Significant genes (padj < {FDR})')
ax.set_title(f'Top-{top_n} perturbations by discovery count')
ax.tick_params(axis='x', rotation=45, labelsize=8)
plt.tight_layout()
plt.show()

6.4 Comparison with a marginal Wilcoxon test

As a reference point, each perturbation is also compared with the pooled non-targeting controls using a Wilcoxon rank-sum test on library-size-normalized, log1p counts (via the optional crispyx package). Both methods are summarized at BH-adjusted p < 0.05.

Wilcoxon is marginal: unlike causarray it adjusts for neither the observed covariates nor the estimated latent factors. Gaps in discovery counts can therefore reflect confounding, model assumptions, filtering, or power, and do not by themselves mark discoveries as true or false.

import crispyx

# Wilcoxon is applied to normalized, log-transformed counts, following the
# same convention as the Adamson tutorial.
norm_path = Path('data/replogle_subset_norm.h5ad')
if not norm_path.exists():
    crispyx.normalize_total_log1p(
        'data/replogle_subset.h5ad', output_path=str(norm_path), verbose=False
    )
# The Wilcoxon result does not depend on causarray; cache it on disk and reuse
# it on later runs (crispyx skips the computation when the output exists).
wc_result = crispyx.wilcoxon_test(
    str(norm_path),
    perturbation_column=PERT_COL,
    control_label=CTRL_LABEL,
    verbose=False,
    output_path='results/replogle_subset_norm_cx_wilcoxon.h5ad',
)

wc_count_map = {
    str(pert): int((np.asarray(wc_result[pert].pvalue_adj) < FDR).sum())
    for pert in wc_result.groups
}
all_perts = pd.Index(A.columns.astype(str), name='Perturbation')
comparison_counts = pd.DataFrame(index=all_perts)
comparison_counts['causarray'] = (
    disc_per_pert.rename(index=str).reindex(all_perts, fill_value=0).astype(int)
)
comparison_counts['Wilcoxon'] = (
    pd.Series(wc_count_map).reindex(all_perts, fill_value=0).astype(int)
)
comparison_counts = comparison_counts.sort_values(
    ['Wilcoxon', 'causarray'], ascending=False
)

print(f'causarray: {comparison_counts["causarray"].sum():,} significant pairs')
print(f'Wilcoxon : {comparison_counts["Wilcoxon"].sum():,} significant pairs')
print('Median discoveries per perturbation: ' +
      f'causarray={comparison_counts["causarray"].median():.0f}, ' +
      f'Wilcoxon={comparison_counts["Wilcoxon"].median():.0f}')
display(comparison_counts.head(10))

causarray: 145,340 significant pairs
Wilcoxon : 425,125 significant pairs
Median discoveries per perturbation: causarray=602, Wilcoxon=2220

	causarray	Wilcoxon
Perturbation
SUPT5H	1591	6842
SUPT6H	1231	6476
MED12	3068	5883
MED30	1983	5458
INTS2	2138	5373
HSPA9	2169	5326
SFPQ	504	5188
DDX41	518	5049
MED19	1920	4946
TSR2	2940	4941

Comparing gene lists

The saved comparison reports more discoveries with Wilcoxon. To explore the difference, align results by perturbation and gene, then inspect direction-specific enrichment and expression support. Neither list size nor enrichment establishes which individual discoveries are correct.

wc_groups = pd.Index([str(group) for group in wc_result.groups], name='trt')
wc_genes = pd.Index(wc_result.genes.astype(str), name='gene_names')
wc_lfc = np.asarray(wc_result.logfoldchanges, dtype=float)
wc_padj = np.asarray(wc_result.pvalues_adj, dtype=float)
wc_pts_treated = np.asarray(wc_result.pts, dtype=float)
wc_pts_control = np.asarray(wc_result.pts_rest, dtype=float)

expected = pd.MultiIndex.from_product([wc_groups, wc_genes])
ca_indexed = (
    df_res.assign(
        trt=df_res['trt'].astype(str),
        gene_names=df_res['gene_names'].astype(str),
    ).set_index(['trt', 'gene_names'])
)
if ca_indexed.index.has_duplicates:
    raise ValueError('causarray contains duplicate perturbation-gene tests')
aligned_ca = ca_indexed.reindex(expected)
if aligned_ca['tau'].isna().all():
    raise ValueError('Wilcoxon and causarray results could not be aligned')

ca_lfc = (
    aligned_ca['tau'].to_numpy(dtype=float).reshape(wc_lfc.shape)
    / np.log(2.0)
)
ca_padj = aligned_ca['padj'].to_numpy(dtype=float).reshape(wc_lfc.shape)
ca_estimable = aligned_ca['estimable'].fillna(False).to_numpy().reshape(wc_lfc.shape)
wc_reject = wc_padj < FDR
ca_reject = ca_padj < FDR
del aligned_ca, ca_indexed, expected

print(f'Aligned {wc_lfc.size:,} perturbation-gene tests.')

Aligned 1,712,600 perturbation-gene tests.

Gene Ontology enrichment for SUPT5H

SUPT5H is the screen’s canonical case: knocking down the SPT5 elongation factor should hit translation and ribosome biogenesis, so its gene list has a known expectation to be read against. Its propensity weights were among the most concentrated before factor selection, so the cell prints that alongside the result rather than leaving it implicit. Both methods are tested on the same gene background.

# SUPT5H is the screen's canonical case: knocking down the SPT5 elongation
# factor should hit translation and ribosome biogenesis. Its propensity
# support is reported alongside the result.
ps_audit = pd.read_csv('results/replogle_propensity_selection.csv').set_index('treatment')
go_focus = 'SUPT5H'
_row = ps_audit.loc[go_focus]
print(f"GO focus: {go_focus} -- factors dropped: {_row['dropped'] if _row['flagged'] else 'none'}, "
      f"treated ESS {_row['ess_treated_before']:.2f} -> {_row['ess_treated_after']:.2f}, "
      f"overlap {_row['overlap_before']:.2f} -> {_row['overlap_after']:.2f}")
print(f"({int((~ps_audit['flagged']).sum())} of {len(ps_audit)} arms kept every factor)")

go_result_path = Path(f'results/replogle_{go_focus.lower()}_go_results.csv')
focus_i = wc_groups.get_loc(go_focus)
go_queries = {
    'causarray up': wc_genes[ca_reject[focus_i] & (ca_lfc[focus_i] > 0)].tolist(),
    'causarray down': wc_genes[ca_reject[focus_i] & (ca_lfc[focus_i] < 0)].tolist(),
    'Wilcoxon up': wc_genes[wc_reject[focus_i] & (wc_lfc[focus_i] > 0)].tolist(),
    'Wilcoxon down': wc_genes[wc_reject[focus_i] & (wc_lfc[focus_i] < 0)].tolist(),
}
if go_result_path.exists():
    go_results = pd.read_csv(go_result_path)
else:
    go_results = sc.queries.enrich(
        go_queries, org='hsapiens',
        gprofiler_kwargs={
            'sources': ['GO:BP'],
            'background': wc_genes.tolist(),
            'no_evidences': False,
        },
    )
    go_results.to_csv(go_result_path, index=False)

go_query_summary = pd.DataFrame({
    'DE_genes': {query: len(genes) for query, genes in go_queries.items()},
    'significant_GO_BP_terms': (
        go_results.loc[go_results['significant']].groupby('query').size()
    ),
}).fillna(0).astype(int)
display(go_query_summary)

go_top = (
    go_results.loc[go_results['significant'] & go_results['term_size'].between(20, 1000)]
    .sort_values(['query', 'p_value'])
    .groupby('query', group_keys=False).head(5)
)
display(go_top[['query', 'name', 'p_value', 'intersection_size', 'query_size']])

if not go_top.empty:
    plotted = go_top.sort_values('p_value', ascending=False).copy()
    plotted['label'] = plotted['name'] + ' — ' + plotted['query']
    colors = plotted['query'].map({
        'causarray up': '#1976d2', 'causarray down': '#64b5f6',
        'Wilcoxon up': '#ef6c00', 'Wilcoxon down': '#ffb74d',
    })
    fig, ax = plt.subplots(figsize=(9, max(4, 0.35 * len(plotted))))
    ax.scatter(
        -np.log10(plotted['p_value'].clip(1e-300)), plotted['label'],
        s=25 + 4 * plotted['intersection_size'], c=colors, alpha=0.8,
    )
    ax.set_xlabel('-log10(g:Profiler adjusted p-value)')
    ax.set_ylabel('')
    ax.set_title(f'{go_focus}: selected GO Biological Process enrichment')
    fig.tight_layout()
    plt.show()

GO focus: SUPT5H -- factors dropped: U10, treated ESS 0.36 -> 0.63, overlap 0.38 -> 0.55
(171 of 200 arms kept every factor)

	DE_genes	significant_GO_BP_terms
Wilcoxon down	5646	22
Wilcoxon up	1196	1
causarray down	642	71
causarray up	949	0

	query	name	p_value	intersection_size	query_size
10	Wilcoxon down	ribosome biogenesis	8.496590e-13	249	5430
13	Wilcoxon down	ribonucleoprotein complex biogenesis	9.065321e-12	351	5430
17	Wilcoxon down	rRNA metabolic process	5.047599e-10	201	5430
19	Wilcoxon down	rRNA processing	3.027228e-09	178	5430
25	Wilcoxon down	translation	2.632631e-07	421	5430
91	Wilcoxon up	proton transmembrane transport	3.728893e-02	37	1124
0	causarray down	translation	1.859947e-46	150	621
1	causarray down	protein biosynthetic process	1.859947e-46	150	621
2	causarray down	cytoplasmic translation	1.099410e-43	85	621
3	causarray down	ribonucleoprotein complex biogenesis	2.308388e-30	112	621
4	causarray down	ribosome biogenesis	1.501356e-23	82	621

In the saved enrichment results, both down-regulated lists highlight ribosome biogenesis and RNA-related processes; the up-regulated lists show less enrichment. This shared theme is more informative than comparing the number of significant terms. It does not establish that one method has fewer false discoveries.

Effect-size agreement and extreme Wilcoxon effects

Near-empty groups can produce large marginal log2 fold changes that depend strongly on the pseudocount. The full-range and zoomed panels show the tail and central agreement separately. Both axes use log2 units; colors indicate rejection at the chosen BH-adjusted threshold.

selected_perts = ['SRRT', 'MED30', 'INTS2']
category_colors = {
    'Neither': '#bdbdbd',
    'Wilcoxon only': '#64b5f6',
    'causarray only': '#ffb74d',
    'Both': '#ef5350',
}

def plot_lfc_comparison(limits=None):
    fig, axes = plt.subplots(1, len(selected_perts), figsize=(15.5, 4.8))
    for ax, perturbation in zip(axes, selected_perts):
        i = wc_groups.get_loc(perturbation)
        x, y = wc_lfc[i], ca_lfc[i]
        finite = np.isfinite(x) & np.isfinite(y)
        categories = np.select(
            [
                wc_reject[i] & ca_reject[i],
                wc_reject[i] & ~ca_reject[i],
                ~wc_reject[i] & ca_reject[i],
            ],
            ['Both', 'Wilcoxon only', 'causarray only'],
            default='Neither',
        )
        for category, color in category_colors.items():
            keep = finite & (categories == category)
            ax.scatter(
                x[keep], y[keep], s=7 if category == 'Neither' else 11,
                alpha=0.22 if category == 'Neither' else 0.55,
                color=color, edgecolors='none', rasterized=True,
                label=f'{category} (n={keep.sum():,})',
            )
        extreme = finite & (x < -5)
        if extreme.any():
            ax.scatter(
                x[extreme], y[extreme], s=55, facecolors='none',
                edgecolors='black', linewidths=0.9, zorder=5,
            )
            for gene, xv, yv in zip(wc_genes[extreme], x[extreme], y[extreme]):
                ax.annotate(gene, (xv, yv), xytext=(4, 4),
                            textcoords='offset points', fontsize=7)
        pearson = np.corrcoef(x[finite], y[finite])[0, 1]
        spearman = pd.Series(x[finite]).corr(pd.Series(y[finite]), method='spearman')
        ax.axhline(0, color='#777777', linewidth=0.7)
        ax.axvline(0, color='#777777', linewidth=0.7)
        ax.axvline(-5, color='black', linestyle='--', linewidth=0.8)
        if limits is not None:
            ax.set_xlim(limits)
            ax.set_ylim(limits)
        ax.set_title(
            f'{perturbation}\nn={finite.sum():,}, r={pearson:.2f}, rho={spearman:.2f}'
        )
        ax.set_xlabel('Wilcoxon log2FC')
        ax.set_ylabel('causarray log2FC')
        ax.legend(fontsize=6.5, frameon=False, loc='best')
    suffix = ' (central range)' if limits is not None else ''
    fig.suptitle('Wilcoxon versus causarray effects' + suffix)
    fig.tight_layout()
    return fig, axes

plot_lfc_comparison()
plt.show()
plot_lfc_comparison(limits=(-3, 3))
plt.show()

Inspecting the negative effect tail

The following diagnostic selects Wilcoxon log2 fold changes below a fixed display threshold and reports detection counts in each group. Direct targets are flagged to distinguish possible on-target depletion from other responses.

group_counts = adata.obs[PERT_COL].astype(str).value_counts()
treated_sizes = np.asarray([group_counts[group] for group in wc_groups])[:, None]
control_size = int(group_counts[CTRL_LABEL])
treated_detected = np.rint(wc_pts_treated * treated_sizes).astype(int)
control_detected = np.rint(wc_pts_control * control_size).astype(int)
extreme_positions = np.argwhere(wc_lfc < -5)

extreme_rows = []
for i, j in extreme_positions:
    if wc_reject[i, j] and ca_reject[i, j]:
        decision = 'Both'
    elif wc_reject[i, j]:
        decision = 'Wilcoxon only'
    elif ca_reject[i, j]:
        decision = 'causarray only'
    else:
        decision = 'Neither'
    extreme_rows.append({
        'perturbation': wc_groups[i],
        'gene': wc_genes[j],
        'direct_target': wc_groups[i] == wc_genes[j],
        'treated_detected': f'{treated_detected[i, j]}/{treated_sizes[i, 0]}',
        'control_detected': f'{control_detected[i, j]}/{control_size}',
        'Wilcoxon_log2FC': wc_lfc[i, j],
        'Wilcoxon_padj': wc_padj[i, j],
        'causarray_log2FC': ca_lfc[i, j],
        'causarray_padj': ca_padj[i, j],
        'causarray_estimable': bool(ca_estimable[i, j]),
        'decision': decision,
    })
extreme_tests = pd.DataFrame(extreme_rows).sort_values('Wilcoxon_log2FC')
display(extreme_tests.round({
    'Wilcoxon_log2FC': 3, 'Wilcoxon_padj': 3,
    'causarray_log2FC': 3, 'causarray_padj': 3,
}))

	perturbation	gene	direct_target	treated_detected	control_detected	Wilcoxon_log2FC	Wilcoxon_padj	causarray_log2FC	causarray_padj	causarray_estimable	decision
2	SRRT	SDSL	False	0/307	153/2000	-25.379	0.0	-2.887	0.199	True	Wilcoxon only
3	MED30	SLCO4A1	False	1/288	237/2000	-5.511	0.0	-2.816	0.361	True	Wilcoxon only
0	SRRT	OPN3	False	1/307	227/2000	-5.470	0.0	-3.459	0.109	True	Wilcoxon only
1	SRRT	C6orf52	False	1/307	246/2000	-5.193	0.0	-3.585	0.094	True	Wilcoxon only
4	INTS2	INTS2	True	3/354	517/2000	-5.048	0.0	-4.757	0.020	True	Both

Candidate support rules

Compare the candidate support rules before adopting one. Requiring detections in both groups removes near-empty comparisons but can also remove genuine complete knockouts. A higher-expression-group rule retains those cases, so it answers a different question. align_test_mask attaches a rule by label without refitting.

negative = wc_lfc < 0
higher_expression_detected = np.where(negative, control_detected, treated_detected)
higher_expression_size = np.where(negative, control_size, treated_sizes)
support_masks = {
    'both groups >= 1 detected':
        (treated_detected >= 1) & (control_detected >= 1),
    'both groups >= 3 detected':
        (treated_detected >= 3) & (control_detected >= 3),
    'both groups >= 5 detected':
        (treated_detected >= 5) & (control_detected >= 5),
    'both groups >= 10 detected':
        (treated_detected >= 10) & (control_detected >= 10),
    'both groups >= max(10, 1%)': (
        (treated_detected >= np.maximum(10, np.ceil(0.01 * treated_sizes)))
        & (control_detected >= max(10, int(np.ceil(0.01 * control_size))))
    ),
    'higher-expression group >= max(10, 1%)': (
        higher_expression_detected
        >= np.maximum(10, np.ceil(0.01 * higher_expression_size))
    ),
}

extreme_mask = wc_lfc < -5
rule_rows = []
for rule, keep in support_masks.items():
    removed = ~keep
    rule_rows.append({
        'rule': rule,
        'tests_removed': int(removed.sum()),
        'tests_removed_percent': 100 * removed.mean(),
        'extreme_tests_removed': int((removed & extreme_mask).sum()),
        'Wilcoxon_discoveries_removed': int((removed & wc_reject).sum()),
        'causarray_discoveries_removed': int((removed & ca_reject).sum()),
    })
support_rule_impact = pd.DataFrame(rule_rows)
display(support_rule_impact.round({'tests_removed_percent': 3}))

chosen_rule = 'higher-expression group >= max(10, 1%)'
keep_in_df_res = align_test_mask(
    df_res, support_masks[chosen_rule],
    treatment_names=wc_groups, gene_names=wc_genes,
)
print(f'{keep_in_df_res.sum():,}/{len(keep_in_df_res):,} causarray rows retained')
# Materialize only if needed; no LFC refit is performed:
# df_res_screened = df_res.loc[keep_in_df_res].copy()

	rule	tests_removed	tests_removed_percent	extreme_tests_removed	Wilcoxon_discoveries_removed	causarray_discoveries_removed
0	both groups >= 1 detected	1	0.000	1	1	0
1	both groups >= 3 detected	12	0.001	4	12	0
2	both groups >= 5 detected	41	0.002	5	37	2
3	both groups >= 10 detected	392	0.023	5	331	19
4	both groups >= max(10, 1%)	392	0.023	5	331	19
5	higher-expression group >= max(10, 1%)	9	0.001	0	0	0

1,712,591/1,712,600 causarray rows retained

Reading the support diagnostic

The saved extreme tail contains comparisons with few treated detections and better control support. This explains why a both-group rule removes more of the tail than a higher-expression-group rule. Inspect effect sizes, uncertainty, and raw support together rather than automatically discarding strong depletion.

These post hoc masks describe how the gene lists respond to a support rule, not a newly validated testing procedure. Keep the original adjusted p-values when displaying retained rows; filtering does not by itself establish FDR control within the selected subset.

7. Summary

Batch fitting makes a large perturbation screen manageable. The saved comparisons show that adjusted and marginal analyses can yield different gene lists while sharing enrichment themes. Expression-support checks help interpret extreme effects, and propensity diagnostics reveal weight concentration that discovery counts alone miss.

Report the selected model, the support rules, and the per-arm factor choices together. Recheck these qualitative patterns whenever the analysis or caches are refreshed.