CIPHER package

Availability

The package is available on PyPI as cipher-perturb (PyPI page) or GitHub as cipher (GitHub page).

How do I use it?

1. Installation

   It can be installed from PyPI using the following:

   pip install cipher-perturb
   

   It can be installed from GitHub using the following:

   pip install -e .
   

   The core install is lightweight (numpy, scipy, pandas, anndata, h5py, tqdm). Optional extras add heavier, feature-specific dependencies:

   pip install -e ".[plot]"    # matplotlib plotting helpers (cipher.plotting)
   pip install -e ".[bayes]"   # PyMC horseshoe reverse (bayesian_reverse)
   pip install -e ".[dev]"     # tests (pytest)
   

   On the lab HPC the full conda environment (matching the exact versions used for the paper) can be built from the pinned spec:

   conda env create -f environment.yaml
   

2. Preparing the input sample sheet.

   The input dataset should be a single-cell RNA-seq dataset with chemical or genetic perturbation in .h5ad format, with the following requirements:

   • Should contain perturbation column in the adata.obs, which indicates the perturbation applied to each cell. The control cells should be labeled as control.

   • Should contain gene_names column in the adata.var, which indicates the gene names corresponding to each gene expression value.

3. Forward Prediction: ForwardResult

   This code provides the predicted transcriptome from the given perturbation set.

   import cipher
   
   res = cipher.forward_prediction(
       "path/to/perturbseq.h5ad",
       normalization="log1p",
       nulls=("meanfield", "shuffled"),    # baseline covariance models
       holdout_frac = 0.0,                 # 0.5 for out-of-sample gene holdout (paper setting)
       max_perturbation = None             #  int for a quick smoke test
   )
   
   res.results # DataFrame per perturbation: r2_uncentered_real, r2_centered_real
               # pearson_real, spearman_real, cosine_real, mse_real, sign_accuracy_real
               # a_hat, n_train_genes, n_test_genes, r2_uncentered_<null>, ...
   
   res.summary      # dict: mean_r2_uncentered_real, mean_pearson_real, mean_r2_uncentered_<>null>
   res.save("path/to/output_dir")
   

   You can also run forward metrics straight from a preprocessed directory (real Sigma only, matching the paper’s final forward recompute): cipher.forward_from_precomputed(“path/to/output_dir”, “log1p”, holdout_frac=0.5).

4. Reverse Prediction: ReverseResult

   This code identifies the causal perturbation from the Perturb-seq transcriptomic profile.

   import cipher
   
   res = cipher.reverse_prediction(
       "path/to/perturbseq.h5ad",
       normalization="log1p",
       method="matched_filter",   # matched_filter (default) | pinv | ridge | lstsq
       top_k=10,
   )
   
   res.summary["mean_auc"]                 # mean one-vs-rest ROC-AUC of the true driver
   res.summary["top10_accuracy"]           # fraction with the target gene in the top-k
   res.results.sort_values("auc").head()   # per-perturbation ranks / AUC
   res.save("path/to/output_dir")          # writes <dataset>_reverse_<norm>_<method>.csv
   

   cipher.reverse_from_precomputed(…) runs the same from a preprocessed directory. With the [bayes] extra, cipher.bayesian_reverse(Sigma, delta_x) fits a horseshoe prior and returns posterior inclusion probabilities.

5. Posterior inverse: InverseResult

   The strongest inverse solver is the fullH_diag posterior: it whitens dx by its per-perturbation sampling covariance (lambda/n0 + projected_var/nu), fits a prior variance by empirical Bayes, and scores each gene by its posterior perturbation strength (or a single-effect PIP). It runs on the artifacts written by preprocess_dataset (needs save_mean_var=True, the default) and is scored with pooled ROC / precision-recall curves.

   import cipher
   
   res = cipher.posterior_inverse_from_precomputed(
       "path/to/output_dir", "log1CP10k",
       method="posterior",         # "posterior" (default) or "pip"
   )
   res.summary["pooled_auc"]                 # pooled one-vs-rest ROC-AUC
   res.summary["mean_per_pert_auc"]          # mean over perturbations
   res.summary["pooled_average_precision"]   # pooled AP (precision-recall)
   res.roc, res.prc                          # (fpr, tpr), (precision, recall) for plotting
   
   # or in memory straight from an .h5ad (moderate datasets):
   res = cipher.posterior_inverse_prediction("path/to/perturbseq.h5ad",
                                           normalization="log1CP10k", method="posterior")
   

   To reproduce the paper’s per-dataset CRISPRi vs CRISPRa ROC/PR figure, run the inverse for each dataset, group with cipher.dataset_group(name), and draw each group with cipher.plotting.plot_inverse_group (thin per-dataset traces + a mean trace):

   import matplotlib.pyplot as plt
   from cipher import plotting
   
   results = {"CRISPRi": [...], "CRISPRa": [...]}   # InverseResult per dataset, by group
   fig, axes = plt.subplots(1, 2, figsize=(12, 6))
   plotting.plot_inverse_group(results["CRISPRi"], curve="roc", ax=axes[0], color="#d62728")
   plotting.plot_inverse_group(results["CRISPRa"], curve="roc", ax=axes[1], color="#1f77b4")
   

6. Condition drivers: DriverResult

   The reverse problem applied outside Perturb-seq: given a control group and a condition group, rank the genes most likely to drive the condition. There is no ground-truth label, so the output is a ranked candidate list.

   import cipher
   
   # From an .h5ad / AnnData with a grouping column in obs:
   res = cipher.condition_drivers(
       "path/to/dataset.h5ad",
       condition_key="stim",       # obs column grouping the cells
       control_value="rest",       # value marking control cells
       condition_value="stim",     # None => every non-control cell
       normalization="log1p",
       method="matched_filter",    # robust default (no matrix inverse)
   )
   res.top(20)                     # top-ranked candidate driver genes
   res.save("path/to/output_dir") # writes <name>_drivers_<norm>_<method>.csv
   
   # Or directly from two raw (cells x genes) matrices sharing a gene axis:
   res = cipher.condition_drivers_from_matrices(
       control_X, condition_X, gene_names,
       normalization="log1p",
       method="matched_filter",
   )
   

7. Lower-level Dataset

   For custom pipelines, load_dataset returns a Dataset that handles perturbation detection, target-gene mapping, and control/perturbation matrices.

   import cipher
   
   ds = cipher.load_dataset(
       "path/to/perturbseq.h5ad",
       expression_threshold=1.0,
       min_samples=100,
   )
   ds.n_genes, ds.n_perturbations
   ds.gene_names                       # np.array[str], length p
   ds.perturbations                    # list[str], non-control perturbation labels
   ds.target_gene_indices              # int64 array parallel to perturbations (-1 if absent)
   
   control = ds.control_matrix(dense=True)          # (n_control, p)
   pert = ds.perturbation_matrix(ds.perturbations[0], dense=True)
   g = ds.gene_index("GATA1")
   
   # Estimate the control covariance and run one perturbation forward by hand:
   Sigma = cipher.compute_covariance(cipher.normalize_matrix(control, "log1p"))
   

Command-line usage

Installing the package registers a cipher console script with four subcommands. Each takes an input .h5ad and an output directory (-o).

# Precompute Sigma + per-perturbation stats for chosen modes.
cipher preprocess path/to/perturbseq.h5ad -o out/ --modes log1p log1CP10k

# Forward prediction (transcriptomic shift) against null covariances.
cipher forward path/to/perturbseq.h5ad -o out/ --normalization log1p \
    --nulls meanfield shuffled

# Reverse prediction (recover the perturbed gene).
cipher reverse path/to/perturbseq.h5ad -o out/ --normalization log1p \
    --method matched_filter --top-k 10

# Condition-driver prediction (control vs condition, no ground truth).
cipher driver path/to/dataset.h5ad -o out/ \
    --condition-key stim --control rest --condition stim \
    --method matched_filter --top 25

Run cipher –version or cipher <command> –help for all options.