Data-driven immune checkpoint blocking therapy response prediction

Abstract

A clinical decision support system and method for predicting a clinical response of a target tumor to an immune checkpoint blocking therapy (ICB) by: obtaining a training dataset, the training dataset comprising training transcriptome records of tumor tissue samples of a plurality of known responders and a plurality of known non-responders to the ICB; deconvolution is carried out on the training data set to identify differential expression genes DEG in the training data set, the DEG is regarded as features of the training data set, and related responder and non-responder states are regarded as labels recorded in the training data set; performing regression analysis on the features to select a feature subset which can predict the responder tag or the non-responder tag and related feature weights; incorporating the feature subset and the associated feature weight into a reaction estimation model; receiving transcriptome data of the target tumor tissue sample; the transcriptome data is processed using the generated response estimation model to generate an estimated response indicator indicative of a possible response of the target tumor to the ICB.

Description

TECHNICAL FIELD

[0001] The present disclosure generally relates to methods and systems for data-driven prediction of response to Immune Checkpoint Blockade (ICB) therapy. BACKGROUND

[0002] The purpose of the background description is to present a framework for understanding the present disclosure. The background description is not admitted to be prior art to the present disclosure.

[0003] Breakthroughs in cancer immunotherapy have revolutionized cancer treatment. Immune checkpoint blockade therapy (ICB) is a promising class of immunotherapy that can induce significant tumor shrinkage and long-term disease control in a subset of cancer patients, enabling them to maintain quality of life and continue to contribute to society.

[0004] Although ICB can produce durable clinical responses for many previously untreatable cancer types, ICB is an expensive therapy and not all patients benefit from ICB.

[0005] Biomarker-guided patient selection can improve the response rate and cost-effectiveness of ICB. However, the existing biomarkers are not sufficiently accurate in prediction to stratify treatment for many cancer types. In addition, the application of existing biomarkers in the clinic is constrained by logistical and technical limitations. Therefore, it remains challenging for clinicians to identify which patients can truly benefit from ICB, highlighting the urgent need to develop new methods to help clinicians stratify cancer patients and prioritize for ICB.

[0006] Known patient selection methods are limited and can rely on a small set of biomarkers, such as PD-L1 expression, microsatellite instability, and tumor mutation burden (TMB) greater than 10 mutations / Mb. Known methods have weak predictive power for treatment response and are often not generalizable to cover a wider range of cancer types. Clinical trials to assess the effectiveness of known methods often produce conflicting results on uncertain outcomes to support effective application in a wider set of cancer types in clinical practice. Therefore, it remains challenging for clinicians to identify which patients are more likely to benefit from ICB.

[0007] It would be desirable to overcome or ameliorate at least one of the above problems, or to at least provide a useful alternative. SUMMARY

[0008] Some embodiments relate to a clinical decision support system for predicting clinical response of a target tumor to an immune checkpoint blockade therapy (ICB), the system comprising one or more processing units configured to:

[0009] generate a response estimation model by:

[0010] obtain a training dataset comprising training transcriptomic data (may include training transcriptomic sequence records) of tumor tissue samples of a plurality of known responders and a plurality of known non-responders to an ICB;

[0011] deconvolute the training dataset to identify differentially expressed genes (DEGs) in the training dataset, wherein the DEGs are treated as features of the training dataset and the associated responder, non-responder status are treated as labels recorded in the training dataset;

[0012] perform regression analysis on the features to select a subset of features that best predict the responder or non-responder labels and obtain associated feature weights;

[0013] incorporate the subset of features and the associated feature weights into the response estimation model;

[0014] receive transcriptomic sequence data of a target tumor tissue sample;

[0015] process the received transcriptomic sequence data using the generated response estimation model to generate an estimated response indicator indicative of a likely response of the target tumor to the ICB.

[0016] Some embodiments relate to a system for generating a response estimation model to predict clinical response of a tumor to an immune checkpoint blockade therapy (ICB), the system comprising one or more processing units configured to:

[0017] obtain a training dataset comprising training transcriptomic data (may include training transcriptomic records) of tumor tissue samples of a plurality of known responders and a plurality of known non-responders to an ICB;

[0018] deconvolute the training dataset to identify differentially expressed genes (DEGs) in the training dataset, wherein the DEGs are treated as features of the training dataset and the associated responder, non-responder status are treated as labels recorded in the training dataset;

[0019] perform regression analysis on the features to select a subset of features that best predict the responder or non-responder labels;

[0020] incorporate the subset of features into the response estimation model.

[0021] Some embodiments relate to a computer-implemented method for predicting clinical response of a tumor to an immune checkpoint blockade therapy (ICB), the method comprising:

[0022] generating a response estimation model by:

[0023] obtaining a training dataset comprising training transcriptomic data (may include training transcriptomic records) of tumor tissue samples of a plurality of known responders and a plurality of known non-responders to an ICB;

[0024] deconvoluting the training dataset to identify differentially expressed genes (DEGs) in the training dataset, wherein the DEGs are treated as features of the training dataset and the associated responder, non-responder status are treated as labels recorded in the training dataset;

[0025] performing a regression analysis on the features to select a subset of features that best predict the responder or non-responder labels and obtain associated feature weights;

[0026] incorporating the subset of features and the associated feature weights into the response estimation model;

[0027] obtaining transcriptomic data of a target tumor tissue sample;

[0028] processing the transcriptomic data using the generated response estimation model to generate an estimated response indicator indicative of a possible response of the target tumor to the ICB.

[0029] Some embodiments relate to a computer-implemented method for generating a response estimation model to predict clinical response of a tumor to an immune checkpoint blockade therapy (ICB), the method comprising:

[0030] obtaining a training dataset comprising training transcriptomic data (may include training transcriptomic records) of tumor tissue samples of a plurality of known responders and a plurality of known non-responders to an ICB;

[0031] deconvoluting the training dataset to identify differentially expressed genes (DEGs) in the training dataset, wherein the DEGs are treated as features of the training dataset and the associated responder, non-responder status are treated as labels recorded in the training dataset;

[0032] performing a regression analysis on the features to select a subset of features that best predict the responder or non-responder labels;

[0033] incorporating the subset of features into the response estimation model. BRIEF DESCRIPTION OF DRAWINGS

[0034] Some embodiments of systems and methods for predicting clinical response of a target tumor to immune checkpoint blockade therapy (ICB) according to the present disclosure will now be described, by way of non-limiting example only, with reference to the accompanying drawings in which:

[0035] FIG. 1 A computer-implemented method for predicting clinical response of a tumor to ICB is shown.

[0036] FIG. 2 A data summary and schematic of the workflow for candidate biomarker discovery is shown.

[0037] Figure 3 shows a schematic of differentially expressed genes and pathways between responders and non-responders in the stroma. FIG. 3A A heatmap of genes differentially expressed between all discovery cohorts in the stroma is shown. Heatmaps are colored by log2 fold change in gene expression. P-values from each cohort are combined using the Fisher’s method and presented as -log10(meta p-value) in a bar plot. Genes with a q-value less than 0.01 and a p-value less than 0.1 in at least one ICB cohort are shown. FIG. 3B A bar plot of deconvoluted stroma expression of top differentially expressed genes is shown. Error bars represent estimated standard error greater than 1. FIG. 3C Pathway enrichment of differentially expressed genes in the stroma is shown. Bar plots show -log10(p-value) of significantly enriched pathways. Similar pathways are grouped by color. FIG. 3D A network representation of enriched pathways is shown. Large nodes represent enriched pathways. Similar pathways are merged and shown in the same color. Small nodes represent genes associated with each pathway.

[0038] Figure 4 shows differentially expressed genes between responders and non-responders in the cancer compartment. FIG. 4A A heatmap of genes differentially expressed between all discovery cohorts in the cancer cells is shown. Heatmaps are colored by log2 fold change in gene expression. P-values from each cohort are combined using the Fisher’s method and presented as -log10(meta p-value) in a bar plot. Genes with a q-value less than 0.01 and a p-value less than 0.1 in at least one ICB cohort are shown. FIG. 4B A bar plot showing the number of important genes found in the stroma and cancer compartment at different q-value cutoffs is shown. In FIG. 4C Pathway analysis of genes with cancer-specific expression differences is shown in FIG. 4E Urothelial carcinoma. FIG. 4C Pathway analysis of genes with cancer-specific expression differences is shown in FIG. 4E Urothelial carcinoma. A bar plot of deconvoluted cancer expression of differentially expressed genes in top enriched pathways is shown. Error bars represent estimated standard error.

[0039] Figure 5 shows analysis of ligand and receptor expression of immunomodulators. FIG. 5A and FIG. 5B Heatmaps showing inferred ligand and receptor expression in cancer and stromal compartments for immune checkpoint proteins and cytokines. Heatmaps are colored by differentially expressed signed -log10(p-value).

[0040] Figure 6 shows a predictive model of response to immune checkpoint inhibition. FIG. 6A shows a schematic of the data and methods used in model training and testing. FIG. 6B shows top ranked features of ICB response selected by LASSO logistic regression. FIG. 6C shows performance comparison of the IOselect model to a 3-feature model with existing biomarkers (TMB + IFNG + FCER1A) in 6 independent test cohorts. Paired t-test was used to calculate the difference between the mean AUC of IOselect and other models. ** p-value < 0.01, * p-value < 0.05,. p-value < 0.1. FIG. 6D shows Kaplan-Meier curves of progression-free survival for patients with high and low (greater or less than median) predicted scores; p-values from log-rank test are shown.

[0041] FIG. 7 shows a schematic summarizing key components associated with ICB response. Red arrows indicate genes positively associated with ICB response. Green arrows indicate genes negatively associated with ICB response.

[0042] FIG. 8 shows a schematic of bulk tumor expression deconvolution and permutation test to identify differentially expressed genes between responders (Microsatellite Instable, MSI) and non-responders (Microsatellite Stable, MSS) in cancer and stromal compartments.

[0043] Figure 9 shows identification of genes consistently differentially expressed between responders and non-responders in the discovery cohort. FIG. 9A and FIG. 9BVolcano plots of differentially expressed genes for stromal and cancer compartments are shown in FIGS. 1A and IB, respectively. The y-axis is the -loglO (meta-p value) of differential expression. The meta-p value was calculated using the Fisher method, combining p-values across all cohorts. The x-axis is the mean log2 (fold change) in gene expression for responders (MSI / Epstein-Barr Virus (EBV)) versus non-responders (MSS). The red line represents the cutoff of q-value less than 0.01. Genes with median absolute log2 (fold change) greater than 0.5, q-value less than 0.01, and p-value less than 0.1 in at least one ICB cohort are colored and labeled. FIG. 9C and FIG. 9D Heat maps of differentially expressed genes for stromal and cancer compartments are shown in FIGS. 2A and 2B, respectively. The heat maps are colored by the signed -loglO (p-value) of differential expression. Genes with median absolute log2 (fold change) greater than 0.5, q-value less than 0.01, and p-value less than 0.1 in at least one ICB cohort are shown.

[0044] FIG. 10 Validation of bulk tumor expression deconvolution is shown in FIG. 3. Expression of known series-specific genes in cancer and stromal compartments was inferred in the ICB cohort.

[0045] FIG. 11 Differential expression of PD1 receptor and its ligands PD-L1 and PD-L2 in stroma is shown in FIG. 4. Bar plot of deconvoluted stromal expression of PD1 (PDCD1) receptor and its ligands PD-L1 (CD274) and PD-L2 (PDCD1LG2). Error bars represent estimated standard error.

[0046] FIG. 12 Differentially enriched immune cell types between responders and non-responders are shown in FIG. 5. Box plots show Cibersort x scores for 4 cell types that are significantly different between responders (MSI) and non-responders (MSS).

[0047] FIG. 13 Correlation between abundance of immune cell subtypes and gene expression of top stromal biomarkers in individual cohorts is shown in FIG. 6.

[0048] FIG. 14 Expression of differentially expressed genes in immune cell subtypes is shown in FIG. 7. Normalized expression of 4 representative stromal biomarkers in different immune cell types. Immune cell expression data derived from single-cell and flow cytometry sorted RNA-seq experiments of Human Protein Atlas.

[0049] FIG. 15 A schematic of IOSelect model training and testing is shown in FIG. 8.

[0050] FIG. 16 Benchmark data for the IOselect model is shown.

[0051] FIG. 17 This is a schematic diagram of the integration of the IOselect model and clinical workflow.

[0052] FIG. 18 Schematic diagram showing the verification and improvement of the IOselect model.

[0053] FIG. 19 A schematic diagram showing the ICB reaction estimation method.

[0054] FIG. 20 Shown are experimental results of tumor transcriptome deconvolution.

[0055] FIG. 21 A schematic diagram illustrating tumor purity estimation according to some embodiments is shown.

[0056] FIG. 22 Figure 2 shows an assessment of cell type-specific expression of top DEGs using single-cell RNAseq data. A bar graph showing the percentage of cells expressing each gene in individual cell types is shown using scRNAseq data from renal cell carcinoma patients treated with ICB (Bi et al., Cancer Cell, 2021).

[0057] FIG. 23 The area under the receiver operating curve (AUC) of the cancer-DEG model on the six test datasets is shown. The black dots show the mean AUC, and the error bars represent the standard error. Using a one-sample t-test, we found that the predictive power of the model was not significantly different from the random expectation value of AUC = 0.5 (p value = 0.34).

[0058] FIG. 24 Differential expression of the PD1 receptor and its ligands PD-L1 and PD-L2 in the stroma is shown. Bar graph of deconvoluted stroma expression of the PD1 (PDCD1) receptor and its ligands PD-L1 (CD274) and PD-L2 (PDCD1LG2). Error bars represent estimated standard errors.

[0059] FIG. 25 The association between model prediction and overall survival is shown. Kaplan-Meier curves for overall disease-free survival are shown for patients with high versus low prediction scores (greater or less than the median); p-values ​​from the log-rank test are shown. DETAILED DESCRIPTION

[0060] The following will refer to FIGS. 1-25Embodiments of computer-implemented systems and methods for predicting clinical response of a tumor to immune checkpoint blockade therapy (ICB) are described. It should be understood that the restriction of the description to preferred embodiments of the invention is merely for the convenience of discussing the invention and that it is contemplated not to depart from the scope of the appended claims.

[0061] The disclosed embodiments relate to a clinical decision support system and computer- implemented method for predicting the clinical response of a target tumor to ICB. These embodiments include a response estimation model (also referred to as the IOselect model / IOSelect clinical platform) for processing sequencing transcriptomic data from cancer tissue samples and generating predictions about the potential impact / benefit of ICB. The data-driven clinical platform IOselect prioritizes cancer patients for immunotherapy based on the molecular profile of their individual tumors. IOSelect was developed based on gene signatures identified from analysis of transcriptomic data from cancer tissues of known responders and non-responders to ICB. IOSelect uses multi-omic data and machine learning to predict the response potential of cancer patients to ICB. IOSelect presents a comprehensive view of the tumor immune context in a user-friendly data dashboard, including personalized scores of predicted response from our multivariate model. Clinicians can make more informed treatment recommendations based on the interpretable machine learning platform (IOSelect) rather than relying on single marker tests that fail to capture the full complexity of the tumor.

[0062] Since transcriptomic differences between MSI and MSS tumors can predict tumor response to ICB, a meta-analysis of ICB response was performed on 1486 tumors. These tumor samples included transcriptomes of 534 pre-treatment tumors from patients who received ICB treatment, including 4 studies of 3 tumor types (gastric cancer, urothelial cancer, melanoma), and 952 tumors of 3 cancer types (colorectal, gastric, and endometrial cancer) with the highest MSI tumor frequency in TCGA, as shown in FIG. 2 Using deconvolution techniques to estimate stroma-specific and cancer-specific gene expression changes, common transcriptomic features of ICB response across cancer types were identified. Machine learning models (IOSelect model / algorithm) were built to predict patient response to ICB based on the identified novel biomarkers.

[0063]

[0064]

[0065] Table 1. Clinical characteristics of study cohorts

[0066] Such machine learning models are embodied by the methods used herein. Such computer- implemented methods 100 are reflected in FIG. 1In particular, the present disclosure relates to methods 100 for predicting clinical response of a tumor to immune checkpoint blockade therapy (ICB). The methods 100 include:

[0067] (a) obtaining a training dataset;

[0068] (b) deconvolving the training dataset to identify differentially expressed genes (DEGs);

[0069] (c) performing regression analysis on the features to select a subset of features that best predict the responder or non-responder label and obtain associated feature weights (step 102c);

[0070] (d) incorporating the subset of features and associated feature weights into a response estimation model (step 102d);

[0071] obtaining transcriptomic data for a target tumor tissue sample (step 104); and

[0072] processing the transcriptomic data using the generated response estimation model to generate an estimated response indicator indicative of a likely response of the target tumor to ICB (step 106).

[0073] Analysis and identification of biomarkers of immune checkpoint blockade therapy (ICB) response

[0074] In most published ICB studies, systematic analysis of potential biomarkers of ICB response is limited to a small number of tumors with molecular profiles. In contrast, through systematic efforts in tumor sequencing over the past few decades, molecular data for a large number of microsatellite instability (MSI) tumors are publicly available. Given the large response difference between microsatellite instability (MSI) and microsatellite stable (MSS) tumors (almost half of MSI tumors respond to immune ICB, while only 5% of MSS tumors respond), the embodiments exploit the transcriptomic differences between MSI and MSS tumors that are likely to be generalizable for predicting tumor response to ICB. The embodiments exploit (1) a large The Cancer Genome Atlas (TCGA) cohort with high MSI prevalence and (2) a comprehensive analysis of tumors from patients treated with ICB to identify features that predict ICB response. In some experiments, transcriptomes of 534 tumors from 4 ICB cohorts from 3 cancer types (gastric cancer, urothelial carcinoma, and melanoma) and transcriptomes of 952 tumors from 3 TCGA cancer types (gastric cancer, colorectal cancer, and endometrial cancer) with the highest MSI tumor frequency are meta-analyzed.

[0075] Tumors are heterogeneous masses of cancer cells and non-cancerous stromal cells, including immune cells, fibroblasts, and endothelial cells. Responses to immunotherapy depend on a multitude of co-stimulatory and co-inhibitory interactions between effector cells, antigen-presenting immune cells (e.g., dendritic cells and macrophages), and cancer cells in the Tumor Micro-Environment (TME). Dysregulation of these ligand-receptor interactions or immune checkpoints can lead to immune suppression in tumors. However, it is a challenge to utilize bulk tumor transcriptomics to distinguish the complex interactions between cancer and immune cells. The heterogeneity of the TME can confound the discovery of molecular signatures. To distinguish signals from cancer and stromal cells, an expression deconvolution method was developed and applied to infer differentially expressed genes (DEGs) from multiple candidate genes. In some embodiments, DEGs are identified separately in cancer and stromal compartments. Experiments show that the most consistent differentially expressed signals come from the stromal rather than the cancer compartment. This suggests that cancer-intrinsic features of ICB responses tend to be cancer-type specific, while stromal-intrinsic features are more applicable to other cancer types, and thus are better biomarkers for ICB response prediction. From this analysis of experimental results, 59 stromal genes were identified that are consistently differentially expressed between responders (MSI) and non-responders (MSS) in each cancer group, as shown in FIG. 3A

[0076] ​To build a predictive model of ICB response, various feature selection strategies and machine learning methods were evaluated, including Least Absolute Shrinkage and Selection Operator (LASSO) regression, stepwise regression, random forest, gradient boosting, and Support Vector Machine (SVM). LASSO regression appeared to be the best performing method among the available data, while other regression methods could also be applicable to different data sets. As described in step 102a, a training data set was generated from the training transcriptomic records of ICB for tumor tissue samples of a plurality of known responders and a plurality of known non-responders. The reaction estimation model or IOSelect model was generated from 534 tumors with clinical annotation of ICB response (responder or non-responder label) (used as training data). A robust selection of the most informative features in ICB response was then identified by experiment using LASSO logistic regression. In some embodiments, the IOSelect model was constructed using 4 selected genes. Other experiments using different training data sets can lead to selection of different sets of genes. Next, experiments were performed to evaluate the model performance on the training cohort (using 10-fold cross-validation) and 2 independent test cohorts. Test cohort 1 consisted of 49 tumors of melanoma patients who received treatment with the PD-L1 inhibitor nivolumab. Test cohort 2 consisted of 87 tumors of 20 solid cancer types who received treatment with mixed ICB drugs. The area under the receiver operating curve (AUC) of the IOselect model was calculated and compared with existing biomarkers of ICB response (i.e., tumor mutation burden (TMB), PD-L1 expression, and T-cell gene expression profile (GEP)), as shown in FIG. 6C and FIG. 16

[0077] ​In most cases, the current response estimation models significantly outperformed existing biomarkers. The 4-gene signature in the IOSelect model is a measure of favorable tumor microenvironment, while TMB is a measure of tumor immunogenicity, where increased mutational burden is associated with higher neoantigen burden and anti-tumor immune response. However, the predictive power of TMB varies across cancer types, making it difficult to define a single TMB threshold across various cancer types. Therefore, they represent different aspects of anti-tumor immune response and thus can be complementary biomarkers. The IOselect multi-omic model of some embodiments combines the 4-gene signature and TMB and shows that the IOselect multi-omic model has the best performance for all cohorts except for cohort 1, where the IOselect RNA-only model performs slightly better than the multi-omic model. Overall, the IOSelect multi-omic model is 15-40% higher than PD-L1 expression, as shown in FIG. 6C and FIG. 16 .

[0078] The clinical software platform / clinical decision support system according to embodiments allows prioritization of treatment for immunotherapy for cancer patients. The platform employs state-of-the-art bioinformatics pipelines that streamline pre-processing, quantification / mutation calling of transcriptomic and genomic data. New patient data containing transcriptomic sequences can be standardized against existing cohorts to remove potential bias from batch effects and technical differences. Then, the IOSelect algorithm outputs the predicted probability of ICB response using the RNA-only model and / or the multi-omic model depending on data availability. To get a more comprehensive view of the patient’s potential response to ICB, IOSelect can calculate the percentile score of published biomarkers (e.g., TMB, PD-L1 expression, and T-cell GEP) for clinical discussion, as shown in FIG. 16 .

[0079] If whole exome sequencing is performed, the response estimation model of some embodiments also demonstrates clinically actionable genomic mutations related to immunotherapy, such as V-Raf Murine Sarcoma Viral Oncogene Homolog B (BRAF) mutations in melanoma and Epidermal Growth Factor Receptor (EGFR) mutations in lung cancer. FIG. 15 Individual molecular features of each tumor are shown, as well as the contribution of each feature to the final score.

[0080] Currently, PD-L1 expression is the most widely adopted biomarker to guide clinical decisions for ICB treatment, but single-marker tests such as PD-L1 immunohistochemistry (IHC) are unlikely to capture the full complexity of the tumor. T-cell GEP is another popular biomarker of ICB response, measuring the extent of T-cell infiltration according to gene expression profiles. The multi-omic approach of the present reaction estimation model combines the identified gene signatures with TMB, which can better capture the multi-faceted tumor immune context, thus outperforming the existing biomarkers. Currently, the clinical application of sequencing-based immunotherapy response assays is limited by the lack of technical solutions that can effectively process, interpret, and present the large data sets to clinicians. The proposed clinical platform aims to address these technical hurdles by simplifying data processing, analysis, and response prediction. The clinical platform presents a comprehensive view of the tumor immune context, including personalized scores of predicted response from multi-variate models as well as other complementary biomarkers. Thus, the reaction estimation model facilitates a fast turnaround time from sample collection to data presentation.

[0081] Tumor transcriptome deconvolution for identifying common features of ICB response in cancer and stromal cells

[0082] FIG. 1 Step 102b includes deconvoluting the training data set to identify differentially expressed genes (DEGs) in the training data set. In doing so, the DEGs are treated as features of the training data set, and the associated responder, non-responder status are treated as labels recorded in the training data set. In some embodiments, transcriptome deconvolution is performed to estimate stromal and cancer cell expression of each gene in tumors from responders (MSI) and non-responders (MSS) separately, as shown in FIG. 2 Since Epstein-Barr virus (EBV)-positive gastric tumors are also enriched in ICB responders, MSI and EBV tumors are grouped together as a potential responder group for gastric cancer. In some experiments, tumor purity is computed for each sample based on transcriptome and / or genomic data, if available. Then, for each response group, a negative least squares regression of gene expression on tumor purity can be performed. This is done to estimate the mean expression of each gene in stromal and cancer cells in step 102b. After quantifying the expression differences of each gene between responders and non-responders in cancer and stromal cells, permutation-based statistics are used to identify differentially expressed genes (DEGs), as shown in FIG. 2 and FIG. 8

[0083] ​In some embodiments, the first step of deconvolution is to determine the proportion of cancer cells in a tumor sample. This can be achieved using a combination of methods that leverage somatic variant allele frequency or B-allele frequency data (obtained from WGS or Exome-seq data) to estimate ploidy and purity of a given tumor sample. Since estimates vary widely between individual algorithms, and algorithms sometimes fail to converge, some embodiments use a combination of missing data imputation and quantile normalization to obtain a robust estimate of tumor purity, as shown in FIG. 21 The second step of some embodiments aims to construct a deconvolution method to separate the transcriptomic profiles of cancer, immune, and other stromal cells. This step aims to robustly estimate the relative proportions of immune and stromal cells in a tumor. Some embodiments achieve this step using transcript signatures to estimate the relative proportions of these two cell type supersets in a tumor. With these three relative cell type proportions estimated from tumor genomic and transcriptomic data, a constrained regression method can be applied to estimate the mean gene expression level of a given gene in each cell type across a cohort of tumors. For a given gene, the constrained regression method can be computed using the following equation:

[0084]

[0085] where p 癌症,i + p 免疫,i + p 基质,i = 1.

[0086] Some embodiments make the simplifying assumption that these (non-negative) mean cell type expression levels are constant across the tumor cohort, and estimate using non-negative least squares regression. Experiments show that the above equation tends to overestimate stromal gene expression for genes where somatic copy number alterations (CNAs) affect gene expression in a subset of samples (e.g. ERBB2 in HER2-positive breast tumors). Some embodiments therefore employ a modified approach for such genes. These embodiments first identify genes with a correlation between CNA and mRNA expression in a given cohort of tumors (compare sample expression with diploid and non-diploid CNAs, Mann-Whitney U test, P < 1e-6 to account for multiple testing), and then estimate gene expression for the cancer and stromal compartments. The estimation is performed using a two-step approach. First, the above method is used to infer stromal compartment mRNA expression, using only samples with diploid copy number for the gene. Then, the inferred mean stromal compartment expression, the measured mean tumor expression, and the mean purity of the tumor samples are used to compute the mean cancer compartment expression using the above equation. Permutation-based statistics are used to identify cancer / stromal cell differentially expressed genes (DEGs) between ICB responders and non-responders, as shown in FIG. 2 and FIG. 8are shown.

[0087] Top signature genes associated with ICB response are enriched in immune pathways

[0088] From this unbiased transcriptomic analysis, 59 consensus DEGs were identified in the stroma associated with ICB response, as shown in Figures FIG. 3A , 3B , 26A and 26C, of which 51 are genes with immune-related functions, including: genes encoding checkpoint proteins (e.g. CD274, LAG3, CTLA4), genes involved in interferon gamma signaling (e.g. IFNG, IRF1, STAT1), genes involved in T cell effector functions (e.g. PRF1, GZMA, GZMB), and genes involved in NK cell-mediated cytotoxicity (KLRD1, KLRC2). The top 2 DEGs overexpressed in ICB responders were chemokines CXCL9 and CXCL13. CXCL9 is a key chemokine responsible for T cell trafficking through binding to the receptor CXCR3 on T cells, which is associated with increased infiltration of T cells to tumors. CXCL13 mediates involvement in B cell and T cell migration through its receptor CXCR5. T follicular helper (TFH) cell production of CXCL13 is critical for the formation of tertiary lymphoid structures at tumor sites and activation of germinal center B cells, which is associated with anti-tumor response of ICT. Several recent studies have suggested that CXCL9 and CXCL13 are predictive factors for improved ICB response. Notably, a recent meta-analysis also highlighted these two chemokines, which evaluated existing biomarkers of ICB response, suggesting that CXCL9 and CXCL13 are indeed the strongest individual transcriptomic predictors of ICB response independent of tumor type.

[0089] Stromal genes negatively associated with ICB response

[0090] In some experiments, a set of immune-related genes that were negatively correlated with ICB outcome were also observed. FCER1A showed the most consistent and largest overexpression in ICB non-responders. FCER1A encodes a subunit of the IgE receptor. While FCER1A is an established marker of basophils, mast cells, and dendritic cells, it was also shown to be expressed on immunosuppressive M2 macrophages and Tumor Associated Macrophages (TAMs). CH25H is a gene encoding cholesterol 25-hydroxylase, which catalyzes the formation of 25-hydroxycholesterol (25-HC), and was also upregulated in ICB non-responders. 25-HC oxysterols are produced by macrophages in response to type I interferon signaling and have diverse immune effects, including B-cell chemotaxis, macrophage differentiation, and modulation of inflammatory responses. Two (2) natriuretic peptides, NPR1 and NPR3, among the genes upregulated in ICB non-responders were also observed. While the primary function of natriuretic peptides is to regulate electrolytes, natriuretic peptides are expressed in immune cells and there is growing evidence that natriuretic peptides play a role in inflammation and immunity.

[0091] Next, the biological pathways enriched in the DEGs associated with ICB response were analyzed. Based on this analysis, 31 enriched gene ontology biological process terms were identified, as shown in Table 3. FIG. 3C As shown, all of these terms are associated with immunity. FIG. 3D These 31 terms can be further grouped into 4 clusters based on similarity, as shown in Table 4.

[0092] In cancer cells, 9 genes differentially expressed between ICB responders and ICB non-responders were identified by the experiments, as shown in Tables 5, FIG. 4A 9B and 9D. However, the differential expression signal was strongest and most consistent among the TCGA cohort and the gastric cancer ICB cohort, where all responders belonged to the MSI subtype. Thus, the DEGs discovered in cancer cells can be related to MSI biology, but are not specific to ICB response. Overall, fewer cancer-specific consensus DEGs were observed compared to stromal-specific consensus DEGs, regardless of the q-value cutoff disclosed in FIG. 5B In addition, no functional enrichment was found even at a lower q-value cutoff of 0.25. The lack of consensus DEGs across cancer types suggests that the ICB response-related intrinsic differences within cancer tend to be cancer type-specific.

[0093] ​Next, pathway enrichment analysis was performed for each cancer type to determine potential cancer-specific processes involved in ICB response. For melanoma, the experiment identified only 22 DEGs common between the two melanoma cohorts, and no enriched pathways were found. In gastric cancer, the experiment found 190 consensus DEGs between the gastric ICB cohort and the TCGA STAD cohort. Interestingly, FIG. 4C Two angiogenesis-related functional clusters, “positive regulation of collagen biosynthetic process” and “regulation of vascular smooth muscle cell proliferation”, were disclosed to be upregulated in ICB non-responders. FIG. 4D The pro-angiogenic genes ENG, F2R, and PDGFRB were shown to be overexpressed in ICB non-responders (MSI / EBV) compared to ICB responders (MSS). This result is in line with the understanding that angiogenesis is associated with immune suppression and suggests that a pro-angiogenic environment is associated with ICB resistance in gastric cancer. In urothelial carcinoma, cancer-specific DEGs were enriched in neuron-related functional terms, as shown in FIG. 4E Neuron genes were overexpressed in ICB responders compared to ICB non-responders, as shown in FIG. 4F This supports previous findings that the neuron TCGA expression subtype of urothelial carcinoma is associated with poorer survival but higher response rate to ICB.

[0094] Single-cell RNA sequencing confirms differential expression of stromal biomarkers in specific immune cell subsets

[0095] A single-cell RNAseq (scRNAseq) dataset of renal cell carcinoma patients treated with ICB was used to investigate the expression of stromal DEGs in specific stromal cell subsets associated with ICB response. Stromal (non-cancerous) cells differentially expressed stromal DEGs between ICB responders and non-responders, in the predicted direction, as shown in FIG. 26A The expression distribution of stromal DEGs in individual cell types was then examined. The study found that CD8+ T cells were enriched in CXCL13 and IFNG, while M1-like tumor-associated macrophages (TAMs) expressed the highest levels of CXCL9 and CD274, as shown in FIG. 22 PID1 and CH25H were enriched in myeloid cells and fibroblasts, while NPR3, NPR1, PREX2, and SDPR were exclusively expressed by endothelial cells. Of the DEGs associated with poor response to ICB, FCER1A was most frequently expressed by dendritic cells, mast cells, and M2-like TAMs, as shown in FIG. 26AThe top DEGs are shown in Figure 2. Thus, the stromal DEGs include several indicators of the tumor microenvironment, such as the abundance of fibroblasts, T cell infiltration, macrophage polarization, and angiogenesis. The cell type-specific expression of the top stromal DEGs was confirmed using orthogonal expression data of immune cell subtypes from flow cytometry and scRNAseq experiments from the Human Protein Atlas, as shown in Figure 3. FIG. 14 In addition, the study found that the increased abundance of specific immune cell subsets contributed to the gene expression signature differences between ICB responders and non-responders. For example, the frequency of CXCL9+Ml-like TAMs and CXCL13+T cells was higher in ICB responders than in non-responders, while the frequency of FCER1A+M2-like TAMs and FCER1A+dendritic cells was higher in non-responders, as shown in Figure 4. These findings imply that some elements of the TME, such as stromal DEGs, are responsible for determining the ICB sensitivity of tumors. FIG. 26B

[0096] Immune cell subtype analysis identifies Ml macrophages as closely associated with ICB response

[0097] The immune cell composition of the TME is a determinant of immunotherapy response. The experiments used Cibersortx to estimate the immune cell composition of each tumor and compared the abundance of immune cell subsets between ICB responders and non-responders. It was found that macrophages Ml, T follicular helper cells (Tfh), and activated CD4 memory T cells were most abundant in ICB responders compared to non-responders. Macrophages Ml are known to have pro-inflammatory and anti-tumor effects. The results further suggest that Ml macrophage infiltration can be a universal feature of ICB response in multiple cancer types. Tfh cells are specialized CD4+T cells that contribute to the formation of germinal centers in Tertiary Lymphoid Structures (TLSs). These cells interact with B cells at the TLS to activate the antibody response of B cells and enhance CD8+T cell effector functions. The presence of Tfh cells has been associated with better survival rates in breast, lung, and colorectal cancers. However, the extent of Tfh cell involvement in ICB response has not been well-established. Our analysis also showed that while resting CD4 memory T cells were enriched in non-responders, activated CD4 memory T cells were enriched in ICB responders. This suggests that pre-existing CD4+T cell immunity is important for effective anti-tumor responses following ICB treatment.

[0098] Next, the top N (N being, for example, 10 or another desired number) DEGs associated with ICB response were analyzed to assess whether these DEGs represent differential infiltration of specific immune cell subsets. The experiments found that the top gene signatures of ICB response were associated with the abundance of Ml macrophages, CD4+and CD8+T lymphocytes, and activated NK cells, as shown in Figure 5. FIG. 26D ​and FIG. 13 In contrast, DEGs associated with poor ICB response, including FCER1A, were associated with the resting state of mast cells, CD4+ memory T cells, and dendritic cells, suggesting that these genes can be markers of a quiescent or immunosuppressive microenvironment, as FIG. 26D and FIG. 13 as shown.

[0099] Intrinsic cancer characteristics of ICB response tend to be tumor type specific

[0100] Previous studies have proposed multiple mechanisms by which cancer cells can directly suppress tumor immunity and hinder ICB therapy. Surprisingly, our meta-analysis across different tumor types showed that there were few DEGs associated with cancer cell ICB response. In cancer cells, we identified 9 genes that were differentially expressed between ICB responders and ICB non-responders in each tumor type, as shown in FIG. 4A and FIG. 26B and 26D However, we observed that the differential expression signal of these 9 genes was mainly driven by the MSI ICB responder enriched cohort. In contrast to the stromal DEGs shown in FIG. 3A the signal of cancer DEGs in different ICB cohorts was less consistent, as shown in FIG. 4A Therefore, it was hypothesized that the cancer DEGs identified by the meta-analysis can be related to the life mechanisms of MSI tumors rather than the general mechanism of ICB response. In fact, when tested on tumor transcriptome data from 6 independent ICB studies in FIG. 23 the 8 cancer DEGs could not predict ICB response. It was also confirmed that the lack of cancer cell DEGs compared to stromal cell DEGs persisted regardless of the significance cutoffs in the meta-analysis, as shown in FIG. 4B and no functional enrichment or relationship between cancer DEGs could be identified.

[0101] Next, the same DEG analysis was performed within individual tumor types to explore the possibility of tissue-dependent processes involved in ICB response. For melanoma, only 22 cancer DEGs were identified that were common between the two melanoma cohorts, and these genes were not enriched in specific pathways. In gastric cancer, 190 cancer ICB response DEGs were identified in both the gastric cancer ICB cohort and the TCGA cohort. Interestingly, genes associated with poor ICB response were enriched for functions related to angiogenesis, as shown in FIG. 4C genes ENG, F2R, and PDGFRB, which are pro-angiogenic, were overexpressed in ICB non-responders (MSI / EBV) compared to ICB responders (MSS), as shown in FIG. 4Das shown. These data are consistent with the understanding that angiogenesis is associated with immune suppression, suggesting that a pro-angiogenic environment hinders the ICB response in gastric cancer. In urothelial carcinoma, we found that cancer DEGs were enriched for neuron-related functional terms, such as FIG. 4E as shown. This observation is consistent with previous studies that associated neuronally expressed subtypes of urothelial carcinoma with poorer overall survival but higher response rates to ICB. FIG. 4F

[0102] Overall, while these results suggest an overall lack of tissue-agnostic, cancer-intrinsic features of ICB response, they also highlight the possibility that tissue- and tumor type-specific expression features in cancer cells can predict ICB response.

[0103] Lack of cancer-intrinsic immune checkpoint regulation across cancer types

[0104] While transcriptome-wide analysis failed to identify conserved cancer-specific DEGs of ICB response, experiments were performed to examine whether crosstalk between immune checkpoint ligands and receptors between cancer and immune cells can contribute to primary ICB resistance. Specifically, the expression patterns of 30 known immune checkpoint ligand-receptor pairs in cancer and stromal components of the TME were examined to identify ligand-receptor pairs with coordinated expression changes between ICB responders and non-responders. In this targeted analysis, all ligand-receptor genes with nominal p-values greater than 0.1 in the differential gene expression permutation test were investigated to find more subtle signals that can have been missed in the transcriptome-wide analysis. The study found that stromal overexpression of 13 checkpoint receptors and 2 checkpoint ligands was associated with ICB response across cancer types, as shown in FIG. 5A as shown. The checkpoint ligands PD-L1 and PD-L2 and their receptor PD-1 were coordinately upregulated in the stroma of ICB responders compared to non-responders. PD-L1 (CD274) and PD-1 were significantly upregulated in the stroma of ICB responders compared to non-responders in 6 out of 7 cohorts studied (p-value < 0.05, permutation test), while PD-L2 was upregulated in 5 out of 7 cohorts as shown in FIG. 11 as shown. Interestingly, no checkpoint ligands were consistently differentially expressed in cancer cells across tumor types and associated with ICB response.

[0105] In addition to immune checkpoints, cytokines also play an important role in tumor immunity. The expression of 13 pairs of cytokines and cytokine receptors known to be involved in tumor immunity was examined with the aim of identifying cytokine signals associated with ICB response, as shown in FIG. 5B ​Similar to the analysis of immune checkpoints, cancer-intrinsic cytokine expression was not associated with ICB response across tumor types. However, six cytokines (IFNG, FASLG, CXCL13, CCL5, CXCL10, and CXCL9) and the corresponding receptors for three cytokines (CCR5 for CCL5, CXCL9, and CXCR3 for CXCL10) were significantly upregulated in the stroma of ICB responders, consistent with previous observations. Overall, our analysis did not find evidence for a universal mechanism by which cancer cells suppress immune cells through checkpoint interactions or chemokine / cytokine signaling.

[0106] Multivariate classifier robustly predicts ICB response across cancer types

[0107] Next, experiments were performed to verify whether the previously identified consensus DEGs could improve the prediction accuracy of ICB response, as FIG. 6A The experiment focuses on FIG. 3A The 59 stromal-specific DEGs are shown because they are more consistent across cancer types than cancer-specific DEGs and are therefore more likely to generalize to other cancer types. To capture the balance of immune activation and inhibitory signals in the TME, the IOselect immune score was developed, which is the ratio of the average expression of the top 10 DEGs associated with ICB response to the average expression of all 8 DEGs associated with ICB resistance.

[0108] While the stromal transcriptome signature is a measure of tumor microenvironment favorability, TMB is a measure of tumor immunogenicity. Because stromal DEGs and TMB represent different aspects of tumor immunity, they may be orthogonal predictors of ICB response. Therefore, a multivariate classifier IOselect, or clinical decision support system, was developed to predict the clinical response of target tumors to immune checkpoint blockade therapy (ICB), using the IOselect immune score and TMB as input features. The IOSelect classifier was trained on 461 samples from the ICB discovery cohort with transcriptome and genomic data. Model performance was evaluated on 6 independent test cohorts receiving ICB treatment, including 232 patients with melanoma, colorectal cancer, lung cancer, and other cancers. The area under the receiver operating curve (AUC) of the model was calculated and compared with established biomarkers of ICB response (i.e., CXCL9+TMB, T cell GEP, VIGex, cytolysis score, TIDE, TMB alone, and PD-L1 expression). The IOSelect classifier achieved an average AUC of 0.72, significantly outperforming existing biomarkers with average AUCs ranging from 0.60 (TMB alone) to 0.68 (T cell GEP), such as FIG. 6A shown.

[0109] Next, an evaluation was performed to determine whether a simple model containing only the top few predictive features was sufficient to robustly predict ICB response. To identify the most predictive features for ICB response, the embodiments used LASSO logistic regression and 10-fold cross-validation to robustly select features from the set of 59 stromal DEGs and TMB, as shown in FIG. 6B In some embodiments, a response estimation model was constructed using the top three features selected by LASSO (TMB, IFNG, and FCER1A expression), as shown in FIG. 6C The simple three-feature model achieved a mean AUC of 0.70 across the 6 test cohorts, only slightly lower than the AUC achieved by the full IOselect model shown in FIG. 6D FIG. 7 In addition, the three-feature classifier could predict Progression Free Survival (PFS) in cohorts with survival data, as shown in

[0110] Tumor transcriptome deconvolution using in silico simulation estimated the mean cancer and stromal expression of each gene for responder / MSI and non-responder / MSS tumors in each cancer cohort, respectively. A set of DEGs that were shared across multiple cancer types was then identified. The experiments found 59 consensus DEGs that were specific to stroma, which were highly enriched in immune-related pathways and processes, as shown in FIG. 7 The top 2 stromal genes overexpressed in ICB responders were CXCL9 and CXCL13. Multiple recent studies have demonstrated the role of CXCL9 and CXCL13 in antitumor immunity and their association with good ICB response.

[0111] Previous studies have suggested that CXCL13 can play a dual role in ICB response. First, CXCL13 is produced by Tfh cells and involved in B cell organization in TLS, the presence of which can serve as a positive predictor of ICB response. Second, CXCL13 was reported to be highly expressed in neoantigen-reactive CD8+ T cells that directly participate in cancer cell killing. Based on these findings, we found that both Tfh and CD8+ T cells were enriched in tumors of ICB responders compared to non-responders across various cancer types (as shown in FIG. 26 and FIG. 26A In addition, the abundance of these two cell types was associated with CXCL13 expression. Consistent with recent studies, we found that CXCL9 expression was closely associated with Ml macrophage infiltration, and Ml macrophages were one of the most differentially enriched immune cell types in ICB responders, as shown in FIG. 26B ​Although myeloid infiltration has traditionally been associated with immunosuppression, our data suggest that M1 macrophage infiltration and macrophage polarization can be important determinants of tumor sensitivity to ICBs.

[0112] Furthermore, a set of immune-related genes that are negatively correlated with ICB response has been identified. Among these genes, FCER1A and CH25H were previously reported to be overexpressed in anti-inflammatory M2 macrophages, whereas FCER1A+ TAMs were reported to promote tumor progression by participating in a positive feedback loop with tumor-initiating cells.

[0113] It has been observed that resting CD4 memory T cells, dendritic cells, and mast cells are enriched in ICB nonresponders. FIG. 7 As shown, the activated state of these cell types is enriched in tumors of ICB responders. Furthermore, although FCER1A is expressed in antigen-presenting cells such as dendritic cells and macrophages, it is noteworthy that FCER1A expression is 25-fold higher in resting dendritic cells than in activated dendritic cells, and 14-fold higher in M2 macrophages than in M1 macrophages. These results suggest that a pre-existing immunosuppressive microenvironment associated with resting immune cells and high FCER1A expression prevents effective antitumor immune responses in ICB non-responders.

[0114] Analysis of immune checkpoint and cytokine ligand–receptor interactions revealed coordinated upregulation of several ligand–receptor pairs in the tumor stroma of ICB responders, including genes encoding PD-L1 / PD-L2 and their receptor PD-1, CXCL9 / CXCL10 and their receptor CXCR3, and CCL5 and its receptor CCR5, as shown in Figures 5 and FIG. 4AIn contrast, consistent changes in expression of checkpoint proteins and cytokines were not found in cancer cells from ICB responders. Among the chemokines upregulated in ICB responders, CXCL9, CXCL10, and CXCL13 promote lymphocyte expansion, while CCL5 is a chemoattractant for T cells and myeloid cells. CCL5 also directs CCR5+ blood monocytes and neutrophil expansion and differentiation into immunosuppressive tumor-associated macrophages and neutrophils. Furthermore, CCL5 has direct pro-tumor effects on cancer cells, promoting cancer proliferation and metastasis. Preclinical and clinical studies have shown that targeting the CCL5-CCR5 axis can enhance anti-tumor responses. Since both CCL5 and CCR5 are overexpressed in ICB responders, the findings suggest that targeting CCL5-CCR5 signaling in combination with ICB can further enhance anti-tumor immune responses. In fact, a recent trial found that the combination of anti-CCR5 and anti-PD-1 had little clinical benefit for patients with microsatellite stable colorectal cancer who received a high amount of prior therapy, and more trials are underway to assess the safety and efficacy of anti-CCR5 and checkpoint inhibitor combinations. Currently, most companion diagnostics for PD-1 / PD-L1 inhibitors only measure PD-L1 expression levels in cancer cells, or measure PD-L1 expression levels in both cancer cells and infiltrating immune cells, but rarely measure PD-L1 expression levels in immune cells alone. A recent trial found that PD-L1 expression on immune cells had prognostic effects for head and neck cancer patients treated with ICB. Our study further suggests that PD-L1 expression on immune cells is a more consistent determinant of ICB response across tumor types.

[0115] Using the stromal consensus DEGs identified from the meta-analysis, the IOselect score was developed, which robustly predicted ICB response in 6 independent test cohorts, outperforming existing gene expression signatures of ICB response. Furthermore, we proposed a simple three-feature model that was able to achieve comparable performance by combining a single positive marker of immune infiltration (IFNG), a single negative marker of immune suppression (FCER1A), and TMB. While an ideal pan-cancer biomarker of ICB response should be robust across various patient cohorts, it was noted that the melanoma cohort of Liu et al. did not share many DEGs with the other cohorts, including other melanoma cohorts. This result highlights the significant heterogeneity between ICB cohorts, even within the same tumor type. It is unclear to what extent this heterogeneity is caused by patient selection, differences in ICB regimens, differences in prior treatments, different sample collection protocols, or differences in gene expression assays. This highlights the need for a large, uniformly selected patient cohort and standardized gene expression profiling to further improve biomarker discovery.

[0116] Since cancer cells are under negative selection pressure from the immune system, several immune escape mechanisms have been proposed for cancer cells: (i) cancer cells downregulate antigen presentation by deleting HLA alleles, (ii) cancer cells deplete neoantigens to avoid detection, (iii) oncogenic signals can upregulate PD-L1 expression and suppress T cell growth. However, the prevalence of these immune evasion mechanisms remains unclear. A recent meta-analysis found no association between loss of HLA allele heterozygosity and ICB response across ICB cohorts. Furthermore, a pan-cancer analysis failed to detect neoantigen depletion signals in most tumor types. Similarly, we found a lack of conserved intra-tumoral immune evasion mechanisms at the transcriptomic level, challenging the prevailing dogma that checkpoint ligand intra-tumoral expression can modulate immune activity and predict ICB response across cancer types. These insights have profound implications for biomarker assays and the development of improved in vitro model systems for ICB drugs. These results suggest that cancer cells employ multiple immune evasion mechanisms, and cancer cells do not employ a universal strategy to avoid immune-mediated killing. The growing number of molecularly characterized ICB cohorts will help further explore the importance and prevalence of different intra-tumoral immune evasion strategies in determining ICB response.

[0117] Gene expression data collection and normalization

[0118] Expression data for TCGA colorectal (COAD, READ), gastric (STAD), and endometrial (UCEC) tumors were downloaded from the UCSC Xena repository. Raw genomic and transcriptomic data were obtained from the European Nucleotide Archive and processed using the bcbio next generation pipeline. Briefly, RNAseq reads were aligned to the human transcriptome using STAR53 and quantified for TPM counts using SALMON54. Transcriptomic data for processed ICB cohorts were also obtained from publicly available sources. Genes with expression of 0 in more than 50% of samples were removed. All gene expression data were log-transformed and upper quantiles were normalized across all cohorts.

[0119] Tumor purity assessment

[0120] Purity estimates for TCGA tumors were obtained from publicly available sources mentioned above, where consensus purity was estimated from mutation allele frequency, copy number variation, and RNA profiles for each tumor. Tumor purity estimates for samples from the ICB cohort were estimated using PUREE, a machine learning method to estimate tumor purity from bulk RNA transcriptomics. For records where genomic data were also available, consensus purity estimates were based on both genomic and transcriptomic data. Purity estimates were computed using 2 RNA-based methods (PUREE and Estimate) and DNA-based methods (Sequenza, PurBayes). The purity estimates from each method were then quantile-quantile standardized, and the mean standardized purity estimate was selected as the consensus purity.

[0121] MSI status and ICB response classification

[0122] MSI status for TCGA tumors was obtained from two prior pan-cancer studies of microsatellite instability. EBV status for TCGA STAD tumors was obtained from a prior TCGA study of gastrointestinal cancers. Clinical annotations of ICB response based on RECIST criteria were obtained from published studies.

[0123] Expression deconvolution

[0124] Bulk tumor expression was modeled as the sum of cancer and stromal expression, weighted by cancer purity. This can be written as:

[0125] y j = c j p + s j (1 - p) + ε j

[0126] where y j is an n x 1 vector of bulk tumor expression of gene j in n samples, c j and s j are the average cancer-specific and stromal-specific expression of gene j, p is an n x 1 vector of cancer purity, and ε j is the residual. Given bulk expression (y j ) and tumor purity (p), non-negative least squares (NNLS) regression can be used to estimate cancer-specific expression (c j ) and stromal-specific expression (s j ) for each gene. c j and s j represent the estimates of average cancer-specific and stromal-specific expression of gene j when p c = 1 and p s = 0, respectively. The standard error of the prediction for gene j is computed as:

[0127]

[0128] where p h = p c = 1 represents cancer-specific expression, p h = p s = 0 represents stroma-specific expression; MSE is the mean square error.

[0129] Permutation statistics to identify differentially expressed genes

[0130] To test for differences in cancer (or stroma) expression between two response groups (e.g., c Rj -c NRj ) for each group and compute the test statistic as:

[0131]

[0132] To assess the significance of the null hypothesis that there is no difference in cancer (stroma) expression between two response groups, a null distribution is generated by permutation, where the group labels are first randomly shuffled 10,000 times and the above procedure is applied to each permuted dataset, thereby generating a null distribution of statistics. The p-value is computed as:

[0133]

[0134] i.e., the probability of observing or more, given that the null hypothesis is true.

[0135] Identification of consensus DEGs across cancer types

[0136] For each gene and each tumor compartment (cancer or stroma), a meta-p-value is computed. This computation can include combining permutation p-values across all cancer cohorts using the Fisher method. The examples use a Bonferroni correction to adjust for multiple testing to maintain an overall a of 0.01. Genes that are differentially expressed in at least 2 cohorts in different directions are removed. Finally, to eliminate DEGs that are involved in MSI biology but not directly related to ICB response, the final consensus DEG list requires differential expression (nominal p-value < 0.1) in at least 1 ICB cohort.

[0137] Pathway enrichment analysis

[0138] To identify whether the identified DEGs are enriched in specific biological processes, pathway enrichment analysis was performed using ClueGo 60 with default parameters. Briefly, the degree of enrichment of a given list of DEGs in GO biological pathway terms was computed using two-sided hypergeometric statistics and the resulting p-values were corrected for multiple testing using the Bonferroni step-down method. To aid interpretation, important biological terms that were highly overlapping (measured by the kappa coefficient) were merged into functional groups.

[0139] To compute the degree of enrichment of stromal DEGs in immune-related functions, immune-related genes were defined as genes in the GO biological process data set in the “immune system process” term and its children. In addition, the curated stromal DEGs were not manually labeled as immune genes by GO and FCER1A, GZMH, FCRL6, SIRPG, JAKMIP1, and WARS were further labeled as immune-related genes based on literature search.

[0140] Analysis of single-cell RNAseq data

[0141] Gene-level normalized counts from Bi et al. were downloaded from the Single Cell Portal. A total of 34,327 cells from 8 renal cell carcinoma tumors were analyzed, of which 17,672 cells from 5 tumors were from patients who received ICB treatment. Cells with low library size and high mitochondrial contamination were removed. For each annotated cell type, the mean expression and the proportion of expressing cells were computed for each DEG separately for cells of ICB responders and non-responders. To generate the dot plots in FIG. 6, the expression counts of each gene were scaled using z-score normalization and the mean scaled expression and the expression prevalence of each gene in all stromal (non-cancer) cell types were computed. ​

[0142] Predictive model for ICB response

[0143] The immune score determined using the present response estimation model was defined as the ratio of the mean expression of the top 10 DEGs upregulated in ICB responders (CXCL9, CXCL13, IFNG, KLRD1, GBP1P1, CRTAM, FASLG, GBP5, CD8A, PRF1) to the mean expression of all 8 DEGs downregulated in ICB responders (FCER1A, PID1, CH25H, NPR3, PREX2, NPR1, FAM134B, SDPR). To establish the response estimation model, RNAseq and TMB from 461 patients from 3 discovery cohorts (Mariathasan et al., Liu et al., Kim et al.) were used as training data.

[0144] ​FASTQ files for whole exome sequencing of the Kim et al. gastric ICB cohort were downloaded from the European Nucleotide Archive and processed using the bcbio next generation pipeline. TMB was calculated as the number of non-synonymous mutations per megabase. For the Mariathasan et al. and Liu et al. cohorts, TMB was obtained from the original publication. All TMB values were standardized using z-score transformation (mean = 0, standard deviation = 1). A logistic regression model was fitted using the immune score generated using the reaction estimation model and TMB as input features.

[0145] Next, experiments were performed to construct more compact models with comparable performance. Using the same training data, LASSO was applied to select top features (“top” is a predetermined number, e.g., 4 or 10) from the 59 stromal DEGs and TMB. LASSO penalizes the sum of the absolute values of the regression coefficients, shrinking some coefficients to zero. Thus, LASSO yields a simpler and more interpretable model, where terms associated with zero coefficients or coefficients with values below a predetermined threshold can be excluded. The regularization parameter, l, was selected by 10-fold cross-validation such that the error of the selected model was within 1 standard deviation from the minimum error. LASSO regression and cross-validation were performed using the “glmnet” package in R. To select the most robust features, 100 samples were bootstrapped using 90% of the data in each bootstrap, and features discovered by at least 75% of the bootstraps were selected for the final model. TMB, IFNG, and FCER1A were found to be the most stably selected features. A logistic regression model was then built using these three features on all 461 training samples.

[0146] Benchmarking of ICB prediction models against existing biomarkers

[0147] T-cell GEP scores were calculated from normalized RNAseq data using weights provided by Ayers et al. in 2018. TIDE scores were calculated from normalized RNAseq data using the TIDE web platform. A CXCL-9+TMB model was trained using logistic regression on the training cohort defined herein. Cytolytic score was calculated as the geometric mean of GZMA and PRF1 expression. Vigex score was calculated as the mean of z-score scaled expression of 12 genes (CXCL9, CXCL10, CXCL11, IFNG, PRF1, IL7R, GZMA, GZMB, PDCD1, CTLA4, CD274, FOXP3).

[0148] Immune cell enrichment analysis

[0149] Immune deconvolution was performed using CibersortX (with LM22 signature matrix) to estimate the absolute abundance of immune cell subsets for each tumor. Wilcoxon rank-sum test was used to identify immune cell subsets that differentially infiltrated between responders and non-responders.

[0150] The reference in this specification to any prior publication (or information derived from it), or to any matter which is known, is not, and should not be taken as, an acknowledgement or admission or any form of suggestion that that prior publication (or information derived from it) or known matter forms part of the common general knowledge in the field of endeavour concerned.

[0151] Throughout this specification and the claims which follow, unless the context requires otherwise, the word "comprise", and variations such as "comprises" and "comprising", will be understood to imply the inclusion of a stated integer or step or group of integers or steps but not the exclusion of any other integer or step or group of integers or steps.

[0152] The scope of the disclosure encompasses all changes, substitutions, variations, alterations, and modifications that a skilled person could understand as being made to the exemplary embodiments described or illustrated herein. The scope of the disclosure is not limited to the exemplary embodiments described or illustrated herein. Furthermore, although the disclosure describes and illustrates various embodiments herein as comprising particular components, elements, features, functions, operations or steps, any one of these embodiments can comprise any combination or permutation of any of the components, elements, features, functions, operations or steps described or illustrated anywhere herein. Although the disclosure describes or illustrates particular embodiments as providing particular advantages, particular embodiments can not provide the advantages described, or may provide some or all of the advantages described.

Claims

1. A clinical decision support system for predicting the clinical response of a target tumor to immune checkpoint blockade therapy (ICB), the system comprising one or more processing units for: (i) Generate a response estimation model by: (a) obtaining a training dataset comprising training transcriptome sequence records of tumor tissue samples from multiple known responders and multiple known non-responders to ICB; (b) performing deconvolution on the training dataset to identify differentially expressed genes (DEGs) in the training dataset, wherein the DEGs are considered as features of the training dataset, and the associated responder and non-responder statuses are considered as labels of the records in the training dataset; (c) performing regression analysis on the features to select a feature subset that best predicts the responder or non-responder label and obtain associated feature weights; (d) incorporating the feature subset and associated feature weights into the response estimation model; (ii) receiving transcriptome sequence data of a target tumor tissue sample; (iii) processing the received transcriptome sequence data using the generated response estimation model to generate an estimated response index indicating the likely response of the target tumor to ICB.

2. The system according to claim 1, wherein: The DEGs were identified based on the p-value of each candidate gene j and calculated as: in is the null distribution estimated by randomly perturbing the responder or non-responder labels in the training dataset to obtain permutations of the training dataset, and calculating for each permutation of the training dataset To estimate c Rj is the estimated expression level of candidate gene j in the responder records of the training dataset, with the corresponding standard error SE Rj , c NRj is the estimated expression level of candidate gene j in the non-responder records of the training dataset, with the corresponding standard error SE RNj , the standard error of candidate gene j is calculated as: Among them, p h =p c =1 represents cancer-specific gene expression, p h =p s = 0 represents matrix-specific gene expression, MSE is the mean squared error, and n is the number of records in the training dataset.

3. The system according to claim 2, wherein: c was estimated by modeling bulk tumor expression as the sum of cancer and stromal expression Rj and c NRj The value of , weighted by cancer purity, is expressed as: y j =c j p+s j (1-p)+ε j , Among them, y j is an n×1 vector of bulk tumor expression of gene j in n samples, c j and s j is the average cancer-specific and stroma-specific expression of gene j, p is an n × 1 vector of cancer purity, ε j is the residual.

4. The system according to claim 1, wherein: The response estimation model also incorporates tumor mutation burden (TMB) signature; and The response estimation model processes the TMB index of the target tumor when generating the estimated response index.

5. The system according to claim 1, wherein: The regression analysis was performed using least absolute shrinkage and the selection operator LASSO.

6. The system according to claim 1, wherein: The one or more processing units are further configured to: receiving exome sequencing data of the target tumor tissue sample; and Clinically actionable genomic mutations associated with ICB were identified by the response estimation model.

7. The system according to claim 1, wherein: The transcriptome data of the target tumor tissue sample includes transcriptome data of stromal cells in the tissue sample.

8. A system for generating a response estimation model to predict the clinical response of a tumor to immune checkpoint blockade therapy (ICB), the system comprising one or more processing units for: Obtaining a training dataset, the training dataset comprising training transcriptome records of tumor tissue samples of multiple known responders and multiple known non-responders to ICB; performing deconvolution on the training dataset to identify differentially expressed genes (DEGs) in the training dataset, wherein the DEGs are considered as features of the training dataset, and the associated responder and non-responder statuses are considered as labels recorded in the training dataset; performing a regression analysis on the features to select a subset of features that best predicts the responder or non-responder label; The feature subset is incorporated into the response estimation model.

9. A computer-implemented method for predicting the clinical response of a tumor to immune checkpoint blockade therapy (ICB), the method comprising: (i) Generate a response estimation model by: (e) obtaining a training dataset comprising training transcriptome records of tumor tissue samples from multiple known responders and multiple known non-responders to ICB; (f) performing deconvolution on the training dataset to identify differentially expressed genes (DEGs) in the training dataset, wherein the DEGs are considered as features of the training dataset, and the associated responder and non-responder statuses are considered as labels of the records in the training dataset; (g) performing regression analysis on the features to select a feature subset that best predicts the responder or non-responder signature and obtain associated feature weights; (h) incorporating the feature subset and associated feature weights into the response estimation model; (ii) obtaining transcriptome data of target tumor tissue samples; (iii) processing the transcriptomic data using the generated response estimation model to generate an estimated response index indicating the likely response of the target tumor to ICB.

10. The method according to claim 9, wherein: The DEGs were identified based on the p-value of each candidate gene j and defined as: in is the null distribution estimated by randomly perturbing the responder or non-responder labels in the training dataset to obtain permutations of the training dataset, and calculating for each permutation of the training dataset To estimate c Rj is the estimated expression level of candidate gene j in the responder records of the training dataset, with the corresponding standard error SE Rj , c NRj is the estimated expression level of candidate gene j in the non-responder records of the training dataset, with the corresponding standard error SE RNj , the standard error of candidate gene j is calculated as: Among them, p h =p c =1 represents cancer-specific gene expression, p h =p s = 0 represents matrix-specific gene expression, MSE is the mean squared error, and n is the number of records in the training dataset.

11. The method according to claim 10, wherein: c was estimated by modeling bulk tumor expression as the sum of cancer and stromal expression Rj and c NRj The value of , weighted by cancer purity, is expressed as: y j =c j p+s j (1-p)+ε j , Among them, y j is an n×1 vector of bulk tumor expression of gene j in n samples, c j and s j is the average cancer-specific and stroma-specific expression of gene j, p is an n × 1 vector of cancer purity, ε j is the residual.

12. The method according to claim 9, wherein The response estimation model also incorporates tumor mutation burden (TMB) signature; and The response estimation model processes the TMB index of the target tumor tissue sample when generating the estimated response index.

13. The method according to claim 9, wherein: The regression analysis was performed using least absolute shrinkage and the selection operator LASSO.

14. The method according to claim 9, further comprising: Receiving exome sequencing data of the target tumor tissue sample; as well as Clinically actionable genomic mutations associated with ICB were identified by the response estimation model.

15. The method according to claim 9, wherein The transcriptome data of the target tumor tissue sample includes transcriptome data of stromal cells in the tissue sample.

16. A computer-implemented method for generating a response estimation model to predict the clinical response of a tumor to immune checkpoint blockade therapy (ICB), the method comprising: Obtaining a training dataset, the training dataset comprising training transcriptome records of tumor tissue samples of multiple known responders and multiple known non-responders to ICB; performing deconvolution on the training dataset to identify differentially expressed genes (DEGs) in the training dataset, wherein the DEGs are considered as features of the training dataset, and the associated responder and non-responder statuses are considered as labels recorded in the training dataset; performing a regression analysis on the features to select a subset of features that best predicts the responder or non-responder label; The feature subset is incorporated into the response estimation model.

Images