Prognostic nomogram based on the lymph node metastasis indicators for patients with bladder cancer: A SEER population‐based study and external validation

Abstract Purpose This study aimed to compare the prognostic value of multiple lymph node metastasis (LNM) indicators and to develop optimal prognostic nomograms for bladder cancer (BC) patients. Methods BC patients were obtained from the Surveillance, Epidemiology, and End Results (SEER) database between 2004 and 2015, and randomly partitioned into training and internal validation cohorts. Genomic and clinical data were collected from The Cancer Genome Atlas (TCGA) as external validation cohort. The predictive efficiency of LNM indicators was compared by constructing multivariate Cox regression models. We constructed nomograms on basis of the optimal models selected for overall survival (OS) and cause‐specific survival (CSS). The performance of nomograms was evaluated with calibration plot, time‐dependent area under the curve (AUC) and decision curve analysis (DCA) in three cohorts. We subsequently estimated the difference of biological function and tumor immunity between two risk groups stratified by nomograms in TCGA cohort. Results Totally, 10,093 and 107 BC patients were screened from the SEER and TCGA databases. N classification, positive lymph nodes (PLNs), lymph node ratio (LNR) and log odds of positive lymph nodes (LODDS) were all independent predictors for OS and CSS. The filtered models containing LODDS had minimal Akaike Information Criterion, maximal concordance indexes and AUCs. Age, LODDS, T and M classification were integrated into nomogram for OS, while nomogram for CSS included gender, tumor grade, LODDS, T and M classification. The nomograms were successfully validated in predictive accuracy and clinical utility in three cohorts. Additionally, the tumor microenvironment was different between two risk groups. Conclusions LODDS demonstrated superior prognostic performance over N classification, PLN and LNR for OS and CSS of BC patients. The nomograms incorporating LODDS provided appropriate prediction of BC, which could contribute to the tumor assessment and clinical decision‐making.


| INTRODUCTION
Bladder cancer (BC) is the most common malignant tumor of urinary tract. 1 Radical cystectomy with pelvic lymph node dissection (PLND) remains recommended treatment for patients with muscle-invasive bladder cancer (MIBC) and high-risk non-muscle invasive bladder cancer. 2,3 The rate of lymph node metastasis (LNM) is up to 25% for MIBC patients, which is significantly associated with tumor recurrence and lower survival rates. 4,5 However, the assessment for LNM is still unsatisfactory which is most widely estimated on the basis of the American Joint Committee on Cancer (AJCC) tumor-node-metastasis (TNM) staging system. 6 The AJCC N classification ignores the different number of dissected regional lymph nodes (LNs), which may decrease the accuracy of prognosis prediction. 7 In addition, not only has the consensus of the optimal extent of PLND not been reached, there are individual differences in the pattern of regional LNM. 8,9 Therefore, novel LNM indicators are urgently needed for better treatment strategies.
In recent years, several LN prognostic factors were proposed to estimate the prognosis of BC patients, including the positive lymph nodes (PLNs) and the lymph node ratio (LNR). 10,11 Furthermore, the log odd of positive lymph nodes (LODDS) was proved to be more reliable for multiple tumor types. [12][13][14] Jin et al. discovered that LODDS promoted better predictive efficiency for MIBC patients than AJCC N classification and LNR. 15 However, this study did not explore the ability for predicting cause-specific survival (CSS) or further establish available prognostic nomogram for clinical usage. To our knowledge, there is no report on the evaluation of novel LNM indicators for predicting both overall survival (OS) and CSS of BC patients.
The present study aimed to compare the prognostic values among different LN status factors of BC patients by analyzing data from the Surveillance, Epidemiology, and End Results (SEER) database. Subsequently, we established novel nomograms incorporating LODDS for predicting OS and CSS, and successfully validated them in internal and external validation cohorts.

| Data source
We collected patients from the SEER database (SEER*Stat version 8.3.9.2) of the National Cancer Institute (NCI) program, which is one of the most representative tumor databases and covers approximately 28% of the US population. 16,17 We obtained permission to access the database using the reference number 15912-Nov2020. In addition, we collected the genomic and clinical information of BC from The Cancer Genome Atlas (TCGA, https://tcga-data. nci.nih.gov/tcga/). Ethics approval for this study was exempt as all data were download from publicly databases.

| Study population
Patients diagnosed with primary bladder cancer (ICD-O-3/ WHO 2008: "Urinary Bladder") between 2004 and 2015 were enrolled into the study. The exclusion criteria for data extraction were (1) patients aged <20 years or ≥ 80 years at diagnosis; (2) patients with no surgery or undergoing local resection; (3) patients with no LN removed or with unclear examined LNs (ELNs) and PLNs; (4) patients with T0/Ta/ Tis classification or with unclear AJCC TNM stage and tumor grade; (5) patients with unclear survival data or survival time less than 1 month. Finally, we included 10,093 BC patients after cystectomy or radical cystectomy from SEER and partitioned them into training cohort (n = 7065, 70%) and internal validation cohort (n = 3028, 30%), using a random sampling method. Additionally, we obtained 107 BC samples from TCGA based on the exclusion criteria as an external validation cohort. Figure 1 displays the selection process of patients.

| Measurement of variables
We collected variables of patients including age at diagnosis, gender, T/N/M classification, tumor grade, the amount of regional ELNs and PLNs, survival time and status. The tumor stage was based on the sixth edition of AJCC staging system, which was adapted to SEER-derived patients diagnosed from 2004 to 2015. Low-and high-grade tumors in TCGA were considered as grade I-II and III-IV in SEER database, respectively. The LNR was defined as the ratio of the amount of PLNs and ELNs. The LODDS was formulated by log [(PLNs +0.5) / (ELNs -PLNs +0.5)]. To avoid division by zero error, we added 0.5 to both numerator and denominator.
The primary and secondary endpoints were OS and CSS, which were shown as "COD to site recode" and "SEER cause-specific death classification" in SEER K E Y W O R D S bladder cancer, LODDS, lymph node metastasis, nomogram, prognosis, SEER database, respectively. Of note, the data from the SEER and TCGA were downloaded on December 9, 2021.

| Independence of LN status indicators
We selected the clinicopathological predictors with univariable Cox regression analyses through "survival" R package for OS and CSS in training cohort. To further assess the predictive values, each LN status factor (including N classification, PLN, LNR and LODDS) was integrated into multivariate regression models together with other risk variables with p < 0.05 in univariate analyses, respectively.

| Comparison of predictive performance among LN status indicators
Backward stepwise selection (via "MASS" R package) was utilized into the above models using Akaike Information Criterion (AIC) as the stopping rule, respectively. 18 The model with minimum AIC was selected as the optimal one. The predictive efficiency of these filtered models incorporating different LN factors was compared using AIC, bootstrapped concordance index (C-index) and area under the curve (AUC) via "riskRegression" R package.

| Construction and validation of nomograms
Variables included in the filtered models with the highest accuracy were integrated to develop nomograms for predicting OS and CSS in training cohort (via "rms" R package), respectively. The efficiency of nomograms was evaluated by bootstrapped C-indexes, time-dependent AUCs and calibration plots in training, internal validation and TCGA cohorts. Additionally, we preformed decision curve analysis (DCA) to assess the net benefit and clinical performance of the nomograms using the "ggDCA" R package. 19

| Survival risk classifiers established by nomograms
The multivariate Cox regression formulas of the nomograms for OS and CSS formed in training cohort were applied into patients in three cohorts using "nomogram-Formula" R package. All patients were divided into high-and low-risk groups according to the total points calculated via "survminer" R package. The Kaplan-Meier (K-M) method was used to assess the survival difference of OS and CSS between two risk groups.

| Functional enrichment
Differential expression genes (DEGs) from TCGA were searched for by comparing high-and low-risk groups with the threshold of |log fold change| >1 and p < 0.05 (via "limma" R package). We performed enrichment analyses for Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) with Metascape (https://metas cape.org/gp/index.html) to further verify the biological function enrichment. 20 The number of min overlaps, min enrichment and P value cutoff were set at 3, 1.5 and 0.01, respectively. Gene Set Enrichment Analysis (GSEA) software (Version 4.2.0) was also used to explore the function difference between two risk groups. 21 Hallmark gene sets were analyzed and the number of permutations was set at 1000. Normalized Enrichment Score (NES) > 1.5 and p < 0.05 were considered statistically significant.

| Estimation of tumor immune infiltration
To estimate the tumor microenvironment (TME), the stromal and immune scores of each sample in TCGA cohort were analyzed through the "estimate" R package. 22 We calculated the relative proportion of 22 human immune cell subsets by the CIBERSORT algorithm through "pre-processCore" R package. 23 Additionally, we also investigated the relationship between the risk classifiers and the genes related to immune checkpoint inhibitors (ICIs) via violin plot visualization.

| Statistical analysis
Continuous variables with non-normal distribution were reported as median with interquartile range (IQR) and categorical variables were presented as frequencies with percentages. Statistical significance was achieved with a two-sided p value less than 0.05. All statistical analyses were performed using R software (version 4.0.5).

| Patient characteristics and survival
The patients' characteristics of training, internal validation and TCGA cohorts are shown in  Figure S1.

| Prognostic analyses for OS and CSS
The detailed results of the univariate Cox regression analyses in training cohort are demonstrated in Table 2 We further performed multivariate analyses and generated prognostic models including different LN indicators, respectively. Briefly, N classification, PLN, LNR and LODDS were all independent risk factors for OS (Table 3) and CSS (Table 4). Besides, age, T and M classification were independent prognostic factors for OS; gender, T and M classification were independent risk factors for CSS.

PLN, LNR and LODDS
The comparison of LN status indicators in training cohort is shown in Table 5. After backward stepwise selection for the above prognostic models, the AICs of the filtered models containing LODDS for OS (65458.57) and CSS (48732.61) were lowest than that of other LN factors (Table S1). The C-indexes of the filtered models containing LODDS were higher over that of N classification, PLN and LNR for OS (0.705) and CSS (0.727). Meanwhile, the 1-, 3-and 5-year AUCs of the selected model including LODDS were higher than that of others except for the 5- . Taken together, the results indicated that the selected models containing LODDS had more enhanced predictivity for OS and CSS, and LODDS might act as the strongest predictive indicator over N classification, PLN and LNR.

| Construction and validation of nomograms
We developed nomograms based on the selected models containing LODDS in training cohort. As results, age, LODDS, T and M classification were incorporated into final nomogram for predicting OS ( Figure 2A); gender, tumor grade, LODDS, T and M classification were included in nomogram for CSS ( Figure 3A). Calibration plots of the three cohorts are displayed in Figure 2B-D and Figure 3B-D, indicating excellent agreement between predictive and actual observation for OS and CSS. The time-dependent AUC values of the nomograms for OS ( Figure 4A) and CSS ( Figure 4E) demonstrated more stable accuracy and better predicting efficiency over time. Besides, the DCA curves of nomograms in three cohorts for OS ( Figure 4B-D) and CSS ( Figure 4F-H) are displayed. X-axis represented threshold probabilities, and net benefit was plotted on the y-axis. The net benefit achieved by making decision on the basis of predictive models balanced the risk of real mortality with the unnecessary intervention of false prediction. The DCA curves indicated that the nomogram added more net benefits compared both the treat-all-patients scheme and treat-none scheme. We simultaneously evaluated detailed C-indexes and 1-, 3-and 5-year  AUCs values of nomograms in each cohort (Table S2). The above results revealed both appreciable reliability and clinical practicality of the prognostic nomograms.

| Survival risk classifiers based on nomograms
To further verify the performance of nomograms, we divided patients into high-and low-risk groups based on the total points calculated by nomograms for OS and CSS.
The cutoff values of the total points were 100.57 for OS and 82.94 for CSS, respectively ( Figure S2). K-M curves showed significant distinction in survival outcomes stratified by the risk classifiers for OS and CSS in three cohorts ( Figure 5).

| Pathway enrichment analyses
After standardization among the RNA-seq in TCGA, we extracted 78 DEGs based on the risk classifier for OS. The representative statistically enriched pathways (both GO and KEGG) were clustered together and shown with a network plot ( Figure 6A). The results showed that the DEGs were associated with immunoglobulin complex, extracellular matrix (ECM) organization, immunoglobulin production, and so on. In additional, GSEA analyses was performed using Hallmark database, illustrating that high-risk group was positively associated with epithelial mesenchymal transition, adipogenesis and protein secretion ( Figure 6B). In contrast, the DNA reparation pathway was enriched in low-risk group. Meanwhile, we also selected 84 DEGs based on the risk classification for CSS, and the results of enrichment analyses are shown in Figure S3A,B.

| The relationship between risk classifiers and tumor immunity
We investigated the TME and calculated stromal and immune scores for each BC sample based on the risk classifier for OS. The higher stromal score was observed in samples of high-risk group ( Figure 6C). We further investigated the differential relative proportion of immune infiltrating cells between two groups. The results showed the highrisk group was more abundant in M0 macrophages and plasma cells, whereas T follicular helper cells had higher infiltration in low-risk patients ( Figure 6E). The correlation between risk classifier and ICI-related genes were subsequently analyzed, and the expression of CTLA4 and PDCD1 were higher in low-risk group ( Figure 6D). The relationship between immune infiltration and the risk classifier based on CSS are presented in Figure S3C-E.

| DISCUSSION
Despite improvement in diagnosis and surgery, bladder cancer remains a heterogeneous malignancy with poor prognosis. 6,24 Individual patients within same AJCC stage vary widely in survival outcomes. Additionally, LNM is of critical importance to the prognosis of BC, while the N classification based on AJCC staging system is insufficient to accurately evaluate LNM. 4,25 Therefore, there is an urgently need to find novel LN status indicators to assess LNs involvement and stratify BC patients for individualized management. Considering this situation, several modified LN status factors has been proposed to predict the survival of BC, including PLN, LNR and LODDS. 10,11,15 Our study demonstrated that all four LN indicators were independent prognostic factors for OS and CSS in BC. We further compared the predictive values among N classification, PLN, LNR and LODDS. As results, selected models containing LODDS with minimal AICs, maximal C-indexes and AUCs were considered as optimal prognostic models for OS and CSS, indicating that LODDS had better discrimination capability in predicting prognosis of BC patients.
The N classification of AJCC staging system evaluates LN status mainly by detecting in certain regional areas, which was unsatisfactory due to the unclear number of negative and removed LNs. 26 an unfavorable risk factor for cancer patients. 28 Although more LN dissection could lead to a smaller probability of undiscovered PLNs, the threshold for LN dissection or the extent of PLND varied considerably and had not achieved consensus at present. 8,29,30 LNR has been shown to reflect the progression of disease and to be a prognostic factor for BC. 4,8,31 However, the association between LNR and tumor prognosis remains controversial when all or none examined LNs exhibited metastasis, especially with insufficient regional LNs removed. 32 Compared with the above, LODDS could reflect both the absolute number of PLNs and negative LNs, indicating that it could stratify patients categorized as N0 classification (no regional LNM) and patients with all examined LNs metastasized. 29,33,34 The current study suggested that the prognostic models with LODDS had best predictive capability for OS and  CSS, which was similar to Jin et al.'s study for OS of BC patients. However, we paid more attention to clinical utilization and constructed two novel nomograms incorporating LODDS for predicting OS and CSS of BC patients for the first time. Then, the nomograms were further validated in internal and external validation cohorts. The calibration curves showed stable linearity and appropriate efficacy of the nomograms, and the calculated C-indexes and timedependent AUCs were most above 0.70 and 0.75 in three cohorts, respectively. In terms of clinical utility, the DCA curves demonstrated that the nomogram showed consistent larger net benefit across a broad range of threshold probabilities. We believed that the nomograms had satisfactory applicability in predicting survival for BC patients. Taken together, better predictive accuracy and clinical validity were verified in our nomograms compared with AJCC and other LN status systems. Furthermore, risk classifiers for OS and CSS were established on basis of the total scores from nomograms, stratifying BC patients into different risk groups. Patients in two high-risk groups had worse survival in each cohort. We subsequently explored the biological function between the two groups. The results showed that there were major differences in immunoglobulin production and secretion, ECM organization, epithelial mesenchymal transition, adipogenesis and DNA reparation. It has been reported that the immunoglobulin, produced by B lymphocytes and plasma cells in response to immunogen, was associated with the prognosis in solid tumors. 35,36 ECM is a fundamental node of cell-extrinsic metabolic regulation and the organization of ECM was related to tumorigenesis and migration. 37,38 The other enrichment pathways were also involved in the progression of multiple tumors. [39][40][41] In addition, the TME and the sensitivity to ICI were evaluated in two groups. Noticeably, high-risk patients had higher stroma scores, consistent with previous studies in BC patients. 42,43 The stromal cell alteration was previously shown to affect the development of organ-specific metastasis. 44 The extracellular vesicles derived from tumor cells could be transferred to adjacent organs and internalized by stromal cells, developing a metastatic-designated microenvironment. Furthermore, specific stromal cells like cancer-associated fibroblasts contribute to the upregulation of multiple cytokines and the remodeling of ECM, which could facilitate the tumor invasion and metastasis. 45 The analysis of immune cell infiltration suggested that the high-risk group had a higher infiltration of M0 macrophages and plasma cells, while T follicular helper

Risk Threshold
Net Benefit cells were more abundant in low-risk group. The high proportion of M0 macrophages, indicating a higher level of inflammatory activation, was correlated with poor prognosis of BC patients. 46,47 It has previously been reported that more recruitment of naïve M0 macrophages is associated with early LNM, which could prevent T cells from attacking tumor cells and secrete growth factors promoting tumor angiogenesis. 48,49 T follicular helper cells might contribute to the maintenance of a protective immune response, which was proved to be associated with better prognosis in BC. 50,51 Furthermore, we observed higher expression levels of CTLA4 and PDCD1 in low-risk group, suggesting BC patients with low risk points based on OS classifier might benefit more from immunotherapy. 52 Collectively, these results explained the possible mechanisms from the perspective of tumor biological and immune function and further corroborated the validity of our risk classifiers for OS and CSS. Indeed, as a retrospective study, we need prospective and multicenter study to verify the prognostic value of LODDS. Secondly, the details of surgical approach, such as the extent of LN dissection at specific nodal level are not recorded in detail, which is worth further investigation. Despite these limitations, our study proves better predictive value of LODDS and firstly incorporates it into prognostic nomograms for both OS and CSS in BC patients.   (43) 380 (33) 290(25) 229 (20) 186 (16) 156 (14) 132 (11) 113(10) 9 3 (8) 68 (6)  4 5 (4)  2 6 (2) 8 (1) 1876 (100)    AUTHOR CONTRIBUTIONS Shuai Li: Conceptualization (lead); data curation (lead); formal analysis (lead); investigation (lead); methodology (lead); project administration (lead); visualization (lead); writing -original draft (lead); writing -review and editing (lead). Yicun Wang: Conceptualization (supporting); data curation (supporting); formal analysis (supporting); investigation (supporting); methodology (supporting); project administration (supporting); validation (supporting); writing -review and editing (supporting). Xiaopeng Hu: Conceptualization (supporting); data curation (supporting); formal analysis (supporting); investigation (supporting); methodology (supporting); project administration (supporting); supervision (supporting); validation (supporting); writing -review and editing (supporting).