Significance of CD8⁺ T cell infiltration-related biomarkers and the corresponding prediction model for the prognosis of kidney renal clear cell carcinoma

Results

De-batch effect of immune cells in GEO chips

To compare immune cell response in KIRC and identify cell type-specific markers and related pathways in the tumor-infiltrating CD8⁺T (TIL-CD8T) cells, transcriptome data was collected from immune cells across multiple studies. The expression results of the immune cell chip data was downloaded from the GEO database (Table 1), and the data underwent an inter-study bias correction by the inSilicoMerging software package. The chip data distribution before and after batch effect calibration, were displayed in (Figure 1A, 1B), respectively.

Table 1. Summary of immune cell types in the chip.

Population	Data set	Cell numbers
B cell activated	GSE28490, GSE49910	9
CD4 T cell activated	GSE28726, GSE49910	7
CD4 T cell resting	GSE28726	4
CD8 T cell activated	GSE49910	6
CD8 T cell resting	GSE49910	4
Dendritic cells activated	GSE59237	4
Dendritic cells resting	GSE59237	6
Eosinophils	GSE28698	3
Immature dendritic cells	GSE23371, GSE6863	6
Mast cells activated	GSE25320	4
Monocytes	GSE49910	6
Myeloid dendritic cells	GSE42058	4
Neutrophils	GSE39889, GSE49910	7
NK activated	GSE27838, GSE8059	11
NK resting	GSE8059	1
NKT activated	GSE28726	6
Plasmacytoid dendritic cells	GSE37750	8
T gamma delta	GSE13906, GSE27291	10
T helper 17	GSE51540	9

Figure 1. Elimination of inter-study bias in immune cell chip data. The batch effect correction of each chip was completed by the inSilicoMerging software package of Bioconductor. (A) Principal component analysis (PCA) analysis of each chip data before batch effect was removed; (B) PCA analysis of each chip data after the batch effect was removed; (C) The differential expression (|logFC|>1 and p<0.05) gene heat map of CD8 T cells and other immune cells.

Screening for specific marker genes of TIL-CD8 T cells

To explore CD8 T cell response upon tumor infiltration and identify a set of genes that were unique markers for TIL-CD8T cells, a systematic transcriptome analysis across multiple immune cells was performed by the Seurat's Findmarkers function. Genes with significant differential expressions (|logFC|>1 and p<0.05) were considered to be specifically expressed by TIL-CD8 T cells. In total of 63 genes were found to be up-regulated in CD8 T cells, while 396 genes were down-regulated (Figure 1C). As the next step, an analysis was put into practice to clarify whether these genes were correlated with patient survival data. Through the pan-cancer analysis, 12 of 63 up-regulated and 162 of 396 down-regulated genes were found to be significantly correlated with the survival prognosis (Supplementary Table 1). Annotations for all cell lines were provided in (Supplementary Material 1), while the details of all cell lines were summarized in (Supplementary Material 2).

GO and KEGG enrichment analyses of specific marker genes of TIL-CD8 T cells

In order to further investigate the biological function of TIL-CD8 T cell-specific marker genes, GO enrichment analysis and KEGG pathway enrichment analysis on the differentially expressed genes were performed. Specifically, the results of the top significant enrichment analyses for molecular function (MF), biological process (BP) and cellular component (CC) (Supplementary Figure 1A–1C) were evaluated, respectively. Neutrophil activationa and degranulation as well as T cell activation were included in the top biological GO categories. The enrichment analyses results of KEGG pathway were shown in (Supplementary Figure 1D), which chemokine signaling pathway was listed as the top hit.

TIL-CD8T-related prognostic factors in KIRC screened by univariate cox regression analysis

Then, the prognostic value of the identified marker genes was assessed. After filtering out patients with no survival information in the TCGA-KIRC cohort, a total of 307 samples with clinical survival data were further analyzed. Through univariate Cox regression analysis, 38 prognostic factors related to KIRC TIL-CD8T were obtained. The maximum selection rank sum statistics of the R software package "maxsat" was used to determine the cutoff value of the prognostic factors.

The Kaplan-Meier curve of top 6 factors (PDK4, MPP1, ASGR1, MS4A14, FCER1A, MX2) was shown in (Figure 2A–2F), and the cutoff values of the prognostic factors were shown in (Supplementary Table 2).

Figure 2. Kaplan-Meier survival curve for TIL-CD8T-related prognostic factors. The abscissa axis represents survival time; the ordinate axis represents survival probability. The survival curves of different colors represent different expression status of related genes. (A) Kaplan-Meier survival curve of PDK4 (P < 0.0001): Hazard Ratio=0.76, 95%CI [0.69, 0.84]; C-index= 0.66. (B) Kaplan-Meier survival curve of MPP1 (P < 0.0001): Hazard Ratio=0.60, 95%CI [0.45, 0.81]; C-index= 0.62. (C) Kaplan-Meier survival curve of ASGR1 (P = 0.00011): Hazard Ratio=1.34, 95%CI [1.17, 1.54]; C-index= 0.62. (D) Kaplan-Meier survival curve of MS4A14 (P = 0.00019): Hazard Ratio=1.33, 95%CI [1.12, 1.57]; C-index= 0.59. (E) Kaplan-Meier survival curve of FCER1A (P < 0.0001): Hazard Ratio=0.81, 95%CI [0.73, 0.89]; C-index= 0.63. (F) Kaplan-Meier survival curve of MX2 (P < 0.0001): Hazard Ratio=1.69, 95%CI [1.36, 2.09]; C-index= 0.61.

Construction of a TIL-CD8T-related tumor gene risk scoring model for KIRC (6 genes)

Through univariate Cox regression analysis, 38 TIL-CD8T-related genes were found to be significantly related to the prognosis of patients with KIRC (Supplementary Figure 2A). Then, the "cv.glmne" function of the R "glmnet" package was used for LASSO regression implement to further screen prognostic-related genes. The family parameter was set as "cox" and nfold=10. According to the lambda, the complexity of the model was optimized by LASSO regression. The "cv.glmne" analysis result contained two optimization models, one was "lambda.min" and the other was “lambda.1se”. The lambda between these two values was considered appropriate. The model constructed by "lambda.1se" was the simplest, which the number of involved genes was few, while the accuracy of "lambda.min" was higher. Here, 13 prognostic-related genes in the LASSO model corresponding to "lambda.min" was collected for follow-up analyses.

After the above analysis, 13 genes related to the prognosis of KIRC patients were identified (Supplementary Figure 2B, 2C). Then, based on these 13 genes, a multivariate Cox regression analysis was used for the TIL-CD8T risk scoring model construction. During the construction process, both the prediction accuracy and the simplicity of the model was taken into account. The "stepAIC" function in the "MASS" software package was used to screen the 13-variable factors in the multivariate Cox regression model. The parameter was set as "direction=both", "both" corresponded to the stepwise regression analysis method, starting from no predictor variables, and then adding the most contributory (calculate the "AIC" value according to the model after adding the variables, and select the smaller "AIC" value Variables) predictor variables in turn (such as forward selection). After adding each new variable, any variables that no longer provided improved model fit (the AIC value of the model could not be reduced after the variable was enrolled) would be deleted. Based on our multi-step iterative analyses, the following risk scoring model was constructed, which involved 6 genes:

TILCD8Sigscore = (-0.121 * PDK4) + (-0.492 * MPP1) + (0.22 * ASGR1) + (0.176 * MS4A14) + (-0.184 * FCER1A) + (0.5 * MX2)

Verification of the TIL-CD8T-related tumor gene risk scoring model by external data

In order to verify the prediction efficiency of the above model in an independent data set, the model was used to calculate the risk scores (TIL-CD8T-Sigscore) in KIRC sample chip data (GSE22541, GSE29609 and TCGA-KIRC). The disease-free survival (DFS) in GSE22541 was used as an index for clinical benefit evaluation. It was found that the survival prognosis of patients with low TIL-CD8T-Sigscore in these data sets was significantly better than those with high risk scores (Figure 3A–3B). The AUC of the 5-year OS, predicted by the TCGA data, was 0.765 (Figure 3C), while the AUC of the 5-year survival time predicted by the GSE22541 data was 0.629 (Figure 3D). In contrast, the survival prognosis of patients with low risk score of GSE29609 was worse than those of the high risk score group (Figure 3F). In order to further assess the correlation between clinical information and the individual expression of prognostic factors, this information was visualized in a heatmap. As shown in (Figure 3E), high risk group corresponded to the high mortality. This trend could also be found in TNM staging (Figure 3E).

Figure 3. Prognostic efficacy evaluation of risk scoring model and external data verification. (A) Kaplan-Meier overall survival curve of patients in the TCGA-KIRC cohort. The abscissa axis represents survival time; the ordinate axis represents survival probability. The survival curves of different colors represent different risk score subgroups. (B) Kaplan-Meier overall survival curve of patients in the GSE22541 cohort. The abscissa axis represents survival time; the ordinate axis represents survival probability. The survival curves of different colors represent different risk score subgroups. (C) The predictive efficiency of the risk scoring model in the TCGA-KIRC cohort (AUC= 0.765). The abscissa axis represents false positive rate; the ordinate axis represents true positive rate. (D) The predictive efficiency of the risk scoring model in the GSE22541 cohort (AUC= 0.629). The abscissa axis represents false positive rate; the ordinate axis represents true positive rate. (E) Heatmap representation of the expression levels of genes included in the KIRC scoring model of the low and the high-risk groups, and the distribution of clinicopathological characteristics in the low and high-risk groups. (F) Kaplan-Meier overall survival curve of patients with KIRC in the GSE29609 cohort. The abscissa axis represents survival time; the ordinate axis represents survival probability. The survival curves of different colors represent different risk score subgroups.

Genes expression profile included in the scoring model across KIRC cell lines

The expression profile of the 6 genes risk scoring model was compared between KIRC cell line and CD8 T cell data set, and then the Wilcoxon rank sum test was used to evaluate the differential expression significance of these genes in the CD8 T and KIRC cell lines. However, no significant expression difference or specificity corresponding to certain KIRC cell lines was found (Supplementary Figure 3A–3F).

Single-sample gene set enrichment analysis (ssGSEA) of immune cell infiltration level in TIL-CD8T Sig score group

Next, the immune infiltration levels in the high and low TIL-CD8T Sig score groups were further evaluated. ssGSEA was used in 19 immune cell subgroups. Dendritic and mast cells were found to be enriched in the low score group, while Act.B, Act.dendritic, Act.CD4 T, Act.CD8 T, Monocyte, Th17, NK.T, NK, Eosinophil, Plas.dendritic and Neutrophil cells were hyper-infiltrated in the high score group (Supplementary Figure 4A).

The landscape analysis of mutation and copy number variation (CNV) in high and low TIL-CD8Sig score groups

In order to clarify the correlation between TIL-CD8 Sig score and gene mutations, the R software package maftools was used to deal with the publicly available mutation annotation format (maf) files (Figure 4A), and evaluate the relationship between different TIL-CD8Sig scores and CNV (Figure 4B). It was found that the gene mutations and CNV status of patients in different TIL-CD8Sig score groups were different, and the CNV frequency was significantly different from each other (Figure 4C).

Figure 4. Genetic landscape analysis of mutation and copy number variation (CNV) in high and low TIL-CD8Sig score groups. (A) Mutation waterfall graph of different TIL-CD8Sig score subgroups in the TCGA-KIRC cohort. (B) CNV spectrum of different TIL-CD8Sig score subgroups in the TCGA-KIRC cohort. Different colors represent different CNV types (gain or loss); the abscissa axis represents chromosome locus; the ordinate axis represents CNV frequency. (C) The difference of mutation frequency in different TIL-CD8 Sig score subgroups. The abscissa axis represents different scoring groups; the ordinate axis represents mutation frequency.

Differences of the immunotherapy efficacy in patients with different TIL-CD8 Sig scores predicted by tumour immune dysfunction and exclusion (TIDE) score

In order to investigate the differences of the immunotherapy efficacy in patients with different TIL-CD8 Sig scores, the TIDE score was used for prediction practice. The TIDE method could integrate the expression characteristics of T cell dysfunction and rejection to simulate tumor immune escape, and then adopt untreated tumor data to predict the clinical response of immune checkpoint blockade (ICB). The TIDE online analysis website (http://tide.dfci.harvard.edu/) was used to calculate the TIDE score, and evaluate the difference among different TIL-CD8Sig score groups. However, no significant difference was found in the two groups (Supplementary Figure 4B).

The drug resistance (Sorafenib) prediction in TIL-CD8Sig score group

Sorafenib represents one of the standard treatment options for patients with metastatic renal cell carcinoma. However, highly treatment resistant is common in KIRC. The mechanism of drug resistance was still not well understood. To reveal the relationship between TIL-CD8 Sig score and drug resistance in patients with KIRC, the R software package pRRophetic was adopted for predicting the sensitivity of Sorafenib. Through analysis, it was found that the sensitivity of Sorafenib (Supplementary Figure 4C) in the low-scoring group was significantly different from that of the high-scoring group.

Independence verification of prognostic factors

To further verify the predictive value of TIL-CD8Sig score as a prognostic factor, Univariate and Multivariate Cox regression analyses were used to deal with the data of TCGA-KIRC and GSE22541. Through these analyses, it was found that the 6-gene TIL-CD8 Sig score model could be used as an independent prognostic factor in both TCGA and GSE22541 (Supplementary Figure 5A–5D).

Construction of nomogram for the prediction of overall survival rate

Nomograms had been previously shown as a reliable alternative tool that could help clinicians make individual predictions for survival. By integrating the TIL-CD8Sig score groups and clinicopathological risk factors, a nomogram for predicting the overall survival rate of TCGA-KIRC was constructed (Figure 5A). Through the standard curve diagram, it was found that the performance of the nomogram was equivalent to that of the ideal model (Figure 5B). Decision curve analysis (DCA) was used to quantify the clinical utility of it. For the overall survival probabilities of 2, 3, and 5 years, the decision curve showed that the nomogram could provide a better net benefit than the alternatives (Figure 5C).

Figure 5. Nomogram and decision analysis curve for predicting the overall survival of KIRC patients. (A) Combining the TIL-CD8Sig score of the TCGA-KIRC cohort data and clinical pathological risk factors to predict the 2-year, 3-year and 5-year overall survival probability. (B) According to the consistency of the prediction and observation results, the correction map of the nomogram was drawn. The performance of the nomogram was shown by the chart relative to the dotted line, which the dotted line indicated a perfect forecast. (C) Nomograph's decision analysis curve. None: Hypothetical events will not occur in any patients (horizontal solid line); All: Hypothetical events will occur in all patients (dotted line), the expected net income based on the nomogram prediction under different threshold probabilities was displayed by it.

Analyses of the combination of TIL-CD8 Sig score and immune check sites expression as the prognostic factors in KIRC

Immune checkpoint molecules, such as PD-1 and PD-L1, had been identified as crucial regulators of the immune response. However, the prognostic significance of these immune checkpoint molecules was still controversial. We wanted to further investigate the relationship between TIL-CD8 Sig score and the PD-1/PD-L1 expression levels. As shown in (Figure 6A), the TIL-CD8 Sig score was significantly correlated with both PD-1 and PD-L1 (PD-1: spearman correlation coefficient rho= 0.441, P <0.001; PD-L1: spearman correlation coefficient rho=0.129, P=0.024). To further evaluate the expression patterns of immune checkpoint genes in different TIL-CD8Sig score groups, it was found that the PD-1 (Wilcoxon rank sum test P<0.001) and PD-L1 (Wilcoxon rank sum test P<0.01) expression levels in the high-risk group were significantly higher than those of the low-risk group (Figure 6B).

Figure 6. The effect of TIL-CD8Sig and immune checkpoint gene expression on patient survival. (A) The relationship between TIL-CD8Sig and the expression of immune checkpoint genes (PD-1 and PD-L1). The abscissa axis represents the TIL-CD8Sig Score, and the ordinate axis represents the expression level of PD-1/PD-L1. (B) Box plot of the expression distribution of immune checkpoint genes (PD-1 and PD-L1) in the high and low risk groups of TIL-CD8Sig. The abscissa axis represents different risk groups; the ordinate axis represents the expression level of PD-1/PD-L1. (C) Kaplan-Meier survival curve of OS in four groups of patients stratified by TIL-CD8Sig and immune checkpoint gene expression. The abscissa axis represents survival time; the ordinate axis represents survival probability.

Finally, to further clarify the impact of the interaction between TIL-CD8Sig score and immune checkpoint genes on survival, the patients were divided into 4 groups based on the combination of TIL-CD8Sig score (high/low) and immune checkpoint genes (high/low). Survival comparisons were put into practice in the 4 groups, and the results were displayed in (Figure 6C). The comparison results showed that the prognosis of patients with the same or similar immune checkpoint gene levels could be distinguished from each other by TIL-CD8Sig score (PD-1: log rank P <0.0001; PD-L1: log rank P = <0.0001) (Figure 6C).

Abbreviations

TIL-CD8TSig: Tumor-Infiltrating CD8 T cell gene Signature; GDSC: Genomics of Drug Sensitivity in Cancer; CCLE: Cancer Cell Line Encyclopedia; GEO: Gene Expression Omnibus; TCGA: The Cancer Genome Atlas; LASSO: Least Absolute Shrinkage and Selection Operator; ssGSEA: Single-sample Gene Set Enrichment Analysis; NES: Normalized Enrichment Score; DCA: Decision Curve Analysis; GO: Gene Ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; TIDE: Tumour Immune Dysfunction and Exclusion; FPKM: Fragments Per Kilobase per Million; ICB: Immune Checkpoint Blockade; BP: Biological Process; MASS: Modern Applied Statistics with S; CNV: Copy Number Variation; AIC: Akaike Information Criterion.

Author Contributions

Yuan Tian, and Alan Huang collected the data; Yuan Tian, Alan Huang, and Qi Dang performed data cleaning and analysis; Yuan Tian drafted the manuscript; Yuan Tian and Yuping Sun reviewed the manuscript for scientific soundness. All authors reviewed the final draft and approved its submission. Yumei Wei, Hongmei Liu, Heli Shang, Yuedong Xu, Tong Wu and Wei Liu took part in the manuscript revisions.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Funding

This study was funded by the Academic Promotion Program of Shandong First Medical University (2019QL025; Yuping Sun), Natural Science Foundation of Shandong Province (ZR2019MH042; Yuping Sun), Jinan Science and Technology Program (201805064; Yuping Sun), and Postdoctoral Innovation Project of Jinan (Yuan Tian).

References

• 1. Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, Bray F. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin. 2021; 71:209–49. https://doi.org/10.3322/caac.21660 [PubMed]
• 2. Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer Statistics, 2021. CA Cancer J Clin. 2021; 71:7–33. https://doi.org/10.3322/caac.21654 [PubMed]
• 3. Cai Z, Liu Q. Understanding the Global Cancer Statistics 2018: implications for cancer control. Sci China Life Sci. 2021; 64:1017–20. https://doi.org/10.1007/s11427-019-9816-1 [PubMed]
• 4. Feng RM, Zong YN, Cao SM, Xu RH. Current cancer situation in China: good or bad news from the 2018 Global Cancer Statistics? Cancer Commun (Lond). 2019; 39:22. https://doi.org/10.1186/s40880-019-0368-6 [PubMed]
• 5. Zhang J, Wu T, Simon J, Takada M, Saito R, Fan C, Liu XD, Jonasch E, Xie L, Chen X, Yao X, Teh BT, Tan P, et al. VHL substrate transcription factor ZHX2 as an oncogenic driver in clear cell renal cell carcinoma. Science. 2018; 361:290–95. https://doi.org/10.1126/science.aap8411 [PubMed]
• 6. Hsieh JJ, Le VH, Oyama T, Ricketts CJ, Ho TH, Cheng EH. Chromosome 3p Loss-Orchestrated VHL, HIF, and Epigenetic Deregulation in Clear Cell Renal Cell Carcinoma. J Clin Oncol. 2018; 36:3533–39. https://doi.org/10.1200/JCO.2018.79.2549 [PubMed]
• 7. Bakouny Z, Barbie DA. TBK1 Activation by VHL Loss in Renal Cell Carcinoma: A Novel HIF-Independent Vulnerability. Cancer Discov. 2020; 10:348–50. https://doi.org/10.1158/2159-8290.CD-19-1525 [PubMed]
• 8. Gao YH, Wu ZX, Xie LQ, Li CX, Mao YQ, Duan YT, Han B, Han SF, Yu Y, Lu HJ, Yang PY, Xu TR, Xia JL, et al. VHL deficiency augments anthracycline sensitivity of clear cell renal cell carcinomas by down-regulating ALDH2. Nat Commun. 2017; 8:15337. https://doi.org/10.1038/ncomms15337 [PubMed]
• 9. Ricketts CJ, Linehan WM. Multi-regional Sequencing Elucidates the Evolution of Clear Cell Renal Cell Carcinoma. Cell. 2018; 173:540–42. https://doi.org/10.1016/j.cell.2018.03.077 [PubMed]
• 10. Rabadán R, Mohamedi Y, Rubin U, Chu T, Alghalith AN, Elliott O, Arnés L, Cal S, Obaya ÁJ, Levine AJ, Cámara PG. Identification of relevant genetic alterations in cancer using topological data analysis. Nat Commun. 2020; 11:3808. https://doi.org/10.1038/s41467-020-17659-7 [PubMed]
• 11. Riaz N, Blecua P, Lim RS, Shen R, Higginson DS, Weinhold N, Norton L, Weigelt B, Powell SN, Reis-Filho JS. Pan-cancer analysis of bi-allelic alterations in homologous recombination DNA repair genes. Nat Commun. 2017; 8:857. https://doi.org/10.1038/s41467-017-00921-w [PubMed]
• 12. Buckley AR, Ideker T, Carter H, Harismendy O, Schork NJ. Exome-wide analysis of bi-allelic alterations identifies a Lynch phenotype in The Cancer Genome Atlas. Genome Med. 2018; 10:69. https://doi.org/10.1186/s13073-018-0579-5 [PubMed]
• 13. Zhao Q, Shi X, Xie Y, Huang J, Shia B, Ma S. Combining multidimensional genomic measurements for predicting cancer prognosis: observations from TCGA. Brief Bioinform. 2015; 16:291–303. https://doi.org/10.1093/bib/bbu003 [PubMed]
• 14. Hao X, Luo H, Krawczyk M, Wei W, Wang W, Wang J, Flagg K, Hou J, Zhang H, Yi S, Jafari M, Lin D, Chung C, et al. DNA methylation markers for diagnosis and prognosis of common cancers. Proc Natl Acad Sci USA. 2017; 114:7414–19. https://doi.org/10.1073/pnas.1703577114 [PubMed]
• 15. Zeng D, Li M, Zhou R, Zhang J, Sun H, Shi M, Bin J, Liao Y, Rao J, Liao W. Tumor Microenvironment Characterization in Gastric Cancer Identifies Prognostic and Immunotherapeutically Relevant Gene Signatures. Cancer Immunol Res. 2019; 7:737–50. https://doi.org/10.1158/2326-6066.CIR-18-0436 [PubMed]
• 16. Chevrier S, Levine JH, Zanotelli VR, Silina K, Schulz D, Bacac M, Ries CH, Ailles L, Jewett MA, Moch H, van den Broek M, Beisel C, Stadler MB, et al. An Immune Atlas of Clear Cell Renal Cell Carcinoma. Cell. 2017; 169:736–49.e18. https://doi.org/10.1016/j.cell.2017.04.016 [PubMed]
• 17. Miao D, Margolis CA, Gao W, Voss MH, Li W, Martini DJ, Norton C, Bossé D, Wankowicz SM, Cullen D, Horak C, Wind-Rotolo M, Tracy A, et al. Genomic correlates of response to immune checkpoint therapies in clear cell renal cell carcinoma. Science. 2018; 359:801–06. https://doi.org/10.1126/science.aan5951 [PubMed]
• 18. Wang W, Green M, Choi JE, Gijón M, Kennedy PD, Johnson JK, Liao P, Lang X, Kryczek I, Sell A, Xia H, Zhou J, Li G, et al. CD8⁺ T cells regulate tumour ferroptosis during cancer immunotherapy. Nature. 2019; 569:270–74. https://doi.org/10.1038/s41586-019-1170-y [PubMed]
• 19. van der Leun AM, Thommen DS, Schumacher TN. CD8⁺ T cell states in human cancer: insights from single-cell analysis. Nat Rev Cancer. 2020; 20:218–32. https://doi.org/10.1038/s41568-019-0235-4 [PubMed]
• 20. Held W, Speiser DE. Not All Tumor-Infiltrating CD8⁺ T Cells Are Created Equal. Cancer Cell. 2021; 39:145–47. https://doi.org/10.1016/j.ccell.2021.01.015 [PubMed]
• 21. Thompson ED, Zahurak M, Murphy A, Cornish T, Cuka N, Abdelfatah E, Yang S, Duncan M, Ahuja N, Taube JM, Anders RA, Kelly RJ. Patterns of PD-L1 expression and CD8 T cell infiltration in gastric adenocarcinomas and associated immune stroma. Gut. 2017; 66:794–801. https://doi.org/10.1136/gutjnl-2015-310839 [PubMed]
• 22. Gómez-Aleza C, Nguyen B, Yoldi G, Ciscar M, Barranco A, Hernández-Jiménez E, Maetens M, Salgado R, Zafeiroglou M, Pellegrini P, Venet D, Garaud S, Trinidad EM, et al. Inhibition of RANK signaling in breast cancer induces an anti-tumor immune response orchestrated by CD8+ T cells. Nat Commun. 2020; 11:6335. https://doi.org/10.1038/s41467-020-20138-8 [PubMed]
• 23. Dong MB, Wang G, Chow RD, Ye L, Zhu L, Dai X, Park JJ, Kim HR, Errami Y, Guzman CD, Zhou X, Chen KY, Renauer PA, et al. Systematic Immunotherapy Target Discovery Using Genome-Scale In Vivo CRISPR Screens in CD8 T Cells. Cell. 2019; 178:1189–204.e23. https://doi.org/10.1016/j.cell.2019.07.044 [PubMed]
• 24. Liu J, Lichtenberg T, Hoadley KA, Poisson LM, Lazar AJ, Cherniack AD, Kovatich AJ, Benz CC, Levine DA, Lee AV, Omberg L, Wolf DM, Shriver CD, et al, and Cancer Genome Atlas Research Network. An Integrated TCGA Pan-Cancer Clinical Data Resource to Drive High-Quality Survival Outcome Analytics. Cell. 2018; 173:400–16.e11. https://doi.org/10.1016/j.cell.2018.02.052 [PubMed]
• 25. Alimohamed S, Geller JI, Aronow B. TCGA Kidney Cell Atlas: A differential gene expression database for the dissection of tumor-specific gene modules. J Clin Oncol. 2020; 38:e17078. https://doi.org/10.1200/jco.2020.38.15_suppl.e17078
• 26. Lin J, Yu M, Xu X, Wang Y, Xing H, An J, Yang J, Tang C, Sun D, Zhu Y. Identification of biomarkers related to CD8⁺ T cell infiltration with gene co-expression network in clear cell renal cell carcinoma. Aging (Albany NY). 2020; 12:3694–712. https://doi.org/10.18632/aging.102841 [PubMed]
• 27. Qi Y, Xia Y, Lin Z, Qu Y, Qi Y, Chen Y, Zhou Q, Zeng H, Wang J, Chang Y, Bai Q, Wang Y, Zhu Y, et al. Tumor-infiltrating CD39⁺ CD8⁺ T cells determine poor prognosis and immune evasion in clear cell renal cell carcinoma patients. Cancer Immunol Immunother. 2020; 69:1565–76. https://doi.org/10.1007/s00262-020-02563-2 [PubMed]
• 28. Shaul ME, Fridlender ZG. Tumour-associated neutrophils in patients with cancer. Nat Rev Clin Oncol. 2019; 16:601–20. https://doi.org/10.1038/s41571-019-0222-4 [PubMed]
• 29. Nicolás-Ávila JÁ, Adrover JM, Hidalgo A. Neutrophils in Homeostasis, Immunity, and Cancer. Immunity. 2017; 46:15–28. https://doi.org/10.1016/j.immuni.2016.12.012 [PubMed]
• 30. Tsukita Y, Fujino N, Miyauchi E, Saito R, Fujishima F, Itakura K, Kyogoku Y, Okutomo K, Yamada M, Okazaki T, Sugiura H, Inoue A, Okada Y, Ichinose M. Axl kinase drives immune checkpoint and chemokine signalling pathways in lung adenocarcinomas. Mol Cancer. 2019; 18:24. https://doi.org/10.1186/s12943-019-0953-y [PubMed]
• 31. Mellado M, Rodríguez-Frade JM, Mañes S, Martínez-A C. Chemokine signaling and functional responses: the role of receptor dimerization and TK pathway activation. Annu Rev Immunol. 2001; 19:397–421. https://doi.org/10.1146/annurev.immunol.19.1.397 [PubMed]
• 32. Villarino AV, Kanno Y, O’Shea JJ. Mechanisms and consequences of Jak-STAT signaling in the immune system. Nat Immunol. 2017; 18:374–84. https://doi.org/10.1038/ni.3691 [PubMed]
• 33. Brooks AJ, Putoczki T. JAK-STAT Signalling Pathway in Cancer. Cancers (Basel). 2020; 12:1971. https://doi.org/10.3390/cancers12071971 [PubMed]
• 34. Sanchez-Vega F, Mina M, Armenia J, Chatila WK, Luna A, La KC, Dimitriadoy S, Liu DL, Kantheti HS, Saghafinia S, Chakravarty D, Daian F, Gao Q, et al, and Cancer Genome Atlas Research Network. Oncogenic Signaling Pathways in The Cancer Genome Atlas. Cell. 2018; 173:321–37.e10. https://doi.org/10.1016/j.cell.2018.03.035 [PubMed]
• 35. Rakoff-Nahoum S, Medzhitov R. Toll-like receptors and cancer. Nat Rev Cancer. 2009; 9:57–63. https://doi.org/10.1038/nrc2541 [PubMed]
• 36. Kuenzi BM, Ideker T. A census of pathway maps in cancer systems biology. Nat Rev Cancer. 2020; 20:233–46. https://doi.org/10.1038/s41568-020-0240-7 [PubMed]
• 37. Nagarsheth N, Wicha MS, Zou W. Chemokines in the cancer microenvironment and their relevance in cancer immunotherapy. Nat Rev Immunol. 2017; 17:559–72. https://doi.org/10.1038/nri.2017.49 [PubMed]
• 38. Davies H, Glodzik D, Morganella S, Yates LR, Staaf J, Zou X, Ramakrishna M, Martin S, Boyault S, Sieuwerts AM, Simpson PT, King TA, Raine K, et al. HRDetect is a predictor of BRCA1 and BRCA2 deficiency based on mutational signatures. Nat Med. 2017; 23:517–25. https://doi.org/10.1038/nm.4292 [PubMed]
• 39. Feng B, Shen Y, Pastor Hostench X, Bieg M, Plath M, Ishaque N, Eils R, Freier K, Weichert W, Zaoui K, Hess J. Integrative Analysis of Multi-omics Data Identified EGFR and PTGS2 as Key Nodes in a Gene Regulatory Network Related to Immune Phenotypes in Head and Neck Cancer. Clin Cancer Res. 2020; 26:3616–28. https://doi.org/10.1158/1078-0432.CCR-19-3997 [PubMed]
• 40. Reichling C, Taieb J, Derangere V, Klopfenstein Q, Le Malicot K, Gornet JM, Becheur H, Fein F, Cojocarasu O, Kaminsky MC, Lagasse JP, Luet D, Nguyen S, et al. Artificial intelligence-guided tissue analysis combined with immune infiltrate assessment predicts stage III colon cancer outcomes in PETACC08 study. Gut. 2020; 69:681–90. https://doi.org/10.1136/gutjnl-2019-319292 [PubMed]
• 41. Li Z, Peng Y, Li J, Chen Z, Chen F, Tu J, Lin S, Wang H. N⁶-methyladenosine regulates glycolysis of cancer cells through PDK4. Nat Commun. 2020; 11:2578. https://doi.org/10.1038/s41467-020-16306-5 [PubMed]
• 42. Wittig R, Nessling M, Will RD, Mollenhauer J, Salowsky R, Münstermann E, Schick M, Helmbach H, Gschwendt B, Korn B, Kioschis P, Lichter P, Schadendorf D, Poustka A. Candidate genes for cross-resistance against DNA-damaging drugs. Cancer Res. 2002; 62:6698–705. [PubMed]
• 43. Chen C, Li K, Jiang H, Song F, Gao H, Pan X, Shi B, Bi Y, Wang H, Wang H, Li Z. Development of T cells carrying two complementary chimeric antigen receptors against glypican-3 and asialoglycoprotein receptor 1 for the treatment of hepatocellular carcinoma. Cancer Immunol Immunother. 2017; 66:475–89. https://doi.org/10.1007/s00262-016-1949-8 [PubMed]
• 44. Choi J, Zhang T, Vu A, Ablain J, Makowski MM, Colli LM, Xu M, Hennessey RC, Yin J, Rothschild H, Gräwe C, Kovacs MA, Funderburk KM, et al. Massively parallel reporter assays of melanoma risk variants identify MX2 as a gene promoting melanoma. Nat Commun. 2020; 11:2718. https://doi.org/10.1038/s41467-020-16590-1 [PubMed]
• 45. Gossage L, Eisen T, Maher ER. VHL, the story of a tumour suppressor gene. Nat Rev Cancer. 2015; 15:55–64. https://doi.org/10.1038/nrc3844 [PubMed]
• 46. Kudo M, Ueshima K, Yokosuka O, Ogasawara S, Obi S, Izumi N, Aikata H, Nagano H, Hatano E, Sasaki Y, Hino K, Kumada T, Yamamoto K, et al, and SILIUS study group. Sorafenib plus low-dose cisplatin and fluorouracil hepatic arterial infusion chemotherapy versus sorafenib alone in patients with advanced hepatocellular carcinoma (SILIUS): a randomised, open label, phase 3 trial. Lancet Gastroenterol Hepatol. 2018; 3:424–32. https://doi.org/10.1016/S2468-1253(18)30078-5 [PubMed]
• 47. Escudier B, Eisen T, Stadler WM, Szczylik C, Oudard S, Siebels M, Negrier S, Chevreau C, Solska E, Desai AA, Rolland F, Demkow T, Hutson TE, et al, and TARGET Study Group. Sorafenib in advanced clear-cell renal-cell carcinoma. N Engl J Med. 2007; 356:125–34. https://doi.org/10.1056/NEJMoa060655 [PubMed]
• 48. Braun DA, Bakouny Z, Hirsch L, Flippot R, Van Allen EM, Wu CJ, Choueiri TK. Beyond conventional immune-checkpoint inhibition - novel immunotherapies for renal cell carcinoma. Nat Rev Clin Oncol. 2021; 18:199–214. https://doi.org/10.1038/s41571-020-00455-z [PubMed]
• 49. Braun DA, Hou Y, Bakouny Z, Ficial M, Sant’ Angelo M, Forman J, Ross-Macdonald P, Berger AC, Jegede OA, Elagina L, Steinharter J, Sun M, Wind-Rotolo M, et al. Interplay of somatic alterations and immune infiltration modulates response to PD-1 blockade in advanced clear cell renal cell carcinoma. Nat Med. 2020; 26:909–18. https://doi.org/10.1038/s41591-020-0839-y [PubMed]
• 50. Dai S, Zeng H, Liu Z, Jin K, Jiang W, Wang Z, Lin Z, Xiong Y, Wang J, Chang Y, Bai Q, Xia Y, Liu L, et al. Intratumoral CXCL13⁺ CD8⁺ T cell infiltration determines poor clinical outcomes and immunoevasive contexture in patients with clear cell renal cell carcinoma. J Immunother Cancer. 2021; 9:e001823. https://doi.org/10.1136/jitc-2020-001823 [PubMed]