Contents

1 What path2path computes

A genome-scale Perturb-seq screen is, in effect, thousands of parallel differential-expression experiments: each perturbed gene yields a signed per-gene statistic across all measured genes, and the results stack into a perturbation-by-gene matrix \(Z\). path2path turns this matrix into a directional pathway-by-pathway interaction map \(\mathbf{M}\) by running preranked GSEA twice (Subramanian et al. 2005; Korotkevich et al. 2019):

Entry \(\mathbf{M}_{jk} > 0\) means that perturbations targeting members of pathway \(j\) are enriched toward higher Stage 1 scores for pathway \(k\). In a loss-of-function screen, when those scores reflect decreased expression of pathway \(k\) after perturbation, a positive entry suggests that \(j\) supports or activates the transcriptional program represented by \(k\); a negative entry suggests suppression. The map therefore identifies candidate activating and suppressive relationships. It is asymmetric because \(\mathbf{M}_{jk}\) and \(\mathbf{M}_{kj}\) summarize different perturbation-readout combinations.

Three safeguards are built in: the perturbed gene is excluded from its own ranked list (removing an on-target tautology that otherwise dominates the self-interaction diagonal), effective-size filters are applied in both stages, and entries whose pathway pair overlaps at Jaccard \(> 0.5\) are masked.

1.1 Relation to the Bioconductor ecosystem

path2path builds directly on Bioconductor infrastructure: enrichment statistics are computed by fgsea, pathway collections are retrieved through msigdbr, parallelization uses BiocParallel backends throughout, and inputs and results interconvert with SummarizedExperiment objects (as_path2path_matrix(), as_path2path_se()), so maps drop into any SummarizedExperiment-based workflow.

Its scope differs from existing enrichment packages. Tools such as fgsea, clusterProfiler, or GSVA answer the question “which pathways respond in this contrast (or sample)?” — one enrichment analysis per ranked list, with pathways as the output rows. path2path asks the next question that a genome-scale Perturb-seq screen makes answerable: “which pathway, when perturbed, moves which other pathway, and in which direction?” Its output is not a pathway list but an asymmetric pathway-by-pathway matrix, where rows are perturbed pathways and columns are transcriptional readouts. To our knowledge no current Bioconductor package derives such directional pathway-level interaction maps from Perturb-seq data; the closest workflows (per-perturbation enrichment followed by manual aggregation) correspond to running only Stage 1 of this package.

2 A five-minute tour on synthetic data

The package ships a small simulated dataset with two planted relations: PATHWAY_A activates PATHWAY_B, and PATHWAY_C suppresses PATHWAY_D.

library(path2path)
data(demo_z)
data(demo_pathways)
dim(demo_z)
#> [1] 120 400

Build the map. The arguments:

res <- path2path(demo_z, demo_pathways,
    stage1 = list(nperm = 200, min_size = 5),
    stage2 = list(min_pu = 5),
    BPPARAM = BiocParallel::SerialParam(RNGseed = 1))
round(res$M, 2)
#>           PATHWAY_A PATHWAY_B PATHWAY_C PATHWAY_D PATHWAY_E PATHWAY_F
#> PATHWAY_A     -1.28      3.46     -1.65      2.14     -1.40     -1.74
#> PATHWAY_B        NA        NA        NA        NA        NA        NA
#> PATHWAY_C      1.50     -1.85      1.36     -3.45      1.73      2.05
#> PATHWAY_D        NA        NA        NA        NA        NA        NA
#> PATHWAY_E     -0.56     -2.28     -0.55      2.00      1.08      0.80
#> PATHWAY_F        NA        NA        NA        NA        NA        NA

The result is a list whose central component M is the pathway-by-pathway map: entry \((j, k)\) is a Stage 2 normalized enrichment score, positive when perturbing members of pathway \(j\) (rows) coherently activates pathway \(k\) (columns), negative for suppression, and NA where the entry is untestable (too few of \(j\)’s members among the perturbations) or masked by the Jaccard overlap filter — in the printed matrix above, the testable rows are exactly A, C, and E, the three pathways whose members were perturbed. Alongside M, the list carries the Stage 1 matrices (nes, pval_stage1), the unmasked map (M_raw), the per-entry nominal p-values (pval), and their Benjamini–Hochberg adjustment over all populated entries (padj):

res$M["PATHWAY_A", "PATHWAY_B"]   # planted activation, > 0
#> [1] 3.463951
res$M["PATHWAY_C", "PATHWAY_D"]   # planted suppression, < 0
#> [1] -3.454238
res$padj["PATHWAY_A", "PATHWAY_B"]
#> [1] 4.332938e-28
res$padj["PATHWAY_C", "PATHWAY_D"]
#> [1] 4.332938e-28

Both planted relations are recovered with the right sign and direction, as the two strongest entries of the map. The reverse entries (B acting on A, D on C) are untestable because B and D members were not perturbed, so those rows fall below the Stage 2 universe filter. Note that in a toy map this small (18 testable entries), the two planted blocks also leave weaker significant echoes in other entries: once A’s members occupy one extreme of a column’s ranking, the remaining pathways’ members are displaced from uniformity too. On realistically sized collections this rank-compression effect is diluted across thousands of perturbations; here it is a useful reminder that entries of \(\mathbf{M}\) are not independent tests, which is exactly why the package’s inference treats the populated map as one BH family and offers whole-pipeline permutation (Route B, below) for entries that will be individually reported.

Plot the map with significance dots:

plot_map(res, padj_cutoff = 0.05)

And inspect why one entry is strong, at the level of the underlying ranking:

plot_entry(res, demo_pathways, "PATHWAY_A", "PATHWAY_B")

3 Real data

For a real screen, the input is the perturbation-by-gene matrix of differential-expression statistics (a DESeq2 Wald \(z\), a \(z\)-normalized pseudobulk difference, or a signed log-fold change), with rows named by the perturbed gene symbols. Such matrices are typically stored as (gzipped) CSV with the perturbation names in the first column; read_input_matrix() reads that format. The package ships the demo matrix in this on-disk form, which stands in for your own file below:

f <- system.file("extdata", "zscore_demo.csv.gz", package = "path2path")
z <- read_input_matrix(f)  # your screen: "zscore_myscreen.csv.gz"
all.equal(z, demo_z)       # the demo matrix, up to CSV round-trip precision
#> [1] TRUE
write_result(res, dir = file.path(tempdir(), "results"),
    tag = "demo", collection = "demo")
list.files(file.path(tempdir(), "results"))
#> [1] "M_matrix_demo_demo_fgsea.csv"             
#> [2] "M_matrix_demo_demo_fgsea_raw.csv"         
#> [3] "nes_matrix_demo_demo_fgsea.csv.gz"        
#> [4] "padj_matrix_demo_demo_fgsea.csv"          
#> [5] "pval_matrix_demo_demo_fgsea.csv"          
#> [6] "pval_matrix_stage1_demo_demo_fgsea.csv.gz"

On real data you would draw pathways from MSigDB instead of the demo list — load_pathways("hallmark") (or "reactome", "kegg_legacy", …) downloads the collection through msigdbr — and run path2path(z, pathways) with the production defaults. write_result() writes each result component (M, pval, padj, the Stage 1 matrices) as a CSV named <component>_matrix_<tag>_<collection>_<method>.csv, as listed above.

If the statistics are stored in a SummarizedExperiment assay with genes in rows and perturbations in columns, use as_path2path_matrix() to transpose and validate the assay. as_path2path_se() converts a completed map back to a SummarizedExperiment whose map and p-value matrices are assays.

Reproducibility: BiocParallel backends draw worker seeds from entropy unless given one explicitly, so set.seed() alone is not enough — pass BPPARAM = MulticoreParam(RNGseed = 1) (or SerialParam(RNGseed = 1)) to make a run exactly repeatable.

Orientation: statistics are negated internally for loss-of-function screens (CRISPRi/knockout, the default); set stage1 = list(negate = FALSE) for gain-of-function screens (CRISPRa, overexpression).

4 Two routes to significance

path2path() attaches Route A p-values: the Stage 2 nominal p-values (fgsea::fgseaMultilevel) with one Benjamini-Hochberg family over all populated entries. On the datasets analyzed in the accompanying manuscript, the realized false-discovery rate of BH screening at 0.05 remained close to the nominal level across datasets and pathway collections, and the null distribution of the nominal p-values was close to uniform over the range that could be assessed by permutation.

For individually examined strong entries, Route B reruns the full pipeline on row/column-permuted inputs, pools the null \(|\mathbf{M}'|\) values, and extrapolates the extreme tail with a generalized-Pareto (GPD) fit (Knijnenburg et al. 2009). The arguments: res is the observed path2path() result (its stored $config is reused so the null runs repeat the observed configuration exactly), demo_z/demo_pathways are the inputs that produced it, B is the number of whole-pipeline permutations (B = 50 is the function default and takes a minute or two on the demo data; production runs in the accompanying manuscript used B = 20–50 on maps with two to nearly four orders of magnitude more populated entries), and seed_base fixes the permutation in null run \(b\) to seed_base + b, making the calibration reproducible (the GSEA side is already seeded here through the RNGseed backend stored in res$config). This is the one deliberately expensive chunk of the vignette — Route B’s point is to pay \(B\) full pipeline runs for a null distribution of the entire construction:

cal <- calibrate_pvalues(res, demo_z, demo_pathways, B = 50, seed_base = 42)
#> calibrate_pvalues: reusing the non-default configuration in `result` for the null runs
unlist(cal$gpd[c("u", "xi", "beta", "n_exc", "q_used")])
#>           u          xi        beta       n_exc      q_used 
#>  1.57453648  0.04866889  0.17344878 45.00000000  0.95000000
cal$pval_gpd["PATHWAY_A", "PATHWAY_B"]
#> [1] 8.001607e-06
cal$pval_pool["PATHWAY_A", "PATHWAY_B"]
#> [1] 0.001109878

Reading the output: cal$gpd reports the tail fit — threshold u (the q_used quantile of the pooled null), shape xi and scale beta of the GPD fitted to the n_exc exceedances. cal$pval_pool is the pooled empirical p-value, floored at \(1/(L+1)\) for a pool of \(L\) null entries; cal$pval_gpd equals it below the threshold and replaces it with the GPD tail extrapolation above, which is why the planted A\(\to\)B entry — stronger than every one of the 900 pooled null values — obtains a p-value well below the pooled floor. cal$padj_gpd applies one BH family over the populated entries of pval_gpd. The calibration plot compares the nominal (Route A) and calibrated (Route B) p-values entry by entry:

plot_calibration(cal)

For long production runs, null_dir persists each completed null map so interrupted runs resume, and the returned cal$pool_abs can be passed back via pool = to recalibrate without recomputation. Route A provides the computationally lighter default for routine screening, while Route B carries uncertainty from both stages through the permutation.

5 Session info

sessionInfo()
#> R version 4.6.1 Patched (2026-06-24 r90190)
#> Platform: x86_64-apple-darwin20
#> Running under: macOS Ventura 13.7.8
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.6-x86_64/Resources/lib/libRblas.0.dylib 
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.6-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
#> 
#> locale:
#> [1] C/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#> 
#> time zone: America/New_York
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] path2path_0.99.4 BiocStyle_2.41.0
#> 
#> loaded via a namespace (and not attached):
#>  [1] sass_0.4.10         generics_0.1.4      fgsea_1.39.2       
#>  [4] lattice_0.23-1      digest_0.6.39       magrittr_2.0.5     
#>  [7] evaluate_1.0.5      grid_4.6.1          RColorBrewer_1.1-3 
#> [10] bookdown_0.48       fastmap_1.2.0       R.oo_1.27.1        
#> [13] jsonlite_2.0.0      Matrix_1.7-6        R.utils_2.13.0     
#> [16] tinytex_0.61        BiocManager_1.30.27 scales_1.4.0       
#> [19] codetools_0.2-20    jquerylib_0.1.4     cli_3.6.6          
#> [22] rlang_1.3.0         R.methodsS3_1.8.2   cowplot_1.2.0      
#> [25] withr_3.0.3         cachem_1.1.0        yaml_2.3.12        
#> [28] otel_0.2.0          tools_4.6.1         parallel_4.6.1     
#> [31] BiocParallel_1.47.0 dplyr_1.2.1         ggplot2_4.0.3      
#> [34] fastmatch_1.1-8     assertthat_0.2.1    vctrs_0.7.3        
#> [37] R6_2.6.1            lifecycle_1.0.5     magick_2.9.1       
#> [40] pkgconfig_2.0.3     bslib_0.12.0        pillar_1.11.1      
#> [43] gtable_0.3.6        data.table_1.18.6.1 glue_1.8.1         
#> [46] Rcpp_1.1.2          xfun_0.62           tibble_3.3.1       
#> [49] tidyselect_1.2.1    knitr_1.52          dichromat_2.0-1    
#> [52] farver_2.1.2        msigdbr_26.1.1      htmltools_0.5.9    
#> [55] labeling_0.4.3      rmarkdown_2.32      compiler_4.6.1     
#> [58] S7_0.2.2

References

Knijnenburg, Theo A., Lodewyk F. A. Wessels, Marcel J. T. Reinders, and Ilya Shmulevich. 2009. “Fewer Permutations, More Accurate p-Values.” Bioinformatics 25 (12): i161–68. https://doi.org/10.1093/bioinformatics/btp211.
Korotkevich, Gennady, Vladimir Sukhov, and Alexey Sergushichev. 2019. “Fast Gene Set Enrichment Analysis.” bioRxiv, ahead of print. https://doi.org/10.1101/060012.
Subramanian, Aravind et al. 2005. “Gene Set Enrichment Analysis: A Knowledge-Based Approach for Interpreting Genome-Wide Expression Profiles.” Proceedings of the National Academy of Sciences 102 (43): 15545–50. https://doi.org/10.1073/pnas.0506580102.