KAMP: Robust, scalable Ripley's K for spatial proteomics

This repository contains the code and analysis scripts to reproduce the results in:

Julia Wrobel* and Hoseung Song*.

A robust, scalable K-statistic for quantifying immune cell clustering in spatial proteomics data.

The method, KAMP (K adjustment by Analytical Moments of the Permutation distribution), corrects Ripley's K for spatial inhomogeneity in spatial proteomics samples. Instead of comparing the observed K to the homogeneous-CSR null (πr²), KAMP uses closed-form first and second moments of the permutation null distribution of Ripley's K, computed from the background cells in each sample. This gives:

• an empirical null E(K) that adjusts the degree-of-clustering statistic K̃(r) = K̂(r) − E(K(r)) for inhomogeneity,
• a within-sample test (Z statistic and approximate p-value) for excess clustering / colocalization, and
• major speedups over explicit permutation (no permutations are actually drawn — the moments are analytic), with a thinned variant, KAMP lite, for very large samples.

Both univariate (clustering of one cell type) and bivariate (colocalization of two cell types) versions are implemented.

The same method is also distributed as part of the spatialTIME and mxfda R packages.

---

Quick start: KAMP on a single sample

The core function is get_permutation_distribution() in source/get_permutation_distribution.R. It takes a spatstat point pattern (ppp) whose marks label cells as "immune" / "background" (univariate) or "immune1" / "immune2" / "background" (bivariate), and returns the KAMP statistics at a given radius — no cluster or external data needed.

library(spatstat.geom)
library(spatstat.random)
library(spatstat.explore)
library(tidyverse)

source("source/get_permutation_distribution.R")

# simulate one sample: immune + background cells in a unit window
set.seed(1)
W  <- owin(c(0, 1), c(0, 1))
pp <- rmpoispp(c(200, 1800), types = c("immune", "background"), win = W)

# univariate KAMP at a single radius
get_permutation_distribution(pp, rvalue = 0.05)

# across several radii
purrr::map_dfr(c(0.05, 0.1, 0.15),
               get_permutation_distribution, ppp_obj = pp)

# bivariate colocalization (marks must be "immune1"/"immune2"/"background")
# get_permutation_distribution(pp_biv, rvalue = 0.05, bivariate = TRUE)

The returned tibble has one row per radius with columns:

Column	Meaning
r	Radius.
khat	Observed Ripley's K (translation edge correction).
expectation	Permutation-null mean of K, i.e. the empirical null E(K(r)).
var	Permutation-null variance of K.
Z	Test statistic (khat − expectation) / sqrt(var).
pvalue	One-sided approximate p-value for excess clustering/colocalization.

The degree of clustering used in downstream models is K̃(r) = khat − expectation.

---

Repository structure

permuted_k/
├── source/                         # methods + simulation/analysis drivers
│   ├── get_permutation_distribution.R   # *** core KAMP function (analytic moments) ***
│   ├── utils_k.R                        # univariate drivers: get_k(), get_kamplite(),
│   │                                    #   get_k_power(), get_k_power_permOnly(), get_coxPH()
│   ├── utils_k_bivariate.R              # bivariate versions of the above
│   ├── utils.R                          # result-merging helpers + get_coxph() for real data
│   ├── simulate_ppp.R                   # homogeneous/inhomogeneous multitype Poisson sims
│   ├── simulate_scSpatialSim.R          # clustered sims with "holes" (via scSpatialSIM)
│   ├── vectraPolarisdata.R / .sh        # real ovarian data: download, process, run KAMP (univariate)
│   ├── vectraPolarisdata_biv.R / .sh    # real data: B cell / macrophage colocalization (bivariate)
│   └── simulations/                     # SLURM array jobs that populate output/
│       ├── k_univariate_expectation.R / .sh
│       ├── k_univariate_variance.R / .sh
│       ├── k_univariate_survival.R / .sh
│       ├── kampLite_univariate.R / .sh
│       ├── k_bivariate_expectation.R / .sh
│       └── k_bivariate_variance.R / .sh
├── analysis/                       # R Markdown that reads output/ and makes paper figures
│   ├── k_expectation.Rmd                # Fig 3 (degree of clustering) + Fig 4 (timing)
│   ├── k_power.Rmd                      # Fig 5 (Type I error & power)
│   ├── k_survival.Rmd                   # survival-bias simulation
│   ├── k_expectation_biv.Rmd            # bivariate expectation (Supplement C)
│   ├── k_power_biv.Rmd                  # bivariate power (Supplement C)
│   ├── kamplite_univariate.Rmd          # KAMP lite thinning (Supplement D)
│   └── vpData_analysis.Rmd              # Figs 1, 6, 7, 8 (real ovarian data analysis)
├── output/                         # simulation results + figures  (git-ignored; regenerated)
└── data/                           # processed data objects         (git-ignored; regenerated)

data/ and output/ are intentionally not tracked in git. You regenerate data/ from the public ovarian dataset (below) and populate output/ by running the simulation and real-data scripts.

---

Dependencies

Methods and simulations use the spatstat family (spatstat.geom, spatstat.random, spatstat.explore), the tidyverse, tictoc, here, survival, broom, and MASS. Simulated clustering with holes uses scSpatialSIM. Figures additionally use patchwork, viridis, ggbreak, ggstance, latex2exp, RESI, and mxfda.

The real data come from the Bioconductor package VectraPolarisData (HumanOvarianCancerVP()).

install.packages(c(
  "spatstat.geom", "spatstat.random", "spatstat.explore", "tidyverse",
  "tictoc", "here", "survival", "broom", "MASS", "scSpatialSIM",
  "patchwork", "viridis", "ggbreak", "ggstance", "latex2exp", "RESI", "mxfda"
))

if (!require("BiocManager")) install.packages("BiocManager")
BiocManager::install("VectraPolarisData")

---

Reproducing the paper

The pipeline has three stages: (1) prepare the data, (2) run the simulations / real-data jobs to fill output/, and (3) knit the analysis .Rmd files to build the figures.

1. Prepare the ovarian data

Running source/vectraPolarisdata.R interactively (it auto-detects a local session) downloads HumanOvarianCancerVP(), processes it, and saves data/processed_ovarian_data.Rda. The survival simulation also expects data/densities_ovarian_data.Rda (densities of KAMP / K−KAMP values derived from the real data), which is created during the real-data analysis.

2. Run the simulations and real-data jobs

Each driver in source/simulations/ is a SLURM array job: scenario i (the array index) selects one row of the parameter grid. The .sh wrappers submit the full grid; results are saved (date-stamped) under output/.

Driver (source/simulations/)	Array size	Output folder	Paper result
k_univariate_expectation.R	1–48	output/univariate_expectation/	Fig 3, Fig 4
k_univariate_variance.R	1–128	output/univariate_variance/	Fig 5
k_univariate_survival.R	1–240	output/univariate_survival/	survival sim
kampLite_univariate.R	1–30	output/kamplite/	Supplement D
k_bivariate_expectation.R	1–48	output/bivariate_expectation/	Supplement C
k_bivariate_variance.R	1–480	output/bivariate_variance/	Supplement C

Real data (in source/):

Driver	Output folder	Paper result
vectraPolarisdata.R / .sh	output/vpData/	Figs 1, 6, 7 (univariate within-sample)
vectraPolarisdata_biv.R / .sh	output/vpData_bivariate/	Figs 7, 8 (bivariate + survival)

Submit a job on a SLURM cluster with, e.g.:

cd source/simulations
sbatch k_univariate_expectation.sh

To run a single scenario locally instead, pass the scenario index directly:

Rscript source/simulations/k_univariate_expectation.R 3

(The driver scripts treat a working directory under /Users as a local run and shrink the workload automatically.)

3. Build the figures

Once output/ is populated, knit the analysis files. Each one reads the relevant output/ subfolder via the merge_files*() helpers in source/utils.R and writes figures to output/:

Analysis file (analysis/)	Reads	Produces
k_expectation.Rmd	output/univariate_expectation/	Fig 3, Fig 4
k_power.Rmd	output/univariate_variance/	Fig 5
k_survival.Rmd	output/univariate_survival/	survival-bias figure
k_expectation_biv.Rmd	output/bivariate_expectation/	Supplement C
k_power_biv.Rmd	output/bivariate_variance/	Supplement C
kamplite_univariate.Rmd	output/kamplite/	Supplement D
vpData_analysis.Rmd	data/processed_ovarian_data.Rda, output/vpData/, output/vpData_bivariate/	Figs 1, 6, 7, 8

---

Notes for reproduction

A few things to keep in mind when running end to end:

• Fresh output/ per run. The analysis files read all result files under each output/ subfolder (recursively). The simulation drivers save results into date-stamped subfolders, so this picks them up automatically regardless of the run date — but if you run the same simulation more than once, clear the old results first so runs are not mixed together.
• Regenerate data/ and output/. Both folders are git-ignored. Build data/processed_ovarian_data.Rda from VectraPolarisData (Section 1) and populate output/ by running the simulation and real-data scripts (Section 2) before knitting the analysis files (Section 3).
• Cluster specifics. The .sh scripts use SLURM directives and site-specific partitions/modules (module load R, partitions encore / wrobel). Adjust these for your environment, or run scenarios locally as shown above.

---

Acknowledgments

Hoseung Song is supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) RS-2025-16066571.