Development of an immune-related gene signature applying Ridge method for improving immunotherapy responses and clinical outcomes in lung adenocarcinoma

The acquisition of LUAD patient datasets

The transcriptional profiles of patients with LUAD and their clinical information were collected from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov) and Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo) under the following criteria: (1) Complete data on overall survival (OS) and survival status, stage, gender, age; (2) IlluminaHiSeq platform or Affymetrix HG-U133_Plus 2.0 platform; and (3) Sample size larger than 200. After filtering, 1,327 LUAD patients from five datasets were included in this study. Specifically, GSE30219 dataset contained 83 patients (Rousseaux et al., 2013); GSE31210 dataset contained 226 patients (Okayama et al., 2012); GSE50081 dataset contained 127 patients from (Der et al., 2014); GSE72094 dataset contained 398 patients (Schabath et al., 2016); 493 patients were collected from TCGA-LUAD (https://xenabrowser.net/datapages/?cohort=GDC%20TCGA%20Lung%20Adenocarcinoma%20(LUAD)&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443). The TCGA-LUAD dataset served as a training dataset for analyzing the features of immune genes, while the other four datasets were independent test datasets for validation. Detailed clinical information of the patients from these five groups was shown in Table 1.

Gene expression profiling data processing

Background correction, quantile normalization, and log-2 transformation of GEO microarray datasets were processed in the R package oligo (Carvalho & Irizarry, 2010), and the Robust Multi-array Average algorithm was used for normalizing all the microarray datasets. According to the annotation files of the microarrays, the median value of multiple probes mapping to the same gene was taken, while a probe was deleted when it was mapped to multiple genes. Subsequently, transcripts per kilobase million (TPM) was converted from the RNA-seq read count in TCGA-LUAD.

Assessment of immune infiltrating cells

Marker genes for 28 types of highly reliable immune-infiltrating cells were collected from a previous study (Charoentong et al., 2017). Single-sample gene set enrichment analysis (ssGSEA) in the R package GSVA (Hänzelmann, Castelo & Guinney, 2013) was conducted for calculating the proportion of immune-infiltrating cells in each patient, followed by verifying the results by ESTIMATE (Yoshihara et al., 2013), CIBERPSORT (Newman et al., 2015), TIMER (Li et al., 2016), EPIC (https://gfellerlab.shinyapps.io/EPIC_1-1/) (Racle et al., 2017), and MCP-counter (Becht et al., 2016) algorithms.

Classification of immune-related molecular subtypes

Patients in the TCGA-LUAD cohort were classified applying a clustering method of unsupervised resampling based on the components of 28 types of immune infiltrating cells. The R package ConsensusClusterPlus (Wilkerson & Hayes, 2010) was used to calculate immune infiltrating cell scores for individual samples as a data matrix. Specifically, 80% of the samples were randomly selected for each iteration, and the distance between samples was measured by Pearson correlation coefficient. This was process was repeated 1,000 times, where samples were partitioned into up to k (maximum k = 10) clusters using the partioning around medoids (PAM) algorithm. Subsequently, the optimal number of clusters was determined according to consensus score matrix, cumulative distribution function (CDF) curve and proportion of ambiguous clustering (Pagès et al., 2018) score. The significance of the clustering results was assessed using sigclust (Harrison et al., 2003).

Screening tumor-infiltrating immune-related gene markers

A novel computational framework for a comprehensive analysis of immune-infiltrating cells in the TME was created to identify immune genes related to tumor infiltration. Multiple machine learning methods were used to establish an immunity prognostic signature (IMMPS). The median absolute deviation (MAD) (Okayama et al., 2012) of all genome-wide protein-coding genes was computed, and potentially deregulated genes were defined as having a MAD greater than top 50%. Weighted gene co-expression network analysis (WGCNA) (Langfelder & Horvath, 2008) was constructed and the appropriate soft threshold β was calculated to satisfy a scale-free network. Additionally, the weighted adjacency matrix was converted to topological overlap matrix (TOM) and modules were identified by the dynamic tree cutting method. Next, immune-related modules were recognized based on the correlation of each module with immune clusters. The relationship between the prognosis of LUAD and genes in the immune-related modules was analyzed by univariate Cox regression analysis to filter significantly prognostically relevant genes as a hub gene set.

In addition, a total of 101 combinations of 10 machine learning algorithms, including Least absolute shrinkage and selection operator (LASSO), CoxBoost, Ridge, generalized boosted regression modeling (GBM), partial least squares regression for Cox (plsRcox), supervised principal components (SuperPC), elastic network (Enet), stepwise Cox, random survival forest (RSF), survival support vector machine (survival-SVM) together with 10-fold cross-validation were introduced to select candidate genes with the highest C-index for IMMPS. The stepCox (both) and Ridge were used to develop the IMMPS model. Specifically, the former was applied to filter the most valuable TIICGs, while the later was used to fit the most reliable model. After the Ridge regression, the IMMPS formula for each variable was as follow:

$$IMMPS=\sum _{k=1}^n{x}^k{a}^k$$where n denoted the number of feature genes, $a^k$ denoted the ridge regression coefficient of the k^th feature gene, and $x^k$ denoted the expression level of the k^th feature gene.

Enrichment analysis

Gene oncology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of specific gene modules were performed using the R package “clusterProfiler” (Yu et al., 2012).

Immunotherapy response prediction

GSE126044 (Cho et al., 2020) and GSE135222 (Kim, Choi & Jung, 2020) were two cohorts containing follow-up information of each NSCLC patient’s response to anti-PD-1/PD-L1 therapy. Subclass mapping was conducted using the SubMap (Hoshida et al., 2007) algorithm to evaluate the response of immune molecular subtypes to immunotherapy, and potential treatment benefit was determined by Tumor Immune Dysfunction and Exclusion (TIDE) algorithm.

Cell culture and transfection

Human normal lung epithelial cells BEAS-2B (BNCC359274) and human LUAD cells lines PC-9 (BNCC340767) and H1395 (BNCC100270) were purchased from Bena Biotechnology Co. (Beijing, China). Penicillin/streptomycin and Dulbecco’s modified Eagle’s medium (A1896701; DMEM, Gibco, Grand Island, NY, USA) with 10% fetal bovine serum (FBS, Gibco, USA) was used for cell culturing at 37 °C in 5% CO₂.

BTNL9 plays a critical role in immune regulation, especially in TME that may affect the activity of immune infiltrating cells (Zhang et al., 2023). In this study, to investigate the regulatory effects of BTNL9 on the expression of the seven immune-related genes, human LUAD cell lines (PC9 and H1395) and normal lung epithelial cells (BEAS-2B) were transfected with BTNL9-specific siRNA (5′-3′ sequence: GCCTCTAACTCCACAACAACACT) using Lipofectamine 3000 (Thermo Fisher, Waltham, MA, USA) following the manufacturer’s protocol. After 48 h, the transfection efficiency was tested by qRT-PCR.

Western blot testing

Cells were lysed using RIPA buffer supplemented with protease inhibitors. Protein samples were separated on 15% SDS-PAGE gel, subsequently transferred onto PVDF membranes (Beyotime, Shanghai, China) and blocked. The membranes were incubated overnight at 4 °C with primary antibodies against BTNL9 (1:1,000) and GAPDH (1:1,000), and then treated with HRP-conjugated secondary antibodies for 1 h at room temperature. Finally, the protein bands were visualized utilizing the ECL system (Amersham Biosciences, Inc, Buckinghamshire, UK) and the band intensities were measured using Quantity One software (BioRad, Hercules, CA, USA).

Quantitative real-time PCR (qRT-PCR)

Total RNA was extracted from BEAS-2B, PC9 and H1395 (Thermo Fisher, Waltham, MA, USA) using TRIzol reagent. The HiScript II SuperMix (Vazyme, China) was employed for isolating cDNA from 500 ng of RNA. QRT-PCR was carried out with the use of the SYBR Green Master Mix in ABI 7500 System (Thermo Fisher Scientific, Waltham, MA, USA). PCR amplification began with 45 cycles a 94 °C for 10 min, at 94 °C for 10 seconds (s), and 60 °C for 45 s. GAPDH was an internal reference. All the primer sequences were shown in Table 2.

Statistical analyses

Heatmap was generated using the R package ComplexHeatmap (Gu, Eils & Schlesner, 2016). Wilcox test and Kruskal-Wallis test were employed to analyze two-group differences and comparisons between multiple groups, respectively. The correlation between two groups of samples was calculated by Pearson’s correlation test. Prognostic differences among patients were assessed using Kaplan-Meier. The R package “timeROC” was used to plot time-dependent receiver operating characteristic (ROC) curves (Blanche, Dartigues & Jacqmin-Gadda, 2013). All the statistical analyses were conducted in R 3.6.4 (R Core Team, 2021). GraphPad Prism V.8.0.2 was applied for the statistical analysis on the experimental data, which were compared using one-way analysis of variance or unpaired t-test. A p < 0.05 indicated a statistical significance.

Development and validation of immune infiltration molecular subtypes

Charoentong et al. (2017), Bindea et al. (2013) revealed that at least 28 different types of immune cell types are crucially involved in tumor infiltration. Immune cell infiltration in the TME of LUAD were analyzed by single-sample gene set enrichment analysis (ssGSEA). Patients were classified by consensus clustering. CDF curves of the consensus score matrix (Fig. 1A) and PAC results (Fig. 1B) showed that the optimal number of clusters was when k = 2 (Fig. 1C). The sigclust analysis detected significant differences (p = 1.141007e−05) between the two consensus clusters of C1 and C2 in terms of immune infiltration, with C1 showing a notably lower overall infiltration than C2 (Fig. 1D, Fig. S2A). In addition, ESTIMATE, CIBERPSORT, TIMER, EPIC, and MCP-counter results showed similar results (Fig. 1E, Figs. S2B–S2F). Hence, C1 and C2 were accordingly defined as an immunosuppressed and immune-activated subtypes. Survival analyses all demonstrated a significantly worse prognosis for C1 than C2 (Fig. 1F).

The expression patterns of immune checkpoint genes in immune-activated subtype and potential benefits from immunotherapy

After literature review, a total of 79 immune checkpoint genes (ICGs) (Table S1) mainly involved in ligand-receptor interactions were selected. These ICGs affect immune activity in different ways, including stimulation, inhibition or a combination of both. Analysis of expression patterns of these genes in the C1 and C2 subtypes demonstrated that almost all the ICGs, such as PDCD1, CD274 and CTLA4 (Fig. 2B), were upregulated in C2 (Fig. 2A). TIDE model could evaluate T cell exclusion and dysfunction signatures across 33,000 samples in 188 tumor cohorts stored in public databases, including TCGA (Weinstein et al., 2013), METABRIC (Curtis et al., 2012), and PRECOG (Gentles et al., 2015). Applying TIDE model, it was observed that the C1 subtype had significantly higher T cell exclusion and the C2 subtype had significantly higher T cell dysfunction, and that C2 subtype might benefit from taking immunotherapy based on combined TIDE scores (Fig. 2C). Recent studies have confirmed immunophenoscore (IPS) as an accurate predictor of ICI response, with a higher IPS predicting greater sensitivity to ICI treatment. Here, we found a significantly higher IPS score of C2 than C1 (Fig. 2D, p = 0.045), which was also validated by the cytotoxic T lymphocyte score (Fig. 2E). Further, the gene expression patterns of immune phenotypes of NSCLC patients who received ICIs were compared by SubMap analysis. Notably, the gene expression profiles of lung cancer patients who responded to anti-PD-1 immunotherapy were highly similar to those with immune-activated subtypes in both the training and validation cohorts (Fig. 2F). These results indicated that patients with immune-activated subtypes could benefit from receiving ICI immunotherapy, especially more from anti-PD-1 therapy.

Identification of immune infiltration-associated gene modules and prognostic features of the genes

The soft threshold β = 8 was selected (Figs. S2A, S2B) to ensure a scale-free nature of the network. A total of 30 gene modules were detected, as indicated by different colors, and the heatmap displayed the characteristic gene neighborhood of the modules (Fig. S2C). Additionally, the modules was weakly correlated with AJCC Stage, Age and Gender (Fig. 3A). Specifically, darkred, orangered4, and green modules were considered as the modules most closely correlated with C2 (immune activation) subgroup. Among the three modules, the correlation coefficient between module membership (MM) and gene significance (GS) of the darkred, orange4, and green modules reached 0.94, 0.94, and 0.71 (Figs. 3B–3D), respectively, indicating that the identification of these gene modules was reliable. Further functional enrichment analysis showed that the darkred module was enriched to pathways of human T−cell leukemia virus 1 infection, PD−L1 expression, PD−1 checkpoint pathway in cancer, etc. (Fig. 3E). The green module was significantly enriched to Cytokine-cytokine receptor interaction, B cell receptor signaling pathway, and other pathways (Fig. 3F). The orangered4 module was significantly enriched to primary immunodeficiency, MAPK signaling pathway and some other pathways (Fig. 3G). These pathways all play critical role in immune regulation. Noticeably, biological processes such as cytokine activity, cytokine receptor activity, cytokine binding, and cytokine receptor binding were enriched multiple times. Under stringent criteria, genes significantly associated with LUAD prognosis were identified by performing one-way survival analysis from at least four cohorts, and finally 26 genes were screened as potential prognostic markers with significant immune relevance in LUAD (Fig. S2E).

Development of an immune-correlated prognostic signature and assessment

Using machine learning-based integration analysis, an immune-related gene signature was developed based on the expression profiles of 26 immune-related genes. In the TCGA-LUAD dataset, 101 predictive models were developed by the global leave-one-out cross-validation (LOOCV) framework and their C-index in all the validation datasets was calculated. Interestingly, the optimal model was the combination of stepCox and Ridge, which showed the highest average C-index (0.697) in all the validation datasets (Fig. 4A). In step cox regression, when the Akaike Information Criterion (Pagès et al., 2018) reached its minimum value, the optimal combination of seven variable genes were identified. Therefore, the final model (IMMPS) was built with these seven genes applying regression analysis with Ridge: IMMPS = 0.093*SEMA7A + 0.145*EFHD2 + 0.108*CHST11 − 0.502*SLC24A4 − 0.089*MAL − 0.074*JCHAIN − 0.128* SCARF1. Each patient was assigned with a risk score by the model and weighted for their regression coefficient. Based on the optimal cut-off value, all the patients were categorized into high- and low- risk groups (Figs. 4B–4F). It was observed that the OS of high-risk patients was significantly worse than the low-risk patients in the TCGA-LUAD training dataset and the five validation datasets (all p < 0.05). ROC analysis showed that the IMMPS reached 1-, 3-, and 5-year area under the curves (AUCs) of 0.69, 0.69, and 0.65 in TCGA-LUAD, respectively. In GSE72094, the AUC values were consistent at 0.63 for all three time points. For GSE31210, 1-, 3-, and 5-year AUCs were 0.74, 0.81, and 0.74, respectively, while in GSE50081, they were 0.74, 0.78, and 0.77, respectively. Finally, the GSE30219 cohort showed AUC values of 0.88, 0.85, and 0.84 at the corresponding time points. These data supported a stable and robust performance of IMMPS across different cohorts. Consistently, multivariate Cox regression analysis revealed that IMMPS was an independent risk factor for LUAD (p < 0.001, Fig. S3).

Comparison of gene expression-based prognostic features for LUAD

Next-generation sequencing and big data technology has promoted the prediction of genes by machine learning-based methods (Ahluwalia, Kolhe & Gahlay, 2021). This study comprehensively compared the performance of published signatures to the IMMPS. A severe lack of miRNA and lncRNA data in the validation dataset excluded the miRNA and lncRNA signatures, resulting in the collection of a total of 154 signatures (Table S2). These signatures were related to different biological processes, such as ferroptosis, WNT, stemness, autophagy, hypoxia, adipogenesis, glycolysis, epigenetics, immune response, vitamin D, aging, epithelial-mesenchymal transition, N6-methyladenosine, Toll-like receptor signaling, and drug sensitivity. The c-index of each signature in all the datasets were calculated to evaluate their predictive performance (Fig. 5). According to the c-index, ranking of the signatures in each dataset showed that one signature (0.6%) outperformed the IMMPS in three datasets. Conversely, 49 signatures (31.8%) exhibited a lower c-index than the IMMPS in three cohorts (60% of the total cohorts), 47 signatures (30.5%) showed a lower c-index than the IMMPS in four cohorts (80% of the total cohorts), and 57 signatures (37%) had a lower c-index than the IMMPS across all cohorts (Fig. S4). These results supported the stability of the IMMPS. We also noted that while most of the models performed well in their own training set and some specific external datasets, their performance was notably poor in other datasets. This discrepancy could be attributed to poor model generalization due to overfitting problem. Next, two machine learning algorithms were employed to reduce the gene numbers in our signature, which thereby manifested a stronger generalization ability. In comparison to the model developed by (Deng et al., 2022), the IMMPS contained only seven genes. Fewer genes with comparable predictive performance could significantly reduce the cost in clinical testing. These results confirmed that IMMPS had the advantage of higher stability but lower detection cost in comparison with previous gene signatures.

The impact of the IMMPS on ICI treatment

Comparison on the differences of the IMMPS in the two subtypes showed a significantly lower IMMPS in C2 patients than in C1 patients (Fig. 6A, p = 0.00026). The patients were then classified into two groups of IMMPS-High and IMMPS-Low using the IMMPS. It can be observed that immune infiltrating cell scores of patients with IMMPS-High were significantly lower (Fig. 6B). A total of 18 types of immune infiltrating cells (64.3%) were negatively correlated with the IMMPS. In particular, eosinophils had the strongest correlation with IMMPS (Fig. 6C), and patients with low IMMPS but more eosinophils had a lower risk of Coronavirus disease 2019 (COVID-19) and better clinical outcomes (Xie et al., 2021). The stability of the results was verified using five other algorithms to avoid bias from different algorithms. IMMPS-High patients showed remarkably lower immune infiltrating cell scores (Figs. S5A–S5E), which was verified by correlation analysis (Fig. 6F). These results validated a negative correlation between immune infiltration and IMMPS (Fig. 6D), specifically, patients who had low IMMPS and higher immune cell infiltration could benefit more from taking immunotherapy. Additionally, analysis on the correlation between immune checkpoint genes and the IMMPS (Fig. 6E) revealed that most of the HLA Class II genes were significantly negatively related to IMMPS.

Experimental verification using cell lines

Considering that BTNL9 was the most relevant immune checkpoint gene, BTNL9 expression was inhibited in BEAS-2B, PC9, and H1395 cell lines and simultaneously the expression of seven genes in the IMMPS model was measured. First, the knockdown efficiency of BTNL9 in the LUAD cell lines was detected based on qRT-PCR and Western blot assays. As shown in Fig. 7A, knockdown of BTNL9 significantly downregulated the mRNA expression level of BTNL9 in human normal lung epithelial cells BEAS-2B and LUAD cell lines PC9 and H1395. Similarly, BTNL9 protein levels were significantly lowered in the si-BTNL9 group in comparison to the blank control group (Ctrl) and the control group transfected with non-targeting siRNA (si-NC) (Figs. 7B,7C). After BTNL9 knockdown, the mRNA expression levels of the seven genes in three cell lines (BEAS-2B, PC9, and H1395) were detected. The results showed that after knockdown of BTNL9, the expression of the seven genes in LUAD cell lines (PC9 and H1395) demonstrated significant differences than before the knockdown. This suggested that BTNL9 may be an upstream regulator of these genes, and that knockdown of BTNL9 changed the expression of these seven genes in LUAD cells, indicating that these genes may be involved in the occurrence and development of LUAD (Fig. 8).