Construction of a Novel Lung Adenocarcinoma Immune-Related lncRNA Pair Signature

Introduction

Lung cancer is a major cause of cancer-related deaths worldwide. The 5-year survival rate for lung cancer patients is below 20% in most countries, accounting for approximately 2.2 million new cases and 1.8 million deaths in 2020.¹ Lung adenocarcinoma (LUAD) is the most common historical subtype of lung cancer - accounting for more than 40% of the cases.² Immunotherapy, especially Immune Checkpoint Inhibitors (ICIs), has been regarded as a promising therapy for LUAD patients, who do not harbor a known mutation, such as in the epidermal growth factor receptor gene. For example, Pembrolizumab and Nivolumab - anti-programmed cell death-1 (PD-1) antibodies - are approved for advanced non-small cell lung cancer treatment.^3,4 Also, another anti-programmed cell death ligand-1 (PD-L1) antibody, Atezolizumab is approved.⁵ Additionally, clinical trials such as CHIO3, to evaluate the response of operable stage III A/B non-small cell lung cancer patients to Durvalumab, are in progress.

Long non-coding RNAs (lncRNAs) are a type of RNA with a length greater than 200 nucleotides. Although they are not translated into proteins, lncRNAs play an essential role in gene transcription and expression. Studies have suggested that lncRNAs are associated with many tumor-related biological processes. For example, overexpression of LINC00472 was proved to inhibit Epithelial-Mesenchymal Transition (EMT) processes of A549 and PC-9 cells via YBX1 modulation, a known cancer-promoting gene.⁶ Additionally, a growing number of studies showed that lncRNAs are related to many aspects of cancer-immune activities including immune infiltrating cells with an identified role in immunotherapy responses.⁷ For example, in patients with non-small cell lung cancer, a higher expression level of PD-1 may be accompanied by CD8⁺TILs depletion.⁸ Therefore, a novel understanding of the relationship between lncRNA and tumor-infiltrating immune cells deserves to be explored.

Several immune-related lncRNAs (irlncRNAs) LUAD signatures have been established and demonstrate potential diagnostic and predictive value. The Areas Under Curve (AUC) values of Receiver Operating Characteristic (ROC) curve calculated for the LUAD dataset of The Cancer Genome Atlas database (TCGA, https://tcga-data.nci.nih.gov/tcga/) ranged from 0.686 to 0.783 in the established signature. All signature were constructed based on lncRNAs’ expression levels,^9–14 and two research analyzed the relationship between infiltration immune cell subtypes and risk score.^13,15

A two-gene combination strategy to seek cancer-related biomarkers from gene expression big data was identified as more sensitive in multiple cancers.¹⁶ Hong et al¹⁷ constructed a novel irlncRNA signature for hepatocellular carcinoma with paired irlncRNAs, and the AUC of the 5-year ROC curve in their model was 0.904, which showed the combination strategy potential to construct prognostic signatures. However, few studies have applied this strategy to LUAD.

Materials and Methods

Identification of Immune-Related lncRNA

The transcriptome profiling (RNAseq) fragments per kilobase million (FPKM) data from the LUAD project were downloaded from the TCGA database. The expression profile data of each sample was combined in an RNAseq expression matrix. Gene Transfer Format (GTF) files, downloaded from Ensemble database (http://asia.ensembl.org) for annotation, were used to isolate lncRNAs from the RNAseq expression matrix. The list of immune genes was obtained from the ImmPort database (http://www.immport.org) and an immune genes-related expression matrix was also extracted from the RNAseq expression matrix. Immune-related lncRNAs (irlncRNAs) were screened by a correlation between immune genes and lncRNAs. Those with correlation coefficients > 0.4 and p < 0.001 were recognized as irlncRNAs. We identified differential expression of irlncRNAs (DEirlncRNAs) using the R limma package and set the log Fold Change (FC) > 1.5 along with a False Discovery Rate (FDR) < 0.05 to screen irlncRNAs. A DEirlncRNAs heatmap was drawn with the R pheatmap package.

Construction of the DEirlncRNA Pairs Matrix

The construction of the DEirlncRNA pairs matrix was performed based on each lncRNAs pair expression level. If the expression level of lncRNA A was higher than that of lncRNA B, it was recorded in the matrix as 1, otherwise, as 0. In the 0-or-1 matrix constructed, pairs without a certain rank were not appropriate for predicting patients’ survival outcomes. Finally, when the number of lncRNA pairs with an expression of 0, or 1, accounted for less than 20% of the total pairs, they were considered an invalid match.

Clinical Data Acquisition

Clinical data of LUAD patients were downloaded from the TCGA-LUAD project. Patients without documented survival, or survival status, were excluded.

Establishment of DEirlncRNAs RiskScore Model

First, we conducted a univariate Cox regression analysis with a p-value of 0.05 to screen DEirlncRNA pairs candidates with prognostic value. Then, we performed a Least Absolute Shrinkage and Selection Operator (LASSO) regression with 10-fold cross-validation. To improve model stability, the LASSO procedure was repeated 1000 times with a random stimulation, and pairs with frequency superior to 100 times were reserved for stepwise Cox proportional hazard regression analysis. A multivariate Cox regression analysis was also conducted. The ROC curve analysis for 1-, 2-, 3-, 4- and 5-years was plotted and the AUC values were calculated. RiskScore calculation for each patient was conducted as follows:

The “βi” and “Si” represent each prognostic DEirlncRNA pair coefficient and expression status, respectively. Patients were divided into high- and low-risk groups according to the optimal Akaike Information Criterion (AIC) values, where the maximum inflection point of the 4-year ROC curve was regarded as the risk cut-off point. The R survival, glmnet, survivalROC, and survminer packages were utilized for these operations.

Constructed Risk Model Clinical Prognostic Signatures

To valid the cut-off point and clarify the survival differences between high- and low-risk groups, we carried out a Kaplan–Meier analysis, and the survival curve was plotted. ROC curves for 4-year as well as for clinical characteristics such as the TNM stage, age, and gender were visualized. Then, we evaluated the association between the risk model and clinical characteristics with a χ² or a Wilcoxon signed-rank test. To better explore if the risk model can be regarded as an independent prognostic predictor, we used univariate and multivariate Cox regression analyses. The following R packages were applied in these operations: survival, glmnet, pbapply, survivalROC, survminer, pHeatmap, and ggupbr.

Tumor-Infiltrating Immune Cells Estimation

We calculated infiltration values of immune cell subtypes for TCGA-LUAD dataset samples based on 7 algorithms: TIMER,^18,19 QUANTISEQ,^20,21 EPIC,²² CIBERSORT,²³ XCELL,^24,25 MCPcounter,²⁶ and CIBERSORT-ABS.²⁷ A Wilcoxon signed-rank test was used to explore infiltration value differences between high- and low-risk groups, and a p < 0.05 was considered to be significant. The relationship of infiltration immune cell subtypes and RiskScore value was assessed by a Spearman correlation analysis. The following R packages were applied in these operations: limma, ggplot2, ggtext, scales, and ggpubr.

Immunosuppressive Molecules Expression

We compared the gene expression levels of CTLA4, PDCD1 that were related to ICIs between high- and low-risk groups in our risk model. We used limma and ggpubr R packages in that case.

Chemotherapeutic and Targeted Therapy Response Prediction

Chemotherapeutic and targeted therapy is commonly used in LUAD, so we decided to use the pRRophetic R package to calculate the half-maximal Inhibitory Concentration (IC₅₀) of anticancer drugs (gefitinib, erlotinib docetaxel, cisplatin, vinorelbine, and paclitaxel) for each TCGA-LUAD dataset sample. Gene expression and drug sensitivity data were downloaded from the Genomics of Drug Sensitivity in Cancer website (GDSC, https://www.cancerrxgene.org). Additionally, we compared the IC₅₀ difference between high- and low-risk groups in our risk model by a Wilcoxon signed-rank test.

Results

Verification of DEirlncRNAs in TCGA-LUAD Database Samples

The research roadmap of our study is represented in Figure 1. A total of 594 transcriptome profiling data - 535 cancer and 59 normal samples - were obtained from the TCGA-LUAD project. A summary list containing 2483 immune genes was downloaded from the ImmPort database. After annotation of the transcriptome profiling data with GTF files, an expression matrix containing 14,087 lncRNAs was extracted. Then, the correlation analysis was conducted between lncRNAs and immune genes, 1201 irlncRNAs were identified according to the pre-defined filtering conditions (Supplementary Table S1), 293 of which were DEirlncRNAs with a log FC > 1.5 and FDR < 0.05, including 249 upregulated and 44 downregulated. The DEirlncRNAs volcano map and heatmap are shown in Figure 2.

Figure 1 Study’s roadmap.

Figure 2 Identification of DEirlncRNAs: (A) heatmap and (B) volcano map of 293 DEirlncRNAs.

RiskScore Model Establishment and Validation

A total of 31,914 DEirlncRNA pairs were obtained by comparing the expression level of each DEirlncRNA pair. After univariate Cox and LASSO regressions, eleven pairs retained prognostic significance were included in the risk model (Figure 3), among which 7 DEirlncRNAs had hazard ratio (HR) values greater than 1, and the remaining 4 including LINC02036|AL365181.2, AC026355|LINC01547, AC104958.2|AC031985.3 and AC007991.2|AL365181.2 were less than 1. In the multivariate Cox regression model, seven DEirlncRNAs pairs were further confirmed to regarded as independent prognostic factors. Then, the ROC curve analysis for 1-, 2-, 3-, 4- and 5-years was plotted and the AUC values were calculated (Figure 4A and B). All AUC values were over 0.7 and the 4-year ROC curve presented the maximum (0.819). Additionally, we compared the 4-year ROC curve with common clinical characteristics ROC curve (Figure 4D). It is apparent that the prognostic value of our RiskScore model is superior to the clinical characteristics, including age, gender, sex, tumor stage and TNM stage, which all have AUC values below 0.7. According to the AIC, we identified 1.003 as the cut-off point (Figure 4C), and patients in TCGA-LUAD with a risk score higher than 1.003 were classified as the high-risk group and those above 1.003 as the low-risk group.

Figure 3 Construction of the irlncRNA signature: (A) partial likelihood deviance was plotted against the lambda logarithm in the 10-fold cross-validation. (B) LASSO coefficient profiles. Forest maps of the 11 DEirlncRNA pairs identified by (C) univariate and (D) multivariate Cox regression analysis.

Figure 4 ROC curves of the irlncRNA signature: (A) ROC curve with the optimal AUC value within 5 years. (B) The 1-, 2-, 3-, 4- and 5-year ROC curves. (C) The cut-off point is calculated by AIC. (D) A comparison between the 4-year and clinical characteristics ROC curves.

Constructed Model Effectiveness Validation

Figure 5A and B shows the distribution of risk scores and survival status of patients in the high- and low-risk groups and suggested that patients in the low-risk have a higher survival rate. A Kaplan–Meier analysis verified that patients in the low-risk group exhibited a better prognosis than those in the high-risk one (p < 0.001) (Figure 5C). The relationships between clinical characteristics (TNM stage, tumor stage, age, and gender) and risk score are presented in Figure 6. TNM and tumor stage were significantly related to the risk score. Patients with advanced lung adenocarcinoma had significantly higher risk score compared to early stage. To further explore the prognostic value of the constructed model, a univariate and multivariate Cox regression analysis were conducted and the results showed that tumor stage (p < 0.001, HR = 1.577, 95% CI [1.348–1.845]), T stage (p < 0.001, HR = 1.579, 95% CI [1.296–1.923]), M stage (p = 0.037, HR = 1.843, 95% CI [1.038–3.272]), N stage (p < 0.001, HR = 1.706, 95% CI [1.405–2.072]) and riskScore (p < 0.001, HR = 1.363, 95% CI [1.291–1.438]) can be regarded as prognostic factors, while the riskScore (p < 0.001, HR = 1.334, 95% CI [1.257–1.417]) was the one that had an independent prognostic role.

Figure 5 Survival characteristics of patients with lung adenocarcinoma: (A) risk scores distribution. (B) Survival outcomes distribution. (C) Survival curve in low- and high-risk groups.

Figure 6 Clinical characteristics evaluation by the risk model: (A) a strip chart aggregating common clinical characteristics. Associations between: (B) T-stage and risk scores; (C) N-stage and risk scores; (D) M-stage and risk scores; (E) clinical stage and risk scores. (F) Univariate and (G) multivariate Cox regression analysis of risk score and clinical characteristics. ** represented p < 0.01 and *** represented p < 0.001.

Comparison with Other Existing Prognostic Signatures for Lung Adenocarcinoma

We have compared our risk model with other existing immune-related lncRNAs prognostic model for lung adenocarcinoma based on TCGA dataset. In a signature consisting of 7 immune-related lncRNAs, the ROC curve for 5-year overall survival (OS) corresponded to an AUC value of 0.702.¹² Time-ROC curve analysis of a 10 immune-related lncRNAs signature in the TCGA dataset in 5-year was 0.783.⁹ The AUC value of ROC curve for a 4 immune-related lncRNAs prognostic signature, consisted of ITGB1-DT, ABALON, TMPO-AS1 and VIM-AS1, plotted by Zhang et al was 0.756¹⁴ and The 5-year AUC value in another 4 immune-related lncRNAs prognostic signature including HSPC078, DRAIC, AP004608.1 and MIR223HG was 0.63.¹³ The 1-year AUC value in a 5 immune-related lncRNAs signature based on TCGA was 0.686.¹¹ The AUC of the 6 immune-related lncRNA signature of Miao et al was 0.761 in training cohort and 0.723 in testing cohort.¹⁵ In our constructed signature, the highest AUC value was 0.819, and the AUC value of the 5-year ROC curve was 0.791. Comparing with other existing prognostic models, our signature has superior predictive sensitivity and specificity.

Constructed Model Immune Landscape

We presented the differences in tumor-infiltrating immune cells between high- and low-risk groups, based on 7 algorithms, with a boxplot (Supplementary Figure S1). All differences were significant and we could infer that B cells, T cells CD8, and T cell CD4 memory resting were tested for high expression in the low-risk group (Supplementary Table S2). Additionally, a Spearman correlation analysis was further conducted (Figure 7A) and the detailed results are listed in Supplementary Table S3. B cell, CD8⁺ T cell, Macrophage M2, and T CD4⁺ T cell [except for T cell CD4⁺ (non-regulatory)], T cell CD4⁺ memory activated, T cell CD4⁺ Th2 were negatively correlated to risk scores, while macrophage M0 and cancer-associated fibroblast were positively correlated. Tow ICIs-related genes were tested in the constructed model (Figure 7B and C) and the results suggested that a higher CTLA4 expression was associated with the low-risk group, while there was not a statistical difference in PDCD1.

Figure 7 Immune landscape of the signature and the prediction of drug response: (A) the correlation of tumor-infiltrating immune cells and risk scores based on 7 known algorithms. (B and C) CTLA4 and PDCD1 expressed levels in high- and low-risk groups. (D) The abilities of the risk model to predict drug sensitivity: gefitinib, erlotinib, docetaxel, cisplatin, vinorelbine, and paclitaxel. *** represented p < 0.001 and ns represented p ≥ 0.05.

Constructed Model Effectiveness for Drug Response Prediction

We evaluated the difference in IC₅₀ between high- and low-risk groups for 7 drugs: 2 EGFR-TKIs and 5 chemotherapeutic agents (Figure 7D). The results showed that the high-risk group is significant correlated with lower IC₅₀.

Discussion

LncRNA possesses many biological functions and has gained increasing attention for its value as a prognostic factor in malignant tumors. Although several prognostic irlncRNA signatures of lung adenocarcinoma based on the TCGA-LUAD dataset have been developed, they are all dependent on specific lncRNA expression levels of transcripts.^9–14 Therefore, this is the first study to apply a two-gene combination strategy to establish lung adenocarcinoma signatures.

First, we obtained DEirlncRNAs by a differential co-expression strategy using an immune gene list from the ImmPort database, while previous studies retrieved them from the Molecular Signature Database or the Gene Set Enrichment Analysis databases.^9,10,12,14 Additionally, Wang et al obtained DEirlncRNAs by differential expression analysis between high- and low-ESTIMATE score groups.¹³ Second, the DEirlncRNAs pairs matrix was arranged based on a cyclic comparison of each lncRNA pair expression level. We constructed the matrix to minimize sample errors derived from expression variations. Third, univariate and LASSO penalized regression analysis was performed to establish the risk signature and a ROC curve was plotted, since its commonly used in bioinformatics research. Fourth, the novel signature prognostic value was evaluated in different situations: various patients’ clinical characteristics, tumor-infiltrating cells, drug response, and ICI-related genes.

A total of 11 irlncRNA pairs were included in our signature, of which LINC02036|AL365181.2, AC026355.2|LIN01547, AC104958.2|AC125807.2, and AC007991.2|AL365181.2 were regarded as risk protective factors. They were identified for the first time, except AC007991.2, which has been verified as a protective autophagy-related lncRNA (HR: 0.609 p = 0.004) in bladder cancer.²⁸ The remaining 7 pairs were considered as risk-related factors. LINC01116, DRAIC, AC022784-1, and LINC01843 were explored in former irlncRNA LUAD signatures.^11–13 Previous studies have implied that overexpressed LINC01116 may promote tumor proliferation, metastasis, and drug resistance in LUAD.^29–31 DRAIC was recognized to have a role in migration and invasion regulation, as well as a predictor in various cancers.³² The ROC curve indicated that the signature was stable and reliable - AUC superior to most existing signatures. We demonstrated that our signature was significantly related to patients’ LUAD prognosis, and proved to be an independent prognostic factor by univariate and multivariate Cox regression analysis. These findings suggested that our signature could potentially help clinicians in the LUAD patients’ treatment.

Tumor-infiltrating immune cells have been identified as a pivotal component of the tumor immune microenvironment and are involved in the response to immunotherapy. Immune scores are commonly applied to guide immunotherapy and evaluate prognosis. As mentioned in the previous reference, lung cancer patients in the high-immune score group had significantly longer disease-free and overall survival than those in the low-immune score group.³³ Due to the complexity and defects of each immune infiltrate algorithm, we performed a correlation analysis between risk scores and immune scores using 7 algorithms. Prior studies have investigated the correlation between TIMER estimated infiltration value and risk scores of a 4-irlncRNA signature including HSPC078, DRAIC, AP004608.1, and MIR223HG in LUAD patients. The results showed that CD4⁺ T cell, CD8+ T cells, neutrophils, macrophages, and myeloid dendritic cells were significantly negatively associated with risk scores.¹³ In the immune landscape that we constructed, B cells (B cell memory and B cell plasma) were significantly negatively correlated to the risk score. T cells - T cell CD4⁺ central memory, T cell CD4⁺ effector memory, T cell CD4⁺ memory resting, T cell CD8⁺, T cell follicular helper, T cell NK, and T cell regulatory (Tregs) - were significantly negatively correlated to risk score, except for T cell CD4⁺ (non-regulatory), T cell CD4⁺ memory activated, and T cell CD4⁺ Th2. These findings suggested that the 7 algorithm combination that we performed presented a more abundant immune landscape compared to individual algorithms. We assessed CTLA4 and PDCD1 expression levels and the results showed that the signature was significantly related to CTLA4. Although CTLA4 and PDCD1 expression have an irreplaceable guiding role in current immunotherapy, this enriches the immune landscape of our signature. Our model suggested that the low-risk group was associated with sensitivity to chemotherapeutics (docetaxel, cisplatin, vinorelbine, and paclitaxel), and targeting-therapy (gefitinib and erlotinib). This indicated that our signature may help clinicians with the LUAD’s treatment.

Finally, we constructed a novel irlncRNA LUAD signature using the TCGA database, based on paired irlncRNAs, and that does not rely on lncRNAs’ specific expression levels. Further analysis confirmed that the prognostic model could be used as an independent prognostic factor for LUAD patients and presented an abundant immune landscape.