The interaction between cancer-associated fibroblasts (CAFs) and the tumor microenvironment (TME) is a key factor for promoting tumor progression. In lung cancer, the crosstalk between CAFs and malignant and immune cells is expected to provide new directions for the development of immunotherapy. In this study, we have systematically analyzed a single-cell dataset and identified interacting genes between CAFs and other cells. Subsequently, a robust fibroblast-related score (FRS) was developed. Kaplan-Meier (KM) and ROC analyses showed its good predictive power for patient prognoses in the training set comprising of specimens from the cancer genome atlas (TCGA) and in three external validation sets from the Gene Expression Omnibus (GEO). Univariate and multivariate Cox regression analyses suggested that FRS was a significant prognostic factor independent of multiple clinical characteristics. Functional enrichment and ssGSEA analyses indicated that patients with a high FRS developed “cold” tumors with active tumor proliferation and immunosuppression capacities. In contrast, those with a low FRS developed “hot” tumors with active immune function and cell killing abilities. Genomic variation analysis showed that the patients with a high FRS possessed a higher somatic mutation burden and copy number alterations and were more sensitive to chemotherapy; patients with a low FRS were more sensitive to immunotherapy, particularly anti-PD1 therapy. Overall, these findings advance the understanding of CAFs in tumor progression and we generated a reliable FRS-based model to assess patient prognoses and guide clinical decision-making.

Introduction

Lung cancer has the highest incidence among all cancer types and is the leading cause of cancer-related deaths (1); lung adenocarcinoma (LUAD) is its most common histological type. Several epidemiological investigations and experimental studies have attributed the onset and progression of LUAD primarily to environmental factors and genetic alterations (2–4). Given a large number of non-smokers with LUAD, previous theories based solely on environmental factors have been disproven and research attention has been re-focused on profound alterations in the genetic content. To date, there are two main genetic factor-related treatment strategies, namely, targeted therapy and immunotherapy (5). However, most patients who receive targeted therapy are prone to resistance, and only a minority of them may benefit from immunotherapy. Therefore, it is crucial to develop robust tools for prognostic prediction and assessment of treatment responses to further facilitate accurate diagnoses and devise individualized treatment strategies.

Tumor microenvironment (TME) is defined as the environment surrounding the tumor, including the extracellular matrix, immune cells, and stromal cells, all of which are closely associated with tumor progression and treatment outcomes (6). Accumulating evidence elucidate the role of TME infiltration in immune therapeutic responses and resistance against different cancer types; these studies have also investigated their impact on patient prognoses (7, 8). Previous studies have focused more on immune cells. However, several findings have now highlighted the importance of stromal cells in tumor progression (9, 10). Cancer-associated fibroblasts (CAFs), a representative component of stromal cells, play crucial roles in cancer genesis, progression, and invasion (11, 12). Recently, the interaction between CAFs and the tumor immune microenvironment (TIME) has been identified as a key factor in promoting tumor progression (13, 14). CAFs interact with immune cells and other immune components within the TIME through various secreted cytokines, growth factors, and chemokines, resulting in an immunosuppressive TME that allows cancer cells to evade the surveillance mechanisms of the immune system (14, 15). Therefore, further investigation into the crosstalk between CAFs and TME is expected to provide new strategies for LUAD treatment, in particular for immunotherapy.

In this study, we used the single-cell dataset, GSE131907, to evaluate the crosstalk between CAFs and other cells. In addition, receptor-ligand pairs were systematically identified for interactions of CAFs with other cells. Based on these receptor-ligand genes, we generated the fibroblast-associated score (FRS) using the LASSO algorithm in the TCGA-LUAD cohort and the GEO meta-cohort to predict patient prognoses and estimate their sensitivity to chemotherapy and immunotherapy. Additionally, the associations among FRS, biological functions, TIME, and genomic alterations were systematically assessed. In summary, our findings are expected to advance the understanding of CAF functions in cancer as we have constructed and described here a robust scoring system to accurately predict patient prognoses and guide clinical decision-making.

Material and methods

Data extraction from online databases

The single-cell transcriptome dataset, GSE131907, was extracted from the GEO database (https://www.ncbi.nlm.nih.gov/geo/), consisting of data of 58 sequences from 44 patients. Next, these data were processed using the 10x Genomics method. Of these, we selected 29 normal lung tissues and early, advanced, and brain-metastasized lung tissues for further analyses. Detailed data processing procedures and ethical approval have been described previously (16).

The data of transcriptome RNA sequencing, Mutect2 mutation, HumanMethylation450 array, copy number variations (CNVs), and the corresponding clinical information were downloaded from TCGA database (https://cancergenome.nih.gov/) using the GDC API. A total of 492 LUAD samples were collected after the exclusion of patients with missed visits and incomplete clinical information. The raw FPKM sequencing data were normalized by TPM and used as the training cohort. Three mature LUAD cohorts were collected from GEO, including dataset GSE30219 from the Affymetrix HG-U133 Plus 2.0 Array platform, dataset GSE72094 from the Rosetta/Merck Human RSTA Custom Affymetrix 2.0 platform, and dataset GSE42127 from the Illumina HumanWG-6 v3.0 expression bead chip. To prevent batch effects on these chips, we merged the three GEO datasets and normalized the data by the log2 transformation using the combat function of the “sva” package (17). Subsequently, LUAD meta-data containing the complete clinical information of 615 individuals were used as the validation cohort. Additionally, we collected the publicly available immunotherapy cohorts with complete clinical information and transcriptomic data. Finally, the information of a cohort of advanced uroepithelial carcinoma treated with anti-PD-L1 immunotherapy (Imvigor210) consisting of 298 patients (8) and a cohort of non-small cell lung cancer (NSCLC) of 27 patients treated with PD1 (GSE135222) was collected.

Single-cell data analysis

The R package, “Seurat”, was used to process the scRNA-seq data. In addition, cells with “min.cells < 3” and “min.features < 200” were excluded. After filtering out the cells with > 60% mitochondrial sequencing count and nFeature_RNA > 7000, a total of 47822 cells were retained for subsequent analyses. The dataset was then normalized using the NormalizeData and ScaleData functions in Seurat. Cell types were identified according to the cell annotations provided in the original article.

To unravel the changes in the cell clusters during tumor progression, we used the R package, ‘monocle’, for the analysis of the single-cell trajectory. Subsequently, single-cell developmental trajectories were identified using the top 1500 variable genes (18). The Python package, “CellphoneDB”, was used to identify receptor-ligand exchanges between cell clusters; the receptor-ligand interactions between eight-cell clusters were thus identified at the molecular level (19). Receptor-ligand pairs with p-values < 0.05 were screened to assess the molecular interaction network among CAFs and other cells. Corresponding interacting genes were identified as fibroblast-related genes (FRGs). Finally, the GGplot2 package was used to visualize these results.

Construction and validation of the FRS model

LUAD-TCGA cohort was used to train the model. Specifically, independent prognostic factors among FRGs were first screened by univariate Cox regression, and genes with P < 0.05 were included for further analysis. Subsequently, a Cox proportional risk model with LASSO penalties was used to identify the best prognostic model. To prevent overfitting, a five-fold cross-validation process was set up. Considering random sampling for cross-validation, 300 iterations were performed to identify the most stable prognostic model. The model with the highest frequency of occurrence in the 300 iterations served as the final prognostic model. Finally, FRS was calculated according to the following equation:

$$FRS=\sum iCoefficient\left(mRN{A}_i\right)\times Expression\left(mRN{A}_i\right)$$

To assess the predictive power of the risk scores in the training and validation sets, the consistency index (C-index) was calculated using the “survcomp” R package, with a larger C-index indicating a more accurate predictive power of model (20). Patients were classified into high- and low-risk groups based on the median FRS. Furthermore, the prognostic value of the risk model was systematically assessed using the Kaplan Meir (KM) survival curves, univariate and multivariate Cox regression analyses, and time-dependent ROC curves.

Functional enrichment and immune infiltration analyses

We performed a single-sample gene set enrichment analysis (ssGSEA) based on the previously published molecular markers using the R package, “gsva”, to assess the activities of the biological pathways for the samples, including angiogenesis, epithelial-mesenchymal transition (EMT), myeloid inflammation, and molecular markers for other immune-related pathways (21–24). Molecular markers for hypoxia were collected from Msigdb (25). Detailed pathway-related gene markers are shown in Table S1. Additionally, GSEA was performed between high- and low-FRS groups, and the significant KEGG pathways were screened using the set threshold of P < 0.05. Moreover, functional enrichment of genes was obtained using the Metascape (www.metascape.org/) database.

The abundances of immune cell infiltrate in tumor samples were estimated using the R package, “CIBERSORT”, to evaluate the degree of infiltration of 22 immune cell types (26). The immune activity and tumor purity of the samples were assessed using the Estimate algorithm (27). The immunophenoscores (IPS) of the samples were calculated based on a previous study, with a higher IPS indicating a stronger immune activity of the sample (28). In short, IPS is calculated on a scale of 0-10 based on the transcriptome of the representative genes of the immunophenotype. Samplewise Z scores were positively weighted according to effective immune cells, negatively weighted according to inhibitive immune cells, and then averaged. Z score ≥ 3 is defined as IPS10, and Z score ≤ 0 is defined as IPS0.

Finally, homologous recombination deficiency (HRD) scores, indel neoantigens, and SNV neoantigens of the samples were obtained from Thorsson et al. (29).

Comparison of genomic variation landscapes between two groups

To compare the differences in mutation burdens between the two groups, the mutation data were processed using the ‘maftools’ package in R. The total number of mutations in the samples was first calculated, and genes with a minimum number of mutations > 30 were identified. The differences in mutation frequencies between the high- and low-FRS groups were then compared using a chi-square test and visualized using maftools (30). CNV data were processed using the GISTIC 2.0 webtool in Genepattern. Subsequently, significantly amplified and missing chromosomal segments were identified and differences in CNVs on the chromosomal arms were assessed. Finally, these CNV results were visualized using the R package, “ggplot2”.

Clinical significance of the risk model

The five most commonly used first-line drugs, including cisplatin, docetaxel, gemcitabine, paclitaxel, and vinorelbine, were selected for the treatment of LUAD. Ridge regression was used to calculate the half-maximal inhibitory concentration (IC50) for each sample, which was then used to assess the sensitivities of patients to chemotherapy in the high- and low-risk groups. Moreover, the accuracy of these predictions was assessed by a 10-fold cross-validation process (31). Furthermore, differentially expressed genes between the two groups were considered as potential therapeutic targets. The CMap database (https://clue.io/) was used to obtain the potential compounds targeting these genes. This database can not only predict drugs based on the gene expression profiles but also elucidate the mode of action (MoA) of these compounds targeting the corresponding molecular pathways. To assess the patient responses to immunotherapy, the TIDE online tool (http://tide.dfci.harvard.edu) was used (32). In addition, the unsupervised subclass mapping algorithm (https://cloud.genepattern.org/gp/) was used to assess the patient responses to anti-PD1 and anti-CTLA-4 immunotherapeutic regimens. Finally, we validated the predictive efficacy of FRS in the immunotherapy cohorts, Imvigor210 and GSE135222.

Clinical specimens

We obtained 50 tissue specimens from patients who received surgical resection of primary LUAD in Shanghai Pulmonary Hospital from September 2015 to April 2016, including 27 males and 23 females, with a mean age of (66.24 ± 7.3) years. And all patients were followed up every three months for five years. Inclusion criteria: 1, all were diagnosed as lung adenocarcinoma by postoperative pathological examination; 2, all did not receive radiotherapy or chemotherapy before surgery; 3, clinical data were complete. Exclusion criteria: 1, combined with chronic systemic diseases; 2, combined with other malignant tumors. Written informed consents of all patients were obtained before the study. The study was approved by Shanghai Pulmonary Hospital Ethics Committee (ethical lot number: K21-111Y).

Immunohistochemistry staining

After obtaining the tumor tissue, the tissue was routinely paraffin-embedded and preserved. In the experiment, tissue sections were dewaxed and rehydrated, and antigen repair was performed by incubating the slides in 10 mmol/L sodium citrate buffer and microwave treating the samples for 20 min. After being closed with 3% H2O2 and 10% normal goat serum (NGS), the slides were incubated with primary antibody at 4°C overnight. The paraffin-embedded LUAD sections were incubated with anti-TNFSF14 (ab115544, ABCAM), anti- JAM2 (ab156586, ABCAM), anti- LIFR (ab202847, ABCAM)、anti-SPN(ab101533, ABCAM)、anti- HGF(ab118871, ABCAM). The slides were then incubated with biotin-coupled anti-rabbit secondary antibody (1:1000, ab205718, Abcam, UK) for 2h at 37°C using the ABC kit from Vector Laboratories (Burlingame, CA, USA).

The sections were then incubated with polymeric HRP reagent and peroxidase activity was observed with diaminobenzidine tetrahydroxyl chloride solution (Vector Laboratories), and the sections were re-stained with hematoxylin. Ouant center is the analysis software of Pannoramic viewer. When the images of the tissue microarrays are scanned, the TMA software in Ouant center sets the numbers that correspond to the arrangement of the tissue sections. Thereafter, the densito quantification software in Quant center automatically identifies all dark brown dots in the microarray tissue as strongly positive, tan dots as moderately positive, light yellow dots as weakly positive, and blue cell nuclei only as negative, and analyzes the percentage of each stained (strong, moderate, weak, and negative) area in pixels, and finally performs an H-Score score.

Bioinformatic and statistical analyses

All statistical analyses and graph plotting were performed using the R software (version: 4.04). Comparisons between the two groups were made using the Wilcoxon test and differences in proportions were compared using the chi-square test. The KM plotter was used to generate the survival curves and statistically significant differences were assessed using the log-rank test. Time-dependent ROC curves (tROC) were plotted using the R package ‘survivalROC’. Univariate and multivariate COX regression analyses were performed using the R package, “survival”. Additionally, ‘rms’ was used to plot nomograms and calibration curves, and decision curve analysis (DCA) was performed using the “DCA” package (33). Unless specified otherwise, two-tailed P < 0.05 denoted statistical significance.

Results

CAF clustering and identification of FRGs

In order to explore the cross-talk between CAF and other cells and identifying the FRGs. We first analyzed the dataset, GSE131907, at single-cell resolution and identified a total of eight-cell clusters according to their original annotation (Figure 1A). The pseudo-time analysis suggested that CAFs were mainly aligned at the beginning of the trajectory (Figure 1B). Subsequently, the communication network between the eight-cell clusters was analyzed (Figure 1C). Specifically, CAFs were found to communicate the most with endothelial cells, followed by myeloid cells (Figure 1D). Significant receptor-ligand pairs were obtained as FRGs for subsequent analysis based on a set threshold of P < 0.05. The Dot plot was used to visualize the top five receptor-ligand pairs between CAFs and other cells (Figure 1E). Detailed results were shown in Table S2.

FIGURE 1

Figure 1 The single-cell profile of CAFs. (A) Eight cell clusters identified by cell clustering; (B) Pseudo-time analysis of the eight-cell clusters. Upper panel: cell cluster distribution, lower panel: pseudo-time distribution; (C) Communication network of the eight-cell clusters; (D) Communication network of CAFs with other cells, wherein the numbers represent the number of receptor-ligand pairs; (E) The top 5 receptor-ligand pairs of CAFs communicating with other cell types.

We then focused on FRGs with independent prognostic values. Univariate Cox analysis showed 33 independent prognostic factors in 127 FRGs. The loop graph in Figure 2A shows the correlation network and hazard ratios (HRs) for these 33 FRGs. Figure 2B displays the mutational landscape of the 37 FRGs. EGFR and COL5A2 were the top two genes with the highest mutation frequencies. The most common mutation was the missense mutation, whereas single nucleotide point mutation was the most common type of mutation, with the most frequently occurring change being from cytosine to adenine. The waterfall plot in Figure 2C shows the mutational landscape of 16 FRGs in patients. The bar chart displays the CNV profile of the 37 FRGs in TCGA-LUAD. Furthermore, FRGs underwent extensive CNV events. LAMC1 and TNFSF14 were the genes that experienced the most amplification and deletion events, respectively (Figure 2D). The loop graph was plotted to visualize the overall CNV profile of the 37 FRGs on the chromosomes (Figure 2E).

FIGURE 2

Figure 2 Genomic profile of FRGs in LUAD. (A) Correlation network of FRGs; (B) Summary of mutational events for FRGs in TCGA-LUAD; (C) Oncoplot showing the mutational mapping of FRGs; (D) Summary of CNV events for FRGs in TCGA-LUAD; (E) A loop graph showing CNV events for FRGs on chromosomes.

Construction of an FRG-related risk model

An FRG-related risk model was constructed using the 37 FRGs with a prognostic value, on which 300 iterations of LASSO regression were performed. Of all the five combinations, the model containing seven genes was found to be the most stable and showed good accuracy in both the training and validation cohorts (TCGA: 0.715; GEO: 0.667) (Figure 3A). This LASSO model was constructed based on the optimal λ value of 0.01608, and FRS was calculated based on the following equation:

$$FRS=\sum iCoefficient\left(mRN{A}_i\right)\times Expression\left(mRN{A}_i\right)$$

FIGURE 3

Figure 3 Construction of the FRG-related risk model. (A) Screening for the best LASSO model. Left panel: frequency of different gene combinations in the LASSO Cox regression model; Right panel: C-index of the best model in both TCGA and GEO cohorts; (B) LASSO coefficients for the 27 model genes; (C) KM survival curves for the high-FRS and low-FRS groups in the TCGA cohort; (D) Survival status and FRS of patients in TCGA cohort; (E) 1-, 3-, 5-, and 8-year ROC curves for FRS in TCGA cohort; (F) tROC curves for FRS and clinical characteristics in TCGA cohort.

shows the LASSO coefficientFigure 3B shows the LASSO coefficients for the model genes, detailed coefficients of 27 FRGs can be found in Table S3. Patients at high- and low-risk were distinguished based on the median FRS. Survival analysis suggested that patients in the high-risk group had significantly lower survival rates relative to those in the low-risk group (Figure 3C; P < 0.0001). Figure 3D shows the distribution of FRS in TCGA cohort and the transcriptional profiles of the model genes. 1, 3, 5, and 8-year AUC values for the model were 0.66, 0.67, 0.68, and 0.70 respectively (Figure 3E). The tROC analysis suggested that FRS and TNM staging were the best predictors (Figure 3F). Subsequently, the predictive efficacy of the model was also assessed in the validation set. The survival analysis suggested that patients in the high-FRS group showed significantly worse survival (Figure S1A; P < 0.0001). The ROC analysis suggested that the model had satisfactory predictive power in the external validation set, with 1, 3, 5, and 8-year AUC values of 0.68, 0.69, 0.69, and 0.71 respectively (Figure S1B). Figure S1C shows the distribution of FRS and model gene expression in the GEO cohort.

Predictive independence of the risk model

We then validate the prognosis value of the FRS model in the TCGA cohort and GEO meta cohort. The relationship between the risk scores and the clinical characteristics of the patients and their prognoses were analyzed using the univariate and multivariate Cox regression analyses. The results of the univariate Cox regression analysis suggested that FRS was an independent prognostic indicator in both the training and validation cohorts (P < 0.0001) (Figure 4A). The results of the multivariate Cox regression analysis showed that FRS remained an independent prognostic factor for overall survival (OS) in both the training and validation cohorts after correcting for other clinical characteristics (P < 0.0001) (Figure 4B). Furthermore, subgroup analysis indicated that FRS remained a reliable prognostic factor in different clinical groups (Figure S2). Therefore, risk scores could serve as a reliable prognostic marker for predicting OS in patients with LUAD. Subsequently, the nomogram was constructed to better assess the risk of patients with LUAD (Figure 4C). The correction curves for the nomogram showed a good 1-, 3-, and 5-year stability and accuracy of the nomogram model (Figure 4D). tROC analysis suggested that the nomogram model was a better predictor relative to the clinical characteristics (Figure 4E). Additionally, a DCA was conducted to assess the decision benefit of the nomogram model. The results showed that this nomogram model was suitable for 1-, 3-, and 5-year risk assessments of patients with LUAD (Figure 4F).

FIGURE 4

Figure 4 Validation of the FRG-related risk model. (A) Univariate Cox regression analysis of OS in TCGA and GEO cohorts; (B) Multivariate Cox regression analysis of OS in TCGA and GEO cohorts. (C) FRS-based nomogram; (D) Calibration curves for the nomogram; (E) tROC curves for the nomogram and clinical characteristics; (F) 1-, 3-, and 5-year DCA curves for the nomogram.

Functional enrichment analysis of FRS

We tried to explain the potential biological logic of the differences in clinical outcomes among high- and low-FRS groups. Therefore, we assessed the correlation between FRS and some typical biological pathways. The heat map was plotted to illustrate the relationship among FRS, biological pathway activities, and clinical characteristics (Figure 5A), and the correlational analysis between FRS and biological pathways is shown on the right panel (Figure 5B). Angiogenesis, myeloid inflammation, and hypoxia were significantly positively correlated with FRS and their levels were markedly higher in the high-FRS group. GO analysis showed that the upregulated genes in the high-FRS group were mainly associated with the cell cycle, mitosis, and cytoskeleton (Figure 5C), whereas those in the low-FRS group were mainly related to antigen presentation and the complement system (Figure 5E). Further, GSEA showed that the cell cycle-related pathways such as the P53 signaling cascade, spliceosome, and DNA repair were significantly enriched in the high-risk group (Figure 5D), whereas antigen presentation, hematopoietic cell lineage, and the JAK-STAT signaling cascade were significantly enriched in the low-risk group (Figure 5F). Thus, these results suggested that tumor angiogenesis and DNA replication were active in the high-FRS group, whereas immune activity and immune cell differentiation were active in the low-FRS group.

FIGURE 5

Figure 5 Functional analysis for FRS. (A) Heat map showing the correlation between FRS, biological pathway activities, and clinical characteristics; (B) Correlational analysis of FRS and biological pathways; (C) Functional enrichment analysis of the upregulated genes in the high group; (D) GSEA enrichment plot showing the five pathways of interest in the high group; (E) Functional enrichment analysis of the upregulated genes in the low group; (F) GSEA enrichment plot showing the five pathways of interest in the low group.

Immune landscape in the risk model

TME plays a dual role in the tumorigenesis and progression of tumor and anti-tumor response. The correlation between FRS and the immune landscape was assessed in further detail. The heat map in Figure 6A demonstrates the relationship of FRS with the Estimate scores, abundances of immune-infiltrating cells, typical immune checkpoints (including CD274, CTLA4, HAVCR2, IDO1, LAG3, and PDCD1), immune active features (including CD8A, CXCL10, CXCL9, GZMA, GZMB, IFNG, PRF1, TBX2, and TNF), and clinical characteristics of the patients. The corresponding correlational analysis is shown on the right side of the heat map (Figure 6B). Tumor purity, M0 macrophages, and T regs were significantly positively correlated with FRS and these levels were significantly elevated in the high-FRS groups. In contrast, the Estimate score, immune score, DC cells, B cells, and monocytes were negatively correlated with FRS and these levels were significantly lowered in the low-FRS groups. Furthermore, the activities of CXCL9, GZMA, IFNG, PRF1, CD8A, CTLA4, TNF, and HAVCR2 were negatively correlated with FRS and enhanced in the low-FRS group. Subsequently, we focused on the four indicators related to tumor-specific antigens, including HRD score (Figure 6C), indel neoantigens (Figure 6D), IP S(Figure 6E), and SNV neoantigens (Figure 6F). The results showed that FRS was negatively correlated with the HRD score, IPS, and SNV neoantigens, and their levels were significantly elevated in the low-FRS group. These results suggested that patients in the low-FRS group experienced more chromosomal instability and had more tumor neoantigens, thereby contributing to a stronger immune system activity. Thus, we inferred that the patients with a low FRS may stand to gain more benefits from immunotherapy (34–36).

FIGURE 6

Figure 6 Immune infiltration analysis for FRS. (A) Heat map showing the correlation of FRS with the Estimate score, immune cell infiltration abundances, immune checkpoint expression, and clinical characteristics; (B) From top to bottom: correlational analysis of FRS with the Estimate score, immune cell infiltration abundances, and immune checkpoint expression. Scatter plot and box plot show the correlation of FRS with (C) HRD score; (D) Indel neoantigens; (E) IPS, and (F) SNV neoantigens.”Red name with * represents upregulated in high-risk score group, and blue name with * represents upregulated in low-risk score group; *p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.001.

Correlation between FRS and somatic variations

Several recent studies indicate that TMB is associated with patient responses to immunotherapy, whereby more somatic mutations may generate more potential mutation-derived antigens that can be recognized by the immune system. Further, the recognition of these antigens with mutant peptides by the immune system can activate immune functions and enhance anti-tumor immunity (37–39). Considering the clinical significance of TMB, we examined the correlation between TMB and FRS. The forest plot showed that the mutational frequencies of ZFHX4, ADAMTS12, TP53, KRAS, TTN, XIRP2, LRP1B, and CSMD3 were significantly greater in the high FRS-group (Figure 7A). The results of the mutation co-occurrence analysis suggested that the mutations in all the eight genes were highly co-occurring (Figure 7B). Correlation analysis showed that all mutation burdens and the non-synonymous mutation burden were significantly positively correlated with FRS and markedly increased in the high-FRS group (Figures 7C, D). Figure 7E details the mutational landscape of the high-frequency mutated genes in patients with LUAD. CNVs cause chromosomal variations differently. Thus, we further evaluated the correlation between FRS and CNV and found an increased frequency of amplifications and deletions in the low-FRS group at the level of the chromosome arm (Figure 7F). The box plots in Figures 7G, H show a significant increase in the chromosomal deletion events and an upward trend in amplification events in the high-FRS group.

FIGURE 7

Figure 7 Genomic variation landscape of FRS. (A)Forest plots showing statistically significant differentially mutated genes between the high- and low-FRS groups; (B) Co-occurrence analysis of differentially mutated genes; (C) Correlation of FRS and all types of mutation burdens; (D) Correlation of FRS with non-synonymous mutation burden; (E) Oncoplot of high-frequency mutated genes between the high- and low-FRS groups; (F) Bar graph showing the CNV events on different chromosome arms in the high- and low-FRS groups; (G) Box plot showing the differences in the number of chromosome amplifications between high- and low-FRS groups; (H) Box plots showing the differences in the number of chromosomal deletions high- and low-FRS groups. *p < 0.05; **p < 0.01.

FRS-related guidance for clinical decision-making

Previous results suggested that patients in different FRS groups have interesting differences in biological function, TME, and genomic variation, which may lead to different responses to chemotherapy and immunotherapy. Differences in patient sensitivities towards chemotherapeutic agents for LUAD were assessed and the results showed that patients in the high-FRS group in TCGA cohort were more sensitive to the commonly used five first-line agents (Figure 8A). The same results were observed in the validation cohort (Figure S1D). Overall, patients in the high-FRS group were more sensitive to chemotherapy. Based on the value of |log2 FC|, the top 300 differential expression genes between the high and low FRS groups were uploaded to the CMap database to search the underlying small molecular drugs. As shown in Figure S3, a total of 47 potential small molecular drugs were identified to target 35 biological process. Differences in the immune landscape and genomic alterations between the two groups suggested that FRS may be associated with immunotherapeutic efficacy. Therefore, we assessed the patient response rates to immunotherapy using the TIDE algorithm. The results showed a higher response rate to immunotherapy in the low-FRS group in TCGA cohort (Figure 8B; P = < 0.001). In the validation cohort, patients in the low-FRS group also responded substantially more to immunotherapy (Figure S1E; P < 0.001). The results of subclass mapping suggested that the patients in the low-FRS group were more sensitive to anti-PD1 therapy in both TCGA and GEO cohorts (TCGA: FDR = 0.011; GEO: FDR = 0.027) (Figure 8C; Figure S1F). Subsequently, we evaluated the prognostic performances of FRS in an immunotherapy cohort of NSCLC. The results showed that patients in the high-FRS group had a worse survival (Figure 8D; P = 0.078). Finally, we evaluated the utility of FRS in a large immunotherapy cohort, which also suggested that the patients in the high-FRS group had a significantly worse survival (Figure 8E; P = 0.00038),. Overall, these results demonstrated that the risk model constructed in this study was a powerful tool to guide decisions related to chemotherapy and immunotherapy for the treatment of LUAD.

FIGURE 8

Figure 8 Role of the FRS-related risk model in guiding clinical treatment decisions. (A) Box plots showing the predicted IC50 values for the five most commonly used drugs in high- and low-FRS groups; (B) Predicted responses to immunotherapy for patients in the high- and low-FRS groups using the TIDE algorithm; (C) Sensitivity of patients in the high- and low-FRS groups to PD1 and CTLA4 treatment regimens predicted by subclass mapping; (D) KM survival curves for patients in the high- and low-FRS groups in GSE135222 cohort; (E) KM survival curves for patients in the high- and low-FRS groups in IMvigor210 cohort.

Validation of key FRGs in the clinical samples

We extracted the most representative top 5 genes according to lasso coefficient for external validation. The staining intensity of TNFSF14, JAM2, HGF, SPN, and LIFR in the tumors of 50 lung adenocarcinoma patients was first analyzed by immunohistochemistry and quantified according to H-scores. Subsequently, they were defined as the high expression group (H scores > median value) and low expression group (H scores < median value) according to the median value of H-scores, respectively. Subsequently, we performed a prognostic analysis of their Kaplan-Meier according to staining intensity, and we found that patients with higher staining intensity in TNFSF14 had a significantly worse prognosis and had shorter survival cycles (Figure 9A). However, patients with higher staining intensity for JAM2, HGF, and LIFR had significantly higher survival cycles than those with lower staining intensity for lung adenocarcinoma (Figures 9B, D, E), but for SPN, there was no significant correlation between their staining intensity and patients’ survival cycles (Figure 9C).

FIGURE 9

Figure 9 Immunohistochemical analysis of the key FRGs in the FRS model. High TNFSF14 expression, low JAM2, HGF, and SPN expression were associated with poor prognosis in patients with LUAD. (A) Immunohistochemical analysis of the intensity of TNFSF14 staining in tumors from 50 patients with lung adenocarcinoma and Kaplan-Meier analysis of the correlation between H-scores of immunohistochemistry for TNFSF14 and the survival cycle of patients with LUAD; (B) Immunohistochemical analysis of the intensity of HGF staining in tumors from 50 patients with lung adenocarcinoma and Kaplan-Meier analysis of the correlation between H-scores of immunohistochemistry for HGF and the survival cycle of patients with LUAD. (C) Immunohistochemical analysis of LFR staining intensity in tumors from 50 lung adenocarcinoma patients, Kaplan-Meier analysis of H-scores of immunohistochemistry for LFR correlated with the survival cycle of LUAD patients; (D) Immunohistochemical analysis of JAM2 staining intensity in tumors from 50 lung adenocarcinoma patients, Kaplan-Meier analysis of JAM2 H-scores of immunohistochemistry correlated with the survival cycle of LUAD patients; (E) immunohistochemistry analyzed the intensity of SPN staining in the tumors of 50 lung adenocarcinoma patients, and Kaplan-Meier analyzed the correlation of H-scores of immunohistochemistry of SPN with the survival cycle of LUAD patients.

Discussion

Considering the complexity of TME in LUAD patients, previous research attention was focused more on the immune cells, however, the crosstalk between CAFs and other cells remains far less understood. In this study, we used single-cell RNA sequencing data to assess the communication between CAFs and other cells and identified the interacting molecules. Subsequently, FRS was constructed from bulk sequencing data based on interacting genes and its significance in prognostic and therapeutic decision-making was determined. Functional enrichment analysis was employed to understand FRS-related biological functions. Additionally, CIBERSORT, ssGSEA, and ESTIMATE algorithms were used to map the TIME landscape and assess the associations between FRS and TIME by analyzing the LUAD-related genomic information. Finally, the intrinsic associations between FRS and genomic alterations were assessed in terms of tumor mutation burden and CNV effects.

CAF is a pro-tumor stromal cell component in most solid tumors. The interaction between CAFs and various cellular components in TME regulates tumor progression and invasion (13, 14, 40). In this study, we first examined the intercommunication between CAFs and other cells. The pseudotime analysis showed that the distribution of CAFs was mainly at the beginning of cell sorting trajectories, thus suggesting that CAFs were involved in the formation of stromal components in the early stages of cancer progression. Cellular communication analysis revealed additional interactions of CAFs with endothelial, malignant, and myeloid cells. These findings demonstrated that CAFs were not only involved in stromal formation and regulation of tumor progression but also interacted extensively with immune cells. Receptor-ligand analysis suggested that the CAFs could regulate immune cells mainly through the TNF signaling pathway. Subsequently, the significant receptor-ligand pairs were identified as FRGs and an FRG-based FRS model was generated using the LASSO algorithm. This model showed excellent predictive performances in both the training and the external validation cohorts and suggested a significant deterioration in survival among the high-risk patients.

Patients with LUAD have a high tumor mutation burden and show strong immunogenicity. Therefore, LUAD is an ideal indication for immunotherapy (41). However, the overall response rate of patients towards immunotherapy is low and only a certain proportion of patients benefit from it (42). Therefore, identifying “hot” tumors in LUAD is expected to enhance the decision and selection of those who would benefit from immunotherapy. In the present study, functional enrichment analysis suggested that a high FRS was associated with hypoxia, angiogenesis, and myeloid immunity; among them, hypoxia is often considered a limiting factor for TME and can lead to treatment resistance (43). Angiogenesis is essential for tumor growth and metastases, and both angiogenesis and myeloid immunity are inhibitors of the immune system functions (44–47). The activities of most anti-tumor immune and antigen-presentation pathways were markedly increased in patients with a low FRS. The interaction between TIME and immune cells is closely related to immunotherapy and patient prognoses (48–50). We further analyzed the abundances of immune infiltrates in TME and the findings suggested that a high FRS was associated with a higher tumor purity and an elevated Treg level, ultimately leading to immunosuppression (51). In contrast, patients in the low FRS group had higher immune scores, an increased proportion of DC cells, and enhanced immune checkpoint activities. Further, the immunophenotype scores were found to be negatively correlated with FRS and were markedly high in the low-FRS group. These results suggested a low FRS-activated immunophenotype (29, 52), consistent with the better survival of the patients in the low-FRS group and the resultant development of “hot” tumors sensitive towards immunotherapy. Additionally, HRD scores, indel neoantigens, and SNV neoantigens were elevated in the low-FRS group. These findings suggested that more tumor-specific neoantigens may be present in patients in the low-FRS group and they may be more likely to benefit from immunotherapy (53–55).

A recent study shows that genomic alterations are closely related to neoantigen formation and immunotherapeutic responses (56). However, the results of this study suggest that patients in the low-FRS group experience less TMB and that the high-frequency mutated genes were mutated more frequently in the high-FRS group. To elucidate this phenomenon, the mutational co-occurrence of the high-frequency mutated genes was examined. As these were all highly co-mutated genes, patients in the high-FRS group showed a higher TMB frequency. Furthermore, patients in the low-FRS group experienced a higher frequency of CNVs in the chromosomal arms but fewer CNV events in total. These results suggested that FRS could better reflect the immune status of the tumor and predict the patient responses to immunotherapy relative to TMB and CNV.

In summary, low FRS resulted in “hot” tumors with an immune-activating phenotype and possibly the production of more tumor neoantigenic peptides. We then systematically evaluated the patient responses to chemotherapy and immunotherapy. Patients with a high FRS were more sensitive to chemotherapy. Previous functional enrichment results suggested that the cell cycle-related pathways, as targets for chemotherapy, were active in patients with a high FRS, thereby leading to better chemotherapeutic benefits. Subsequently, the TIDE and subclass mapping algorithms predicted a higher patient sensitivity towards anti-PD1 therapy in those with a low FRS, consistent with our previous findings. Moreover, we observed better survival in patients with a low FRS in both the external NSCLC immunotherapy and the large immunotherapy cohorts. Overall, these results demonstrated that the FRS model is a powerful tool that can guide the treatment-decision making for patients with LUAD. Patients with a high FRS are better suited for chemotherapy, whereas those with a low FRS are more likely to benefit from immunotherapy.

However, the present study has some limitations. First, the similarity of expression profiles of CAFs and vascular cells may confound our analysis due to the lack of finer cell classification. Second, bulk sequencing only reflects inter-patient heterogeneity and not intra-tumoral heterogeneity. Third, although we have employed several algorithms to assess the accuracy of this FRS model for predicting patient sensitivity towards chemotherapy and immunotherapy, further validation of these findings by prospective cohort studies and clinical data is required. Finally, additional in vivo and in vitro experiments should be performed to confirm the specific mechanisms underlying the crosstalk of FRGs with other cells in CAFs, which are expected to contribute to the further understanding of the functions of these CAFs.

In summary, this study contributed towards the understanding of cellular interactions in CAFs and TME, and we developed a novel, FRS-based model. This model allowed for the systematic quantification of “cold” and “hot” tumor patterns from multiple perspectives, including function, immune infiltration, and genomic alterations. Moreover, it can also facilitate the quantitative estimation of patient prognoses and guide the clinical decision-making for chemotherapy and immunotherapy.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.

Ethics statement

The study was approved by Shanghai Pulmonary Hospital Ethics Committee (K21-111Y). The patients consented to participate.

Author contributions

SW analyzed the data and wrote the manuscript. XG carried out data interpretations and helped data discussion. WZ conceived and designed the whole project and drafted the manuscript. All authors read and approved the final manuscript.

Funding

This work was supported by the Youth Found of National Natural Science Foundation of China (82103309).

Acknowledgments

The authors hereby express their gratitude to all participants who supported the study.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2022.905212/full#supplementary-material

Supplementary Figure 1 | External validation of FRS (A) KM survival curves for patients in the high- and low-FRS groups in the GEO cohort; (B) 1-, 3-, 5-, and 8-year ROC curves for FRS in the GEO cohort; (C) Survival status and FRS of patients in the GEO cohort; (D) Box plots showing the predicted IC50 values of the five most commonly used drugs in high- and low- FRS groups in the GEO cohort; (E) Immunotherapeutic responses of patients in the high- and low-FRS groups in the GEO cohort predicted using the TIDE algorithm; (F) Sensitivity of the patients in the high- and low-FRS groups to PD1 and CTLA4 treatment regimens in the GEO cohort predicted using the subclass mapping algorithm.

Supplementary Figure 2 | Subgroup Cox analysis of FRS Subgroup Cox regression analysis of FRS in TCGA (A) and GEO (B) cohorts.

Supplementary Figure 3 | Prediction of FRS-related small molecule compounds.

References

1. Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer statistics, 2021. CA: Cancer J Clin (2021) 71(1):7–33. doi: 10.3322/caac.21654

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Ortega-Gómez A, Rangel-Escareño C, Molina-Romero C, Macedo-Pérez EO, Avilés-Salas A, Lara-García A, et al. Gene-expression profiles in lung adenocarcinomas related to chronic wood smoke or tobacco exposure. Respir Res (2016) 17:42. doi: 10.1186/s12931-016-0346-3

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Song MA, Benowitz NL, Berman M, Brasky TM, Cummings KM, Hatsukami DK, et al. Cigarette filter ventilation and its relationship to increasing rates of lung adenocarcinoma. J Natl Cancer Instit (2017) 109(12):djx075. doi: 10.1093/jnci/djx075

CrossRef Full Text | Google Scholar

4. Vaz M, Hwang SY, Kagiampakis I, Phallen J, Patil A, O'Hagan HM, et al. Chronic cigarette smoke-induced epigenomic changes precede sensitization of bronchial epithelial cells to single-step transformation by KRAS mutations. Cancer Cell (2017) 32(3):360–76.e6. doi: 10.1016/j.ccell.2017.08.006

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Herbst RS, Morgensztern D, Boshoff C. The biology and management of non-small cell lung cancer. Nature (2018) 553(7689):446–54. doi: 10.1038/nature25183

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Binnewies M, Roberts EW, Kersten K, Chan V, Fearon DF, Merad M, et al. Understanding the tumor immune microenvironment (TIME) for effective therapy. Nat Med (2018) 24(5):541–50. doi: 10.1038/s41591-018-0014-x

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Gentles AJ, Bratman SV, Lee LJ, Harris JP, Feng W, Nair RV, et al. Integrating tumor and stromal gene expression signatures with clinical indices for survival stratification of early-stage non-small cell lung cancer. J Natl Cancer Instit (2015) 107(10):djv211. doi: 10.1093/jnci/djv211

CrossRef Full Text | Google Scholar

8. Mariathasan S, Turley SJ, Nickles D, Castiglioni A, Yuen K, Wang Y, et al. TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature (2018) 554(7693):544–8. doi: 10.1038/nature25501

PubMed Abstract | CrossRef Full Text | Google Scholar

9. 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(8):1277–89. doi: 10.1038/s41591-018-0096-5

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Auciello FR, Bulusu V, Oon C, Tait-Mulder J, Berry M, Bhattacharyya S, et al. A stromal lysolipid-autotaxin signaling axis promotes pancreatic tumor progression. Cancer Discover (2019) 9(5):617–27. doi: 10.1158/2159-8290.Cd-18-1212

CrossRef Full Text | Google Scholar

11. Sahai E, Astsaturov I, Cukierman E, DeNardo DG, Egeblad M, Evans RM, et al. A framework for advancing our understanding of cancer-associated fibroblasts. Nat Rev Canc (2020) 20(3):174–86. doi: 10.1038/s41568-019-0238-1

CrossRef Full Text | Google Scholar

12. Boyle ST, Poltavets V, Kular J, Pyne NT, Sandow JJ, Lewis AC, et al. ROCK-mediated selective activation of PERK signalling causes fibroblast reprogramming and tumour progression through a CRELD2-dependent mechanism. Nat Cell Biol (2020) 22(7):882–95. doi: 10.1038/s41556-020-0523-y

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Barrett RL, Puré E. Cancer-associated fibroblasts and their influence on tumor immunity and immunotherapy. eLife (2020) 9:e57243. doi: 10.7554/eLife.57243

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Mao X, Xu J, Wang W, Liang C, Hua J, Liu J, et al. Crosstalk between cancer-associated fibroblasts and immune cells in the tumor microenvironment: new findings and future perspectives. Mol Canc (2021) 20(1):131. doi: 10.1186/s12943-021-01428-1

CrossRef Full Text | Google Scholar

15. Erez N, Truitt M, Olson P, Arron ST, Hanahan D. Cancer-associated fibroblasts are activated in incipient neoplasia to orchestrate tumor-promoting inflammation in an NF-kappaB-Dependent manner. Cancer Cell (2010) 17(2):135–47. doi: 10.1016/j.ccr.2009.12.041

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Kim N, Kim HK, Lee K, Hong Y, Cho JH, Choi JW, et al. Single-cell RNA sequencing demonstrates the molecular and cellular reprogramming of metastatic lung adenocarcinoma. Nat Commun (2020) 11(1):2285. doi: 10.1038/s41467-020-16164-1

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Leek JT, Johnson WE, Parker HS, Jaffe AE, Storey JD. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinf (Oxford England) (2012) 28(6):882–3. doi: 10.1093/bioinformatics/bts034

CrossRef Full Text | Google Scholar

18. Qiu X, Mao Q, Tang Y, Wang L, Chawla R, Pliner HA, et al. Reversed graph embedding resolves complex single-cell trajectories. Nat Methods (2017) 14(10):979–82. doi: 10.1038/nmeth.4402

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Efremova M, Vento-Tormo M, Teichmann SA, Vento-Tormo R. CellPhoneDB: inferring cell-cell communication from combined expression of multi-subunit ligand-receptor complexes. Nat Protoc (2020) 15(4):1484–506. doi: 10.1038/s41596-020-0292-x

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Schröder MS, Culhane AC, Quackenbush J, Haibe-Kains B. Survcomp: an R/Bioconductor package for performance assessment and comparison of survival models. Bioinf (Oxford England) (2011) 27(22):3206–8. doi: 10.1093/bioinformatics/btr511

CrossRef Full Text | Google Scholar

21. McDermott DF, Huseni MA, Atkins MB, Motzer RJ, Rini BI, Escudier B, et al. Clinical activity and molecular correlates of response to atezolizumab alone or in combination with bevacizumab versus sunitinib in renal cell carcinoma. Nat Med (2018) 24(6):749–57. doi: 10.1038/s41591-018-0053-3

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Gibbons DL, Creighton CJ. Pan-cancer survey of epithelial-mesenchymal transition markers across the cancer genome atlas. Dev Dyn (2018) 247(3):555–64. doi: 10.1002/dvdy.24485

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Ayers M, Lunceford J, Nebozhyn M, Murphy E, Loboda A, Kaufman DR, et al. IFN-γ-related mRNA profile predicts clinical response to PD-1 blockade. J Clin Invest (2017) 127(8):2930–40. doi: 10.1172/jci91190

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Liang JY, Wang DS, Lin HC, Chen XX, Yang H, Zheng Y, et al. A novel ferroptosis-related gene signature for overall survival prediction in patients with hepatocellular carcinoma. Int J Biol Sci (2020) 16(13):2430–41. doi: 10.7150/ijbs.45050

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Liberzon A, Subramanian A, Pinchback R, Thorvaldsdóttir H, Tamayo P, Mesirov JP. Molecular signatures database (MSigDB) 3.0. Bioinf (Oxford England) (2011) 27(12):1739–40. doi: 10.1093/bioinformatics/btr260

CrossRef Full Text | Google Scholar

26. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods (2015) 12(5):453–7. doi: 10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Yoshihara K, Shahmoradgoli M, Martínez 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. doi: 10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, et al. Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade. Cell Rep (2017) 18(1):248–62. doi: 10.1016/j.celrep.2016.12.019

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Thorsson V, Gibbs DL, Brown SD, Wolf D, Bortone DS, Ou Yang TH, et al. The immune landscape of cancer. Immunity (2018) 48(4):812–30.e14. doi: 10.1016/j.immuni.2018.03.023

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Mayakonda A, Lin DC, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res (2018) 28(11):1747–56. doi: 10.1101/gr.239244.118

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Geeleher P, Cox N, Huang RS. pRRophetic: an r package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PloS One (2014) 9(9):e107468. doi: 10.1371/journal.pone.0107468

PubMed Abstract | CrossRef Full Text | Google Scholar

32. 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(10):1550–8. doi: 10.1038/s41591-018-0136-1

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Vickers AJ, Cronin AM, Elkin EB, Gonen M. Extensions to decision curve analysis, a novel method for evaluating diagnostic tests, prediction models and molecular markers. BMC Med Inf Decis Makin (2008) 8:53. doi: 10.1186/1472-6947-8-53

CrossRef Full Text | Google Scholar

34. van Wilpe S, Tolmeijer SH, Koornstra RHT, de Vries IJM, Gerritsen WR, Ligtenberg M, et al. Homologous recombination repair deficiency and implications for tumor immunogenicity. Cancers (2021) 13(9):13092249. doi: 10.3390/cancers13092249

CrossRef Full Text | Google Scholar

35. Desrichard A, Snyder A, Chan TA. Cancer neoantigens and applications for immunotherapy. Clin Cancer Res (2016) 22(4):807–12. doi: 10.1158/1078-0432.Ccr-14-3175

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Hu Z, Ott PA, Wu CJ. Towards personalized, tumour-specific, therapeutic vaccines for cancer. Nat Rev Immunol (2018) 18(3):168–82. doi: 10.1038/nri.2017.131

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Chan TA, Yarchoan M, Jaffee E, Swanton C, Quezada SA, Stenzinger A, et al. Development of tumor mutation burden as an immunotherapy biomarker: utility for the oncology clinic. Ann Oncol (2019) 30(1):44–56. doi: 10.1093/annonc/mdy495

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Matsushita H, Vesely MD, Koboldt DC, Rickert CG, Uppaluri R, Magrini VJ, et al. Cancer exome analysis reveals a T-cell-dependent mechanism of cancer immunoediting. Nature (2012) 482(7385):400–4. doi: 10.1038/nature10755

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Rizvi NA, Hellmann MD, Snyder A, Kvistborg P, Makarov V, Havel JJ, et al. Cancer immunology. mutational landscape determines sensitivity to PD-1 blockade in non-small cell lung cancer. Sci (New York NY) (2015) 348(6230):124–8. doi: 10.1126/science.aaa1348

CrossRef Full Text | Google Scholar

40. Gagliano T, Shah K, Gargani S, Lao L, Alsaleem M, Chen J, et al. PIK3Cδ expression by fibroblasts promotes triple-negative breast cancer progression. J Clin Invest (2020) 130(6):3188–204. doi: 10.1172/jci128313

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Li X, Shao C, Shi Y, Han W. Lessons learned from the blockade of immune checkpoints in cancer immunotherapy. J Hematol Oncol (2018) 11(1):31. doi: 10.1186/s13045-018-0578-4

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Forde PM, Chaft JE, Smith KN, Anagnostou V, Cottrell TR, Hellmann MD, et al. Neoadjuvant PD-1 blockade in resectable lung cancer. New Engl J Med (2018) 378(21):1976–86. doi: 10.1056/NEJMoa1716078

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Riera-Domingo C, Audigé A, Granja S, Cheng WC, Ho PC, Baltazar F, et al. Immunity, hypoxia, and metabolism-the ménage à trois of cancer: Implications for immunotherapy. Physiol Rev (2020) 100(1):1–102. doi: 10.1152/physrev.00018.2019

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Sharma RA, Harris AL, Dalgleish AG, Steward WP, O'Byrne KJ. Angiogenesis as a biomarker and target in cancer chemoprevention. Lancet Oncol (2001) 2(12):726–32. doi: 10.1016/s1470-2045(01)00586-1

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Motz GT, Coukos G. The parallel lives of angiogenesis and immunosuppression: cancer and other tales. Nat Rev Immunol (2011) 11(10):702–11. doi: 10.1038/nri3064

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Welti J, Loges S, Dimmeler S, Carmeliet P. Recent molecular discoveries in angiogenesis and antiangiogenic therapies in cancer. J Clin Invest (2013) 123(8):3190–200. doi: 10.1172/jci70212

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Gabrilovich DI, Nagaraj S. Myeloid-derived suppressor cells as regulators of the immune system. Nat Rev Immunol (2009) 9(3):162–74. doi: 10.1038/nri2506

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Zhang Y, Zhang Z. The history and advances in cancer immunotherapy: understanding the characteristics of tumor-infiltrating immune cells and their therapeutic implications. Cell Mol Immunol (2020) 17(8):807–21. doi: 10.1038/s41423-020-0488-6

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Duan Q, Zhang H, Zheng J, Zhang L. Turning cold into hot: Firing up the tumor microenvironment. Trends Canc (2020) 6(7):605–18. doi: 10.1016/j.trecan.2020.02.022

CrossRef Full Text | Google Scholar

50. Fu T, Dai LJ, Wu SY, Xiao Y, Ma D, Jiang YZ, et al. Spatial architecture of the immune microenvironment orchestrates tumor immunity and therapeutic response. J Hematol Oncol (2021) 14(1):98. doi: 10.1186/s13045-021-01103-4

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Lim SA, Wei J, Nguyen TM, Shi H, Su W, Palacios G, et al. Lipid signalling enforces functional specialization of t(reg) cells in tumours. Nature (2021) 591(7849):306–11. doi: 10.1038/s41586-021-03235-6

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Hamanishi J, Mandai M, Iwasaki M, Okazaki T, Tanaka Y, Yamaguchi K, et al. Programmed cell death 1 ligand 1 and tumor-infiltrating CD8+ T lymphocytes are prognostic factors of human ovarian cancer. Proc Natl Acad Sci USA (2007) 104(9):3360–5. doi: 10.1073/pnas.0611533104

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Shi Z, Zhao Q, Lv B, Qu X, Han X, Wang H, et al. Identification of biomarkers complementary to homologous recombination deficiency for improving the clinical outcome of ovarian serous cystadenocarcinoma. Clin Trans Med (2021) 11(5):e399. doi: 10.1002/ctm2.399

CrossRef Full Text | Google Scholar

54. Ganesh K, Stadler ZK, Cercek A, Mendelsohn RB, Shia J, Segal NH, et al. Immunotherapy in colorectal cancer: rationale, challenges and potential. Nat Rev Gastroenterol Hepatol (2019) 16(6):361–75. doi: 10.1038/s41575-019-0126-x

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Eso Y, Shimizu T, Takeda H, Takai A, Marusawa H. Microsatellite instability and immune checkpoint inhibitors: toward precision medicine against gastrointestinal and hepatobiliary cancers. J Gastroenterol (2020) 55(1):15–26. doi: 10.1007/s00535-019-01620-7

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Anagnostou V, Smith KN, Forde PM, Niknafs N, Bhattacharya R, White J, et al. Evolution of neoantigen landscape during immune checkpoint blockade in non-small cell lung cancer. Cancer Discov (2017) 7(3):264–76. doi: 10.1158/2159-8290.Cd-16-0828

PubMed Abstract | CrossRef Full Text | Google Scholar