Ewing sarcoma (ES) is a rare, high-grade tumor affecting primarily young adults and adolescents, with an incidence of 0.13 per 100,000. It is genetically characterized by chromosomal translocations involving the EWSR1 gene, most commonly EWS-FLI1 t(11;22), which is confirmed by in situ hybridization with EWSR1 break-apart probes. Despite combined-modality treatment incorporating dose-intensive chemotherapy, surgery, and radiotherapy, only 60-70% of patients achieve 5-year survival. There are currently no prospectively validated tissue biomarkers in clinical use for prognostic stratification.
Tumor heterogeneity is the biological phenomenon where cells within a single tumor differ in genomic, epigenomic, and proteomic profiles. Intra-tumoral heterogeneity (variation between cells within a tumor) is mechanistically linked to the rapid emergence of treatment resistance. Despite this, most current biomarker assessments in clinical pathology remain semi-quantitative (scored as +1, +2, +3), lack cellular-level resolution, and fail to quantify heterogeneity as a feature. Ki67 proliferative index is the only biomarker showing potential prognostic significance in Ewing sarcoma in large series, but it has been assessed using subjective IHC scoring.
The study was motivated by the need for unbiased, high-dimensional, cell-level quantification of biomarkers in tumors. The authors developed an integrated pipeline combining: automated fluorescent confocal image acquisition, novel image segmentation algorithms to define individual cell nucleus and cytoplasm boundaries, probability density function-based feature summarization for each patient, and Random Survival Forest (RSF) machine learning to identify prognostic classifiers without requiring arbitrary threshold setting.
Ewing sarcoma tissue microarrays (TMAs) were assembled from three European cohorts sourced from University of Valencia (Spain), University of Dusseldorf (Germany), and Leiden University Medical Centre (Netherlands), all participating in the EuroBoNeT biobank. All patients were treated with standardized chemotherapy, radiotherapy, and surgical protocols (EICESS92/EE99 CESS81). The three cohorts were characterized by different TMA core sizes: Cohort a had one 2mm core per case, Cohort b had three 2mm cores, and Cohort c had two 1mm cores.
Of 524 initially identified patients, 472 had TMA cores. After excluding 52 cases with no cores listed on obtained TMAs and 77 with missing cores on slides (non-adherent), 314 potential biopsies were imaged. Following strict image quality control, 67 cases (36%) were eliminated due to poor DAPI localization or inadequate biomarker preservation, leaving 118 informative cases (Cohort a: 43, Cohort b: 16, Cohort c: 59) with 185 total cases having both visible cores and outcome data for combined RSF analysis.
Clinical characteristics of the 118 imaged patients: mean age 18-19 years across cohorts, male-to-female ratio approximately 1.2-3:1, survival rates varying by cohort (Cohort a: 70% alive, Cohort b: 25% alive, Cohort c: 54% alive). Cohorts differed in tumor site distribution, tumor volume, and metastasis rates. Importantly, Kaplan-Meier analysis confirmed that the 118 imaged cases were representative of the original 185 combined cohort in terms of survival distribution, validating that quality control filtering did not introduce survival bias.
Tissue sections (5 mm thick) were processed for simultaneous multiplex immunofluorescence with five markers: DAPI (nuclear DNA marker), CD99 (Ewing sarcoma cytoplasmic marker), Ki67 (proliferative index), and signaling biomarkers EGR1, Foxo3a, pS6 (phospho-ribosomal protein S6, downstream of PI3K), and pMAPK (phospho-MAPK). All antibodies were monoclonal rabbit primary antibodies from New England Biolabs, detected with secondary antibodies conjugated to Alexa 488, 555, and 647 fluorophores.
Antibody specificity was validated in Ewing sarcoma cell lines using pharmacological inhibitors: rapamycin (10 nM) for pS6 signal reduction, LY294002 (10 mM) for Foxo3a nuclear localization, and U0126 (10 mM) for EGR1 cytoplasmic localization. Antibodies were validated by titration to confirm signal above background. Images were acquired with an Olympus Fluoview FV1000 confocal microscope at 60x oil objective (NA 1.35), capturing 2048x2048-pixel images at 211x211 mm field with 1.292 mm slice thickness. Between 1 and 6 images per patient were captured depending on core size. Total processing time was 2-3 days per TMA slide.
The signal-to-noise quality of confocal images is critical for automated segmentation. Median filtering and Gaussian smoothing were applied to all images to reduce noise. An automated quality control classifier using 12 image features (including focus score, inverse coefficient of variation, image correlation from GLCM, and power log-log slope across 8 spatial frequency ranges) was built to objectively distinguish good and poor quality images, though its ROC AUC of 0.72 in both classifier variants indicates that fully automated QC remains an open problem.
The segmentation pipeline (OxBioPathv1, written in MATLAB) operates in two stages. Nuclear segmentation uses Otsu's method to initialize a level set algorithm, followed by a hybrid geodesic region-based level set to estimate nuclear location and shape, then a watershed tesselation to separate touching or overlapping nuclei. Cytoplasmic segmentation begins with a Voronoi tesselation for equidistant partitioning between neighboring cells, followed by iterative marker-controlled watershed using the intensity gradient of the CD99 cytoplasmic marker to define cell and cytoplasmic boundaries.
Segmentation accuracy was validated against manual delineations by three independent experts using Bland-Altman plots (comparing segmented nucleus areas to expert areas) and Hausdorff distance (measuring the maximum boundary discrepancy). Validation was performed in three cell lines (A673, SK-N-MC, RD-ES) on coverslips and in agar pellet cell cores, and then in actual TMA tissue sections. The resulting cell count was 50-500 cells per image depending on cell density, with a total of 113,201 individual cells segmented across all 476 images from 118 patients, yielding up to 4 biomarker measurements per cell and 104 features per patient.
Patient-level features were derived from single-cell measurements by kernel density estimation to approximate the probability density function (PDF) of each feature across all cells from all images of that patient. Each PDF was evaluated at 100 equally spaced points, converting heterogeneous single-cell data into a patient-level distribution representation. This approach captures the full shape of the biomarker distribution, including bimodality and skewness, rather than collapsing it to a single mean or median.
Random Survival Forest (RSF) is a supervised machine learning variant designed specifically for time-to-event (survival) data. It builds an ensemble of 1,000 decision trees using bootstrap sampling (bagging), with each split testing a random subset of N/3 features (default). Out-of-bag samples from each tree provide unbiased error estimates using Harrell's concordance index (C-index), where 0 indicates perfect prediction and 0.5 is expected by chance. Unlike Cox proportional hazards models, RSF requires no threshold specification and handles non-linear feature interactions.
Nine separate RSFs were trained on different combinations of biomarker distributions: DAPI alone, CD99-negative cells DAPI, CD99-positive cells DAPI, Ki67 alone, CD99-negative Ki67, CD99-positive Ki67, combined CD99/DAPI/Ki67, combined CD99/DAPI/Ki67/EGR1/Foxo3a/pS6, and combined CD99/DAPI/Ki67/EGR1/Foxo3a/pS6/pMAPK. The number of patients per analysis ranged from 50 (with pMAPK, requiring complete 7-marker staining) to 118 (for analyses requiring only DAPI, CD99, and Ki67). A variable hunting algorithm with 80/20 train-test splits across 100 independent iterations was used for feature selection.
Internal cross-validation was performed with 50 repetitions of 2/3 train and 1/3 test partitioning to characterize the distribution of error rates across resamples rather than relying on a single estimate. This approach is particularly important given the limited sample size of 118 cases, where a single train-test split could give misleadingly good or poor results by chance. The 50-repeat cross-validation distribution provides a realistic view of expected performance under varying data compositions.
The lowest cross-validation error rates (below 0.40) were achieved by RSFs incorporating Ki67, Ki67 combined with DAPI, and Ki67 combined with EGR1/Foxo3a/pS6. The CD99-negative Ki67 RSF performed among the best, with variable hunting cross-validation mean error of 0.38 (standard deviation 0.064), compared to the all-features baseline error of 0.36 (standard deviation 0.085) for the same marker combination. The full 7-marker RSF (including pMAPK) had worse mean error (0.48) due to its smaller sample size (n=50) limiting statistical power.
The most frequently selected feature in the CD99-negative Ki67 RSF was the 'CD99-negative Ki67 mean nuclear/cytoplasmic ratio' at the 63rd centile, selected 73 times out of 100 variable hunting iterations. This feature identifies a subpopulation of cells that are Ki67-positive (actively proliferating) but have relatively low cytoplasmic CD99 labeling (low nuclear/cytoplasmic CD99 ratio below 1). This CD99-low, Ki67-high subpopulation is prognostically informative and would have been invisible to single-biomarker analysis.
Conventional threshold-based analysis of Ki67 on the same data showed that systematic threshold optimization and independent observer scoring produced qualitatively similar Kaplan-Meier curves, with optimal threshold giving Cox regression p = 1.6 x 10^-6 (log-rank p = 2.6 x 10^-7) and observer scoring giving p = 2.5 x 10^-4 (log-rank p = 1 x 10^-4). While both approaches identified Ki67 as prognostic, they are subject to threshold selection bias and the RSF framework avoids this by analyzing the full biomarker distribution simultaneously across multiple features.
The prognostically informative feature revealed by RSF represents a specific subpopulation of cells that are Ki67-positive but CD99-low, defined as cells where the mean nuclear CD99 exceeds mean cytoplasmic CD99 (nuclear/cytoplasmic ratio greater than 1 for Ki67 but less than 1 for CD99). This subpopulation would be systematically missed or diluted by standard Ki67 index calculations that consider all CD99-positive cells. The biological identity of this subpopulation is unknown and may represent an undifferentiated progenitor or cancer stem cell population.
The hypothesis of a cancer stem cell connection is supported by existing literature: CD133-positive Ewing sarcoma cells have been isolated with stem cell properties, capable of sustaining tumor growth through serial transplantation and differentiating into adipocytes and osteogenic lineages. Additionally, EWS-FLI1 regulates mir145 and SOX2, and TARBP2-dependent miRNA maturation controls cancer stem cell identity in ES. Whether the CD99-low proliferative cells identified here correspond to CD133-positive stem cells requires further experimental investigation.
Individual patient survival curves were predicted by the RSF classifier, enabling a personalized survival prediction approach rather than grouping patients into two binary risk categories. Cross-validation showed that when the error rate was low (near 0.22), the predicted high and low mortality groups showed clearly divergent survival curves. When the error rate was high (near 0.53), separation was poor, illustrating the inherent variability from small sample sizes and highlighting the need for larger prospective cohorts.
The most critical limitation is the small effective sample size of 118 informative cases after quality control from an initial 524. The high attrition rate (36% discarded for poor image quality, plus 52 and 77 cases lost earlier) means that the majority of available material was unusable. The authors explicitly note that a conventional separate training and validation approach is not feasible at this sample size, even with the combined three-cohort dataset. Cross-validation is used as a substitute but introduces uncertainty that a separate validation cohort would resolve.
Pre-analytical processing variability is identified as a major cause of the 36% image failure rate. Samples from three different European pathology laboratories were collected and stored under local guidelines with no standardized ischemic time or fixation conditions. The authors demonstrate that automated QC classification of good vs poor images achieved only AUC 0.72, insufficient for reliable automated exclusion. Addressing pre-analytical variability through prospective collection following standardized protocols (e.g., BRISQ guidelines) is identified as essential for future studies.
The three cohorts differed significantly in survival outcomes (Cohort a: 70% alive, Cohort b: only 25% alive, Cohort c: 54% alive), tumor site distribution, and metastasis rates. These clinical differences could introduce confounding in the combined RSF analysis. PCA of the image features showed cohort-specific clustering, but TMA cores from the same patient imaged separately clustered together, demonstrating intra-patient consistency. The authors chose combined analysis because separate cohort analyses would have insufficient power.
The authors are applying this methodology in the LINES trial (EuroSarc), a Ewing sarcoma exploratory biomarker Phase II study in patients with advanced disease testing linsitinib, a dual IGF1R/IR-A kinase inhibitor. The RSF pipeline is being used to identify biomarker distributions in core biopsies before and after linsitinib exposure, with the goal of identifying cell populations resistant to target inhibition. This represents the translational application of the discovery methodology to drug development.
Future prospective validation requires standardized pre-analytical tissue processing as a prerequisite. The authors recommend mandating rapid tissue ischemic time and combined paraffin-coated nitrogen storage to maximize biomarker preservation. Without standardized sample handling, even well-designed computational pipelines cannot overcome the signal-to-noise limitations introduced by poor tissue preservation, as evidenced by the 36% exclusion rate in this study.
Emerging dimensionality reduction strategies and genomic enrichment classification of cancer subtypes will likely reshape biomarker study design, reducing reliance on rigid training-validation cohort separation as established in REMARK guidelines. Integration of the image-based heterogeneity quantification with genomic profiling data (e.g., EWSR1 translocation type, 1q gain, CDT2 overexpression) could generate multi-modal prognostic classifiers combining structural tumor biology with molecular genotype, which is the likely direction of next-generation biomarker research in Ewing sarcoma.