Integrative analysis from multi‐centre studies identifies a function‐derived personalized multi‐gene signature of outcome in colorectal cancer

Abstract Colorectal cancer (CRC) is highly heterogeneous leading to variable prognosis and treatment responses. Therefore, it is necessary to explore novel personalized and reproducible prognostic signatures to aid clinical decision‐making. The present study combined large‐scale gene expression profiles and clinical data of 1828 patients with CRC from multi‐centre studies and identified a personalized gene prognostic signature consisting of 46 unique genes (called function‐derived personalized gene signature [FunPGS]) from an integrated statistics and function‐derived perspective. In the meta‐training and multiple independent validation cohorts, the FunPGS effectively discriminated patients with CRC with significantly different prognosis at the individual level and remained as an independent factor upon adjusting for clinical covariates in multivariate analysis. Furthermore, the FunPGS demonstrated superior performance for risk stratification with respect to other recently reported signatures and clinical factors. The complementary value of the molecular signature and clinical factors was further explored, and it was observed that the composite signature called IMCPS greatly improved the predictive performance of survival estimation relative to molecular signatures or clinical factors alone. With further prospective validation in clinical trials, the FunPGS may become a promising and powerful personalized prognostic tool for stratifying patients with CRC in order to achieve an optimal systemic therapy.

staging system. However, the TNM staging system is not sufficient for treatment decisions and prognosis prediction of patients with CRC. [2][3][4] For example, patients with stage IIB CRC tend to show poor prognosis, with a 5-year relative survival rate of 46%-61% compared with those with stage IIIA (~70%). 5,6 This limitation of TNM staging indicates an increasing and urgent need for identifying novel biomarkers to improve the outcome and treatment of patients with CRC.
With the application of muti-omics technologies to study CRC, it was demonstrated that CRC is of high heterogeneity at the intertumoral and intratumoral levels. Patients with CRC often have variable prognosis and treatment responses, even in tumours that are histologically identical. 7 Advances in molecular profiling have enabled a better understanding of CRC development and provide additional clinically relevant prognostic information beyond the current classic staging system. 8 During the past years, considerable substantial efforts have been made to identify gene expression-based biomarkers for predicting prognosis in patients with CRC. Although these efforts were an important step to guide clinical decision-making, few recently proposed signatures have been incorporated into clinical practice. A previous systematic review performed comparison analysis and revealed limited overlap across different signatures and poor prognostic performance across different independent datasets. 9 The potential issues preventing the translation of in silico data into clinical practice include (a) the fact that these existing signatures were usually generated from a small sample size or a single dataset; (b) did not account for biological heterogeneity and technical biases across different datasets, which led to overfitting and concentrating on the discovery dataset; and (c) insufficient independent validation. These limitations highlighted the requirement for adequate sample size and multi-institutional patient cohorts for sufficient statistical power when trying to identify a robust prognostic signature. 10 Another major concern for low reproducibility is that these signatures only focus on statistical values and often fail to incorporate the biological rationale, thus leading to the incorporation of unrelated genes, which are correlative rather than causative. 11,12 Therefore, it is necessary to explore novel personalized and reproducible prognostic signatures based on multi-institutional cohorts of patients with CRC of sufficient size to aid clinical decision-making and improve the outcomes of patients with CRC.
The present study combined large-scale gene expression profiles and clinical data of patients with CRC from multi-centre studies and developed a novel computational framework to identify a function-derived personalized gene signature (FunPGS) for improving outcome. This prognostic signature was validated in multiple independent datasets across different technology platforms, and its performance was assessed in comparison to recently proposed signatures. TA B L E 1 CRC patient datasets enrolled in the study three patient datasets profiled with different platforms and outcome measure were utilized as an independent testing cohort for validating the prognostic value of the signature. Detailed information about these 10 CRC datasets is shown in Table 1 and Table S1.

| Pre-processing of profiling data
Raw microarray datasets (.CEL files) from the GEO database were pre-processed and normalized using the Robust Multi-array Average (RMA) algorithm for background correction, log 2 -transformation and quantile normalization using the R package 'affy'. All microarray probes were mapped to Entrez Gene ID, and the mean value of multiple probes mapping to the same gene ID was used to represent the gene expression level using the R package 'limma'. To avoid systematic measurement bias, each microarray datum was normal-

| Gene set function analysis
Gene Ontology (GO) function enrichment analysis of the prognostic gene sets was performed using the R package 'clusterProfiler'. 13 GO terms with P < 0.01 were considered to be significantly different and were selected for further analysis. Significantly enriched GO terms were clustered and visualized using the Enrichment Map plugin in Cytoscape. 14 Functional similarity among different sets of enriched GO terms was computed using the R package 'GOSemSim'. 15

| Development of a function-derived personalized gene signature
The workflow describing the development of the FunPGS is illustrated in Figure 1. The FunPGS was developed based on an improved version of the single-sample gene set enrichment analysis (ssGSEA) scoring method, 16,17 as follows: In the above equations, G represents the prognostic gene set, G protective represents the protective prognostic gene set, G risky represents the risk-associated prognostic gene set, S represents the single patient sample, r g j represents the rank of gene j in the prognostic gene set and N G is the number of prognostic genes. Patients with a high FunPGS displayed better outcomes than those with a low FunPGS.

| Statistical analysis
Survival analysis, including univariate and multivariate analyses with Cox proportional hazards regression and Kaplan-Meier analysis with log-rank test was performed using the R package 'survival'. Hazard ratios (HRs) and 95% confidence intervals (CIs) were calculated.
Harrell's concordance index (C-index) was calculated for each dataset to evaluate its prognostic performance using the R package 'survcomp'. Time-dependent ROC curves and AUC at 3-and 5-years were calculated to assess the predictive performance of molecular signature in comparison with clinical prognostic factors (stage and age) using the R package 'timeROC'. Meta-analysis was performed using Schematic representation of the computational workflow to derive and validate a function-derived gene signature as a personalized prognostic predictor of outcome in patients with colorectal cancer the R package 'meta'. Because of potential heterogeneity among clinical samples and studies, heterogeneity among studies was assessed using the Higgins' I 2 and Q statistics. When heterogeneity existed among studies (P < 0.1 and I 2 > 50%), random-effect models were used for the meta-analyses. Otherwise, the fixed-effect model was employed. Similarity of gene membership between each prognostic signature pair was assessed by the Jaccard index. All statistical analyses were performed using R (v3.3.3) and Bioconductor.

| Consistency evaluation of prognostic gene sets from different patient datasets
The association of genes with survival was first assessed in each dataset of the meta-training cohort using univariate Cox regression analysis followed by multivariate analysis adjusted by clinical variables including stage, gender and age. This resulted in 0. The blue diamond shows the fix-effects meta-analysis summary of HRs over the seven datasets (HR = 0.67, 95% CI = 0.59-0.76, P < 0.001). C-index, Harrell's concordance index; CI, confidence interval, HR, hazard ratio low reproducibility of prognostic genes from a single patient dataset, as shown previously. 9,18 Then, GO enrichment analysis was performed for each prognostic gene set, and the functional similarity among enriched GO terms was calculated for each prognostic gene set pair. Compared with random signatures, these prognostic gene sets derived from a single patient dataset exhibited significantly high functional consistency ( Figure 2B).

| External validation of the FunPGS on three independent datasets across different technology platforms
The ability of the FunPGS for stratifying patients with CRC and different prognosis was first confirmed in two validation datasets with microarray platforms. For the GSE14333 dataset, disease-free sur-  Figure 4E).
The C-index for the FunPGS was calculated in each of the three independent datasets. As shown in Figure 4F, the FunPGS yielded a significant prognostic value, exhibiting a high C-index in all the independent datasets (GSE14333, C-index = 0.59, P = 0.043; GSE33113, C-index = 0.77, P = 0.004; and TCGA, C-index = 0.76, P < 0.001) ( Figure 4F).

| Independence of the FunPGS from other clinicopathological factors
To evaluate whether the FunPGS was independent from other clinical or pathological factors, multivariate Cox regression analysis was conducted in the meta-training cohort and three independent datasets using the following factors as categorical variables: FunPGS

| Prognostic performance of the FunPGS in comparison with five previously published gene signatures
In the present study, five gene expression-based signatures associated with outcome of patients with CRC were collected

| Integrated prognostic signature obtained by combining the FunPGS with clinical factors
The AJCC TNM staging system and age are well-known important prognostic factors. 24,25 In the multivariate analysis conducted in the present study, TNM stage and age remained significantly associated with survival besides the FunPGS in the meta-discovery and TCGA cohorts, indicating the independent and complementary value of TNM stage and age. Therefore, the FunPGS, TNM stage and age were combined to construct an integrated molecular and clinical factors-based prognostic signature (IMCPS), and its prognostic value was tested. The IMCPS was quantified by subjecting the FunPGS, TNM stage and age to a multivariate Cox regression model in the meta-discovery cohort as follows: (−0.0001 × FunPGS) + (0.7792 × TNM stage) + (0.0236 × age). The median score of the IMCPS derived from the meta-discovery cohort was used as a cut-off value to stratify patients into low-or high-risk groups. As shown in Figure 7, the IMCPS had higher separation than the FunPGS alone between the high and low-risk groups in the meta-training cohort and the  Figure 7A and 7). The present study also compared the IMCPS with TNM stage and age by time-dependent ROC analysis and observed that the AUC of the IMCPS was significantly higher than TNM stage and age at 3-year and 5-year OS in both the meta-training cohort and the independent TCGA dataset.
Thus, the IMCPS greatly improved the predictive performance of survival estimation (Figure 7C and 7). This may be partly because of heterogeneous prognosis and treatment responses, which may arise from molecular heterogeneity. 8,26 Therefore, reliable and independent molecular biomarkers capa- patients with CRC at the individual level. Furthermore, the FunPGS represents differential activity levels of key biological processes and pathways rather than simply different expression levels of individual genes. 16,17 The FunPGS was constructed in a meta-training cohort comprising 1192 patients from seven datasets, which overcomes the weakness of small sample size and a single simple source. Limitations in the application of previously proposed signatures have demonstrated the importance of rigorous validation and reproducibility for medical applications. 28 The FunPGS effectively discriminated patients with CRC and significantly different survival and was successfully validated in a completely independent TGCA RNA-seq dataset, suggesting that the FunPGS is robust and insensitive to technical biases that are inherent to measurements by different platforms. It is well known that the outcome of patients with CRC has improved by the combined treatment of surgical resection and adjuvant chemotherapy. However, because of molecular heterogeneity, a group of patients with a relatively low incidence of recurrence may be (WIBEZD2017009-05). The funders had no roles in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

CO N FLI C T O F I NTE R E S T
The authors declare that they have no competing interests.