Identification of a prognostic signature in colorectal cancer using combinatorial algorithm‐driven analysis

Abstract Colorectal carcinoma is one of the most common types of malignancy and a leading cause of cancer‐related death. Although clinicopathological parameters provide invaluable prognostic information, the accuracy of prognosis can be improved by using molecular biomarker signatures. Using a large dataset of immunohistochemistry‐based biomarkers (n = 66), this study has developed an effective methodology for identifying optimal biomarker combinations as a prognostic tool. Biomarkers were screened and assigned to related subsets before being analysed using an iterative algorithm customised for evaluating combinatorial interactions between biomarkers based on their combined statistical power. A signature consisting of six biomarkers was identified as the best combination in terms of prognostic power. The combination of biomarkers (STAT1, UCP1, p‐cofilin, LIMK2, FOXP3, and ICOS) was significantly associated with overall survival when computed as a linear variable (χ 2 = 53.183, p < 0.001) and as a cluster variable (χ 2 = 67.625, p < 0.001). This signature was also significantly independent of age, extramural vascular invasion, tumour stage, and lymph node metastasis (Wald = 32.898, p < 0.001). Assessment of the results in an external cohort showed that the signature was significantly associated with prognosis (χ 2 = 14.217, p = 0.007). This study developed and optimised an innovative discovery approach which could be adapted for the discovery of biomarkers and molecular interactions in a range of biological and clinical studies. Furthermore, this study identified a protein signature that can be utilised as an independent prognostic method and for potential therapeutic interventions.

. Frequency distribution of the immunostaining scores in primary colorectal cancer and normal colon mucosa for antibodies that have not previously been published Figure S2. Treemap representation of the main pathways associated with 29 biomarker targets included in the combinatorial analysis Figure S3. Correlation matrix of biomarkers included in combinatorial analysis. Figure S4. Kaplan-Meier survival analysis of the biomarker signature in the validation cohort Figure S5. Distribution of the pathological risk groups across the prognostic groups of biomarker signature Figure S6. Distribution of UICC stage, T stage, lymph node stage and MMR status across the prognostic groups of the signature Table S1. Clinicopathological characteristics of patients in the discovery cohort (Grampian cohort), their tumours and the relationship of each variable with overall survival Table S2. Clinicopathological characteristics of patients in the validation cohort, their tumours and the relationship of each variable with overall survival Table S3. Biomarkers screened on the discovery cohort of colorectal cancer Table S4. Peptide sequences used as immunogens to generate monoclonal antibodies which have not been previously tested on the Grampian patient cohort Table S5. The numbers of normal and tumour tissue samples in the multi-tissue microarray Table S6. Dilutions, antigen retrieval conditions, subcellular localisation of antibodies which have not previously been published Table S7. Comparison of the expression of each protein in normal colonic mucosa and different UICC stages in primary colorectal cancer using antibodies that have not previously been published on the discovery cohort Table S8. Starting list of biomarkers included in the biomarker combinatorial analysis Table S9. Associations of the biomarker signature and clinicopathological variables. Survival analysis of the biomarker signature was stratified by clinicopathological variables Table S10. Comparisons of the clinicopathological characteristics of cases for which current pathological risk classifications differ from the risk groups of the biomarker signature Remark checklist

Supplementary materials and methods
The histopathological processing of colorectal tissue specimens All the colorectal cancer were received fresh immediately after resection and were dissected, fresh tumour and normal samples taken for the Grampian Biorepository, fixed and sampled for histopathology using a standardised protocol. All the specimens used to construct the tissue microarray were processed and fixed in 10% neutral buffered formalin for at least 48 hours at room temperature. Representative tissue blocks were embedded in wax before tissue slides were cut and stained with haematoxylin and eosin (H&E) for histopathological diagnosis. After examining each H&E stained slide an expert gastro-intestinal pathologist (GIM) selected the appropriate areas of tissue to be sampled for constructing the tumour microarray.

Colorectal cancer tissue microarray
The tissue microarray was constructed using 50 normal colon samples (distance of 10 cm from the tumour) and 650 primary colorectal tumours. The tumour resections were from patients who had undergone elective surgery with curative intent for primary colorectal cancer in the following periods: 1994-1998 (n=99), 1999-2003 (n=198) and 2004-2009 (n=353). After identifying suitable areas to be sampled from representative blocks, two cores (1mm in diameter) were taken from each sample and placed in the tissue microarray block.

Immunohistochemistry procedure
The tissue slides were firstly processed as follow: dewaxed in xylene for a minimum of 10 minutes and rehydrated for two minutes by immersion in decreasing ethanol concentrations.
For antibodies that did not need antigen retrieval, the slides were then placed in water ready for immunostaining. For the remaining slides, antigen retrieval was performed by microwaving (800W) for 20 minutes while fully immersed in citrate buffer (pH 6.0) or in EDTA (pH 7.8). The primary antibodies were applied on the slides and incubated for 60 minutes at room temperature. Blocking with peroxidase was applied for 7 minutes after the primary antibody was washed twice using Dako buffer washes. The blocking buffer was washed twice before applying the secondary antibody (peroxidase-polymer labelled goat antirabbit/mouse, Envision, Dako) for 30 minutes at room temperature. This was followed by further two washes before adding the diaminobenzidine substrate which was incubated for 7 minutes and washed off with water. The slides were counterstained with haematoxylin for 10 seconds and washed using tap water. The slides were then dehydrated in increasing ethanol concentrations for couple of minutes before being immersed in a xylene and mounted.
Antibody diluent was used instead of primary antibodies as a negative control. Appropriate normal tissue known to express the relevant protein was used as a positive control.

Immunohistochemistry assessment
The immunostaining results were assessed by either a semi-quantitative scoring system or a quantitative semi-automated image analysis. The semi-quantitative scoring system was conducted using light microscopy by two independent observers unaware of clinical data.
Tumour cores were scored based on the intensity of immunostaining of none (score 0), weak (score 1), moderate (score 2) and strong (score 3).
The automated quantitative image analysis was completed using QuPath, Slidepath or Aperio image software (supplementary material, Table S3). The scores (continuous) were based on H-score which is calculated using the following formula: 3x% strongly stained cells+ 2x% moderately stained cells + 1x% weakly stained. Therefore, to allow for meaningful comparisons with the semi-quantitative categorised data, the quantitative scores were dichotomised into four categories (scores 0, 1, 2 and 3) based on 3 cut-off points (quartiles). Categorical score 0 represented cases with H-score which fall under the first quartile, categorical score 1: between the first and second quartiles, categorical score 2: between the second and third quartiles and score 3: higher than the upper quartile.
The immunostaining was assessed in the tumour epithelial compartment (PTEN and VAV1 were assessed in both the stromal and epithelial tumours), except for the immune markers which were assessed in the tumour stroma (supplementary material, Tables S7 &   S8). The scores for each biomarker were recorded in Excel before being combined into one master spreadsheet and transferred to SPSS for analysis.

Development of monoclonal antibodies
Antibodies were produced using hybridomas technology and ELISA screening [15]. Short amino acid sequences (8-12 in length) were used as immunogens and were chosen using a range of opensource bioinformatic software [15]. The sequences of peptides and their location on each protein are specified in supplementary material, Table S3. The specificity and sensitivity of all antibodies, developed in our lab, were evaluated by immunoblotting as previously described [15]. Antibodies were also validated by immunohistochemistry using comprehensive set of multi-tissue microarray (supplementary material, Table S4). Each antibody was tested without antigen retrieval and with antigen retrieval solution. Antigen retrieval was performed by microwaving the tissue slides for 20 minutes while fully immersed in 10mM citrate buffer (pH 6.0) or in EDTA solution (pH=7.8). The immunostaining profile for each antibody was assessed by an expert consultant pathologist (GIM) and interpreted in the context of existing knowledge of the expression of corresponding protein.

Selection criteria of prognostic biomarker signatures
After running the iteration algorithm, solutions were manually evaluated based on the following criteria: Kaplan-Meier plot, median survival of prognostic groups, chi-square value, p-value, hazard ratio, the number of patients in each prognostic group, and finally evaluate the prognostic significance of the biomarker combination against clinical prognostic parameters), the number of biomarkers in the solution combination (i.e., a smaller number is preferred) and the concordance index using internal and external validation sets. The combinatorial programme by default evaluates the hazard ratio as a fitness score. To evaluate other fitness scores, the function "compute fitness" need to be adjusted. This is done by changing cph1.hazard_ratios_ to cph1. hazards_ (coefficients) or to cph1.summary.

Data processing and statistical analysis
The following Python libraries and modules were used: NumPy, Pandas, lifelines, scikitlearn, scikit-survival, SciPy, random, collection, itertools, os and logging.
All quantitative scores were dichotomised into four categories based on three cut-off points (i.e., first, second and third quartiles) to allow meaningful comparison with the remaining semi-quantitative data. There was no other transformation of data or any computations of missing values. The following statistical tests were used: Mann-Whitney U test, Wilcoxon signed rank test, chi-squared test, Spearman correlation, Kaplan-Meier survival analysis (log-rank test), Cox's proportional hazard model (univariate and multivariate) and hierarchical cluster analysis (Ward's method). Cox regression was used to calculate the hazard ratio and estimated coefficient of biomarkers and to compute the concordance index score of fitted models. Train-test split option was implemented for evaluating the performance of prognostic models.

Immunostaining profiles of new biomarkers
All antibodies showed immunoreactivity in primary colorectal cancer. Of those tested on normal colonic mucosa, only CK7, SPATA2 and cytoplasmic BCLAF were negative in all normal samples. The majority of biomarkers were localised to the tumour cell cytoplasm except for SPATA2L, SPATA25, BCLAF1, HMGB1, SNTB2 and PRDM16 which showed both nuclear and cytoplasmic immunoreactivities. HER2 was localised to the cell membrane while CDX2 was nuclear. The expression of PTEN, VAV1 and y-H2AX was recorded in both epithelial and stromal compartments, whereas the immunostaining of the remaining biomarkers was limited to tumour cells. There was no evidence of tumour heterogeneity between duplicate cores as shown in whole section immunohistochemistry of a sub-set of tumours [10].
Trends of increased and decreased expression were both observed when comparing primary tumour to normal colonic mucosa and when comparing different tumour stages (supplementary material, Figure S1 and Table S8). The immunostaining was significantly higher in primary colorectal cancer compared with normal colonic mucosa for BAG3, BCLAF1c, DRAM1, DRAM2, ERO1lb, PRDM16c, SLA2, CYP2J2, GCSP, yH2, SNTB2c, HMGB1c, EphB6, SPATA2L, SPATA2, TNFRSF18s and UCP2. Whereas the expression of P53, P21, p-Akt, NOTCH, PRDM16n, PTEN, SNTB2n, HMGB1n, SPTA25n, SPATA25n, STAT1 and VAV1 was significantly less in primary tumour samples. The majority of primary tumours (63%) were CK20+/CK7-/CDX2+, with 35% expressing non-classic phenotype in at least one of the three biomarkers and 2% of tumours showing the opposing immunophenotype of CK20-/CK7+/CDX2-( Figure S1). A proportion of 6% of tumour were HER2 positive. The frequency distributions of each immunostaining category for biomarkers scored using semi-quantitative method are shown in supplementary material, Figure S1. Figure S1. Frequency distribution of the immunostaining scores in primary colorectal cancer and normal colon mucosa for antibodies that have not previously been published. The frequency of negative, weak, moderate and strong immunostaining for biomarkers in normal colonic mucosa and colorectal cancers in the discovery cohort.
8 Figure S2. Treemap representation of the main pathways associated with 29 biomarker targets included in the combinatorial analysis.

Figure S3. Correlation matrix of biomarkers included in combinatorial analysis.
Chi-Square correlation was used to compare relationships between biomarkers. The figure shows four strongly-correlated biomarker subsets; cytochromes P450, immune markers, metabolic and cytoskeleton/signalling.

Figure S4. Kaplan-Meier survival analysis of the biomarker signature in the validation cohort.
Survival plots of prognostic groups identified through the biomarker signature in the validation cohort. Group 1 was used as the reference group for comparison with the other groups (groups 2, 3, 4 and 5). Pathological risk classification: low risk (pT1-pT3 and (G1 or G2) and V0 and pN0) vs high risk ((pN0 and pT4 or G3 or V1) or either pN1 or pN2).       Pathological risk classification: low risk (pT1-pT3 and (G1 or G2) and V0 and pN0) Vs high risk ((pN0 and pT4 or G3 or V1) or any pN1 or pN2). * Fisher exact test expected value of cells< 5). Significant values are highlighted in bold. Pathological risk classification: low risk (pT1-pT3 and (G1 or G2) and V0 and pN0) Vs high risk ((pN0 and pT4 or G3 or V1) or any pN1 or pN2).

MATERIALS AND METHODS Patients
2 Describe the characteristics (e.g., disease stage or co-morbidities) of the study patients, including their source and inclusion and exclusion criteria. 6,7 and 8 (Table S1) 3 Describe treatments received and how chosen (e.g., randomised or rule-based). 6 ,7 and 8 Specimen characteristics 4 Describe type of biological material used (including control samples) and methods of preservation and storage. 6 ,7 and 8, Materials and methods S1 Assay methods 5 Specify the assay method used and provide (or reference) a detailed protocol, including specific reagents or kits used, quality control procedures, reproducibility assessments, quantitation methods, and scoring and reporting protocols. Specify whether and how assays were performed blinded to the study endpoint.
6, Materials and methods S1, Table  S3-S6 Study design 6 State the method of case selection, including whether prospective or retrospective and whether stratification or matching (e.g., by stage of disease or age) was used. Specify the time period from which cases were taken, the end of the followup period, and the median follow-up time.
6-8 7 Precisely define all clinical endpoints examined. 6 8 List all candidate variables initially examined or considered for inclusion in models. Table S3, S7 9 Give rationale for sample size; if the study was designed to detect a specified effect size, give the target power and effect size.

Statistical analysis methods
10 Specify all statistical methods, including details of any variable selection procedures and other model-building issues, how model assumptions were verified, and how missing data were handled.

9-1
11 Clarify how marker values were handled in the analyses; if relevant, describe methods used for cutpoint determination. 9-11 RESULTS Data 12 Describe the flow of patients through the study, including the number of patients included in each stage of the analysis (a diagram may be helpful) and reasons for dropout. Specifically, both overall and for each subgroup extensively examined report the numbers of patients and the number of events.

27
13 Report distributions of basic demographic characteristics (at least age and sex), standard (disease-specific) prognostic variables, and tumour marker, including numbers of missing values.

Table S1
Analysis and presentation 14 Show the relation of the marker to standard prognostic variables. Table 1&2 15 Present univariable analyses showing the relation between the marker and outcome, with the estimated effect (e.g., hazard ratio and survival probability). Preferably provide similar analyses for all other variables being analysed. For the effect of a tumour marker on a time-to-event outcome, a Kaplan-Meier plot is recommended.

Figure 2 &3
16 For key multivariable analyses, report estimated effects (e.g., hazard ratio) with confidence intervals for the marker and, at least for the final model, all other variables in the model. Table 1 &2 17 Among reported results, provide estimated effects with confidence intervals from an analysis in which the marker and standard prognostic variables are included, regardless of their statistical significance. Table 1 &2 18 If done, report results of further investigations, such as checking assumptions, sensitivity analyses, and internal validation.

DISCUSSION
19 Interpret the results in the context of the pre-specified hypotheses and other relevant studies; include a discussion of limitations of the study.