Discovery of tumor immune infiltration-related snoRNAs for predicting tumor immune microenvironment status and prognosis in lung adenocarcinoma.

Lung adenocarcinoma (LUAD) has a high mortality rate and is difficult to diagnose and treat in its early stage. Previous studies have demonstrated that small nucleolar RNAs (snoRNAs) play a critical role in tumor immune infiltration and the development of a variety of solid tumors. However, there have been no studies on the correlation between tumor-infiltrating immune-related snoRNAs (TIISRs) and LUAD. In this study, we filtered six immune-related snoRNAs based on the tissue specificity index (TSI) and expression profile of all snoRNAs between all LUAD cell lines from the Cancer Cell Line Encyclopedia and 21 types of immune cells from the Gene Expression Omnibus database. Further, we performed real-time quantitative polymerase chain reaction (RT-qPCR) to validate the expression status of these snoRNAs on peripheral blood mononuclear cells (PBMCs) and lung cancer cell lines. Next, we developed a TIISR signature based on the expression profiles of snoRNAs from 479 LUAD patients filtered by the random survival forest algorithm. We then analyzed the value of this TIISR signature (TIISR risk score) for assessing tumor immune infiltration, immune checkpoint inhibitor (ICI) treatment response, and the prognosis of LUAD between groups with high and low TIISR risk score. Further, we found that the TIISR risk score groups showed significant differences in biological characteristics and that the risk score could be used to assess the level of tumor immune cell infiltration, thereby predicting prognosis and responsiveness to immunotherapy in LUAD patients.

Introduction

Lung cancer is one of the most common malignancies worldwide and is the leading cause of cancer-related mortality [1]. Among them, lung adenocarcinoma (LUAD) is the most common shared histological subtype, accounting for approximately 50% of the primary lung cancer. LUAD is usually diagnosed in the advanced stage and has high mortality and poor prognosis [2], [3]. Although the application of therapies such as surgery, chemotherapy, and molecular targeted therapy has improved LUAD patient survival in the past decade, the poor prognosis caused by its high invasiveness, easy metastasis, and recurrence has not been fully resolved. Immune checkpoint inhibitors (ICIs) are by far the most widely used immunotherapeutic agents for LUAD [4]. However, in clinical practice, ICIs are effective only in a minority of patients [5]. Therefore, in order to effectively guide personalized treatment, it is important to identify biomarkers to predict the treatment response and prognosis of LUAD patients.

With the development of high-throughput and accurate detecting technologies, the tumor immune microenvironment (TIME) is considered to have a significant impact on the clinical treatment response and prognosis of patients with tumors [6], [7]. Rachel and other scholars have demonstrated that immune cell infiltration in tumors is positively correlated with LUAD prognosis [8], [9], [10]. Based on this evidence, we speculate that the TIME status is a promising indicator of LUAD treatment response and prognosis.

Noncoding RNA (ncRNA) has gradually become important for gene regulation and plays a key role in the biological processes of different tumors [11]. Small nucleolar RNAs (snoRNAs) comprise the largest group of small ncRNAs in mammalian cells. They consist of 60–300 nucleotides and regulate gene expression at multiple levels, including involvement in chromatin architecture, RNA editing, and translation [12], [13], [14]. Since Chang et al. [14] first reported the relationship between snoRNA and tumors, a new understanding of the role of snoRNA in tumors has emerged. Growing evidence shows that snoRNA plays a paramount role in the progression of a variety of tumors, including gastric and colon cancers [15], [16]. A recent study suggested that snoRNA was significantly elevated in the tumor tissue and plasma of non-small cell lung cancer patients and was related to patient survival rate [14], [17]. Meanwhile, snoRNA was also associated with TIME-related features, such as CD8+ T cell infiltration, and had the potential to become a noncoding molecular marker of TIME [18]. In addition, the TIME status is actively involved in the occurrence and development of tumors. For instance, the intertumoral infiltration of cytotoxic lymphocytes is positively correlated with better response to treatment in LUAD patients [19], [20]. Accordingly, we conclude that snoRNA is associated with the immune status of tumors and may predict the treatment response and prognosis of tumor patients through the TIME, thereby accurately guiding personalized therapy.

Studies have shown that snoRNA plays a key role in tumorigenesis and tumor development; however, the role of immune cell-derived snoRNA in the TIME remains unknown. In this study, we searched for immune-related snoRNAs by comparing LUAD tumor cell lines from the Cancer Cell Line Encyclopedia (CCLE) database with the human immune cell data from the Gene Expression Omnibus (GEO) database. We then constructed a prediction model by identifying the housekeeping snoRNAs from immune cells (Fig. 1), in order to predicting the immune infiltration status, prognosis, and even the immunotherapeutic effect in LUAD patients.

Fig. 1

Strategy for identifying the tumor-infiltrating immune-related snoRNA (TIISR) signature in this study. snoRNAs were sorted in immune cells and cancer cell lines in descending order of expression. The snoRNAs that were most expressed in immune cells and least expressed in lung adenocarcinoma (LUAD) cell lines were screened, and the specificity of candidate snoRNAs expressed against different types of immune cells was calculated using the tissue specificity index (TSI). Six snoRNAs were finally included. Furthermore, the TIISR signature was identified by the RSF model; accordingly, LUAD patients were divided into high-risk and low-risk groups based on median risk score to explore the characteristics in terms of biological processes, immune infiltration status, immunotherapy response, and prognosis.

Methods

Data source

Publicly available RNA sequencing (RNA-seq) datasets were gathered from the GEO database (https://www.ncbi.nlm.nih.gov/geo), The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/), and the CCLE project (https://portals.broadinstitute.org/ccle). In detail, the transcriptional profiles of patients with LUAD (GSE81089 and TCGA-LUAD) were obtained from the GEO and TCGA databases using GEOquery [21], TCGAbiolinks [22], or manually for further analysis (Table S1). The IMvigor210 [23] dataset was downloaded from http://research-pub.gene.com/IMvigor210CoreBiologies/. 479 LUAD samples with complete clinical and expressional information from TCGA-LUAD were used for the training set. The data of 106 LUAD samples with survival information were extracted from GSE81089 dataset was used as validation set. Moreover, 348 bladder cancer patients who received ICIs treatment form IMvigor210 were used as the validation set because of lacking eligible LUAD datasets with ICIs treatment patients. Transcriptome sequencing data of purified human immune cells (GSE114765, GSE133145, GSE135635 and GSE107011) were downloaded from the GEO database using GEOquery [21] or manually (Table S2). RNA-Seq data of LUAD cell lines were obtained from the CCLE database based on histology subtype.

Processing of the transcriptomic sequencing data

All counts or fragments per kilobase of transcript per million (FPKM) mapped reads values were transformed to transcripts per kilobase million (TPM) values under GENCODE annotation (https://www.gencodegenes.org/) version 34 for further analysis. Samples with only TPM values were directly used without datatype transformation. Due to the complexity of the purified immune cell data, we excluded samples from patients with specific disease (e.g., primary Sjogren's syndrome) and then used the xCell [24] algorithm to further confirm immune cell types using the maximum enrichment score. Finally, a total of 188 samples from 21 small or 7 main categories of immune cells were enrolled in this study (Table S2). Combat from the sva package (3.35.2) was used to remove batch effects across immune cell samples.

Identifying tumor-infiltrating immune-related snoRNAs and developing a related signature

Tumor-infiltrating immune-related snoRNAs (TIISRs) were identified based on the tissue specificity index (TSI) and the expression profile between LUAD cell lines and immune cells. The TSI was developed by Yanai et al. in 2005, as follows [25]:where N is the total number of immune cell samples and is defined as the expression intensity calculated by maximal component normalization.

A higher TSI represents higher immune cell-type specificity, while a lower TSI indicates more general expression across immune cell types. We selected snoRNAs that are more commonly expressed among immune cells (TSI ≤ 0.4) but are not expressed in LUAD cell lines (housekeeping) as TIISRs. Due to unremovable batch effects between immune cells and cancer cell lines, we sorted the snoRNA expression profile by the mean value of each snoRNA. We then selected the intersection between the high expression snoRNAs among immune cells and the low expression snoRNAs among LUAD cell lines as TIISRs. Further, we excluded TIISRs that were not expressed on more than half of the samples from the training set. Then, six TIISRs (SNORD14A, SNORD59A, SNORD99, SNORD100, SNORD63 and SNORD19) were included for random survival forest (RSF) model construction. This is an ensemble-tree method that adapts the random forests to right-censored data and survival analysis [26], [27].

The log-rank score test for splitting survival trees is described by Hothorn and Lausen (2003) [28]. According to this test, we assume that the -variable has been ordered so that . The “ranks” were then calculated for each survival time where . This gives the following equation:where . Note that is the index of the order for . The log-rank score test is defined as follows:where and are the sample mean and sample variance of . Log-rank score splitting defines the measure of node separation by . Maximizing this value over and yields the best split.

The comprehensive working framework for identifying marker tumor-infiltrating snoRNAs and developing a TIISR signature is shown in Fig. 1.

Gene set variation analysis (GSVA), single-sample gene set enrichment analysis (ssGSEA), and TIME cell infiltration estimation

The “c2.cp.hallmark.v7.1.symbols” and“c5.all.v7.1.symbols” gene sets were downloaded from MSigDB and used for hallmark and Gene Ontology (GO) enrichment by GSVA and ssGSEA from “GSVA” R package [29]. The characteristics of the TIME include immune cell infiltration, anticancer immunity cycle activation, and immune checkpoint expression. TIME cell infiltration estimation was performed using the xCell [24] and CIBERSORT-Abs [30], [31] algorithms. Seven steps involved in anticancer immunity cycle activation were scored by ssGSEA based on the gene expression of each sample, and the score of each step reflected the activation of antitumor immunity [32], [33]. The inhibitory immune checkpoints were obtained from studies of Auslander et al. and Hu et al. [34], [35].

Separating peripheral blood mononuclear cells (PBMCs)

Blood samples were collected from 10 health volunteers and processed within 2 h. The whole separating was followed by product information sheet of LymphoPrep (a density gradient medium for the isolation of mononuclear cells from STEMCELL technologies, cat. 07851) and PBMCs were washed by phosphate buffer saline (PBS) and centrifuged by 400g at 4 ℃ for 5 mins. Then the cells were ready for RNA extract.

Lung cancer cell lines and cell culture

H1299 and A549 cell lines were provided by Stem Cell Bank, Chinese Academy of Sciences (Beijing, China). H1975, and PC-9 cells were kindly provided by Professor Chun Fang Zhang (Department of Thoracic Surgery, Xiangya Hospital). Cells were cultured in RPMI 1640 medium (Gibco, Waltham, MA, USA) supplemented with 10% fetal bovine serum (Gibco), 100 U/mL penicillin, and 100 μg/mL streptomycin (Gibco) and incubated at 37 °C and 5% CO2.

RNA extract and real-time quantitative polymerase chain reaction (RT-qPCR)

RNA was extracted, reverse transcripted to cDNA and RT-qPCR were all followed by the guide of miDETCT A track miRNA RT-qPCR Starter kit (RiboBio, cat. C10721-1) and performed RT-qPCR on a QuantStudio 5 real time PCR system (Applied Biosystems). The primers for SNORD63 (miRA100955), SNORD100 (miRA100954), SNORD99 (miRA100953), SNORD19 (miRA100957), SNORD59A (miRA100951), SNORD14A (miRA101316) and U6 (miRAN0002-1-100) were also provided by RiboBio. The U6 was used as loading control.

Statistical analysis

Spearman correlation was used for all correlation analyses. Kruskal–Wallis tests were used for comparing numeric values between two or more groups, except as otherwise noted. The “limma” package and the empirical Bayesian approach were used for performing differential analysis using the enrichment analysis results. The medians were considered cutoff points for the risk score and other numeric variants; cutoffs were also used for subgrouping. Survival analyses were based on the Kaplan–Meier method and log-rank tests. Univariate and multivariable Cox regression models were used for calculating hazard ratios (HRs) and identifying independent prognostic factors. Receiver operating characteristic (ROC) curves were used to assess the specificity and sensitivity of risk score based on the pROC package. P values or adjusted P-values <0.05 were considered statistically significant, while pathways with adjusted P < 0.01 were used for visualization. Data processing, analysis, and visualization were performed using R software (3.6.3).

Results

snoRNA expression in immune cells and cancer cell lines

After excluding immune cell samples that collected from patients with specific diseases, the identities of all 188 samples were further confirmed by the xCell algorithm. Seven main types of immune cells were involved: T cells, B cells, neutrophils, hematopoietic stem cells, eosinophils, monocytes/macrophages, natural killer (NK) cells, and 21 subtypes. Among them, there were eleven types of T cells, five types of monocytes/macrophages, and one type each of B cells, neutrophils, hematopoietic stem cells, eosinophils, and NK cells (Fig. 2A). To investigate the expression patterns of snoRNAs in different immune cells, we analyzed the expression levels of snoRNAs in 21 types of immune cells (Fig. 2B). The expression data of cancer cell lines from the CCLE database are listed in Fig. 2C. To take into consideration the inevitable batch effects between immune cells and cancer cell lines, we ranked the snoRNA expression profiles of immune cells and LUAD cell lines from high to low/ low to high mean expression of each snoRNA (Fig. 2B and C). The intersection between the most highly expressed snoRNAs among immune cells and the least expressed snoRNAs among LUAD cell lines were selected to be TIISRs, which included nine snoRNAs (SNORD14A, SNORD59A, SNORD63B, SNORD99, SNORD100, SNORD63, SNORD19, SNORA31, and SNORD13E). They were considered to be ubiquitously expressed in human immune cells, but not in LUAD. Further, we excluded snoRNAs what were not expressed in more than half of the samples from TCGA-LUAD cohort to improve the performance and stability of the RSF model.

Fig. 2

The expression of snoRNAs in immune and cancer cell lines. A. The human immune cell types screened by xCell. Color represents different types of immune cells. B. The expression of snoRNAs in immune cells from the Gene Expression Omnibus (GEO) database. From top to bottom, the expression decreases. Color represents different types of immune cells. C. The expression of snoRNA in cancer cell lines from the Cancer Cell Line Encyclopedia (CCLE) database. snoRNA expression increases from top to bottom. Color represents different types of cancer cell lines.

Finally, six snoRNAs (SNORD14A, SNORD59A, SNORD99, SNORD100, SNORD63, and SNORD19) were screened out as TIISRs. We separately identified the expression levels of these six TIISRs in different immune cells. We found that they were expressed in 21 types of purified human immune cells (Fig. S1A, S1B), suggesting that they were housekeeping for almost all kinds of immune cells.

Also, RNA from PBMC and four lung cancer cell lines were prepared for RT-qPCR assays to verify findings from bioinformatic analysis. As shown on Fig. 3, all six snoRNAs were expressed much higher in PBMCs than PC-9, A549, H1299 and H1975 cell lines.

Fig. 3

The expression of six snoRNAs from PBMCs and four lung cancer cell lines. A-F. The expression of SNORA14A, SNORD19, SNORD59A, SNORD63, SNORD99 and SNORD100 in PBMCs and PC9, A549, H1299, H1975 lung cancer cell lines. (Compared by t tests.)

Consequently, these six TIISRs were used for the subsequent analysis.

Identification and biological characteristics of TIISR signature

Using TCGA database, 479 samples containing complete clinical information and transcriptomic data of LUAD were collected as a training dataset. Then, based on the expression profile of these six TIISRs and the overall survival information (Fig. 4A), we identified the TIISR signature by RSF using TCGA-LUAD dataset. When the TIISR signature was applied to the training dataset, 479 LUAD patients were divided into the high-risk (n = 240) and low-risk groups (n = 239) by the median value of the risk score (0.5335).

Fig. 4

Identification and biological characteristics of the tumor-infiltrating immune cell-related snoRNA (TIISR) signature in The Cancer Genome Atlas (TCGA)-lung adenocarcinoma (LUAD) dataset. A. The distribution of the TIISR, overall survival (OS) status of patients, and snoRNA expression pattern. B. Gene set variation analysis (GSVA) enrichment analysis showed activated (red) or inhibited (blue) pathways between high and low risk score. Different colored squares represent different pathways. C. Survival analysis for risk score of TCGA-LUAD dataset. Blue and red lines represent low or high risk score, respectively. D. The predictive value of the TIISR signature in LUAD patients among TCGA-LUAD dataset. (area under the curve [AUC] 0.83, 3-year OS, and AUC 0.82, 5-year OS). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Firstly, GSVA enrichment analysis was performed to explore the biological characteristics of TIISR signature under different risk score. As shown in Fig. 4B and Table S3, the high-risk group was remarkably enriched in G2M checkpoint, E2F targets, mTORC1 signaling, glycolysis, MYC targets V1 and V2, protein secretion pathways, and biological states. Meanwhile, the low-risk group was mainly enriched in cholesterol homeostasis, PI3K/AKT/mTOR signaling pathways, and biological states. These results demonstrate that there were significant differences in the pathways and the biological status between two risk score groups in LUAD patients.

Furthermore, we extracted the expression of these six TIISRs under different risk score. As shown in Fig. 4A, the heatmap revealed that SNORD14A, SNORD63, SNORD59A, and SNORD100 had high expression in the low-risk group and low expression in the high-risk group. In addition, the high-risk individuals had poor prognosis, indicating that the TIISR signature might play a harmful role in LUAD. Thus, we further explored whether the TIISR signature could reflect the probability of patient survival. Survival analysis of the high/low-risk groups in TCGA-LUAD dataset (n = 479) demonstrated that the patients who had low risk score had a longer survival time than those with high risk score (Fig. 4C). Furthermore, the ROC of the 3-year and 5-year survival correlation of the TIISR signature was analyzed, and the area under the curve (AUC) was calculated. The AUC of the TIISR signature was 0.83 and 0.82 at 3 and 5 years of survival, respectively (Fig. 4D); Also, the sensitivity and specificity were 82.7%, 62.4% for 3-year-OS and 77.2%, 64.4% for 5-year-OS with median risk score as cutoff. These results demonstrated the effectiveness of the signature as a prognostic biomarker for predicting the survival status of LUAD patients.

In order to further validate the role of the TIISR signature in LUAD patient prognosis, multivariate Cox analysis was applied. The results are presented in Table 1. The HR of the high risk score versus low risk score for survival probability was 4.854 (95% CI 3.445–6.84, P < 0.001), suggesting that the LUAD patients who have a higher risk score are more likely to have a poor prognosis.

Table 1

Univariate and multivariate regression analysis of disease-free survival in each data set.

	Univariate analysis			Multivariate analysis
Variables	HR	95% CI	P Value	HR	95% CI	P Value
TCGA-LUAD
Age	1.155	0.855–1.56	0.348	1.227	0.907–1.659	0.185
Gender (Male/Female)	1.053	0.782–1.416	0.735	1.076	0.799–1.45	0.629
Stage (High/Low)	2.522	1.843–3.451	<0.001	2.147	1.564–2.947	<0.001
Risk Score (High/Low)	4.854	3.445–6.84	<0.001	4.605	3.259–6.508	<0.001
GSE81089
Age	1.07	0.582–1.965	0.829	1.025	0.553–1.9	0.938
Gender (Male/Female)	1.231	0.696–2.178	0.475	1.304	0.729–2.332	0.37
Stage (High/Low)	3.208	1.799–5.72	<0.001	3.159	1.771–5.634	<0.001
Risk Score (High/Low)	2.016	1.125–3.611	0.018	2.044	1.134–3.684	0.017

TCGA, the cancer genome atlas.

LUAD, lung adenocarcinoma.

HR, hazard ratio.

CI, confidence interval.

The above results revealed that the risk score groups showed remarkably different in multiple biological features. Furthermore, the TIISR signature can also predict the prognosis of LUAD patients.

Immunological characteristics of different groups based on the TIISR signature

To explore the potential biological processes related to immune and influenced by the TIISR signature, we first conducted GO enrichment analysis by performing ssGSEA algorithm on TCGA-LUAD cohort. The results demonstrated that a variety of immune-related biological processes were different between two risk score groups, especially in the somatic diversification of immune receptors via somatic mutation (Fig. 5A and Table S4). Also, we found that the risk score was statistically different among the various immune subtypes (Fig. 5B). Then, we conducted xCell and CIBERSORT-Abs analysis based on multiple immune subpopulations. By assessing the immune infiltration stage of the high-risk and low-risk groups, we aimed to further verify the specificity of the TIISR signature in response to tumor cell immune infiltration. There were significant differences in the immune infiltration status estimated by xCell between high-risk score and low-risk score groups. Twenty-one subsets were enriched in the low-risk group, while only T cell CD4 + Th2 and T cell CD4 + Th1 were enriched in the high-risk group. Ten subsets, including B cells, monocytes, and neutrophils, were not statistically different among the different risk groups (Fig. 5C). A similar result could be found at estimating infiltrated level of various kinds of immune cells by CIBERSORT-Abs algorithm (Fig. S2). These findings suggest that a higher risk score is correlated with less immune cell infiltration.

Fig. 5

The immunological characteristics of the tumor-infiltrating immune cell-related snoRNA (TIISR) signature. A. Gene Ontology enrichment analysis showed a significant difference in immune-related biological processes between high- and low-risk score groups. B. The differences of risk score in five immune subtypes from The Cancer Genome Atlas (TCGA)-lung adenocarcinoma (LUAD) cohort. C. Volcano plots for the enrichment of immune cell types for tumors with high and low risk score, calculated based on the NES score from the xCell algorithm. NES, enrichment of immune subpopulations. D. Correlation analysis between risk score and immune score. The two were linearly negatively correlated. Higher immune score was correlated with lower risk score. E. Immune score in the two risk score groups. Kruskal–Wallis tests were applied for testing statistical differences.

In addition, the xCell algorithm was used to estimate the immune cell infiltration score and examine the correlation between the immune and risk score. The results showed that the two were negatively correlated with each other. A higher immune score was correlated with a lower risk score, and vice versa (Fig. 5D, E). Therefore, we inferred that the TIISR signature might indicate the degree of immune cell infiltration.

Taken together, the above results demonstrated that the TIISR signature could greatly assess the degree of tumor immune cell infiltration, and a higher risk score indicated less immune cell infiltration.

The correlation between the TIISR signature and anti-tumor immunity or immunotherapy response

First, we analyzed the correlation between the risk score and the degree of immune cell infiltration and found that 16 types of immune cells had significant correlation with the risk score. Moreover, while CD4 + Th1 and Th2 cells were positively correlated with the risk score, the rest were negatively correlated with it. In addition, we analyzed the correlation between the expression levels of 18 immune checkpoints and the risk score. This analysis demonstrated a significant negative correlation between the risk score and the expression levels of nine immune checkpoints, including ADORA2A, BTLA, CD200R1, CD80 and other 5 checkpoints (Fig. 6A and Table S5).

Fig. 6

The correlation between risk score and immune status and immune checkpoints. A. Spearman correlation analysis between risk score with immune cells (left) and immune checkpoints (right). B. The difference in antitumor immune steps between the high- and low-risk groups among various cancer immunity cycles.

Then, we analyzed the antitumor immunity of different risk score groups by tumor immune cycle. The results are displayed in Fig. 6B. Among the various anticancer immunity steps, such as those of T cell recruitment, dendritic cell recruitment, Treg cell recruitment, and tumor immune cell infiltration, the anticancer immunity status was different between the high- and low-risk score groups. In particular, during tumor immune cell infiltration, the antitumor immunity of the low-risk group was more pronounced compared with that of the high-risk group. Accordingly, we speculate that a lower risk score indicates stronger antitumor immunity.

To further investigate the potential role of the TIISR signature in predicting ICI treatment response, we examined the expression levels of 18 potential tumor immune checkpoints between two risk score groups. Previous studies have shown that CTLA4, BTLA, CD200R1, and CD80 played important roles in the development and immunotherapy of solid tumors [36], [37], [38], [39]. As shown in Fig. 7A, patients with low risk score were more likely to express CTLA4 than those with high risk score, and a similar trend was observed for BTLA, CD200R1, CD80 and other four immune checkpoints (Fig. S3A). While there were no differences on the expression of other 10 immune checkpoints between the high- and low-risk score groups (Fig. S3B)

Fig. 7

The role of the tumor-infiltrating immune cell-related snoRNA (TIISR) signature in immune checkpoint inhibitor (ICI) treatment. A. The relative expression of four immune checkpoints between the high-risk score and low-risk score groups (Kruskal–Wallis test). Left to right, CTLA4, BTLA, CD200R1, and CD80. B. Kaplan–Meier survival curves of overall survival status among four patient groups stratified by the TIISR signature and different immune checkpoints. From left to right, CTLA4, BTLA, CD200R1, and CD80.

We further determined whether the expression of CTLA4, BTLA, CD200R1, or CD80 had an impact on the survival of patients. The 479 samples from TCGA-LUAD cohort were divided into 4 subgroups according to the risk score and immune checkpoint expression levels: the high-risk score and high immune checkpoint expression group, the high-risk score and low immune checkpoint expression group, the low risk score and high immune checkpoint expression group, and the low risk score and low immune checkpoint expression group. Then, the survival probability of LUAD patients was compared among these four subgroups with different risk score and immune checkpoint expression levels. The results were shown in Fig. 7B. When the survival rate of LUAD patients was higher in low risk score group, which was consistent with our previous results. These results demonstrated that the TIISR signature had discriminatory power in patients with similar levels of immune checkpoint expression.

Thus, the TIISR signature might be a potential positive marker for representing anti-tumor immunity status and predicting the response to ICI treatment.

Validation in independent datasets

To assess the effectiveness of the TIISR-based RSF model in other datasets and prevent model overfitting, we used an independent dataset (GSE81089) as a test set to validate its ability to evaluate LUAD prognosis. One hundred and eight LUAD patients were divided into low-risk (n = 54) and high-risk groups (n = 54). The results showed that the survival rate of the low-risk group was significantly higher than that of the high-risk group (HR 2.044, 95% CI 1.134–3.684, P = 0.017), and that the discrepancy between the long-term survival rates of the two groups widened as survival time increased. The 5-year survival rate of the high-risk group was only about two-thirds that of the low-risk group (Fig. 8A). Likewise, the AUC of the TIISR signature was 0.73 and 0.71 at 3 and 5 years of survival, respectively (Fig. 8B), supporting the notion that the TIISR signature is a powerful predictor of LUAD prognosis.

Fig. 8

Validation of the tumor-infiltrating immune cell-related snoRNA (TIISR) signature on an independent dataset (GSE81089). A. Survival analysis for risk score from the GSE81089 dataset. Blue and red lines represent low and high risk score, respectively. B. The predictive value of the TIISR signature in lung adenocarcinoma (LUAD) patients among the GSE81089 dataset (area under the curve [AUC] 0.73, 3-year overall survival [OS], and AUC 0.71, 5-year OS). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

As the LUAD cohort with ICIs treatment lacks the expression data of snoRNAs. We employed a dataset of bladder cancer from IMvigor210 cohort which contains expression information of snoRNAs and survival data. Patients on this cohort were received ICIs treatment. We further applied the RSF model and found that the high risk score group had worse prognosis (Fig. S4A, HR 1.313, 95% CI 1.093–1.577, P = 0.004) and response to ICIs treatment (Fig. S4B, P = 0.01549).

Discussion

LUAD is the most common histological type of lung cancer [40]; its degree of malignancy, tumor progression, and prognosis depend not only on the intrinsic characteristics of tumor cells, but also on tumor immune infiltration [41]. Recently, molecular map-based methods and signatures, such as DNA methylation [42] and lncRNA [43], have been applied to the inferences of tumor immune infiltration, and the role of snoRNA in controlling cell fate and tumorigenesis has been widely considered. The potential roles of different types of snoRNAs in tumors are gradually being revealed. For example, Grigory et al. [44] demonstrated that immune response-related genes were activated in cells transfected with snoRNA; meanwhile, SNORA50A/B might inhibit tumor progression [18]. In our study, six housekeeping snoRNAs were screened out: SNORA14A, SNORD63, SNORD59A, SNORD100, SNORD19, and SNORD99. Among these, SNORA59A has been revealed to play a role in prostate cancer [45], and SNORD63 has been considered as a noninvasive diagnostic biomarker for clear cell renal cell carcinoma [46]. Meanwhile, there were no relevant studies that demonstrated the application of the remaining four snoRNAs in cancers. Accordingly, based on the expression and TSI of snoRNAs in human immune cells and cancer cell lines, we herein identified the TIISR signature through RSF and then classified LUAD samples into high-risk and low-risk groups.

Studies have found that MYC and KRAS gene mutations and ALK gene rearrangements are pathophysiological alterations in lung cancer [47]. Among them, KRAS and MYC mutation-driven tumors are characterized by altered metabolic pathways, including enhanced nutrient intake, protein, and glycolysis [48], while the ALK gene activates mTORC1. Mutations in this pathway are extremely common in primary invasive cancers [47], [49]. In this study, mTORC1 signaling, glycolysis, MYC targets, protein secretion pathways, and biological status were significantly enriched in the high-risk group. Furthermore, our study also demonstrated that snoRNA expression was higher in the low-risk group compared to that in the high-risk group, which is consistent with the results of previous studies showing that snoRA50A/B may exert oncogenic effects by suppressing KRAS oncogenes [18]. The prognostic predictive roles of SNORD93 in breast cancer, SNORA42 in lung cancer, and SNORA21 in colon cancer have been confirmed and reported [12], [18]. Therefore, we further performed survival analysis and multifactorial Cox analysis in different risk score groups among the training set. We found that the low-risk group had longer survival time and better long-term prognosis. Further analysis revealed that the TIISR signature could be a valid marker for predicting LUAD patient prognosis; patients with a higher risk score tended to have worse prognosis.

To further confirm whether the predictive effect of the TIISR signature for LUAD prognosis is TIME-related, we performed an immune status analysis. The results showed that differential genes in dissimilar risk score groups were significantly enriched in a diverse range of immune-related biological processes. Immune cells such as NK cells, B cells, and macrophages have been demonstrated to have a positive or negative effect on cancer prognosis [50], [51], [52]. We therefore hypothesized that TIISRs might predict LUAD prognosis and reflect the degree of immune infiltration. Generally, active immunity is essential for the clearance of damaged cells. After stimulation by various cytokines or tumor cells, immune cells can target tumor cells and induce their death, thereby eliminating malignant cells and generating further immune responses [53]. Moreover, studies have indicated that the presence of infiltrating immune cells is associated with longer survival and better prognosis [54]. Similarly, we further evaluated the degree of immune infiltration under different risk score and found a negative correlation between them; greater immune infiltration was accompanied by a lower risk score and better prognosis. These results strongly demonstrate that the TIISR signature can predict the immune infiltration status and the prognosis of LUAD. More importantly, we performed validation on overall survival by using an independent dataset. The results were consistent with those displayed in the training set, suggesting that the prediction model had good generalizability in LUAD.

ICI treatment has led to substantial progress in the clinical treatment of a variety of solid tumors, including LUAD. Studies have reported that CTLA4, BTLA, CD200R1, and CD80 are important immune checkpoints for LUAD. CTLA4 is a member of the immunoglobulin superfamily, which encodes a transmembrane protein that negatively regulates T cell function and thereby acts as a tumor suppressor in LUAD [55]. BTLA, CD200R1, and CD80 are also considered promising targets for ICI treatment [38], [39], [56], [57]. However, as we all known, the expression of PD-1/PD-L1 were the most important biomarker for predicting response of ICIs treatment [55], [58]. BTLA, CD200R1, and CD80 are also considered promising targets for ICI treatment [38], [39], [56], [57]. However, as we all known, the expression of PD-1/PD-L1 were the most important biomarker for predicting response of ICIs treatment [58]. In this study, the TIISR signature showed no correlation with them which may considered less clinical significance nowadays. But, as shown in this study and many other papers, anti-tumor immunity was a complex system, many studies have revealed that PD-1/PD-L1 were not a perfect predictor and finding better predictors were never stopped in the past few years [59], [60]. In this study, we not only evaluated the relationship between TIISR signature and immune checkpoints, but also relationship between tumor microenvironment and cancer immunity circle. So, we think that our findings may become a possible solution for evaluate whether a patient with lung adenocarcinoma could benefit from a lot of developed or potential ICIs and offer a novel sight into evaluating anti-tumor immunity.

However, this study was far away from perfect. As the raw reads or eligible small RNA-Seq data were hardly obtained or not provided by GEO, CCLE and TCGA-data base and limited by our time and fundings. We downloaded level 3 data generated by common total RNA-Seq procedure for this study. Limited by theory of paired-end RNA-Seq pipeline, though part of snoRNAs could be detected in this type of data and finding could be validated by further biological experiments like this study and other works. The values of mature snoRNAs may be under represented and some related information will be lost by performing analysis using based on level 3 data directly. So, this study may be updated and extended when more appropriate data provided by databases or generated by researchers in the future.

Consistent with the results found in our study, a lower risk score was accompanied by a higher level of immune checkpoint expression. Subgroup analysis further suggested that LUAD patients with a low risk score and high immune checkpoint expression have better survival, also indicating that patients with a low risk score may benefit more from ICI treatment. Moreover, studies have demonstrated that ICI treatment is only effective in 10–20% of tumor patients [61]. This study showed that the TIISR signature can provide an important tool for screening for LUAD patients that are most likely to benefit from ICI treatment.

Conclusion

Taken together, this study screened immune-related snoRNAs by comparing the expression and TSI of immune cells and cancer cell lines. We then identified the TIISR signature by RSF model as a way to predict the immune infiltration of tumors and LUAD prognosis. Our findings show that the TIISR signature can provide an effective tool for identifying patients who might benefit from ICI treatment.

Funding

This work was supported by grants from the National Natural Science Foundation of China (81873406, 82170041 and 82100099), the Natural Science Foundation of Hunan Province (2019JJ50939) and Innovative Research Platform of Hunan Development and Reform Commission (2021-212).

Availability of data and materials

All data used in this work can be acquired from the Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) and the GDC portal (https://portal.gdc.cancer.gov/). Readers can download them directly by using GEOquery and TCGAbiolinks in R console or via following hyperlinks:

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE81089

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE114765

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE133145

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE135635

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE107011

https://portal.gdc.cancer.gov/projects/TCGA-LUAD

Also, the data from IMvigor210 cohort can be found at:

http://research-pub.gene.com/IMvigor210CoreBiologies/

And the source file of CCLE dataset can be downloaded via:

https://depmap.org/portal/download/

Ethics approval and consent to participate

The patient data used this study were acquired from publicly available datasets. The ethical approval for acquiring peripheral blood from volunteers was provided by Ethics Committee of the Xiangya Hospital of Central South University.

CRediT authorship contribution statement

Rongjun Wan: Methodology, Software, Data curation, Formal analysis, Writing – review & editing, Visualization. Lu Bai: Investigation, Data curation, Writing – original draft. Changjing Cai: Conceptualization, Methodology. Wang Ya: Methodology. Juan Jiang: Writing – review & editing. Chengping Hu: Supervision, Funding acquisition. Qiong Chen: Supervision. Bingrong Zhao: Project administration. Yuanyuan Li: Conceptualization, Supervision, Funding acquisition.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.