Abstract

Introduction

Lung adenocarcinoma (LUAD) is the most predominant type of lung cancer, with a high invasion and mortality rate [1] . Presently, with a thorough understanding of LUAD development, earlier diagnosis, earlier detection, and diverse treatments have been implemented for patients. However, dismal median overall survival (OS) and 5-year survival rates of LUAD patients persist [2] . Moreover, substantial prognostic differences exist in LUAD patients with similar clinical characteristics, indicating the unsatisfactory risk stratification ability of traditional clinical classification strategies based on clinicopathological features [3–5] . Given these challenges, it is imperative for clinicians to elucidate tumor heterogeneity and develop novel stratification approaches to improve prognosis and treatment efficacy.

The circadian rhythm is a vital biological mechanism in almost all organisms, which coordinates multiple physiological and behavioral processes through the construction of the circadian clock [6,7] . In recent years, mounting research demonstrates that circadian rhythm disruption (CRD) and altered hub circadian gene expression are linked to abnormalities in cell metabolism, cell proliferation, tumor microenvironment (TME), and intratumoral heterogeneity, which contribute to cancer development and progression [8,10] . For instance, enhanced stemness of tumor cells and an immunosuppressive TME were detected in breast cancer mice with chronic circadian disruption [11] . A previous study revealed the existence of intratumoral heterogeneity and related resistance to anti-cancer treatments in LUAD patients with CRD based on single-cell RNA-seq analysis. Besides, Ruan et al. indicated that CRD was a potential target to facilitate anti-tumor therapeutic efficacy [12] . With the continuous advancements in circadian rhythm research, several CRD-related genes have been detected, such as BMAL1, CLOCK, PER, and CRY [13,14] . Nevertheless, the relationship among CRD status, molecular characteristics, and clinical outcomes in LUAD remains to be elucidated.

Besides, with the advancement of tumor research, plentiful effective therapies have been developed (e.g., chemotherapy, radiation therapy, targeted therapy, and immunotherapy) [5] . Due to the provision of diverse treatment options, the need for individualized treatment and precision medicine ensues [15,16] . Obviously, a traditional therapy strategy with an insufficient understanding of molecular characteristics was powerless for this requirement [17,19] . Thus, individualized comprehensive treatment, which included novel approaches such as immunotherapy and targeted therapy, was brought to the forefront [20,21] . However, an effective and rational classification is the essential prerequisite to determining the appropriate treatments. Therefore, it is warranted to identify CRD status heterogeneity and propose new insights for molecular classification, which could offer proper clinical management and precision medicine to LUAD patients.

In our present study, we aimed to address a significant clinical gap in the stratification of LUAD patients by identifying three heterogeneous subtypes based on the expression of CRD-related genes. By utilizing six independent databases, we validated the robustness of the CRD-related subtypes, demonstrating consistent relative fractions, gene expression profiles, and prognostic outcomes across cohorts. Furthermore, we explored the distinct differences among the three subtypes from multiple perspectives, including biological functions, genomic variations, and tumor microenvironment (TME) characteristics. This comprehensive approach allowed us to illustrate how CRD-based classification uniquely adds clinical value compared to existing biomarkers. We also assessed patients' sensitivity to common clinical therapies, which is essential for advancing personalized medicine strategies. Overall, our findings suggest that LUAD patients may benefit from more efficient and precise management if our promising CRD-based stratification platform is implemented in clinical practice.

Materials and Methods

Data acquisition and procession

In this study, LUAD cohorts were collected from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) dataset according to the following inclusion criteria: 1) Patients were primary lung adenocarcinoma; 2) The number of patients is more than 100; 3) Patients in the dataset have complete gene expression profiles and corresponding survival information; 4) The probe and gene ID were clearly labeled, with more than 20000 genes. Finally, 2090 patients from TCGA-LUAD (n = 497), GSE72094 (n = 442), GSE68465 (n = 462), GSE50081 (n = 127), GSE42127 (n = 133), GSE41271 (n = 183) and GSE31210 (n = 246) were included. The expression profile and clinical information of the TCGA-LUAD cohort were from UCSC Xena portal. The expression data were converted from fragments per kilobase of million mapped reads (FPKM) transcripts to trans per million (TPM) format and log2 transformed. The remaining cohorts were retrieved from GEO, normalized, and processed using Affy and Lumi packages based on different platforms subsequently [22] . Furthermore, TCGA-LUAD somatic mutation and segmented copy number variation data were received from the TCGA portal. Across all cohorts, the expression of each gene was converted into a Z-score value before model construction.

Development and validation of CRD-related subtypes

With the help of previous studies, we retrieved a total of 2091 circadian rhythm disruption-related genes from CircaDB and MSigDB for the development of CRD-related subtypes (Table S1) [23] . Nonnegative matrix factorization (NMF) algorithm was applied to identify the optimal number of consensus clusters from the TCGA-LUAD cohort through the NMF package. After decomposing the nonnegative matrix of CRD-related genes and iterating, the cophenetic coefficient was executed to determine the optimal factorization rank. The identification criteria were as follows: possible factorization ranks = 2-7, number of iterations = 100, and method = "lee". In general, the rank before the most obvious decrease of the cophenetic coefficient value was considered to be the optimal rank [24] . Then, Partial Least Squares Discriminant Analysis (PLS-DA) was applied to evaluate the separation of three subtypes. PLS-DA was a linear classification method that uses partial least squares regression to identify latent variables that maximize class separation. It is particularly effective in high-dimensional data settings. The limma package was used to decipher the differences among distinct CRD-related subtypes and obtain signature genes for each subtype (log2 fold change (log2 FC) > 1 and adjust P value < 0.05). A flexible technique, the nearest template prediction (NTP), was a helpful tool to assess class prediction confidence for single patient [25,26] . Using signature genes, we implemented the NTP algorithm with the CMScaller package to evaluate the stability and robustness of clusters across multiple GEO validation cohorts from different platforms. An FDR threshold of less than 0.05 was applied to establish appropriate classification confidence thresholds.

Exploring specific biological characteristics of three subtypes

To further explore the biological functional heterogeneity among the three subtypes, we conducted the gene set variation analysis (GSVA), which was widely applied to evaluate the activity of pathways [16,27] . GSVA was a method used to assess variability in gene sets across samples or conditions. GSVA transforms gene expression data into enrichment scores for gene sets, helping to reveal changes in the activity of biological processes or signaling pathways. Differentially expressed genes (DEGs) were analyzed by limma package and the expression matrix was obtained by arranging all genes in descending order according to log2 FC. Based on this matrix and gene sets from Gene Ontology (GO), the Kyoto Encyclopedia of Genes and Genomes (KEGG), and HALLMARK, the GSVA package was implemented to determine specific biological characteristics between each subtype and the others. Meanwhile, the gene set enrichment analysis (GSEA) algorithm was performed to exhibit CRD-related pathway activities through the clusterProfiler package.

Somatic mutation and copy number variation analysis

The maftools package was applied for processing and visualizing the genomic alteration data. Based on CNV data obtained from GISTIC 2.0 pipeline, the burden of copy number alteration, including amplification and deletion, was quantified at focal and arm levels. We also calculated the fraction of genome alteration (FGA), fraction of genomic gained (FGG), and fraction of genome lost (FGL) to evaluate genetic changes in three CRD-related clusters.

Depicting distinct immune landscape and evaluating Immunotherapy

For deciphering features of the tumor immune infiltration, single-sample gene set enrichment analysis (ssGSEA) was exploited for the quantification of 28 immune cell subsets [28] . Besides, GSVA package was conducted to evaluate the relative infiltration of 24 TIME cells. Six other algorithms including ESTIMATE, TIMER, quantTIseq, MCP counter, EPIC, and xCell were further implemented to verify the stability and reliability of the results. Meanwhile, the assessment of immunogenicity and immunosuppression status was accomplished by computing Immunophenoscore (IPS) and novel S score respectively (Table S2) [29] . The IPS of TCGA-LUAD patients was acquired from the cancer-immune group atlas (TCIA, https://tcia.at/home). Human leukocyte antigen (HLA) molecule expression was compared to assess the antigen presentation ability of three subtypes [30] . A range of immune escape-related signatures was collected and estimated to reveal underlying distinct immune escape mechanisms among three subtypes [31] .

For distinct subtypes, immune checkpoint molecules (ICM) expression, T cell inflammatory signature (TIS), and subclass mapping (Submap) algorithm were employed to deduce immunotherapeutic efficacy. A total of 27 ICM were enrolled, including the B7-CD28 superfamily, TNF superfamily, and eight other molecules [32] . As a signature obtained based on ssGSEA algorithm for 18 inflammatory genes, TIS could play a predictor of the response to PD-1 inhibitors. GSEA was an extension of traditional gene set enrichment analysis that evaluates the enrichment of gene sets in individual samples. Unlike GSEA, ssGSEA provides specific scores for each sample, allowing for the analysis of subtle differences between samples [33-35] . A higher TIS score represents a better response to PD-1 inhibitors. Moreover, the Submap algorithm was utilized to estimate the similarity between the three phenotypes and the patients with different immunotherapy responses from two independent immunotherapy cohorts.

Personalized management for distinct subtypes

Among the considerable number of published LUAD signatures, we aimed to identify an optimal signature tailored for distinct subtypes to facilitate personalized management. In our research, we retrieved a total of 151 published signatures that were based on various biological processes [36] . To strengthen the rationale for their use, we employed a systematic approach that included univariate Cox regression analysis, followed by a comparison of the concordance index (C-index) across these candidates. This allowed us to identify and prioritize the superior prognostic signatures for each subtype.

Statistical analysis

All data processing, plotting, and statistical analysis were conducted in R 4.1.2. Cox regression and Kaplan-Meier analyses were performed via the survival package. The comparison of the survival of categorical variables was completed through the log-rank test. Kruskal-Wallis test was applied to compare the difference among the three clusters. A two-sided P < 0.05 was considered statistically significant for all statistical tests.

Result

Identification of three CRD-related subtypes

Based on expression profiles of 2091 CRD-related genes, the NMF approach was applied to decode heterogeneous phenotypes. We selected three subtypes as the optimal rank based on the cophenetic coefficient score and consensus matrix (Figure S1A, 1A). Meanwhile, PLS-DA exhibited the obvious separation of three distinct subtypes (Figure 1B). We further explored the prognostic value of CRD-related subtypes to boost their practice in the clinic. For three subtypes called C1, C2, and C3, C2 displayed better overall survival (OS), while C1 harbored a dismal prognosis (Figure 1C).

Validation of CRD-related subtypes

To further demonstrate the stability and reliability of CRD-related subtypes, NTP analysis was performed in five independent cohorts, including GSE72094, GSE68465, GSE41271, GSE50081, GSE42127, and GSE31210 (Figure 1D-F, S1B-D). Signature genes for NTP were defined as specific upregulated DEGs in the individual subtype. Corresponding with a previous study, patients with a false discovery rate (FDR) of more than 0.05 were eliminated in subsequential analysis. Then, a resemblant proportion of three clusters in divergent datasets was exhibited, which hinted stability of CRD-related subtypes (Figure 2A). Moreover, we evaluated the difference in clinical outcomes for distinct subtypes through Kaplan-Meier curves in validation cohorts (Figure 2B-F, S2A). Relatively speaking, C2 subtypes possessed the most favorable OS (P < 0.05), whereas C1 exhibited the most frustrating OS, which was consistent with the above findings. Taken together, CRD-related subtypes were stable and robust in LUAD patients and revealed discrepant clinical prognoses.

Underlying biological functional processes linked to three subtypes

Based on GO, KEGG, and HALLMARK, enrichment analysis was performed to decode specific biological functions of LUAD patients in distinct subtypes. As illustrated in Figure 3A, significantly activated proliferation pathways were revealed in C1, including cell cycle, DNA replication, and E2F targets. C2 possessed a high correlation with lipid metabolisms such as bile acid metabolism and fatty acid metabolism. Notably, positive regulation of the circadian sleep-wake cycle was detected in C2 patients (Figure 3B), which demonstrated a healthy circadian rhythm. For tumor patients with CRD, clock ablation supports glycolysis and fatty acid synthesis, which was the signature of proliferative metabolism. Consistent with that, patients in C1 displayed negative regulation of the circadian sleep-wake cycle (Figure S2B), which signified its CRD status. However, elevated lipid metabolism activity impeded such a metabolism change, which coincided with depressed proliferation pathways in C2. Meanwhile, we observed C3 mainly harbored enhanced immune response and up-regulated signaling pathways of immune factors including cytokine, chemokine, and complement. Therefore, C3 was characterized as immune LUAD. Analogously, we defined C1 as proliferative LUAD and C2 as lipid-metabolic LUAD.

Genomic alterations of three distinct subtypes

By displaying and comparing the mutation frequency of the top 20 mutant genes in three subtypes, we detected C1 patients possessed the highest mutated frequency, especially TP53 and TTN mutations (Figure 4A). Meanwhile, an overview of single-nucleotide polymorphism (SNP), insertion and deletion (INDEL), and tumor mutation burden (TMB) were displayed to describe somatic variants comprehensively. C1 demonstrated the richest somatic mutations, which was consistent with the result of mutated frequency (Figure 4A). As a general mutation in LUAD, TP53 mutation conferred more vigorous malignant proliferation and poorer prognosis to patients [37-39] . In line with that, a significantly higher stemness index of C1 was indicated, which represented powerful malignant proliferation ability (Figure S2C). C2 had lower somatic mutation compared to other subtypes, implying better outcomes. A comparison of CNV in three CRD-related subtypes was also performed to further decode genomic alterations. Strikingly, C1 harbored obvious CNV at bases, fragments, and chromosome levels, which suggested a higher likelihood of cell proliferation and immune escape (Figure 4B) [40] . Overall, prominent genomic alterations were revealed for patients in C1, suggesting a subtype with high genomic instability.

The depiction of immune infiltration and immune escape landscape

Immune infiltration and immune escape play an essential role in tumorigenesis and the development and prognosis of patients. Therefore, we depicted the landscape of immune infiltration and immune escape in three heterogeneous subtypes. Firstly, we quantified the relative abundance of 28 immune cells in three subtypes. Consistent with our results through multiple algorithms, C3 showed a more abundant immune-cell infiltration than the other two subtypes (Figure 5A, Figure S2D). However, IPS and the S score hinted lowest immunity and conspicuous immunosuppression of C3 patients, which might be associated with elevated infiltration of immunosuppressive cells, such as MDSC, Treg, and Th17. Conversely, C2 possessed vigorous immunity and the least immunosuppression (Figure 5B-C). Besides, we explored underlying mechanisms of immune escape for the three subtypes. Among the three subtypes, C3 exhibited the highest expression of HLA molecules, which represented superior power for antigen presentation, while C1 displayed the deficient capability to present antigen (Figure 5D). Furthermore, a spectrum of immunogenicity indicators was evaluated, including neoantigen load (single nucleotide variant (SNV) and indel neoantigens), cancer/testis-antigens score (CTA Score), and genomic instability-related indicators. As illustrated in Figure 5E, C1 owned a high level of immunogenicity but insufficient antigen processing and presentation capabilities, including TCR Richness and Shannon. Meanwhile, C2 and C3 showed lower immunogenicity, and C3 showed the highest intratumor heterogeneity (ITH) specifically. Taken together, inadequate immune cell infiltration and deficient antigen procession and presentation capacity were the main immune escape mechanism for the C1 subtype. The absence of immunogenicity might be responsible for immune escape in C2. For C3, copious infiltration of immunosuppressive cells, high level of ICI expression, and ITH contributed to immune escape.

The assessment of response to immunotherapy

Immunotherapy is recommended for clinical treatment because of its good efficacy and fewer side effects. Therefore, we evaluated the benefit of immunotherapy in CRD-related subtypes to further guide the clinical application of immunotherapy approaches. Notably, C3 harbored the highest level of immune checkpoints (Figure 6A), which suggested the benefit of immune checkpoint inhibitors (ICIs). We further enrolled TIS and Submap algorithm to evaluate the efficacy of immunotherapy. As expected, C3 had the highest TIS score, indicating that it was more likely to benefit from ICI treatment (Figure 6B). In line, in two independent immunotherapy cohorts, C3 showed an expression profile that was more similar to that of patients who responded to anti-PD-1 treatment (both Bonferroni corrected and Nominal P value < 0.05) (Figure 6C). In conclusion, immunotherapy was a recommended treatment modality for patients in C3.

Potential drug development

To facilitate the personalized treatment of each subtype, we used the pRRophetic package to evaluate the half-maximal inhibitory concentration (IC50) of potential sensitive drugs for C1 and C2 subtypes based on drug sensitivity data from CTRP and PRISM databases. For common clinical chemotherapeutic agents such as docetaxel and paclitaxel, C1 patients displayed better sensitivity, suggesting an effect of chemotherapy (Figure 6D, Figure S2E). Meanwhile, C2 patients may benefit more from corticosteroids, including dexamethasone and prednisone (Figure 6E, Figure S2F).

Technologies of next-generation sequencing and machine learning have flourished, which provided a powerful basement for the development of prognostic tools. Several prognostic signatures were constructed and validated to aid clinical decisions. After a systematic search of the published literature, a total of 151 LUAD prognostic signatures based on ample biological processes and robust machine-learning algorithms were collected. In the TCGA-LUAD cohort, the risk score of each model was recalculated for all patients and appraised their efficiency via the C-index (Figure 6F). Interestingly, Chen et al. (2017) possessed the most accurate discrimination for both C1 (C-index = 0.741) and C2 (C-index = 0.698) patients [41] . Therefore, the model was suitable to predict their prognosis to optimize clinical management. Similarly, Ma et al. (2020) was the optimum signature for patients in C3 (C-index = 0.796) [42] .

Discussion

The heterogeneity of LUAD, which contributes to elusive prognosis and treatment sensitivity [43] , continues to perplex clinicians and researchers. Considering the critical role of CRD in tumorigenesis and progression [44] , it is essential for LUAD patients to accept proper stratification strategies based on CRD status, which helps with delicacy management and personalized treatment. As far as we know, deficiency has hitherto existed in CRD molecular heterogeneity research for LUAD. In this study, heterogeneous CRD subtypes were identified and systematically analyzed specific characteristics from multiple perspectives, including biological function, immune landscape, and genomic alteration. These results could improve our understanding of CRD and refine clinical management and personalized treatment.

We identified three CRD-related subtypes through the NMF algorithm. The reliability and stability of the three subtypes were validated in multiple ways. As an efficient tool, the NTP algorithm was applied to assess the stability of subtypes through the specific DEGs expression profile. Finally, similar specific gene expression profiles, proportions, and prognoses of each subtype were demonstrated in six independent GEO cohorts, suggesting the rationality of CRD-related subtypes. As displayed, C1 possessed the worst prognosis, C2 owned the best outcome, and the OS of C3 was between C1 and C2.

As described, there was distinct heterogeneity of biological functions in three subtypes. C1 was depicted by activated proliferation pathways, C2 was characterized by enrichment of lipid metabolic, whereas C3 was distinguished by prominent association to immune-related function. Besides, C1 displayed negative regulation of the circadian sleep-wake cycle. In contrast, C2 was enriched in the regulation of the circadian sleep-wake cycle. Notably, the CRD status of the three subtypes revealed similar trends with their prognosis. Besides, genomic characteristics of diverse genomic characteristics. C1 possessed the most profound genomic instability from analyses of both somatic mutation and CNV. Previous research revealed that TP53 mutations were associated with active DNA damage repair (DDR) and elevated cell proliferation levels [45,46] . Besides, an increased risk of immune escape and dismal prognosis appeared with additional TP53 mutation [47] . Correspondingly, with the highest TP53 mutation, C1 presented enrichment in proliferation pathways and a poorer prognosis. Besides, we explored the immune landscape to further depicted the heterogeneous immune microenvironment of three subtypes. As illustrated, C3 possessed conspicuous immune cell infiltration and antigen presentation ability, which was consistent with the "immune-hot" phenotype. Moreover, the assessment of immune escape mechanisms demonstrated that C1 may achieve immune escape primarily due to insufficient TCR richness and TCR Shannon diversity; C2 appears to rely on deficient immunogenicity; and C3 likely utilizes a mechanism driven by high intratumor heterogeneity.

As described above, C1 displayed negative regulation of the circadian sleep-wake cycle, suggesting CRD status. For LUAD patients, systemic and somatic disruption of circadian rhythms contribute to enhanced proliferation and metabolic dysregulation, which resulting in cancer progression and related poor prognosis [48] . Besides, recent research provided evidence that CRD promoted genomic instability in LUAD patients [49] , which was consistent with our findings in the C1 subtype. Nevertheless, the specific regulatory relationships among CRD, genomic instability and TP53 mutation still need to be further studied. Notably, CRD-related lipid metabolic dysregulation was a significant factor for tumor invasion and TME remodeling [50] . However, for C2 with relative healthy circadian rhythms, undisturbed lipid metabolism made a contribution for relatively good prognosis.

As is well acknowledged, it's crucial and essential to provide individualized treatment. Development and clinical application of tumor immunotherapy, especially immune checkpoint inhibitor (ICI) therapy, has revolutionized treatment patterns in LUAD. As humanized monoclonal antibodies for blocking immune checkpoints (such as CTLA-4 and PD-1), ICI works by restoring effective immune cell function [51,52] . C3, the "immune-hot" subtype, harbored abundant immune cell infiltration and higher expression of immune checkpoints, such as CD8+ T cell, CD4+ T cell, CD274 (PD-L1), PDCD1, and CTLA-4 molecular. Meanwhile, together with the highest TIS score and similar characteristics with patients who respond to ICIs, C3 was thought to derive potential benefits from ICIs treatment. Due to the heterogeneity of LUAD, ill-fitted chemotherapeutic caused additional side effects, which was severe challenges for LUAD patients. Our study suggested that C1 patients with significant TP53 mutation and increased proliferative activity may have better efficacy with two kinds of important chemotherapeutic drugs, docetaxel and paclitaxel. In addition, for C2 patients, glucocorticoids were recommended for their potential effect on adverse events for LUAD. Finally, we included 151 LUAD prognostic signatures to boost clinical management. The 21-gene model of Chen et al. possessed the highest accuracy for both C1 and C2 and the 16-gene model of Ma et al. exhibited the best discrimination for C3 patients, implying their excellent ability to be applied for prognostic management to individual subtypes.

Several limitations need to be acknowledged in the present study. Firstly, the clinical information was incomplete for some patients in public datasets, which contribute to potential bias. Then, we focused on inter-tumor heterogeneity from the bulk RNA sequencing data, further intra-tumor heterogeneity studies haven't been considered from the single cell level. Additional experimental and validation data for the biological pathway characterizations associated with each subtype were meaningful. At last, further clinical validation for sensitivity to immunotherapy and specific drugs predicted by machine learning algorithms was still needed.

In conclusion, we revealed three heterogeneous CRD subtypes in LUAD. Among the three subtypes, different survival times, biological features, genomic alterations, immune landscape, and treatment responses were spotted. Our work provided a promising classification platform and individualized treatment strategies, which would be helpful for clinical management. Based on RNA sequencing results, clinicians can gain insights into the characteristics of different subtypes, enabling accurate classification and targeted treatment.

Abbreviations

CI: Confidence Interva; CNV: Copy Number Variation; CRD: Circadian Rhythm Disruption; FGA: Fraction of Genome Alteration; FGG: Fraction of Genomic Gained; FGL: Fraction of Genome Lost; FPKM: Fragments Per Kilobase of Million Mapped Reads; GEO: Gene Expression Omnibus; GO: Gene Ontology; GSEA: Gene Set Enrichment Analysis; HR: Hazard Ratio; KEGG: Kyoto Encyclopedia of Genes and Genomes; LUAD: Lung Adenocarcinoma; NMF: Non-negative Matrix Factorization; TCGA: The Cancer Genome Atlas; TCR: T Cell Receptor; TP53: Tumor Protein p53; TPM: Transcripts Per Million

Supplementary Material

Supplementary methods, results, spectra, figures.

Declarations

Acknowledgements

We acknowledge assistance with the access of analytic instruments from Translational Medical Center at The First Affiliated Hospital of Zhengzhou University. We acknowledge the support from Henan Key Laboratory for Pharmacology of Liver Diseases.

Author contributions

R.W contributed study design and paper revisiting. K.K and X.S contributed project oversight. P.Y, Y.Y, and P.X contributed data analysis and visualization. X.H, Z.S, and S.G contributed paper writing. H.Z, S.Z, C.Y, T.L, and G.Z contributed paper revisiting. All authors approved this manuscript.

Ethics approval and consent to participate

Not applicable.

Funding information

Not applicable.

Competing Interests

The authors declare that they have no existing or potential commercial or financial relationships that could create a conflict of interest at the time of conducting this study.

Data Availability

The datasets generated during and/or analyzed during the current study are available from the corresponding author upon reasonable request. Publicly available datasets analyzed during the current study are available in the TCGA database and GEO database under accession codes. TCGA (https://www.cancer.gov/tcga), GSE72094 (https://www.ncbi. nlm.nih.gov/geo/query/acc.cgi?acc=GSE72094)(54), GSE68465 (https://www.ncbi. nlm.nih.gov/geo/query/acc.cgi?acc=GSE68465)(55), GSE50081 (https://www.ncbi. nlm.nih.gov/geo/query/acc.cgi?acc=GSE50081)(56), GSE42127 (https://www.ncbi. nlm.nih.gov/geo/query/acc.cgi?acc=GSE42127)(57), GSE41271 (https://www.ncbi. nlm.nih.gov/geo/query/acc.cgi?acc=GSE41271)(58), GSE31210 (https://www.ncbi. nlm.nih.gov/geo/query/acc.cgi?acc=GSE31210)(59), GSE93157 (https://www.ncbi. nlm.nih.gov/geo/query/acc.cgi?acc=GSE93157)(60), and a comprehensive immunotherapy annotations cohorts(61). Additional data related to this paper may be requested from the authors.

References