Abstract

(1) Background

   In lung cancer, the use of small-molecule inhibitors, chemotherapy and
   immunotherapy has led to unprecedented survival benefits in selected
   patients. Considering most patients will experience a relapse within a
   short period of time due to single drug resistance, combination therapy
   is also particularly important to improve patient prognosis. Therefore,
   more robust biomarkers to predict responses to immunotherapy, targeted
   therapy, chemotherapy and rationally drug combination therapies may be
   helpful in clinical treatment choices.

(2) Methods

   We defined tumor-specific T cells (TSTs) and their features (TSTGs) by
   single-cell RNA sequencing. We applied LASSO regression to filter out
   the most survival-relevant TSTGs to form the Tumor-specific T cell
   score (TSTS). Immunological characteristics, enriched pathways, and
   mutation were evaluated in high- and low TSTS groups.

(3) Results

   We identified six clusters of T cells as TSTs in lung cancer, and four
   most robust genes from 9 feature genes expressed only on tumor-specific
   T cells were screened to construct a tumor-specific T cells score
   (TSTS). TSTS was positively correlated with immune infiltration and
   angiogenesis and negatively correlated with malignant cell
   proliferation. Moreover, potential vascular-immune crosstalk in lung
   cancer provides the theoretical basis for combined anti-angiogenic and
   immunotherapy. Noticeable, patients in high TSTS had better response to
   ICB and targeted therapy and patients in the low TSTS group often
   benefit from chemotherapy.

(4) Conclusion

   The proposed TSTS is a promising indicator to predict immunotherapy,
   targeted therapy and chemotherapy responses in lung cancer patients for
   helping clinical treatment choices.

   Keywords: Lung cancer, Tumor-specific T cell, Combination therapies,
   Immune infiltration, Immunotherapy

1. Introduction

   Lung cancer comprises the majority, around 25% of total cancer deaths
   worldwide and the 5-years overall survival (OS) has achieved very
   little progress in recent years [[37]1]. Despite advances in our
   understanding of new pathogenesis and treatment options, unfortunately,
   lung cancer is still one of the highly aggressive and rapidly fatal
   tumor subtypes [[38]2].New insights into the causal roles in lung
   cancer spurred the latest belief that tumorigenesis is entirely driven
   not only by cancer-cell-intrinsic genetic aberrations, but also is
   affected by the tumor microenvironment (TME) [[39]3]. Lung cancer
   develops in complex TMEs consisting of endothelial cells,
   tumor-infiltrating immune cells (TIICs) and fibroblast cells [[40]4].

   Immune checkpoint inhibitors (ICIs) which enhance immunization and
   relieve immune suppression and anti-angiogenesis therapy which inhibits
   tumor growth through inhibiting angiogenesis both had a significant
   effect on lung cancer. In ICI, a number of patients do not benefit from
   the initial therapy (primary resistance) and some responders recur
   after a period of response (acquired resistance) [[41]5]. The
   suppressive tumor immune microenvironment (TIME) is the main barrier to
   ICI efficacy [[42]6]. T cells with exhausted or regulatory phenotype
   are correlated with failure in anti-cancer immunity [[43]7]. Tumor
   cell-extrinsic mechanisms that induce ICI resistance involve components
   within the TME, including regulatory T cells (Tregs), myeloid-derived
   suppressor cells (MDSCs), and other inhibitory immune cells [[44]8].In
   anti-angiogenesis therapy, the majority of agents have shown some
   evidence of activity, but increased toxicities have been observed,
   including an increased risk of death for some agents, limiting their
   development [[45]9]. It is important to further understand the
   influence of immune cells in the TME on the proliferation of vascular
   endothelial cells. Pre-existing tumor-infiltrating lymphocyte (TIL) is
   the necessary precondition for potent tumor regression. Based on the
   status of pre-existing TIL, tumor microenvironments are classified into
   three types, in excluded infiltration type, where abnormal angiogenesis
   and immunosuppressive reactive stroma prevent the infiltration of T
   cells [[46]10]. Neo-vasculatures upregulate some inhibitory signals for
   an anti-tumor immune response such as PD-L1, indoleamine 2,
   3-dioxygenase (IDO), interleukin-6 (IL-6), and interleukin-10 (IL-10)
   [[47]11]. The interaction between immunity and angiogenesis renders
   tumor immune escape and treatment resistance. Many clinical studies
   have been conducted to investigate the synergistic effect of ICI plus
   anti-angiogenesis therapy in patients and have achieved good efficacy
   [[48]12].

   Recent advances in single-cell RNA sequencing (scRNA-seq) analysis have
   shed new light on the heterogeneity of T cell clusters in TME at
   single-cell resolution. Nevertheless, the tumor-specific T cells (TSTs)
   and how they affect the formation of the immune microenvironment,
   remains however poorly elucidated. How endothelial cells recruit TSTs
   to communicate is also unknown in lung cancer. The microbial abundance
   between cold and hot tumors has been even less discussed in lung
   cancer. Furthermore, combination strategies were identified in melanoma
   to overcome immunotherapy resistance [[49]13], thus, further studies
   are required to identify predictive biomarkers for combination
   approaches in lung cancer.

   In the present study, we aimed to uncover the characterization of the
   tumor-specific microenvironment by single-cell technology to guide the
   selection of therapeutic strategies. First, we identify TSTs with
   tumor-specific genes (TSTGs) for the tumor environment in lung cancer.
   Additionally, we applied Least absolute shrinkage and selection
   operator (LASSO) regression to filter out the most survival-relevant
   TSTGs to form the Tumor-Specific T Cell Score (TSTS). And we also
   defined these TSTs can receive the inhibitory signal from tip cells, a
   key phenotype in sprouting angiogenesis through specific crosstalk
   pairs. Finally, the TSTS was confirmed that can precisely predict the
   efficacy of combination strategies, immunotherapy, chemotherapy, and
   targeted therapy respectively in the 12 clinical treatment cohorts.

   Thus, these findings suggest that selecting treatment strategies based
   on the characteristics of the tumor microenvironment is a more accurate
   and reliable approach and that TSTS can as an indicator to measure the
   choice of clinical treatment in lung cancer.

2. Results

2.1. ScRNA-seq identified tumor-specific T cell clusters and tumor-specific T
cell features

   To investigate the heterogeneity of tumor-infiltrated immune cells, we
   obtained scRNA-Seq data from patients diagnosed with lung cancer. After
   normalization and principal component analysis, we applied graph-based
   clustering ([50]Fig. S1A, [51]Table S1). These clusters could be
   assigned to eight known cell lineages through marker genes ([52]Fig.
   S1B). The specific multicell group, characterized by high expression of
   MKI67, is a group of proliferating cells with active mitosis ([53]Fig.
   S1C). We divided T cells into CD4^+ and CD8^+ T cells by CD4, CD8A and
   CD8B expression ([54]Fig. S1D). Nine clusters of CD4^+ T cells and
   eleven clusters of CD8^+ T cells were defined ([55]Fig. 1A–D). The
   expression of signature genes and known functional markers suggested
   the presence of conventional naïve (CD4_CCR7, CD4_IL7R, CD8_CCR7),
   effector (CD4_GZMA, CD4_FGFBP2, CD8_GZMK and CD8_TYROBP), memory
   (CD4_FOS, CD8_TXNIP, CD8_ZNF683), exhausted clusters (CD4_CXCL13,
   CD8_CXCL13, CD8_TIGIT), CD4^+ Tregs (CD4_FOXP3), proliferative clusters
   (CD4_STMN1, CD8_STMN1) ([56]Fig. 1B–E [57]Table S2).To determine
   tumor-specific T cells, we compared the observed and expected cell
   numbers of each cluster (Ro/e) ([58]Fig. 1C–F) and the density of T
   cells between tumor and normal ([59]Fig. S1E). Combining the above two
   conditions, the exhausted clusters (CD4_CXCL13 and CD8_CXCL13), ISG
   gene-related subcluster (CD4_ISG15, CD8_ISG15), regulatory clusters
   (CD4_FOXP3), and pre-exhausted (CD8_ZNF683) enriched in the tumor,
   which can be defined as tumor-specific T cells (TSTs) which highly
   expressed immune checkpoint inhibitors like PD1, CTLA4 and LAG3
   ([60]Fig. S1F). Through pseudotime analysis, we determined that CD4^+
   exhausted T cells were transformed from effector T cells and
   CD8 + exhausted cells were transformed from effector memory T cells
   ([61]Fig. 1G and H). These T cells often lose their effector function
   and infiltrate in the microenvironment in tumor in the late stage of T
   cell transformation. We found the percentage of four representative
   TSTs was positively with pre-exhausted state cells which have effector
   function, which demonstrated that the more cells in the effector state,
   the more cells will transition to the exhausted state ([62]Fig. S1G).
   Cells in the ISG gene related subcluster (CD4_ISG15, CD8_ISG15) had
   highly expressed ISG-related genes like ISG15, IFITM1 and IFITM3, but
   we found that these genes were expressed not only in T cells, but in
   myeloid cells can not represent TSTs ([63]Fig. S1H). Using MCP_Counter
   deconvolution analysis with highly expressed genes in each cluster, we
   determined TSTs clusters-CD4_CXCL13, CD8_CXCL13, CD4_FOXP3 and
   CD8_ZNF683 had higher infiltration in tumor samples from TCGA compared
   with normal samples from GTEX ([64]Fig. 1I; [65]Table S3). Compared to
   normal tissues, nine representative genes -CXCL13, ZNF683, TOX2, TOX,
   LAIR2, CXCR6, LAYN, ACP5 and MAGEH1 were specially expressed in the
   tumor tissue in TSTs clusters ([66]Fig. 1J, [67]Fig. S1I) and these
   genes exclusively expressed in the TSTs clusters-CD4_CXCL13, CD4_FOXP3,
   CD8_CXCL13, CD8_ZNF683-were recognized as tumor-specific T cell genes
   (TSTGs) ([68]Fig. 1K). These genes were expressed only on T cells and
   not on other cells in the TME ([69]Fig. S1H). It is recognized that
   LAYN is closely related to tumorigenesis, and ZNF683 is chiefly
   implicated in the differentiation of T cells [[70]14]. Both TOX2 and
   TOX are transcription factors, which correlate with T cell exhaustion
   and tumor immunosuppressive microenvironment [[71]15].

Fig. 1.

   [72]Fig. 1
   [73]Open in a new tab

   ScRNAseq identified tumor-specific T cell clusters and tumor-specific
   genes in lung cancer (A) UMAP projection of CD4^+ T cells. Cells were
   color-coded for distinct clusters. (B) Dotplot exhibited the expression
   of marker genes (Top three) in nine subclusters of CD4^+ T cells. (C)
   Tissue prevalence estimated by Ro/e score in CD4^+ T cells. Cells
   derived from tumor, normal tissue, normal lymph node (NLN), metastases
   in lymph nodes (MLN), brain metastases (MBrain), and peritoneum (PE).
   (D) UMAP projection of CD8^+ T cells. (E) Dotplot exhibited the
   expression of marker genes in eleven subclusters of CD8^+ T cells. (F)
   Tissue prevalence estimated by Ro/e score in CD8^+ T cells. (G) UMAP
   visualization in the trajectory of defined CD4 T cell types marked by
   colors, the dark line denoted cellular evolutionary directions (left).
   Visualization of single-lineage cells marked with inferred pseudotime
   by Monocle3 (using UMAP). The black curve denotes the inferred lineage
   (right). (H) UMAP visualization in the trajectory of defined CD8 T cell
   types marked by colors, the dark line denoted cellular evolutionary
   directions (left). Visualization of single-lineage cells marked with
   inferred pseudotime by Monocle3 (using UMAP). The black curve denotes
   the inferred lineage (right). (I) Boxplots representing the MCP_counter
   score of each subclusters of CD4 T cells (upper) and CD8 T cells
   (lower) in tumor versus normal. (J) Dotplot denotes the expression of
   each gene (row) in tumor versus normal. (K) Dotplot indicating the
   expression of tumor-specific T cell genes in each cluster of T cells.

2.2. Construction of tumor-specific T cell gene signatures

   To avoid heterogeneity bias from single-cell data, we validated the
   correlation between TSTGs and immune-related genes using bulk RNA-seq
   in TCGA. The result demonstrated that the TSTGs (CXCL13, ZNF683, TOX2,
   TOX, LAIR2, CXCR6, LAYN, ACP5 and MAGEH1) and the exhausted T cells
   associated genes (PDCD1, CTLA4, LAG3, HAVCR2, TIGIT), the effector T
   cells associated genes (GZMA, GZMB), the Tregs-related genes (FOXP3,
   IL2RA) were positively correlated ([74]Fig. 2A). Subsequently,
   single-sample Gene Set Enrichment Analysis (ssGSEA) was applied in TCGA
   to infer from bulk-RNA seq to reveal the relative level of immune cell
   infiltration. Positive correlation between gene expression of TSTGs and
   immune cell infiltration score was proved that TSTGs can reflect the
   state of immune infiltration in bulk RNA-seq ([75]Fig. 2B–[76]Tables
   S4–5). The immune microenvironment is often related to over survival.
   The expression of these genes was further assessed using univariate Cox
   regression analysis. The results indicated that 4 of 9 candidate genes
   were significantly correlated with OS of lung cancer patients ([77]Fig.
   2C, p < 0.05). In order to identify the optimal prognostic biomarker,
   the least absolute shrinkage and selection operator (LASSO) logistic
   regression was applied to screen the tumor-specific features and
   recruited to construct a multi-gene prognostic signature ([78]Figs. S2A
   and B). The four genes including LAYN, TOX, MAGEH1 and CXCR6 are
   considered to be the optimal biomarkers with the best features that
   build tumor-specific T cell score (TSTS) ([79]Fig. 2D). In the TCGA
   training set, patients were divided into two groups according to the
   TSTS score. The patients with lower TSTS score had a worse overall
   prognosis compared to the patients with higher TSTS score ([80]Fig.
   2E), which was validated in five independent validation cohorts
   ([81]Fig. 2F). Further survival analysis demonstrated that the high
   TSTS score group showed a decreased disease-free survival than the low
   TSTS score group ([82]Fig. S2C). To test the independence of our
   scoring system, a multivariate Cox regression analysis was then
   performed. After adjusting by common clinical parameters such as age,
   gender, stage and smoking status, the TSTS score still displayed robust
   assessment ability ([83]Fig. S2D). In clinical staging, we also noted
   that the TSTS score represents early-clinical stages in TCGA dataset
   ([84]Fig. S2E). Early-stage lung cancer patients had higher TSTS score,
   suggesting that there may be higher tumor-specific T cells infiltrating
   in the TME at this time. To further investigate the clinical
   potentiality of the TSTS score of lung cancer, stratified analysis
   based on these clinical characteristics was implemented. The results
   indicated that the TSTS score seemed to predict OS of lung cancer
   patients in subgroup of T1, T2, T3, T4, N0, N1, N2, N3, M0 and M1
   ([85]Fig. S2F). Patients in low groups had significantly poorer
   clinical outcomes than those in high groups. The prognostic value of
   the TSTS score for the survival outcome of patients was verified in
   pan-cancer types using TCGA datasets, TSTS may be used as a biomarker
   of pan-cancer to predict patients’ survival ([86]Fig. S3A).

Fig. 2.

   [87]Fig. 2
   [88]Open in a new tab

   Construction of tumor specific T cell score with survival outcomes (A)
   The correlation between the expression of TSTGs (CXCL13, ZNF683, TOX2,
   TOX, LAIR2, CXCR6, LAYN, ACP5 and MAGEH1) and the expression of immune
   checkpoint genes (PDCD1, CTLA4, LAG3, HAVCR2, TIGIT and BTLA), the
   effector T cells associated genes (GZMA, GZMB and PRF1) and the
   Tregs-related genes (FOXP3 and IL2RA). P-value <0.001, Pearson
   correlation analysis for the datasets that meet the normal
   distribution, Spearman correlation analysis for those that are not
   normally distributed (detailed statistical information of these dataset
   is attached to [89]Table S11). (B) The correlation between the
   expression of TSTGs and ssGSEA score which represents the infiltration
   of immune cells. (C) The p-value and HR of selected genes in
   univariable Cox regression analysis. (D) The bar plot shows the lasso
   cox coefficient of four optimal genes (LAYN, TOX, MAGEH1, and CXCR6)
   associated with prognosis. (E) Kaplan-Meier (K–M) curves demonstrated
   low TSTS group (orange) had noticeably poorer over survival than
   patients in the high TSTS group (blue). (F) Kaplan-Meier curve of
   validation cohorts ([90]GSE68465, [91]GSE14814, [92]GSE30219,
   [93]GSE31210, [94]GSE37745) showing low TSTS group (orange) had
   noticeably poorer survival than patients in the high TSTS group.
   P-value was calculated by the log-rank test.

2.3. Immune infiltration was enhanced in the high TSTS group

   Furthermore, in the lung cancer samples from TCGA, firstly the high and
   low groups are divided at the median, then it is intriguing to
   investigate that the high TSTS group was featured by high immune
   infiltration, which enriched genes were immune-related genes such as
   CXCL9, IFNG, IL2, MS4A1, CLEC6A, and KLRC4 ([95]Fig. 3A, [96]Table S6).
   The elevated gene of the high TSTS group was enriched in immunological
   enhancement pathways, including immunoglobulin production, positive
   regulation of lymphocyte activation, T cell proliferation, cell
   killing, and Type I interferon production ([97]Fig. 3B). Similarly, the
   Hallmarks pathway in the high TSTS group also revealed highly
   significant enrichment in immunological activity including
   interferon-alpha response, interferon-and gamma response, inflammatory
   response ([98]Fig. 3C–[99]Table S7). Five immunomodulators (chemokine,
   immunostimulator, MHC, and chemokine receptors, inhibitory immune
   checkpoint) and the infiltration levels of five types of TIICs (CD8^+ T
   cells, DC, NK cells, macrophages, and Th1 cells) were enriched in high
   TSTS ([100]Fig. 3D). MCP counter, XCELL, EPIC, quanTIseq and TIMER
   algorithms were used to calculate an immune score and estimate the
   abundance of various types of immune cells. We found higher proportion
   of CD4^+ T cells, CD8+T cells, NK cells and B cells in lung cancer with
   higher TSTS ([101]Fig. 3E-F, [102]Tables S8–9). Consistent with the
   above studies, TSTS closely related to immunopotentiation validated by
   multiple deconvolution algorithms ([103]Fig. S4A). Meanwhile, TSTS
   highly positively correlated with cytotoxic lymphocytes ([104]Fig.
   3G),. To our surprise, there were a large number of endothelial cells
   in the high TSTS group. ([105]Fig. 3F, [106]Fig. S4A). In addition, the
   higher number and diversity of T-cell receptors (TCRs) in the high TSTS
   group than that in the low TSTS group (P < 0.01) suggested higher tumor
   antigen stimulation and TCR expansion occurs in the high TSTS group
   ([107]Fig. 3H). Combined cell compositions and pathway of the high TSTS
   group, more favorable immune environment consisting of immune cells and
   immune components exist in higher groups. Meanwhile, the most
   ICI-related genes used as therapy targets were the highly positively
   correlated with TSTS, a direct clue that TSTS may be suitable. for ICI
   therapy ([108]Fig. 3I).

Fig. 3.

   [109]Fig. 3
   [110]Open in a new tab

   Enrichment of immune infiltration in high TSTS group (A) Volcano plot
   demonstrating the DE genes between high TSTS group and low TSTS group.
   (B) Enriched Gene Ontology (GO) pathways of high TSTS group shows
   immunological enhancement in high TSTS group. Statistical analysis was
   performed using GSEA analysis. The adjust P value labeled to legend.
   (C) Enriched hallmark pathways of the high TSTS group. Statistical
   analysis was performed using GSEA analysis. The adjust P value labeled
   to legend. (D) Heatmap shows immune-related genes (chemokine,
   immunostimulator, MHC, chemokine receptors and inhibitory immune
   checkpoint) and canonical genes of five TIICs (CD8^+ T cells, DC, NK
   cells, macrophages, and Th1 cells) were highly expressed in high TSTS
   group. (E) By deconvolution method TIMER, boxplot demonstrated the six
   immune cells' score in the high TSTS versus the low TSTS group. (F) By
   deconvolution method MCP_counter analysis, boxplot demonstrated the ten
   types of cells' score in the high TSTS versus the low TSTS group. (G)
   The correlation between the TST scores and tertiary lymphoid structure
   (TLS) score. TLS score was calculated from the average expression of
   CCL19, CCL21, CXCL13, CCR7, SELL, LAMP3, CXCR4, CD86 and BCL6.
   Statistical significance was determined by Spearman's rank correlation
   analysis (p value < 0.001). (H) TCR clone number (left) and TCR clone
   diversity (right) was significantly elevated in high TSTS group. (I)
   Correlation between TSTS score and expression of known immune
   checkpoint genes. ns P > 0.05; *P < 0.05; ***P < 0.001; ****P < 0.0001.
   Statistical significance was determined by Spearman's rank correlation
   analysis (p value < 0.001).

2.4. Angiogenesis and tumor-specific T cell infiltration have synergistic
effects

   In consideration of the infiltration of endothelial cells in the high
   TSTS group, we evaluated how endothelial cells are involved in the
   formation of tumor microenvironment in the high TSTS group. We found
   TSTS score was positively correlated with MCP_counter score of
   endothelial cells ([111]Fig. 4A). We also found TSTS score had positive
   correlation with migration (NRP1), and vessel formation (FLT1, KDR,
   NRP1, FLT4, NRP2 and ENG), adhesion (CDH5 and PECAM1) related gene
   expression level ([112]Fig. 4B). By immunohistochemical analysis, we
   found high expression of TOX, LAYN, and MAGEH1 proteins that make up
   the TSTS fraction in the perivascular area. In contrast,
   immunohistochemical samples in which no obvious blood vessels could be
   found often showed negative protein staining ([113]Fig. 4C). With the
   utilization of the MCP_counter de-convolution algorithm, T cells and
   endothelial cells had higher infiltration level in the high TSTS group
   in five independent validation bulk datasets ([114]Fig. 4D). We used
   spatial transcriptome from a patient with lung cancer to demonstrate
   that the emergence of TSTS was often accompanied by vascular
   infiltration ([115]Fig. S5A). In order to better understand the
   relationship between immune infiltration and angiogenesis, we divided
   all endothelial cells into six clusters ([116]Fig. 4E; [117]Table S2).
   By integrating the function-related gene signatures in previous study
   [[118]16], we determined the arteries, capillaries, tips, and veins of
   the blood vasculature, as well as the lymphatic vessels ([119]Fig.
   S5B). Remarkably, cluster Tip_COL4A1, which enriched in tumor tissue
   the key phenotype in sprouting angiogenesis, expressed a signature of
   genes involved in angiogenesis and endothelial cell growth (FLT1, KDR,
   PGF, PIGF and SPARC), cell adhesion, migration (LAMC1, LAMA4, LAMB1 and
   SPARC), facilitating their navigating role in the process of vessel
   sprouting ([120]Fig. 4F and G). With the MCP_counter de-convolution
   algorithm, we found that six endothelial cell clusters had high
   infiltration level in the high TSTS group in TCGA ([121]Fig. 4H). These
   observations suggested that angiogenesis serves as an important
   characteristic in the high TSTS group, and the tip EC phenotype is
   likely to be a promising target of anti-angiogenic therapy (AAT),
   especially anti-vascular endothelial growth factor receptor (VEGFR)
   antibodies in lung cancer. To dissect the complex activities of tip
   cells in TME, we subsequently used CellPhoneDB to investigate the
   molecular communication between T cells and tip cells. In this
   analysis, we observed widespread interactions between Tip_COL4A1 and
   TSTs cells, which expressed the CD94:NKG2A receptors on surface and
   received the inhibitory signal HLA-E from Tip_COL4A1 ([122]Fig.
   4I–[123]Table S10). As a result, CD8^+ T cell activation is suppressed,
   it was proved a mechanism of immune escape in various tumor types such
   as glioma, melanoma, or colorectal cancer [[124]17], it may hint the
   presence of tip cells may contributed to immune escape in lung cancer.
   Taken together, these observations reflected the therapeutic potential
   for AAT alone or in combination with ICB in lung cancer.

Fig. 4.

   [125]Fig. 4
   [126]Open in a new tab

   Angiogenesis and tumor-specific T cell infiltration have synergistic
   effects (A) Positive correlation between the TSTS and the endothelial
   cells' score (calculated by MCP-counter). Statistical significance was
   determined by Spearman's rank correlation analysis (p value < 0.001).
   (B) Positive correlation between TSTS and the expression of vascular
   related genes. P-value <0.001, Pearson correlation analysis for the
   datasets that meet the normal distribution, Spearman correlation
   analysis for those that are not normally distributed (detailed
   statistical information of these dataset is attached to [127]Table
   S11). (C) Immunostaining of representative genes (TOX, LAYN, and
   MAGEH1). Orange arrows and dashed boxes indicate blood vessels. (D)
   Utilization of MCP_counter de-convolution algorithm to compare the
   infiltration of endothelial cells between TSTS low and high groups in
   five independent cohorts. (E) UMAP analysis identified six clusters
   that named arteries, veins, capillaries and tips of the blood
   vasculature, as well as the lymphatic vessels. (F) Density plot of UMAP
   representation comparing normal and tumor samples. Solid line
   highlights the different infiltration level between normal and tumor.
   (G) Dotplot shows the expression of angiogenesis and endothelial cell
   growth genes in Tip_COL4A1 cluster. (H) By MCP_counter analysis,
   infiltration score of six clusters endothelial cells presents different
   state in the high and low TSTS group. (I) Dotplot showing the
   statistically significant ligand-receptor interactions between
   Tip_COL4A1 cells and subclusters of T cells, inferred by CellPhoneDB in
   which cluster of tumor specific T cells highlighted by red. P values
   are indicated by circle size. ns P > 0.05; *P < 0.05; ***P < 0.001;
   ****P < 0.0001.

2.5. Elevated malignant cell proliferation status in the low TSTS group

   Of note, In single cell analysis, the proportion of epithelial cells
   was negatively correlated with tumor-specific T cells (TSTs) which also
   found in immune cell deconvolution of TCGA dataset ([128]Fig. S4A),
   suggesting that more 8epithelial cells were clustered in the low TSTS
   group ([129]Fig. 5A). Moreover, genes involved in proliferation (MK167
   and TUBB) and synthesis of important nucleotides necessary for cell
   growth (TYMS, RRM1and ATIC) had high expression in low TSTS groups in
   five validation cohorts ([130]Fig. 5B). And we found that the TSTs
   score was positively correlated with the expression of pyroptosis and
   ferroptosis sensitive genes and negatively correlated with the
   expression of ferroptosis resistance genes which suggests that there is
   an anti-apoptotic mechanism in the low TSTS group ([131]Fig. 5C).

Fig. 5.

   [132]Fig. 5
   [133]Open in a new tab

   The low TSTS group showed a state of rapid proliferation of malignant
   cells (A) The correlation analysis between the proportion of Tregs
   cells (CD4_FOXP3), exhausted T cell cells (CD4_CXCL13 and CD8_CXCL13),
   pre-exhausted cells (CD8_ZNF683) and the proportion of epithelial
   cells. Each dot represents a patient. In one patient, the ratio was
   calculated as the proportion of the number of cells in the appointed
   group to the total number of cells. Statistical significance was
   determined by Spearman's rank correlation analysis (p value < 0.001).
   (B) Boxplot shows the different expression level of
   proliferation-related genes in high TSTS and TSTS group in five
   independent cohorts. (C) In the TCGA dataset, the correlation between
   the TSTS and the expression of pyroptosis sensitive genes (left) and
   ferroptosis resistance or sensitive genes (right). Statistical
   significance was determined by Spearman's rank correlation analysis (p
   value < 0.001).

2.6. Associated oncogenic event reveals different TSTS states

   Subsequently, we probed into the deep-seated reasons at the genome
   level underlying the distinct features between these two groups. The
   comparison of the proportions of tumor mutation burden (TMB) indicated
   no discernible difference ([134]Fig. S6A, p > 0.05). Nevertheless, the
   mutant landscape of the top 10 most frequently mutated genes was not
   significantly different between the two groups ([135]Fig. S6B).
   Strikingly, KEAP1, DGKB, PSG8, and ITPR2, were identified as
   statistically different in odds ratios for low versus high groups
   ([136]Fig. S6C). The Low TSTS group displayed enriched gene mutation of
   KEAP1, which accelerate proliferation and oncogenic transformation
   [[137]18]. We found that there were the same or different mutational
   signatures in the high and low TSTS groups. Tobacco signature, APOBEC
   signature and defective DNA mismatch repair signature can explain most
   of the mutation signatures. The aging-associated signature SBS5 only
   appeared in the high TSTS group, and SBS40 only appeared in the low
   TSTS group. This may be an important reason for two kinds of
   microenvironment infiltration in lung cancer ([138]Figs. S6D–E). What's
   more, we do further analysis to investigate whether mutation and TSTs
   can provide a more effective prediction of tumor progression. The
   results indicate that among the common mutations (TP53, KRAS, EGFR and
   ALK) in lung adenocarcinoma, the survival outcomes predicted by
   environmental factors represented by TSTs and mutations are fairly
   consistent, with no statistically significant differences ([139]Fig.
   S6F). We found that when patients with STK11 mutation and lower TSTS
   score, the patient had a worse prognosis, indicating that STK11
   mutations could make this indicator more clinically meaningful.

2.7. TSTS was identified as predictor of immunotherapy and targeted therapy

   Taken together, high TSTS represents high immune environment with more
   exhausted T cells infiltration combining abundant tertiary lymphoid
   structure and angiogenesis while low TSTS exhibited state of high
   malignant cell proliferation, ([140]Fig. S5B). Further in the clinical
   treatment cohort, in the NSCLC immunotherapy cohort, the high TSTS
   group was significantly associated with better survival outcome,
   suggesting that TSTS can predict immunotherapy efficacy, and this
   ability was also validated in the immunotherapy cohort of the CCRCC,
   BLCA and SKCM ([141]Fig. 6A). The prognosis of the high TSTS group was
   better due to higher ICB response rates ([142]Fig. 6B). Furthermore, we
   examined a single-cell dataset ([143]GSE207422) of individuals
   undergoing anti-PD-1 therapy and assessed the levels of TSTS gene
   expression across groups with different treatment outcomes. The results
   revealed that TSTS genes were notably enriched in the CD4_CXCL13,
   CD4_FOXP3, and CD8_CXCL13 cell subclusters ([144]Figs. S7A–B).
   Importantly, the degree of enrichment within these three cell subgroups
   was significantly higher in the partial response (PR) group compared to
   the stable disease (SD) group. This provides compelling evidence that
   TSTS is indeed a relevant indicator associated with the efficacy of
   immunotherapy ([145]Fig. S7C).

Fig. 6.

   [146]Fig. 6
   [147]Open in a new tab

   TSTS identified as a predictor of immunotherapy, chemotherapy and
   targeted/anti-angiogenic therapy (A) Kaplan-Meier (K–M) curves
   demonstrated low TSTS group (red) had noticeably poorer survival than
   patients in high TSTS group (green) in immunotherapy treatment cohorts
   (B) The distribution of TSTS between responders (sTable. disease,
   partial response and complete response) and non-responders (progressive
   disease) in immunotherapy treatment cohorts. (C) Heatmap shows the
   different expression level of NSCLC-related drug-target genes between
   low TSTS and high TSTS (P < 0.05) (D) The distribution of TSTS between
   responders (sTable. disease, partial response and complete response)
   and non-responders (progressive disease) in anti-angiogenic treatment
   cohorts. (E) Kaplan-Meier (K–M) curves demonstrated low TSTS group
   (red) had noticeably poorer survival than patients in high TSTS group
   (green) in anti-angiogenic treatment cohorts. (F) Targeted-therapeutic
   sensitivity of six common drugs were estimated and compared in TCGA
   cohort. (G) Chemotherapeutic sensitivity of four common drugs were
   estimated and compared in TCGA cohort. (H) The distribution of TSTS
   between responders (sTable. disease, partial response and complete
   response) and non-responders (progressive disease) in TCGA chemotherapy
   treatment cohorts. ns P > 0.05; *P < 0.05; ***P < 0.001;
   ****P < 0.0001.

   Notably, comparing TSTS score to other well-recognized immunotherapy
   predictive signature [[148]19,[149]20] in the NSCLC immunotherapy
   cohort, Cytotoxic T cells (CTL), IFNG, T cell inflamed GEP score were
   statistically significant in predicting the efficacy of immunotherapy,
   furthermore, the Hazard Rates (HR) of TSTS was the smallest of all
   indicators (TSTS HR = 0.20, CTL HR = 0.45, IFNG HR = 0.36, T cell
   inflamed GEP = 0.37, CD8 HR = 0.40), implying that TSTS was the
   strongest protective factor for immunotherapy ([150]Fig. S8A).

   Patients with DGKB, PSG8, and ITPR2 mutations correlated with better
   survival than those with no mutation ([151]Fig. S8B). Nevertheless, the
   patients carrying KEAP1 mutation had a poor survival outcome ([152]Fig.
   S8C). Thus, the gene mutations mentioned above may serve as potential
   biomarkers and predictors of therapeutic efficacy in lung cancer
   patients.

   Additionally, we searched the Drugbank database for therapeutic drugs
   and drug targets that have currently been approved in non–small-cell
   lung cancer (NSCLC). Genes that are the targets of chemotherapy drugs
   (vinorelbine, Gemcitabine and Pemetrexed) were highly expressed in the
   low TSTS group. Genes that are the targets of targeted therapy drugs
   (ALK inhibitor, EGFR inhibitor and antiangiogenic inhibitor) were
   highly expressed in the high TSTS group ([153]Fig. 6C). Strikingly, the
   result of the NSCLC cohort provided unique insights into the better
   synergistic effect of Erlotinib (EGFR-TKI) combined with Bevacizumab
   (an antiangiogenic drug) in high TSTS group (p < 0.05). By analysis Two
   antiangiogenic therapy (Bevacizumab) cohorts, we also found patients
   with better responses had higher TSTS ([154]Fig. 6D). Obviously, in
   antiangiogenic therapy cohorts, patients with high TSTS may have a
   better prognosis outcome than those with low TSTS ([155]Fig. 6E). The
   estimated IC50 values of targeted therapy drugs were significantly
   elevated in TSTS-low samples of TCGA cohort ([156]Fig. 6F); however,
   the estimated IC50 values of chemotherapy drugs were significantly
   elevated in TSTS-high samples ([157]Fig. 6G). And in the chemotherapy
   cohort from the TCGA database, the low TSTS group had a better response
   efficacy ([158]Fig. 6H).

   Taken together, patients of advanced lung cancer with low TSTS were
   preferable for chemotherapy, whereas patients of early lung cancer with
   high TSTS fit targeted therapy (antiangiogenic inhibitor and small
   molecule inhibitors) and/or in combination with immunotherapy.

3. Discussion

   TILs especially T lymphocytes, which are the main anti-tumor cells.
   There is ample evidence to demonstrate that the presence of TILs is
   associated with enhanced host response to the tumor and therefore is
   associated with improved patient outcomes [[159]21].

   Previous studies revealed that TILs consist of not only those T cells
   specific for TSAs, but also the ones that recognize epitopes unrelated
   to the tumor, such as virus-specific antigens, known as bystander T
   cells [[160]22]. In contrast to bystander T cells, although exhibiting
   a dysfunctional phenotype, tumor-specific T cells possess significant
   capabilities of tumor-specific killing [[161]23]. Therefore, we tried
   to find and use tumor-specific T cells related index to measure the
   level of immune infiltration and to precisely predict the response to
   ICB. Meanwhile, the strategy of combining ICB with ramucirumab
   demonstrated favorable effects in clinical trials [[162][24],
   [163][25], [164][26]]. Indicators urgently needed to guide combination
   therapy offer hope for the treatment of advanced lung cancer.

   Through the single-cell sequencing technology, we first identified
   unique TSTs. Intriguingly, our findings with the gene expression
   profiles in TSTs further identified the tumor-specific T cell features,
   consisting of CXCL13, ZNF683, TOX2, TOX, LAIR2, CXCR6, LAYN, ACP5 and
   MAGEH1. LASSO Cox regression model was successively used to filter out
   the most robust markers and tumor-specific T cell-related prognostic
   score (TSTS) was established. Therefore, a TSTS model consisting of
   LAYN, TOX, MAGEH1 and CXCR6, classifies NSCLC patients into high TSTS
   group with better prognosis and low TSTS group with poor prognosis.
   Strikingly, TOX was reported as a critical regulator of tumor-specific
   T cell differentiation, and highly expressed in dysfunctional TST cells
   from tumors [[165]27]. Remarkably, LAYN has been corroborated to
   suppress CD8^+ T cells, and high expression of LAYN, MAGEH1 in Treg
   cells correlated with poor prognosis of NSCLC and colorectal cancer
   (CRC) [[166]28]. Also, we found that chemokine receptor CXCR6 was
   enriched in CD4^+ Treg and CD8^+ exhausted T cell clusters. Knockdown
   of CXCR6 inhibited cell proliferation migration and invasion of
   hepatocellular carcinoma (HCC) cells, followed by the deregulated
   expression of vascular endothelial growth factor (VEGF) [[167]29].
   CXCR6 is vital for effector-like activate cytotoxic lymphocytes (CTLs)
   in the TME to maximize their antitumor efficacy before progressing to
   irreversible dysfunction [[168]30].

   Higher tumor antigen stimulation and TCR expansion occurs in the high
   TSTS group than that in the low TSTS group suggested higher ability to
   recognize tumor-specific antigen that ultimatly kill tumor cells.
   Meanwhile the exhaustion of effector cells indicates that the more
   effector T cells become toxic before moved towards terminal
   differentiation. By using bulk RNA sequencing, we identified that
   immune-related pathways and cellular adhesion pathways were
   significantly elevated in the high TSTS group, immune-related genes,
   such as PD1, TIGIT and CTLA4 were highly expressed in the high TSTS
   group. The Low TSTS group exhibits low immune infiltration and high
   proliferative properties. Therefore, converse to the immunologically
   hot (T cells inflamed) TME in high TSTS group, cold (non-T cells
   inflamed) TME have been related with low TSTS, indicating that TSTS
   defines hot or cold TME. Ultimately, this study validated the ability
   of TSTS as a predictor of immunotherapy efficacy in pan-cancer in the
   SKCM, BTCC, CCRCC and NSCLC immunotherapy datasets and found strong
   predictive power in NSCLC, SKCM, CCRCC, and BTCC. We believe that this
   can be explained by that T cell subgroups are relatively stable. and do
   not tend to be heterogeneous across cancer types, resulting in the
   presence of tumor-specific T cells that are likely to be shared across
   cancer types. This phenomenon is also consistent with that reported by
   Zheng, L. et al. [[169]31].

   Strikingly, a strong positive correlation between TSTS and angiogenesis
   was identified in our present work. Sustained angiogenesis is
   considered to be an important hallmark of NSCLC [[170]32]. Angiogenesis
   contributes to NSCLC development by providing oxygen and nutrients and
   facilitating tumor growth and metastasis [[171]33]. We found that the
   tumor angiogenesis was predominantly featured by tip cell
   infiltration-the main cells at the tips of vascular sprouts integrating
   multiple activities during angiogenesis [[172]34,[173]35]. The
   sprouting angiogenesis promotes the creation of emerging blood vessels
   from pre-existing vessels, along with the tip endothelial cells, which
   may significantly affect immunity. Our research confirmed angiogenesis
   high enriched in high TSTS group. Subsequently, we further identified
   the comprehensive intercellular communication between tip cells and
   tumor-specific T cells. Notably, HLA-E/β2m-CD94/NKG2A, a new druggable
   inhibitory immune checkpoint, was found to be a highly communicative
   ligand pair between tip cells and CD8^+ exhausted T cells, which
   contributes to the functional exhaustion of cytotoxic NK cells and CD8
   T cells [[174]36,[175]37]. In addition to inhibiting the induction of
   local hypoxia, promoting the aggregation of MDSCs, Treg and leading to
   a suppressive tumor immune microenvironment, inhibition of
   neovascularization can reduce the expression of immunosuppressive check
   pairs and provide new evidence for synergistic effects of immunotherapy
   and anti-angiogenesis therapy. High TSTS group, enriched in tip cells,
   were proved better response in Bevacizumab (AAT)-therapy cohort. Thus,
   it is plausible to make recommendations on the combination of
   anti-angiogenic and immunotherapy for patients with high TSTS.

   Overall, low TSTS patients had cold TME, characterized by decreased
   immune infiltration, decreased cells apoptosis and increased cells
   proliferation. These patients represent a phenotype of “immune desert”,
   which has been associated with the poor outcome to ICB
   [[176]38,[177]39]. Surprisingly, not only were epithelial cells
   enriched in the low TSTS group but also the TSTS score is negatively
   correlated with the proliferative genes and Ferroptosis resistance
   genes and is positively correlated with pyroptosis sensitive genes
   ferroptosis sensitive genes, implying that proliferation and growth of
   cells in the low TSTS group would be more sensitive, meaning that
   patients in the low TSTS group tended to have a better response to
   chemotherapy. This is one of the unexpected findings of this study.

   Meanwhile, also, we found that TMB as an immunotherapy-guided biomarker
   did not seem to differ between high and low TSTS groups. This may
   explain why previous research has reported that TMB has a limited
   predictive effect on non-small cell lung cancer [[178]40]. It should be
   noted that TSTs, as an indicator of immune inflammation, is higher in
   patients with early-stage lung cancer, while patients with advanced
   lung cancer tend to be low TSTS with “immune desert”. This is in line
   with one previous clinical trial that earlier use of immunotherapy will
   bring more significant benefits to NSCLC patients [[179]41].

   The treatment of advanced NSCLC is mostly based on a combination of
   chemotherapy or targeted therapy. In recent years, immunotherapy
   targeting programmed cell death protein 1 (PD-1) and its ligand (PD-L1)
   or cytotoxic T-lymphocyte-associated antigen 4 (CTLA-4) has made
   significant progress in the treatment of NSCLC. Anti-CTLA-4 inhibitor
   has achieved breakthroughs in the treatment of NSCLC, demonstrating
   significant clinical efficacy, prolonging patient survival and
   dramatically changing the treatment landscape of NSCLC. Thus, such a
   biomarker, which can evaluate the efficacy of chemotherapy, targeted
   therapy and immunotherapy, brings a new breakthrough in the
   stratification of patients for comprehensive drug use.

   The limitation of our approach should also be taken into consideration.
   Although RNA sequencing, whole exon sequencing, and single-cell RNA
   sequencing were used and validated in as many independent cohorts as
   possible, all conclusions were based on bioinformatics analysis. Our
   results warrant further validation in a large number of real-world
   cohorts.

4. Methods and materials

4.1. Data acquisition

   Single-cell RNA-sequencing (ScRNA-seq) data (accession number:
   [180]GSE131907, [181]GSE207422) of lung adenocarcinoma samples from the
   initial publication was downloaded and reanalyzed for this manuscript.
   Spatial transcriptomic data for one patient with NSCLC were downloaded
   from 10X
   ([182]https://www.10xgenomics.com/resources/datasets/human-lung-cancer-
   11-mm-capture-area-ffpe-2-standard).The Cancer Genome Atlas (TCGA)
   data: For comparison within a single TCGA-LUAD cohort, lung
   adenocarcinoma (LUAD) RNA Sequencing data (TPM), somatic mutation data
   and survival information were downloaded from the UCSC Xena data portal
   ([183]https://xenabrowser.net/datapages/). Raw RNA-seq was annotated by
   Gencode V23 and quantificated by RSEM (v1.3.3), tpm values are
   extracted, log2 (tpm+0.001) transformed, and combined all sample. Gene
   Expression Omnibus (GEO): RNA sequencing data (number: [184]GSE68465,
   [185]GSE14814, [186]GSE30219, [187]GSE31210 and [188]GSE37745) of lung
   cancer with detailed survival data were downloaded as validation
   queues. RNA sequencing of immunotherapy cohort of non -small-cell lung
   cancer ([189]GSE135222, [190]GSE207422) and melanoma ([191]GSE78220 and
   PRJEB23709) were downloaded from GEO. RNA sequencing of the
   immunotherapy cohort of clear cell renal cell carcinoma (PMID32895571
   and PMID29301960) was downloaded from supplementary materials. POPLAR
   immunotherapy cohort of non -small-cell lung cancer was applied from
   European Genome-Phenome Archive
   ([192]https://ega-archive.org/studies/EGAS00001004343). RNA sequencing
   of the immunotherapy cohort of bladder cancer was downloaded from R
   package Mvigor210CoreBiologies [[193]42]. MSK-IMPACT assay of
   immunotherapy cohort was downloaded from cbioportal
   ([194]http://www.cbioportal.org/). RNA sequencing of patients with
   antiangiogenic treatment in non-small-cell lung cancer ([195]GSE61676),
   –triple-negative breast cancer ([196]GSE103668), colorectal cancer
   ([197]GSE60331 and [198]GSE53127) was downloaded from GEO. All public
   datasets are downloads of the processed datasets, so background
   correction, normalization and expression value calculation follow the
   original authors' processing.

4.2. Tumor mutational burden (TMB) estimates

   TMB was defined as the number of somatic non-synonymous single
   nucleotide variants. For the TCGA-LUAD dataset, somatic mutation data
   were retrieved from the UCSC Xena multi-omics database platform
   ([199]https://xenabrowser.net/datapages/?dataset=TCGA.LUAD.sampleMap%2F
   mutation_broad&host=https%3A%2F%2Ftcga.xenahubs.net&removeHub=https%3A%
   2F%2Fxena.treehouse.gi.ucsc.edu%3A443) and preprocessed at the Broad
   Institute Genome Sequencing Center [[200]43]. WES data were generated
   on the Illumina platform. Mutation calls were calculated using the
   MuTect method [[201]44], and only calls with variant allele frequency
   (VAF) > 4.0% were included. The R package “maftools” [[202]45] was then
   used to calculate the total number of somatic non-synonymous point
   mutations within each sample. The high and low TMB groups are cut off
   by the median. The copy number variation data were calculated by
   GISTIC2 [[203]46].

4.3. Immunohistochemical staining

   Six samples of IHC was available from the Human Protein Atlas (HPA,
   [204]https://www.proteinatlas.org/). This database contains the
   immunohistochemical (IHC) staining profiles of numerous proteins in a
   variety of cancerous and non-cancerous tissues based on more than 8800
   antibodies. All samples detailed information was attached in [205]Table
   S12.

4.4. Single-cell transcriptome sequencing and data preprocessing

   Seurat (v3.1.3) R toolkit was used to analyze the single-cell
   transcriptome sequencing data [[206]47]. Firstly, cells with low
   quality were filtered out [[207]48]. Briefly, the dead or dying cells
   with more than 10% mitochondrial RNA content were removed, and the
   cells with a too low number (less than 200) were also removed. Cell
   doublets were predicted using DoubletFinder with default parameters.
   Then, the filtered gene expression matrix for each sample was
   normalized using the “NormalizeData” function in Seurat with default
   parameters which employs a global-scaling normalization method
   “LogNormalize” then multiplies by a scale factor (10,000 by default)
   and log-transforms. Only highly variable genes remained using the
   “FindVariableFeatures” function in Seurat. Next, “Runharmony” functions
   in harmony were used to integrate the gene expression matrices of all
   samples, where batch effects between different samples have been
   adjusted [[208]49]. Next, the “RunPCA” function was used to perform the
   principal component analysis (PCA) and the “FindNeighbors” function was
   used to construct a K-nearest-neighbor graph. Next, the most
   representative principal components (PCs) selected based on PCA were
   used for clustering analysis with the “FindCluster” function to
   determine different cell types. Lastly, UMAP was used to visualize the
   different cell types.

   We annotated the cell types using the following rules: Based on the
   most ten differentially expressed genes that were derived using the
   “wilcoxauc” function in presto, genes such as CD3D, CD3G were used as T
   cell markers, SDC1, TNFRSF17 were used as plasma cell markers, CD68,
   CD14 were used as myeloid cell markers, MS4A2, CPA3 were used as mast
   cell markers, ACTA2, COL6A1 were used as fibroblast/smc cell markers,
   KRT18, EPCAM were used as epithelial cell markers, VWF, PECAM1 were
   used as endothelial cell markers, CD19, MS4A1 were used as B cell
   markers. CD4 and CD8 gene expression were used to differentiate CD4^+
   and CD8^+ T cells. Sub-cluster of CD4^+ T cells, CD8+T cells and
   endothelial were named by the first marker gene.

4.5. Pathway analysis of single-cell

   The pathway analysis was performed using the Correlation Adjusted MEan
   RAnk gene set test (CAMERA)that has been implemented in the
   SingleSeqGset (version 0.1.2) R package. Lists of metabolic genes and
   pathways were obtained from the KEGG database ([209]http://www.kegg.jp)
   were used for the CAMERA analysis [[210]48].

4.6. Tissue distribution of clusters

   We used to observe and expected cells in each cluster
   Ro/e = observed/expected, and expected cell numbers for each
   combination of cell clusters and tissues are obtained from the
   chi-square test. We defined R o/e > 1 as specific tissue cluster
   [[211]50].

4.7. Trajectory analysis

   Monocle 3 was used to infer the cell differentiation trajectories
   [[212]51]. These algorithms place the cells along with a trajectory
   corresponding to a biological process by taking advantage of an
   individual cell's asynchronous progression under an unsupervised
   framework.

4.8. Cell-cell communication analysis

   To investigate the potential cell-cell communications between T cells
   and endothelial cells from the tumor, we performed ligand-receptor
   analyses using CellPhoneDB [[213]52] (version 2.0.6). CellPhoneDB
   applies an algorithm that considers only receptors and ligands with
   broad expression among the tested cell types, followed by calculating
   the likelihood of cell-type specificity of a given receptor-ligand
   complex with a sufficient number of permutations.

4.9. Selection of candidate genes and establishment of tumor-specific T cell
signatures

   Tumor-specific T cell clusters can be defined by Roe. Genes expressed
   in these clusters and did not express in any other clusters were
   involved in the establishment of the TSTS score. We chose these genes
   for more functional investigation and developed a possible risk score
   using the lasso cox regression algorithm [[214]53].Finally, minimum
   criteria defined six genes and their constants, choosing the perfect
   penalty parameter λ related to the min 20fold-cross validation inside
   the training set. The EIC score was calculated:
   [MATH: <mrow><munderover><mo>∑</mo><mrow><mi>i</mi><mo
   linebreak="badbreak">=</mo><mn>1</mn></mrow><mi>n</mi></munderover><mro
   w><mo
   stretchy="true">(</mo><mrow><msub><mrow><mi>C</mi><mi>o</mi><mi>e</mi><
   mi>f</mi></mrow><mi>i</mi></msub><mo>*</mo><msub><mi>x</mi><mi>i</mi></
   msub></mrow><mo stretchy="true">)</mo></mrow></mrow> :MATH]

   Coefi is the coefficient and xi is the value of chosen genes.

4.10. Functional and pathway enrichment analysis of bulk RNA-Seq

   Limma R package was used to identify the signaling pathways that were
   differentially 473 activated between the low TSTS score and high TSTS
   score groups [[215]54]. Signal pathways specific to the tumor
   microenvironment (TME) phenotype between the high TSTS and low TSTS
   groups were explored by GSEA, Gene Set Enrichment Analysis (GSEA) was
   used with adjusted p < 0.05 using the clusterfiler R package [[216]55].

4.11. Immunological characteristics of the TME in lung cancer

   Evaluation of Immune Cell Infiltration The proportions of immune cell
   types (i.e., TICs) in each sample (with immune infiltration scores)
   were calculated using the Timer [[217]56]and MCP_counter algorithm
   [[218]57],EPIC, Xcell, quanTIseq. Dividing between the high and low
   groups for comparing the proportions of immune cell types in the text
   are thresholded at the median.

4.12. Deciphering mutational signature operative in the genome

   We used R packages Maftools to extract mutational signatures from the
   aggregated WES genomic data of NSCLC cohorts [[219]58]. The
   ExtractSignatures function from Maftools factorized the mutation
   portrait matrix into two nonnegative matrices, “signatures” and
   “contributions,” where “signatures” represented mutational processes
   and “contributions” represented the corresponding mutational
   activities. Mutational signatures were annotated by calculating cosine
   similarity against 30 validated mutational signatures in COSMIC,
   Version 2 [[220]59].

4.13. Spatial transcriptome analysis

   To infer the spatial organization of certain cell subtypes, we used
   Seurat R package to integrate spatial and single-cell data. Raw UMI
   counts were normalized by “SCTransform” function. Dimensionality
   reduction and clustering were performed as before. Cell subtypes
   distributions were visualized in spatial context over H&E images.

4.14. Statistical analysis

   Dividing between the high and low groups for all survival curves in the
   text are thresholded at the median and the sample sizes between high
   and low groups are indicated in the pictures. The Kaplan-Meier method
   was used to estimate OS or FPS, the log-rank test was used to compare
   the Kaplan-Meier curves. A two-sided P value of less than 0.05 was
   considered significant. All sample sizes were large enough to ensure
   proper statistical analysis. Pearson correlation analysis for the
   datasets that meet the normal distribution, Spearman correlation
   analysis for those that are not normally distributed (detailed
   statistical information of these dataset is attached to [221]Table
   S11). Statistical analyses were performed using GraphPad Prism
   (GraphPad Software, Inc.). P values < 0.05 were considered
   statistically significant. All t-test analyses were one-tailed t-tests
   (paired or unpaired depending on the experiments).

5. Conclusions

   In summary, we applied single-cell RNA sequencing to lung cancer
   patients and identified tumor-specific T cells and their features. The
   model constructed by LAYN, TOX, MAGEH1 and CXCR6 exclusively expressed
   in the tumor-specific T cells, which score represented how extend of
   the immune-inflammatory state and angiogenesis state, dividing lung
   cancer patients into two types. Patients with high TSTS are sensitive
   to immunotherapy and targeted therapy, while patients with low TSTS are
   sensitive to chemotherapy. Collectively, this article provides a basic
   explanation of the environmental context in which TSTS is used as a
   treatment decision-making strategy. And the limitation of our approach
   should also be taken into consideration. Although RNA sequencing, whole
   exon sequencing, and single-cell RNA sequencing were used and validated
   in as many independent cohorts as possible, all conclusions were based
   on bioinformatics analysis. Our results warrant further validation in a
   large number of real-world cohorts. The underlying mechanism should be
   further investigated.

Source(s) of support

   The study was supported by National Natural Science Foundation of China
   (81800050), and Natural science fund of Yangzhou City (YZ2017119).

Ethics approval

   Not applicable.

Consent to publish

   Not applicable.

Consent to participate

   Not applicable.

Data availability statement

   The processed Single-cell RNA-sequencing (ScRNA-seq) data of lung
   adenocarcinoma samples are available in GEO (accession number:
   [222]GSE131907, [223]GSE207422). The processed RNA sequencing data of
   lung cancer (number: [224]GSE68465, [225]GSE14814, [226]GSE30219,
   [227]GSE31210 and [228]GSE37745), non-small-cell lung cancer
   ([229]GSE135222, [230]GSE207422) and melanoma ([231]GSE78220 and
   PRJEB23709), and cohorts with antiangiogenic treatment in
   non-small-cell lung cancer ([232]GSE61676), –triple-negative breast
   cancer ([233]GSE103668), colorectal cancer ([234]GSE60331 and
   [235]GSE53127) are available in GEO. The processed RNA sequencing of
   the immunotherapy cohort of bladder cancer are available in R package
   Mvigor210CoreBiologies [[236]42]. POPLAR immunotherapy cohort of non
   -small-cell lung cancer are available in European Genome-Phenome
   Archive.

   ([237]https://ega-archive.org/studies/EGAS00001004343). MSK-IMPACT
   assay of immunotherapy cohort are available in cbioportal
   ([238]http://www.cbioportal.org/). The TCGA-LUAD RNA Sequencing data
   (TPM), somatic mutation data and survival information are available in
   the UCSC Xena data portal ([239]https://xenabrowser.net/datapages/).
   The processed RNA sequencing of the immunotherapy cohort of clear cell
   renal cell carcinoma (PMID32895571 and PMID29301960) are available in
   the original authors' supplementary materials. Spatial transcriptomic
   data for one patient with NSCLC are available in 10X
   ([240]https://www.10xgenomics.com/resources/datasets/human-lung-cancer-
   11-mm-capture-area-ffpe-2-standard).

   The source code for TSTS is available at
   [241]https://github.com/liuxuefei99/LungTSTS.

CRediT authorship contribution statement

   Ziwei Luo: Writing – review & editing, Writing – original draft,
   Visualization, Validation, Software, Project administration,
   Methodology, Investigation, Formal analysis, Data curation,
   Conceptualization. Xuefei Liu: Writing – review & editing,
   Visualization, Software, Project administration, Methodology, Data
   curation, Conceptualization. Ying Chen: Visualization, Software,
   Methodology, Data curation. Lize Shen: Visualization, Software, Data
   curation. Hui Qin: Visualization, Software, Methodology, Data curation.
   Qiongfang Zha: Validation, Software, Data curation. Feng Hu:
   Visualization, Methodology, Data curation. Yali Wang: Writing – review
   & editing, Writing – original draft, Funding acquisition,
   Conceptualization.

Declaration of competing interest

   The authors declare the following financial interests/personal
   relationships which may be considered as potential competing interests:
   Yali Wang reports financial support was provided by National Natural
   Science Foundation of China. Yali Wang reports financial support was
   provided by Natural science fund of Yangzhou City. If there are other
   authors, they declare that they have no known competing financial
   interests or personal relationships that could have appeared to
   influence the work reported in this paper.

Acknowledgements