PI / Ti Clinical Association Index (cell-subset ratio vs clinical change)

Overview

Both indices reduce a cell-subset frequency (ratio) to a single signed number that
measures how strongly that subset tracks a clinical score change. The statistic is
index = (-slope / |slope|) * R2 from lm(clin_value ~ ratio):
the sign of the slope (a decrease in the subset with a better clinical outcome
gives a positive index) multiplied by the regression R2 (effect strength).
Range is [-1, 1]; 0 and NA slopes become NA.

PI is horizontal (between patients at baseline); Ti is longitudinal (within patient
across treatment). Use PI when only baseline samples exist, Ti when both time points
are available for enough patients.

Canonical script

• scripts/pi_ti_clinical_index.R.txt — CLI-parameterized version of the author's
  original PI_TI.R. Statistical definition is byte-identical in behaviour
  (verified: base-R recomputation on mock data reproduced all slope / R2 / Pi
  values exactly).
• Original (read-only, do not modify): the author's PI_TI.R, with README.md in the same directory.
• The author's script hardcoded absolute project paths and the input rds path; those paths do not exist in this container.

Inputs

An .rds of a Seurat object (only meta.data is used), or an .rds holding a plain
data.frame of meta.data (useful for testing). One row per cell. Required columns
are checked explicitly; a missing one aborts with the list of available columns.

The clinical change column is usually derived, e.g. the author commented
meta$NPS.change <- meta$NPS.post - meta$NPS.pre. The script never computes it —
add it to meta.data yourself before running.

Steps

PI (pre-treatment, cross-sectional)
1. Subset meta to time_col == pre_level (aborts if the level is absent, listing the
observed levels).
2. group_by(sample_col, cluster_col) → n(); then per sample ratio = n / sum(n)
→ one row per sample × subset.
3. Build a per-sample clinical table (distinct(sample_col, .keep_all = TRUE) →
sample, clin_value = clin_col_pi, treatment = group_col) and left_join it
onto the ratio table.
4. Per treatment × subCluster1: keep only n() >= min_n and sd(ratio) > 0.
5. lm(clin_value ~ ratio) → slope, R2, Pi = (-slope/|slope|) * R2;
sort by treatment, desc(Pi).
6. One lollipop panel per treatment, arranged side by side.

Ti (paired pre/post, longitudinal)
1. group_by(sample_col_ti, time_col, cluster_col) → n(); ratio = n / sum(n)
within each patient × time point.
2. Join the same per-patient clinical table (by patient).
3. pivot_wider(names_from = time, values_from = ratio, names_prefix = "ratio_")
→ columns ratio_pre / ratio_post (names follow --pre-level / --post-level).
4. delta_ratio = ratio_post - ratio_pre; drop rows with NA delta — i.e. patients
with only one time point (the pairing requirement; count is logged).
5. Per treatment × subCluster1: keep only n() >= min_n and sd(delta_ratio) > 0.
6. lm(clin_value ~ delta_ratio) → slope, R2, Ti = (-slope/|slope|) * R2;
sort by treatment, desc(Ti).
7. One lollipop panel per treatment.

--mode PI|Ti|both selects which blocks run; both share the same palette and output
directory.

Usage

S=scripts/pi_ti_clinical_index.R.txt
Rscript $S --help

# author's own design (CD4T object: PI on baseline samples, Ti on paired patients)
Rscript $S \
  --input-rds /path/CD4T.seu.globalC.rds \
  --output-dir ./PI_TI_results \
  --time-col time --pre-level pre --post-level post \
  --group-col treatment \
  --sample-col sample --sample-col-ti patient \
  --cluster-col subCluster1 \
  --clin-col-pi NPS.change --clin-col-ti SNOT.change \
  --mode both --min-n 3

Parameters

Built-in palette (the author's CD4 map, used when --palette-csv is absent) — the
names must be edited for other datasets:

CD4_Tcm_IL7R #9AA7B1   CD4_Tem_CCL4 #BA5B49    CD4_Th17_CCR6 #D3CCD6
CD4_Th2_IL13 #a382bd   CD4_Tn_CCR7  #BC836B    CD4_Treg_FOXP3 #B1D5C8
CD4_Treg_TNFRSF4 #EE7959   CD4_Ttrm_TOX2 #E5A84B

Outputs

<output-dir>/

The lollipop geometry is the author's: grey stem from 0 (geom_segment, grey70,
linewidth 0.8) plus a grey90 halo point (shape 21, size 6, stroke 1.2) and the filled
point on top (shape 21, size 5, stroke 0.3), theme_classic(base_size = 14), x labels
at 45°, no legend. Subsets are ordered by descending index, so sign and magnitude are
read straight off the y-axis.

Pitfalls

• Column names must exist in meta.data exactly. The script aborts with the full
  available-column list rather than guessing; check colnames(obj@meta.data) first.
• library(SingleCellExperiment) silently masks dplyr::rename. S4Vectors
  exports its own rename, so a later library(SingleCellExperiment) breaks
  rename(sample = !!col) with 参数类型无效 / invalid argument type (base ! gets
  evaluated instead of tidyselect injection). Fix in the script: re-attach dplyr
  after all libraries and call dplyr::rename explicitly. The author's script was
  immune only because its library(dplyr) line came after
  library(SingleCellExperiment). Same class of trap: matrixStats::count,
  IRanges::slice (ggpubr::mutate is identical to dplyr::mutate, so it is safe).
• n() >= 3 and sd(x) > 0 silently drop whole groups. These filters remove
  any treatment × subset combination with fewer than 3 samples (PI: samples; Ti:
  paired patients) and any combination where the ratio/delta is constant (R2 is
  undefined). Small subsets therefore vanish entirely from the table and the plot —
  the run log prints before -> after counts per block, so check it before
  concluding a subset has no signal. Loosening --min-n trades coverage for
  unreliable slopes.
• Ti needs true pairing. delta_ratio is NA whenever a patient has only one
  time point (or one of the two, e.g. post missing) and those rows are dropped —
  a patient contributes to Ti only if both --pre-level and --post-level rows
  exist for the same --sample-col-ti value. The ratio is also computed per patient
  × time point, so an unbalanced pre/post cell count is fine.
• Palette names must match your subset names. scale_fill_manual(values = pal, limits = names(pal)) means any subset absent from the palette is set to NA and
  the point is dropped from the panel (ggplot warns "Removed k rows…"). The script
  pre-empts this with its own warning listing the offending subsets; pass
  --palette-csv (or edit default_palette) with your real subset names.
• PI and Ti use different sample/clinical columns in the author's design
  (sample + NPS.change vs patient + SNOT.change). Keep --sample-col-ti and
  the two --clin-col-* flags straight; --clin-col is a shortcut that sets both.
• --sample-col-ti must be constant within a patient across pre/post, otherwise the
  pairing breaks and every delta_ratio is NA.
• The clinical column has to be the change score. Passing a raw value silently
  produces a meaningless index (no error) — PI/Ti assume clin_value already is a
  delta or a change-from-baseline score.
• The library(Matrix/scales/matrixStats/slingshot/SingleCellExperiment) lines are
  carried over from the author's script; only Seurat, tidyverse/dplyr/tidyr,
  ggplot2, ggpubr, RColorBrewer are actually needed. The optional ones are
  soft-checked (a message instead of an error if missing) so the script still runs.

Environment

• Author's environment (README): a Singularity/Apptainer base image,
  conda activate R.4.4.3, start-rstudio.
• Abio container: Rscript is /opt/conda/bin/Rscript (R 4.5.3). .libPaths() =
  /opt/user_R_libs → /opt/public_R_libs → /opt/conda/lib/R/library.
  The env exports a stale R_HOME=/usr/local/lib/R; R 4.5.3 ignores it and prints
  WARNING: ignoring environment value of R_HOME — harmless.
• optparse is NOT installed in the Abio container, so the script parses
  commandArgs(trailingOnly = TRUE) in base R (--flag value, --flag = TRUE).
  Do not add an optparse dependency unless the user asks for it.
• Verified present (container, R 4.5.3): Matrix, scales, matrixStats,
  slingshot, Seurat, SingleCellExperiment, tidyverse, RColorBrewer,
  ggplot2, dplyr, ggpubr, tidyr.
• Write outputs under <output-dir>/{Img,Stat,Rdata} (the Abio convention), and keep
  anything user-visible in the container /output tree so the frontend can preview it.
• Validation actually run (not the real object): a synthetic meta.data with 2
  treatments × 5 subsets × 8 samples (one unpaired post sample, one 2-patient arm)
  went through --mode both — 15 → 10 groups kept, 1 unpaired patient dropped,
  PDFs written, and all slope / R2 / Pi values matched an independent base-R
  recomputation exactly. Rscript --help exits 0 and
  invisible(parse(file = ...)) is clean.

> Packaging note. Scripts and tables in this skill carry a .txt suffix on their
> original extension (e.g. .R.txt, .tsv.txt) because the SkillHub upload allowlist
> does not accept .R or .tsv. This is cosmetic: Rscript is extension-agnostic, so
> Rscript scripts/<name>.R.txt ... behaves exactly like .R, and tab-separated data
> parses identically from .tsv.txt. Keep the suffix when repackaging.