DOI: 10.5061/dryad.s4mw6m9jh
title: Single-cell spatial transcriptomics uncovers niches that govern response to PD-1/PD-L1 blockade in cutaneous squamous cell carcinoma

description: Neoadjuvant PD-1/PD-L1 blockade demonstrates robust efficacy in cutaneous
 squamous cell carcinoma, yet a significant proportion of patients do not
 achieve a major pathological response for unclear reasons. Here, we
 characterize spatially resolved gene expression profiles before, during,
 and after treatment in cutaneous squamous cell carcinoma (cSCC) patients
 enrolled on two pivotal phase II clinical trials of neoadjuvant PD-1/PD-L1
 blockade. Using single-cell RNA-based spatial molecular imaging, we
 profiled more than 56 million transcripts from over 250,000 cells across
 27 patients. Traditional cell type clustering uncovered five broad cell
 types (keratinocytes, fibroblasts, B-lineage cells, myeloid cells, and a
 mixed dendritic/T/NK compartment). Genes associated with interferon
 signaling and metabolic pathways were significantly upregulated in
 responders, but bulk gene expression alone did not fully distinguish
 responders from non-responders. We then applied a spatial clustering
 approach, which revealed six distinct spatial niches defined by their
 differential gene expression signatures. Notably, these spatial niches
 were differentially enriched in responders versus non-responders,
 suggesting niche-specific mechanisms of immunotherapy sensitivity and
 resistance. In-depth analyses of individual patients demonstrated that
 these spatial niches correspond to diverse mechanisms of immune evasion,
 including impaired interferon-stimulated antigen presentation,
 immunosuppressive myeloid infiltration, or epithelial–mesenchymal
 transition. Taken together, our study provides a comprehensive single-cell
 spatial transcriptomic atlas of cSCC evolution during PD-1/PD-L1 blockade,
 highlighting the heterogeneity of resistance mechanisms. These findings
 illuminate potential targets for rational combination therapies and
 underscore the importance of spatial context in understanding and
 overcoming immunotherapy resistance in cSCC.
description: # Single-cell spatial transcriptomics uncovers niches that govern response
 to PD-1/PD-L1 blockade in cutaneous squamous cell carcinoma Dataset DOI:
 [10.5061/dryad.s4mw6m9jh](10.5061/dryad.s4mw6m9jh) ## Description of the
 data and file structure This repository contains a Seurat object (named
 “seurat_object”) generated from CosMx SMI using the manufacturer's
 1000 plex panel. The object enables researchers to reproduce or extend the
 analyses performed in our study. All cell metadata has been de-identified
 according to our IRB-approved protocol, and no personally identifiable
 information (PII) is included. The data comes from 2 tissue microarrays of
 formalin fixed paraffin embedded tissue of human cutaneous SCC specimens
 taken before, during, or after immunotherapy.  ### Files and variables
 #### Seurat_Object_cSCC.RDS ``` – A serialized R object in Seurat format,
 which includes: – Filtered and normalized gene expression data. –
 Cell-level metadata (excluding direct identifiers). – Cluster assignments
 and other annotations used in the study. ``` **Data De-identification** •
 In compliance with applicable regulations and IRB guidelines, all
 patient‐specific identifiers have been removed or replaced with hashed or
 random codes. • Only variables relevant to the analyses are retained
 within the Seurat object. **Using the data** To load the Seurat object in
 R: ``` library(Seurat) # Replace "path/to/" with the local path
 where you placed the .rds file: seurat_object &lt;-
 readRDS("path/to/seurat_object.rds")` ``` After loading, you can
 explore gene expression levels, run Seurat’s built-in visualization tools,
 or incorporate the object into your own pipelines. **Requirements** • R
 (version 4.0 or higher recommended) • Seurat (version 4 or higher
 recommended) ## Code/software **Creating_object_and_clustering_cells.R**
 This code demonstrates how to preprocess, transform, and cluster
 single‐cell (or spatial) transcriptomic data using Seurat. It includes
 filtering low‐count cells, performing SCTransform, detecting/removing
 doublets, and using PCA/UMAP for dimension reduction. The pipeline also
 applies InSituType for immunofluorescence‐based classification and
 identifies top differentially expressed markers in each cluster for
 downstream biological insights. **GSEA.R** This code demonstrates a
 pipeline for performing spatially resolved GSEA, starting from normalized
 gene expression data and metadata about cell locations. It uses GSVA to
 compute per‐cell and summarized enrichment scores for hallmark gene sets
 and then compares them across samples, clusters, and time points. Finally,
 the code generates heatmaps, plots, and statistical outputs to help
 interpret how these gene sets vary within and across patients and biopsy
 timepoints. **InSituCor_spatial_gene_modules.R** This code uses InSituCor
 to identify gene modules that are spatially co‐expressed across the
 tissue. After normalizing and preparing the data (including single‐cell
 metadata and neighboring information), it fits correlation modules and
 calculates both per‐cell and environment scores. Finally, it visualizes
 the modules via heatmaps and correlation networks, providing a functional
 look at how gene expression groups are spatially organized.
 **Moran_s_I.R** This code calculates bivariate Moran’s I statistics to
 measure spatial autocorrelation between ligand–receptor pairs (e.g., IFN-γ
 pathways) across multiple fields of view (FOVs). It subsets metadata and
 gene expression, identifies nearest neighbors for each cell, and then
 computes Moran’s I along with bootstrap confidence intervals. Finally, the
 script merges these results with patient metadata for comparison and
 visualization of spatial ligand–receptor interaction patterns by sample or
 clinical grouping. **ROC_analysis.R** This code performs ROC (receiver
 operating characteristic) analysis to assess the ability of various
 features—namely the proportion of “favorable niches” (spatial cluster
 composition) and PD‐L1 expression levels—to distinguish between responders
 and nonresponders. It merges and filters relevant metadata, calculates AUC
 values via the pROC package, identifies optimal cutoffs
 (sensitivity/specificity), and uses statistical tests (e.g., Wilcoxon) to
 compare distributions of those features by response category, ultimately
 creating ROC plots to visualize classifier performance.
 **Spatial_Clustering.R** This code demonstrates an approach for clustering
 cells based on their local neighborhood gene‐expression profiles. It first
 computes a nearest‐neighbor graph, extracts neighborhood‐level expression
 values (e.g., average per gene in each cell’s local neighborhood), and
 then performs PCA on these neighborhood features. Next, it uses Mclust to
 define six spatial clusters, which are integrated into the Seurat
 metadata. Finally, the script generates bar plots and summary tables
 examining how these clusters vary across biopsy timepoints, patients, and
 clinical response groups, and identifies spatial‐cluster‐specific marker
 genes. **Transcriptome_Correlation.R** This code computes an inter-patient
 Pearson correlation matrix of average transcriptome profiles.
 Specifically, it calculates gene‐expression averages for each patient and
 then generates a correlation heatmap, splitting rows and columns by
 response type. This helps visualize which patients have similar global
 expression patterns and whether these clusters align with clinical
 outcomes. ## Human subjects data All patient data used in this study were
 de-identified in accordance with institutional review board (IRB)
 requirements and relevant data privacy regulations (e.g., HIPAA).
 Specifically, the research team removed or replaced all direct and
 indirect identifiers, including names, dates of birth, addresses, and
 medical record numbers, with coded pseudonyms or aggregated values. No
 link to identifiable personal information was retained, ensuring that the
 data cannot be traced back to individual participants. All de-identified
 data were stored in secure, access-controlled repositories with strict
 oversight, and were only used for the purposes of this study.