The tumor microenvironment (TME) is a complex mixture of cell types whose interactions affect tumor growth and clinical outcome. To discover such interactions, we developed CODEFACS (COnfident DEconvolution For All Cell Subsets), a tool deconvolving cell type–specific gene expression in each sample from bulk expression, and LIRICS (Ligand–Receptor Interactions between Cell Subsets), a statistical framework prioritizing clinically relevant ligand–receptor interactions between cell types from the deconvolved data. We first demonstrate the superiority of CODEFACS versus the state-of-the-art deconvolution method CIBERSORTx. Second, analyzing The Cancer Genome Atlas, we uncover cell type–specific ligand–receptor interactions uniquely associated with mismatch-repair deficiency across different cancer types, providing additional insights into their enhanced sensitivity to anti–programmed cell death protein 1 (PD-1) therapy compared with other tumors with high neoantigen burden. Finally, we identify a subset of cell type–specific ligand–receptor interactions in the melanoma TME that stratify survival of patients receiving anti–PD-1 therapy better than some recently published bulk transcriptomics-based methods.

Significance:

This work presents two new computational methods that can deconvolve a large collection of bulk tumor gene expression profiles into their respective cell type–specific gene expression profiles and identify cell type–specific ligand–receptor interactions predictive of response to immune-checkpoint blockade therapy.

This article is highlighted in the In This Issue feature, p. 873

The importance of the tumor microenvironment (TME) in cancer has been recognized since the late 1800s (1). The recent success of immune-checkpoint blockade (ICB) has further sparked interest in studying the various cellular states of cell types and their interactions in the TME, aiming to find biomarkers of treatment response and new treatment opportunities (2). One key step in elucidating the different cellular states of cell types de novo is the characterization of the molecular profiles of each cell type in each patient's tumor sample. FACS and single-cell RNA sequencing (scRNA-seq) have emerged as effective tools to address this challenge (3). The use of FACS is limited because it can measure expression of only a few protein markers simultaneously. The use of scRNA-seq has been limited due to its cost and the scarcity of fresh tumor biopsies. As bulk tumor gene expression from preserved biopsies accompanied by clinical outcome metadata is abundant, computational methods that can effectively extract cell type–specific expression from such data in each individual sample could be very helpful. If successful, such deconvolution methods could be used to identify clinically relevant cellular states of cell types and prioritize cell–cell interactions associated with patient survival and response to treatment. Furthermore, they may be readily applied to interrogate the troves of large bulk expression data sets that are in the public domain, to further advance what we can learn from them.

Recovery of cell type–specific gene expression profiles from bulk gene expression is a computationally challenging problem because the standard formulation requires one to solve a large system of underdetermined linear equations. This problem is closely related to the problem of compressed sensing in signal processing (4). Several recent studies have developed a variety of algorithms to address this challenge. DeMix (5) was designed to estimate individual-specific expression for three cell components when provided prior reference samples of two of these cell components. ISOpure (6) has aimed to derive sample-specific tumor cell expression only, with the assumption that the observed bulk gene expression profile is a mixture of predefined stromal and immune cell expression profiles that are shared across all the samples. Building on this work, Fox and colleagues extended ISOpure to predict individual-specific nontumor cell expression by subtracting cancer cell expression profiles from the bulk mixtures in a two-cell-type model (7). More recently, Newman and colleagues (8) developed CIBERSORTx, the first approach that aims to predict the sample-specific gene expression of all cell types composing it by using a set of novel deconvolution heuristics. As a proof of concept, Newman and colleagues showed that CIBERSORTx can accurately reconstruct the cell type–specific expression of genes in each input sample under certain modeling assumptions. This groundbreaking work has, however, some notable limitations: (i) The number of genes whose cell type–specific expression can be reconstructed in each sample is quite small, especially for low-abundance cell types, and (ii) CIBERSORTx does not provide confidence estimates of the predictions made. Confidence estimates on the expression predictions could be useful in most deconvolution applications because of the absence of ground-truth data.

Starting from the work of Newman and colleagues, we introduce a new deconvolution algorithm and software, CODEFACS (Confident DEconvolution For All Cell Subsets), that markedly advances our ability to successfully reconstruct cell type–specific gene expression of each bulk sample. CODEFACS receives as input bulk gene expression profiles of tumor samples and either precomputed estimates of abundance of expected tumor, immunologic and stromal cell types in each sample, or their prototypical molecular signatures, which serve as seeds for estimating the abundance of each cell type in each sample. CODEFACS then predicts the cell type–specific gene expression profiles in each sample. It is a heuristic approach aimed at maximizing the number of genes in each cell type whose expression across the samples can be estimated as confidently predicted. Using 15 benchmark data sets where the ground truth is known, we show that CODEFACS robustly improves over CIBERSORTx, in terms of both gene coverage and the individual gene expression estimation accuracy. We focus our comparisons on CIBERSORTx because it is widely used (cited more than 350 times as measured by Google Scholar since its publication in 2019), and, to our knowledge, it is the only widely used package that makes sample-specific estimates of gene expression for each cell type. The essential methodologic contributions of CODEFACS that distinguish CODEFACS from CIBERSORTx are: (i) a principled, quantitative definition of the confidence with which the expression of a gene or cell type has been predicted (which is highly important, if not essential, in any application) and (ii) two new methods of iterative refinement to leverage the high-confidence predictions to improve predictions for genes predicted with low confidence at an earlier stage, implemented with a modular software engineering design.

We additionally introduce LIRICS (LIgand–Receptor Interactions between Cell Subsets), which integrates the output of CODEFACS with a database of prior immunologic knowledge that we curated to infer the cell type–specific ligand–receptor pairs that are likely to interact in each sample. These data can then be analyzed in conjunction with any sample-associated clinical annotations (e.g., response to treatment) to systematically prioritize clinically relevant immune-related interactions between any pair of cell types in each patient's cancer cohort. This opens the possibility of identifying phenotypic/clinical associations of individual interactions in large clinical tumor expression cohorts.

Building on the enhanced coverage and accuracy of CODEFACS, we applied it to reconstruct the cell type–specific transcriptomes of approximately 8,000 tumor samples from 21 cancer types in The Cancer Genome Atlas (TCGA). Analyzing these fully deconvolved TCGA expression data sets using LIRICS, we find a shared repertoire of cell type–specific ligand–receptor interactions unique to the TME of mismatch repair–deficient (MMRD) tumors from different tissues of origin. These interactions provide additional insights into the very high response rates of MMRD tumors to anti–programmed cell death protein 1 (PD-1) treatment compared with other highly mutated tumors. Finally, using machine learning techniques, we identify a subset of intercellular TME interactions that stratify survival outcomes of patients with melanoma receiving anti–PD-1 therapy better than some recently published transcriptomics-based methods.

In summary, CODEFACS and LIRICS effectively build upon the statistical power of large bulk RNA-seq data sets and prior immunologic knowledge to characterize sample-specific cell type–specific gene expression and learn more about the association of different tumor–immune interactions with different clinical measures. Notably, the potential scope of applications of both CODEFACS and LIRICS goes beyond studying the TME, as these tools can be applied to study any disease of interest given bulk gene expression data and relevant reference signatures of cell types involved.

Overview of CODEFACS and LIRICS

CODEFACS is designed to characterize the TME by reconstructing the cell type–specific transcriptomes of each sample from bulk expression. It takes as input the bulk RNA-seq expression values of a cohort of tumor samples and either the estimations of the cell fractions of a predefined set of cell types in each sample or their cell type–specific molecular signature profiles, derived based on reference data sets or from the literature.

CODEFACS then employs a heuristic approach that sequentially executes three modules: (module 1) high-resolution deconvolution, (module 2) hierarchical deconvolution, and (module 3) imputation (see Fig. 1A and Supplementary Note 1). In module 1, we perform a high-resolution deconvolution, which extends the CIBERSORTx algorithm. In module 2 (hierarchical deconvolution), bulk expression is modeled as a mixture of two components: a specific cell type of interest and all the remaining cell types. The expression for the cell type of interest is predicted by removing the estimated expression in the second component (using high-resolution deconvolution from module 1) from the bulk mixture. In module 3—imputation-based deconvolution—we impute the cell type–specific expression of a specific gene based on the predicted cell type–specific expression of other high-confidence genes that are coexpressed with that gene in the bulk mixture. Each module is designed to overcome the shortcomings of its predecessor based on their respective modeling assumptions.

Figure 1.

Overview of CODEFACS and LIRICS. A, CODEFACS takes bulk gene expression profiles and prior knowledge of the cellular composition of each sample and executes a heuristic three-step procedure to infer the deconvolved gene expression in each sample, as described earlier and in more detail in Supplementary Note 1. B, LIRICS takes the output of CODEFACS and processes it in three steps, as described in the main text and in detail in Supplementary Note 2.

Figure 1.

Overview of CODEFACS and LIRICS. A, CODEFACS takes bulk gene expression profiles and prior knowledge of the cellular composition of each sample and executes a heuristic three-step procedure to infer the deconvolved gene expression in each sample, as described earlier and in more detail in Supplementary Note 1. B, LIRICS takes the output of CODEFACS and processes it in three steps, as described in the main text and in detail in Supplementary Note 2.

Close modal

A key component of CODEFACS is its confidence ranking system, which receives cell type–specific expression predictions from the different modules and labels them as high- or low-confidence estimations. (For an overview of the rationale behind the confidence ranking system, see Methods.) Genes whose expression is determined with high confidence in each module are added to the output set, whereas low-confidence predictions are continued to be processed in subsequent modules. The final output of CODEFACS consists of two items: (i) a three-dimensional gene expression matrix, where each entry represents the predicted gene expression of a gene in a given cell type in a specific sample, and (ii) a two-dimensional matrix of confidence scores ranging from 0–1 representing which gene–cell type pairs have the most confident predictions (Fig. 1A, Output).

Given fully deconvolved gene expression data from CODEFACS, one can use LIRICS (Fig. 1B) to infer which cell type–specific ligand–receptor pairs are likely to interact in each sample. Specifically, LIRICS takes the output of CODEFACS and processes it in three steps. Step 1: The first step queries a database of all plausible ligand–receptor interactions between any two cell types A and B, which we have systematically assembled and curated from the literature. This database is publicly available as part of LIRICS (see Supplementary Note S2; Supplementary Tables S1–S3) and serves as a repository of all prior immunologic knowledge. Step 2: In the second step, given the deconvolved expression profiles of cell types A and B in each bulk tumor sample, LIRICS denotes as “active” or “likely to occur” (“1”) the interactions where both the ligand and receptor are overexpressed in the relevant cell types in that sample, or otherwise “inactive” (“0”). A ligand or receptor is overexpressed in each cell type if its expression exceeds the median expression in that cell type (Supplementary Note 2). Step 3: In the third step, a Fisher enrichment analysis is performed to test the association of the activity of specific ligand–receptor interactions with any relevant phenotypes of interest (e.g., treatment response and mutational subtype; Fig. 1B).

Evaluating the Performance of CODEFACS versus Numerous Different Benchmarks

To assess the accuracy of CODEFACS, we generated 15 benchmark data sets (see Methods) by merging publicly available scRNA-seq (9, 10) and FACS-purified RNA-seq (11). Thereafter, we applied CODEFACS to deconvolve these generated bulk data sets and define the accuracy of its predictions by computing the Kendall correlation between the predicted and ground-truth expression in each cell type across individual samples (the Kendall correlation provides a less inflated measure of accuracy by accounting for ties in the data). In the main text, we show the results obtained on three of these benchmark bulk data sets: one derived from FACS lung cancer data (11), one from a single-cell melanoma RNA-seq data set (9), and one from a single-cell glioblastoma (GBM) RNA-seq data set (10), respectively. Each bulk sample from these data sets represents a biopsy from a patient. We show that CODEFACS predicts the cell type–specific expression of many more genes than CIBERSORTx (with Kendall correlation ≥0.3; Fig. 2AC), and its predictions are overall more accurate (Fig. 2DF). Furthermore, the estimated cell type–specific expression profiles from CODEFACS recapitulate known tumor cell subtypes derived from single-cell data (Supplementary Fig. S1 demonstrates with the GBM data set as an example). The results for all the remaining 12 benchmark data sets, created via artificial mixing of single-cell profiles and simulation of batch effects between bulk and single-cell expression data (see Methods and Supplementary Table S4 for more details), further show the superiority of CODEFACS (Supplementary Figs. S2A–S2L and S3A–S3L). The comparisons based on different cutoffs (the baseline of 0.3, Kendall correlation ≥0.4, and Kendall correlation ≥0.5) also support our conclusion (Supplementary Figs. S2A–S2L, S3A–S3L, S4A–S4O, and S5A–S5O). Overall, we observe that the more abundant the cell type is, the better CODEFACS can predict its cell type–specific gene expression (Supplementary Fig. S6).

Figure 2.

Evaluating the performance of CODEFACS. A–C, Bar plots depicting the number of genes with a prediction accuracy (Kendall correlation) ≥0.3 with the ground truth for each cell type, estimated from bulk-generated samples of lung cancer [non–small cell lung cancer (NSCLC) data set; sample size = 26], melanoma (SKCM data set 1; sample size = 28) and glioblastoma (GBM data set 1; sample size = 24) benchmark data sets, as estimated by CODEFACS (yellow bars) and CIBERSORTx (blue bars). D–F, Box plots depicting prediction accuracy distributions of all genes across different cell types in the lung cancer (NSCLC with sample size = 26), melanoma (SKCM with sample size = 28) and glioblastoma (GBM with sample size = 24) benchmark data sets, using CODEFACS (yellow) and CIBERSORTx (blue). A two-sided Wilcoxon signed rank test was performed to compare the prediction accuracies of CODEFACS and that of CIBERSORTx for each cell type in each data set. ***, P < 0.001. G, Spearman correlations between prediction accuracies and confidence scores among cell types in the lung cancer benchmark data set (NSCLC data set; sample size = 26). The y-axis indicates the Spearman correlation coefficient value and the x-axis indicates the cell type. H, AUCs obtained in classifying informative and uninformative predictions among cell types in the lung cancer benchmark data set (NSCLC data set; sample size = 26). I, Bar plots depicting the Spearman correlations between mean deconvolved cell type–specific expression in TCGA LUAD, and mean cell type–specific expression derived from publicly available single-cell data sets of LUAD. The overlapping gene set between high-confidence deconvolved genes in TCGA LUAD and genes in the single-cell data for each cell type are used for computing the Spearman correlation. The sizes of the gene sets for the nine cell types are listed as follows according to the cell-type order in the figure: 11,765, 8,342, 7,881, 8,790, 9,132, 6,743, 7,448, 10,379, and 7,267. The y-axis indicates the Spearman correlation coefficient value and the x-axis indicates the cell type. Details of all scRNA-seq–based pseudobulk data sets can be found in Supplementary Table S4. In addition, Supplementary Fig. S25 shows a scatter plot of measured mean gene expression in malignant LUAD single cells versus CODEFACS-predicted gene expression for malignant cells in bulk TCGA LUAD samples.

Figure 2.

Evaluating the performance of CODEFACS. A–C, Bar plots depicting the number of genes with a prediction accuracy (Kendall correlation) ≥0.3 with the ground truth for each cell type, estimated from bulk-generated samples of lung cancer [non–small cell lung cancer (NSCLC) data set; sample size = 26], melanoma (SKCM data set 1; sample size = 28) and glioblastoma (GBM data set 1; sample size = 24) benchmark data sets, as estimated by CODEFACS (yellow bars) and CIBERSORTx (blue bars). D–F, Box plots depicting prediction accuracy distributions of all genes across different cell types in the lung cancer (NSCLC with sample size = 26), melanoma (SKCM with sample size = 28) and glioblastoma (GBM with sample size = 24) benchmark data sets, using CODEFACS (yellow) and CIBERSORTx (blue). A two-sided Wilcoxon signed rank test was performed to compare the prediction accuracies of CODEFACS and that of CIBERSORTx for each cell type in each data set. ***, P < 0.001. G, Spearman correlations between prediction accuracies and confidence scores among cell types in the lung cancer benchmark data set (NSCLC data set; sample size = 26). The y-axis indicates the Spearman correlation coefficient value and the x-axis indicates the cell type. H, AUCs obtained in classifying informative and uninformative predictions among cell types in the lung cancer benchmark data set (NSCLC data set; sample size = 26). I, Bar plots depicting the Spearman correlations between mean deconvolved cell type–specific expression in TCGA LUAD, and mean cell type–specific expression derived from publicly available single-cell data sets of LUAD. The overlapping gene set between high-confidence deconvolved genes in TCGA LUAD and genes in the single-cell data for each cell type are used for computing the Spearman correlation. The sizes of the gene sets for the nine cell types are listed as follows according to the cell-type order in the figure: 11,765, 8,342, 7,881, 8,790, 9,132, 6,743, 7,448, 10,379, and 7,267. The y-axis indicates the Spearman correlation coefficient value and the x-axis indicates the cell type. Details of all scRNA-seq–based pseudobulk data sets can be found in Supplementary Table S4. In addition, Supplementary Fig. S25 shows a scatter plot of measured mean gene expression in malignant LUAD single cells versus CODEFACS-predicted gene expression for malignant cells in bulk TCGA LUAD samples.

Close modal

To test whether CODEFACS confidence scores can filter out potentially noisy predictions, we quantified how well these scores align with the Kendall scores that measure the prediction accuracy with the ground truth for each (gene, cell type) pair, as described above. We quantified this using two metrics: Spearman correlation and a classification area under the ROC curve (AUC). As evident, this correlation is strong and positive (Fig. 2G depicts the results for the FACS lung cancer benchmark data set, and Supplementary Fig. S7A–S7N depicts the results for the remaining benchmark data sets). To perform a classification-based quantification, we grouped the genes in each cell type into two classes based on the correlation between their predicted and actual expression, as informative (Kendall score ≥0.1 and P value ≤ 0.05) and uninformative (Kendall score <0.1 or P > 0.05). We then tested whether the confidence scores could be used to classify genes into these two classes for each cell type. We found that the confidence score could effectively filter out uninformative predictions (Fig. 2H depicts the results for the FACS-sorted lung cancer benchmark data set, and Supplementary Fig. S8A–S8N depicts the results for the remaining benchmark data sets). Supplementary Fig. S9A–S9O visualizes the correlation between predicted and measured gene expression. This figure reveals that the predictions are generally less accurate for low-expressed genes compared with genes with high expression. Supplementary Fig. S10A–S10O shows scatter plots of the confidence score (x-axis) and the measured expression (y-axis). Confidence levels mildly increase with expression levels; the Spearman correlation coefficients have min = 0.07, median = 0.23, max = 0.75) across the 15 data sets.

To illustrate the biological information that may be gleaned from further analysis of the results using the confidence scores, we performed a KEGG pathway enrichment analysis over the highest-confidence genes (confidence score = 1) in TCD8 cells of the SKCM (skin cutaneous melanoma) benchmark data set 1 and found 58 significantly (q-value <0.1) enriched biological pathways, including the T-cell receptor signaling pathway, Th1 and Th2 cell differentiation, PD-L1 expression and PD-1 checkpoint pathway in cancer, chemokine signaling pathway, among others. However, after permuting the confidence score profile across genes, none of those findings recur. The rationale to identify the enriched pathways in high-confidence genes is that high-confidence genes are expected to be substantially expressed on specific types of cells, and thus cell type–specific pathways are expected to be enriched in those genes.

We used the package clusterProfiler (12) to identify KEGG pathways that are significantly (q-value <0.1) enriched for high-confidence genes. We found that the high-confidence genes involved in significantly enriched pathways exhibit higher expression. Supplementary Fig. S11 shows that the difference in expression comparing the median expression values of all relevant high-confidence genes to the median expression values of all other genes is significant across all cancer types (two-sided Wilcoxon test, P < 1e−14, Bonferroni correction).

To test the robustness of CODEFACS to variations in cell fraction estimation, we performed deconvolution on two benchmark data sets (SKCM data set 1 and GBM data set 1) by providing the true cell fractions instead of predicted cell fractions as input to CODEFACS. The performance of CODEFACS for the two scenarios is not substantially different, as Supplementary Fig. S12A and S12B shows. Additionally, we revealed that the overall performance of CODEFACS is not strongly affected by random perturbations in cell fractions (Supplementary Fig. S13A and S13B).

In addition to testing on real data sets, we simulated 3,000 cells using R package “dyngen” (13) and performed cell clustering using R package “Rtsne.” Three cell clusters (see Supplementary Fig. S14A and S14B) were defined as proxy cell types. After one fifth of random cells from each cell type were selected for cell-type signature derivation, all other cells were randomly assigned to 16 proxy individuals to generate a bulk mixture data set (1,035 genes and 16 individual samples) with ground truth cell type–specific expression (1,035 genes and 3 cell types) and cell fractions (16 individual samples and 3 cell types). Thereafter the performance of CODEFACS and CIBERSORTx was compared based on this benchmark data set. We conclude that CODEFACS is superior to CIBERSORTx based on this simulated benchmark data set (see Supplementary Figs. S15 and S16).

Finally, to evaluate CODEFACS on real bulk tumor data, we applied it to deconvolve bulk expression data from 21 cancer types (∼8,000 RNA-seq samples) in TCGA. To infer the cellular abundance of each cell type in each sample that is required as input for CODEFACS, we made use of matched bulk methylation data available for these samples and methylation-based reference signature profiles of distinct cell types. These include 11 cell-type signatures [macrophages/dendritic cells—CD14+; B cells—CD19+; CD4+ T cells, CD8+ T cells, T regulatory cells; natural killer (NK) cells—CD56+; endothelial cells, fibroblasts, neutrophils, basophils, eosinophils, and tissue-specific tumor cells] obtained from MethylCIBERSORT (14). Reassuringly, we found high Spearman correlations between the resulting predicted tumor cell fraction and the tumor purity estimates derived from matched mutation and copy-number data (based on ABSOLUTE; ref. 15) for the same samples in each of 10 different cancer types (Spearman correlation: min = 0.72, max = 0.88, avg = 0.8; Supplementary Fig. S17A–S17J). The correlations are consistently high among the comparisons with other reported tumor purity estimates (based on ESTIMATE, LUMP, IHC, and CPE, refs. 16 and 17; Supplementary Figs. S18A–S18N, S19A–S19O, S20A–S20O, and S21A–S21O).

We then asked if CODEFACS can recover the cell type–specific gene expression signature of different cell types in a given cancer type. To this end, we computed the Spearman correlation between (i) the mean deconvolved gene expression of the top confidently deconvolved genes in a given cell type (confidence score ≥0.95) and (ii) the mean expression of these genes, which we derived from completely independent single-cell expression data of the same cancer type (Methods). We find that (i) and (ii) are substantially correlated [Fig. 2I depicts results for the TCGA LUAD (lung adenocarcinoma) data set as an example and Supplementary Fig. S22A–S22H depicts results for the remaining cancer types that have publicly available scRNA-seq data]. The concordance level is higher for cell types that are abundant (e.g., tumor cells and fibroblasts) and decreases for less abundant cell types. Additionally, we observed that tumor cells have the largest fraction of genes whose expression is predicted with high confidence, with the highest in thyroid cancer (∼67.4% of all genes; Supplementary Fig. S23A). Furthermore, seven KEGG pathways are significantly enriched (adjusted P < 0.1) with highly confident genes in tumor cells (confidence score ≥0.95) across the 21 cancer types (Supplementary Fig. S23B). Those pathways mostly involve mRNA surveillance, spliceosome, DNA replication, and mismatch repair.

We additionally assessed the correlation between mean confidence levels for tumor cell expression and the mean aneuploidy scores (data borrowed from ref. 18) across cancer types in TCGA (Supplementary Fig. S24A–S24C), finding that they are not significantly correlated. As an additional quality control test of our pseudobulk mixtures, Supplementary Fig. S25 shows a scatter plot of measured mean gene expression in malignant LUAD single cells versus CODEFACS-predicted gene expression for malignant cells in bulk TCGA LUAD samples.

Inference of Cellular Cross-talk Linked with Enhanced Sensitivity of Mismatch-Repair–Deficient Tumors to Anti–PD-1 Therapy

In normal cells, DNA is constantly repaired in response to DNA damage or DNA replication errors (19). However, defects in specific DNA repair pathways in cancer cells may result in the accumulation of many somatic mutations, resulting in hypermutated tumors [tumor mutation burden (TMB) >10 nonsynonymous mutation/Mb; refs. 20–22]. One cause of hypermutability is MMRD that leads to the accumulation of insertions and deletion (indel) mutations in microsatellite regions of the genome due to uncorrected DNA replication polymerase slippage events. This is known as microsatellite instability (MSI; refs. 23, 24).

Solid tumors with MMRD were shown to be sensitive to ICB therapy irrespective of tumor type, leading the FDA to approve MSI as the first cancer type–agnostic biomarker for patients receiving anti–PD-1 treatment (25). Prior work has led to the prevailing hypothesis that elevated TMB in mismatch-repair deficient tumors leads to more neoantigens, and thus is more likely to activate a host immune response against tumor cells (23, 26–28).

However, not all tumor types with elevated TMB have similar response rates to anti–PD-1 (29, 30), and recent studies have revealed that T cells recognize and respond to only a few neoantigens per tumor (31–34). In the TCGA collection, there are three solid tumor types with a significant association between hypermutability and survival benefit in patients. These three are the MSI-enriched solid tumor types (Supplementary Fig. S26A and S26B, borrowed from refs. 18, 35), MSI tumors highlighted as red dots in Supplementary Fig. S26A (Supplementary Fig. S26C log-rank test P = 0.00084), but not in other tumor types (Supplementary Fig. S26D, log-rank test P = 0.4). Most importantly, immunogenic indel mutations generated by MSI tumors need to escape nonsense-mediated mRNA decay to be presented by the antigen presentation system and activate immune responses (36, 37).

Despite these barriers, MMRD solid tumors have very high response rates to anti–PD-1 therapy. We therefore searched for other tumor–immune interactions that are possibly unique to the TME of MMRD tumors and could contribute to their strong immune responses. One key lead supporting this possibility comes from the work of McGrail and colleagues (38), which shows that MMRD tumors have elevated unfolded protein response (UPR). Recent work has shown that tumor cells with elevated UPR can initiate intercellular communication to spread the UPR to other cells in the TME, including immune cells (39).

To identify cell–cell interactions that are differentially active between MMRD tumors and non-MMRD tumors (red vs. black dots in Supplementary Fig. S26A; data borrowed from ref. 35), we applied CODEFACS to deconvolve the bulk gene expression of all solid tumors from TCGA and integrated their predicted cell type–specific gene expression levels with LIRICS. The top 50 interactions likely to occur in MMRD solid tumors (ordered by FDR adjusted P value) based on LIRICS are shown in a network in Fig. 3A. The complete enrichment table listing results for all plausible ligand–receptor interactions is provided in Supplementary Table S5. An odds ratio >1 implies the interaction is more likely to occur in MMRD solid tumors versus non-MMRD solid tumors; odds ratio <1 implies enrichment in the opposite direction. The top 50 interactions shown in Fig. 3A were tumor type–agnostic and less likely to occur in other hypermutated (but not MMRD) tumors as depicted in Fig. 3B and Supplementary Fig. S27A–S27F.

Figure 3.

A, Interaction network consisting of the top 50 interactions most highly activated in the TME of tumors with DNA MMRD. Interactions highlighted in green represent costimulatory interactions/having an activating effect on the target cell. Interactions highlighted in red represent checkpoint interactions/having an inhibitory effect on the target cell. Interactions highlighted in black represent proinflammatory/chemotaxis interactions involved in inflammatory response and immune cell trafficking to tumor sites. The ligand-expressing cell type is the sender of the signal, and receptor expressing cell type is the receiver. The functional effect of the interaction typically takes place at the receiver's end as shown in immunologic literature. Hence, the functional annotation is indicated with colored sectors on the cell type at the signal receiving end. B, A volcano plot depicting on the x-axis the log2 fold change in the frequency of occurrence of each cell–cell interaction in the TME of hypermutated tumors with an underlying DNA MMRD versus other hypermutated tumors. The y-axis indicates the −log10 FDR-adjusted P value of the observed enrichment. Highlighted in red in the scatter plot are the top 50 interactions that are most differentially active between all MMRD versus non-MMRD tumors (shown in A). C, Scatter plot depicting the relationship between the MMRD cellular cross-talk scores (y-axis) and the nonsynonymous mutation burden (x-axis) across all solid tumor samples in the TCGA. D, Scatter plot depicting the relationship between the MMRD cellular cross-talk scores (y-axis) and the fraction of indel mutations (x-axis) across all solid tumor samples in the TCGA. E, Violin plot depicting the distribution of MMRD cellular cross-talk scores in all MMRD and non-MMRD tumors in the TCGA. F, The estimated effect sizes of each feature when fitting a multivariable logistic regression model to MMRD status data in the TCGA. The x-axis represents the estimated logistic regression coefficient for each feature.

Figure 3.

A, Interaction network consisting of the top 50 interactions most highly activated in the TME of tumors with DNA MMRD. Interactions highlighted in green represent costimulatory interactions/having an activating effect on the target cell. Interactions highlighted in red represent checkpoint interactions/having an inhibitory effect on the target cell. Interactions highlighted in black represent proinflammatory/chemotaxis interactions involved in inflammatory response and immune cell trafficking to tumor sites. The ligand-expressing cell type is the sender of the signal, and receptor expressing cell type is the receiver. The functional effect of the interaction typically takes place at the receiver's end as shown in immunologic literature. Hence, the functional annotation is indicated with colored sectors on the cell type at the signal receiving end. B, A volcano plot depicting on the x-axis the log2 fold change in the frequency of occurrence of each cell–cell interaction in the TME of hypermutated tumors with an underlying DNA MMRD versus other hypermutated tumors. The y-axis indicates the −log10 FDR-adjusted P value of the observed enrichment. Highlighted in red in the scatter plot are the top 50 interactions that are most differentially active between all MMRD versus non-MMRD tumors (shown in A). C, Scatter plot depicting the relationship between the MMRD cellular cross-talk scores (y-axis) and the nonsynonymous mutation burden (x-axis) across all solid tumor samples in the TCGA. D, Scatter plot depicting the relationship between the MMRD cellular cross-talk scores (y-axis) and the fraction of indel mutations (x-axis) across all solid tumor samples in the TCGA. E, Violin plot depicting the distribution of MMRD cellular cross-talk scores in all MMRD and non-MMRD tumors in the TCGA. F, The estimated effect sizes of each feature when fitting a multivariable logistic regression model to MMRD status data in the TCGA. The x-axis represents the estimated logistic regression coefficient for each feature.

Close modal

For each deconvolved solid tumor sample in the TCGA, we additionally computed a sample-specific FDR-weighted interaction enrichment score linked with MMRD cross-talk score (Methods), its nonsynonymous TMB, and fraction of indel mutations, with the goal of comparing their values between MMRD samples (shown as red dots in Fig. 3C and D and Supplementary Fig. S26A) and non-MMRD samples (shown as black dots in Fig. 3C and D and Supplementary Fig. S26A; Methods). The x-axis in Fig. 3C records the nonsynonymous mutation, whereas the x-axis in Fig. 3D records the fraction of indel mutations. The y-axis records those MMRD cross-talk scores. Overall, as expected, the MMRD cross-talk scores are strongly associated with MMRD status (Fig. 3E; Wilcoxon test P < 2.2E−16). However, the associations of these scores with the nonsynonymous TMB or fraction of indel mutations are weak (Fig. 3C, Spearman Rho: 0.06, P = 2.4E−7; Fig. 3D, Spearman Rho: −0.009, P = 0.49).

Finally, we fit a multivariable logistic regression model to the MMRD status data across TCGA using as features the sample-specific MMRD cross-talk scores, the nonsynonymous mutation burden, and the fraction of frameshifted/indel mutations. Figure 3F plots the estimated regression coefficients after transforming the values taken by each feature into Z scores. Overall, we see that the sample-specific MMRD score is significantly associated with MMRD status, independent of other features, thus suggesting that the inferred ligand–receptor interaction repertoire is an additional unique feature of MMRD and could potentially be contributing to the exceptional response rates of MMRD solid tumors to anti–PD-1 therapy.

The top 50 interactions from the differential analysis (Fig. 3A) include the well-established PD-L1–PD-1 checkpoint interaction between tumor cells and CD8+ T cells, but additionally, numerous T-cell costimulatory interactions such as the 41BBL–41BB interaction between tumor cells and CD8+ T cells and ULBP2–NKG2D interaction between tumor cells and CD4+ T cells. They also encompass several chemotaxis-related interactions involved in the trafficking of lymphocytes in and around the tumor mass, such as the CXCL9–CXCR3 chemokine interaction between macrophages and CD4+ T cells and CCL3/4/5–CCR5 interactions between various immune and stromal cell types. These observations are supported by recently emerging results analyzing differences in the TME of 34 MMRD versus 28 mismatch repair–proficient colorectal cancers at a single-cell resolution (40). In this study, the analysis of core cell type–specific gene expression programs upregulated specifically in MMRD tumors revealed elevated expression levels of genes such as PD1, CXCL13, GZMB, and PRF1 in T cells, indicative of chronic T-cell stimulation (41–43) and elevated expression levels of T cell–attracting chemokines in malignant and myeloid cells. For instance, the T cell–attracting chemokine CXCL9 is among the top upregulated genes in cells of myeloid lineage, and Pelka and colleagues (40) suggested that expression of CXCL9, CXCL10, and CXCL11, all of which bind to CXCR3 on T cells, creates a positive feedback loop for increased T-cell activity. Increased CXCR3 binding to its ligands has also been previously associated with response to anti–PD-1 treatment (44). These findings overlap with our differential analysis using LIRICS, wherein we observe a significant enrichment of CXCL9–CXCR3 chemokine interaction between macrophages and CD4+ T cells in MMRD solid tumors across the TCGA collection (fold change: 2.58, FDR-adjusted P = 2.38 × 10−10, displayed in Fig. 3A among the top 50 most significantly enriched ligand–receptor interactions), showing that this phenomenon is not specific to colorectal MMRD tumors.

Machine Learning–Guided Discovery of Cell Type–Specific Ligand–Receptor Interactions that Predict Patient Response to ICB Therapy

It has been previously demonstrated in pan-cancer clinical trials that differences in mutation burden result in differences in survival upon receiving ICB treatments, leading to the recent FDA approval of TMB as a biomarker for anti–PD-1 therapy (45). However, open questions remain on the contribution of additional factors within individual cancer types (46, 47). Hence, to find additional tumor–immune interactions predictive of response to ICB therapy, we chose to focus on studying TME differences between hypermutated tumors (>10 nonsynonymous mutations/Mb) and nonhypermutated tumors in TCGA data. We focus on melanoma, which has one of the best response rates to ICB, and where there are many independent publicly available bulk expression data sets of patients receiving anti–PD-1 treatment. Starting from the deconvolved TCGA–SKCM data set as our training set (N = 469), we trained a genetic algorithm to find cell type–specific ligand–receptor interactions that are predictive of hypermutation in melanomas (see Methods; Fig. 4A). For all analyses in this study, the definition of hypermutated tumor was set to any tumor with >10 nonsynonymous mutations per megabase based on previous papers (20, 21, 48). In melanoma, these excess mutations primarily arise from UV radiation–induced DNA damage and not replication slippage and MMRD (21). We term the interactions identified by the algorithm in this pretext task as melanoma hypermutation-specific functional interactions (MSFI), and the network formed by these interactions is displayed in Fig. 4B.

Figure 4.

A, Overview of the machine learning analysis used to identify cell type–specific interactions that are predictive of response to ICB therapy. B, A chord diagram of the resulting MSFI network. Each individual interaction is represented by a link from the source cell type (ligand-expressing cell type) to the target cell type (receptor-expressing cell type) and the color of the link represents the color of the source cell type. The ligand-expressing cell type is the sender of the signal and the receptor-expressing cell type is the receiver. The functional effect of the interaction typically takes place at the receiver's end as shown in immunologic literature. Hence, the functional annotation is indicated with colored sectors on the cell type at the signal receiving end. For interactions that are activating/costimulatory, the sector in the corresponding target cell type is highlighted in green. For inhibitory/checkpoint interactions, the sector in the target cell type is highlighted in red. Interactions involved in chemotaxis are highlighted in black, those mediating a proinflammatory response are highlighted in blue, and cell-adhesion interactions are highlighted in gray. C, Kaplan–Meier plot depicting the progression-free survival of the combined set of patients with melanoma receiving ICB (N = 244). The patients are stratified into low-risk/high-risk groups based on the median value of MSFI score. D, Kaplan–Meier plots depicting the overall survival of all patients with melanoma receiving ICB (N = 244). The patients are stratified into low-risk/high-risk groups based on the median value of MSFI score. E and F, Estimated hazard ratios of MSFI scores and TMB on progression-free survival and overall survival based on the multivariable Cox proportional hazards model. E, Estimates of Cox proportional hazards ratios for MSFI score and TMB for a multivariable model fitted to progression-free survival outcome data. F, Estimates of Cox proportional hazards ratios for MSFI score and TMB for a multivariable model fitted to overall survival outcome data.

Figure 4.

A, Overview of the machine learning analysis used to identify cell type–specific interactions that are predictive of response to ICB therapy. B, A chord diagram of the resulting MSFI network. Each individual interaction is represented by a link from the source cell type (ligand-expressing cell type) to the target cell type (receptor-expressing cell type) and the color of the link represents the color of the source cell type. The ligand-expressing cell type is the sender of the signal and the receptor-expressing cell type is the receiver. The functional effect of the interaction typically takes place at the receiver's end as shown in immunologic literature. Hence, the functional annotation is indicated with colored sectors on the cell type at the signal receiving end. For interactions that are activating/costimulatory, the sector in the corresponding target cell type is highlighted in green. For inhibitory/checkpoint interactions, the sector in the target cell type is highlighted in red. Interactions involved in chemotaxis are highlighted in black, those mediating a proinflammatory response are highlighted in blue, and cell-adhesion interactions are highlighted in gray. C, Kaplan–Meier plot depicting the progression-free survival of the combined set of patients with melanoma receiving ICB (N = 244). The patients are stratified into low-risk/high-risk groups based on the median value of MSFI score. D, Kaplan–Meier plots depicting the overall survival of all patients with melanoma receiving ICB (N = 244). The patients are stratified into low-risk/high-risk groups based on the median value of MSFI score. E and F, Estimated hazard ratios of MSFI scores and TMB on progression-free survival and overall survival based on the multivariable Cox proportional hazards model. E, Estimates of Cox proportional hazards ratios for MSFI score and TMB for a multivariable model fitted to progression-free survival outcome data. F, Estimates of Cox proportional hazards ratios for MSFI score and TMB for a multivariable model fitted to overall survival outcome data.

Close modal

Having identified the MSFIs, we applied CODEFACS to deconvolve the bulk expression data of pretreatment samples from the three largest publicly available melanoma data sets where patients received anti–PD-1 treatment [either monotherapy or in combination with anti–cytotoxic T-lymphocyte associated protein 4 (CTLA4); Methods; refs. 49–51]. We then used LIRICS (steps 1 and 2) to the respective deconvolved expression of each of these checkpoint data sets and, without any additional training, simply quantified the number of MSFIs that are active in each of these patients’ tumor samples, which we denote as the tumor's MSFI score. Remarkably, we find that the MSFI score of each sample can robustly stratify patients into those who are likely to respond to ICB versus those who are unlikely to respond. Figure 4C plots progression-free survival differences between the low-risk group (MSFI score > median value; N = 120) and the high-risk group (MSFI score ≤ median value; N = 124); log-rank test P = 0.00057. Figure 4D plots overall survival differences between the same two groups; log-rank test P = 0.0031. Supplementary Fig. S28A–S28H plots the progression-free survival and overall survival differences between the same two groups for each data set separately. Importantly, a multivariable Cox proportional hazards model with MSFI scores and TMB of each patient receiving anti–PD-1 (wherever TMB data were available) shows that the MSFI score remains significantly associated with improved progression-free survival (P = 0.013) and overall survival (P = 0.0258) after accounting for differences in mutation burden (Fig. 4E and F), suggesting that the genetic algorithm learns how to predict immunotherapy outcomes without overfitting to TMB.

In addition, we plotted Kaplan–Meier survival curves depicting the stratification performance of seven other previously published bulk transcriptomics-based scores (52–60) besides MSFI, based on their respective median values (Figs. 5AH and 6A–H). For TIDE (Tumor Immune Dysfunction and Exclusion; ref. 59) and MPS (Melanocytic Plasticity Signature; ref. 60) scores, a patient was put into the high-risk group if that patient's TIDE and MPS scores were greater than respective median values. This is because these scores were previously linked with immune resistance. Overall, for patients receiving anti–PD-1 monotherapy, the MSFI score achieved the best survival stratification results, whereas for patients who are receiving or have previously received anti-CTLA4 in addition to anti–PD-1, the immune signature score (58) performed best.

Figure 5.

Survival stratification performance of MSFI score versus seven other previously published bulk transcriptomics-based scores for patients with melanoma receiving anti–PD-1 only. A, Kaplan–Meier plots showing progression-free survival and overall survival differences between the low-risk groups defined by the median value of the MSFI score. B–H, Kaplan–Meier plots showing progression-free and overall survival differences between low- and high-risk groups defined by the median value of previously published bulk transcriptomics-based signatures. The significance of survival differences was estimated using the log-rank test. Time on the x-axis is measured in days. The signatures evaluated in each panel are: B, IMPRES (8); C, melanocytic plasticity signature (9); D, TIDE (10); E, IFNG signature (11); F, T cell–inflamed GEP (12); G, cytotoxic signature (13); H, immune signature (14).

Figure 5.

Survival stratification performance of MSFI score versus seven other previously published bulk transcriptomics-based scores for patients with melanoma receiving anti–PD-1 only. A, Kaplan–Meier plots showing progression-free survival and overall survival differences between the low-risk groups defined by the median value of the MSFI score. B–H, Kaplan–Meier plots showing progression-free and overall survival differences between low- and high-risk groups defined by the median value of previously published bulk transcriptomics-based signatures. The significance of survival differences was estimated using the log-rank test. Time on the x-axis is measured in days. The signatures evaluated in each panel are: B, IMPRES (8); C, melanocytic plasticity signature (9); D, TIDE (10); E, IFNG signature (11); F, T cell–inflamed GEP (12); G, cytotoxic signature (13); H, immune signature (14).

Close modal
Figure 6.

Survival stratification performance of MSFI score versus seven other previously published bulk transcriptomics-based markers for patients with melanoma receiving anti-CTLA4 + anti–PD-1. A, Kaplan–Meier plots showing progression-free survival and overall survival differences between the low-risk groups defined by the median value of the MSFI score. B–H, Kaplan–Meier plots showing progression-free and overall survival differences between low- and high-risk groups defined by the median value of previously published bulk transcriptomics-based signatures. The significance of survival differences was estimated using the log-rank test. Time on the x-axis is measured in days. The signatures evaluated in each panel are: B, IMPRES (8); C, melanocytic plasticity signature (9); D, TIDE (10); E, IFNG signature (11); F, T cell–inflamed GEP (12); G, cytotoxic signature (13); H, immune signature (14).

Figure 6.

Survival stratification performance of MSFI score versus seven other previously published bulk transcriptomics-based markers for patients with melanoma receiving anti-CTLA4 + anti–PD-1. A, Kaplan–Meier plots showing progression-free survival and overall survival differences between the low-risk groups defined by the median value of the MSFI score. B–H, Kaplan–Meier plots showing progression-free and overall survival differences between low- and high-risk groups defined by the median value of previously published bulk transcriptomics-based signatures. The significance of survival differences was estimated using the log-rank test. Time on the x-axis is measured in days. The signatures evaluated in each panel are: B, IMPRES (8); C, melanocytic plasticity signature (9); D, TIDE (10); E, IFNG signature (11); F, T cell–inflamed GEP (12); G, cytotoxic signature (13); H, immune signature (14).

Close modal

We further report the performance of the MSFI score in classifying (partial/complete responses vs. stable/progressive) disease as determined by posttreatment radiologic review. The MSFI score achieves an average AUC of 0.63 for this task across different data sets and treatment groups (the AUCs for each data set and treatment group are reported in Supplementary Table S6). A similar performance could not be achieved if the placement of the ligand and receptor between interacting cell types in the MSFI network was swapped (average AUC, 0.58) or by randomly shuffling the interaction activity profiles (average AUC ∼0.5), testifying that the selected cell–cell interactions are best predictive of response to ICB. AUCs for the seven other scores we benchmarked against are also provided in Supplementary Table S6. Overall, for this task, cytotoxic signature score (55) achieved the best average performance (average AUC: 0.68).

We note that the performance levels of some previously published scores on more recent larger bulk expression data sets, where the original RNA-seq reads have been uniformly aligned and normalized as described in the Methods, are different from those expected from the original publications (Supplementary Table S6), pointing to the growing awareness of the sensitivity of expression-based predictors to batch effects, the normalization and alignment methods used, highlighting the need to preprocess the raw sequencing data uniformly in making meaningful comparisons (see Discussion).

We next examined the interactions in the MSFI network (Fig. 4B). Supplementary Table S7 reports the number of times each interaction from the MSFI network was selected by the genetic algorithm over multiple training runs on the TCGA data (Methods). Overall, we find an overrepresentation of cell type–specific costimulatory interactions from immunologic literature (hypergeometric test P < 0.05; see Methods; refs. 61–71), highlighting the importance of costimulation in mediating successful antitumor–immune responses. Additionally, the MSFI network includes various cytokine–chemokine interactions involved in the trafficking of specific NK, T-, and B-cell subsets to the tumor site.

Taken together, these results suggest that although recognition of melanoma tumor neoantigens serves as an initiating trigger, the activation of additional, cell type–specific, costimulatory signals in the TME is required to further enhance an effective host immune response upon ICB treatments.

This study presents a new computational tool, CODEFACS, and a supporting immune interactions framework, LIRICS, that enable an (averaged) “virtual single-cell” characterization of the TME from bulk tumor expression data. Applying these tools across the TCGA, we systematically identify, for the first time, cell type–specific ligand–receptor pairs that are likely to interact directly in the TME of specific patient populations and prioritize those that are associated with patients’ response to ICB. Specifically, we identified a shared core of intercellular TME interactions unique to DNA MMRD solid tumors that are associated with improved patient survival and enhanced sensitivity to ICB therapy versus other highly mutated tumors. Finally, focusing on melanoma, we show that one can bootstrap on the large deconvolved data resource from TCGA using machine learning techniques to discover cell–cell interactions within the TME that successfully stratify patients’ responses to ICB.

The heightened cellular cross-talk unique to the TME of MMRD tumors suggests that T cells are being actively recruited to the tumor mass and activated by multiple costimulatory signals, in addition to immunogenic neoantigens, only to be kept in balance by other immunoregulatory mechanisms such as the PD-1–PD-L1 checkpoint interaction between CD8+ T cells and tumor cells. This suggests that when this interaction is blocked by anti–PD-1 treatment, additional costimulatory interactions such as the 41BBL–41BB or ULBP2–NKG2D interaction between tumor and T cells are likely to lead to the observed enhanced response of MMRD tumors to ICB therapy. These results further support the investigation of the UPR in tumor cells or infiltration of specific immune cell subsets such as NKG2D+, 41BB+, and PD-1+ T cells in the TME of patients. Notably, recent preclinical studies have shown that combination therapies dually targeting specific immune-checkpoint and costimulatory interactions rescue antitumor–immune responses in low TMB and highly immunosuppressive settings (72–77). Currently, several clinical trials to assess the safety and efficacy of such combinations in humans are in progress (78–81).

Although we have provided a tool kit to prioritize clinically relevant cell–cell interactions from bulk tumor expression, there are limitations that should be noted and potentially improved upon in the future. First, several recent publications have reported discrepancies between different RNA-seq expression quantification pipelines based on the reference transcriptome version used, choice of method, and normalization (e.g., alignment-based vs. alignment-free; refs. 82–86). This can potentially affect the reproducibility of findings. Hence, whenever possible, we recommend that all bulk RNA-seq data sets are homogeneously preprocessed using the same RNA-seq quantification method and reference transcriptome before the application of our tools. Indeed, following this notion, in this study we have preprocessed all bulk RNA-seq data sets in a uniform manner, using STAR (v2.7.6a) + RSEM (v1.3.3) and GENCODE v23 human genome annotation (Methods).

Second, CODEFACS itself has several limitations. (i) It requires prior information about the cell-type composition of the input tumors, or, alternatively, knowledge of the pertaining cell types’ gene expression or methylation signatures that can be used to infer their abundances, and its accuracy depends on the accuracy of the latter. (ii) It is predictive typically for cell types with higher abundance, and its performance deteriorates for low-abundance cell types. However, the confidence scores provided help by allowing the user to rank genes in each cell type by the quality of predictions. During the revision of our manuscript, we became aware of another software package called BayesPrism (87) that solves the problem of estimating gene expression in malignant cells taking as input read counts from a bulk mixture of malignant and nonmalignant cells. BayesPrism uses a Bayesian approach and Markov chain Monte Carlo methods, which are completely different from the methods in CODEFACS. We found that the results of BayesPrism are comparable to CODEFACS, but very dependent on the input sample size. Because the problems solved are different and because of the BayesPrism dependence on sample size, we elected not to present a full head-to-head comparison.

Third, LIRICS is currently restricted to well-defined protein–protein interactions and does not consider the spatial localization of cells in the TME. The inclusion of the latter with the advent of forthcoming spatial transcriptomics data is likely to lead to considerably more informative interaction inference approaches. In Supplementary Note 2, we have described how LIRICS fundamentally differs from previously published single cell–based approaches in terms of both its inputs and the problem being solved. We have further provided a brief walkthrough of how interested users can use LIRICS to conduct a meta-enrichment analysis of their own (as we demonstrated with the TCGA) with the example of the three large melanoma ICB data sets studied in this work. In case one would like to do a more exploratory analysis, we provide a walkthrough of how this can be achieved via a joint analysis of single-cell and corresponding deconvolved bulk data.

Although this work focuses on studying the TME, the tools presented here can be applied to prioritize important cell–cell interactions in noncancerous tissues under a variety of normal and disease states. One interesting application that we envision is the characterization of clinically relevant intercellular interactions occurring at the maternal–fetal interface using corresponding bulk gene expression data and pregnancy outcome information, whose elucidation may help treat and mitigate preeclampsia and other pregnancy-related complications. One can also use our tools to study bulk gene expression data from premalignant tissue samples and compare them against malignant samples to elucidate cell–cell interaction dynamics on the way to malignancy. Finally, one can deconvolve expression data from autoimmune disorders to learn more about the underlying immune interactions. In sum, the computational tools developed and presented here offer a cost-effective way to study immune responses at a cell type–specific resolution from the widely abundant bulk gene expression in both health and disease, complementing first-line single-cell technologies in situations where the latter are less applicable.

Data Collection and Preprocessing

scRNA-seq data sets

To benchmark the performance of CODEFACS, we first set out to obtain publicly available scRNA-seq data sets where both tumor and nontumor cells were successfully isolated. This search led us to the identification of nine such scRNA-seq data sets from the literature, each from a different cancer type. Collection of additional single-cell data sets was frozen after December 2019. Details of the single-cell data sets that we curated for this study are in Supplementary Table S8. For each data set sequenced on the SmartSeq2 platform, the log-normalized transcript counts for each gene in each sequenced cell were made publicly available by the original authors. For the application of deconvolution, these counts were transformed back to the transcripts per million (TPM) scale. For data sets sequenced on the 10X platform, unique molecular identifier counts for each gene were made publicly available and were scaled by the library size of each cell and multiplied by a factor of 1 million to get expression values in the TPM scale.

Summary of the CODEFACS algorithm and confidence ranking system rationale

The CODEFACS algorithm is a greedy algorithm that is executed in three stages: (i) high-resolution deconvolution with recursive splitting (module 1), (ii) hierarchical deconvolution (module 2), and (iii) imputation (module 3). Each of these stages aims to accurately estimate a gene's expression in a given cell type and sample based on some assumptions.

The fundamental rationale behind the confidence ranking system is that we expect to have high confidence in the estimated expression of a gene in a cell type if its expression distribution across samples satisfies as many stage-specific modeling assumptions as possible. This intuition is captured by computing multiple rankings of gene–cell type pairs based on different stage-specific confidence features that have been engineered (see Supplementary Table S9), and subsequently averaging over these rankings.

The aim of the confidence ranking cutoffs defined at each stage was to help the CODEFACS algorithm heuristically converge to a solution by forcing a distinction between gene–cell type pairs whose estimated expression levels are satisfactory (i.e., belonging to a high-confidence set) or require further refinement (i.e., belonging to a low-confidence set). The “greediness” of the algorithm arises from the fact that we only focus on refining the expression of gene–cell type pairs belonging to the low-confidence set and that we choose cutoffs as conservatively as possible. After execution of all stages, a single confidence score is assigned to each gene–cell type pair. The sequential confidence-building design of our approach is inspired, for example, by the quality scores that are assigned to reads generated from sequencing machines. The final confidence score (0–1) indicates the confidence level of each predicted gene in each cell type. We tend to assign higher confidence scores to those genes in a specific cell type whose predicted expression levels are highly correlated with the high-confidence gene set. A detailed description is provided in Supplementary Note 1.

Guide to Some of the Supplementary Information Concerning Methods of CODEFACS

The methods for CODEFACS and LIRICS are complex and have many stages. At the suggestion of a reviewer, we provide here some cross-reference to the Supplementary Information that should aid the reader to find specific information about different stages of the algorithms and some additional evaluations of performance not covered in Supplementary Figs. S1–S28. Supplementary Note 1 describes in detail the algorithms in CODEFACS. A schematic description of CODEFACS module 1 can be found in Supplementary Fig. S29. An evaluation of how confidence of gene expression predictions changes from module 1 to module 2 to module 3 can be found in Supplementary Fig. S30A–S30O. Supplementary Fig. S31 shows a schematic of the optional cell fraction estimation step in CODEFACS module 1. Supplementary Fig. S32 shows a schematic of the optional batch correction step in CODEFACS module 1. Supplementary Fig. S33 shows a schematic of the recursive splitting method in CODEFACS module 1.

Supplementary Note 1, part 3, delves further into the confidence ranking system, which relies on the assumption that the expression levels of functionally related genes in the same cell type are correlated and that this correlation can be used incrementally to boost confidence in the prediction of gene expression levels within a cell type. The correlation of gene expression levels of different genes within a cell type is quantified in Supplementary Fig. S34A–S34C.

Supplementary Note 1, part 4, describes CODEFACS module 2, which is illustrated schematically in Supplementary Fig. S35. Part 5 describes how confidence values are changed in module 2. Part 6 describes CODEFACS module 3, which is illustrated schematically in Supplementary Fig. S36. Part 7 describes how confidence values are changed in CODEFACS module 3. Part 8 summarizes the three modules together and includes Supplementary Fig. S37A–S37F that evaluates some of the algorithmic choices in the recursive splitting and sliding window steps of module 1. Part 8 also includes Supplementary Fig. S38A–S38O that quantifies how gene expression prediction performance improves from module 1 to module 2 to module 3.

Bulk RNA-seq data sets

Gene expression and matching bulk tumor methylation data from fresh-frozen tumor biopsies in TCGA were downloaded from http://xena.ucsc.edu/ (88). In addition, publicly available bulk expression data from formalin-fixed paraffin-embedded tumor biopsies of patients with melanoma receiving ICB treatment were downloaded from refs. 49–51. All bulk RNA-seq data sets were collected such that they have a sufficiently large sample size to reliably perform complete deconvolution of expression profiles (>4 times the number of cell types involved; ref. 8). Collection of data sets was frozen after December 2019. Details on bulk RNA-seq data sets deconvolved in this study can be found in Supplementary Table S10. To maintain consistency with the pipeline used for preprocessing TCGA data, bulk gene expression levels in ICB data sets were requantified using STAR v2.7.6a and RSEM v1.3.3 (89) with GENCODE v23 human genome annotation (90). Furthermore, to mitigate technical biases, between-sample scaling factors were estimated using the TMM method implemented in edgeR (91), and TPM values in each sample were further rescaled by these scaling factors (92).

Bulk RNA-seq data sets Used for benchmarking

To evaluate the performance of CODEFACS, we generated different pseudobulk RNA-seq data sets from mixing experiments with single-cell data. Each sample in each benchmark data set has matching cell type–specific gene expression profiles derived from averaging scRNA-seq profiles of individual cells from the same sample and cell type. For each pseudo-bulk sample, to derive the ground-truth cell fraction, we normalized the number of cells for each cell type by dividing by the total number of cells. These profiles serve as the ground truth for the evaluation of deconvolution performance. To avoid any circularity in our validations, for each of the single-cell data sets involved, single-cell data from four randomly chosen patients were separated from the rest. These data were used to derive reference gene expression signatures for each cell type. The mixing experiments were then performed on single-cell data of the remaining patients who were hidden from the reference signature derivation process. The data sets and the mixing experiment used to generate the pseudobulk samples are listed in Supplementary Table S4. In addition, we obtained a FACS lung cancer data set, which includes purified RNA-seq for four cell types (11) and matched bulk RNA-seq for each sample.

To inject technical noise (e.g., PCR amplification bias or degradation of mRNA in preserved biopsies; ref. 92) into the pseudobulk data sets, the following general procedure was followed:

For each bulk sample,

  1. Multiply the TPM count values of a random subset of genes with a random nonnegative scalar value.

  2. Renormalize the count values of all genes in that sample to add up to a million (thereby satisfying the definition of TPM).

With the general procedure described above, here is how we generated each noisy replicate data set with batch effects:

  • For each sample in SKCM data set 2, SKCM data set 5, GBM data set 2, and GBM data set 5: artificially perturb the TPM values of 500 randomly chosen genes with randomly chosen multiplicative scalars in the range 50–100 to generate SKCM data set 3, SKCM data set 6, GBM data set 3, and GBM data set 6, respectively.

  • For each sample in SKCM data set 2, SKCM data set 5, GBM data set 2, and GBM data set 5: artificially perturb the TPM values of 3,000 randomly chosen genes with randomly chosen multiplicative scalars in the range 0–100 to generate SKCM data set 4, SKCM data set 7, GBM data set 4, and GBM data set 7, respectively.

For reproducibility, we fixed the seed in our implementations. The codes to generate these data sets are made available in the code repository. Additionally, we also generated data sets with noise injected into the cell fraction estimates and found that the overall performance of CODEFACS is not strongly affected by random perturbations in cell fractions (Supplementary Fig. S13A and S13B).

Deconvolving bulk RNA-seq data sets Used for benchmarking with CIBERSORTx

To perform high-resolution deconvolution of bulk expression data sets used in benchmarking experiments using CIBERSORTx, we relied on their online portal (https://cibersortx.stanford.edu/; ref. 8). We ran their impute cell expression module in high-resolution analysis mode. To run this analysis, three inputs are required: bulk gene expression profiles (Mixture file), the cell type–specific signature matrix (Signature Matrix file), and the names of genes whose expression needed to be deconvolved (we fed the names of all genes). In addition, batch correction mode (B-mode) was checkmarked to adjust for any batch effects between the signature expression matrix and bulk expression profiles. The outputs obtained are high-resolution cell type–specific expression values of each gene of interest in each cell type and sample.

Deconvolving ICB and TCGA bulk RNA-seq data sets Using CODEFACS

For the application of CODEFACS, molecular profiles of signature genes of each cell type of interest are needed to estimate the relative cell fractions in the bulk. We used single-cell expression–derived signatures as priors to deconvolve the melanoma ICB data sets. To derive these signatures from single-cell data, we first obtained the class labels of each cell type of interest. These data are publicly available for each single-cell data set we collected (9, 10, 93–98). Hence, we primarily used these labels in our study (unless further refinement of labels into specific cell subtypes of interest was needed for a specific usage). With a collection of single-cell expression profiles and matching cell-type labels as input, we used the CIBERSORT online tool to derive a cell type–specific signature matrix. Thereafter, we applied CODEFACS to ICB data sets with default parameter settings and batch correction requirement specified. For TCGA deconvolution, we first estimated cell fractions based on bulk methylation and then applied CODEFACS to corresponding bulk gene expression for the 21 cancer types that had both types of data available. The main motivation to use DNA methylation to generate initial cell fraction estimates in TCGA is that the methylation signatures from MethylCIBERSORT (14) cover exactly the same set of cell types for all 21 cancer types in TCGA we analyzed (see Supplementary Table S11). In contrast, cell types covered by expression-based signatures (derived from publicly available single-cell data) across cancer types are very heterogeneous due to design differences of the original studies. In addition, DNA methylation–based signatures are considered to be more stable marks of cellular identity compared with dynamic RNA expression–derived signatures (99). The methylation-based cell-type signatures were obtained from MethylCIBERSORT (14). We applied CODEFACS to TCGA data sets with default parameter settings and without batch correction option specified. For the details of CODEFACS algorithm and hyperparameter settings, see Supplementary Note 1.

Guide to some of the supplementary information concerning methods of LIRICS

Supplementary Note 2 describes the methods of LIRICS and includes a tutorial related to the MMRD analysis in Results above. We begin with a description of the literature search and curation of the ligand–receptor data used in LIRICS. Second, we summarize LIRCS steps 1, 2, and 3. Third, we compare the problems that LIRICS solves to what can be done with other methods such as CellPhoneDb. Fourth and finally, we give a detailed tutorial whose two key outputs of predicted ligand–receptor interactions are illustrated in Supplementary Figs. S39A, S39B, S40A, and S40B.

Defining sample-specific MMRD cross-talk scores

The MMRD-specific ligand–receptor interaction enrichment score for a given tumor sample j was defined as

formula

where 1ij is an indicator variable that takes value 1 if ligand–receptor interaction i tends to occur with greater frequency in MMRD tumors and is observed as “activated” in the given sample j (0 otherwise; see LIRICS step 2).

Machine Learning for Discovery of Cell–Cell Interactions Predictive of Clinical Responses to ICB

The genetic algorithm is a randomized heuristic search algorithm designed to select optimal features for a prediction task given some user-defined fitness function for training (100). For the genetic algorithm, the training set is from TCGA, and the test sets are the three ICB data sets. Overall, the features of interest are ligand–receptor interactions between cell types, and the prediction task is predicting response to ICB treatment. The fitness function guides the search for predictive cell–cell interaction features, which predict the hypermutation status of a sample (TMB >10 mutations per Mb) based on the total number of these features that are observed as activated in a given (nontreated, TCGA melanoma) sample. Prediction accuracy is quantified by the AUC. The overall feature space consists of 4,288 possible ligand–receptor interactions between 10 cell types in the melanoma TME. To minimize the risk of overfitting to the TCGA data set and to aid in faster convergence of the genetic algorithm during training, the size of the feature space is first reduced by filtering out complex and harder-to-measure interactions involving receptors/ligands encoded by multiple genes, interactions among deconvolved cell types that are not shared between train and test data sets, and interactions with fold change in frequency of activation between hypermutated versus nonhypermutated TCGA melanoma samples ≤1. A total of 769 ligand–receptor interactions passed these filters and compose the search space of the genetic algorithm.

The genetic algorithm then proceeds by randomly selecting subsets of features from the set of 769 interactions, defined as the seed population. These subsets are iteratively modified to maximize the fitness function specified above, following standard evolutionary algorithms practices, adding or deleting features aiming to gradually increase the fitness function until a (local) optimum is converged upon. Because the fitness function landscape is often nonconvex and the training process is stochastic, we repeat the training process 500 times, each with a randomly chosen seed population. We then chose the most frequently selected features over all 500 solutions to approximate a solution that was likely to be closer to the global optimum solution and less likely to overfit. For this aggregate feature selection, one can estimate the empirical P value associated with each feature being selected more than the observed number of times by random chance. This is done by repeatedly shuffling 1,000 times the values of the binary solution vector returned by the genetic algorithm for each training run. Given 1,000 random permutations of 500 solutions, one can then derive a null distribution representing the number of times a feature is selected across all 500 solutions by random chance. The empirical P value for each feature j is then estimated as

formula

We found that features with empirical P < 0.01 were selected >100 times by the genetic algorithm. Overall, 103 interactions were selected >100 times and constitute the final MSFI network displayed in Fig. 4B. The downstream results are qualitatively similar for more stringent thresholds. The genetic algorithm was implemented in R using the genalg package using the following default hyperparameter settings: (population size = 200, number of training epochs = 100, mutation probability = 0.001, elitism = true). A more detailed description of the pros and cons of this approach versus other feature selection approaches has been provided in Supplementary Note 2.

Hypergeometric Test for Overrepresentation of Costimulatory Interactions in MSFI Network

The overrepresentation P value for costimulatory interactions is calculated as follows: Let S = 103 be the number of MSFI interactions selected by the genetic algorithm. The foreground F = 103 is the number of all plausible interactions between deconvolved cell types common to both train and test sets and costimulatory. The background B = 2004 − 103 is the number of all plausible interactions between deconvolved cell types common to both train and test sets but not costimulatory. The overlap O between the foreground and the selection is = 11. The hypergeometric enrichment P value for having ≥O costimulatory interactions selected by random chance is then calculated as 1 – phyper(O −1, F, B, S). The phyper function is implemented in the R stats package.

Data and Code Availability

All deconvolved expression data generated in this study, codes, and tools (CODEFACS and LIRICS) are available via the ZENODO repository: https://zenodo.org/record/5790343.

K. Wang, S. Patkar, J.S. Lee, A.A. Schäffer, and E. Ruppin report a patent pending for methods of CODEFACS and LIRICS. No disclosures were reported by the other authors.

K. Wang: Conceptualization, data curation, software, formal analysis, visualization, methodology, writing–original draft. S. Patkar: Conceptualization, data curation, software, formal analysis, visualization, methodology, writing–original draft. J.S. Lee: Data curation, formal analysis, writing–review and editing. E. Gertz: Software, formal analysis, methodology, writing–review and editing. W. Robinson: Data curation. F. Schischlik: Visualization, writing–review and editing. D.R. Crawford: Data curation. A.A. Schäffer: Conceptualization, formal analysis, supervision, methodology, writing–original draft, project administration, writing–review and editing. E. Ruppin: Conceptualization, supervision, funding acquisition, methodology, writing–original draft, project administration, writing–review and editing.

This research is supported in part by the Intramural Research Program of the NIH, NCI, Center for Cancer Research. This work utilized the computational resources of the NIH HPC Biowulf cluster. The results here are in part based upon data generated by the TCGA Research Network: https://www.cancer.gov/tcga. We would additionally like to acknowledge Dr. Noam Auslander, Dr. Chi-Ping Day, Dr. Di Wu, and other members of the Cancer Data Science Laboratory for their helpful feedback on this work.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

Note: Supplementary data for this article are available at Cancer Discovery Online (http://cancerdiscovery.aacrjournals.org/).

1.
Paget
S
.
The distribution of secondary growths in cancer of the breast
.
Lancet
1889
;
133
:
571
3
.
2.
Smyth
MJ
,
Ngiow
SF
,
Ribas
A
,
Teng
MW
.
Combination cancer immunotherapies tailored to the tumour microenvironment
.
Nat Rev Clin Oncol
2016
;
13
:
143
58
.
3.
Wagner
A
,
Regev
A
,
Yosef
N
.
Revealing the vectors of cellular identity with single-cell genomics
.
Nat Biotechnol
2016
;
34
:
1145
60
.
4.
Candès
EJ
,
Wakin
MB
.
An introduction to compressive sampling
.
IEEE Sig Proc Mag
2008
;
25
:
21
30
.
5.
Ahn
J
,
Yuan
Y
,
Parmigiani
G
,
Suraokar
MB
,
Diao
LX
,
Wistuba
II
et al
.
DeMix: deconvolution for mixed cancer transcriptomes using raw measured data
.
Bioinformatics
2013
;
29
:
1865
71
.
6.
Quon
G
,
Haider
S
,
Deshwar
AG
,
Cui
A
,
Boutros
PC
,
Morris
Q
.
Computational purification of individual tumor gene expression profiles leads to significant improvements in prognostic prediction
.
Genome Med
2013
;
5
:
29
.
7.
Fox
NS
,
Haider
S
,
Harris
AL
,
Boutros
PC
.
Landscape of transcriptomic interactions between breast cancer and its microenvironment
.
Nat Commun
2019
;
10
:
3116
.
8.
Newman
AM
,
Steen
CB
,
Liu
CL
,
Gentles
AJ
,
Chaudhuri
AA
,
Scherer
F
et al
.
Determining cell type abundance and expression from bulk tissues with digital cytometry
.
Nat Biotechnol
2019
;
37
:
773
82
.
9.
Jerby-Arnon
L
,
Shah
P
,
Cuoco
MS
,
Rodman
C
,
Su
MJ
,
Melms
JC
et al
.
A cancer cell program promotes T cell exclusion and resistance to checkpoint blockade
.
Cell
2018
;
175
:
984
97
.
10.
Neftel
C
,
Laffy
J
,
Filbin
MG
,
Hara
T
,
Shore
ME
,
Rahme
GJ
et al
.
An integrative model of cellular states, plasticity, and genetics for glioblastoma
.
Cell
2019
;
178
:
835
49
.
11.
Gentles
AJ
,
Hui
ABY
,
Feng
WG
,
Azizi
A
,
Nair
RV
,
Bouchard
G
et al
.
A human lung tumor microenvironment interactome identifies clinically relevant cell-type cross-talk
.
Genome Biol
2020
;
21
:
107
.
12.
Yu
GC
,
Wang
LG
,
Han
YY
,
He
QY
.
clusterProfiler: an R package for comparing biological themes among gene clusters
.
OMICS
2012
;
16
:
284
7
.
13.
Cannoodt
R
,
Saelens
W
,
Deconinck
L
,
Saeys
Y
.
Spearheading future omics analyses using dyngen, a multi-modal simulator of single cells
.
Nat Commun
2021
;
12
:
3942
.
14.
Chakrayarthy
A
,
Furness
A
,
Joshi
K
,
Ghorani
E
,
Ford
K
,
Ward
MJ
et al
.
Pan-cancer deconvolution of tumour composition using DNA methylation
.
Nat Commun
2018
;
9
:
3220
.
15.
Carter
SL
,
Cibulskis
K
,
Helman
E
,
McKenna
A
,
Shen
H
,
Zack
T
et al
.
Absolute quantification of somatic DNA alterations in human cancer
.
Nat Biotechnol
2012
;
30
:
413
21
.
16.
Yoshihara
K
,
Shahmoradgoli
M
,
Martinez
E
,
Vegesna
R
,
Kim
H
,
Torres-Garcia
W
et al
.
Inferring tumour purity and stromal and immune cell admixture from expression data
.
Nat Commun
2013
;
4
:
2612
.
17.
Aran
D
,
Sirota
M
,
Butte
AJ
.
Systematic pan-cancer analysis of tumour purity
.
Nat Commun
2015
;
6
:
8971
.
18.
Taylor
AM
,
Shih
J
,
Ha
G
,
Gao
GF
,
Zhang
X
,
Berger
AC
et al
.
Genomic and functional approaches to understanding cancer aneuploidy
.
Cancer Cell
2018
;
33
:
676
89
.
19.
Drake
JW
,
Charlesworth
B
,
Charlesworth
D
,
Crow
JF
.
Rates of spontaneous mutation
.
Genetics
1998
;
148
:
1667
86
.
20.
Loeb
LA
.
Mutator phenotype may be required for multistage carcinogenesis
.
Cancer Res
1991
;
51
:
3075
9
.
21.
Alexandrov
LB
,
Nik-Zainal
S
,
Wedge
DC
,
Aparicio
SAJR
,
Behjati
S
,
Biankin
AV
et al
.
Signatures of mutational processes in human cancer
.
Nature
2013
;
500
:
415
21
.
22.
Muzny
DM
,
Bainbridge
MN
,
Chang
K
,
Dinh
HH
,
Drummond
JA
,
Fowler
G
et al
.
Comprehensive molecular characterization of human colon and rectal cancer
.
Nature
2012
;
487
:
330
7
.
23.
Funkhouser
WK
,
Lubin
IM
,
Monzon
FA
,
Zehnbauer
BA
,
Evans
JP
,
Ogino
S
et al
.
Relevance, pathogenesis, and testing algorithm for mismatch repair-defective colorectal carcinomas a report of the association for molecular pathology
.
J Mol Diagn
2012
;
14
:
91
103
.
24.
Cortes-Ciriano
I
,
Lee
S
,
Park
WY
,
Kim
TM
,
Park
PJ
.
A molecular portrait of microsatellite instability across multiple cancers
.
Nat Commun
2017
;
8
:
15180
.
25.
Boyiadzis
MM
,
Kirkwood
JM
,
Marshall
JL
,
Pritchard
CC
,
Azad
NS
,
Gulley
JL
.
Significance and implications of FDA approval of pembrolizumab for biomarker-defined disease
.
J Immunother Cancer
2018
;
6
:
137
.
26.
Galon
J
,
Costes
A
,
Sanchez-Cabo
F
,
Kirilovsky
A
,
Mlecnik
B
,
Lagorce-Pages
C
et al
.
Type, density, and location of immune cells within human colorectal tumors predict clinical outcome
.
Science
2006
;
313
:
1960
4
.
27.
Llosa
NJ
,
Cruise
M
,
Tam
A
,
Wicks
EC
,
Hechenbleikner
EM
,
Taube
JM
et al
.
The vigorous immune microenvironment of microsatellite instable colon cancer is balanced by multiple counter-inhibitory checkpoints
.
Cancer Discov
2015
;
5
:
43
51
.
28.
Le
DT
,
Uram
JN
,
Wang
H
,
Bartlett
BR
,
Kemberling
H
,
Eyring
AD
et al
.
PD-1 blockade in tumors with mismatch-repair deficiency
.
New Engl J Med
2015
;
372
:
2509
20
.
29.
Gurjao
C
,
Tsukrov
D
,
Imakaev
M
,
Luquette
LJ
,
Mirny
LA
.
Limited evidence of tumour mutational burden as a biomarker of response to immunotherapy
.
bioRxiv
2020
:
2020.09.03.260265
.
30.
Yarchoan
M
,
Hopkins
A
,
Jaffee
EM
.
Tumor mutational burden and response rate to PD-1 inhibition
.
New Engl J Med
2017
;
377
:
2500
1
.
31.
Robbins
PF
,
Lu
YC
,
El-Gamil
M
,
Li
YF
,
Gross
C
,
Gartner
J
et al
.
Mining exomic sequencing data to identify mutated antigens recognized by adoptively transferred tumor-reactive T cells
.
Nat Med
2013
;
19
:
747
52
.
32.
Bassani-Sternberg
M
,
Braunlein
E
,
Klar
R
,
Engleitner
T
,
Sinitcyn
P
,
Audehm
S
et al
.
Direct identification of clinically relevant neoepitopes presented on native human melanoma tissue by mass spectrometry
.
Nat Commun
2016
;
7
:
13404
.
33.
Roh
W
,
Chen
PL
,
Reuben
A
,
Spencer
CN
,
Prieto
PA
,
Miller
JP
et al
.
Integrated molecular analysis of tumor biopsies on sequential CTLA-4 and PD-1 blockade reveals markers of response and resistance
.
Sci Transl Med
2017
;
9
:
eaah3560
.
34.
Kalaora
S
,
Wolf
Y
,
Feferman
T
,
Barnea
E
,
Greenstein
E
,
Reshef
D
et al
.
Combined analysis of antigen presentation and T-cell recognition reveals restricted immune responses in melanoma
.
Cancer Discov
2018
;
8
:
1366
75
.
35.
Bonneville
R
,
Krook
MA
,
Kautto
EA
,
Miya
J
,
Wing
MR
,
Chen
HZ
et al
.
Landscape of microsatellite instability across 39 cancer types
.
JCO Precis Oncol
2017
;
1
:
PO.17.00073
.
36.
Lindeboom
RGH
,
Vermeulen
M
,
Lehner
B
,
Supek
F
.
The impact of nonsense-mediated mRNA decay on genetic disease, gene editing and cancer immunotherapy
.
Nat Genet
2019
;
51
:
1645
51
.
37.
Litchfield
K
,
Reading
JL
,
Lim
EL
,
Xu
H
,
Liu
P
,
Al-Bakir
M
et al
.
Escape from nonsense-mediated decay associates with anti-tumor immunogenicity
.
Nat Commun
2020
;
11
:
3800
.
38.
McGrail
DJ
,
Garnett
J
,
Yin
J
,
Dai
H
,
Shih
DJH
,
Lam
TNA
et al
.
Proteome instability is a therapeutic vulnerability in mismatch repair-deficient cancer
.
Cancer Cell
2020
;
37
:
371
86
.
39.
Rodvold
JJ
,
Chiu
KT
,
Hiramatsu
N
,
Nussbacher
JK
,
Galimberti
V
,
Mahadevan
NR
et al
.
Intercellular transmission of the unfolded protein response promotes survival and drug resistance in cancer cells
.
Sci Signal
2017
;
10
:
eaah7177
.
40.
Pelka
K
,
Hofree
M
,
Chen
JH
,
Sarkizova
S
,
Pirl
JD
,
Jorgji
V
et al
.
Spatially organized multicellular immune hubs in human colorectal cancer
.
Cell
2021
;
184
:
4734
52
.
41.
Thommen
DS
,
Koelzer
VH
,
Herzig
P
,
Roller
A
,
Trefny
M
,
Dimeloe
S
et al
.
A transcriptionally and functionally distinct PD-1+ CD8+ T cell pool with predictive potential in non-small-cell lung cancer treated with PD-1 blockade
.
Nat Med
2018
;
24
:
994
1004
.
42.
Gros
A
,
Robbins
PF
,
Yao
X
,
Li
YF
,
Turcotte
S
,
Tran
E
et al
.
PD-1 identifies the patient-specific CD8(+) tumor-reactive repertoire infiltrating human tumors
.
J Clin Invest
2014
;
124
:
2246
59
.
43.
Simoni
Y
,
Becht
E
,
Fehlings
M
,
Loh
CY
,
Koo
SL
,
Teng
KWW
et al
.
Bystander CD8+ T cells are abundant and phenotypically distinct in human tumour infiltrates
.
Nature
2018
;
557
:
575
9
.
44.
Chow
MT
,
Ozga
AJ
,
Servis
RL
,
Frederick
DT
,
Lo
JA
,
Fisher
DE
et al
.
Intratumoral activity of the CXCR3 chemokine system is required for the efficacy of anti–PD-1 therapy
.
Immunity
2019
;
50
:
1498
512
.
45.
Marabelle
A
,
Fakih
M
,
Lopez
J
,
Shah
M
,
Shapira-Frommer
R
,
Nakagawa
K
et al
.
Association of tumour mutational burden with outcomes in patients with advanced solid tumours treated with pembrolizumab: prospective biomarker analysis of the multicohort, open-label, phase 2 KEYNOTE-158 study
.
Lancet Oncol
2020
;
21
:
1353
65
.
46.
Subbiah
V
,
Solit
DB
,
Chan
TA
,
Kurzrock
R
.
The FDA approval of pembrolizumab for adult and pediatric patients with tumor mutational burden (TMB) >= 10: a decision centered on empowering patients and their physicians
.
Ann Oncol
2020
;
31
:
1115
8
.
47.
Strickler
JH
,
Hanks
BA
,
Khasraw
M
.
Tumor mutational burden as a predictor of immunotherapy response: is more always better?
Clin Cancer Res
2021
;
27
:
1236
41
.
48.
Cancer Genome Atlas Research Network
,
Weinstein
JN
,
Collisson
EA
,
Mills
GB
,
Shaw
KR
,
Ozenberger
BA
et al
.
The Cancer Genome Atlas pan-cancer analysis project
.
Nat Genet
2013
;
45
:
1113
20
.
49.
Riaz
N
,
Havel
JJ
,
Makarov
V
,
Desrichard
A
,
Urba
WJ
,
Sims
JS
et al
.
Tumor and microenvironment evolution during immunotherapy with Nivolumab
.
Cell
2017
;
171
:
934
49
.
50.
Gide
TN
,
Quek
C
,
Menzies
AM
,
Tasker
AT
,
Shang
P
,
Holst
J
et al
.
Distinct immune cell populations define response to anti–PD-1 monotherapy and anti–PD-1/Anti-CTLA-4 combined therapy
.
Cancer Cell
2019
;
35
:
238
55
.
51.
Liu
D
,
Schilling
B
,
Liu
D
,
Sucker
A
,
Livingstone
E
,
Jerby-Arnon
L
et al
.
Integrative molecular and clinical modeling of clinical outcomes to PD1 blockade in patients with metastatic melanoma
.
Nat Med
2019
;
25
:
1916
27
.
52.
Auslander
N
,
Zhang
G
,
Lee
JS
,
Frederick
DT
,
Miao
B
,
Moll
T
et al
.
Robust prediction of response to immune checkpoint blockade therapy in metastatic melanoma
.
Nat Med
2018
;
24
:
1545
9
.
53.
Ayers
M
,
Lunceford
J
,
Nebozhyn
M
,
Murphy
E
,
Loboda
A
,
Kaufman
DR
et al
.
IFN-gamma-related mRNA profile predicts clinical response to PD-1 blockade
.
J Clin Invest
2017
;
127
:
2930
40
.
54.
Cristescu
R
,
Mogg
R
,
Ayers
M
,
Albright
A
,
Murphy
E
,
Yearley
J
et al
.
Pan-tumor genomic biomarkers for PD-1 checkpoint blockade-based immunotherapy
.
Science
2018
;
362
:
eaar3593
.
55.
Davoli
T
,
Uno
H
,
Wooten
EC
,
Elledge
SJ
.
Tumor aneuploidy correlates with markers of immune evasion and with reduced response to immunotherapy
.
Science
2017
;
355
:
eaaf8399
.
56.
Steiniche
T
,
Rha
SY
,
Chung
HC
,
Georgsen
JB
,
Ladekarl
M
,
Nordsmark
M
et al
.
T-cell–inflamed gene expression profile (GEP) and PD-L1 expression in patients (pts) with esophageal cancer (EC)
.
J Clin Oncol
2019
;
37
:
26
.
57.
Fehrenbacher
L
,
Spira
A
,
Ballinger
M
,
Kowanetz
M
,
Vansteenkiste
J
,
Mazieres
J
et al
.
Atezolizumab versus docetaxel for patients with previously treated non-small-cell lung cancer (POPLAR): a multicentre, open-label, phase 2 randomised controlled trial
.
Lancet
2016
;
387
:
1837
46
.
58.
Ock
CY
,
Hwang
JE
,
Keam
B
,
Kim
SB
,
Shim
JJ
,
Jang
HJ
et al
.
Genomic landscape associated with potential response to anti-CTLA-4 treatment in cancers
.
Nat Commun
2017
;
8
:
1050
.
59.
Jiang
P
,
Gu
S
,
Pan
D
,
Fu
J
,
Sahu
A
,
Hu
X
et al
.
Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response
.
Nat Med
2018
;
24
:
1550
8
.
60.
Perez-Guijarro
E
,
Yang
HH
,
Araya
RE
,
El Meskini
R
,
Michael
HT
,
Vodnala
SK
et al
.
Multimodel preclinical platform predicts clinical response of melanoma to immunotherapy
.
Nat Med
2020
;
26
:
781
91
.
61.
Billadeau
DD
,
Leibson
PJ
.
ITAMs versus ITIMs: striking a balance during cell regulation
.
J Clin Invest
2002
;
109
:
161
8
.
62.
Staub
E
,
Rosenthal
A
,
Hinzmann
B
.
Systematic identification of immunoreceptor tyrosine-based inhibitory motifs in the human proteome
.
Cell Signal
2004
;
16
:
435
56
.
63.
Varin
MM
,
Le Pottier
L
,
Youinou
P
,
Saulep
D
,
Mackay
F
,
Pers
JO
.
B-cell tolerance breakdown in Sjogren's syndrome: focus on BAFF
.
Autoimmun Rev
2010
;
9
:
604
8
.
64.
Murphy
K
,
Weaver
C
.
Janeway's immunobiology
.
New York, NY
:
Garland Science
;
2016
.
65.
Chen
L
,
Flies
DB
.
Molecular mechanisms of T cell co-stimulation and co-inhibition
.
Nat Rev Immunol
2013
;
13
:
227
42
.
66.
Campbell
KS
,
Purdy
AK
.
Structure/function of human killer cell immunoglobulin-like receptors: lessons from polymorphisms, evolution, crystal structures and mutations
.
Immunology
2011
;
132
:
315
25
.
67.
Pende
D
,
Falco
M
,
Vitale
M
,
Cantoni
C
,
Vitale
C
,
Munari
E
et al
.
Killer Ig-like receptors (KIRs): their role in NK cell modulation and developments leading to their clinical exploitation
.
Front Immunol
2019
;
10
:
1179
.
68.
Ward-Kavanagh
LK
,
Lin
WW
,
Sedy
JR
,
Ware
CF
.
The TNF receptor superfamily in co-stimulating and co-inhibitory responses
.
Immunity
2016
;
44
:
1005
19
.
69.
Gonçalves
CM
,
Henriques
SN
,
Santos
RF
,
Carmo
AM
.
CD6, a rheostat-type signalosome that tunes T cell activation
.
Front Immunol
2018
;
9
:
2994
.
70.
Steri
M
,
Orru
V
,
Idda
ML
,
Pitzalis
M
,
Pala
M
,
Zara
I
et al
.
Overexpression of the cytokine BAFF and autoimmunity risk
.
N Engl J Med
2017
;
376
:
1615
26
.
71.
Chen
M
,
Lin
X
,
Liu
Y
,
Li
Q
,
Deng
Y
,
Liu
Z
et al
.
The function of BAFF on T helper cells in autoimmunity
.
Cytokine Growth Factor Rev
2014
;
25
:
301
5
.
72.
Chen
S
,
Lee
LF
,
Fisher
TS
,
Jessen
B
,
Elliott
M
,
Evering
W
et al
.
Combination of 4-1BB agonist and PD-1 antagonist promotes antitumor effector/memory CD8 T cells in a poorly immunogenic tumor model
.
Cancer Immunol Res
2015
;
3
:
149
60
.
73.
Ma
HS
,
Poudel
B
,
Torres
ER
,
Sidhom
JW
,
Robinson
TM
,
Christmas
B
et al
.
A CD40 agonist and PD-1 antagonist antibody reprogram the microenvironment of nonimmunogenic tumors to allow T-cell-mediated anticancer activity
.
Cancer Immunol Res
2019
;
7
:
428
42
.
74.
Barber
A
.
Costimulation of effector CD8+ T cells: which receptor is optimal for immunotherapy?
MOJ Immunol
2014
;
1
:
00011
.
75.
Melero
I
,
Hirschhorn-Cymerman
D
,
Morales-Kastresana
A
,
Sanmamed
MF
,
Wolchok
JD
.
Agonist antibodies to TNFR molecules that costimulate T and NK cells
.
Clin Cancer Res
2013
;
19
:
1044
53
.
76.
Chester
C
,
Sanmamed
MF
,
Wang
J
,
Melero
I
.
Immunotherapy targeting 4-1BB: mechanistic rationale, clinical results, and future strategies
.
Blood
2018
;
131
:
49
57
.
77.
Philipson
BI
,
O'Connor
RS
,
May
MJ
,
June
CH
,
Albelda
SM
,
Milone
MC
.
4-1BB costimulation promotes CAR T cell survival through noncanonical NF-kappaB signaling
.
Sci Signal
2020
;
13
:
eaay8248
.
78.
Siu
LL
,
Steeghs
N
,
Meniawy
T
,
Joerger
M
,
Spratlin
JL
,
Rottey
S
et al
.
Preliminary results of a phase I/IIa study of BMS-986156 (glucocorticoid-induced tumor necrosis factor receptor–related gene [GITR] agonist), alone and in combination with nivolumab in pts with advanced solid tumors
.
J Clin Oncol
2017
;
35
:
104
.
79.
Tolcher
AW
,
Sznol
M
,
Hu-Lieskovan
S
,
Papadopoulos
KP
,
Patnaik
A
,
Rasco
DW
et al
.
Phase Ib study of utomilumab (PF-05082566), a 4-1BB/CD137 agonist, in combination with Pembrolizumab (MK-3475) in patients with advanced solid tumors
.
Clin Cancer Res
2017
;
23
:
5349
57
.
80.
Cohen
EEW
,
Pishvaian
MJ
,
Shepard
DR
,
Wang
D
,
Weiss
J
,
Johnson
ML
et al
.
A phase Ib study of utomilumab (PF-05082566) in combination with mogamulizumab in patients with advanced solid tumors
.
J Immunother Cancer
2019
;
7
:
342
.
81.
Choi
Y
,
Shi
Y
,
Haymaker
CL
,
Naing
A
,
Ciliberto
G
,
Hajjar
J
.
T-cell agonists in cancer immunotherapy
.
J Immunother Cancer
2020
;
8
:
e000966
.
82.
Wu
DC
,
Yao
J
,
Ho
KS
,
Lambowitz
AM
,
Wilke
CO
.
Limitations of alignment-free tools in total RNA-seq quantification
.
BMC Genomics
2018
;
19
:
510
.
83.
Sahraeian
SME
,
Mohiyuddin
M
,
Sebra
R
,
Tilgner
H
,
Afshar
PT
,
Au
KF
et al
.
Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis
.
Nat Commun
2017
;
8
:
59
.
84.
Everaert
C
,
Luypaert
M
,
Maag
JLV
,
Cheng
QX
,
Dinger
ME
,
Hellemans
J
et al
.
Benchmarking of RNA-sequencing analysis workflows using whole-transcriptome RT-qPCR expression data
.
Sci Rep
2017
;
7
:
1559
.
85.
Teng
M
,
Love
MI
,
Davis
CA
,
Djebali
S
,
Dobin
A
,
Graveley
BR
et al
.
A benchmark for RNA-seq quantification pipelines
.
Genome Biol
2016
;
17
:
74
.
86.
Robert
C
,
Watson
M
.
Errors in RNA-seq quantification affect genes of relevance to human disease
.
Genome Biol
2015
;
16
:
177
.
87.
Chu
T
,
Danko
C
.
Bayesian inference of cell composition and gene expression reveals tumor-microenvironment interactions
.
bioRxiv
2020
.
88.
Goldman
MJ
,
Craft
B
,
Hastie
M
,
Repecka
K
,
McDade
F
,
Kamath
A
et al
.
Visualizing and interpreting cancer genomics data via the Xena platform
.
Nat Biotechnol
2020
;
38
:
675
8
.
89.
Bray
NL
,
Pimentel
H
,
Melsted
P
,
Pachter
L
.
Near-optimal probabilistic RNA-seq quantification
.
Nat Biotechnol
2016
;
34
:
525
7
.
90.
Harrow
J
,
Frankish
A
,
Gonzalez
JM
,
Tapanari
E
,
Diekhans
M
,
Kokocinski
F
et al
.
GENCODE: the reference human genome annotation for The ENCODE Project
.
Genome Res
2012
;
22
:
1760
74
.
91.
Robinson
MD
,
McCarthy
DJ
,
Smyth
GK
.
edgeR: a Bioconductor package for differential expression analysis of digital gene expression data
.
Bioinformatics
2010
;
26
:
139
40
.
92.
Robinson
MD
,
Oshlack
A
.
A scaling normalization method for differential expression analysis of RNA-seq data
.
Genome Biol
2010
;
11
:
R25
.
93.
Puram
SV
,
Tirosh
I
,
Parikh
AS
,
Patel
AP
,
Yizhak
K
,
Gillespie
S
et al
.
Single-cell transcriptomic analysis of primary and metastatic tumor ecosystems in head and neck Cancer
.
Cell
2017
;
171
:
1611
24
.
94.
Peng
J
,
Sun
BF
,
Chen
CY
,
Zhou
JY
,
Chen
YS
,
Chen
H
et al
.
Single-cell RNA-seq highlights intra-tumoral heterogeneity and malignant progression in pancreatic ductal adenocarcinoma
.
Cell Res
2019
;
29
:
725
38
.
95.
Ma
L
,
Hernandez
MO
,
Zhao
Y
,
Mehta
M
,
Tran
B
,
Kelly
M
et al
.
Tumor cell biodiversity drives microenvironmental reprogramming in liver cancer
.
Cancer Cell
2019
;
36
:
418
30
.
96.
Li
H
,
Courtois
ET
,
Sengupta
D
,
Tan
Y
,
Chen
KH
,
Goh
JJL
et al
.
Reference component analysis of single-cell transcriptomes elucidates cellular heterogeneity in human colorectal tumors
.
Nat Genet
2017
;
49
:
708
18
.
97.
Lambrechts
D
,
Wauters
E
,
Boeckx
B
,
Aibar
S
,
Nittner
D
,
Burton
O
et al
.
Phenotype molding of stromal cells in the lung tumor microenvironment
.
Nat Med
2018
;
24
:
1277
89
.
98.
Karaayvaz
M
,
Cristea
S
,
Gillespie
SM
,
Patel
AP
,
Mylvaganam
R
,
Luo
CC
et al
.
Unravelling subclonal heterogeneity and aggressive disease states in TNBC through single-cell RNA-seq
.
Nat Commun
2018
;
9
:
3588
.
99.
Bird
A
.
DNA methylation patterns and epigenetic memory
.
Genes Dev
2002
;
16
:
6
21
.
100.
Mitchell
M
.
Genetic algorithms: an overview
.
Citeseer
1995
;
31
39
.

Supplementary data