Identification of P2RY13 as an immune-related prognostic biomarker in lung adenocarcinoma: A public database-based retrospective study

Data source and processing

The fragments per kilobase of transcript per million mapped reads level of gene-expression matrixes of patients with LUAD were obtained from TCGA. The raw data of the mRNA expression matrix of GSE68465 were collected from the GEO database and normalized with the ‘affy’ package in R 3.6.3 (https://www.r-project.org). The clinicopathological parameters of each patient were also downloaded. Cases lacking pathological diagnosis were excluded.

ESTIMATE algorithm-derived immune and stromal scores

As described previously (Chen et al., 2019b; Chen et al., 2019c), the immune and stromal scores of each sample were calculated with the package ‘ESTIMATE’ in R. Next, the scores of patients in TCGA dataset were evaluated with the package ‘survminer’ to infer the optimal cut-off value, which divided patients into high- or low-immune/stromal score groups. Kaplan–Meier plot and log-rank test were performed to construct a survival curve to illustrate the association of immune/stromal scores and overall survival of patients with LUAD. Patients in the GEO dataset were classified into high- or low-immune/stromal score groups according to the optimal cut-off value, and survival analysis was performed.

Expression analysis of DEGs

Using | log fold-change (FC) | >1.2 and false discovery rate (FDR) <0.05 as the criteria, the Bioconductor package ‘edgeR’ was used to determine DEGs between high- and low-immune/stromal score groups in TCGA dataset. The overlapping DEGs were utilized for further analysis.

Functional enrichment analysis and protein-protein interaction (PPI) network construction

All the overlapping DEGs were used for Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) analyses. FDR < 0.05 was set as the threshold. In addition, a PPI network of all the overlapping DEGs was obtained from Search Tool for the Retrieval of Interacting Genes/Proteins (https://string-db.org) with a confidence >0.9 as the threshold, and was reconstructed with Cytoscape version 3.6 (https://cytoscape.org).

Identification of prognostic immune-related biomarkers

All the overlapping DEGs in TCGA dataset were subjected to univariate Cox regression analysis to identify prognosis-related genes. Survival analysis was applied to compare the survival difference between high- and low-expression of DEGs. Similarly, univariate Cox regression and survival analyses of DEGs were also performed in the GEO dataset to validate the result in TCGA dataset. Genes meeting the criteria (univariate Cox regression analysis, P < 0.01; survival analysis, P < 0.05) in both datasets were selected as prognostic biomarkers for further research.

Group comparisons of the expression levels of biomarkers among different clinical characteristics were performed with Student’s t-test. Univariate and multivariate Cox regression models were conducted to identify whether the expression of biomarkers was an independent prognostic factor for LUAD.

Single-sample gene set enrichment analysis (ssGSEA)

ssGSEA was performed to assess the level of tumor-infiltrating immune cells (TIICs) and immune functions with the BiocManager package ‘Gene Set Variation Analysis’. The method for ssGSEA was based on a rank value of each gene, which defined a score representing the degree of absolute enrichment of a particular gene set in each sample. In total, 28 specific gene sets (15 immune cell gene sets and 13 immune function gene sets) were acquired from other studies (Charoentong et al., 2017; Shi et al., 2020; Zuo et al., 2020). The association of biomarkers with TIICs and immune functions were investigated with Spearman correlation analysis.

Association of biomarkers with immunomodulators and patients’ response to immunotherapy

In the present study, several key immunomodulators [cytotoxic T-lymphocyte-associated protein 4 (CTLA-4), intercellular adhesion molecule 1, inducible T-cell costimulatory, interferon- γ, lymphocyte activating gene 3, T cell immunoreceptor with Ig and ITIM domains, natural killer gene 2A, programmed cell death protein 1 (PD-1), programmed death-ligand 1 (PD-L1), T-cell immunoglobulin and mucin-domain containing-3 and V-domain Ig suppressor of T cell activation] were quantified. Student’s t-test and Spearman correlation analysis were performed to determine the association of immunomodulators and immune/stromal scores as well as prognostic biomarkers.

The Cancer Immunome Atlas (https://tcia.at/) is a public database, which analyzes next-generation sequencing data to present immune landscapes and anti-genomes of 20 solid tumors, and calculates the immunophenoscore (IPS) (Charoentong et al., 2017). The IPS value is ranked from 0 to 10, and is positively correlated with tumor immunogenicity. Furthermore, the IPS could reflect the response to immune checkpoint inhibitors treatment. The present study analyzed two types of IPS values (IPS: PD-1/PD-L1/PD-L2 blocker and IPS: CTLA-4 blocker) to investigate the different responses to anti-PD-1/PD-L1 and anti-CTLA-4 treatment between patients with low- and high-stromal score/immune score/P2Y purinoceptor 13 (P2RY13) level with Wilcoxon signed-rank test.

GSEA

To determine the potential signaling pathways altered by prognostic biomarkers, GSEA was performed with FDR<0.05 as the threshold.

Patient cohorts

A total of 896 LUAD samples (465 samples in TCGA dataset and 431 samples in the GEO dataset) were identified. The detailed demographic and baseline characteristics of these 896 patients with LUAD are presented in Table 1. The flow diagram of the present study is shown in Fig. 1.

Evaluation of immune and stromal scores

Using the ESTIMATE algorithm, the present study calculated the immune and stromal scores of patients in both the GEO and TCGA datasets. The optimal cut-off value for immune score and stromal score was 1,795.5 and 13.5, respectively, which divided patients into high- and low-immune/stromal score groups. Survival analysis in TCGA dataset demonstrated that the prognosis of patients with LUAD with high-stromal/immune scores was better than that of patients with low-stromal/immune scores (Figs. 2A and 2C), which was similar to the results obtained in the GEO dataset (Figs. 2B and 2D).

In addition, the immune scores were significantly different among different tumor-node-metastasis (TNM) stages and tumor sizes in both TCGA dataset (TNM stage, P = 0.039; tumor sizes, P < 0.001; Fig. S1) and the GEO dataset (TNM stage, P = 0.046; tumor sizes, P = 0.045; Fig. S2). The stromal scores were significantly different among tumor sizes in the GEO dataset (P = 0.014; Fig. S2), whereas in TCGA dataset there were no differences (P = 0.210; Fig. S1). There was no statistically significant difference in immune or stromal scores between patients who were < 65 and ≥65 years of age; between female and male patients; between patients with and without lymph node metastasis; between patients with and without distant metastasis; or between patients with and without epidermal growth factor receptor (EGFR)/KRAS proto-oncogene, GTPase (KRAS) mutation (data not shown). Furthermore, no difference in stromal scores was observed among different TNM stages.

Functional enrichment analysis and PPI network

A total of 195 downregulated and 108 upregulated overlapping DEGs were identified between the high- and low-immune/stromal score groups in TCGA dataset (Fig. 3A). Next, KEGG and GO analyses were performed to illustrate the role of 303 overlapping DEGs. KEGG analysis demonstrated that 17 pathways were enriched in these DEGs, including ‘Cytokine-cytokine receptor interaction’, ‘Chemokine signaling pathway’ and ‘B cell receptor signaling pathway’, which were highly associated with the immune system (Fig. 3B). In GO analysis, 34 terms (four terms of molecular function, five terms of cellular component and 25 terms of biological processes) were identified (Table S1).

For exploring the interactions among 303 overlapping DEGs, a PPI network was constructed, which consisted of 98 nodes (90 upregulated and 8 downregulated DEGs) and 395 edges (Fig. 3C).

P2RY13 is an immune-related prognostic biomarker

Univariate Cox regression and survival analyses were applied to determine the best prognosis-related genes. The results of the univariate Cox regression analysis showed that, among 303 DEGs, 24 genes in TCGA dataset and 20 genes in the GEO dataset exhibited P < 0.001 (Fig. 4A and 4B). In addition, survival analysis demonstrated that 30 genes in TCGA dataset and 21 genes in the GEO dataset had an important effect on the prognosis of LUAD (data not shown). Subsequently, immune-related prognostic biomarkers were identified with the criteria univariate Cox regression analysis, P < 0.01 and survival analysis, P < 0.05 in both sets, and it was found that only one gene (P2RY13) met the aforementioned criteria: Univariate Cox regression analysis: TCGA dataset, hazard ratio (HR) = 0.736, 95% confidence interval (CI) = 0.601–0.900, P = 0.003 (Fig. 4A) and GEO dataset, HR = 0.474, 95% CI [0.278–0.811], P = 0.006 (Fig. 4B); and survival analysis: TCGA dataset, P = 0.006 (Fig. 4C) and GEO dataset, P < 0.001 (Fig. 4D). Therefore, P2RY13 was selected as a prognostic biomarker for further investigation.

Correlation between P2RY13 and clinical characteristics

As shown in Figs. 5A and 5B, the expression of P2RY13 in patients with high stromal/immune score was significantly upregulated in TCGA dataset (stromal score, P < 0.001; immune score, P < 0.001), which was in agreement with the results obtained in the GEO dataset (stromal score, P < 0.001; immune score, P < 0.001). In both datasets, patients with stage I/II had higher P2RY13 expression level (TCGA dataset, P = 0.045; GEO dataset, P = 0.039). Similar results were observed in patients with T1/T2. In addition, the P2RY13 expression level of patients with lymph node metastasis in TCGA dataset was significantly decreased (P = 0.037). However, in the GEO dataset, there was no difference between patients with or without lymph node metastasis (P = 0.293). In TCGA dataset, no difference was observed between patients with or without EGFR mutation (Fig. S3A). A similar phenomenon was found in patients with or without KRAS mutation (Fig. S3A).

Univariable Cox regression analysis in two datasets revealed that the expression level of P2RY13 was a meaningful factor influencing the prognosis of patients with LUAD (TCGA dataset, HR = 0.615, 95% CI [0.433–0.873], P = 0.006; GEO dataset, HR=0.631, 95% CI=0.487–0.819, P < 0.001) (Table 2). In addition, multivariable Cox regression analysis indicated that P2RY13 was an independent prognosis-related factor (TCGA dataset, HR = 0.601, 95% CI [0.423–0.856], P = 0.005; GEO dataset, HR = 0.760, 95% CI [0.582–0.994], P = 0.045) (Table 2).

Next, the expression of P2RY13 at the protein level was investigated with the online database The Human Protein Atlas (https://www.proteinatlas.org). The results of immunohistochemistry showed that the representative protein expression of P2RY13 in LUAD tissues was downregulated (Figs. S3C and S3D). The prognostic value of P2RY13 at the protein level was explored with Clinical Proteomic Tumor Analysis Consortium (National Cancer Institute; https://proteomics.cancer.gov/programs/cptac), and the results revealed that a high P2RY13 protein level predicted improved prognosis (P = 0.021; Fig. S3B).

P2RY13 is associated with TIICs and immune functions

ssGSEA was performed to evaluate the level of TIICs and immune functions. The distribution of 28 gene sets, including 15 immune cell gene sets and 13 immune function gene sets, is presented in Fig. 6A. Spearman correlation analysis revealed a strong positive correlation between the expression level of P2RY13 and 28 gene sets (Fig. 6B). In addition, the association between the immune and stromal scores and these 28 gene sets was statistically significant (Fig. 6B).

Association of P2RY13 with immunomodulators and patients’ response to immunotherapy

The present study quantified 11 immunomodulators, all of which were significantly upregulated in the high-stromal/immune score group (Fig. S4A and S4B). Furthermore, P2RY13 was positively correlated with all the 11 immunomodulators (Fig. 6C).

Next, the responses to anti-PD-1/PD-L1 and anti-CTLA-4 treatment among different groups were explored. Both the IPS: PD-1/PD-L1/ PD-L2 blocker and the IPS: CTLA-4 blocker were higher in patients with LUAD with high P2RY13 expression level (Figs. 6D and 6E), indicating that the relative probabilities to respond to anti-PD-1/PD-L1 and anti-CTLA-4 treatment were higher in patients with high P2RY13 expression level. Similar results were observed in patients with high-stromal/immune score (Figs. S4C and S4D).

GSEA

For exploring the changes in KEGG pathways between the high- and low-expression level of P2RY13 groups, GSEA analysis was performed. The result indicated that 26 pathways were enriched (Fig. 7A). Of note, in patients with high P2RY13 expression, a few immune-related pathways were identified, including ‘B cell receptor signaling pathway’, ‘T cell receptor signaling pathway’, ‘Intestinal immune network for IgA production’, ‘Natural killer cell-mediated cytotoxicity’, ‘Primary immunodeficiency’, ‘Chemokine signaling pathway’ and ‘Cytokine-Cytokine receptor interaction’ (Fig. 7B).