Glycosylation of plasma IgG in colorectal cancer prognosis

In this study we demonstrate the potential value of Immunoglobulin G (IgG) glycosylation as a novel prognostic biomarker of colorectal cancer (CRC). We analysed plasma IgG glycans in 1229 CRC patients and correlated with survival outcomes. We assessed the predictive value of clinical algorithms and compared this to algorithms that also included glycan predictors. Decreased galactosylation, decreased sialylation (of fucosylated IgG glycan structures) and increased bisecting GlcNAc in IgG glycan structures were strongly associated with all-cause (q < 0.01) and CRC mortality (q = 0.04 for galactosylation and sialylation). Clinical algorithms showed good prediction of all-cause and CRC mortality (Harrell’s C: 0.73, 0.77; AUC: 0.75, 0.79, IDI: 0.02, 0.04 respectively). The inclusion of IgG glycan data did not lead to any statistically significant improvements overall, but it improved the prediction over clinical models for stage 4 patients with the shortest follow-up time until death, with the median gain in the test AUC of 0.08. These glycan differences are consistent with significantly increased IgG pro-inflammatory activity being associated with poorer CRC prognosis, especially in late stage CRC. In the absence of validated biomarkers to improve upon prognostic information from existing clinicopathological factors, the potential of these novel IgG glycan biomarkers merits further investigation.


Results
IgG glycan measurements. The IgG glycan analysis resulted in 23 directly measured IgG glycans structures, and 54 derived traits that represent common features shared among several measured glycans (galactosylation, sialylation, core fucosylation and the incidence of bisecting GlcNAc; Fig. 1, Supplementary Table 1) 22,24 . We restricted our survival analysis to those IgG glycan traits that were found to be robustly analysed. Robustness was calculated as follows. On each plate from the CRC cohort we put 3 standards that were biologically identical. Therefore, differences between measurements of standards are consequence of only experimental noise. We then Each IgG contains one conserved N-glycosylation site on Asn197 of its heavy chain. Different glycans can be attached to this site and the process seems to be highly regulated. UPLC analysis can reveal composition of the glycome attached to a population of IgG molecules by separating total IgG N-glycome into 24 chromatographic glycan peaks (GP1-GP24), mostly corresponding to individual glycan structures.
Scientific RepoRts | 6:28098 | DOI: 10.1038/srep28098 calculated the variance of standards only and the variance in the whole CRC population. "Robustness" is defined as the ratio of those two variances (Var(Stand)/Var(CRC)) * 100 (i.e. lower values indicate higher robustness) and represents the contribution of experimental variation in total variation. Thirty nine of the 77 glycan traits whose percentage of experimental variation was below 20% were included in the analysis (Supplementary Table 1). The correlation coefficients of the robust measured and derived glycan traits are presented in Supplementary Fig. 1.
The univariate glycan HRs for the whole sample complete case analysis are presented in Table 2 and  Supplementary Table 2 for all-cause mortality and in Table 2 and Supplementary Table 3 for CRC specific mortality. IgG glycans linked to mainly galactosylation were strongly associated with all-cause mortality and CRC mortality. In particular an increase in the percentage of agalactosylated structures (G0 n ) and a decrease in mono-and di-galactosylated structures (G1 n , G2 n ) was associated with poorer all-cause and CRC-specific mortality. Statistically significant associations were also observed for decreased sialylation and increase in the incidence of bisecting GlcNAc (Table 2). Results were similar when AJCC stage 4 patients were excluded from the analysis (Supplementary Tables 4 and 5). The minus logarithm of the q-values (FDR corrected p-values) of all 39 glycan traits for all-cause mortality and CRC-specific model III are presented in a Manhattan-like plot (Fig. 2). Stratified analysis by stage for all-cause and CRC-specific mortality is presented in Supplementary Tables 6 and 7. An increase in the percentage of agalactosylated structures (G0 n ) and a decrease in mono-and di-galactosylated structures (G1 n , G2 n ) was associated with poorer all-cause and CRC-specific mortality in stages 1, 2 and 3 (p-values from all-cause mortality models for G0 n : stage 1: 0.05, stage 2: 0.009 and stage 3: 0.01) but not in stage 4 (p-value for G0 n : 0.38). In contrast, decrease in sialylation and increase in incidence of bisecting GlcNAc were statistically significantly associated with all-cause and CRC-specific mortality only in stage 4 (p-values from all-cause mortality models for stage 4 for FGS/(FG + FGS): 0.003; FGS/(F + FG + FGS): 0.01; FG2S1/(FG2 + FG2S1 + FG2S2): 0.002; FBG2S1/(FBG2 + FBG2S1 + FBG2S2): 0.008; FBS1/FS1: 0.008). Finally, only in stage 2 disease IgG glycans linked to core fucosylation were associated with all-cause and CRC-specific mortality (Supplementary Tables 6 and 7).
Multivariate Cox regression clinical algorithms (including all the covariates of model III) showed good prediction of subsequent all cause (Harrell's C = 0.73, AUC = 0.75, IDI = 0.02 [as compared to model II that included AJCC stage, age and sex]) and CRC-mortality (Harrell's C = 0.77, AUC = 0.79, IDI = 0.04 [as compared to model II that included AJCC stage, age and sex]). Using glycans in addition to the clinical factors (that were selected by generalised boosted regression) did not lead to any statistically significant improvements for the whole sample analysis (Table 3) or after stage stratification (Supplementary Tables 8 and 9). This was reconfirmed by using Cox regression with L1 (LASSO) penalties on model parameters 25 , as there were no significant differences in the validation deviances of models with and without glycans both for the whole sample and stage-stratified designs. Similarly, predictions of the 5-year risk of death using the clinical factors stage, age, sex, BMI and CRP (e.g.

Figure 2. Manhattan-type plot of the association FDR corrected p values (q-values) of all 39 glycan variables for all-cause (red dots) and CRC-specific (blue dots) mortality adjusted for AJCC CRC stage, age, sex, time between sample and surgery, operation type, CRP and BMI (Model II).
Analysed glycans are plotted on the X-axis. Y-axis plots the -logarithm of the q values. The names of glycans, which associations achieved statistical significance after correction for multiple testing (q < 0.05), are named in the plot. Details of the association analysis results are presented in Table 2  When we stratified by stage, adding glycans to the clinical variables improved the prediction results, did not change them, or made them worse, depending on the stage of cancer and on the chosen models. We tested whether independently of the choice of a model class, adding glycans to clinical covariates would improve predictions of a model of the same class estimated on independent test data using cross-validation. We performed two instances of the paired Wilcoxon sign-rank test comparing models with and without glycans, including all the considered models (W), or including only the models of disparate classes (Wd) as discussed in Methods. We showed that there was no significant improvement in the prediction of the rapid progressors using glycans (in addition to the clinical factors) for stage 2 (p W ~ 0.99, p Wd ~ 0.98) as measured by cumulative (merged) AUC on the validation data (Supplementary Table 12). Similarly, for stage 3 the impact of the glycans was not consistent across the models, varied depending on the modelling assumptions, and was not significant overall (p W ~ 0.75, p Wd ~ 0.58; Supplementary Table 13).
On the other hand, there was a significant improvement in the prediction of the rapid progressors using glycans for stage 4 (Supplementary Table 14), with p W ~ 0.01, p Wd ~ 0.04, leading to the median gain in the test AUC of 0.08. Importantly, the inclusion of glycans in the models consistently resulted in the improved quality of predictions across the range of the considered models, and independently of whether the restricted or extended sets of clinical variables were used in the adjustments. The results were qualitatively similar for multiple repetitions of 10-fold cross-validation with random partitions into non-overlapping test folds, and independently of whether 10-fold or two-fold cross-validation was used to estimate the AUC on test data for the considered models. The best extended clinical model had the test AUC of 0.58, with the PPV of 0.35. The best model augmented with unfiltered log-transformed glycans had the test AUC of 0.66, with the PPV of 0.62. (Note: we acknowledge that since the choice of these best models uses the validation data, new external validations are needed to confirm the differences).

Discussion
In this study, we investigated the relationship between the IgG glycome composition in plasma from of CRC patients with survival outcomes. We applied univariate and multivariate statistical models to examine the associations between specific glycan changes and CRC-specific or all-cause mortality. IgG glycans linked to galactosylation, sialylation and bisecting GlcNAc were strongly associated with all-cause mortality and CRC mortality.
Since the method used to analyse glycans normalises measured data to the total glycome, this effectively measures glycome content per molecule of IgG and is not sensitive to changes in total IgG concentration. Multivariate Cox regression clinical algorithms showed good prediction of outcome for all cause and CRC-mortality, but using glycans in addition to the clinical factors did not lead to any statistically significant improvements. However, when we investigated the prediction of rapid progressors within each AJCC stage, there was an improvement in the prediction of the rapid progressors using glycans for stage 4. It is well established that glycosylation changes are involved in the aetiology of cancer, and specifically mark tumour proliferation and metastasis 26 . IgG is produced and secreted by CRC cells and the expression levels of CRC-tumour derived IgG correlated with clinical and pathological characteristics of the tumour (including stage) 27 . In particular it has been shown that expression of IgG was stronger in CRC tissues with TNM stage III-IV, than in those with TNM I-II. Similarly, in this study we observe different changes in IgG glycosylation status (levels of sialylation and incidence of bisecting GlcNAC) in late-stage disease and we see an improvement in the prediction algorithms using glycans in addition to clinical factors in AJCC stage 4.
IgG biological activity is influenced by its Fc-glycosylation. Each heavy chain of IgG carries a single covalently attached bi-antennary N-glycan at the highly conserved asparagine 297 residue in each of the CH2 domains of the Fc region of the molecule. The attached oligosaccharides are structurally important for the stability of the antibody and its effector functions 28 . In addition, 15-20% of normal IgG molecules also bear complex bi-antennary oligosaccharides attached to the variable regions of the light chain, heavy chain or both 18,29 .
Changes in IgG galactosylation, sialylation, bisecting GlcNaC and fucosylation have been previously reported in cancer studies (Supplementary Table 15). In particular, a decline in plasma IgG galactosylation has been observed in multiple myeloma, prostate, gastric, lung and ovarian cancers [30][31][32][33][34][35][36][37] . Previous studies investigated small sample sizes (< 100 cancer cases). This is the first time that similar changes in IgG galactosylation have been shown to be associated with CRC prognosis in a study with > 1000 CRC patients (Supplementary Table 15). It has been shown that the immune system can identify and destroy new tumour cells through cancer immunosurveillance, which functions as an important defence against cancer. A recent review on the natural innate and adaptive immunity to cancer has presented evidence from mouse models that B cells (which create and release IgG) are important in the surveillance of CRC 38 . Furthermore, it has been hypothesised that decreased IgG galactosylation leads to a greater pro-inflammatory antibody response 39,40 , which might influence cancer survival after diagnosis. Interestingly, galactosylation of IgG increases during pregnancy and reverts to normal levels after delivery, indicating that IgG glycome composition is prone to natural modifications 41 .
Similar to galactosylation, decreased sialylation of IgG also results in a pro-inflammatory IgG phenotype 40 . In this study we found that decreased sialylation was also linked to poorer prognosis, which replicated findings of two small studies on ovarian 37 and gastric cancer 33 . However, there was conflicting evidence in a study of multiple myeloma 42 , where increased sialylation of IgG was linked to higher risk of multiple myeloma. Therefore, through both decreased galactosylation and decreased sialylation, IgG in CRC patients with poorer prognosis had significantly greater pro-inflammatory properties (decreased galactosylation and sialylation) than CRC patients with better prognosis. Furthermore elevated occurrence of bisecting GlcNAc and lack of core fucose results in increased ADCC activity. In our study we found higher occurrence of bisecting GlcNAc in CRC patients of poorer prognosis, but IgG core fucosylation changes were associated with all-cause or CRC-specific mortality only in stage 2 CRC patients.
To our knowledge, this is the largest study to date examining the complexities of IgG glycan structure in relation to CRC prognosis. It includes prospective CRC cases from almost all hospitals in Scotland therefore is broadly representative of the colorectal cancer population. Cases were recruited as soon as possible after diagnosis to limit survival bias among those recruited and maximize the person-years of follow up. In addition, data relevant to the survival analysis were obtained from the Scottish registries General Register Office and the Scottish Cancer Registry (which are known to have high levels of data quality and data completeness) after linkage of our participants with their databases using the Community Health Index number. In addition, care was taken to determine the AJCC stage by experienced study clinicians reviewing individual computerized tomography scans and staging and metastasis information for every case. Information on date of diagnosis, date of sample taken and date of operation were collected and used in the univariate and multivariate models. Finally, we applied the state-of-the-art technology to measure the IgG glycosylation status in collaboration with the most experienced glyco-analytical laboratories in Europe.
This study has some limitations. One single measurement of IgG glycosylation status was taken after colorectal cancer diagnosis. Protein glycosylation status changes through time and is affected by factors such as age, operation and obesity. We have adjusted our analyses for all these factors, but it remains possible that glycosylation changes observed at the time of diagnosis do not reflect the glycosylation status at the time of death. This might also explain why IgG glycosylation status had more predictive power for late stage disease (since measurements were taken close to the end-event). However, the aim of this current work was to identify glycosylation changes that occur at colorectal cancer diagnosis and also investigate whether a measurement at time of diagnosis can act as a prognosis biomarker and therefore any biases due to single measurements are not relevant for these hypotheses. The stratified analysis leads to a reduction in the sample size, which may adversely affect the quality of the constructed predictive models. Finally, since our analysis included CRC patients of Scottish origin only, it may be possible that these findings will not be applicable to other ethnic groups. However, we have recently analysed the IgG glycosylation in 3 independent cohorts of Systemic Lupus Erythematosus (SLE) patients of different ethnicity and we observed very similar changes in IgG glycome composition in African Caribbean, Han Chinese, and Latin American Mestizo SLE patients despite known differences in SLE manifestation in different ethnic groups 43 . Therefore, there is a good reason to believe that the IgG glycosylation changes we observed in relation to CRC prognosis would be relevant to other ethnic populations too.
Some of the classification methods considered in this paper can be improved for handling imbalanced distributions of class labels. We are currently developing such extensions based on the advancements in imbalanced classification 44 , with the specific focus on predicting the prognosis for patients diagnosed at stages 1 and 2. The predictions may potentially be improved further by adapting recent approaches based on interaction networks 45 of glycans and clinical factors. The definition of rapid progressors may require a refinement, especially at the earlier stages of colon cancer. For example, at stage 1, rather than defining a rapid progressor to be someone who dies of CRC during the lowest tertile of the follow-up, it could be useful to consider the speed of progression to the more advanced stages. We are currently investigating a range of alternative definitions. However, the key finding of IgG gyclans as a predictor of rapid CRC progression at the later stages remains potentially of great importance and should ideally be replicated in an independent study population.
The plasma IgG glycan differences which we observed at the time of CRC diagnosis are consistent with significantly increased IgG pro-inflammatory activity being associated with poorer CRC prognosis, especially in late stage (stages 3 and 4) CRC. In the absence of validated biomarkers to improve upon prognostic information from existing clinicopathological factors the potential of these novel IgG glycan biomarkers merits further investigation. In particular, the improved predictive power in models including glycan factors in stage 4 patients is interesting. Currently, there are various strategies that are employed when using chemotherapy 46 in stage 4 disease. Therefore having a novel biomarker or prediction model that could help select patient groups that may have a better prognosis and a more indolent disease course would be useful as these patients could perhaps be offered chemotherapy agents with lower toxicity when compared to a more aggressive combination strategy. Furthermore, there is a great interest in novel immunotherapies in cancer and therefore it will be of useful to identify a more 'immunogenic' tumour phenotype based on identified IgG glyco-markers with a particular response to immunotherapy. Certainly to date the most encouraging results for immunotherapies (including PD-1 inhibitors) have been in tumours such as melanoma that are thought to be highly immunogenic and there is further interest in investigating mismatch repair deficient colon cancers which are often associated histologically with a heavy immune infiltrate. Finally, recent studies 47 demonstrated that IgG glycosylation changes are very dynamic and variable between individuals, thus longitudinal studies are needed to fully investigate the prognostic potential of IgG glycosylation changes in CRC.

Methods
Study population. The SOCCS study (1999-2006) is a case-control study designed to identify genetic and environmental factors associated with non-hereditary colorectal cancer risk and survival outcome. Approval for the study was obtained from the MultiCentre Research Ethics Committee for Scotland and Local Research Ethics committee, and all participants gave written informed consent. The study has been described in detail elsewhere 48 .
The present study comprises 1229 patients with pathologically confirmed colorectal adenocarcinoma, in whom we assayed IgG glycan levels after CRC diagnosis. Participants completed a detailed lifestyle questionnaire and a semi-quantitative food frequency and supplements questionnaire (http://www.foodfrequency.org). Blood was collected and transferred to the research centre within 72 h of sampling. Plasma was prepared by gentle ficoll hypaque gradient centrifugation of sodium EDTA tubes and 1.5 mL of each participant's plasma was stored at − 80 °C.
IgG glycans measurement and normalisation. All methods for IgG glycans measurement and normalisation were carried out in accordance with the approved guidelines as described as described previously 24 .
Purification of IgG: The IgG was isolated from plasma samples using 96-well protein G monolithic plates (BIA Separations, Ajdovščina, Slovenia). Briefly, 50-90 μ l of plasma was diluted 10× with 1× PBS, pH 7.4 and applied to the protein G plate. IgG was eluted with 0.1 M formic acid (v = 1 mL; Merck, Darmstadt, Germany) and neutralized with 1 M ammonium bicarbonate (Merck).
Release and labelling of IgG glycans was performed as described previously 22 . IgG was first denatured with the addition of 30 μ L 1.33% SDS (w/v) (Invitrogen, Carlsbad, CA, USA) and 10 min incubation at 65 °C. Subsequently, 10 μ L of 4% Igepal-CA630 (Sigma-Aldrich, St. Louis, MO, USA) and 1.25 mU of PNGase F (ProZyme, Hayward, CA, USA) in 10 μ L 5× PBS were added to the samples and incubated overnight at 37 °C to release N-glyans. The released N-glycans were labelled with 2-aminobenzamide (2-AB). The labelling mixture was freshly prepared by dissolving 2-AB (Sigma-Aldrich) in DMSO (Sigma-Aldrich) and glacial acetic acid (Merck) mixture (85:15, v/v) to a final concentration of 48 mg/mL. A volume of 25 μ L of labelling mixture was added to each N-glycan sample in the 96-well plate. Also, 25 μ L of freshly prepared reducing agent solution (106.96 mg/ml 2-picoline borane (Sigma-Aldrich) was added and the plate was sealed using adhesive tape. Mixing was achieved by shaking for 10 min, followed by 2 hour incubation at 65 °C. Samples (in a volume of 100 μ L) were brought to 80% ACN (v/v) by adding 400 μ L of ACN (J.T. Baker, Phillipsburg, NJ, USA). Free label and reducing agent were removed from the samples using HILIC-SPE. An amount of 200 μ L of 0.1 g/mL suspension of microcrystalline cellulose (Merck) in water was applied to each well of a 0.45 μ m GHP filter plate (Pall Corporation, Ann Arbor, MI, USA). Solvent was removed by application of vacuum using a vacuum manifold (Millipore Corporation, Billerica, MA, USA). All wells were prewashed using 5× 200 μ L water, followed by equilibration using 3× 200 μ L acetonitrile/ water (80:20, v/v). The samples were loaded to the wells. The wells were subsequently washed 7× using 200 μ L acetonitrile/water (80:20, v/v). Glycans were eluted with 2× 100 μ L of water and combined eluates were stored at −20 °C until usage. Hydrophilic Interaction Chromatography (HILIC)-UPLC: Fluorescently labelled N-glycans were separated by HILIC on a Waters Acquity UPLC instrument (Milford, MA, USA) with fluorescence detector set with excitation and emission wavelengths of 330 and 420 nm, respectively. The instrument was under the control of Empower 2 software, build 2145 (Waters, Milford, MA, USA). Labelled N-glycans were separated on a Waters BEH Glycan chromatography column, 100× 2.1 mm i.d., 1.7 μ m BEH particles, with 100 mM ammonium formate, pH 4.4, as solvent A and acetonitrile as solvent B. Separation method used linear gradient of 75-62% acetonitrile (v/v) at flow rate of 0.4 ml/min in a 25 min analytical run. Samples were maintained at 5 °C before injection, and the separation temperature was 60 °C. The system was calibrated using an external standard of hydrolyzed and 2-AB labelled glucose oligomers from which the retention times for the individual glycans were converted to glucose units. Data processing was performed using an automatic processing method with a traditional integration algorithm after which each chromatogram was manually corrected to maintain the same intervals of integration for all the samples. The chromatograms were all separated in the same manner into 24 peaks. In addition to 24 directly measured glycan structures, 53 derived traits were calculated as described previously 22 . These derived traits average particular glycosylation features (galactosylation, fucosylation, sialylation) across different individual glycan structures. Consequently, they are more closely related to individual enzymatic activities, and underlying genetic polymorphisms.
Glycan normalisation: Normalization of glycan measurements across samples, was performed by total area of chromatograms, where peak area of each glycan was divided by total area of corresponding chromatogram.
Survival and risk related parameters. Mortality outcomes were ascertained through linkage with the National Records of Scotland. Primary cause of death ("CRC" or "other") was assigned from death certificates separately by two researchers (concordance was > 99%). Time to event was measured from the date of diagnosis. Survival follow-up was censored at the date of death or at January, 31 2013, for participants who were not known to have died. Clinicopathological staging data was collected where possible (e.g. TNM is not feasible in patients who did not undergo surgery). Clinical records were reviewed and tumour site and multiplicity were determined from clinical and pathological records. Pre-operative staging imaging was collected through participating centres. Using the collated pathology, imaging and clinical data, tumour stage was assigned according TNM staging system and mapped onto the American Joint Committee on Cancer (AJCC) tumour-node-metastasis system (AJCC 1-4).
Blood was collected after CRC diagnosis (and after surgery). Median time to sampling was 5.4 months after the diagnosis (interquartile range, IQR: 3.2 to 8.3 months). Since illness and treatment may acutely affect IgG glycan levels and confound the analysis, a variable describing time from operation to blood collection and a variable determining the type of operation were created. Statistical analysis. Data were analysed using STATA (version 12.0) and R. We initially examined the association between IgG glycan levels (continuous and rank transformed) and CRC/all-cause mortality using Cox proportional hazards models. The non-CRC cause of death was right censored in CRC related mortality analysis. Tests of the proportional hazards assumptions and linearity on the log hazard rate scale were performed prior the analysis. Deviation from the proportional hazard assumption was noted for stage of disease variable in the overall analysis. No violation of assumption for any of the covariates was observed after stage stratification. Three models were applied in the whole dataset, in a dataset excluding stage 4 disease (to evaluate the glycan associations among patients with no-metastatic disease -AJCC stages 1-3): a crude model (Model I), a model where hazard ratios (HR) were adjusted for age at diagnosis, sex and stage of disease (Model II) and a model where HRs were adjusted for age at diagnosis, sex, stage of disease, body mass index (BMI), time from operation to blood collection, type of operation and CRP (Model III). P-values were adjusted for multiple testing using false discovery rate method (Benjamini-Hochberg procedure). Then we investigated whether the profiles of biomarkers differed across the stages of CRC by performing AJCC stage stratification.
We estimated the predictive value of a clinical-only Cox-regression algorithm for model II (which included age, sex, disease stage) and model III (which was adjusted for age, sex, disease stage, BMI and CRP level), by calculating the Harrell's C concordance coefficient, the time-dependant cumulative/dynamic Area under the ROC curve (AUC) as suggested by Chambless, L. E. and G. Diao 49 and the Integrated Discrimination Index (IDI defined as a difference in discrimination slopes) and compared this to an algorithm that also included glycan predictors. We ran this analysis in the whole data set and after AJCC stage stratification. The glycan variables included in the final model were selected by applying generalised boosted regression, which orders the variables by their relative importance, in 1000 bootstrap samples 50,51 , over the 10 inner training folds, and forward selection of ranked glycans by applying log-likelihood ratio test. The predictive value of the models was evaluated on independent samples using 10-fold cross-validation for all models except for AJCC stage 1 strata, where cross-validation was not possible due to the small number of events.
Finally, we performed classification analyses to a) predict whether a patient will die of CRC within 5-years from the time of diagnosis and b) to predict the rapid progressors within each stage. The analysis of the rapid progressors was pre-planned and motivated by the successful development of molecular signatures of prognostic biomarkers in other areas 52 ; our goal here was to evaluate the potential of igG glycosylation as prognostic markers for CRC. A rapid progressor was defined as someone who died of CRC and whose follow-up time was in the lower 1/3 rd of the patients to die of CRC in that stage of cancer, with the cut-off thresholds at 2.9, 2.4, and 1.3 years for stages 2-4 respectively. We applied several families of classification models (LASSO, nearest neighbours, sparse shrunken centroids (PAM), Support Vector Machines, naive Bayes, Decision Trees, and boosted stump classifiers), with and without stratification, with and without initial filtering on the training data, with and without log transformations of glycan expressions and clinical factors. The choice of the models was influenced by their popularity in biomarker studies, and their ability to address high-dimensional (large-p, small-n) problems via regularization or an explicit control for model complexity. More information about these estimators is presented in Supplementary Box 1, and the motivations for considering multiple classifiers for this problem are discussed in Supplementary Section on Model Comparison. All the results for this analysis were aggregated over 10 runs using 10-fold cross validation, where the validation folds were used neither for filtering nor for estimation of model parameters. Where necessary and appropriate, we also used 10 inner folds to estimate the stopping criteria or optimal value of hyperparameters (such as the regularization parameter for LASSO).
Then we estimated whether adding glycans to clinical covariates would improve the predictive performance of a model of the same class on independent test data; that is, we compared LASSO using clinical variables with LASSO using clinical variables and glycans, DTs using clinical variables with DTs using clinical variables and glycans, etc. This task is different from the association analysis, or from identifying specific glyco-clinical models outperforming a known baseline, where corrections for multiple tests are needed to control the probability of false discoveries. We applied the paired Wilcoxon sign-rank test comparing models with and without glycans, testing whether the difference in the cross-validated AUC of the clinical and glycol-clinical models is significantly different from zero. One limitation of this approach is the assumption of independence of the paired observations. (Whereas the choice of the models is independent by construction, it is difficult to ensure independence of the estimates of the performance on new, previously unused test data without constructing an extensive simulation study mimicking the joint distribution of the glycans and clinical variables, which goes beyond the scope of this paper). To address this limitation, we performed the test twice, first by considering all the models, and then by considering models of disparate classes, where we retained one model of each class and discarded SVM-linear, SVM-quadratic, and SVM-cubic differing in the degree of the polynomial used for the kernel construction. We performed the tests multiple times for different training/test fold partitions of the 10-fold and 2-fold cross-validation, and found that the results were qualitatively similar irrespectively of the specifics of cross-validation.