Gene expression signatures as candidate biomarkers of response to PD-1 blockade in non-small cell lung cancers

Abstract

Although anti-PD-1/PD-L1 monotherapy has achieved clinical success in non-small cell lung cancer (NSCLC), definitive predictive biomarkers remain to be elucidated. In this study, we performed whole-transcriptome sequencing of pretreatment tumor tissue samples and pretreatment and on-treatment whole blood samples (WB) samples obtained from a clinically annotated cohort of NSCLC patients (n = 40) treated with nivolumab (anti-PD-1) monotherapy. Using a single-sample gene set enrichment scoring method, we found that the tumors of responders with lung adenocarcinoma (LUAD, n = 20) are inherently immunogenic to promote antitumor immunity, whereas those with lung squamous cell carcinoma (LUSC, n = 18) have a less immunosuppressive tumor microenvironment. These findings suggested that nivolumab may function as a molecular targeted agent in LUAD and as an immunomodulating agent in LUSC. In addition, our study explains why the reliability of PD-L1 expression on tumor cells as a predictive biomarker for the response to nivolumab monotherapy is quite different between LUAD and LUSC.

About the Authors:

Tomoiki Aiba

Roles Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Visualization, Writing – original draft

* E-mail: [email protected]

Affiliation: Department of Pulmonary Medicine, Sendai Kousei Hospital, Sendai, Japan

ORCID logo https://orcid.org/0000-0001-8991-237X

Chieko Hattori

Roles Investigation, Resources

Affiliation: Department of Pulmonary Medicine, Sendai Kousei Hospital, Sendai, Japan

Jun Sugisaka

Roles Investigation

Affiliation: Department of Pulmonary Medicine, Sendai Kousei Hospital, Sendai, Japan

ORCID logo https://orcid.org/0000-0002-1000-0077

Hisashi Shimizu

Roles Investigation

Affiliation: Department of Pulmonary Medicine, Sendai Kousei Hospital, Sendai, Japan

Hirotaka Ono

Roles Investigation

Affiliation: Department of Pulmonary Medicine, Sendai Kousei Hospital, Sendai, Japan

Yutaka Domeki

Roles Investigation

Affiliation: Department of Pulmonary Medicine, Sendai Kousei Hospital, Sendai, Japan

Ryohei Saito

Roles Investigation

Affiliation: Department of Pulmonary Medicine, Sendai Kousei Hospital, Sendai, Japan

Sachiko Kawana

Roles Investigation

Affiliation: Department of Pulmonary Medicine, Sendai Kousei Hospital, Sendai, Japan

Yosuke Kawashima

Roles Investigation

Affiliation: Department of Pulmonary Medicine, Sendai Kousei Hospital, Sendai, Japan

Keisuke Terayama

Roles Investigation

Affiliation: Department of Pulmonary Medicine, Sendai Kousei Hospital, Sendai, Japan

Yukihiro Toi

Roles Investigation

Affiliation: Department of Pulmonary Medicine, Sendai Kousei Hospital, Sendai, Japan

Atsushi Nakamura

Roles Investigation

Affiliation: Department of Pulmonary Medicine, Sendai Kousei Hospital, Sendai, Japan

Shinsuke Yamanda

Roles Investigation

Affiliation: Department of Pulmonary Medicine, Sendai Kousei Hospital, Sendai, Japan

Yuichiro Kimura

Roles Investigation

Affiliation: Department of Pulmonary Medicine, Sendai Kousei Hospital, Sendai, Japan

Yutaka Suzuki

Roles Investigation, Resources

Affiliation: Department of Computational Biology and Medical Sciences, Graduate School of Frontier Sciences, The University of Tokyo, Chiba, Japan

Atsushi Niida

Roles Data curation, Resources, Supervision

Affiliation: Human Genome Center, The Institute of Medical Science, The University of Tokyo, Tokyo, Japan

ORCID logo https://orcid.org/0000-0002-6851-8004

Shunichi Sugawara

Roles Funding acquisition, Investigation, Supervision

Affiliation: Department of Pulmonary Medicine, Sendai Kousei Hospital, Sendai, Japan

ORCID logo https://orcid.org/0000-0002-3427-4558

Introduction

For the last several decades, there have been remarkable advances in cancer immunology and cancer immunotherapy. The clinical efficacy of immuno-oncology (IO) agents that inhibit the programmed death 1 (PD-1)-programmed death ligand 1 (PD-L1) signaling pathway in advanced NSCLC has been described. In particular, nivolumab, a fully human IgG4 anti-PD-1 monoclonal antibody, is the first approved IO agent for use in patients with previously treated advanced NSCLC. Unfortunately, nivolumab monotherapy is not effective in all patients with advanced NSCLC, with an objective response rate (ORR) of no higher than 20% [1, 2].

Preclinical and clinical studies have revealed that various factors can affect the clinical outcomes of anti-PD-1/PD-L1 monotherapy in patients with advanced NSCLC; these factors include PD-L1 expression, the presence of tumor-infiltrating lymphocytes (TILs), the tumor mutation burden (TMB), human leukocyte antigen class I (HLA-I) genotype, T-cell repertoire diversity, the gene expression profile and the gut microbiota [3–7]. Among these factors, PD-L1 expression, which is determined by immunohistochemical assays, is the only biomarker currently approved as a companion or complementary diagnostic biomarker for the response to anti-PD-1/PD-L1 monotherapy. In general, PD-L1 expression on tumor cells is positively correlated with the clinical response to anti-PD-1/PD-L1 monotherapy in patients with nonsquamous NSCLC [2, 8]. Importantly, however, some patients with a positive PD-L1 expression rate of less than 1% also benefit from anti-PD-1/PD-L1 monotherapy, with an ORR of approximately 10% [9]. Unlike in patients with nonsquamous NSCLC, PD-L1 expression has no predictive value in patients with squamous NSCLC (1). PD-L1 expression alone remains an imperfect biomarker, and none of the abovementioned factors has yet proven to be sufficiently robust for clinical use. Therefore, extensive efforts have been devoted to exploring predictive biomarkers for the response to anti-PD-1/PD-L1 monotherapy.

The recent development of next-generation sequencing (NGS) technology and computational genomic tools has enabled deep analysis of the so-called omics data, including genomic, transcriptomic, proteomic, and epigenomic data [10]. To date, biomarker studies in NSCLC have mainly focused on the histological or genomic analyses, and there have been very few transcriptomic analyses. In metastatic melanoma and other types of cancer, multiple comprehensive transcriptomic studies have identified immune-related gene expression signatures that are positively or negatively associated with the clinical response to anti-PD-1 and anti-CTLA4 monotherapy [11, 12]. We thus assumed that by combining gene expression signatures obtained from pretreatment and on-treatment transcriptomic profiles with clinical data of patients, we could identify a novel promising biomarker to predict the response to anti-PD-1/PD-L1 monotherapy in NSCLC patients.

The two main histological subtypes of NSCLC are adenocarcinoma and squamous cell carcinoma. In biomarker studies, these subtypes have often been combined and analyzed in the same way; however, each exhibits distinct mutational and genomic profiles. In this cohort, we analyzed the patient’s transcriptomic features separately to identify histology-specific gene expression signatures that are associated with the clinical response to nivolumab monotherapy. Characterization of these signatures will help us to decipher the complexity of tumor-immune interactions and deepen the understanding of the tumor microenvironment (TME) that favors a better clinical response to nivolumab monotherapy.

Materials and methods

Ethical statement

All clinical data and patient samples were collected following approval by the Sendai Kousei Hospital Institutional Review Board (IRB) (IRB number: 29–4). The study period is from 18 Jul 2017 to 31 Dec 2020. In all cases, written informed consent was obtained from the patients.

Patient characteristics and sample collection

A total of 40 patients with advanced NSCLC were enrolled in this cohort. All enrolled patients were administered at least one dose of nivolumab. Tumor tissue samples were collected before the first dose of nivolumab (pretreatment tumor tissue). Immediately after the biopsy procedure (three tumor tissue samples of approximately 1.5−3.0 mm in diameter per patient), the obtained tumor tissues were suspended in RNAlater RNA Stabilization Reagent (QIAGEN, Hilden, Germany) and stored at −80°C for batched RNA extraction. Whole blood (WB) samples were collected before the first dose and after the fourth or fifth dose of nivolumab (pretreatment WB and on-treatment WB, respectively). To obtain whole blood samples, 2.5 mL of blood was collected into PAXgene Blood RNA Tubes (BD, Franklin Lakes, NJ) and stored at −80°C for batched RNA extraction.

For 28 of the 40 patients (15 LUAD and 13 LUSC), pretreatment tumor tissue samples were available and used in subsequent analyses; the other 12 patients were excluded because of an inadequate quantity of tumor cells. Pretreatment WB was obtained from 39 patients (20 LUAD and 17 LUSC); for one patient, isolation of WB was unsuccessful due to a blood sampling error. On-treatment WB was obtained from 32 patients (15 LUAD and 15 LUSC); for the other 8 patients, blood sampling was not performed due to death, early disease progression or early treatment discontinuation. In this cohort, we employed progression-free survival (PFS) time as an outcome measure of treatment efficacy. The PFS time was defined as the length of time from the start of nivolumab monotherapy until progression. Patients with a PFS time equal to or more than 6 months were defined as responders; those with a PFS time less than 6 months, non-responders.

The PD-L1 tumor proportion score (TPS) was determined by a commercial PD-L1 IHC assay with PD-L1 IHC 22C3 pharmDx (Dako, Carpinteria, CA). During the enrollment and follow-up period of this cohort, four genetic screens of driver mutations, including EGFR and BRAF mutations and ALK and ROS1 fusions, had been approved and were commercially available for advanced NSCLC in Japan. Other genetic aberrations, such as KRAS mutations, were screened in the LC-SCRUM-Asia (formerly LC-SCRUM-Japan) consortium [13], which employed an amplicon-based next-generation sequencing (NGS) panel, Oncomine Comprehensive Assay (Thermo Fisher Scientific, Waltham, MA). Tumor responses to nivolumab monotherapy were assessed by the investigators according to RECIST (Response Evaluation Criteria in Solid Tumors) version 1.1.

RNA extraction and whole transcriptome sequencing (RNA-seq)

Total RNA was extracted from tumor tissue and WB samples using an RNeasy Mini Kit (QIAGEN) and a PAXgene Blood RNA Kit (QIAGEN), respectively, according to the manufacturer’s protocol. Homogenization of tumor tissue samples was carried out using a QIAshredder homogenizer (QIAGEN). DNase I was used during processing with the manufacturer’s protocol to minimize DNA contamination.

The RNA-seq libraries were prepared from the total RNA extracts using an Illumina TruSeq Stranded mRNA Library Kit (Illumina, San Diego, CA), following the manufacturer’s protocol. Paired-end 100 bp sequencing of these libraries was performed on an Illumina HiSeq3000 platform. We utilized HISAT2 [14] to align the RNA-seq reads to the human genome assembly GRCh38. The raw read counts were generated with featureCounts in the Rsubreads Bioconductor package (version 2.4.2) [15] and normalized by conversion to TPM (transcripts per million) [16].

Identification of differentially expressed genes (DEGs)

The RNA-seq raw read counts were employed for differential expression analysis. As a preprocessing step, filtering was applied to exclude genes with low expression, retaining genes with a raw read counts > 1 in at least n samples (where n is the number of samples in each analysis). We fitted a generalized linear model to the preprocessed count data, which is generally assumed to follow a negative binomial distribution, using the DESeq2 Bioconductor package (version 1.29.6) [17]. The statistical significance of differences between responders and nonresponders was assessed with the Wald test. DEGs were defined as genes with |Log2(fold change)| ≥ 1 and an adjusted p-value < 0.1. MA and volcano plots were generated using the ggplot2 R package (version 3.3.2). Hierarchical clustering and heatmap representation of the DEGs were implemented using the ComplexHeatmap Bioconductor package (version 2.5.3) [18]. DESeq2-normalized counts were converted to log10 values, followed by normalization to the z-score values for all genes to reduce expression variance for the genes expressed at different levels.

Fast gene set enrichment analysis (FGSEA)

In classical GSEA [19, 20], correction for multiple hypothesis testing is performed using permutation tests, where independent random gene sets are generated for each permutation and each gene set. Accordingly, standard GSEA is relatively slow because of the huge computational burden imposed by the permutation tests. FGSEA with the fgsea Bioconductor package (version 1.15.0) reuses sampling for different query gene sets to reduce the computational burden, which enables quick and accurate performance of enrichment analysis [21]. As a preprocessing step, we first generated a preranked gene list from the DESeq2 results. This preranked gene list was used for FGSEA with 7,573 gene sets in the Gene Ontology (GO) biological process ontology from MSigDB v7.2 (http://www.gsea-msigdb.org/gsea/msigdb/). The minimal and maximal thresholds of the gene sets were set to 10 and 500 genes, respectively. Statistical significance was calculated with 1,000 permutation tests. Significantly enriched gene sets were defined as those with an adjusted p-value < 0.1. Lollipop plots of the normalized enrichment scores (NESs) obtained from FGSEA were constructed using ggplot2. Running enrichment scores were plotted using the fgsea package. The enrichment maps were generated using Cytoscape App (version 3.8.2: https://cytoscape.org).

Singscore

The TPM normalized counts were used in single-sample gene set scoring with the singscore Bioconductor package (version 1.9.0) [22, 23]. Genes with low counts were filtered based on the CPM (counts per million) value in the edgeR Bioconductor package (version 3.31.4) [24] to avoid rank duplication: we retained genes with a CPM > 1 across more than 50% of the samples. Genes were ranked based on count in increasing order. For each gene set of interest, the mean rank was calculated and normalized to the theoretical minimum and maximum values, centered on zero and then summed to provide a single-sample enrichment score, which ranged between −1 and 1. As with FGSEA, single-sample enrichment scores were generated using the singscore method for gene sets in the GO biological process ontology from MSigDB.

The distributions of the single-sample enrichment scores across all the gene sets were visualized by stacked histograms. The single-sample enrichment scores of the gene sets of interest were shown by scatter plots or box plots and statistically evaluated with the Wilcoxon rank-sum test. In box plots, the lower and upper box hinges indicate the 25th and 75th percentiles, respectively; the central bold line indicates the median; and the whiskers extend to the largest and smallest scores within no more than 1.5× the interquartile range. All the statistical tests were performed using the R program (version 4.0.2; https://www.r-project.org/). Stacked histograms and box plots were generated with ggplot2.

Correlation analysis and regression model fitting using cubic regression splines

Correlation coefficients and p-values between PFS time and the single-sample enrichment scores of the gene sets of interest were calculated by Spearman rank correlation analysis using R. The correlation matrix obtained was visualized with ggplot2.

Nonparametric regression models using cubic regression splines with the R and mgcv packages (version 1.8.31) were fitted to the calculated relationship between PFS time and the single-sample enrichment scores of the gene sets of interest, from which we estimated the predicted PFS time from the scores [25]. Assuming that Log[observed PFS] is normally distributed, we employed the Gaussian family and identity link function. In addition, we used REML (restricted maximum likelihood) as a smoothing parameter estimation method and performed model selection based on AICc (second-order Akaike information criterion) with the MuMIn package (version 1.43.17) [26]. The regression splines were plotted using ggplot2.

Survival analysis

Bubble plots displaying the relationships between PFS time and two of the gene sets of interest were generated using ggplot2. All Kaplan-Meier curves were visualized with the survminer package (version 0.4.9). With the survival package (version 3.2.3), Kaplan-Meier estimates of PFS time for two independent groups were assessed using the two-sided log-rank test, and those for four independent groups were assessed using the two-sided Wald test based on the multivariable Cox proportional hazards model. When comparing the PFS times of patients with high scores to those of patients with low scores, we set the threshold as the median value of the score. Both the waterfall plots describing changes in tumor size and swimmer plots showing patient responses were generated using ggplot2. Cox proportional hazard models were built using the “coxph” function from the survival package and visualized as forest plots using ggplot2.

Results and discussion

From clinically annotated NSCLC patients treated with nivolumab monotherapy in the second- or later-line settings, we prospectively collected tumor tissue and whole blood (WB) samples before the first dose of nivolumab (pretreatment tumor tissues and WB), and WB samples after the fourth or fifth dose of nivolumab (on-treatment WB) (S1 Fig). Patient characteristics and clinical courses are summarized in S2 and S3 Figs, Table 1 and S1 Table. All tumor tissue and WB samples obtained were subjected to whole transcriptome sequencing (RNA-seq). We extracted transcriptomic datasets of LUAD and LUSC from the results, analyzed each histological subtype separately, and explored whether transcriptomic features could differentiate between responders and nonresponders to nivolumab monotherapy.

[Figure omitted. See PDF.]

Table 1. Baseline characteristics of patients in this cohort.

https://doi.org/10.1371/journal.pone.0260500.t001

A previous study has provided a metric for immune cytolytic activity based on gene expression in tumors, where immune cytolytic activity was estimated by the average expression level of CD8A, CD8B, GZMA, GZMB and PRF [27]. Using this metric, we assessed immune cytolytic activity in tumor tissues from patients in our cohort. As the result, we found no significant association between the immune cytolytic activity and clinical outcomes (S4 Fig).

A cohort of patients with LUAD

In the LUAD cohort (n = 20), differential expression analysis between responders and nonresponders was performed using the RNA-seq datasets of pretreatment tumor tissues (n = 15), pretreatment WB (n = 20) and on-treatment WB (n = 15). A total of 15, 68 and 160 differentially expressed genes (DEGs) were identified in pretreatment tumor tissues, pretreatment WB and on-treatment WB, respectively. Of the 15 DEGs in pretreatment tumor tissues, six genes were upregulated and nine genes were downregulated in responders (n = 4) compared to nonresponders (n = 11) (S5A and S5B Fig and S2 Table). Of the 68 DEGs in pretreatment WB, 27 genes were upregulated and 41 genes were down-regulated in responders (n = 5) compared to nonresponders (n = 15) (S5C and S5D Fig and S3 Table). In addition, of the 160 DEGs in on-treatment WB, 20 genes were upregulated and 140 genes were down-regulated in responders (n = 4) compared to nonresponders (n = 11) (S5E and S5F Fig and S4 Table). Heatmaps of DEGs identified in the three datasets demonstrate the relationships between the expression levels of DEGs and clinical factors, including progression-free survival (PFS), PD-L1 tumor proportion score (TPS), Brinkman index and driver mutation status (S6 Fig).

To identify biological processes associated with clinical outcomes in LUAD patients treated with nivolumab monotherapy, we performed gene set enrichment analysis (GSEA) (Fig 1, S7 Fig and S5–S7 Tables). In pretreatment tumor tissues of responders, gene sets related to B cell function (e.g., ‘B CELL-MEDIATED IMMUNITY’ [padj = 8.424 × 10−20, NES = 2.738], ‘B CELL RECEPTOR SIGNALING PATHWAY’ [padj = 1.692 × 10−11, NES = 2.625], and ‘POSITIVE REGULATION OF B CELL ACTIVATION’ [padj = 3.990 × 10−7, NES = 2.225]) and humoral immunity (e.g., ‘HUMORAL IMMUNE RESPONSE MEDIATED BY CIRCULATING IMMUNOGLOBULIN’ [padj = 1.534 × 10−20, NES = 3.003], ‘COMPLEMENT ACTIVATION’ [padj = 1.317 × 10−19, NES = 2.896], REGULATION OF HUMORAL IMMUNE RESPONSE’ [padj = 1.234 × 10−15, NES = 2.853], ‘HUMORAL IMMUNE RESPONSE’ [padj = 1.648 × 10−23, NES = 2.664], and ‘IMMUNOGLOBULIN PRODUCTION’ [padj = 4.344 × 10−16, NES = 2.626]) were significantly enriched (Fig 1A and S5 Table). Recently, three studies demonstrated that B cells and tertiary lymphoid structures (TLSs) in the TME are associated with favorable clinical outcomes to immunotherapy in patients with melanoma, renal cell carcinoma and sarcoma [28–30]. It has also been reported that the presence of B cells and TLSs could be a potential prognostic marker for NSCLC [31, 32]. Thus, these findings indicate that the involvement of B cells and TLSs in the TME seems to facilitate a better response to immunotherapy in LUAD.

[Figure omitted. See PDF.]

Fig 1. Classical gene set enrichment analysis (GSEA) in LUAD.

A–F, Lollipop plots depicting the GSEA results in the following samples: pretreatment tumor tissues (Pre-tissue) of responders (A) and nonresponders (B), pretreatment WB (Pre-WB) of responders (C) and nonresponders (D), and on-treatment WB (On-WB) of responders (E) and nonresponders (F). The X-axes show the normalized enrichment score (NES), and the Y-axes show gene sets ranked among the top 30 enriched gene sets with adjusted p-value < 0.10 (in descending order of NES). The dot size is proportional to the size of the corresponding gene set. The dot color indicates the adjusted p-value.

https://doi.org/10.1371/journal.pone.0260500.g001

Although gene sets enriched in nonresponders were extremely similar between pretreatment tumor tissues and pretreatment WB, gene sets enriched in responders were quite different between these sources (Fig 1A–1D and S5 and S6 Tables). Notably, in pretreatment WB, gene sets significantly enriched in responders were those related to type I/II interferon (IFN) signaling (e.g., ‘RESPONSE TO TYPE I INTERFERON’ [padj = 2.383 × 10−14, NES = 3.049], ‘INTERFERON-GAMMA-MEDIATED SIGNALING PATHWAY’ [padj = 6.655 × 10−14, NES = 2.946], ‘RESPONSE TO INTERFERON-GAMMA’ [padj = 2.968 × 10−18, NES = 2.939], ‘RESPONSE TO INTERFERON-BETA’ [padj = 2.708 × 10−7, NES = 2.720], and ‘CELLULAR RESPONSE TO INTERFERON-BETA’ [padj = 3.193 × 10−5, NES = 2.506]) and the process of antigen processing and presentation (e.g., ‘ANTIGEN PROCESSING AND PRESENTATION OF EXOGENOUS PEPTIDE ANTIGEN VIA MHC CLASS I’ [padj = 1.604 × 10−10, NES = 2.755], ‘ANTIGEN PROCESSING AND PRESENTATION’ [padj = 5.540 × 10−17, NES = 2.741], ‘ANTIGEN PROCESSING AND PRESENTATION OF PEPTIDE ANTIGEN’ [padj = 1.532 × 10−15, NES = 2.728], and ‘ANTIGEN PROCESSING AND PRESENTATION OF PEPTIDE ANTIGEN VIA MHC CLASS I’ [padj = 4.537 × 10−11, NES = 2.706]) (Fig 1C and S6 Table). Gene sets related to host defense against viral infection (e.g., ‘DEFENSE RESPONSE TO VIRUS’ [padj = 1.799 × 10−17, NES = 2.750], ‘RESPONSE TO VIRUS’ [padj = 3.428 × 10−18, NES = 2.583], ‘NEGATIVE REGULATION OF VIRAL LIFE CYCLE’ [padj = 1.072 × 10−7, NES = 2.578], ‘NEGATIVE REGULATION OF VIRAL GENOME REPLICATION’ [padj = 4.863 × 10−7, NES = 2.546], ‘REGULATION OF VIRAL LIFE CYCLE’ [padj = 9.7864 × 10−10, NES = 2.464], and ‘NEGATIVE REGULATION OF VIRAL PROCESS’ [padj = 1.642 × 10−7, NES = 2.392]), where type I IFN signaling plays a key role, were also enriched in responders (Fig 1C, S7 Fig and S6 Table). Type I IFNs, such as IFN-α and IFN-β, indirectly elicit antitumor immune responses in the TME by stimulating the maturation of dendritic cells, increasing the expression of perforin and granzymes in cytotoxic T cells, promoting the survival of memory T cells, and inactivating the suppressive function of regulatory T (Treg) cells. In addition, type I IFNs enhance antitumor immunity directly through the inhibition of tumor cell proliferation and the acceleration of senescence and apoptosis [33, 34]. Type II IFN (IFN-γ) also supports antitumor immunity by augmenting the function of tumor-infiltrating immune cells, inactivating suppressive Treg cells, and modulating stromal cell function to alter metabolism and inhibit angiogenesis [35]. Both type I and type II IFN signaling can promote the process of antigen processing and presentation, which means that gene sets related to antigen processing and presentation (hereinafter referred to as ‘APP signatures’) are closely linked to those related to IFN signaling (hereinafter referred to as ‘IFN signatures’). Moreover, type II IFN can induce the expression of PD-L1 on tumor cells and tumor-associated macrophages (TAMs) in the TME, which is an ingenious strategy that tumor cells employ to evade the host immune system [35]. In contrast to pretreatment WB of responders, pretreatment tumor tissues of responders did not show strong enrichment of IFN and APP signatures (Fig 1A and S5 Table). We presume that this difference may arise from the presence of immunosuppressive conditions in the TME mainly due to PD-1/PD-L1 signaling; that is, tumor-infiltrating immune cells are inactivated in the immunosuppressive TME of responders, while WB of responders is entirely unaffected by the TME.

In on-treatment WB, APP signatures, especially gene sets related to the process of antigen processing and presentation via MHC class I molecules (e.g., ‘ANTIGEN PROCESSING AND PRESENTATION OF EXOGENOUS PEPTIDE ANTIGEN VIA MHC CLASS I’ [padj = 3.337 × 10−7, NES = 2.565], and ‘ANTIGEN PROCESSING AND PRESENTATION OF PEPTIDE ANTIGEN VIA MHC CLASS I’ [padj = 3.795 × 10−7, NES = 2.453]), were significantly enriched in responders (Fig 1E and S7 Table). This enrichment of APP signatures in responders seems to reflect a durable antitumor immune response elicited by nivolumab monotherapy. Interestingly, gene sets related to mitochondrial metabolism (e.g., ‘ATP SYNTHESIS COUPLED ELECTRON TRANSPORT’ [padj = 4.044 × 10−6, NES = 2.442], ‘OXIDATIVE PHOSPHORYLATION’ [padj = 1.482 × 10−6, NES = 2.302], ‘RESPIRATORY ELECTRON TRANSPORT CHAIN’ [padj = 3.639 × 10−5, NES = 2.219], ‘ATP SYNTHESIS COUPLED PROTON TRANSPORT’ [padj = 5.748 × 10−3, NES = 2.216], and ‘MITOCHONDRIAL ELECTRON TRANSPORT, NADH TO UBIQUINONE’ [padj = 4.266 × 10−3, NES = 2.042]) were significantly enriched in responders (Fig 1E and S7 Table). Mitochondrial metabolism supports proinflammatory signaling; in addition, the proinflammatory milieu can reprogram mitochondrial metabolism. The electron transport chain produces adenosine triphosphate (ATP) and reactive oxygen species (ROS) via coupling with oxidative phosphorylation (OXPHOS), which can drive the differentiation and activation of T cells [36, 37]. These findings indicate that mitochondrial metabolism may be persistently enhanced to effectively support antitumor immunity in on-treatment WB of responders.

To verify the association between PFS and the gene signatures enriched in responders irrespective of phenotypic information (e.g., responders vs. nonresponders), we used an unsupervised single-sample gene set enrichment scoring approach. Among the so-called unsupervised, nonparametric methods, ssGSEA (single-sample gene set enrichment analysis) [38] and GSVA (gene set variation analysis) [39] have been described as the most common methods. Both ssGSEA and GSVA, however, need expression data and phenotypic information for all samples—the former to normalize enrichment scores and the latter to conduct kernel density estimation to approximate the cumulative density function. Thus, we employed an alternative unsupervised method, singscore, which is a truly single-sample gene set enrichment scoring method [22].

Using the singscore method, we computed and evaluated single-sample gene set enrichment scores of the enriched gene sets identified through the supervised GSEA method above (hereinafter referred to as ‘GSEA gene sets’; Fig 2A, S8–S10 Figs and S8–S10 Tables). The single-sample gene set enrichment scores of GSEA gene sets in pretreatment tumor tissues were not significantly different between responders and nonresponders (S8 Fig). In pretreatment WB of responders, the single-sample enrichment scores showed that IFN and APP signatures were strongly enriched with high reproducibility, indicating that they were robustly enriched regardless of whether the enrichment analysis was supervised or unsupervised (Fig 2A). The single-sample enrichment scores of APP signatures were also significantly higher in on-treatment WB of responders (S10 Fig). In nonresponders, the enrichment scores of type I IFN signaling in on-treatment WB were significantly increased compared with pretreatment WB (‘IFN_I’, p = 0.0128; ‘IFNB1’, p = 0.0108; ‘IFNB2’, p = 0.0025) (Fig 2B). In responders, by contrast, the enrichment scores of type I IFN signaling in on-treatment WB showed no significant differences compared with pretreatment WB (Fig 2B). It has been reported that sustained type I IFN signaling contributes to resistance to anti-PD-1 monotherapy [34, 40]. Hence, these findings suggest that nivolumab-induced activation of type I IFN signaling may be a predictive biomarker for worse clinical outcomes in LUAD patients treated with nivolumab monotherapy.

[Figure omitted. See PDF.]

Fig 2. Single-sample enrichment analysis (singscore) in LUAD.

A, For representative GSEA gene sets that were enriched in pretreatment WB of responders, the single-sample gene set enrichment scores and dispersions were calculated using singscore and visualized on scatter plots. Red circles denote responders (n = 5); cyan circles, nonresponders (n = 15). The enrichment scores were analyzed using the Wilcoxon rank sum test. All calculated p-values are shown on the plots (*p < 0.05, **p < 0.01, and ***p < 0.001). B, Box plots indicating the single-sample enrichment scores of IFN and APP signatures calculated from four different groups with LUAD: pretreatment WB of nonresponders (Pre-Tx/NR, n = 15), on-treatment WB of nonresponders (On-Tx/NR, n = 11), pretreatment WB of responders (Pre-Tx/R, n = 5) and on-treatment WB of responders (On-Tx/R, n = 4). Red dots denote responders; cyan dots, nonresponders. The differences in the single-sample enrichment scores between the groups (i.e., ‘Pre-Tx/NR vs. On-Tx/NR’, ‘Pre-Tx/NR vs. On-Tx/R’, ‘Pre-Tx/R vs. On-Tx/R’ and ‘On-Tx/NR vs. On-Tx/R’) were evaluated by the Wilcoxon rank sum test. Only significant p-values (p < 0.05) are shown on the plots (*p < 0.05, **p < 0.01, and ***p < 0.001).

https://doi.org/10.1371/journal.pone.0260500.g002

Moreover, we noted that the enrichment scores of APP signatures exhibited no significant differences between pretreatment and on-treatment WB. As neoantigens are generated from mutations, the higher the TMB, the greater the chance that some of the neoantigens presented by MHC proteins will be immunogenic and hence enable the induction of anti-tumor immune response. In fact, accumulating evidence has indicated that high TMB is correlated with better clinical outcomes in NSCLC patients with anti-PD-1/PD-L1 therapy [41]. Additionally, several studies have reported that tumor mutation burden (TMB) in blood (or circulating tumor DNA) correlates with TMB in tumor tissue, and that high TMB in blood may serve as a potential biomarker of clinical benefit in NSCLC patients with anti-PD-1/PD-L1 therapy [4, 42]. Given that TMB in blood can be a surrogate for TMB in tumor tissue, it is rational to suggest that TMB in blood reflect the immunogenicity of tumor cells. Therefore, the comparable levels of APP signatures before and after nivolumab monotherapy indicate that nivolumab monotherapy has no significant impact on the immunogenicity of the tumor itself. Based on these findings, we presume that the tumors of responders are inherently sufficiently immunogenic to effectively elicit antigen processing and presentation for antitumor immunity. Clearly, it seems that elevated IFN signaling reveals the magnitude of antitumor immunity, which in turn upregulates PD-L1 expression on tumor cells to induce PD-1/PD-L1 signaling-mediated immune evasion. Thus, the PD-1/PD-L1 immune evasion axis can be one of the primary targets of nivolumab monotherapy in LUAD patients.

For each single-sample enrichment score, we next calculated the Spearman correlation coefficient (ρ) to estimate the correlation between the PFS time and the enrichment score (Fig 3A and S11A Fig). Consequently, we observed that both IFN and APP signatures were significantly correlated with PFS time (e.g., IFN_I, ρ = 0.590, p = 0.0061; IFNB1, ρ = 0.546, p = 0.0127; IFNG1, ρ = 0.526, p = 0.0173; APP1, ρ = 0.511, p = 0.0214). Gene sets related to host defense against viral infection (hereinafter, referred to as ‘VIRUS signatures’), which have extremely strong correlations with IFN signatures, also exhibited significant correlations with PFS time (e.g., VIRUS1, ρ = 0.651, p = 0.0019; VIRUS2, ρ = 0.624, p = 0.0033). Nonparametric regression models using cubic regression splines indicated that the enrichment scores of a single gene set in the IFN and APP signatures in pretreatment WB could be candidate biomarkers to predict PFS time (e.g., IFNB1, R-squared [R-sq] = 0.5861, AICc = 39.0470, p = 0.0007; APP1, R-sq = 0.4353, AICc = 41.9426, p = 0.0010) (Fig 3B and S11B Fig). When the enrichment scores for two gene sets in the IFN and APP signatures were combined, LUAD patients with high values for both enrichment scores tended to have longer PFS time (Fig 4A and S12 Fig). By stratifying the LUAD patients into subsets with high and low enrichment scores (stratification by the median), we found that there were quite large differences in PFS time between the patients with high enrichment of signatures and those with low enrichment of signatures (Fig 4B and S13 Fig). Multivariate Cox regression analysis identified that many combinations of two of the IFN and APP signatures remained to be independent predictive factors of PFS time (Fig 4C and S14 Fig). Collectively, these findings suggest that IFN and APP signatures in pretreatment WB may have potential predictive performance in LUAD patients treated with nivolumab monotherapy.

[Figure omitted. See PDF.]

Fig 3. IFN and APP signatures predict the nivolumab response in LUAD.

A, Spearman correlation matrix between PFS and single-sample gene set enrichment scores of IFN and APP signatures in pretreatment WB. The upper triangular region shows the values of Spearman’s ρ correlation coefficients (significant correlations are in bold; *p < 0.05, **p < 0.01, and ***p < 0.001). In the lower triangular region, positive correlations are visualized in red and negative correlations in blue. The color intensity is proportional to the value of Spearman’s ρ, and the size of the circle to the p-value. B, Scatter plots showing the relationships between PFS and the enrichment score for a single gene set in the IFN and APP signatures in pretreatment WB, with a fitted line representing the regression model using a cubic spline and 95% confidence interval. The accuracy of the fit was assessed by calculating the adjusted R-squared (R-sq) and p-values (*p < 0.05, **p < 0.01, and ***p < 0.001).

https://doi.org/10.1371/journal.pone.0260500.g003

[Figure omitted. See PDF.]

Fig 4. Two gene sets in the IFN and APP signatures as candidate biomarkers for the nivolumab response in LUAD.

A, Bubble plots showing the relationships between PFS and the enrichment scores for two gene sets in the IFN and APP signatures in pretreatment WB. Each bubble represents a patient, and the size of each bubble is proportional to the PFS time. On a gradient color scale based on the PFS time, bubbles representing responders were assigned colors ranging from white to dark red; nonresponders, ranging from white to lavender. B, Kaplan-Meier PFS curves for patients stratified by the enrichment score (‘Low’ vs. ‘High’) for two gene sets in the IFN and APP signatures. Patients with both scores above the median were defined as ‘High’ and the others as ‘Low’. The p-values were calculated by the two-sided log-rank test (*p < 0.05, **p < 0.01, and ***p < 0.001). C, Forest plot showing multivariate Cox regression analysis for potential factors associated with prolonged or shortened PFS time. Squares represent estimated hazard ratios and whiskers represent the 95% confidence intervals. Hazard ratios less than 1 indicate improved PFS time (*p < 0.05).

https://doi.org/10.1371/journal.pone.0260500.g004

A cohort of patients with LUSC

In the LUSC cohort (n = 18), differential expression analysis between responders and nonresponders was performed using the RNA-seq datasets of pretreatment tumor tissues (n = 13), pretreatment WB (n = 17) and on-treatment WB (n = 15). A total of 424 DEGs were identified in pretreatment tumor tissues. Of these 424 DEGs, 36 genes were upregulated and 388 genes were downregulated in responders (n = 3) compared to nonresponders (n = 10) (S15A and S15B Fig and S11 Table). In marked contrast to LUAD, LUSC had many more DEGs in tumor tissues than in WB. In fact, only three DEGs were identified in on-treatment WB, among which two genes were upregulated and one gene was down-regulated in responders (n = 2) compared to nonresponders (n = 13) (S15E and S15F Fig and S12 Table). We found no DEGs in pretreatment WB (S15C and S15D Fig). A heatmap of DEGs obtained from pretreatment tumor tissues shows the relations between the expression levels of DEGs and clinical factors such as PFS, PD-L1 TPS and Brinkman index (S16 Fig).

Surprisingly, the GSEA results indicated that the gene sets enriched in responders or nonresponders were fairly similar between pretreatment and on-treatment WB in LUSC (Fig 5C–5F, S17 Fig and S14 and S15 Tables). In WB of patients with LUSC, we identified no or few DEGs between responders and nonresponders and almost no changes in the enrichment pattern across the gene sets between the pre- and post-nivolumab monotherapy settings. Nivolumab monotherapy seemed to have only limited impacts on WB, indicating that a transcriptome analysis of WB may be unable to elucidate systemic effects and predict the clinical response to nivolumab monotherapy in LUSC patients. We therefore focused on pretreatment tumor tissues of LUSC for the following analysis (Fig 5A, 5B, S18 Fig and S14 Table).

[Figure omitted. See PDF.]

Fig 5. Classical GSEA in LUSC.

A–F, Lollipop plots depicting the GSEA results in the following samples: pretreatment tumor tissues (Pre-tissue) of responders (A) and nonresponders (B), pretreatment WB (Pre-WB) of responders (C) and nonresponders (D), and on-treatment WB (On-WB) of responders (E) and nonresponders (F). The X-axes show the normalized enrichment score (NES); the Y-axes, gene sets ranked among the top 30 enriched gene sets with adjusted p-value < 0.10 (in descending order of NES). The dot size is proportional to the size of the corresponding gene set. The dot color indicates the adjusted p-value.

https://doi.org/10.1371/journal.pone.0260500.g005

In pretreatment tumor tissues, gene sets related to mitochondrial metabolism (e.g., ‘MITOCHONDRIAL RESPIRATORY CHAIN COMPLEX ASSEMBLY’ [padj = 4.225 × 10−14, NES = 2.919], ‘ATP SYNTHESIS COUPLED ELECTRON TRANSPORT’ [padj = 1.102 × 10−12, NES = 2.852], ‘NADH DEHYDROGENASE COMPLEX ASSEMBLY’ [padj = 2.407 × 10−10, NES = 2.850], ‘OXIDATIVE PHOSPHORYLATION’ [padj = 2.145 × 10−15, NES = 2.838], ‘MITOCHONDRIAL GENE EXPRESSION’ [padj = 5.865 × 10−17, NES = 2.803], ‘INNER MITOCHONDRIAL MEMBRANE ORGANIZATION’ [padj = 1.819 × 10−9, NES = 2.792], ‘MITOCHONDRIAL TRANSLATION’ [padj = 4.873 × 10−14, NES = 2.750], ‘CRISTAE FORMATION’ [padj = 6.394 × 10−8, NES = 2.711], ‘RESPIRATORY ELECTRON TRANSPORT CHAIN’ [padj = 5.114 × 10−11, NES = 2.694], and ‘MITOCHONDRIAL ELECTRON TRANSPORT, NADH TO UBIQUINONE’ [padj = 2.189 × 10−7, NES = 2.596]) were significantly enriched in responders (Fig 5A, S18A Fig and S13 Table). Additionally, gene sets related to organization of the tumor microenvironment (e.g., ‘EXTRACELLULAR STRUCTURE ORGANIZATION’ [padj = 9.555 × 10−31, NES = −2.611], ‘COLLAGEN FIBRIL ORGANIZATION’ [padj = 2.039 × 10−8, NES = −2.508], ‘COLLAGEN CATABOLIC PROCESS’ [padj = 4.344 × 10−6, NES = −2.291], ‘REGULATION OF CELLULAR EXTRAVASATION’ [padj = 1.381 × 10−5, NES = −2.261], ‘COLLAGEN METABOLIC PROCESS’ [padj = 1.235 × 10−7, NES = −2.229], ‘CELLULAR EXTRAVASATION’ [padj = 2.989 × 10−6, NES = −2.226], ‘MUCOPOLYSACCHARIDE METABOLIC PROCESS’ [padj = 9.872 × 10−8, NES = −2.221], ‘POSITIVE REGULATION OF CELLULAR EXTRAVASATION’ [padj = 3.971 × 10−5, NES = −2.219], ‘EXTRACELLULAR MATRIX DISASSEMBLY’ [padj = 1.433 × 10−6, NES = −2.199], ‘NEGATIVE REGULATION OF HOMOTYPIC CELL-CELL ADHESION’ [padj = 1.000 × 10−4, NES = −2.199], and ‘TISSUE REMODELING’ [padj = 8.091 × 10−8, NES = −2.114]) were significantly enriched in nonresponders (Fig 5B, S18B Fig and S13 Table) [43–45].

Among these gene sets, singscore reproducibly identified three gene sets related to mitochondrial functions, namely, ‘MITOCHONDRIAL GENE EXPRESSION’, ‘INNER MITOCHONDRIAL MEMBRANE ORGANIZATION’ and ‘CRISTAE FORMATION’ (hereinafter referred to as ‘mitochondrial signatures’), as significantly enriched in responders (S19 Fig and S16 Table). In addition, seven gene sets related to the regulation of the TME, namely, ‘COLLAGEN CATABOLIC PROCESS’, ‘COLLAGEN METABOLIC PROCESS’, ‘MUCOPOLYSACCHARIDE METABOLIC PROCESS’, ‘POSITIVE REGULATION OF CELLULAR EXTRAVASATION’, ‘EXTRACELLULAR MATRIX DISASSEMBLY’, ‘NEGATIVE REGULATION OF HOMOTYPIC CELL-CELL ADHESION’ and ‘TISSUE REMODELING’ (hereinafter referred to as ‘TME signatures’), were significantly enriched in nonresponders (Fig 6A and S16 Table). It has been well documented that the TME signatures identified contribute to the structural organization and metabolic regulation of the extracellular matrix (ECM), which primarily consists of collagens, mucopolysaccharides and proteoglycans [43–45]. ECM remodeling in the TME plays important roles in various biological processes, including proliferation, adhesion, angiogenesis and metastasis, to affect tumor progression [46]. In many solid tumors, the TME exhibits excessive deposition of collagen and other ECM components. Dense collagen and other ECM components give rise to an immunosuppressive and hypoxic microenvironment. Furthermore, the hypoxic TME can not only accelerate tumor proliferation and metastasis but also promote the development of immunosuppressive conditions favoring tumor immune evasion [47]. Thus, it is conceivable that the enrichment of TME signatures in nonresponders may reflect the more immunosuppressive TME. Conversely, the enrichment of mitochondrial signatures in responders probably means a relatively oxygen-rich and less immunosuppressive TME.

[Figure omitted. See PDF.]

Fig 6. Singscore in LUSC.

A, For representative GSEA gene sets that were enriched in the pretreatment tumor tissues of responders or nonresponders, the single-sample gene set enrichment scores and dispersions were calculated using singscore and visualized on scatter plots. Red circles denote responders (n = 3); cyan circles, nonresponders (n = 10). The enrichment scores were analyzed using the Wilcoxon rank sum test. All calculated p-values are shown on the plots (*p < 0.05). B, Box plots indicating the single-sample enrichment scores of IFN and APP signatures calculated from four different groups with LUSC: pretreatment WB of nonresponders (Pre-Tx/NR, n = 13), on-treatment WB of nonresponders (On-Tx/NR, n = 13), pretreatment WB of responders (Pre-Tx/R, n = 4) and on-treatment WB of responders (On-Tx/R, n = 2). Red dots denote responders; cyan dots, nonresponders. The differences in the single-sample enrichment scores between the groups (i.e., ‘Pre-Tx/NR vs. On-Tx/NR’, ‘Pre-Tx/NR vs. On-Tx/R’, ‘Pre-Tx/R vs. On-Tx/R’ and ‘On-Tx/NR vs. On-Tx/R’) were evaluated by the Wilcoxon rank sum test. No significant differences were observed between the groups.

https://doi.org/10.1371/journal.pone.0260500.g006

We next investigated whether the enrichment scores of the IFN and APP signatures, which were candidate predictive biomarkers in LUAD patients treated with nivolumab monotherapy, also dynamically changed between responders and nonresponders or between the pre- and post-nivolumab monotherapy settings in WB of LUSC patients (Fig 6B). We observed no significant differences in the enrichment scores of the IFN and APP signatures between responders and nonresponders in pretreatment WB, indicating that both responders and nonresponders have similar levels of preexisting antitumor immunity. The increase in type I IFN signaling after nivolumab monotherapy observed in WB of nonresponders with LUAD was not observed in those with LUSC (Fig 6B). Collectively, these findings suggest that nivolumab success in LUSC depends entirely on the extent of the immunosuppressive TME, not on the inherent immunogenicity of the tumor itself.

Spearman rank correlation analysis demonstrated significant correlations between PFS time and the enrichment scores of the mitochondrial and TME signatures (Fig 7A). Notably, the TME signatures were negatively correlated with the mitochondrial signatures, supporting the hypothesis that the oxygenation and immunomodulation status of the TME can explain the enrichment status of both the mitochondrial and TME signatures. Among the TME signature scores, the ADHES score was most strongly correlated with PFS time (ρ = −0.786, p = 0.0015). Using a cubic regression spline, we found that the ADHES score is a potential predictive biomarker for PFS in LUSC patients treated with nivolumab monotherapy (R-sq = 0.6009, AICc = 28.1764, p = 0.0009) (Fig 7B). Gene sets related to the regulation of platelet activation (PLAT1 and PLAT2; hereinafter referred to as ‘platelet signatures’) were negatively correlated with PFS time (PLAT1, ρ = −0.687, p = 0.0095; PLAT2, ρ = −0.659, p = 0.0142) and positively correlated with TME signatures (Fig 7A). Platelets can suppress T cell antitumor responses through the production and activation of immunosuppressive factors [48, 49]. Hence, platelet signatures are inferred to reflect the immunosuppressive conditions in the TME, similar to TME signatures. The fitted models of the enrichment scores of platelet signatures exhibited a relatively low but statistically significant predictive power for PFS time (PLAT1, R-sq = 0.3135, AICc = 35.2646, p = 0.0289; PLAT2, R-sq = 0.4101, AICc = 34.4125, p = 0.0116) (Fig 7B and S22 Fig). When two gene sets in the TME and platelet signatures were combined, the LUSC patients with low values for both scores tended to have longer PFS times (Fig 8A and S23 Fig). When the LUSC patients were stratified into subsets with high and low scores (stratification by the median), there were quite large differences in PFS time between the patients with high enrichment of signatures and those with low enrichment of signatures (Fig 8B and S24 Fig). In multivariate Cox regression analysis revealed that the COL1/TISSUE signature remained an independent predictive factor of PFS time (Fig 8). Other combinations of two of the TME signatures seemed to be slightly better predictive factors compared to clinical factors (age, ECOG PS, Brinkman index, PD-L1 TPS and CNS metastasis), but without statistical significance (S25 Fig). These findings indicate that TME signatures may have potential predictive performance for PFS in LUSC patients treated with nivolumab monotherapy.

[Figure omitted. See PDF.]

Fig 7. TME signatures predict the nivolumab response in LUSC.

A, Spearman correlation matrix between PFS and single-sample gene set enrichment scores of MIT and TME signatures in pretreatment tumor tissues. The upper triangular region shows the values of Spearman’s ρ correlation coefficients (significant correlations are in bold; *p < 0.05, **p < 0.01 and ***p < 0.001). In the lower triangular region, positive correlations are visualized in red and negative correlations in blue. The color intensity is proportional to the value of Spearman’s ρ and the size of the circle to the p-value. B, Scatter plots showing the relationships between PFS and the enrichment score for single gene sets in the TME signature in pretreatment tumor tissues, with a fitted line representing the regression model using a cubic spline and 95% confidence interval. The accuracy of the fit was assessed by calculating the adjusted R-squared (R-sq) and p-values (*p < 0.05, **p < 0.01 and ***p < 0.001).

https://doi.org/10.1371/journal.pone.0260500.g007

[Figure omitted. See PDF.]

Fig 8. Two gene sets in the TME signatures as candidate biomarkers for the nivolumab response in LUSC.

A, Bubble plots showing the relationships between PFS and the enrichment scores for two gene sets in the TME signature in pretreatment tumor tissues. Each bubble represents a patient, and the size of the bubble is proportional to the PFS time. On a gradient color scale based on the PFS time, bubbles representing responders were assigned colors ranging from white to dark red; nonresponders, ranging from white to lavender. B, Kaplan-Meier PFS curves for patients stratified by the enrichment score (‘Low’ vs. ‘High’) of two gene sets in the TME signature. Patients with both scores below the median are defined as ‘Low’; the others, as ‘High’. The p-values were calculated by the two-sided log-rank test (*p < 0.05). C, Forest plot showing multivariate Cox regression analysis for potential factors associated with prolonged or shortened PFS time. Squares represent estimated hazard ratios and whiskers represent the 95% confidence intervals. Hazard ratios less than 1 indicate improved PFS time (*p < 0.05).

https://doi.org/10.1371/journal.pone.0260500.g008

Discussion

PD-L1 expression on tumor cells has been widely accepted as a predictive biomarker for therapeutic decision making in NSCLC, although its accuracy is limited and it has virtually no predictive value in patients with LUSC. To the best of our knowledge, the cause of the discrepancy in the reliability of its predictive power between LUAD and LUSC remains to be clarified. In this study, we employed whole-transcriptome sequencing and a single-sample enrichment analysis method—singscore—and found that the discrepancy can be explained by the difference in the immunogenicity of the tumor itself and the immunosuppressive conditions in the TME.

In LUAD, we observed that the IFN and APP signatures, which are closely related to each other and functionally cooperate to activate the antitumor immune response, were significantly enriched in pretreatment WB of responders (Fig 1 and S6 Table). This enrichment of the IFN and APP signatures suggests that responders may have preexisting antitumor immunity prior to nivolumab monotherapy. In contrast, transcriptomic data of pretreatment tumor tissues showed no significant enrichment of either IFN or APP signatures, suggesting that antitumor immunity is locally disturbed in the TME to support tumor progression (Fig 1 and S5 Table). These findings highlighted the common features in responders as follows: First, the tumors possess sufficiently high immunogenicity to induce effective antigen presentation to host immune components (S26 Fig). Second, the host immune system activates IFN signaling to activate tumor immunosurveillance, which is disrupted in the TME so that tumor immune evasion (i.e., adaptive immune resistance) is established locally. Given that activated type II IFN signaling can upregulate PD-L1 expression on tumor cells [35], it is highly likely that the PD-1/PD-L1 axis is responsible for the establishment of adaptive immune resistance. Finally, nivolumab monotherapy restores antitumor activity by inhibiting the PD-1/PD-L1 axis. In this regard, it is possible to say that nivolumab functions just as a molecular targeted agent in LUAD patients and that PD-L1 expression on tumor cells helps to predict the efficacy of nivolumab monotherapy.

In LUSC, TME signatures were significantly enriched in nonresponders, and this enrichment indicates the immunological condition within the TME. We found that the enrichment scores of TME signatures were negatively correlated with PFS time (Figs 7 and 8), indicating that patients with tumors strongly protected by the immunosuppressive TME are unlikely to benefit from nivolumab monotherapy. More importantly, IFN and APP signatures were strongly enriched in pretreatment WB of responders with LUAD, but this relationship was not observed in those with LUSC (Fig 6B). This observation raises the possibility that there may be no difference in the level of preexisting anti-umor immunity or, alternatively, in the immunogenicity of tumors between responders and nonresponders with LUSC. In support of this hypothesis, no significant enrichment of IFN and APP signatures was observed in on-treatment WB of responders (S15 Table). In fact, we identified the similarity in the patterns of enriched gene sets between pretreatment and on-treatment WB of patients with LUSC. Thus, nivolumab monotherapy has just a local impact in LUSC patients, a striking contrast to its systemic impact in LUAD patients treated with nivolumab monotherapy. Specifically, the clinical response to nivolumab monotherapy is actually determined by the extent of the immunosuppressive TME, where immunosuppressive factors other than the PD-1/PD-L1 axis may be considered to be critical, because the preexisting IFN and APP signatures exhibited no differences between responders and nonresponders, as mentioned above. This hypothesis is in line with the well-known finding that PD-L1 expression does not correlate significantly with clinical outcomes in LUSC patients treated with anti-PD-1 monotherapy [1]. Collectively, these findings prompt us to note that nivolumab monotherapy functions just as an immunomodulating agent and cannot overcome the highly immunosuppressive TME alone (S26 Fig).

One limitation of our study is the lack of validation in an independent cohort. However, by using a true single-sample enrichment approach, singscore, we have devised a new workflow for identifying gene expression signatures to predict a patient’s response to immunotherapy and to gain a deeper understanding of cancer biology. For example, combination strategies to enhance the immunogenicity of the tumor itself (e.g., cancer vaccines and CAR-T therapy) [4, 41–52] can be expected to improve the clinical response to nivolumab monotherapy in patients with LUAD, whereas combinational strategies to overcome the immunosuppressive TME are needed in LUSC [53, 54]. We envision that future studies will provide a blueprint for innovating combination immunotherapy approaches and optimizing patient selection and treatment strategies to improve long-term clinical outcomes in NSCLC.

Supporting information

S1 Fig. A schematic illustration of patient enrollment and sample collection.

https://doi.org/10.1371/journal.pone.0260500.s001

(TIFF)

S2 Fig. Visualization of tumor response.

A, Waterfall plot of the best percentage change from baseline during nivolumab monotherapy according to RECIST v1.1. B, Swimmer plot of all 40 patients treated with nivolumab monotherapy. PD, progressive disease.

https://doi.org/10.1371/journal.pone.0260500.s002

(TIFF)

S3 Fig. Kaplan-Meier estimates of progression-free survival (PFS) and overall survival (OS) of patients according to different stratification schemes.

A–B, the total patient cohort and C–D, LUAD or E–F, LUSC patients with PD-L1 TPS ≥ 1% versus < 1%. The p-values were calculated by the two-sided log-rank test.

https://doi.org/10.1371/journal.pone.0260500.s003

(TIFF)

S4 Fig. The association between the immune cytolytic activity in tumor tissues and clinical outcomes of patients in this cohort.

For each tumor sample, the immune cytolytic activity was estimated as the average expression level of the marker genes (CD8A, CD8B, GZMA, GZMB and PRF1). Patients with all expression levels above the average are defined as ‘High’; the others, as ‘Low’. In the upper panel, above-average values are shown in red. The lower panel illustrates Kaplan-Meier estimates of PFS and OS of patients stratified by the estimated immune cytolytic activity. The p-values were calculated by the two-sided log-rank test.

https://doi.org/10.1371/journal.pone.0260500.s004

(TIFF)

S5 Fig. Differential gene expression analysis in LUAD.

A–B, MA plot (A) and volcano plot (B) of DEGs in pretreatment tumor tissues. C–D, MA plot (C) and volcano plot (D) of DEGs in pretreatment WB. E–F, MA plot (E) and volcano plot (F) of DEGs in on-treatment WB. Red dots represent DEGs [adjusted p-value < 0.10 and |log2(fold change)| ≥ 1]. Triangles and diamonds represent genes with log2(fold change) and normalized counts, respectively, out of the plot scale. The horizontal lines in the MA plots and vertical lines in the volcano plots indicate the thresholds log2(fold change) = 1 or −1. The horizontal lines in the volcano plots indicate the threshold −log10(adjusted p-value) = 1.

https://doi.org/10.1371/journal.pone.0260500.s005

(TIFF)

S6 Fig. Heatmap representation of DEGs in LUAD.

Heatmaps of DEGs between responders and nonresponders with hierarchical clustering of samples: A, from pretreatment tumor tissues (n = 15), B, pretreatment WB (n = 20), and C, on-treatment WB (n = 15). The DEGs clearly differentiated between responders and nonresponders in all three datasets.

https://doi.org/10.1371/journal.pone.0260500.s006

(TIFF)

S7 Fig. Enrichment maps for representative gene sets significantly enriched in pretreatment WB of responders with LUAD.

Each node denotes a distinct gene set, and the size of the node is proportional to the number of genes in the set. The thickness of the edges (pale blue lines) represents the degree of overlap between the two connected gene sets.

https://doi.org/10.1371/journal.pone.0260500.s007

(TIFF)

S8 Fig. Scatter plots of single-sample gene set enrichment scores in pretreatment tumor tissues.

For each gene set in the top 20 GSEA gene sets that were enriched in pretreatment tumor tissues of responders (A) and nonresponders (B) with LUAD, the enrichment scores and dispersions were calculated using singscore. Red circles denote responders (n = 4); cyan circles, nonresponders (n = 11). The enrichment scores were analyzed using the Wilcoxon rank sum test.

https://doi.org/10.1371/journal.pone.0260500.s008

(TIFF)

S9 Fig. Scatter plots of single-sample gene set enrichment scores in pretreatment WB.

For each gene set in the top 20 GSEA gene sets that were enriched in pretreatment WB of responders (A, except for those shown in Fig 2A) and nonresponders (B) with LUAD, the enrichment scores and dispersions were calculated using singscore. Red circles denote responders (n = 5); cyan circles, nonresponders (n = 15). The enrichment scores were analyzed using the Wilcoxon rank sum test (*p < 0.05 and **p < 0.01).

https://doi.org/10.1371/journal.pone.0260500.s009

(TIFF)

S10 Fig. Scatter plots of single-sample gene set enrichment scores in on-treatment WB.

For each gene set in the top 20 GSEA gene sets that were enriched in on-treatment WB of responders (A) and nonresponders (B) with LUAD, the enrichment scores and dispersions were calculated using singscore. Red circles denote responders (n = 4); cyan circles, nonresponders (n = 11). The enrichment scores were analyzed using the Wilcoxon rank sum test (*p < 0.05 and **p < 0.01).

https://doi.org/10.1371/journal.pone.0260500.s010

(TIFF)

S11 Fig. Relationships between PFS and single-sample enrichment scores for gene sets significantly enriched in responders with LUAD.

A, Spearman correlation matrix between PFS and the enrichment scores in pretreatment WB. The upper triangular region shows the values of the Spearman’s ρ correlation coefficients (significant correlations are in bold; *p < 0.05, **p < 0.01 and ***p < 0.001). In the lower triangular region, positive correlations are visualized in red and negative correlations in blue. The color intensity is proportional to the value of Spearman’s ρ and the size of the circle to the p-value. B, Scatter plots showing the relationships between PFS and the enrichment score of the above gene sets in pretreatment WB, with a fitted line representing the regression model using a cubic spline and 95% confidence interval. The accuracy of the fit was assessed by calculating the adjusted R-squared (R-sq) and p-values (*p < 0.05 and **p < 0.01).

https://doi.org/10.1371/journal.pone.0260500.s011

(TIFF)

S12 Fig. Bubble plots showing the relationships between PFS and two single-sample enrichment scores in pretreatment WB of LUAD patients.

Each bubble represents a patient, and the size of the bubble is proportional to the PFS time. On a gradient color scale based on the PFS time, bubbles representing responders were assigned colors ranging from white to dark red; nonresponders, ranging from white to lavender.

https://doi.org/10.1371/journal.pone.0260500.s012

(TIFF)

S13 Fig. Kaplan-Meier PFS curves for patients stratified by the single-sample enrichment score (‘Low’ vs. ‘High’) for two gene sets in the IFN and APP signatures.

LUAD patients with both scores above the median are defined as ‘High’; the others, as ‘Low’. The p-values were calculated by the two-sided log-rank test (*p < 0.05, **p < 0.01, and ***p < 0.001).

https://doi.org/10.1371/journal.pone.0260500.s013

(TIFF)

S14 Fig. Forest plot showing multivariate Cox regression analysis for potential factors associated with prolonged or shortened PFS time in LUAD patients.

Squares represent estimated hazard ratios and whiskers represent the 95% confidence intervals. Hazard ratios less than 1 indicate improved PFS time (*p < 0.05).

https://doi.org/10.1371/journal.pone.0260500.s014

(TIFF)

S15 Fig. Differential gene expression analysis in LUSC.

A–B, MA plot (A) and volcano plot (B) of DEGs in pretreatment tumor tissues. C–D, MA plot (C) and volcano plot (D) of DEGs in pretreatment WB. E–F, MA plot (E) and volcano plot (F) of DEGs in on-treatment WB. Red dots represent DEGs [adjusted p-value < 0.10 and |log2(fold change)| ≥ 1]. Triangles and diamonds represent genes with log2(fold change) and normalized counts, respectively, out of the plot scale. The horizontal lines in the MA plots and vertical lines in the volcano plots indicate the thresholds log2(fold change) = 1 or −1. The horizontal lines in the volcano plots indicate the threshold −log10(adjusted p-value) = 1.

https://doi.org/10.1371/journal.pone.0260500.s015

(TIFF)

S16 Fig. Heatmap representation of DEGs in LUSC.

Heatmaps of DEGs between responders and nonresponders with hierarchical clustering of samples from pretreatment tumor tissues (n = 13).

https://doi.org/10.1371/journal.pone.0260500.s016

(TIFF)

S17 Fig. The lists of the top 25 gene sets significantly enriched in WB of responders and nonresponders with LUSC.

Top 25 enriched gene sets in responders (A) and in nonresponders (B). The left panel shows the list for pretreatment WB; the right panel, on-treatment WB. The enriched gene sets common between pretreatment and on-treatment WB are shown in bold and connected by solid lines.

https://doi.org/10.1371/journal.pone.0260500.s017

(TIFF)

S18 Fig.

Enrichment maps for representative gene sets significantly enriched in pretreatment tumor tissues of (A) responders and (B) nonresponders with LUSC. Each node denotes a distinct gene set, and the size of the node is proportional to the number of genes in the set. The thickness of the edges (pale blue lines) represents the degree of overlap between the two connected gene sets.

https://doi.org/10.1371/journal.pone.0260500.s018

(TIFF)

S19 Fig. Scatter plots of single-sample gene set enrichment scores in pretreatment tumor tissues.

For each gene set in the top 20 GSEA gene sets that were enriched in pretreatment tumor tissues of responders (A) and nonresponders (B) (except for those shown in Fig 5A) with LUSC, the enrichment scores and dispersions were calculated using singscore. Red circles denote responders (n = 3); cyan circles, nonresponders (n = 10). The enrichment scores were analyzed using the Wilcoxon rank sum test (*p < 0.05).

https://doi.org/10.1371/journal.pone.0260500.s019

(TIFF)

S20 Fig. Scatter plots of single-sample gene set enrichment scores in pretreatment WB.

For each gene set in the top 20 GSEA gene sets that were enriched in pretreatment WB of responders (A) and nonresponders (B) with LUSC, the enrichment scores and dispersions were calculated using singscore. Red circles denote responders (n = 4); cyan circles, nonresponders (n = 13). The enrichment scores were analyzed using the Wilcoxon rank sum test (*p < 0.05).

https://doi.org/10.1371/journal.pone.0260500.s020

(TIFF)

S21 Fig. Scatter plots of single-sample gene set enrichment scores in on-treatment WB.

For each gene set in the top 20 GSEA gene sets that were enriched in on-treatment WB of responders (A) and nonresponders (B) with LUSC, the enrichment scores and dispersions were calculated using singscore. Red circles denote responders (n = 2); cyan circles, nonresponders (n = 13). The enrichment scores were analyzed using the Wilcoxon rank sum test.

https://doi.org/10.1371/journal.pone.0260500.s021

(TIFF)

S22 Fig. Scatter plots showing the relationships between PFS and the single-sample enrichment score for gene sets significantly enriched in pretreatment tumor tissues of responders with LUSC.

The fitted line on each scatter plot represents the regression model using a cubic spline and 95% confidence interval. The accuracy of the fit was assessed by calculating the adjusted R-squared (R-sq) and p-values (*p < 0.05 and **p < 0.01).

https://doi.org/10.1371/journal.pone.0260500.s022

(TIFF)

S23 Fig. Bubble plots showing the relationships between PFS and two single-sample enrichment scores in the pretreatment tumor tissues of LUSC patients.

Each bubble represents a patient, and the size of the bubble is proportional to the PFS time. On a gradient color scale based on the PFS time, bubbles representing responders were assigned colors ranging from white to dark red; nonresponders, ranging from white to lavender.

https://doi.org/10.1371/journal.pone.0260500.s023

(TIFF)

S24 Fig. Kaplan-Meier PFS curves for patients stratified by the single-sample enrichment score (‘Low’ vs. ‘High’) for two gene sets in the selected TME signatures.

LUSC patients with both scores below the median are defined as ‘Low’; the others, as ‘High’. The p-values were calculated by the two-sided log-rank test (*p < 0.05 and **p < 0.01).

https://doi.org/10.1371/journal.pone.0260500.s024

(TIFF)

S25 Fig. Forest plot showing multivariate Cox regression analysis for potential factors associated with prolonged or shortened PFS time in LUSC patients.

Squares represent estimated hazard ratios and whiskers represent the 95% confidence intervals. Hazard ratios less than 1 indicate improved PFS time.

https://doi.org/10.1371/journal.pone.0260500.s025

(TIFF)

S26 Fig. Schematic depicting the different mechanisms of action of nivolumab monotherapy between LUAD and LUSC.

The success of nivolumab monotherapy depends on the inherent immunogenicity of the tumor itself in LUAD and the preexisting TME favoring an antitumor immune response in LUSC.

https://doi.org/10.1371/journal.pone.0260500.s026

(TIFF)

S1 Table. Assessment of tumor response to nivolumab monotherapy according to RECIST v1.1.

https://doi.org/10.1371/journal.pone.0260500.s027

(DOCX)

S2 Table. The lists of DEGs identified from pretreatment tumor tissue of LUAD patients.

https://doi.org/10.1371/journal.pone.0260500.s028

(XLSX)

S3 Table. The lists of DEGs identified from pretreatment whole blood of LUAD patients.

https://doi.org/10.1371/journal.pone.0260500.s029

(XLSX)

S4 Table. The lists of DEGs identified from on-treatment whole blood of LUAD patients.

https://doi.org/10.1371/journal.pone.0260500.s030

(XLSX)

S5 Table. The results from classical GSEA (FGSEA) of pretreatment tumor tissue of LUAD patients.

https://doi.org/10.1371/journal.pone.0260500.s031

(XLSX)

S6 Table. The results from classical GSEA (FGSEA) of pretreatment whole blood of LUAD patients.

https://doi.org/10.1371/journal.pone.0260500.s032

(XLSX)

S7 Table. The results from classical GSEA (FGSEA) of on-treatment whole blood of LUAD patients.

https://doi.org/10.1371/journal.pone.0260500.s033

(XLSX)

S8 Table. The results from single-sample enrichment analysis (singscore) of pretreatment tumor tissue of LUAD patients.

https://doi.org/10.1371/journal.pone.0260500.s034

(XLSX)

S9 Table. The results from single-sample enrichment analysis (singscore) of pretreatment whole blood of LUAD patients.

https://doi.org/10.1371/journal.pone.0260500.s035

(XLSX)

S10 Table. The results from single-sample enrichment analysis (singscore) of on-treatment whole blood of LUAD patients.

https://doi.org/10.1371/journal.pone.0260500.s036

(XLSX)

S11 Table. The lists of DEGs identified from pretreatment tumor tissue of LUSC patients.

https://doi.org/10.1371/journal.pone.0260500.s037

(XLSX)

S12 Table. The lists of DEGs identified from on-treatment whole blood of LUSC patients.

https://doi.org/10.1371/journal.pone.0260500.s038

(XLSX)

S13 Table. The results from classical GSEA (FGSEA) of pretreatment tumor tissue of LUSC patients.

https://doi.org/10.1371/journal.pone.0260500.s039

(XLSX)

S14 Table. The results from classical GSEA (FGSEA) of pretreatment whole blood of LUSC patients.

https://doi.org/10.1371/journal.pone.0260500.s040

(XLSX)

S15 Table. The results from classical GSEA (FGSEA) of on-treatment whole blood of LUSC patients.

https://doi.org/10.1371/journal.pone.0260500.s041

(XLSX)

S16 Table. The results from single-sample enrichment analysis (singscore) of pretreatment tumor tissue of LUSC patients.

https://doi.org/10.1371/journal.pone.0260500.s042

(XLSX)

S17 Table. The results from single-sample enrichment analysis (singscore) of pretreatment whole blood of LUSC patients.

https://doi.org/10.1371/journal.pone.0260500.s043

(XLSX)

S18 Table. The results from single-sample enrichment analysis (singscore) of on-treatment whole blood of LUSC patients.

https://doi.org/10.1371/journal.pone.0260500.s044

(XLSX)

S19 Table. The raw count data from RNS-seq.

https://doi.org/10.1371/journal.pone.0260500.s045

(CSV)

Citation: Aiba T, Hattori C, Sugisaka J, Shimizu H, Ono H, Domeki Y, et al. (2021) Gene expression signatures as candidate biomarkers of response to PD-1 blockade in non-small cell lung cancers. PLoS ONE 16(11): e0260500. https://doi.org/10.1371/journal.pone.0260500

References

1. Brahmer J, Reckamp KL, Baas P, Crinò L, Eberhardt WE, Poddubskaya E, et al. Nivolumab versus docetaxel in advanced squamous-cell non-small-cell lung cancer. N Engl J Med 2015;373: 123–135. pmid:26028407

2. Borghaei H, Paz-Ares L, Horn L, Spigel DR, Steins M, Ready NE, et al. Nivolumab versus docetaxel in advanced nonsquamous non-small-cell lung cancer. N Engl J Med 2015; 373: 1627–1639. pmid:26412456

3. Camidge DR, Doebele RC, Kerr KM. Comparing and contrasting predictive biomarkers for immunotherapy and targeted therapy of NSCLC. Nat Rev Clin Oncol 2019;16: 341–355. pmid:30718843

4. Wang Z, Duan J, Cai S, Han M, Dong H, Zhao J, et al. Assessment of blood tumor mutational burden as a potential biomarker for immunotherapy in patients with non-small cell lung cancer with use of a next-generation sequencing cancer gene panel. JAMA Oncol 2019;5: 696–702. pmid:30816954

5. Chowell D, Krishna C, Pierini F, Makarov V, Rizvi NA, Kuo F, et al. Evolutionary divergence of HLA class I genotype impacts efficacy of cancer immunotherapy. Nat Med 2019;25: 1715–1720. pmid:31700181

6. Hwang S, Kwon AY, Jeong JY, Kim S, Kang H, Park J, et al. Immune gene signatures for predicting durable clinical benefit of anti-PD-1 immunotherapy in patients with non-small cell lung cancer. Sci Rep 2020;10: 643. pmid:31959763

7. Routy B, Le Chatelier E, Derosa L, Duong CPM, Alou MT, Daillère R, et al. Gut microbiome influences efficacy of PD-1-based immunotherapy against epithelial tumors. Science 2018;359: 91–97. pmid:29097494

8. Garon EB, Rizvi NA, Hui R, Leighl N, Balmanoukian AS, Eder JP, et al. Pembrolizumab for the treatment of non-small-cell lung cancer. N Engl J Med 2015;372: 2018–2028. pmid:25891174

9. Abdel-Rahman O. Correlation between PD-L1 expression and outcome of NSCLC patients treated with anti-PD-1/PD-L1 agents: A meta-analysis. Crit Rev Oncol Hematol 2016;101: 75–85. pmid:26969107

10. Hackl H, Charoentong P, Finotello F, Trajanoski Z. Computational genomics tools for dissecting tumour-immune cell interactions. Nat Rev Genet 2016;17: 441–458. pmid:27376489

11. Hugo W, Zaretsky JM, Sun L, Song C, Moreno BH, Hu-Lieskovan S, et al. Genomic and transcriptomic features of response to anti-PD-1 therapy in metastatic melanoma. Cell 2016;165: 35–44. pmid:26997480

12. Auslander N, Zhang G, Lee JS, Frederick DT, Miao B, Moll T, et al. Robust prediction of response to immune checkpoint blockade therapy in metastatic melanoma. Nat Med 2018;24: 1545–1549. pmid:30127394

13. Furuya N, Matsumoto S, Kakinuma K, Morikawa K, Inoue T, Saji H, et al. Suitability of transbronchial brushing cytology specimens for next-generation sequencing in peripheral lung cancer. Cancer Sci 2021;112: 380–387. pmid:33124129

14. Kim D, Paggi JM, Park C, Bennett C, Salzberg SL. Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nat Biotechnol 2019;37: 907–915. pmid:31375807

15. Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 2014;30: 923–930. pmid:24227677

16. Wagner GP, Kin K, Lynch VJ. Measurement of mRNA abundance using RNA-seq data: RPKM measure is inconsistent among samples. Theory Biosci 2012;131: 281–285. pmid:22872506

17. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 2014;15: 550. pmid:25516281

18. Gu Z, Eils R, Schlesner M. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics 2016;32: 2847–2849. pmid:27207943

19. Mootha VK, Lindgren CM, Eriksson KF, Subramanian A, Sihag S, Lehar J, et al. PGC-1alpha-responsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes. Nat Genet 2003;342: 67–73. pmid:12808457

20. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci USA 2005;102: 15545–15550. pmid:16199517

21. Korotkevich G, Sukhov V, Sergushichev A. Fast gene set enrichment analysis. BioRxiv [Preprint]. 2021 bioRxiv 060012 [posted 2021 Feb 1]: Available from: http://biorxiv.org/content/early/2016/06/20/060012

22. Foroutan M, Bhuva DD, Lyu R, Horan K, Cursons J, Davis MJ. Single sample scoring of molecular phenotypes. BMC Bioinformatics 2018;19: 404. pmid:30400809

23. Bhuva DD, Foroutan M, Xie Y, Lyu R, Cursons J, Davis MJ. Using singscore to predict mutation status in acute myeloid leukemia from transcriptomic signatures. F1000Res 2019;8: 776. pmid:31723419

24. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 2010;26: 139–140. pmid:19910308

25. Wood SN. Generalized Additive Models: An Introduction with R, 2nd ed. London: Chapman & Hall/CRC; 2017.

26. Burnham KP, Anderson DR. Model selection and multimodel inference: a practical information-theoretic approach, 2nd ed. New York: Springer-Verlag; 2002.

27. Rooney MS, Shukla SA, Wu CJ, Getz G, Hacohen N. Molecular and genetic properties of tumors associated with local immune cytolytic activity. Cell 2015;160: 48–61. pmid:25594174

28. Helmink BA, Reddy SM, Gao J, Zhang S, Basar R, Thakur R, et al. B cells and tertiary lymphoid structures promote immunotherapy response. Nature 2020;577: 549–555. pmid:31942075

29. Petitprez F, de Reyniès A, Keung EZ, Chen TW, Sun CM, Calderaro J, et al. B cells are associated with survival and immunotherapy response in sarcoma. Nature 2020;577: 556–560. pmid:31942077

30. Cabrita R, Lauss M, Sanna A, Donia M, Skaarup Larsen M, Mitra S, et al. Tertiary lymphoid structures improve immunotherapy and survival in melanoma. Nature 2020;577: 561–565. pmid:31942071

31. Germain C, Gnjatic S, Tamzalit F, Knockaert S, Remark R, Goc J, et al. Presence of B cells in tertiary lymphoid structures is associated with a protective immunity in patients with lung cancer. Am J Respir Crit Care Med 2014;189: 832–844. pmid:24484236

32. Tang J, Ramis-Cabrer D, Curull V, Wang X, Mateu-Jiménez M, Pijuan L, et al. B cells and tertiary lymphoid structures influence survival in lung cancer patients with resectable tumors. Cancers (Basel) 2020;12: 2644. pmid:32947928

33. Zitvogel L, Galluzzi L, Kepp O, Smyth MJ, Kroemer G. Type I interferons in anticancer immunity. Nat Rev Immunol 2015;15: 405–414. pmid:26027717

34. Budhwani M, Mazzieri R, Dolcetti R. Plasticity of type I interferon-mediated responses in cancer therapy: from anti-tumor immunity to resistance. Front Oncol 2018;8: 322. pmid:30186768

35. Ivashkiv LB. IFNγ: signalling, epigenetics and roles in immunity, metabolism, disease and cancer immunotherapy. Nat Rev Immunol 2018;18: 545–558. pmid:29921905

36. Weinberg SE, Sena LA, Chandel NS. Mitochondria in the regulation of innate and adaptive immunity. Immunity 2015;42: 406–417. pmid:25786173

37. Franchina DG, Dostert C, Brenner D. Reactive oxygen species involvement in T cell signaling and metabolism. Trends Immunol 2018;39: 489–502. pmid:29452982

38. Barbie DA, Tamayo P, Boehm JS, Kim SY, Moody SE, Dunn IF, et al. Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1. Nature 2009;462: 108–112. pmid:19847166

39. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 2013;14: 7. pmid:23323831

40. Jacquelot N, Yamazaki T, Roberti MP, Duong CPM, Andrews MC, Verlingue L, et al. Sustained Type I interferon signaling as a mechanism of resistance to PD-1 blockade. Cell Res 2019;29: 846–861. pmid:31481761

41. Sha D, Jin Z, Budczies J, Kluck K, Stenzinger A, Sinicrope FA. Tumor mutational burden as a predictive biomarker in solid tumors. Cancer Discov 2020;10: 1808–1825. pmid:33139244

42. Gandara DR, Paul SM, Kowanetz M, Schleifman E, Zou W, Li Y, et al. Blood-based tumor mutational burden as a predictor of clinical benefit in non-small-cell lung cancer patients treated with atezolizumab. Nat Med 2018;24: 1441–1448. pmid:30082870

43. Xu S, Xu H, Wang W, Li S, Li H, Li T, et al. The role of collagen in cancer: from bench to bedside. J Transl Med 2019;17: 309. pmid:31521169

44. Henke E, Nandigama R, Ergün S. Extracellular matrix in the tumor microenvironment and its impact on cancer therapy. Front Mol Biosci 2020;6: 160. pmid:32118030

45. Janiszewska M, Primi MC, Izard T. Cell adhesion in cancer: beyond the migration of single cells. J Biol Chem 2020;295: 2495–2505. pmid:31937589

46. Winkler J, Abisoye-Ogunniyan A, Metcalf KJ, Werb Z. Concepts of extracellular matrix remodelling in tumour progression and metastasis. Nat Commun 2020;11: 5120. pmid:33037194

47. Rankin EB, Giaccia AJ. Hypoxic control of metastasis. Science 2016;352: 175–180. pmid:27124451

48. Rachidi S, Metelli A, Riesenberg B, Wu BX, Nelson MH, Wallace C, et al. Platelets subvert T cell immunity against cancer via GARP-TGFβ axis. Sci Immunol 2017;2: eaai7911. pmid:28763790

49. Gaertner F, Massberg S. Patrolling the vascular borders: platelets in immunity to infection and cancer. Nat Rev Immunol 2019;19: 747–760. pmid:31409920

50. Takahashi H, Shimodaira S, Ogasawara M, Ota S, Kobayashi M, Abe H, et al. Lung adenocarcinoma may be a more susceptive subtype to a dendritic cell-based cancer vaccine than other subtypes of non-small cell lung cancers: a multicenter retrospective analysis. Cancer Immunol Immunother 2016;65: 1099–1111. pmid:27448677

51. Mao Y, Fan W, Hu H, Zhang L, Michel J, Wu Y, et al. MAGE-A1 in lung adenocarcinoma as a promising target of chimeric antigen receptor T cells. J Hematol Oncol 2019;12: 106. pmid:31640756

52. Srivastava S, Furlan SN, Jaeger-Ruckstuhl CA, Sarvothama M, Berger C, Smythe KS, et al. Immunogenic chemotherapy enhances recruitment of CAR-T cells to lung tumors and improves antitumor efficacy when combined with checkpoint blockade. Cancer Cell 2021;39: 193–208. pmid:33357452

53. Altorki NK, Markowitz GJ, Gao D, Port JL, Saxena A, Stiles B, et al. The lung microenvironment: an important regulator of tumour growth and metastasis. Nat Rev Cancer 2019;19: 9–31. pmid:30532012

54. Tang T, Huang X, Zhang G, Hong Z, Bai X, Liang T. Advantages of targeting the tumor immune microenvironment over blocking immune checkpoint in cancer immunotherapy. Signal Transduct Target Ther 2021;6: 72. pmid:33608497