Predicting lung adenocarcinoma prognosis with a novel risk scoring based on platelet-related gene expression

Identifying differentially expressed genes (TCGA)

The EdgeR algorithm can identify 188 differentially expressed genes, 83 up-regulated genes and 105 down-regulated genes (Figure 1A). As expected, gene function enrichment analysis showed that the platelet pathway was the most common. Platelet activation is one of the most common biological terms in biological processes (Figure 1B). And platelet activation is most often enriched by differentially expressed genes. (Figure 1C).

Identifying differential genes related to survival

After a single factor analysis, we found that 18 differences are related to lung adenocarcinoma patients’ overall survival rates, and these results are statistically significant (Figure 1D). The forest map shows that most of the identified central genes in the lung adenocarcinoma samples are up-regulated. The forest map of hazard ratios shows that most of these genes are risk factors. Nine genes were selected by multifactor regression, but because the validation group did not have the expression value of TREML1, we finally included only the other eight genes into the study. And we constructed a prognostic marker and divided lung adenocarcinoma patients into two groups. The calculation formula is as follows: Platelet-related risk score (PRS) = RAB27B *0.017537331+CCNA2*0.020495528+PCDH7 *0.038831123+HAS2*0.062731768+P2RX1*0.173961877+PFKP*0.006157281+ENTPD2*0.065271557+TIMP1*0.000897402

Platelet-related receptor genes based on multivariate Cox may be an important tool for diagnosing lung adenocarcinoma patients according to potential, and the ROC curve of CCNA2 is 0.970. (Figure 1E).

Survival analysis showed that the high-risk group’s survival rate was lower than that of the low-risk group, and these results were statistically significant. Figure 2A represents a risk map, including signature-based, individual survival status between groups and the expression levels and distributions of the included modeling genes. It shows a clear difference in survival status between the risk groups, with red dots indicating death and green dots indicating survival. Many deaths occurred in the high-risk group, while most patients in the low-risk group survived through follow-up. In addition, PRS scores are verified in three other independent queues (GSE13213, GSE68465 and GSE72094), and the results are consistent with those of TCGA (Figure 2B–2D). Multivariate Cox regression analysis showed that after adjusting for other parameters like stage and T, the prognostic index can become an independent predictor (Figure 2E). Moreover, multiple Cox regression analysis showed that in GSE13213, GSE68465 and GSE72094 data sets, the prognosis index could also become an independent prediction index (Figure 2F–2H).

The ROC curve was 0.732, indicating that there may be a medium potential for PRS based prognosis in survival monitoring. (Figure 3A) There are still survival differences between the high-risk and low-risk groups in the other three datasets (GSE13213, GSE68465 and GSE72094), and the ROC curve values formed by PRS scores are 0.791, 0.737 and 0.606, respectively (Figure 3B–3D). It was also found that for different lung cancer patient subtypes, the platelet-associated gene risk score was a moderately viable indicator for prognosis. (Supplementary Figure 1).

Using the above eight genes and clinical factors for lung cancer, we developed a nomogram to predict the 1-year, 3-year, and 5-year overall survival. The ROC analysis showed sufficient discrimination with an AUC as 0.746, 0.746 and 0.695, indicating that the nomogram has medium prediction performance (Figure 3E and Figure 3F). Besides, the calibration chart showed the best prediction accuracy, and the predicted survival rate was equal to the actual survival rate (Figure 3G–3I). Moreover, we constructed the nomogram by combining PRS with clinical factors, and provided a free nomogram prediction web tool for clinicians (https://wsxzaq.shinyapps.io/wsxzaq_nomogram/) (Supplementary Figure 2).

Correlation analysis

Correlation analysis between PRS and immune cell infiltration showed that PRS and B_cell or CD4 T_cell (Cor = –0.229, P = 4.952e-07 and Cor = –0.102, P = 0.026) were related, and correlation analysis between PRS and immune checkpoints showed that PRS and PD-L1 (Cor = 0.183, P = 5.575e-05) were related. PRS was negatively correlated with ESTIMATE score (Cor = –0.112, P = 0.014), indicating that PRS was positively correlated with tumor purity of LUAD. PRS was positively correlated with the MKI67 score, which proves that the high-risk group may be closely related to the regulation of cancer cell cycle and cancer cell proliferation (Cor = 0.295, P = 4.865e-11). (Figure 4).

Functional annotation of established eight gene signatures

GSEA was performed in TCGA high-risk populations to study key biological and cellular processes related to prognosis. There were clear abundant GO pathways in high-risk patients, including cell cycle, glycolysis and platelet-derived related factors (Figure 5A). Most of these pathways were involved in cell proliferation and migration and may lead to lung cancer metastasis and recurrence. Previous studies have shown that platelet-derived growth factors regulate glycolysis and therefore affect the cell cycle of cancer cells. This may also be a potential mechanism for high-risk PRS, suggesting poor LUAD prognosis.

Analyzing mutations and CNV in patients with high and low PRS expression

CBioportal database (http://www.cbioportal.org/) shows that among the 8 genes, HAS2 has the highest mutation rate of 10%. (Figure 5B) TTN and MUC16 were detected with high mutation frequencies in tumors expressed in the high-risk group, and they are well-known tumor-associated genes (Figure 5D). This demonstrates that PRS distinguishes mutant genes with poor prognosis. Figure 5C shows the differential copy number in PRS high and low risk, which get copy number in chromosomes 10, 11, 12 and 20, and mitigated in 8, 10, 11 and 12. Among them, the copy numbers of 36 and 20 genes in the high and low-risk groups, respectively, also changed (Figure 5E). The differential expression genes of PRS between high and low expression levels are affected by copy number variation, which indicates that there is no independence between the abnormal copy number expression genes and the differential expression genes of PRS.

The correlation between miRNA and risk score differential expression

MiRNA is a small non-coding RNA, which can be involved in targeted mRNA cleavage and post-transcriptional regulation. Previous studies have shown that miRNA expression imbalance is a key step in lung adenocarcinoma carcinogenesis. Our work evaluated the relationship between miRNA target-controlled down-regulated genes and risk score differentially expressed genes. Therefore, we used TargetScan to predict the relationship between miRNA and its target genes. A total of 120 up-regulated miRNA and 14 down-regulated miRNA were detected in the high-risk group (the low-risk group served as a reference). We identified 233 pairs of down/miRNA-mRNA interactions by TargetScan (i.e. the miRNA and target control gene scores were greater than or equal to 0.6), of which 10 pairs overlapped with risk score differential expression (Figure 5F). Overlapping genes, including XCT and UCK2, are adverse prognostic factors of lung adenocarcinoma, and the mechanism may be related to regulating the cell cycle of cancer cells. These results indicate that miRNA in the down-regulated group has feedback loops in the carcinogenesis processes of the high-risk and low-risk groups.

External data support

By analyzing the HPA (http://www.proteinatlas.org) and TCPA (http://tcpaportal.org) protein databases, we found that the protein expression of related genes constructed by PRS was also consistent with the previous results. (Figure 6A and 6B) Finally, we verified the relationship between PRS and TP53 mutations and KRAS mutations on GSE13213 data sets and found that the risk score was not affected by these two gene mutations. (Figure 6C and 6D).

In addition, among the other 6 cancers, the high-risk scoring group of PRS also has a significantly worse prognosis than the low-risk group, and the difference is statistically significant (P < 0.05). (Supplementary Figure 3).

Clinical samples and data collection

RNA sequencing and clinical data were available from The Cancer Genome Atlas (TCGA: https://portal.gdc.cancer.gov/) and the Gene Expression Omnibus (GEO: https://www.ncbi.nlm.nih.gov/geo/) data portals. We adopted three datasets (GSE13213, GSE72094 and GSE68465) from the GEO database. We also populated a list of platelet-associated receptor genes from the GO Consortium database (http://geneontology.org/). The mutation spectrum and copy number variation (CNV) data for the entire TCGA dataset study population were also downloaded. Samples with less than 30 days of survival and samples with 0 expression values higher than 10% were removed.

Analysis of differentially expressed genes

We performed differential gene analysis on all transcription data, setting the false discovery rate (FDR) to < 0.05, and |fold change| > 1.5 as the cutoff value [31]. The differentially expressed genes were then extracted from all the differentially expressed genes. Functional enrichment analysis was performed through the GO and KEGG pathways.

Survival analysis

Survival related genes were selected by univariate and multivariate COX analysis, which was performed with R’s survival software package.

Building a prognostic index based on platelet-related genes

We performed a multivariate analysis of core platelet-related genes to construct a platelet risk score (PRS). This was based on the expression data multiplied by the Cox regression coefficient. Patients were divided into high-risk and low-risk groups based on the median. The prognostic value of PRS was evaluated in patients with different subtypes. The TIMER database was used to analyze the abundance of tumor-infiltrating immune cells. The association was then calculated between PRS and LUAD immune cell infiltration.

Gene set enrichment analysis

Gene set enrichment analysis (GSEA) is a computational method to identify over-represented biological processes and pathways.

Mutation and copy number variation analysis

The genes were analyzed with different mutations in different risk score subtypes (high and low risks). The mutation spectrum of the entire study population from the TCGA dataset was determined as previously described. Copy number variation (CNV) data were downloaded from GDAC Firehose and divided into datasets according to the risk score expression.

Gene expression and miRNA integration

To study microRNAs’ (miRNAs’) potential gene regulation, we focused on the abnormally expressed miRNAs (adjusted p-value < 0.05, log2 | multiple changes | > 1) and the RNAseq genes with statistically significant differences selected from high and low-risk patients. The correlation between miRNA and regulatory genes was obtained through TargetScan analysis (http://www.targetscan.org/vert_72/).

Statistical analysis

The AUC for the survival ROC curve was calculated by the survivalROC R software package to verify the prognostic signature’s performance. Differences between clinical parameters were tested using independent t-tests. A p-value of less than 0.05 was considered statistically significant.

Ethics approval and consent to participate

Our study did not require the approval of an ethics committee, as it was a secondary analysis of a public database.

Data availability statement

The data that support the findings of this study are openly available in TCGA at https://cancergenome.nih.gov/, and GEO at https://www.ncbi.nlm.nih.gov/geo/, reference number: GSE13213, GSE68465 and GSE72094.