Identification of immune-related long noncoding RNAs (irlncRNAs)
The process flow of this study is shown in Figure 1. First, we obtained the transcriptome profile from high-throughput sequencing (HTseq) count data (npatient=176, nnormal=4) and the clinical information of patients from The Cancer Genome Atlas pancreatic cancer dataset (TCGA-PAAD). Second, we retrieved the gene transfer files (GTFs) from Ensembl to annotate the lncRNAs from the expression matrix. Third, immune-related genes (irgenes) were downloaded from the ImmPort database. Pearson correlation analysis between the lncRNAs and irgenes was completed, and 1903 irlncRNAs were identified (R>0.5, p<0.001) (Table S1).
Establishment of risk model by random survival forest (RSF) analysis
RSF was applied constructed to determine the irlncRNAs that are of most significance to the OS of pancreatic patients. The samples were randomly divided into a training set (n=123) and a test set (n=53) at a ratio of 7:3. In the survival tree analysis, the 9 most significant variables (CTC-529P8.1, RP11-706C16.8, LINC01493, LINC01887, LINC00462, LINC01510, LINC02205, RP11-1082L8.2 and RP11-402N8.1) were selected according to the variable importance (VIMP) and the minimal depth. The subsequent Cox regression analysis identified 3 irlncRNAs (LINC00462, LINC01887, RP11-706C16.8), with a coefficient index and risk score for each sample in the training set calculated. To validate this model, the receiver operating characteristic (ROC) curves were drawn, and the areas under the curves (AUCs) for 36 months, 30 months, 24 months, 18 months, 12 months and 6 months were 0.778, 0.774, 0.751, 0.753, 0.780 and 0.756, respectively (Figure 2). In addition, the C-index for this risk model was 0.696 (p<0.001).
Clinical evaluation of the risk model
Instead of the median risk score value, the best cutoff value calculated by X-Tile software (https://medicine.yale.edu/lab/rimm/research/software/) was utilized to divide the samples into high-risk and low-risk groups in the training set and validation set. It was 1.44 in training set and 1.40 in test set. To clinically evaluate the risk model, several analytical methods were applied. First, Kaplan-Meier analysis showed that the OS of patients in the high-risk group was significantly lower than that of patients in the low-risk group (p<0.001). In the high-risk group, the 5-year OS rate was 5.66%, and the 95% confidence interval (CI) was 0.88 to 361; the 5-year OS rate was 30.9%, and the 95% CI was 18.28 to 52.4 in the low-risk group. In the test set, Kaplan-Meier analysis showed a significant difference in OS between the two groups (p<0.001) (Figure 3).
The risk model is an independent prognostic factor for pancreatic cancer
To construct a more accurate clinical prognostic risk model, the AUC of each ROC for every clinical characteristic and the risk score was calculated; the AUCs of the risk score, age, sex, T stage, N stage, M stage and stage were 0.778, 0.654, 0.548, 0.564, 0.686, 0.462 and 0.557, respectively, with the risk score AUC being the only one above 0.7. After univariate Cox regression analysis, the risk score (p<0.001), N (p=0.006) and T (p=0.034) were selected for multivariate Cox regression analysis, which revealed that the risk score (p<0.001) was the only independent prognostic factor for pancreatic cancer (Table 1). Moreover, the addition of the risk score to the clinical model can raise the C-index from 0.599 to 0.682, indicating that it greatly contributes to prognosis prediction. A nomogram and the related calibration plots were established to visualize the specific method, calculate the risk scores and show the ability of the model to predict OS at 6 months, 12 months and 36 months (Figure 4).
Exploring the correlation between immune cell infiltration and the risk model
The irlncRNAs identified by correlation analysis of irgenes and lncRNAs at the beginning of this study may exert an influence on immune status, including immune cell infiltration. We employed the Tumor IMmune Estimation Resource (TIMER), CIBERSORT, QUANTISEQ, MCPcounter and EPIC resources to delineate the immune cell infiltration of the training set (Table S2). Wilcoxon signed-rank analysis between the risk score and immune cell infiltration revealed that high risk was associated with greater infiltration of cancer-associated fibroblasts, follicular helper T cells, CD4+ T cells, M0 macrophages and M1 macrophages, while low risk was correlated with greater infiltration of B cells and M2 macrophages. Pearson correlation analysis showed the correlation between the risk score and the infiltration of 3 types of immune cells (M0 macrophages, CD4+ T cells and M1 macrophages) (p<0.1) (Table S3). The correlation results were expressed in a lollipop graph. However, the Wilcoxon signed-rank test comparing risk and immune checkpoint inhibitor (ICI) biomarkers, including CTLA4, LAG3, IDO1, PDCD1 and ICOS, showed no significant association, which is reflective of the poor effect of ICIs in clinical trials (Figure 5).