A gradient-aware application of phenotype-to-cell mapping in normal → peri-lesional → lesional skin
Psoriatic skin is usually studied as a two-state contrast (lesional vs. normal), which discards the clinically meaningful peri-lesional margin where disease is actively expanding. We reframe the problem as an ordinal gradient — uninvolved (NN) < peri-lesional (PN) < lesional (PP) — and ask which single cells co-vary with that gradient. Using Scissor, we project an ordinal bulk RNA-seq phenotype (93 biopsies spanning all three tiers, F=254 for the tier axis on bulk PC1) onto a 89,058-cell single-cell reference that contains the peri-lesional compartment. Because the phenotype label is a clinical biopsy site rather than a molecular signature derived from the same cells, the mapping is non-circular by construction. Selected cells are monotonic on the gradient (mean tier: Scissor− 0.79 < background 1.38 < Scissor+ 1.43) and the gradient-tracking (Scissor+) fraction peaks at the peri-lesional tier. Both a reliability test (p = 0.000) and a selection permutation null (p = 0.000) confirm the signal is driven by real phenotype structure. Endothelial cells dominate the gradient-tracking population (5.2× enriched, OR 11.3), and the associated 1,861-gene program is vascular-led; STAT3 is a significant, if modest, member (log2FC 0.43, padj 0.018). An orthogonal bulk deconvolution independently confirms the compositional trend for the three strongest lineages. This document reports the validated laptop-scale backbone; the full-census run and additional robustness tiers are staged for cluster execution.
Standard psoriasis transcriptomics contrasts lesional (PP) against uninvolved (NN) skin. That design is blind to the peri-lesional zone (PN) — the advancing margin a few millimetres outside the visible plaque, where the transition from health to disease is presumably underway. If a distinct cellular program initiates psoriatic conversion, the peri-lesional compartment is where it should be visible, and a two-state design cannot see it.
We therefore treat biopsy site as an ordinal phenotype (NN < PN < PP) and use Scissor (Sun, Guan, … Xia, Nat Biotechnol 2022) to identify the single cells whose expression co-varies with that gradient. Scissor correlates every single cell against a bulk cohort carrying a phenotype label, then solves a network-regularized regression to select the cells most consistent with the phenotype. Our central design choice is that the phenotype is a clinical biopsy-site label, not a molecular score computed from the reference cells — so a cell being flagged cannot be a restatement of how it was labelled. STAT3, a longstanding candidate in psoriasis, is examined as a member of the resulting program rather than as an assumed driver.
We assembled the Ma et al. (2023, Nat Commun) psoriasis atlas from 33 per-sample count matrices: 96,088 cells × 33,538 genes across 22 donors, spanning all three tiers (NN 13,534; PN 35,518; PP 47,036), with 11 donors contributing paired PN+PP biopsies — the design feature that makes the peri-lesional compartment usable.
After Ensembl→symbol mapping (24,185 symbols) and QC (nFeature 200–6000, nCount ≥500, percent-MT <20), 89,058 cells (92.7%) were retained. Standard normalization, HVG selection (2,000), PCA (30), neighbour graph (dims 1:20), Louvain clustering (resolution 0.5) and UMAP produced 23 clusters, annotated to nine broad lineages by canonical markers.
Among five recount3 psoriasis studies carrying tier labels, only SRP165679 (Tsoi et al. 2019) balances all three tiers at usable depth: NN=38, PN=27, PP=28 (93 classified samples). Using a single study removes cross-study confounding. edgeR logCPM (filterByExpr) and Ensembl→symbol collapse gave 24,533 symbols; the ordinal response is y ∈ {0,1,2}. Non-degeneracy checks pass: minimum tier size 27, and bulk PC1 (27.8% variance) tracks the tier axis at F=253.9, p ≈ 10−38.
| Check | Value |
|---|---|
| Tiers / min samples per tier | 3 / 27 |
| Bulk PC1 variance explained | 27.8% |
| Bulk PC1 ∼ tier | F=253.9, p≈10−38 |
| Cross-study confound | none (single study) |
| Phenotype circularity | none (clinical biopsy-site label) |
Scissor's selection step is a graph-regularized elastic net solved by a compiled routine (APML1). That component could not be built in our environment, so we reimplemented the identical objective in pure R on top of glmnet. The network penalty is the symmetric-normalized graph Laplacian, L_sym = I − D−1/2 A D−1/2; at ~105 cells a dense Laplacian is intractable, so we encode it as a sparse edge-difference augmentation — one sparse row per graph edge realizing βTLsymβ = Σedges(βi/√di − βj/√dj)2. Sparsity is controlled by walking the glmnet λ-path to a target selected fraction (the operational form of Scissor's cutoff). On synthetic two-community data the port recovered 100/100 true-positive cells (cor(β,βtrue)=0.72). This is a validated equivalent, not the authors' canonical solver; a cross-check against compiled Scissor is staged for the cluster.
For tractable graph regression the backbone runs on a stratified 20,023-cell subset (preserving cell-type×tier proportions; full census staged for cluster). Correlating 93 bulk samples against 20,023 cells over 1,664 shared HVGs and tuning the graph-smoothing parameter gave alpha=0.40, 14.84% of cells selected (1,574 Scissor+, 1,397 Scissor−), under the 20% cutoff.
The selected cells are monotonic on the gradient: mean tier (NN=0, PN=1, PP=2) rises Scissor− 0.79 < background 1.38 < Scissor+ 1.43. Critically, the gradient-tracking Scissor+ fraction peaks at the peri-lesional tier (9.4% of PN cells vs. 8.0% of PP), the direct signature of an intermediate progression state that a two-state design would miss.
Two independent controls confirm the result. The reliability test (100 label permutations) gives real cross-validated MSE 0.147 vs. a null mean of 0.779 (p = 0.000, 0/100 permutations lower). The selection permutation null (30 reruns with shuffled tier labels, full re-selection) collapses the positive-minus-negative tier gap from the real 0.638 to a null mean of −0.305 (p = 0.000). The selection reflects real phenotype structure, not graph geometry or overfitting.
Endothelial cells overwhelmingly dominate the gradient-tracking fraction, with fibroblasts and melanocytes marking the opposite (normal-associated) pole.
| Cell type | Fold in Scissor+ | Odds ratio | p | Interpretation |
|---|---|---|---|---|
| Endothelial | 5.18× | 11.26 | 6×10−245 | lesional-tracking |
| DC | 1.40× | 1.47 | 0.014 | weakly lesional |
| Keratinocyte | 0.87× | 0.69 | 8×10−12 | near-uniform |
| NK | 0.40× | 0.36 | 2×10−12 | depleted from Scissor+ |
| Fibroblast | 0.18× | 0.15 | 5×10−53 | normal-tracking |
| Melanocyte | 0.04× | 0.04 | 3×10−20 | normal-tracking |
Differential expression of Scissor+ vs. background (Wilcoxon) yields 1,861 genes at padj < 0.05. The top up-regulated genes are vascular/endothelial — CCL14, ACKR1, RAMP3, PLVAP, APLNR, CYTL1, SPNS2 — consistent with dermal angiogenesis as a hallmark of psoriatic progression.
STAT3 is a significant member of the program (log2FC 0.43, padj 0.018; 48.5% vs. 44.7% of cells expressing). The direction matches the STAT3 hypothesis, but the effect is modest and the program is not STAT3-led — the dominant axis is vascular.
As an independent check by a different method, we deconvolved the real SRP165679 bulk (NNLS against a Ma-reference signature) and tested each cell type for a monotonic proportion trend across NN→PN→PP. Scissor per-cell direction is taken from the Fisher enrichment among Scissor+ cells; a cell type is concordant only when its bulk proportion trend is significant and matches.
| Cell type | Scissor direction | Bulk proportion trend | padj | Status |
|---|---|---|---|---|
| Endothelial | lesional-tracking | rises 0.0→0.2→3.8% | 1×10−19 | ✓ concordant |
| Fibroblast | normal-tracking | falls 12.3→8.2→1.4% | 1×10−11 | ✓ concordant |
| Melanocyte | normal-tracking | falls 7.8→7.3→2.2% | 1×10−14 | ✓ concordant |
| NK | normal-tracking | rises to 7.9% (PP) | 2×10−11 | ✗ discordant |
| Keratinocyte | normal-tracking | rises 79→83→83% | 5×10−3 | ✗ discordant |
| DC | lesional-tracking | flat | 0.20 | trend n.s. |
The three strongest-signal lineages are concordant: the endothelial gradient signal reflects a real compositional increase in bulk tissue, not solely a state change. The discordant cases are informative rather than failures — NK cells infiltrate lesional tissue (bulk proportion rises) yet individual NK cells are depleted from the per-cell gradient-tracking set, and keratinocytes dominate every tier so per-cell direction and bulk fraction are not expected to align. Per-cell tracking and bulk composition are genuinely different measurements; they converge where the biology is a clean compositional shift and diverge otherwise.
Reframing psoriasis as an ordinal gradient and mapping it with a non-circular phenotype yields a coherent, reproducible signal: a gradient-tracking cell population that is monotonic on the NN→PN→PP axis, peaks at the peri-lesional margin, survives two orthogonal significance controls, and is independently corroborated by bulk deconvolution. The dominant biology is vascular — endothelial expansion and an angiogenic gene program — which fits the histology of psoriatic progression and points attention to the dermal vasculature at the advancing edge, not only the epidermal keratinocyte compartment that two-state designs emphasize. STAT3 participates in the program in the expected direction but is not its driver here; its modest single-cell effect suggests any STAT3 role is embedded in a broader vascular/inflammatory circuit rather than acting alone.
Single-cell reference: GSE173706 (Ma et al. 2023). Bulk anchor: SRP165679 (Tsoi et al. 2019, via recount3). All code, figures, result tables, the environment specification, and this white paper are in the repository github.com/soahum-b/psoriasis-1-bulk; the full-census cluster script and reproduce/scale instructions are in HANDOFF.md. Heavy inputs are regenerated by code/00_download_data.R.
Method reference. Scissor: Sun D, Guan X, Moran AE, et al. Identifying phenotype-associated subpopulations by integrating bulk and single-cell sequencing data. Nature Biotechnology 40, 527–538 (2022; published online 11 Nov 2021). doi:10.1038/s41587-021-01091-3.
Working draft — internal white paper for the peri-lesional psoriasis Scissor project. All quantitative values are computed from the saved analysis artifacts. Clinical interpretation is preliminary and not intended to guide patient care.