.. _package-info:

===============================
CIPHER package
===============================

Availability
-----------------
The package is available on PyPI as ``cipher-perturb`` (`PyPI page <https://pypi.org/project/cipher-perturb/>`_) or GitHub as ``cipher`` (`GitHub page <https://github.com/GoyalLab/CIPHER>`_).

How do I use it?
------------------

    #. Installation
        It can be installed from PyPI using the following:
        
        .. code-block:: python

            pip install cipher-perturb


        It can be installed from GitHub using the following:

        .. code-block:: python

            pip install -e .

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

            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:

        .. code-block:: python

            conda env create -f environment.yaml


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

    #. Forward Prediction: ``ForwardResult``

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

        .. code-block:: python

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

    #. Reverse Prediction: ``ReverseResult``
    
        This code identifies the causal perturbation from the Perturb-seq transcriptomic profile.

        .. code-block:: python

            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.

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

        .. code-block:: python

            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):

        .. code-block:: python

            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")

    #. 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.
        
        .. code-block:: python

            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",
            )

    #. Lower-level Dataset 
        
        For custom pipelines, load_dataset returns a Dataset that handles perturbation detection, target-gene mapping, and control/perturbation matrices.
        
        .. code-block:: python

            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).
	    
    .. code-block:: bash

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