py-Pando

Pure-Python Pando: gene regulatory network inference from paired single-cell RNA + ATAC data, AnnData-native, parity-tested against the R package.

Pando (Fleck et al., Nature 2023, doi:10.1038/s41586-022-05459-w) fits, for every target gene, a regression of its expression on TF expression x accessibility of TF-binding peaks near the gene, and reads significant TF–peak terms as regulatory links. py-Pando ports it following the omicverse-rebuildr protocol: the R source is the specification, and pre-registered parity gates (data/manifest.yaml) were committed before any algorithmic code.

Install

pip install pypando            # core
pip install "pypando[xgb]"     # + xgboost models
pip install "pypando[umap]"    # + UMAP layout of the module graph

Quick start

import pypando as pp

# rna:  AnnData, cells x genes, log-normalised (Seurat LogNormalize)
# atac: AnnData, cells x peaks, TF-IDF (Signac RunTFIDF), var_names "chr-start-end"
annotation = pp.annotation_from_gtf("gencode.v41.basic.annotation.gff3.gz")

grn = pp.initiate_grn(rna, atac, annotation)                # peaks minus exons
grn = pp.find_motifs(grn, genome="hg38.fa")                 # Pando's 1,590 motifs, MOODS
grn = pp.infer_grn(grn, genes=targets, method="glm")        # one model per target gene
grn = pp.find_modules(grn)                                  # TF modules

pp.coef(grn)            # TF / target / region / estimate / std_err / pval / padj / corr
pp.gof(grn)             # R^2 and number of variables per target gene
pp.NetworkModules(grn)  # module table + genes_pos / genes_neg / regions / peaks

Within Epione: epione.grn.pando(rna, atac, genome="hg38", genes=targets).

Function map

Full numbers: RECONSTRUCTION_REPORT.md and examples/compare_R_vs_Python.ipynb.

What makes the glm network identical

R's glm uses LINPACK's dqrdc2 QR with limited pivoting: a column whose remaining norm falls below 1e-11 of its original norm is moved to the end and its coefficient becomes NA. Pando's per-gene designs (hundreds to thousands of peak x TF products, often nearly collinear) sit right at that threshold, so the pivot decisions depend on the last bits of BLAS level-1 results. py-Pando's QR (pypando._blas, numba) calls the same BLAS routines (dnrm2, ddot, daxpy, dscal) in the same order as R's Fortran, multiplies with dgemv like R's %*%, sums in long double like R's sum(), and inverts with LAPACK dpotri like chol2inv. With the same BLAS build (here OpenBLAS 0.3.34) and one BLAS thread (infer_grn(blas_threads=1), the default), the glm network is bit-identical to R's, including genes where R's IRLS does not converge.

Notes

• cv.glmnet folds come from an emulation of R's set.seed() + sample() (pypando.rrng), so infer_grn(method="cv.glmnet", seed=42) uses R's folds after set.seed(42). Fold fits follow glmnet 5.0: each fold on its own lambda path, predictions interpolated onto the full path (glmnet 4.x refitted folds on the full path; cv_glmnet_gaussian(fold_lambda="full")).
• glmnet: py-Pando's coordinate descent follows glmnet's path rules (standardisation, lambda sequence, fdev / devmax stopping) but not its strong-rule screening, so for genes with hundreds of variables both stop within glmnet's 1e-7 tolerance at slightly different coefficients. fit_glmnet ports family = "gaussian" given as a string, which is how infer_grn calls it; R's fit_glmnet called directly with its default family = gaussian (a function) goes to glmnet's generic glmnet.path solver, which is not ported.
• xgb: Pando calls xgboost::xgboost(data, label, params = list(max_depth = 3, eta = 0.01), ...). xgboost >= 3 changed that function's interface and silently ignores params, training with depth 6 and learning rate 0.3. py-Pando uses Pando's parameters; with them honoured in R (xgb.train), importances match to 1e-15.
• bagging_ridge draws from NumPy's global RNG, as Pando's scikit-learn call does from R; infer_grn(seed=...) seeds it.
• New in Python (not in R): region_background, regions_near_genes and subset_regions scan only the regions a target gene can reach, against the background of all regions — the network is identical to a full scan.

Reproducing the R comparison

python tests/make_fixture.py pbmc_granulocyte_sorted_10k_filtered_feature_bc_matrix.h5 data/fixture
Rscript tests/r_reference_driver.R data/fixture ref                     # R 4.4.3, Pando 1.1.1
python tests/_run_candidate.py ref cand hg38.fa
python tests/parity_check.py ref cand results.json
PYPANDO_REF_DIR=ref PYPANDO_FASTA=hg38.fa PYPANDO_CAND_DIR=cand pytest tests/test_exact_match.py

Citation

Fleck, J.S., Jansen, S.M.J., Wollny, D. et al. Inferring and perturbing cell fate regulomes in human brain organoids. Nature 621, 365–372 (2023).

License

MIT, as the upstream R package (see LICENSE).