Using ESTIMATE algorithm to establish an 8-mRNA signature prognosis prediction system and identify immunocyte infiltration-related genes in Pancreatic adenocarcinoma

The ESTIMATE algorithm identifies DEGs associated with stroma and immune scores

We downloaded the transcriptome and clinical data of 182 PAAD program cases from the TCGA database. Based on the ESTIMATE database results and the median values of the stromal and immune scores (excluding the normal pancreas cases), we divided these patients into high and low stromal score groups and high and low immune score groups. We displayed distinct these clinicopathological features, including age, sex, N stage, M stage, histologic grade, T stage, and primary tumor sites, and mRNA expression forms between high and low stromal (Figure 1A) and between high and low immune (Figure 1B) score groups in the form of heatmaps. We found statistically significant results only for sex and M stage between high and low stromal scores (P = 0.010 and P = 0.010, respectively), whereas only for N stage between high and low immune scores (P = 0.049) as shown in Table 1. And the age, M stage, histologic grade, T stage, and primary tumor sites were not significantly different between two groups. Low stromal scores group was characterized by a higher number of male patients and poorly M stages compared with high stromal scores group. Low immune scores group included more patients with better N stages compared with high immune scores group. Based on our threshold value (log fold change ≥1.5 and P <0.05), there were 333 DEGs between high and low stromal groups and 314 DEGs between high and low immune score groups (Supplementary Tables 5 and 6). The Venn diagram (Figure 1C) used to find the co-DEGs in the two groups showed that 203 genes, which accounted for the majority of DEGs, were differentially expressed simultaneously in two score-related comparisons. Moreover, these simultaneously upregulated and downregulated differentially expressed DEGs were 196 and 7, respectively, with the same alternative trend between the two groups (Figure 1D and 1E). As shown in Figure 1F and 1G, CIBERSORT results revealed that CD4+ activated memory T cells (P = 0.005), M2 macrophages (P = 0.012), and neutrophils (P<0.001) were upregulated, while activated NK cells (P = 0.011) and M0 macrophages (P = 0.001) were downregulated in the high stromal-score group. In the high immune-score group, CD8 T cells (P = 0.001), gamma delta T cell (P = 0.041), resting NK cells (P=0.026), monocytes (P = 0.002), and neutrophils (P = 0.003) were upregulated, while plasma cells (P = 0.043) and M0 macrophages (P = 0.001) were downregulated. Hence, we used the ESTIMATE algorithm to identify DEGs associated with stroma and immune scores successfully.

KEGG and GO functional enrichment analysis for DEGs

To identify differentially regulated genes on immune cell infiltration, we used Metascape tools to conduct the GO and KEGG analyses of the 203 simultaneously expressed DEGs. The Heatmap in Figure 2A exhibits the top 15 enriched GO terms across the DEGs based on Biological Process (BP), Cellular Component (CC), and Molecular Function (MF), and the heatmap in Figure 2B also exhibits the top 15 enriched KEGG pathways across the DEGs. Among them, most terms or pathways were associated with biological processes and lymphocyte activation. We also applied protein-protein interaction networks from the DEGs, which were coloured by different cluster IDs (Figure 2C) and P-values (Figure 2D), and the networks contained more genes that tended to have more significant p-values. Since most networks from the DEGs were related to biological processes and lymphocyte activation, the result confirmed in reverse that the DEGs identified based on stroma and immune scores were accurate. Consequently, combined with the result showed in Figure 1, our study identified that most of the DEGs identified based on the immune and stromal scores were correlated considerably with the immune responses or immune cell infiltration in PAAD.

An 8-mRNA signature system was established to predict the overall survival of PAAD patients

Based on the 203 DEGs obtained from the mRNA difference analysis (Figure 1C), we first conducted Cox univariate regression for initial screening to remove interferences from excessive confounding genes and acquire genes with the most significant impact on prognosis. As shown in Supplementary Table 1, to avoid exclusion of important prognostic genes, 14 of the mRNA with P <0.1 were then moved into LASSO regression [16, 17]. We presented the LASSO coefficient profiles of the 14 mRNA (Figure 3A) and produced 10-fold cross-validation results that identified optimal values of the penalty parameter λ (Figure 3B).

The forest plot demonstrated the association between each gene and overall survival (Figure 3C). According to these results, we identified an 8-mRNA signature to evaluate the overall survival time of PAAD patients based on the expression of the 8 mRNAs and their regression coefficients as follows: Risk score = (-0.01579 × expression level of ADH1B) + (0.07838 × expression level of CA9) + (-0.76361 × expression level of CDHR3) + (0.33812 × expression level of CXCL9) + (-0.06606 × expression level of GIMAP7) + (-0.70384 × expression level of ICAM3) + (-0.11457 × expression level of LDLRAD1) + (-0.24179 × expression level of P2RY8). Patients in the TCGA cohort were divided into a low-risk group (N=88) and a high-risk group (N=88) utilizing the median risk score as the cut-off value.

The Kaplan-Meier curves showed that high-risk patients had significantly worse survival rates in the training set (P<0.001) (Figure 3D and Supplementary Figure 2A). Moreover, the multivariate analysis revealed that this prognostic signature was an independent prognostic factor in PAAD patients (Table 2).

We conducted a time-dependent ROC analysis at one, two, and three years to assess the prognostic accuracy of the risk score (Figure 3F and Supplementary Figure 2C). We used 69 patients from ICGC PACA-Australia as the validation cohort and calculated the risk score of each case based on our 8-mRNA signature. Per the Kaplan-Meier curves, the high-risk group patients had worse outcomes (Figure 3E and Supplementary Figure 2B); the time-dependent ROC analysis verified that the risk score had a good long-term prognostic accuracy (Figure 3G and Supplementary Figure 2D). We, therefore, successfully established an 8-mRNA signature to assess the prognosis of PAAD patients.

CA9, CXCL9, and GIMAP7 have a close relationship with the immune infiltration of the tumour microenvironment in PAAD

To further explore the regulatory mechanisms for the prognosis of PAAD based on the eight differentially expressed genes in our signature, we searched the GEPIA web tool [18] for an association between the expression of genes at the mRNA level and Overall Survival (OS). The result showed that carbonic anhydrase 9 (CA9), C-X-C motif chemokine ligand 9 (CXCL9), and GTPase, IMAP family member 7 (GIMAP7) were independently associated with the overall survival (OS) rate of PAAD (Figure 4A–4C) but not with the other five factors (Supplementary Figure 1A–1E). So CA9, CXCL9, and GIMAP7 were included for further exploration.

As shown in Figure 1F and 1G, the DEGs were identified based on the immune microenvironment scores, and the GO and KEGG analysis showed that the DEGs correlated significantly with immune reactions or immune cell infiltration in PAAD. So we tried to determine the relationship between the expression levels of these three factors and the immune infiltration level in the tumour microenvironment. The infiltration levels of CD8+ T cells, macrophage cells, neutrophil cells, and dendritic cells were down-regulated by CA9 (Figure 4D). The evaluation of CXCL9 in PAAD showed a positive correlation with infiltration levels of B cells, CD8+ T cells, CD4+ T cells, macrophage cells, neutrophil cells, and dendritic cells (Figure 4E). The analysis of GIAMP7 revealed that the infiltration levels of B cells, CD8+ T cells, CD4+ T cells, macrophage cells, neutrophil cells, and dendritic cells were all up-regulated with the increased expression level of GIMAP7 (Figure 4F). When taken together, these data suggest that CA9, CXCL9, and GIMAP7 contribute to the immune infiltration of the tumour microenvironment in PAAD

CA9, CXCL9, and GIMAP7 regulate the expression level of FOXO1 in PAAD

Reportedly, the Forkhead box O (FOXOs) family intrinsically influences the anti-tumour immune response and the infiltration levels of immune cells, including CD4+ T cells and CD8+ T cells, B cells, neutrophil cells, macrophage cells, and dendritic cells [19]. Consistent with this report, our data also showed a significant correlation between the expression level of FOXO1 or FOXO3, but not FOXO4, and the infiltration levels of these immune cells (Figure 5A–5C). Furthermore, by searching the GEPIA web tool, we found a negative correlation between the mRNA expression level of CA9 and the FOXOs family in PAAD patients (Figure 5D). On the contrary, the mRNA expression level of CXCL9 correlated positively with FOXOs (Figure 5E), and so did mRNA expression level of GIMAP7 with FOXOs in PAAD (Figure 5F). Consequently, we hypothesised that CA9, CXCL9, and GIMAP7 regulated the infiltration level of immune cells by modifying FOXOs in PAAD.

To determine whether CA9, CXCL9, and GIMAP7 regulated the immune cell infiltration via FOXOs, we first analyzed the relationship between FOXO3 and FOXO4 and these three genes. We knocked down each of CA9, CXCL9, or GIMAP7 using a gene-specific siRNA and detected the changes in the mRNA level of FOXO3 or FOXO4. We found that these three genes (CA9, CXCL9, or GIMAP7) had apparent effects on FOXO3 and FOXO4 in both PANC-1 and SW 1190 cells (Figure 6A–6C). In contrast, the examination of the correlation between FOXO1 with these three genes indicated that repressing the level of CA9 increased the protein and mRNA levels of FOXO1 in PAAD cells (Figure 6D), whereas, the inhibition of CXCL9 or GIMAP7 by siRNA down-regulated the levels of FOXO1 in PANC-1 and SW 1990 cells (Figure 6E and 6F).

We next examined the downstream inflammatory factors of FOXO1, such as IL-8, IL-10, VEGF, and PD-L1 [20–23] and found that knocking down CA9 significantly decreased the mRNA levels of IL-8, IL-10, VEGF, and PD-L1, but the knockdown of CXCL9 and GIMAP7 substantially increased IL-8, IL-10, VEGF, and PD-L1 expression levels in PAAD cells (Figure 6H and Supplementary Figure 1F). Therefore, our results indicate that CA9 transcriptionally inhibits, while CXCL9 and GIMAP7 transcriptionally promote FOXO1 expression in PAAD.

CA9, CXCL9, and GIMAP7 correlate with FOXO1 in PAAD patient specimens

To determine the correlation between these three proteins and FOXO1 in PAAD specimens, we examined the expression of all four proteins by immunohistochemistry (IHC) on a TMA containing a PAAD samples cohort (n = 31). Representative images of the high and low/no stainings of CA9, CXCL9, GIMAP7, and FOXO1 are shown in Figure 7A, 7D, and 7G. The IHC score was calculated by multiplying the percentage of positively stained cells with the staining intensity. The expression of CA9 correlated negatively with FOXO1 (Spearman’s product-moment correlation coefficient r = -0.5392, P = 0.0017) (Figure 7B, 7C), which was consistent with the protein and mRNA level changes in pancreatic cell lines reported above. Moreover, our data showed that there was a positive relationship between the expression levels of CXCL9 or GINAP7 and FOXO1 in PAAD patient specimens (Spearman’s product-moment correlation coefficient r = 0.7416, P < 0.001 for CXCL7 and FOXO1, Spearman’s product-moment correlation coefficient r = 0.4853, P = 0.0053 for GIMAP7 and FOXO1) (Figure 7D–7G). Consequently, our results suggest that the expression of CA9, CXCL9, and GIMAP7 all correlated with the expression of FOXO1 in PAAD patient specimens.

Downloading transcriptome datasets and the stromal and immune scores of PAAD

TCGA PAAD transcriptome FPKM data were downloaded from the GDC Data Portal (https://portal.gdc.cancer.gov/). Log2-transformed normalized ICGC PACA-Australia data were obtained from UCSC Xena (https://xenabrowser.net/). The included criteria for the analysis were: (a) pathological type was pancreatic carcinoma; (b) overall survival (OS) data were available; (c) raw count or normalized gene expression data were available. Per these criteria, we enrolled 176 and 69 pancreatic carcinoma cases for TCGA and UCSC, respectively [29]. Clinical data, such as the age, gender, race, survival time, and status of the selected patients, were also obtained from the two databases. We got the stromal and immune scores of TCGA PAAD patients from the ESTIMATE database (https://bioinformatics.mdanderson.org/estimate/index.html). The TCGA PAAD cohort was set as the training group for prognostic mRNA signature exploration, and the ICGC PACA group was used for validation. The mRNA names of probes were mapped according to their ensemble id, and the median value was used when several inquiries represented the same mRNA.

Identifying differentially expressed mRNAs based on stromal and immune scores

We separated the TCGA cohort based on the median values of the stromal and immune scores of samples and obtained two groups of DEGs through the linear models for microarray data (LIMMA) R package. Fold change ≥2 and P <0.05 were set as the cut-off values for DEGs screening. Heatmaps were drawn using the heatmap R package. We also performed gene ontology (GO) enrichment analysis, KEGG pathway analysis, and protein-protein interaction enrichment analysis for DEGs using the Metascape (http://metascape.org) website tool [45].

Estimating immune cell type proportions

The CIBERSORT R package with an LM22 signature and 1000 permutations was used to calculate the percentages of immune cells in PAAD cases. Per the description of a previous study [46], we inferred the proportions of immune cell types in the transcriptome data of cases and differences in immune cell type composition between different analysis groups (high immune score group vs low immune score group; high stromal score group vs low stromal score group) using the Microenvironment Cell Populations-counter method. To directly visualize the estimation result, we applied the “barplot” and “vioplot” R packages.

Identifying prognosis-related genes

We selected overall survival (OS) as the endpoint and made verifications using the model for disease-free survival analysis. First, we used the univariate COX regression for the initial mRNA screening. The least absolute shrinkage and selection operator (LASSO) Cox regression model was then employed for the further selection of prognostic mRNAs [47]. We used the “Survival” and “glmnet” R packages for analyses and obtained the final model.

Based on the LASSO Cox regression result, we got a group of mRNAs and built an mRNA-based prognostic signature for the training cohort. Each patient’s risk score was calculated using a combination of the expression levels of mRNAs and LASSO-Cox regression coefficients. We used the median risk score as the cut-off value to divide the training set patients into a high-risk group and low-risk group, and the Kaplan-Meier curves and log-rank analysis were used to identify survival differences. The time-dependent receiver operating characteristic (ROC) curves of the risk scores were generated using the time ROC R package to explore the prognostic accuracy of risk scores. Univariate and multivariable Cox analyses were used to study the independent prognostic values of risk scores compared with other clinical characteristics. Lastly, the ICGC-AU data were used as the validation cohort to verify the effect of the risk-score model.

Screening DEGs for further exploration

We further screened the genes in the LASSO Cox regression model for their biological regulatory mechanisms in PAAD. We used the webserver GEPIA (http://gepia.cancer-pku.cn/) to identify genes that could independently predict overall survival in patients with PAAD [18]. Considering the correlation between these genes and immune components, we also analyzed the relevance of interferon-gamma (as an essential immunoregulatory factor) to these genes. The web server TIMER (https://cistrome.shinyapps.io/timer/) was used to investigate the correlation between these genes and tumour-infiltrating immune cells [48].

Immunohistochemistry (IHC)

To study altered protein expression, an IHC analysis was performed on tissue microarray (TMA) slides purchased from Outdo Biobank (Cat No. XT14-029) (Shanghai, China), using GIMAP7 antibody (Sigma, Cat# HPA020268,1:50), anti-CXCL9 antibody (Abcam, Cat# ab9720, 1:400), anti-CA9 antibody (Proteintech, Cat# 11071-1-AP, 1:200), and anti-FOXO1 antibody (Proteintech, Cat# 18592-1-AP, 1:1000). Semi-quantitative scoring was performed based upon the staining intensity (negative = 0; weak = 1; moderate = 2; and intense = 3). Two senior pathologists rated the degree of immunostaining of formalin-fixed, paraffin-embedded sections independently in a blinded manner. The scores were determined by the percentage of positive cells multiplied by the staining intensity. The clinical information on Tissue microarray (TMA) slides is shown in Supplementary Table 4.

Statistical analysis

The R package was used to perform the differential expression genes (DEGs) analysis to obtain the genomic profile between high and low immune/stromal groups. The LASSO Cox regression model was employed for the further selection of prognostic mRNAs. The Chi-square test was used to calculate differences in clinicopathological variables between groups. The Kaplan-Meier method was used to determine disease-specific survival (DSS) and OS, and the differences between the study groups were compared using the log-rank test. Univariate and multivariate Cox proportional hazard regression models were used to calculate hazard ratios for OS and DSS. All the P values were adjusted using the FDR. A p-value less than 0.05 was considered significant. The described statistical analysis was performed with the IBM SPSS Statistics software version 24 for Windows (IBM Corporation, Armonk, NY, USA).

For RT-PCR and IHC analysis, statistical analysis was carried out using GraphPad Prism 6.0 (GraphPad Software Inc., San Diego, CA). Statistical significance was assessed using the unpaired two-tailed Student t-test between two groups or the one-way ANOVA with Tukey post hoc test for multiple comparisons. P values less than 0.05 were considered significant. All the values are expressed as the mean ± SD. Asterisks used to indicate significance correspond with *, P < 0.05; **, P < 0.01; ***, P <0.001.