Bulk RNA-seq datasetsWe collected five publicly available bulk RNA-seq datasets from formalin-fixed-paraffin-embedded (FFPE) and fresh-frozen tumor tissue biopsies, accompanied by clinical information of melanoma patients undergoing immune checkpoint blockade (ICB) therapy. These datasets encompass Auslander et al. (n = 37), Gide et al. (n = 91), Riaz et al. (n = 98), Liu et al. (n = 119), and PUCH (n = 55) datasets9,16,17,18,19. The collection of datasets was frozen by Mar 2023. Each bulk RNA-seq dataset comprises a minimum of 30 samples, ensuring reliable deconvolution for ten distinct cell types. To ensure consistency for deconvolution using CODEFACS we normalized gene expression values to transcripts per million (TPM) for each dataset, following the guidance of Wang et al.15. Notably, all patients across these clinical studies received anti-PD1 monotherapy or a combination of anti-PD1 with anti-CTLA4 except Auslander et al.9, which includes samples (n = 6) treated exclusively with anti-CTLA4 monotherapy. Across all cohorts, we binarized patients’ ICB response to either responder (RECIST criteria: complete/partial response) or non-responder (RECIST criteria: progressive/stable disease). For Auslander et al., Gide et al., and Riaz et al., both pre- and post/on-treatment samples were collected. For Gide et al., Riaz et al., and Liu et al., TPM expression data, along with CODEFACS’ deconvolved expression and cell fraction data for ten cell types were acquired from Wang et al.15. For Auslander et al., expression data was downloaded from GEO (GSE115821), while expression and clinical data for PUCH was downloaded from (https://github.com/xmuyulab/ims_gene_signature). In the PUCH dataset, we converted the overall survival duration into days by multiplying the reported months by 30.In addition, we obtained the deconvolved expression and cell fraction data utilizing CODEFACS for skin cutaneous melanoma patients in the The Cancer Genome Atlas (TCGA-SKCM) from Wang et al.15. Overall survival and progression free interval (synonymous to progression free survival) for TCGA-SKCM patients was downloaded from the UCSC Xena browser53 (https://xenabrowser.net). The pathology classification (brisk or non-brisk) of TCGA-SKCM samples was acquired from Saltz et al.33. Indeterminate (n = 8) and “none” (n = 5) pathology classifications were excluded from the analyses.Deconvolution of bulk ICB RNA-seq datasets using CODEFACSTo deconvolve newly acquired ICB datasets, we utilized cell-type-specific signatures for ten distinct cell types sourced from Wang et al.15. These signatures are tailored to melanoma tumors, encompassing key cell types that best represent the melanoma TME, comprising malignant cells (Mal), skin dendritic cells (skinDC), plasmacytoid dendritic cells (pDC), CD8 + T lymphocytes (TCD8), CD4 + T lymphocytes (TCD4), macrophages, natural killer cells (NK), B lymphocytes (Bcell), endothelial cells (Endo), and cancer associated fibroblasts (CAF). Leveraging the provided signature and bulk datasets, CODEFACS15 estimated the cell-type-specific expression and cell fractions for the identical set of ten cell types across all five deconvolved melanoma ICB cohorts.Inferring “activated” and “inactivated” cell-type-specific ligand-receptor interactions using LIRICSDeconvolved expression data from the five melanoma ICB cohorts and TCGA-SKCM were inputted into LIRICS15, which infers cell-type-specific ligand-receptor interactions, classified as “activated”1 or “inactivated”(0) for each tumor sample within a given cohort. This procedure was conducted for pre-treatment samples and post-treatment samples separately. The scope of these inferred interactions is confined to those cell-type-specific ligand-receptor interactions present within the original curated LIRICS database (n = 3776).Inferring activated cell-type-specific ligand-receptor interaction from single-cell transcriptomics using SOCIALWe developed an R code, SOCIAL (Single-cell transcriptOmics Cell-cell Interaction ALgorithm), to identify significant ligand-receptor interactions between two specific cell types, drawing upon insights from Kumar et al.‘s20, Vento-Tormo et al.‘s21, and our own LIRICS framework15. Our decision to create our own code stemmed from four primary motivations: 1. Leveraging the strengths of previous methods: By combining aspects of the three approaches, we aimed to maximize the accuracy and robustness of our ligand-receptor interaction predictions. 2. Implementing an R-based solution: While the first method lacked publicly accessible code and the second was in Python, we sought to create an R-based solution for accessibility and ease of use. 3. Incorporating our comprehensive database: Our ligand-receptor interaction database (LIRICS) provided rich and informative annotations, enhancing the depth of our analysis. 4. Accommodating variations in ligand-receptor interaction activity observed across patients.SOCIAL (see Supplementary Fig. 11) comprises three main steps: 1. Querying the LIRICS database: Initially, we queried the LIRICS database to identify plausible ligand-receptor interactions; 2. Computing interaction scores: Next, we computed the ligand-receptor interaction score by multiplying the average expression levels of the ligand and receptor complexes for each interaction pair and cell type. 3. Permutation testing: Following that, we performed permutation tests (utilizing 100 iterations in our study) by randomly shuffling cell type labels. This allowed us to derive empirical p-values by calculating the fraction of permutation tests resulting in a higher interaction score than the foreground score determined in step 2. A lower p-value suggests a higher likelihood of the interaction occurring. 4. Optionally, ligand-receptor interactions can be further denoted as significantly activated if the average expression level of both the ligand and receptor genes is greater than the median across all samples.Validate resistant downregulated interactions using single-cell RNA-seq datasetWe obtained single-cell RNA-seq data (scRNA-seq) in TPM encompassing both naïve/untreated (n = 16) and post-ICB resistant (n = 15) resected melanoma tumors from Jerby-Arnon et al.12. Expression data was downloaded from GEO (GSE115978). Unfortunately, the post-ICB responder sample was omitted from our analyzes due to its limited patient sample size (n = 1). To ensure consistency across the bulk datasets and incorporate dendritic cells (DCs), which play a crucial role in immune response through antigen presentation modulation, we reannotated subsets of macrophages as DCs. Utilizing K-means clustering and PCA analysis on twelve classical DC marker genes54,55, this expanded the number of cell types from eight to ten. The reannotation encompassed skin resident DCs (a blend of dermal and epidermal DCs) along with plasmacytoid DCs. Additionally, T cells that did not align with CD4+ or CD8+ subtypes were excluded. Given the variances in available cell types across each tumor sample, we initially consolidated all cells within each treatment group (naïve or post-ICB resistant). Subsequently, we down sampled 40% of cells for each cell type to create “pseudo-samples” corresponding to their respective treatment group. This procedure was repeated through 200 iterations for each treatment group, yielding 200 naïve and 200 post-ICB resistant pseudo tumor samples. To discern the “activated” or “inactivated” status of curated LIRICS interactions using scRNA-seq data, we employed our proprietary R tool, SOCIAL, as discussed earlier, to identify activated ligand-receptor interactions. Interactions meeting the criteria of an empirical p-value < 0.05 and both ligand and receptor mean expression levels exceeding the median across all samples (similar to LIRICS methodology) were classified as “activated” through this procedure.Summary of IRIS IRIS (Immunotherapy Resistance cell-cell Interaction Scanner) operates as a supervised machine learning algorithm, following a two-step process: 1. Differential activation analysis between pre-treatment and post-treatment non-responder group to extract a set of cell-type-specific ligand-receptor interactions, and 2. Employ a hill climbing aggregative feature selection algorithm to further optimize the selected interactions from step 1 by maximizing their classification power in distinguishing responders from non-responders in the pre-treatment phase. Together, these steps aim to infer either a set of “resistance downregulated interactions” (RDI) or “resistance upregulated interactions” (RUI). When multiple training cohorts are available, the training cohorts used in step 1 and step 2 are mutually exclusive and iteratively exchanged resulting in an ensemble model with multiple sets of inferred RDIs or RUIs for a given testing-cohort.Differential activation analysis in step 1The first step of IRIS aims to determine interactions that are either “upregulated” or “downregulated” following the emergence of immune checkpoint blockade (ICB) resistance. This involves using cell-type-specific ligand-receptor interactions as features, and the objective is to identify differentially activated interactions between pre-treatment and non-responder (RECIST criteria: stable/progressive disease) post-treatment samples. This step both establishes the directionality of the interactions following ICB resistance and reduces the search space for the subsequent step 2.Training involves utilizing the cell-type-specific ligand-receptor interaction activity profile, a direct output of LIRICS. When multiple training cohorts are available, interaction activity profiles are merged, including the union of inferred interactions within each cohort. The initial feature-space can consist of up to 3776 possible ligand-receptor interactions between ten cell types within the melanoma TME, sourced from the LIRICS-curated database. A fisher’s exact test is employed to identify interactions that are significantly activated (FDR < 0.2 per cell type pair) in either the post-treatment non-responder or the pre-treatment groups, categorizing the selected interactions as either resistance “upregulated” or resistance “downregulated” respectively. The output of step 1 is a ranked list (by FDR values) of interactions that are categorized as either RUIs or RDIs. These interactions compose the search space for step 2.Hill climbing aggregative feature selection for interactions that maximize classification of patient response to ICB in step 2The second step of IRIS aims to identify an optimal set of interactions that are functionally relevant in dictating ICB response. This involves starting from the list of interactions (RDI or RUI) identified in the previous step as features, and the objective is to select for an optimal set of interactions that maximizes the classification power in distinguishing responders (RECIST criteria: partial/complete response) and non-responders (RECIST criteria: stable/progressive disease) in the pre-treatment phase. The method of step 2 was developed by amalgamating components of the feature selection methodologies outlined in our previous work from Wang et al15. and Auslander et al.9.The algorithm consists of an iterative procedure comprising three sequential stages. First, the pre-treatment samples in the training cohort are randomly split into three folds: two for training and one for testing. Second, starting from an empty set, the procedure greedily adds the most effective discriminating interactions (ranked by FDR by cell-pair from step 1) in a step-by-step manner until no further improvement in the classification accuracy (in terms of AUC) of ICB response is achieved within the training-folds. Considering the directional nature of individual features established in stage 1, we anticipate that for RUIs from step 1, less activations correspond to a response. Conversely, for RDIs, higher activation levels correlate with a response. Third, the selected set of interactions from the training folds together is evaluated in terms of classification accuracy and scored on the testing-fold. A test AUC > = 0.6 results in a reward score of 1 for all selected set of interactions, while a test AUC < 0.4 incurs a penalty score of -1. AUCs in-between and all unselected features received a score of 0. This process is repeated for 500 iterations.Following 500 iterations, we sum the scores for each individual interaction across all iterations, termed feature score. We identify the features with the greatest feature score from the entire 500 solutions, which approximates a solution that minimizes the risk of overfitting. During the feature selection, it is possible to estimate the empirical p-value linked to each feature having a higher feature score than the feature score obtained by random chance. This estimation is accomplished by shuffling the values of each solution derived from the hill climbing feature selection, repeating this procedure 1000 times for each iteration. By performing 1000 random permutations of 500 solutions, a null distribution representing the feature score via random chance is constructed. We selected features with an empirical p-value < 0.05 as either our RDIs or RUIs depending on the input features directionality. The empirical p-value for each feature j is estimated as follow:$${p}_{j}=\frac{1}{1000}{\sum }_{i=1}^{1000}{1}_{{{{\rm{null}}}}\, {{{\rm{feature}}}}\, {{{\rm{score}}}}_{{ij}}\, > \, {{{\rm{observed}}}}\, {{{\rm{feature}}}}\, {{{\rm{score}}}}_{{ij}}}$$
(1)
In this study, the pre-treatment cohort from Auslander et al. was excluded during step 2 feature selection, owing to the limited number of responders, which hindered the reliable determination of prediction accuracy.Calculating RDS and RUSFollowing the identification of RDI or RUI in melanoma, the calculation of each tumor samples’ resistance-downregulated score (RDS) and resistance-upregulated score (RUS) for a given tumor sample is computed as followed:$${\mbox{RDS or RUS}}=\frac{{\sum }_{i=1}^{n}\frac{{{{\rm{number}}}} \, {{{\rm{of}}}}\, {{{\rm{active}}}}\, {{{\rm{inferred}}}}\, {{{\rm{RDI}}}} \, {{{\rm{or}}}}\, {{{\rm{RAI}}}}_{i}}{{{{\rm{total}}}}\, {{{\rm{number}}}}\, {{{\rm{of}}}}\, {{{\rm{inferred}}}}\, {{{\rm{RDI}}}}\, {{{\rm{or}}}}\, {{RAI}}_{{{\rm{i}}}}}}{n}$$
(2)
First calculate the fraction of activated inferred interactions derived from each set of RDIs or RUIs inferred (i) within the ensemble model. Then average the fraction of activated inferred interactions across all (n) sets of RDIs or RUIs within the ensemble model. The ultimate RDS and RUS scores are scaled for cross-cohort comparisons.A merged ensemble model is also derived by combining all individual RDI ensemble models originally inferred for each cohort: Gide et al., Liu et al., Auslander et al., Riaz et al., and PUCH. This unified merged model is applied to the two independent patient cohorts: TCGA-SKCM and Jerby-Arnon et al.Calculating odds ratioTo establish the proper classification threshold for distinguishing true responder and true non-responder from false responder and non-responder for each testing cohort, we considered both the set of RDIs inferred, and tumor samples extracted from each pre-treatment training cohort used in step 2 are considered. For each training cohort, the proportion of activated inferred interactions in each sample are calculated and scaled across all tumor samples. Using the scaled score in conjunction with response information for each tumor sample in the training cohort, we compute the optimal maximizing cut point using the cutpointr() function from the cutpointr package (v1.1.2) in R56. The final threshold is obtained by taking the mean of the individual thresholds computed using each training cohort.Benchmarking with other state-of-the-art ICB predictorsTo benchmark IRIS’ performance, we computed the accompanying score for each state-of-the-art ICB predictor in each tumor sample following the original method. The T cell exhaustion (Tex) signature was computed by averaging the expression of available genes within the signature, including LAG3, TIGIT, PDCD1, PD1, HAVCR2, TIM3, and CTLA4. The TIDE score was computed using the online tool (http://tide.dfci.harvard.edu). The resF signature was computed using the author’s original published code available from (https://github.com/livnatje/ImmuneResistance). For Melanocytic plasticity score (MPS)10, IFNG signature57, Cytotoxic signature11, T cell inflamed GEP (Tin)58, and IMPRES9, we computed the signature scores following the exact methods detailed in the original publication.Spatial RNA-seq datasetWe collected two publicly available metastatic melanoma spatial transcriptomics datasets from Thrane et al37. and Biermann et al.36, acquired from the legacy and SlideSeqV2 platforms respectively. Thrane et al. includes slides of treatment naïve melanoma lymph node metastasis (n = 8). Biermann et al. includes slides of treatment naïve melanoma brain metastasis (n = 9), treatment naïve extracranial melanoma metastasis (n = 2), and anti-CTLA4 post-treatment extracranial melanoma metastasis (n = 2). The two post-treatment samples were derived from non-responder (RECIST criteria: progressive disease) patients to the ICB regimen. Slides MBM07.1, MBM13.1, and ECM08.1 from Biermann et al. were excluded from our analyzes due to insufficient spatial coordinate coverage or prior adoptive T cell therapy. Additionally, we collected all single-nuclei RNA-seq (snRNA-seq) data from Biermann et al. of profiled metastatic melanoma samples (n = 28). Legacy data from Thrane et al. was downloaded from (https://www.spatialresearch.org/resources-published-datasets/), while SlideSeqV2 and snRNA-seq data for Biermann et al. was downloaded from GEO (GSE185386).Aligning single-cell to spatial transcriptomics using CytoSPACETo infer the activity of ligand-receptor interactions in spatial data using SPECIAL (see below), we first deconvolved spatial transcriptomics (ST) data to single-cell resolution utilizing CytoSPACE (v1.0.6)38 (https://github.com/digitalcytometry/cytospace). CytoSPACE requires two inputs: 1. count matrices for a spatial transcriptomic slide, and 2. count matrices for a reference single-cell transcriptomics dataset with cell type annotations. It outputs deconvolved spatial transcriptomics data, assigning spatial coordinates for selected individual single-cells from the references scRNA-seq dataset. The cell abundance within each coordinate of the spatial transcriptomic slide for each cell type is computed as the fraction of cells within that specific region.For bulk ST data from Thrane et al. (using the legacy platform), we assigned spatial spots with scRNA-seq data (in counts) of treatment naïve samples (n = 16) from Jerby-Arnon et al.12 as the reference. The cell types were limited to the ten specified cell types mentioned earlier in our single-cell analysis. To optimize CytoSPACE performance for smart-seq counts data from Jerby-Arnon et al., we estimated cell fraction composition for each slide using Spatial Seuart59 as user-provided cell fraction input for CytoSPACE. Additionally, we set the geometry as square, down-sampling off, and mean cell numbers to 20 following the guidance of Vahid et al. when applying CytoSPACE to legacy spatial platforms and smart-seq reference data38.For SlideSeqV2 platform (i.e., single-cell resolution ST) from Biermann et al., we first assigned each spatial spots with individual cell types using the RCTD60 pipeline in doublet mode from the spacexr package in R. Each SlideSeqV2 slide was annotated using all metastatic melanoma snRNA-seq data (n = 28) with cell type annotations. The cell types from the snRNA-seq data were manually reannotated as shown in Supplementary Table 1. Doublets, “low-quality” cells, and cycling cells were excluded. The output of RCTD included the first cell type annotation for all spatial spots and the second cell type annotation only if RCTD inferred the spatial spot as a “doublet certain”. The RCTD framework was adopted from Biermann et al.’s methods36. Next, we executed CytoSPACE (single-cell mode) twice for each SlideSeqV2 slide from Biermann et al. using only the matched patient snRNA-seq data as a reference. The first run incorporated all spatial spots with their respective first cell type annotation from RCTD as user-provided input, while the second run was limited to “doublet certain” spatial spots with their respective second cell type annotation. Cell types absent in either the matched patient spatial transcriptomics slide or snRNA-seq were excluded before running CytoSPACE.Inferring activated cell-type-specific ligand-receptor interaction from single-cell transcriptomics aligned to spatial transcriptomics using SPECIALTo quantify the activity of cell-type-specific ligand-receptor interactions within each spatial transcriptomics slide, we further developed our in-house single-cell ligand-receptor inference tool called SOCIAL, into SPECIAL (SPatial transcriptomics cEll-Cell Interaction ALgorithm; see Supplementary Fig. 12). This novel iteration is customized specifically for spatial transcriptomics with aligned single-cell transcriptomes, the direct output of CytoSPACE (see above). The three major steps of SPECIAL are outlined below (see Supplementary Fig. 12).In the first step, spatial coordinates for each slide are divided into spatial “regions” with a diameter set to approximately 250 μm, a distance conducive to paracrine ligand-receptor interactions based on literature39, via two possible methods. For bulk ST platforms (i.e., Legacy platform or Visium 10X), we developed a sliding window approach. This approach consolidates cells within a specified radius at each labeled spatial coordinate. The number of regions is contingent upon the number of spatial coordinates sequenced within a given slide, and each region’s spatial coverage can overlap with multiple other regions. For SlideSeqV2 platform (or similar single-cell resolution ST platforms with a circular sequenced area or “puck”) we developed a K-means clustering-based approach. This method clusters individual cells based on their x and y coordinates into non-overlapping circular regions based on a specified diameter. Here, the number of regions is determined by the number of circle regions with diameter i that can fit within a circular sequenced area of diameter j. For both approaches, regions containing a single cell type are omitted from downstream analysis.In this study, we employed the sliding window approach to Thrane et al.’s cohort (Legacy platform), consolidating cells within a 1-unit radius (maximum 300 μm diameter region). While for Biermann et al.’s cohort (SlideSeqV2 platform), we utilized the K-means clustering approach, clustering each slide into 144 regions. We chose 144 regions based on the number of 250 μm diameter circle regions that can fit within the 3000 μm diameter circle “puck” sequenced by the SlideSeqV2 platform.In the second step, to discern the “activated” or “inactivated” status of curated LIRICS interactions within each spatial region, we employed the same proprietary R tool, SOCIAL, described earlier. Specifically, SOCIAL is independently applied to each region, encompassing the respective single-cell transcriptomes.In the third step, interactions with an empirical p-value < 0.05 and both ligand and receptor mean expression greater than the median across all regions (analogous to LIRICS) are categorized as “activated” through this process.Reporting summaryFurther information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
A machine learning model reveals expansive downregulation of ligand-receptor interactions that enhance lymphocyte infiltration in melanoma with developed resistance to immune checkpoint blockade
