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

Stage 1 scores every row of $Z$ (one perturbation) against a pathway collection, giving each perturbation a pathway-response profile;
Stage 2 ranks the perturbations along each pathway column of the Stage 1 matrix and scores, for every pathway $j$, whether perturbations of $j$’s members concentrate at either end.

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.

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.

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:

demo_z: the perturbation-by-gene numeric matrix of signed differential-expression statistics (rownames = perturbed gene symbols).
demo_pathways: a named list of character vectors of gene symbols.
stage1: overrides forwarded to stage1_nes(). Here nperm = 200 sets the per-perturbation permutation count (the default 1000 is more than this toy dataset needs) and min_size = 5 lowers the minimum effective gene-set size from its production default of 15, because the demo pathways have 30 members each.
stage2: overrides forwarded to stage2_matrix(). min_pu = 5 lowers the minimum number of pathway members that must be present among the perturbations (default 10) to match the small demo.
BPPARAM: any BiocParallel backend; a serial backend with a fixed RNGseed keeps this vignette deterministic. On real data use the default BiocParallel::bpparam().

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

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_raw.csv"         
#> [2] "M_matrix_demo_demo_fgsea.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).

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.

Session info

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> 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] fastmap_1.2.0       R.oo_1.27.1         jsonlite_2.0.0     
#> [13] Matrix_1.7-6        R.utils_2.13.0      BiocManager_1.30.27
#> [16] scales_1.4.0        codetools_0.2-20    jquerylib_0.1.4    
#> [19] cli_3.6.6           rlang_1.3.0         R.methodsS3_1.8.2  
#> [22] cowplot_1.2.0       withr_3.0.3         cachem_1.1.0       
#> [25] yaml_2.3.12         otel_0.2.0          tools_4.6.1        
#> [28] parallel_4.6.1      BiocParallel_1.47.0 dplyr_1.2.1        
#> [31] ggplot2_4.0.3       fastmatch_1.1-8     assertthat_0.2.1   
#> [34] buildtools_1.0.0    vctrs_0.7.3         R6_2.6.1           
#> [37] lifecycle_1.0.5     pkgconfig_2.0.3     bslib_0.12.0       
#> [40] pillar_1.11.1       gtable_0.3.6        data.table_1.18.6.1
#> [43] glue_1.8.1          Rcpp_1.1.2          xfun_0.61          
#> [46] tibble_3.3.1        tidyselect_1.2.1    sys_3.4.3          
#> [49] knitr_1.52          farver_2.1.2        msigdbr_26.1.1     
#> [52] htmltools_0.5.9     rmarkdown_2.32      maketools_1.3.2    
#> [55] labeling_0.4.3      compiler_4.6.1      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.