{'fraction_non_alphanumeric': 0.03469673464754625, 'fraction_numerical': 0.031669756706648, 'mean_word_length': 5.067493112947658, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 2, 'https://': 3, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Screening and discrimination of optimal prognostic genes for pancreatic cancer based on a prognostic prediction model

Abstract The prognosis of pancreatic cancer is poor because patients are usually asymptomatic in the early stage and the early diagnostic rate is low. Therefore, in this study, we aimed to identify potential prognosis-related genes in pancreatic cancer to improve diagnosis and the outcome of patients. The mRNA expression profile data from The Cancer Genome Atlas database and GSE79668, GSE62452, and GSE28735 datasets from Gene Expression Omnibus were downloaded. The prognosis-relevant genes and clinical factors were analyzed using Cox regression analysis and the optimal gene sets were screened using the Cox proportional model. Next, the Kaplan-Meier survival analysis was used to evaluate the relationship between risk grouping and patient prognosis. Finally, an optimal gene-based prognosis prediction model was constructed and validated using a test dataset to discriminate the model accuracy and reliability. The results showed that 325 expression variable genes were identified, and 48 prognosis-relevant genes and three clinical factors, including lymph node stage (pathologic N), new tumor, and targeted molecular therapy were preliminarily obtained. In addition, a gene set containing 16 optimal genes was identified and included FABP6, MAL, KIF19, and REG4, which were significantly associated with the prognosis of pancreatic cancer. Moreover, a prognosis prediction model was constructed and validated to be relatively accurate and reliable. In conclusion, a gene set consisting of 16 prognosis-related genes was identified and a prognosis prediction model was constructed, which is expected to be applicable in the clinical diagnosis and treatment guidance of pancreatic cancer in the future.


Introduction
Pancreatic cancer, which is the most fatal malignancy of the digestive system, is highly aggressive with a poor prognosis, leading to considerable morbidity and mortality worldwide and it continues to be a major health challenge (Forster et al. 2020;Nevala-Plagemann et al. 2020). The 5-year survival rate of pancreatic cancer is <5%, and 85% of patients die within 12 months of diagnosis because the asymptomatic early stages generally lead to late detection (Qian et al. 2019). Conventional risk factors such as advanced age, alcohol consumption, tobacco use, obesity, and history of diabetes mellitus and chronic pancreatitis are associated with the etiology of pancreatic cancer (Rawla et al. 2019). In addition, molecular biomarkers play an important role in the development and prognosis of pancreatic cancer. The discovery of effective biomarkers is essential for the detection of pancreatic cancer at an early stage, which may contribute to improving the prognosis and developing new therapeutic strategies.
Recently, numerous abnormally expressed genes have been identified in pancreatic cancer based on a microarray technology, which can detect subsets of genes that could be potential biomarkers for cancer diagnosis (Nonogaki et al. 2010). For instance, AKT serine/threonine kinase 1 (AKT 1) regulates growth factor-induced cell survival and its phosphorylated form (p-Akt1) is involved in the carcinogenesis of pancreatic cancer (Liu et al. 2010). This molecule has been detected at high levels in patients and is correlated with a lower primary tumor size and extent tumor (T) stage, which may be a favorable prognostic factor for pancreatic cancer (Liu et al. 2010). Song et al. (2020) found that protein arginine methyltransferase 1 (PRMT1) promotes pancreatic cancer growth and is predictive of poor prognosis. Saif et al. (2007) reported that P16 deletion may potentially inhibit cyclin D, CDK4, and CDK6 function and regulate cell-cycle progression. Furthermore, P16 deletion is significantly associated with shorter average survival times and acts as a predictive marker for poor prognosis in patients with pancreatic ductal adenocarcinoma (Luo et al. 2013). In addition, many studies suggest that the abnormal expression of growth factors or their receptors may affect cellular functions and the tumorigenicity of pancreatic cancer (Ebert et al. 1995;Gnatenko et al. 2018). It has been reported that high plasma levels of the transforming growth factor (TGF)-b1 in pancreatic cancer patients is associated with advanced stages of the disease and significantly shorter survival times (Javle et al. 2014).
Genes with variable expressions that are related to the prognosis of pancreatic cancer have been identified in many studies. However, screening for new indicators and the construction of prognostic models based on additional prognosis-related genes, appears to be a more sensitive strategy than using signal genes for the early detection of pancreatic cancer. Moreover, the predictive accuracy of these potential prognostic biomarkers has been largely limited by the sample size in previous studies. Therefore, in the present investigation, the mRNA expression profiles extracted from The Cancer Genome Atlas (TCGA) database were used to identify significant prognosis-related genes and clinical factors. Then, a risk prediction model for prognosis was constructed based on optimal genes and clinical factors, and three datasets (GSE79668, GSE62452, and GSE28735) from the Gene Expression Omnibus (GEO) database were downloaded to validate the correctness and reliability of the prognostic prediction models.

Data sources and data preprocessing
The mRNA-seq expression profile data of pancreatic cancer were downloaded from the TCGA database (https://gdc-portal.nci.nih. gov/), which is based on the Illumina HiSeq 2000 RNA sequencing platform and used as a training dataset in this analysis. A total of 168 samples, including 164 pancreatic cancer tumors and four normal samples, were included in this dataset. Gene expression profiles of GSE79668 (Kirby et al. 2016), GSE62452 (Yang et al. 2016), and GSE28735 (Zhang et al. 2012) were obtained from the National Center for Biotechnology Information (NCBI) GEO database (http://www.ncbi.nlm.nih.gov/geo/) (Barrett et al. 2005). The GSE79668, GSE62452, and GSE28735 profiles were used as three independent validation datasets and included information on 51, 69, and 45 pancreatic cancer tumors, respectively. The distributions of the clinical information for the training and validation datasets are listed in Table 1. The CEL files of the GSE62452 and GSE28375 datasets produced from Affymetrix platform were downloaded and the raw data were processed using the oligo package (version 1.44.0) in R3.4.1 (http://www.bioconductor.org/ packages/release/bioc/html/oligo.html) (Parrish and Spencer 2004). This process included background correction using the MAS method, supplementation of missing values using the median method, and quantile normalization (Irizarry et al. 2003). Furthermore, the gene count data of GSE79668 were acquired from the GEO database based on the Illumina HiSeq 2000 RNA sequencing platform (Supplementary File S1, https://www.ncbi. nlm.nih.gov/geo/query/acc.cgi?acc¼GSE79668), and the gene transcript sequence counts for all genes across all tissues using GENCODE v. 9. Table 1 Distributions of clinical information for training set and testing sets "-" indicates relevant information is missing.

Preliminary screening of expression variable genes
The expression variable genes were preliminarily identified based on the coefficient of variation (CV) value, which is a random variable defined as a ratio obtained by comparing the full range, mean deviation, standard deviation, and mean values (Koopmans et al. 1964). Both the dispersion and average value of the variable can affect the CV, which reflects the fluctuations in data size (Koopmans et al. 1964). The mRNAs were obtained after filtering genes with expression median values <1 in the samples retrieved from the TCGA database. The CV value of the gene expression of samples in the TCGA training dataset was calculated using the genefilter package (version 1.58.1) (Gentleman et al. 2009) in R3.4.1 (https://bioconductor.org/packages/release/bioc/ html/genefilter. html), based on a threshold CV value >0.7.

Prognosis-related gene and clinical factor identification
After screening for expression variable genes and their corresponding samples, the clinical information was extracted for follow-up analysis. The univariate and multivariate Cox regression analyses were performed to identify significant prognosisrelated genes and clinical factors using the R3.4.1 language survival package (version 2.41.3, http://bioconductor.org/packages/ survival/) . Statistical significance was set at P < 0.05. Following the extraction of the expression values of the variable genes from the TCGA database, a coupled two-way clustering analysis based on centered Pearson's analysis (Eisen et al. 1998) was conducted to identify similarities among these genes using the pheatmap package (version 1.0.8) in R3.4.1

Optimal prognostic gene identification
After inputting the matrix of the gene expression values obtained from previous screening steps, the associations between optimal genes and the prognosis of pancreatic cancer were analyzed using the Cox proportional hazards (PHs) model based on L1-penalized regularization in the penalized package (version 0.9.50, http://bioconductor.org/packages/penalized/) in the R3.4.1 language (Tibshirani 1997;Goeman 2010). The optimal parameter "lambda" in the screening model was obtained using 1000 crossvalidation likelihood (cvl) cycle calculations.

Risk prediction model construction based on optimal genes
The risk prediction model of the optimal genes was constructed based on the prognostic coefficient of the Cox-PH generated using L1-penalized regularization and the gene expression levels (clinical factors) in the TCGA training set. The prognosis index (PI) of each sample was calculated using the following formula: PI score ¼ P Coef gene/clinical Â Exp gene/clinical , where Coef gene/clinical represents the coefficient of the identified genes (clinical factors), Exp gene/clinical represents the expression level of the genes (clinical factors) in the TCGA training set, and the median PI score was used as the boundary to divide the samples into high-and lowrisk groups in the training dataset.  The Kaplan-Meier (KM) survival analysis in the R3.4.1 language survival package (version 2.41.3, http://bioconductor.org/ packages/survival/) was used to evaluate the relationship between risk grouping and the overall survival (OS) time of pancreatic cancer patients (Goeman 2010). Finally, three datasets were used to validate the correctness and reliability of the prognostic prediction models using survival and receiver operating characteristic (ROC) curve analyses (Hajian-Tilaki 2013).

Risk prediction model construction based on optimal genes and clinical factors
To identify the significant clinical factors, the Cox-PH model was also constructed again using the clinical information of the samples provided by the TCGA training dataset. We used the same method to construct a risk prediction model for optimal genes and clinical factors and calculated the PI score of each sample. Then, we divided the samples into high-and low-risk groups by setting the median PI score as the dividing point. Similarly, the association between risk grouping and the survival ratio of patients was analyzed again using the survival analysis, and the accuracy of prognosis of the training dataset was evaluated using an ROC curve analysis.

Preliminary screening of variable mRNA
Expression values for 13,238 mRNAs were retained after filtering out genes with median expression values <1 in the TCGA database. Then, 325 expression variable genes from 13,238 mRNAs were obtained using a CV > 0.7 as the threshold.

Screening for prognosis-related genes and clinical factors
The Cox univariate regression analysis identified 129 prognosisrelated genes from 325 expression variable genes. A set of 48 genes significantly associated with the prognosis of pancreatic cancer were identified using the Cox multivariate analysis (Table 2). Similarly, the following three independent prognostic clinical factors were identified: pathologic N (P ¼ 0.05), new tumor (P ¼ 0.01), and targeted molecular therapy (P ¼ 0.00) ( Table 3). The survival curve analysis showed that there were statistically significant differences in OS time between patients with and without targeted molecular therapy (logRank P ¼ 0.00, Figure 1B), pathologic N0 compared with pathologic N1 (logRank P ¼ 0.01, Figure 1C), and between patients with and without new tumors (log-rank P ¼ 0.01, Figure 1D).
Hierarchy cluster analysis was conducted according to the similarity of the expression values of 48 independent prognosisrelated genes. The results showed that samples in the TCGA database were, clearly, divided into clusters 1 and 2, which contained 97 and 67 pancreatic cancer samples, respectively ( Figure 1A). In addition, the distribution of clinical information of the two clusters was determined using the chi-square test, and the results are shown in Table 4. A statistically significant association between new tumors and the formation of clusters was observed (P ¼ 0.02). However, there was no significant correlation between the formation of clusters and any other clinical factors (pathologic N and targeted molecular therapy).

Screening of optimal prognosis genes
A Cox-PH model was constructed to identify the optimal genes related to the prognosis of pancreatic cancer, and the lambda score (19.33) was obtained from 1000 cvl cycle calculations ( Figure 2A). As shown in Table 5, of these 48 genes, 16 optimal genes (CCK, FABP6, FER1L6, G6PC2, GSTM1, HEPACAM2, KIF19, KLK8, MAL, MUC17, NLRP2, PI16, REG4, SFTPA2, TNNT2, and UPK1B) associated with the prognosis of pancreatic cancer were identified (Table 5). Furthermore, the distribution of the gene coefficients is shown in Figure 2, B and C.

Risk prediction model construction based on optimal genes
The training dataset included 164 patients who were divided into the high-and low-risk groups according to the median PI score (0.12). The mean OS time was 21.3 6 17.26 months and 12.42 6 10.73 months for the low-(82 of 164 samples) and highrisk (82 of 164 samples) groups, respectively according to clinical information. The survival curve analysis revealed a significant difference in OS time between the high-and low-risk groups (P < 0.01, Figure 3A). The area under the ROC (AUROC) curve analysis value was 0.99, indicating that our risk prediction model for prognosis might be relatively accurate and in agreement with real-life conditions in estimating the prognosis based on the samples ( Figure 3E).
Then, the risk model of the optimal genes was further validated using three validation datasets. The patient samples in the GSE62452 dataset were divided into high-and low-risk groups, and 32 patients in the low-risk group showed a longer mean OS time (23.98 6 16.02 months) than that of the 33 in the high-risk group (16.55 6 16.73 months). The patients in the GSE79668 dataset were also divided into high-and low-risk groups according to the median PI score (0.81), and the mean OS time of the 25 lowrisk patients was higher (34.04 6 29.63 months) than that of the 26 high-risk patients (20.33 6 21.06 months). There were significant differences in survival rates between the high-and low-risk pancreatic cancer groups represented in the GSE62452 ( Figure 3B) and GSE79668 datasets ( Figure 3C, both P ¼ 0.03). In addition, the AUROC curve analysis values for accurate discrimination of the prognosis of pancreatic cancer patients in the GSE62452 and GSE79668 datasets were 0.96 and 0.94, respectively.
The patients in the GSE28735 dataset were also divided into highand low-risk groups (21 samples each) according to the median PI score (0.44), and they showed mean OS times of 13.146 10.44 and 21.66 12.31 months, respectively. However, no significant correlation between risk grouping and the OS time of the pancreatic cancer patients was observed according to the survival curve analysis of this dataset (P ¼ 0.08, Figure 3D). The AUROC curve value for accurate analysis of the GSE28735 dataset was 0.91.

Risk model construction based on optimal genes and clinical factors
Overall, the following three clinical factors related to the prognosis of pancreatic cancer were identified from the TCGA training dataset: pathologic N, new tumor, and targeted molecular therapy based on the Cox-PH model ( Figure 4B), and the results were consistent with those of the Cox univariate and multivariate analyses. The lambda score was 0.55 after 1000 cvl cycle calculations ( Figure 4A). The result of the discriminant analysis of the two risk groups (82 samples each) using the clinical information showed that the mean OS time was longer for the low-risk group (19.02 6 15.02 months) than it was for the high-risk group, which showed a poor prognosis (14.93 6 11.68 months). The results of the survival curve analysis suggested that risk grouping was significantly associated with the prognosis of pancreatic cancer (P < 0.01, Figure 5A).
To further estimate the prognosis of pancreatic cancer, a risk model was constructed using the optimal genes and three clinical factors. The mean OS time of the low-risk group was higher (21.81 6 14.97 months) than that of the high-risk group  GSE28735 validation datasets assigned to low-risk group (black) compared with high-risk group (red). (E) ROC curve for all datasets based on optimal gene combination. Black, red, green, and purple curves represent TCGA training, GSE62452, GSE79668, and GSE28735 validation datasets, respectively.
(13.78 6 11.5 months, 82 samples each) and the KM curve analysis showed that there was significant difference in survival rate between low and high-risk groups (P < 0.01; Figure 5B), indicating that degrees of risk for these genes and factors were significantly related to the prognosis of pancreatic cancer. In addition, the AUROC value was 0.99 ( Figure 5C), reflecting the accuracy and reliability of the prognostic prediction system.

Discussion
A total of 325 genes with high expression variability were preliminarily identified between tumor and normal tissue samples. Furthermore, 48 genes and the following three clinical factors, pathologic N, new tumor, and targeted molecular therapy relevant to prognosis, were screened using data from the TCGA database. In addition, a risk prediction model containing 16 optimal genes related to prognosis was constructed and validated to be relatively accurate and reliable according to the results of the KM curve analysis of the GSE79668 and GSE62452 validation datasets. Among these genes, FABP6, MAL, KIF19, and REG4 were significantly correlated with the prognosis of pancreatic cancer and are expected to be further applied in the clinical diagnosis of pancreatic cancer and treatment guidance. The most significant prognosis-related gene in our data was the fatty acid binding protein 6 (FABP6), which is a highly conserved cytoplasmic protein that binds bile acids and participates in the metabolism of enterohepatic bile acid (Fisher et al. 2009). Bile acids have been reported to regulate cell growth and proliferation, and abnormal levels are implicated in hepatic inflammation and tumorigenesis (Li and Apte 2015). Thus, genes associated with bile acids are also likely to be involved in cancer pathogenesis. Ohmachi et al. (2006) demonstrated higher expression levels of FABP6 mRNA in colorectal cancer tissues than in adenomas and metastatic lymph nodes, and the tumor size was smaller with high expression levels of FABP6. In addition, colorectal cancer cells transfected with FABP6 showed weaker invasiveness and lower levels of apoptosis than mock cells did (Ohmachi et al. 2006). Fang et al. (2007) reported that a novel transcript of FABP6 may protect colon cancer cells from apoptosis through activation of the nuclear factor (NF)-jB pathway. In this study, we found that FABP6 may act as a favorable biomarker for the prognosis of pancreatic cancer according to the multivariate Cox PH regression analysis. Further comprehensive investigations of the mechanism mediating the  involvement of FABP6 in pancreatic cancer development are needed and will be included in our future work.
The Mal, T cell differentiation protein (MAL) is mainly localized in the endoplasmic reticulum and plasma membrane of T cells, and regulates T cell differentiation (Rancaño et al. 1994). MAL has been reported as a tumor suppressor gene, and its increased expression may be implicated in the reducing the pathogenesis of head and neck squamous cell carcinoma (HNSCC) via suppression of cell proliferation, invasiveness, and tumor growth (Cao et al. 2010). In addition, MAL expression is correlated with prolonged disease-free survival in breast and gastric cancers (Buffart et al. 2008;Horne et al. 2009). However, the predictive accuracy of MAL in cancer prognosis remains controversial. For instance, upregulated MAL mRNA levels were significantly associated with shorter OS and progression-free survival in a cohort of high-grade serious epithelial ovarian carcinoma patients (Zanotti et al. 2017). This previous finding is consistent with our observation that MAL, which has high expression variability, was identified as one of the optimal genes related to the poor prognosis of pancreatic cancer. Our observation suggests that MAL might play an important role in the etiology of pancreatic cancer.
Kinesin superfamily proteins (KIFs) facilitate the transport of mRNAs, protein complexes, and organelles in an ATP-and microtubule-dependent manner (Brendza et al. 2000). KIFs are essential for mitosis and meiosis (Vicente and Wordeman 2015) and any abnormalities in mitosis cause cell death, gene deletion, and even induce carcinogenesis (Yu and Feng 2010;Zhu et al. 2014). Thus, the detection of abnormal kinesin protein expression could be utilized as a biomarker for early tumor diagnosis or to predict the survival of patients with tumors. Numerous studies have reported that altered expression of kinesins is related to the development and progression of various human tumors (Corson et al. 2007;Taniwaki et al. 2007;Lukong and Richard 2008). For instance, altered KIF14 mRNA expression has been shown to be a prognostic indicator for patients with breast and lung cancers (Corson and Gallie 2006). KIFC1 is essential for the viability of certain extra-centrosome-containing cancer cells (Cai et al. 2009). Chen et al. (2017) reported that KIF19 is significantly associated with the prognosis of hepatocellular carcinoma. The results showed that KIF19 expression was significantly correlated with the prognosis of pancreatic cancer, which implies it has a role in carcinogenesis and may have clinical value as a biomarker.
Another significantly poor prognosis-related gene identified in our data was regenerating family member 4 (REG4), a member of the regenerating (REG) islet-derived family of proteins. REG4 has been demonstrated to be highly expressed in ovarian (Lehtinen et al. 2016), colorectal (Zhu et al. 2015), rectal (He et al. 2014), and gastric (Miyagawa et al. 2008) cancers, and acts as a predictive biomarker for early diagnosis and prognosis. High serum REG4 levels in pancreatic ductal adenocarcinoma at tumor stages IA-IIA are associated with a worse survival rate than that of patients with grade I tumors (Saukkonen et al. 2018). The REG4-induced epidermal growth factor receptor (EGFR)/AKT/cAMP-response element binding protein (CREB) signaling pathway is involved in macrophage polarization to the M2 phenotype, which may promote pancreatic cancer cell proliferation and invasion, and regulate the invasion of extra-pancreatic and lymph vessels (Ma et al. 2016). In addition, a study reported that REG4 promotes the proliferation and invasiveness of pancreatic cancer cells by upregulating the invasion-related genes matrix metalloproteinase 7 (MMP-7) and MMP-9 . These findings suggest that REG4 may serve as a useful biomarker for the prognosis of pancreatic cancer.
In addition, the three clinical factors, pathologic N, new tumor, and targeted molecular therapy relevant to the prognosis of pancreatic cancer were screened in this study. In addition, a prognosis prediction model based on optimal genes was constructed, which was verified in the independent validation datasets GSE79668 and GSE62452. However, the model was not verified in the GSE28735 validation dataset (P ¼ 0.08), which might have been due to the limited dataset sample size. Furthermore, a prediction model based on optimal genes and clinical factors was also constructed, and a significant difference was observed in the survival ratio between the high-and low-risk groups, and the AUROC value was 0.99, indicating the correctness and reliability of our prognostic prediction models.
This study used three datasets to verify the prediction model. However, this study had some limitations that are worth mentioning. First, although the trend of the result is consistent, no significant correlation between risk grouping and the OS time of pancreatic cancer was observed according to KM survival curve analysis in GSE28735 dataset. In addition, the comprehensive model was not validated because the clinical information on the pathologic N, new tumor, and targeted molecular therapy were not included in the three validation datasets at the same time. Therefore, more comprehensive data validation needs to be conducted in future studies. In addition, relevant experiments need to be performed to verify the multiple candidate targets identified in this study.

Conclusions
In summary, in this study, 325 expression variable genes and a gene set consisting of 16 prognosis-related genes were identified from data retrieved from the TCGA database. The risk prediction model for prognosis based on optimal genes was validated to be relatively accurate and reliable in the GSE79668 and GSE62452 validation datasets according to the results of the KM curve analysis. Among these genes, FABP6, MAL, KIF19, and REG4 were significantly associated with the prognosis of pancreatic cancer and can expected to be further applied in the clinical diagnosis of pancreatic cancer and in the guidance of its treatment.