Development and validation of a risk prediction score for patients with nasopharyngeal carcinoma

Background We aimed to develop and validate a predictive model for the overall survival (OS) of patients with nasopharyngeal carcinoma (NPC). Methods Overall, 519 patients were retrospectively reviewed in this study. In addition, a random forest model was used to identify significant prognostic factors for OS among NPC patients. Then, calibration plot and concordance index (C-index) were utilized to evaluate the predictive accuracy of the nomogram model. Results We used a random forest model to select the three most important features, dNLR, HGB and EBV DNA, which were significantly associated with the OS of NPC patients. Furthermore, the C-index of our model for OS were 0.733 (95% CI 0.673 ~ 0.793) and 0.772 (95% CI 0.691 ~ 0.853) in the two cohorts, which was significantly higher than that of the TNM stage, treatment, and EBV DNA. Based on the model risk score, patients were divided into two groups, associated with low-risk and high-risk. Kaplan–Meier curves demonstrated that the two subgroups were significantly associated with OS in the primary cohort, as well as in the validation cohort. The nomogram for OS was established using the risk score, TNM stage and EBV DNA in the two cohorts. The nomogram achieved a higher C-index of 0.783 (95% CI 0.730 ~ 0.836) than that of the risk score model 0.733 (95% CI 0.673 ~ 0.793) in the primary cohort (P = 0.005). Conclusions The established risk score model and nomogram resulted in more accurate prognostic prediction for individual patient with NPC.

(EBV) DNA, serum lactate dehydrogenase (LDH) [7], globulin (GLOB), hs-CRP [8], and neutrophil/lymphocyte ratio (NLR) [9]. The infection of EBV is virtually 100% associated with NPC in endemic areas. In addition, the plasma EBV DNA has gradually been used in clinical applications and is considered to be the most attractive potential biomarker to complement TNM staging system [10]. The derived neutrophils to leukocytes ratio (dNLR) has been linked to the inflammatory status and clinical outcomes among several cancers, including NPC [11]. This indicates that dNLR likely has prognostic value and also has the advantage of being inexpensive and easy to calculate. Hemoglobin (HGB) levels have been regarded as important determinants of outcome for a number of cancers treated with radiotherapy, especially gynecological tumors and NPC [12]. However, it remains a challenge to screen and incorporate biomarkers into a new staging system for NPC patients. A random forest model is an effective classifier as it can predict class of input, and select its most important features [13]. Nomograms have currently been proven to be an effective tool to predict prognosis of patients with cancers, including lung cancer [14], rectal cancer [15], and gastric cancer [16]. In this study, we utilized a random forest model to screen factors related to prognosis of NPC, and incorporate them into a new staging system by establishing a nomogram model.

Patients and clinical characteristic
We retrospectively reviewed the patients with first diagnosed NPC at Sun Yat-Sen University Cancer Center (SYSUCC) between January 2009 and December 2011. Patients included criteria: (1) All patients were pathologically diagnosed as NPC in SYSUCC for the first time.
(2) Patients were not any malignancies besides NPC. (3) Patients did not undergo anti-tumor therapy, and data was collected prior to anti-tumor therapy. (4) No patients were infected with hepatitis B virus or hepatitis C virus. (5) Complete records in the database of medical information, and follow-up data. All the patients were classified based on the 8th edition of the AJCC TNM staging guidelines.

Statistical analysis
Data analyses was carried out using SPSS standard version 20.0 (SPSS, Chicago, USA) and R software version 3.6.1 (http:// www.R-proje ct. org). The cut-off values were determined using the R package "survival" and "survminer". The Kaplan-Meier method was utilized to estimate OS of patients in high-risk and low-risk groups. The concordance index (C-index) was utilized to assess discriminative ability and predictive accuracy of established random forest model and nomogram. The C-index was calculated and compared using the "survcomp" package. The area under the curve (AUC) was computed using the "survivalROC" package. Calibration of the nomogram for 1-, 3-, and 5-year OS was executed via comparison of the predicted survival and observed survival. All statistical tests were two-tailed, and a P value less than 0.05 was considered to be statistically significant.

Patients and clinical characteristics
In total, 519 NPC patients were enrolled in this study. All patients were randomly divided into either the primary cohort (n = 363) or the validation cohort (n = 156). The patients' demographic data and clinical characteristics in both the primary cohort and validation cohort are described in Table 1. In the primary cohort, 209 (57.57%) patients were males and 154 (42.43%) were females. The mean age (SD) of patients in the primary cohort was 46.05 (10.87) years, and the median OS was 51.0 months (interquartile range [IQR]: 42.3-66.7 months). In the validation cohort, 92 (58.97%) patients were males and 64 (41.03%) were females. The mean age (SD) was 46.87 (11.58) years, while the median OS was 50.4 months (IQR: 41.7-66.0 months). In the primary cohort, the 1-, 3-, and 5-year OS rates were 95.0%, 84.0% and 46.8%, respectively. In the validation cohort, the 1-, 3-, and 5-year OS rates were 98.7%, 84.0% and 45.5%, respectively.

Model construction based on clinical characteristics
In the primary cohort, we utilized a random forest model to select the most significant variables from 34 different clinical features. In addition, we used the sliding windows sequential forward feature selection method (SWSFS) to identify important variables by minimizing the 'out of bag (OOB)' error rate (Fig. 1A). According to the minimum OOB, three variables including dNLR (HR = 1.14; 95% CI 1.05-1.23; P = 9.14 × 10 -4 ), HGB (HR = 0.98; 95% CI 0.97-0.99; P = 5.24 × 10 -3 ) and EBV DNA (HR = 1.59, 95% CI 1.32-1.93, P = 1.22 × 10 -6 ) were found to be significantly associated with OS among NPC patients (Fig. 1B). Finally, we constructed a risk score model that included dNLR, HGB and EBV DNA. The computational formula of the risk score was as follows: risk sco re = (0.466 × DNA) + (0.129 × dNLR) − (0.02 × HGB). The heatmap of NPC samples among the two cohorts are shown in Fig. 2, in which red represents upregulated imaging features and blue represents downregulated imaging features. The three feature clusters (C1-C3) were identified within the heatmap using unsupervised hierarchical clustering of 3 imaging features.

Model evaluation
ROCs were utilized to evaluate the accuracy of the established risk score model, TNM stage, treatment, and EBV DNA. In the primary cohort, for the 1-year OS (Fig. 3A), the AUC of TNM stage, treatment, EBV DNA and our established model were 0.748, 0.591, 0.751 and 0.797, respectively. Moreover, our model achieved a higher AUC than the TNM stage, treatment, EBV DNA for both the 3-year OS (Fig. 3B) and 5-year OS (Fig. 3C). In the validation cohort, for the 1-year OS (Fig. 3D) (Fig. 3E, F). The results of a time-dependent ROC  Table 2).

Performance of the established model in stratifying risk
Based on the computational formula of risk score (0.46 6 × DNA + 0.129 × dNLR − 0.02 × HGB), NPC patients were subdivided into high risk (risk score ≤ −0. 16) and low risk (risk score > −0.16). We used the R package "survival" and "survminer" in order to determine the cut-off values. The optimum cut-off of our model was − 1.46. The results demonstrated that patients with a high-risk score had significantly shorter OS than patients with a low-risk score in the primary cohort (P < 0.001; Fig. 5A) and the validation cohort (P < 0.001; Fig. 5E). In the primary cohort, the 1-, 3-, and 5-year survival probabilities of the high-risk were 90.0%, 71.3% and 37.3%, respectively. On the other hand, for the lowrisk patients, the 1-, 3-, and 5-year survival probabilities were 99.0%, 93.0% and 53.5%, respectively. Meanwhile, in the validation cohort, the 1-, 3-, and 5-year survival probabilities of the high-risk and low-risk patients were 97.2%, 70.4%, 35.2% and 98.8%, 95.3%, 54.1%, respectively (Table 3). Moreover, Kaplan-Meier curves demonstrated that high-risk and low-risk subgroups were significantly correlated with OS outcomes in the primary cohort (Fig. 5C, P < 0.001; Fig. 5D, P = 0.011), as well as in the validation cohort (Fig. 5G, P = 0.015,  Fig. 5H, P = 0.021) with a respective stage III, and stage IV, with an exception for patients in stage I/II (Fig. 5B,  F).

The nomogram for the prediction of OS
We established a nomogram for OS, which included risk score, TNM stage and EBV DNA in the two cohorts. In the primary cohort, the nomogram model achieved a C-index of 0.783 (95% CI 0.730 ~ 0.836), which was significantly higher than that of the prognostic model 0.733 (95% CI 0.673-0.793, P = 0.005) (Figs. 6A, 7A). On the other hand, in the validation cohort, the nomogram model achieved a C-index of 0.776 (95% CI 0.709 ~ 0.844), which was much higher than that of the prognostic model of 0.772 (95% CI 0.691 ~ 0.853, P = 0.455) (Figs. 6E, 7B). Calibration curves for probability of survival at 1-, 3-, and 5-years showed optimal agreement between prediction established in the nomogram and actual observation in the primary (Fig. 6 B-D) and validation cohort (Fig. 6F-H). Furthermore, RMS curves demonstrated a larger slope in the primary cohort for nomogram, which indicates superior estimation of survival with nomogram (Fig. 7A, B).

The correlations among the variables in the nomogram model
The relationship among variables of nomogram model were shown in Fig. 8. In this figure, blue indicated positive correlations, while red indicated negative corrections. Moreover, the correlation coefficients were proportional to color intensity and circle size. In the primary cohort, there was a highly significant between EBV DNA and risk score (Fig. 8A). Meanwhile, treatment was moderately correlated to TNM stage. In the meantime, we were able to get consistent results in the validation group (Fig. 8B).  Table 2 The C-index of the prognostic model, TNM staging, Treatment, and EBV DNA for prediction of OS in the training cohort and validation cohort C-index: concordance index; CI: confidence interval; P values are calculated based on normal approximation using function rcorrp.cens in Hmisc package

Discussion
The TNM stage is commonly utilized to predict prognosis and guide clinical therapeutic regimen across many cancers. This system for NPC was updated and refined to the 8 th edition in 2016 [2]. However, this system has several controversies as it is completely based on anatomical extent of cancer, and neglects the biological heterogeneity of NPC patients. Many other important risk factors need to be considered in the current staging systems.
In the present study, we used a random forest model to investigate prognostic value of many clinical factors and selected the most significant ones. We revealed that EBV DNA, dNLR and HGB levels could be used for prediction of NPC prognosis. The established risk score model, which included EBV DNA, dNLR, and HGB have a higher AUC and C-index than TNM stage, treatment, and EBV DNA model in both the primary and validation cohort. Based on risk score, we stratified NPC patients into two subgroups, including high-risk and lowrisk, which were significant in OS outcomes. Moreover, according to results from random forest model analysis, we established nomograms that can help predict OS in NPC patients, which integrated risk score, TNM stage, and treatment. The infection of EBV is common in NPC in endemic areas. Levels of plasma EBV DNA have been shown to be the most attractive potential biomarker to predict prognosis and provide accurate risk stratification in NPC [20]. Intriguingly, a recent prospective screening study that involved 20,174 participants showed that plasma EBV DNA detection was useful to screen for NPC, as it was associated with 97.1% sensitivity and 98.6% specificity. NPC was found to be detected significantly earlier by EBV DNA, with a significantly higher proportion of stage I or II disease than in a historical cohort (71% vs 20%), and had superior 3-year progression-free survival (97% vs. 70%; hazard ratio, 0.10) [21]. However, EBV DNA alone for prognosis has limitations, as the methodology of EBV DNA measurement is not globally standardized and measurement is not routinely available in many medical institutions. Moreover, there is accumulating evidence indicating that inflammation plays an important role in carcinogenesis and tumor proliferation [22]. Studies have reported that inflammation-based markers can be used as potential prognostic factors for many cancers, including CRP, neutrophils, lymphocytes, dNLR and ALB [23][24][25]. Neutrophils secrete proangiogenic cytokines, which include IL-8, MMP-9, MMP-8, and VEGF, which are known to contribute to tumor angiogenesis and progression [26,27]. Lymphocytes, especially the CD8+ T cells, which mediated immune response and increased OS of patients with gallbladder cancer [28]. Furthermore, low HGB is a risk factor in cancer patient survival, and HGB level is an important predictor in evaluation and treatment anemia [29,30]. In our study, we used a random forest model to identify that EBV DNA, dNLR and HGB levels can be used to predict NPC prognosis.
There are several limitations to this study. First, this is a retrospective study, and there may be a selection bias during data collection. Secondly, this is a singlecenter study with a limited number of NPC patients. Third, this study established models for predicting OS among patients with NPC. However, models of diseasefree survival are unknown. Therefore, our next aim is to validate our models on a large-scale with a multi-center study.