Identification of tumor immune infiltration-associated lncRNAs for improving prognosis and immunotherapy response of patients with non-small cell lung cancer.

BACKGROUND: Increasing evidence has demonstrated the functional relevance of long non-coding RNAs (lncRNAs) to immunity regulation and the tumor microenvironment in non-small cell lung cancer (NSCLC). However, tumor immune infiltration-associated lncRNAs and their value in improving clinical outcomes and immunotherapy remain largely unexplored.
METHODS: We developed a computational approach to identify an lncRNA signature (TILSig) as an indicator of immune cell infiltration in patients with NSCLC through integrative analysis for lncRNA, immune and clinical profiles of 115 immune cell lines, 187 NSCLC cell lines and 1533 patients with NSCLC. Then the influence of the TILSig on the prognosis and immunotherapy in NSCLC was comprehensively investigated.
RESULTS: Computational immune and lncRNA profiling analysis identified an lncRNA signature (TILSig) consisting of seven lncRNAs associated with tumor immune infiltration. The TILSig significantly stratified patients into the immune-cold group and immune-hot group in both training and validation cohorts. These immune-hot patients exhibit significantly improved survival outcome and greater immune cell infiltration compared with immune-cold patients. Multivariate analysis revealed that the TILSig is an independent predictive factor after adjusting for other clinical factors. Further analysis accounting for TILSig and immune checkpoint gene revealed that the TILSig has a discriminatory power in patients with similar expression levels of immune checkpoint genes and significantly prolonged survival was observed for patients with low TILSig and low immune checkpoint gene expression implying a better response to immune checkpoint inhibitor (ICI) immunotherapy.
CONCLUSIONS: Our finding demonstrated the importance and value of lncRNAs in evaluating the immune infiltrate of the tumor and highlighted the potential of lncRNA coupled with specific immune checkpoint factors as predictive biomarkers of ICI response to enable a more precise selection of patients. © Author(s) (or their employer(s)) 2020. Re-use permitted under CC BY-NC. No commercial re-use. See rights and permissions. Published by BMJ.

Introduction

Lung cancer is one of the most common cancers diagnosed in men and women and accounts for one-quarter of all cancer deaths.1 Non-small cell lung cancer (NSCLC) is the most frequent histological subtype and accounts for approximately 85% of all lung cancers.2 3 NSCLC has the lowest 5-year relative survival rate with approximately 19% partly because of diagnosed at a distant stage and a paucity of late-stage treatments.1 Therefore, continued research into molecular biomarkers and novel therapies is critical to predict prognosis and determine the personalized treatment.

The growing research on tumor microenvironment (TME) has indicated that tumor-infiltrating immune cells play a critical role in cancer progression and aggressiveness.4–6 There is evidence that the microenvironment of NSCLC is rich in different types of immune cells which are associated with clinical outcomes.7 8 Thus, the quantitative molecular signature of tumor-infiltrating immune cells is increasingly recognized as predictive biomarkers to enable personalized treatment selection and improve patient management. Long non-coding RNAs (lncRNAs), a type of non-coding RNAs longer than 200 nucleotides in length, have been shown to be involved in a wide range of biological and cellular functions.9–12 Recent studies have demonstrated that lncRNAs are emerging as critical regulatory elements in the immune system and play important roles in the development and differentiation of different immune cell lineages.13–15 A growing body of lncRNAs has been detected and identified in various immune cells.14 16 For example, Hu has identified 1524 lncRNAs in T cell populations.17 Another RNA-seq analysis of human lymphocyte subsets discovered more than 500 expressed lincRNAs and identified several lymphocyte subset-specific lincRNA signatures.18 A recent study also highlighted the functional relevance of lncRNAs in cancer immunity regulation and TME which contributes the progression and clinical outcomes of multiple cancers.19 Whole transcriptome RNA sequencing in diverse immune cell types of patients with melanoma identified 27,625 lncRNAs, some of which are significantly associated with melanoma status.20 Together these findings demonstrated the roles of lncRNAs in cancer immunology. Furthermore, the correlation between lncRNAs and immune cell infiltrate has also been observed in several cancers,21 22 which implied the potential of lncRNAs in evaluating the immune cell infiltrate of tumor.

In this study, we aimed to develop a novel computational frame for identifying tumor-infiltrating immune-related lncRNAs (TILncRNA) by integrative analysis for molecular profiling of purified immune cells, cancer cell lines and bulk cancer tissues and to explore their potential importance as predictive biomarkers for prognosis and immunotherapy in NSCLC tumors.

Materials and methods

Patient and tumor cell line cohorts

Clinical information and transcriptional profiles of patients with NSCLC were obtained from the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo) and The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov) according to following selection criteria: (1) have basic clinical information of stage, age, gender, overall survival (OS) and survival status; (2) based on the Affymetrix HG-U133_Plus 2.0 platform or IlluminaHiSeq platform; and (3) have larger sample size (>200). After filtering, a total of 1533 patients with NSCLC from three data sets were enrolled in this study, including 293 patients in GSE30219 from Rousseaux et al’s study (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE30219),23 226 patients in GSE31210 from Okayama et al’s study (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE31210)24 and 1014 patients from TCGA (https://portal.gdc.cancer.gov). GSE30219 data set was used as a training data set for discovering lncRNA signature, and other two data sets (GSE31210 and TCGA) were used as independent testing data sets for validating the lncRNA signature. The detailed clinical information of three patient sets was shown in table 1. Transcriptional profiles of 187 NSCLC cell lines based on Affymetrix HG-U133_Plus 2.0 platform were obtained from Cancer Cell Line Encyclopedia project (https://portals.broadinstitute.org/ccle).

Table 1

Clinical characteristics of patients with NSCLC in each data set

Characteristics	Training series	Testing series
	GSE30219	GSE31210	TCGA
Platform	HG-U133_Plus_2	HG-U133_Plus_2	IlluminaHiSeq
Patients, n	293	226	1014
>60	162	118	721
≤60	130	108	265
NA	1	0	28
Gender
Female	43	121	406
Male	250	105	608
Stage
I	152	168	518
II	77	58	283
III	56	0	168
IV	7	0	33
NA	1	0	12
Survival
Dead	200	35	284
Alive	93	191	730

NA, not applicable; NSCLC, non-small cell lung cancer; TCGA, The Cancer Genome Atlas.

Purified immune cell data

Transcriptional profiles of 115 purified cell lines of 19 immune cell types based on Affymetrix HG-U133_Plus 2.0 platform were obtained from GEO database including GSE13906, GSE23371, GSE25320, GSE27291, GSE27838, GSE28490, GSE28698, GSE28726, GSE37750, GSE39889, GSE42058, GSE49910, GSE51540, GSE59237, GSE6863 and GSE8059. Marker genes of tumor-infiltrating immune cell types that are representative and are not expressed in cancer cells or in normal tissues were obtained from Charoentong et al’s study25 which analyzed 366 microarrays of immune cells collected from 37 studies.

Preprocessing of transcriptional profiles

All raw data (.cel files) of microarray data sets profiled by Affymetrix HG-U133_Plus 2.0 platform were downloaded from the GEO database and were processed using the Robust Multi-array Average algorithm with R ‘affy’ packages from Bioconductor for background correction, quantile normalization, and log2 transformation.26 lncRNA expression profiles were obtained by reannotating the probes from the Affymetrix HG-U133_Plus 2.0 microarray data sets. By matching the NetAffx annotation files (HG-U133 Plus 2.0 Annotations, CSV format, release 36, July 12, 2016) of the probe sets and the annotation files of Refseq (release 79) and GENCODE (release 25), we extracted probe sets with Refseq IDs which were labeled as ‘NR_’ and annotated with ‘long non-coding RNA’ in Refseq database and with Ensembl gene IDs which were annotated as ‘long non-coding RNA’ in GENCODE project. Finally, we obtained 2466 unique lncRNAs corresponding to 3431 probe sets in microarray data sets for further analysis. Sequencing lncRNA expression profiles (IlluminaHiSeq platform) were obtained from TCGA (https://portal.gdc.cancer.gov).

Development of tumor-infiltrating immune-related lncRNA signature

We developed a novel computational frame for identifying tumor-infiltrating immune-related lncRNA signature (TILSig) by integrative immune and lncRNA profiling analysis of purified immune cells, cancer cell lines and bulk cancer tissues as follows (figure 1): (1) For each immune cell line, top 5% expressed lncRNAs were obtained as candidate immune-related lncRNAs. (2) The expression specificity of a candidate immune-related lncRNA with respect to different immune cell types was calculated using tissue specificity index (TSI) developed by Yanai 27 as follows:

Figure 1

Strategy for identifying tumor-infiltrating immune-related lncRNA signature (TILSig) in this study. Top 5% expressed long non-coding RNAs (lncRNAs) were obtained as candidate immune-related lncRNAs for each immune cell line. The specificity of expression of a candidate immune-related lncRNA with respect to different immune cell types was calculated using tissue specificity index (TSI). Those housekeeping lncRNAs (hklncRNAs) which are upregulated in immune cell lines and downregulated in non-small cell lung cancer (NSCLC) cell lines were selected as tumor-infiltrating immune-related lncRNAs. A prognostic signature by focusing tumor-infiltrating immune-related lncRNAs (TILSig) was constructed using the linear combination of the expression values of the prognostic tumor-infiltrating immune-related lncRNAs, weighted by their estimated regression coefficients in the multivariate Cox regression analysis.

where N is the total number of immune cell types and is the expression intensity of immune cell normalized by the maximal expression of any immune cell types for lncRNA. The TSI ranges from 0 to 1. The lncRNA is a general immune lncRNA when the value is 0, while the lncRNA is one immune cell-specific lncRNA when the value is 1. Those lncRNAs which are universally highly expressed in all immune cell types were defined as immune-related housekeeping lncRNAs (hklncRNA). (3) Those hklncRNAs which are upregulated in immune cell lines and downregulated in NSCLC cell lines were selected as TILncRNAs. (4) A prognostic signature by focusing TILncRNAs (TILSig) was constructed using the linear combination of the expression values of the prognostic TILncRNAs, weighted by their estimated regression coefficients in the multivariate Cox regression analysis.2 28

Statistical analysis

Differentially expressed lncRNAs between immune cell lines and lung cancer cell lines were determined using significance analysis of microarrays method with R packages ‘samr’. The prognostic value of NSCLC TILncRNAs was evaluated by univariate Cox proportional hazards regression analysis in the training data set. The Kaplan-Meier survival plots were used to estimate OS and disease-free survival (DFS), and the difference in OS or DFS between the high-risk group and the low-risk group was determined using log-rank tests with R package ‘survival’. Univariate and multivariate analyses with Cox proportional hazards regression for OS and DFS were performed on the individual clinical variables with and without TILSig in each data set. HRs and 95% CIs were calculated. The prognostic performance of TILSig was measured using time-dependent receiver operating characteristic (ROC) curves with R packages ‘survivalROC’. The predictive performance of TILSig for recurrence status was measured using ROC curves. The Kruskal-Wallis test was used for comparisons among multiple groups. All the statistical analyses were performed in R V.3.1.3 with additional Bioconductor packages

Gene set enrichment analysis

To evaluate the infiltration of immune cells, gene set enrichment analysis (GSEA) was performed using marker gene set of each immune cell with Bioconductor package ‘clusterProfiler’,29 which can identify immune cell types that are over-represented in high-risk group or low-risk group.25

Results

The landscape of lncRNA expression in human immune cells

To characterize the expression pattern of lncRNAs in various human immune cells, all lncRNAs for each immune type were first ranked based on their expression levels, and those RNAs whose expression levels were ranked top 5% were considered as candidate immune-related lncRNAs. Of them, 91 lncRNAs were commonly highly expressed in all 19 immune cell types and 117 lncRNAs were highly expressed in exclusively one type of 19 immune cell types (online supplementary table S1 and online supplementary figure S1). Then we calculated the TSI score for each of those 208 candidate immune-related lncRNAs to measure the expression specificity with respect to different immune cell types. Finally, we identified nine immune cell type-specific lncRNAs (icsLncRNA) which were expressed in only one type of immune cell and had a high score of cell-type specificity (>0.5) (online supplementary table S2 and online supplementary figure S2). Of them, three icsLncRNAs (LOC101926943, AP003774.1 and LINC00996) were found to be plasmacytoid dendritic cell specific and three icsLncRNAs (LINC01234, TLR8-AS1 and LINC01296) were significantly higher expressed in mast cells. For natural killer T cells and dendritic cells, we detected two specific icsLncRNAs (LINC00892 and LINC00515) and one icsLncRNA (LINC00158), respectively. In contrast, we also identified 57 hklncRNAs which are highly expressed in all immune cell types and have a lower score of cell-type specificity (<0.2), demonstrating their essential for basic immune cellular functions (online supplementary table S3).

Identification of NSCLC TILncRNAs

To identify NSCLC TILncRNAs, we only focused on these 57 hklncRNAs which are universally highly expressed in all immune cell types and are critical for the maintenance of basic immune cellular functions. We performed differential expression analysis of 57 hklncRNAs between immune cell lines and NSCLC cell lines, and identified 17 hklncRNAs which are upregulated in immune cell lines and downregulated in NSCLC cell lines demonstrating their expression specificity to immune cells rather than tumor cells. These 17 hklncRNAs were considered as NSCLC TILncRNAs.

Development of TILSig

We performed univariate Cox proportional hazards regression analysis to investigate the association between expressions of 17 TILncRNAs and patients’ OS time in the training data set. A set of seven TILncRNAs (HCG26, PSMB8-AS1, TNRC6C-AS1, CARD8-AS1, HCP5, LOC286437 and LINC02256) was significantly correlated with NSCLC patients’ OS and therefore was considered as prognostic immune lncRNAs. Then we develop an lncRNA-based prognostic signature indicative of immune infiltration (TILSig) using the expression of seven TILncRNAs weighted by the multivariate Cox regression coefficient as follows: TILSig=(−0.1323*expression value of HCG26)+(−0.2323*expression value of PSMB8-AS1)+(0.0009*expression value of TNRC6C-AS1)+(0.2472*expression value of CARD8-AS1)+(0.1190*expression value of HCP5)+(−0.3019−*expression value of LOC286437)+(−0.0572*expression value of LINC02256).

When TILSig was applied to the training data set, 293 patients were classified into the high-risk group (n=146) and low-risk group (n=147) using the median value (0.0429) of TILSig as risk cut-off. As shown in figure 2A, patients in the low-risk group had significantly longer OS time than those in the high-risk group (HR=1.981, 95% CI 1.491 to 2.633, log-rank p<0.001) (figure 2A). The 5-year survival rate of the low-risk group was 60.4% which is significantly higher than that of the high-risk group (37.9%). The area under the curve (AUC) of the TILSig was 0.665 and 0.646 at 5 and 3 years of OS (figure 2B). The distribution of TILSig and expression pattern of lncRNAs in the TILSig was revealed in figure 2C. To further examine whether the TILSig is highly reflective of immune infiltrates, we estimated the infiltration of 19 immune subpopulations for patients of high-risk and low-risk groups using single-sample GSEA analysis as previously described.25 We found that patient risk groups stratified by TILSig showed distinct immune infiltrate patterns. As shown in figure 2D, patients in the low-risk group were enriched with 10 immune subpopulations, while only four immune subpopulations were enriched in patients with high risk. These results suggested that the higher score of TILSig corresponded to less immune cell infiltration and poor outcome, the lower score of TILSig corresponded to greater immune cell infiltration and better outcome.

Figure 2

The tumor-infiltrating immune-related lncRNA signature (TILSig) is associated with outcome and immune cell infiltrates. (A) Kaplan-Meier survival curves of overall survival between patients with a higher score of TILSig and with the lower score of TILSig. (B) Time-dependent receiver operating characteristic (ROC) curve at 5 and 3 years of overall survival (OS). (C) The distribution of TILSig, patients’ survival status and long non-coding RNA (lncRNA) expression pattern. (D) Volcano plots for the enrichment of immune cell types for tumors with high TILSig and low TILSig calculated based on the NES score from the gene set enrichment analysis (GSEA). AUC, area under the curve; NES, normalized enrichment scores; NK, natural killer; NKT, natural killer T cell.

Confirmation of the TILSig in two independent data sets with different platform

To evaluate the robustness of the TILSig, we test its prognostic power using lncRNA expression and clinical data from an independent cohort of 226 patients profiled by microarray platform. The TILSig could stratify patients into low risk (n=127) and high risk (n=99) with significantly different OS (HR=2.422, 95% CI 1.219 to 4.812, log-rank p=0.009) using the same risk score derived from the training data set (figure 3A). The 5-year survival rate of the low-risk group is 90.5%, whereas the corresponding rate in the high-risk group is 75.5%.

Figure 3

Validation of the tumor-infiltrating immune-related lncRNA signature (TILSig) in independent cohorts. Kaplan-Meier survival curves of overall survival between patients with a higher score of TILSig and with the lower score of TILSig in the GSE31210 data set (A) and The Cancer Genome Atlas (TCGA) data set (B). (C) Distribution and comparison of the TILSig among five immune subtypes. (D) Expression pattern of four tumor-infiltrating immune-related lncRNAs (TILncRNAs) among five immune subtypes. FPKM, fragments per kilobase of exon model per million reads mapped; IFN-γ, interferon gamma; TGF-β, transforming growth factor beta.

The prognostic value of the TILSig was further tested using another completely independent TCGA cohort of 979 patients profiled by the IlluminaHiSeq platform. Due to different platforms, only four of seven lncRNAs in the TILSig were covered by the IlluminaHiSeq platform. Therefore, the TILSig only based on four lncRNAs (PSMB8-AS1, TNRC6C-AS1, HCP5 and CARD8-AS1) without re-estimating parameters was applied to TCGA data set. The median risk score cut-off point obtained from TCGA data set classified 979 patients into the high-risk group (n=490) and low-risk group (n=489). Survival analysis demonstrated that the high and low-risk groups are significantly different in OS (HR=1.305, 95% CI 1.028 to 1.656, log-rank p=0.028) (figure 3B). The 5-year survival rate of the low-risk group is 43.8%, whereas the corresponding rate in the high-risk group is 37.1%. We further examined the distribution of the TILSig among five immune subtypes reported by a recent study30 and found that there is significant difference in TILSig among five immune subtypes (Kruskal-Wallis test p<0.001) (figure 3C), suggesting that the TILSig is closely associated with tumor immune microenvironment. Furthermore, further analysis revealed that expression levels of four TILncRNAs in the TILSig are also significantly different among five immune subtypes (figure 3D). As shown in figure 3D, lncRNAs TNRC6C-AS1 and CARD8-AS1 tended to be overexpressed in transforming growth factor beta (TGF-β) dominant immune subtype than other four immune subtypes, while lncRNA HCP5 has higher expression levels in interferon gamma (IFN-γ) dominant compared with other four immune subtypes. LncRNA PSMB8-AS1 tended to be overexpressed in three of five immune subtypes, including IFN-γ dominant, inflammatory and TGF-β dominant. We further examine the association of the TILSig with OS in other 17 solid cancers using univariate Cox analysis. Pan-cancer analysis of the TILSig revealed significant association between the TILSig and OS in other two cancers: colon adenocarcinoma (COAD) and kidney renal clear cell carcinoma (KIRC) (online supplementary figure S3A). The TILSig could stratify patients into two groups with significantly different OS in COAD and KIRC (log-rank p=0.019 for COAD and p=0.0006 for KIRC) (online supplementary figure S3B).

Independence of the TILSig from other clinical factors

To examine whether the prognostic performance of the TILSig is independent of other clinical factors, we used multivariate Cox analysis to test the performance of the TILSig after adjusted by other clinical factors, including age, gender and stage. In multivariate analysis, the HRs of high TILSig versus low TILSig for OS were 1.584 (p=0.002; 95% CI 1.184 to 2.118) in the training set, 2.322 (p=0.019; 95% CI 1.147 to 4.703) in the GSE31210 test set and 1.338 (p=0.017; 95% CI 1.054 to 1.698) in the TCGA test set (table 2), respectively, showing that TILSig is still significantly correlated with unfavorable OS in the training set and other two independent test sets after adjusting for other clinical factors. The results of the multivariable analysis thus indicated that the prognostic performance of TILSig is independent of other clinical factors for OS prediction.

Table 2

Univariate and multivariate Cox regression analysis of overall survival in each data set

Variables	Univariate analysis			Multivariate analysis
	HR	95% CI of HR	P value	HR	95% CI of HR	P value
GSE30219
TILSig (high/low)	1.981	1.491 to 2.633	<0.001	1.584	1.184 to 2.118	0.002
Age	1.038	1.024 to 1.053	<0.001	1.037	1.022 to 1.052	<0.001
Gender (male/female)	1.701	1.082 to 2.674	0.022	1.398	0.886 to 2.206	0.150
Stage (III/IV vs I/II)	3.036	2.234 to 4.127	<0.001	2.810	2.059 to 3.835	<0.001
GSE31210
TILSig (high/low)	2.422	1.219 to 4.812	0.012	2.322	1.147 to 4.703	0.019
Age	1.025	0.978 to 1.075	0.306	1.035	0.988 to 1.085	0.147
Gender (male/female)	1.519	0.780 to 2.955	0.219	1.100	0.549 to 2.201	0.789
Stage (II /I)	4.232	2.175 to 8.236	<0.001	4.290	2.164 to 8.502	<0.001
TCGA
TILSig (high/low)	1.305	1.028 to 1.656	0.028	1.338	1.054 to 1.698	0.017
Age	1.017	1.003 to 1.030	0.014	1.020	1.006 to 1.034	0.004
Gender (male/female)	0.995	0.780 to 1.270	0.970	0.966	0.756 to 1.234	0.782
Stage (III/IV vs I/II)	1.901	1.476 to 2.45	<0.001	2.008	1.555 to 2.592	<0.001

TCGA, The Cancer Genome Atlas; TILSig, tumor-infiltrating immune-related lncRNA signature.

The TILSig is associated with tumor recurrence

To examine whether the TILSig is also related to tumor recurrence, we first used DFS time as recurrence measures to investigate the effectiveness of the TILSig for recurrence risk prediction in the training set and independent GSE31210 test set. As shown in figure 4A, C, survival analyses demonstrated that the DFS distribution from the high-risk group was significantly different from that of the low-risk group (log-rank p=0.0001 for the training set and p=0.01 for GSE31210 test set). In the high-risk group, the proportions of the disease-free patient were 47.6% and 60.3% at 5 years in the training and GSE31210 sets, respectively, which were also significantly lower than those in the low-risk group (69.3% and 77.1%) (figure 4A, C). Next, we used TILSig to predict patients’ recurrence status, which achieved better predictive performance with AUC of 0.607 in the training set and 0.618 in the GSE31210 test set (figure 4B, D). Moreover, increasing TILSig score is associated with a greater probability of disease recurrence. We also performed multivariate Cox analysis of DFS for each data set and found that the TILSig still maintained an independent correlation with DFS (p=0.024, HR=1.570, 95% CI 1.062 to 2.320 for training set and p=0.034, HR=1.728, 95 % CI 1.043 to 2.863) (table 3).

Figure 4

The tumor-infiltrating immune-related lncRNA signature (TILSig) is associated with tumor recurrence. Kaplan-Meier survival curves of disease-free survival between patients with a higher score of TILSig and with the lower score of TILSig in the training data set (A) and GSE31210 data set (C). Receiver operating characteristic (ROC) curve of the TILSig for predicting tumor recurrence in the training data set (B) and GSE31210 data set (D). AUC, area under the curve.

Table 3

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

Variables	Univariate analysis			Multivariate analysis
	HR	95% CI of HR	P value	HR	95% CI of HR	P value
GSE30219
TILSig (high/low)	2.000	1.373 to 2.914	<0.001	1.570	1.062 to 2.320	0.024
Age	1.021	1.004 to 1.038	0.017	1.019	1.000 to 1.037	0.045
Gender (male/female)	1.355	0.787 to 2.333	0.273	1.216	0.704 to 2.100	0.484
Stage (III/IV vs I/II)	3.959	2.685 to 5.837	<0.001	3.583	2.411 to 5.326	<0.001
GSE31210
TILSig (high/low)	1.863	1.136 to 3.053	0.014	1.728	1.043 to 2.863	0.034
Age	1.034	0.997 to 1.072	0.074	1.039	1.003 to 1.077	0.0342
Gender (male/female)	1.271	0.778 to 2.075	0.338	1.012	0.609 to 1.682	0.964
Stage (II /I)	3.164	1.921 to 5.213	<0.001	3.280	1.968 to 5.466	<0.001

TILSig, tumor-infiltrating immune-related lncRNA signature.

Potential of TILSig as an indicator of immunotherapy response in patients with NSCLC

The previous study has revealed that immune infiltration can be modulated by the immune checkpoint inhibitor (ICI) genes. To further investigate the complex crosstalk between immune infiltration and ICI genes, we first compared the expression pattern of ICI genes (PD1, PD-L1 and CTLA-4) between different patient groups stratified by the TILSig. As shown in figure 5A, patients with high TILSig tend to express high ICI genes compared with those with low TILSig in the GSE20319 data set (Mann-Whitney U test p=0.006, p<0.001, p=0.007 for PD1, PD-L1 and CTLA-4, respectively; figure 5A). Similar trends also were observed both in GSE31210 (Mann-Whitney U test p=0.068, p<0.001, p<0.001 for PD1, PD-L1 and CTLA-4, respectively; figure 5B) and TCGA data sets (Mann-Whitney U test p<0.001, p<0.001, p<0.001 for PD1, PD-L1 and CTLA-4, respectively; figure 5C). This trend is concordant with previous observations linking higher expression of immune checkpoint genes to poor outcomes. We further examined whether the immune infiltration has an impact on clinical outcomes in patients with similar expression levels of immune checkpoint genes. Survival distribution of four patient groups stratified by the TILSig and high/low immune checkpoint gene expression was compared. As shown in figure 5D, patients with low TILSig and high PD-1 have significantly better survival than those with high TILSig and high PD-1 (log-rank p<0.001), and patients with low TILSig and low PD-1 also are associated with prolonged survival relative to those with high TILSig and low PD-1 (log-rank p<0.001) (figure 5D). A similar analysis was repeated using TILSig and PD-L1 or CTLA-4, and revealed that patients stratified by TILSig and PD-L1 or CTLA-4 exhibited survival patterns similar to those of PD-1 (figure 5E, F). When stratifying patients with low TILSig based on immune checkpoint gene expression, we also observed that immune checkpoint gene expressions are associated with noticeable survival differences in patients with low TILSig. However, survival differences were not present for patients with high TILSig when stratifying patients with TILSig-based immune checkpoint gene expression. Furthermore, patients with low TILSig and low immune checkpoint gene expression tend to have significantly better survival than other three patient groups. These results indicated that the TILSig may be a potential predictive biomarker of treatment response to ICI immunotherapy.

Figure 5

Impact of immune checkpoint gene expression and tumor-infiltrating immune-related lncRNA signature (TILSig) on clinical outcome. Comparison of the expression pattern of immune checkpoint genes (PD1, PD-L1, and CTLA-4) between patients with a higher score of TILSig and with lower score of TILSig in the training data set (A), GSE31210 data set (B) and TCGA data set (C). Kaplan-Meier survival curves of overall survival among four patient groups stratified by the TILSig and PD-1 (D), PD-L1 (E) and CTLA-4 (F). TCGA, The Cancer Genome Atlas.

Discussion

The immune infiltrates in the TME are increasingly recognized to be associated with patient prognosis in NSCLC and other cancers.5 31 32 For example, the immunoscore was proposed based on the mean of four density percentiles including two markers (CD3+ and CD8+) and two regions (tumor and invasive margin regions), and has been internationally validated as a risk assessment tool in colon cancer,33 highlighting the potential importance of evaluating the immune infiltrate of tumor in guiding decision-making in the clinic. Therefore, quantitative evaluation of the tumor immune infiltrates in the TME is a major challenge. Considerable efforts have been devoted to quantifying the composition of immune cell infiltration. However, the traditional immunohistochemistry immunoscoring approach remains suboptimal because of lacking a consistent standard and limited number of biomarkers assessed simultaneously.34 Recently, molecular profiling-based approaches and signatures have been used to infer immune infiltration, such as mRNAs,31 microRNA35 and DNA methylation.36 lncRNA, a recently identified key players of genome regulatory network, has also been found to play critical roles in the development and activation of immune cells, and may, therefore, serve as a specific molecular marker for tumor-infiltrating immune cells.14 15 However, genome-wide screening of tumor immune infiltration-associated lncRNA and their value in evaluating immune infiltrate of tumor and clinical outcome has barely been explored.

In this study, we first reannotated the expression profiles of 19 immune cell types and uncovered previously uncharacterized immune-related lncRNAs including nine icsLncRNAs and 57 hklncRNAs. Gene Ontology (GO) enrichment analysis suggested that mRNAs coexpressed with hklncRNAs were significantly enriched in various immune-related biological processes (online supplementary figure S4). Based on the assumption that lncRNAs that have universal expression across different immune cell types may be essential for basic immune cellular functions, we developed a novel computational frame by integrative analysis of molecular profiling of purified immune cells, cancer cell lines, and bulk cancer tissues to identify TILncRNAs for evaluating immune cell infiltration in TME. By applying this computational frame to NSCLC cell line and patient data sets, we identified the lncRNA signature, named TILSig, as an indicator of immune cell infiltration in patients with NSCLC. This computational frame for identifying TILSig may be referenced for other cancer types. The ability of the TILSig to evaluate immune cell infiltration is based on highly differential expression of hklncRNAs in immune cells versus lung cancer cell lines. In particular, some components of the TILSig have recently been associated with immunological functions and cancer prognosis. LncRNA HCG26 is located in the human leucocyte antigen (HLA) class I region which is a crucial determinant of the adaptive immune response, and recent study also linked HCG26 with bone metastasis of NSCLC.37 CARD8-AS1 has recently been shown to be involved in the regulation of immune response and epithelial-mesenchymal transition.38 HCP5, a hybrid HLA class I endogenous retroviral gene, is known for its functional roles in adaptive and innate immune responses, and newly is recognized as SMAD3-responsive lncRNA which can promote lung adenocarcinoma metastasis via miR-203/SNAI axis.39 These published experimental efforts provide further evidence supporting the association between the TILSig and immune responses.

To further understand the prognostic role of the TILSig, we then correlate the TILSig with patient survival, resulting in significant risk stratification of survival in patients with NSCLC (figure 2A). Using the single-sample GSEA analysis to measure the infiltration of 19 immune subpopulations within NSCLC tumors of different risk groups stratified by the TILSig, respectively, we found significantly different immune contexture between different risk groups classified by the TILSig. The high-risk group defined by the TILSig seems to be an immune-cold patient group with less immune cell infiltration and low-risk group defined by the TILSig is an immune-hot patient group with greater immune cell infiltration. This observation is concordant with previous findings that there are two distinct immunophenotypes in patients with NSCLC.40 These results showed that the TILSig is highly reflective of immune cell infiltrates resulting in a better prognosis prediction of patients with NSCLC. To confirm that these associations were not specific to one data set, we applied the TILSig to other two additional larger NSCLC patient data sets from GEO and TCGA. The TILSig is reproducible across two independent patient data sets with microarray or RNA-seq platforms and was capable of distinguishing patients with a different outcome in terms of survival and recurrence.

Over the past several years, checkpoint inhibitor-based immunotherapies have revealed substantial advancement in clinical care for many cancer types, including NSCLC.41 42 An early assessment for ICI response by predictive biomarkers is crucial for the selection of patients who are most likely to benefit from ICI. The expression of immune checkpoint genes such as PD-L1 is currently available biomarkers in clinical practice, but not a sufficient independent predictor of ICI response.43 44 Patients with NSCLC still have heterogeneous response rates to immunotherapy despite high tumor PD-L1/PD-1 levels and only a small subset of patients with NSCLC are responders for immune checkpoint blockade therapy.45 The immune infiltrates in the TME are increasingly recognized to be associated with immunotherapy response. By comparing the survival distributions of patients stratified by the TILSig and immune checkpoint gene expression, we found that the TILSig has a discriminatory power in patients with similar expression levels of immune checkpoint genes, and the prognosis of patients with similar TILSig could be affected by immune checkpoint gene expression. These findings suggest that the complex interplay between immune infiltrate and immune checkpoint gene in the tumor immune microenvironment has an impact on patient survival in NSCLC which is in line with previous study in immunomodulatory interactions between immune cell infiltration and checkpoint gene expression.31 Furthermore, significantly prolonged survival was observed for patients with low TILSig and low immune checkpoint gene expression, implying that these tumors with low TILSig may be associated with better response to ICI therapy.

Conclusions

We performed a genome-wide screening by integrative analysis of molecular profiling of immune cells, cancer cell lines and bulk cancer tissues to identify lncRNAs associated with tumor immune cell infiltration and highlighted the importance and value of lncRNAs in evaluating the immune infiltrate of the tumor. Furthermore, for the first time, our study identified and validated an lncRNA signature (TILSig) that is based on seven TILncRNAs, as an indicator of immune cell infiltration in the TME and has independent prognostic significance for patients with NSCLC. Lastly, we characterized complex interplay between immune infiltrate and immune checkpoint gene in patient’s outcome, and suggested the potential of TILSig coupled with specific immune checkpoint factors as predictive biomarkers of ICI response to enable a more precise selection of patients who will benefit from checkpoint inhibitor immunotherapy which need further validating with continuously releasing immunotherapy data sets.