Ewing sarcoma (ES) is the second most common malignant bone tumor in children and adolescents, yet its clinical behavior remains highly variable and difficult to predict. Patients with the same apparent stage and histopathological features can have drastically different outcomes under identical treatment regimens combining surgical resection, chemotherapy, and radiotherapy. This heterogeneity drives a clinical need for reliable biomarkers that can identify high-risk patients early and guide individualized treatment decisions before disease progression becomes apparent.
The immune microenvironment in ES: Evidence has accumulated that immune cell infiltration plays a central role in shaping ES outcomes. Flow cytometry studies have documented elevated CTLA-4+ T-cell infiltration and CD14+ HLA-DRlo/neg monocytes relative to CD4+ T cells. NK cells and T lymphocytes exert direct cytotoxic effects on ES tumor cells through NKG2D and DNAM-1 receptor-mediated killing. Lymphocyte counts serve as a prognostic indicator for recovery following chemotherapy, and tumor-associated macrophage (TAM) infiltration, particularly CD68+ TAM, correlates with poor overall survival. Despite these individual findings, no unified immunogenomic prognostic model had been developed for ES prior to this work.
Limitations of single-biomarker approaches: Prior ES prognostic efforts focused on individual molecular markers, including SOX2, MTAP, Connexin 43, and Nucleophosmin. Known inflammatory blood biomarkers such as CRP (AUC = 0.676), the Glasgow Prognostic Score (AUC = 0.634), neutrophil-to-lymphocyte ratio (AUC = 0.618), and circulating cell-free mitochondrial DNA (AUC = 0.708) provide only modest predictive power. Single markers inevitably omit prognostic information; conversely, using too many independently selected markers introduces redundancy. A multi-gene machine learning approach addresses both problems simultaneously.
This study from researchers at Xining First People's Hospital and Lanzhou University constructed and validated an 11-gene immune-related gene (IRG) signature for ES prognosis by applying iterative Lasso Cox regression to gene expression data, then characterized the relationship between those genes and immune cell infiltration patterns across the tumor microenvironment.
The study used two independent datasets drawn from major public genomic repositories. The training set consisted of dataset GSE17679 from the Gene Expression Omnibus (GEO) database, comprising 88 ES tumor samples and 18 matched normal skeletal muscle samples profiled on the Affymetrix GPL570 HG-U133_Plus_2 array. The external validation set came from the International Cancer Genome Consortium (ICGC) database and contained 57 ES samples with RNA-seq data and paired clinical outcome information. This two-cohort design, using entirely separate samples for training and independent validation, is critical for assessing whether a prognostic signature will generalize beyond the institution where it was developed.
Patient characteristics: In the GEO training cohort (n = 88), patients were split evenly between those aged 16 and under (50%) and those over 16 (50%), with 63.6% male. The vital status showed 42% alive and 58% deceased at follow-up, reflecting the aggressive nature of the disease in this cohort. In the ICGC validation cohort (n = 57), 63.2% of patients were age 16 or under and 54.4% were male. Importantly, the two cohorts were generated on different technology platforms (microarray versus RNA-seq), making successful cross-platform validation a particularly stringent test of the signature's robustness.
Gene expression preprocessing: The GEO microarray data were subjected to background correction and normalization using the R affy package (version 3.6.3). Quality was verified visually via box plots and principal component analysis (PCA) before and after normalization, confirming that pre-processing produced consistent mean expression values across samples and a reliable data source for downstream analysis. Differential expression analysis used the limma package with strict criteria: absolute log2 fold change greater than 1 and adjusted p-value less than 0.05.
Immune-related gene list: The curated IRG list was sourced from the ImmPort database, a publicly available immunology repository maintained by NIAID. DEGs meeting the fold-change and significance thresholds were intersected with this IRG list to isolate differentially expressed immune-related genes (DEIRGs) as the candidate feature pool for model construction.
Applying differential expression analysis to GSE17679 identified a total of 249 differentially expressed immune-related genes (DEIRGs), visualized in a volcano plot where red nodes represent upregulated genes and green nodes represent downregulated genes relative to normal skeletal muscle. This pool of 249 candidates provided the biological starting material from which the eventual prognostic signature would be distilled.
Gene Ontology (GO) enrichment analysis: GO analysis of the 249 DEIRGs using the R clusterProfiler package revealed biologically coherent enrichment patterns across all three GO domains. In biological processes (BP), the most significant changes involved leukocyte migration, leukocyte chemotaxis, cell chemotaxis, response to steroid hormones, peptidyl-tyrosine phosphorylation, and regulation of chemotaxis, consistent with active immune remodeling of the tumor microenvironment. Cellular component (CC) enrichments were concentrated on the external side of the plasma membrane, vesicle lumen, secretory granule lumen, collagen-containing extracellular matrix, and membrane rafts, pointing to secretory and cell-surface immune activity. Molecular function (MF) enrichments included cytokine activity, growth factor activity, cytokine receptor binding, steroid hormone receptor activity, and transcription factor activity with direct ligand regulation.
Protein-protein interaction (PPI) network: STRING (version 11.0) was used to construct a PPI interaction network for all 249 DEIRGs, which was then imported into Cytoscape for visualization. The cytoHubba plugin, using the MCC (Matthews Correlation Coefficient) algorithm, identified the ten hub genes with the highest network connectivity scores: VEGFA, MMP9, CXCL8, EGF, IGF1, CXCR4, TGFB1, EGFR, SPP1, and ICAM1. These hubs represent the most densely interconnected nodes in the immune-related interaction network, and several have established roles in ES biology. VEGFA and IGF1 promote tumor angiogenesis; CXCR4 and MMP9 drive ES invasion and metastasis; EGF-EGFR signaling stimulates ES cell proliferation through EWS-FLI1 and neurotrophin pathways.
Functional similarity analysis: The R GOSemSim package was used to compute semantic similarity between the GO annotations of hub genes. This revealed strong biological functional similarity between EGFR and SPP1, and between CXCR4 and MMP9, suggesting that SPP1 may share ES cell proliferation functions with EGFR, and that CXCL8, TGFB1, SPP1, and ICAM1 may represent previously unexplored ES immunotherapy targets.
Before applying Lasso regression, the authors performed univariate Cox regression on all 249 DEIRGs and retained only those with p < 0.05 as candidate genes. This step reduced the feature space to a manageable set of prognostically relevant IRGs, filtering out genes with no detectable association with relapse-free survival (RFS) in the training cohort. The use of Lasso (Least Absolute Shrinkage and Selection Operator) regression, which applies an L1 norm penalty that shrinks many regression coefficients to exactly zero, is well-suited for high-dimensional data where many candidate features are correlated with each other.
Iterative approach with 1000 Lasso models: Rather than relying on a single Lasso fit, the authors ran 1000 independent Lasso Cox regression iterations, each producing a different gene combination. Gene combinations appearing in more than 50 of the 1000 iterations were retained as robust candidates. The final signature was selected by identifying the combination with the largest area under the ROC curve and the strongest overall correlation with prognosis. This iterative cross-validation approach reduces the risk of selecting an unstable gene set that performs well by chance in one particular split of the data.
The 11-gene immune signature: The resulting optimal prognostic signature consists of 11 immune-related genes: CRLF3 (cytokine receptor-like factor 3), ECD (ecdysoneless cell cycle regulator), FABP4 (fatty acid binding protein 4), FGF6 (fibroblast growth factor 6), GNRH2 (gonadotropin-releasing hormone 2), NDRG1 (N-myc downstream regulated gene 1), PAK2 (p21-activated kinase 2), PLTP (phospholipid transfer protein), PTGDS (prostaglandin D2 synthase), RBP1 (retinol binding protein 1), and ZC3HAV1 (zinc finger CCCH-type containing, antiviral 1). These genes span diverse biological functions including lipid metabolism, cytokine signaling, antiviral immune response, and cell cycle regulation, reflecting the multidimensional nature of immune-tumor interactions in ES.
Risk scoring: A patient-level risk score was computed by combining the expression levels of the 11 genes weighted by their respective Cox regression coefficients. Patients were divided into high-risk and low-risk groups based on the median risk score in the training set. Kaplan-Meier survival curves showed a highly significant separation between groups (log-rank p < 0.001), with high-risk patients experiencing substantially shorter relapse-free survival.
The 11-gene signature was rigorously tested in the independent ICGC external cohort of 57 ES patients profiled by RNA-seq, a fundamentally different measurement technology than the Affymetrix microarrays used for training. Cross-platform validation is one of the most stringent benchmarks for a genomic prognostic tool because systematic technical differences between microarray and RNA-seq platforms can confound gene expression values and artificially inflate or deflate model performance.
Time-point AUC analysis: The authors constructed ROC curves at 1-year, 3-year, and 5-year endpoints using the external ICGC cohort. The 1-year AUC was 1.00, indicating perfect discrimination between patients who experienced events within one year and those who did not. The 3-year AUC was 0.92 and the 5-year AUC was 0.97, both substantially above the performance of any previously reported individual ES biomarker. These results confirm that the signature is not an artifact of the training data or the microarray platform.
Comparison with established ES biomarkers: The authors performed head-to-head ROC comparisons between the 11-gene signature and four known prognostic markers in the external cohort: SOX2, STAG2, TP53, and CXCR4. The 11-gene signature outperformed all four at every time point tested. In contrast, the inflammatory blood biomarkers reported by Li et al. (CRP AUC = 0.676, GPS AUC = 0.634, NLR AUC = 0.618, PLR AUC = 0.613, LMR AUC = 0.574, NPC AUC = 0.571) and cell-free mitochondrial DNA (AUC = 0.708) all fall far below the 5-year AUC of 0.97 achieved by the 11-gene signature. These comparisons position the signature as the most accurate prognostic tool reported for ES at the time of publication.
Survival stratification in the external cohort: Kaplan-Meier analysis of the ICGC cohort using the risk score derived from the training set showed a highly significant survival difference between high-risk and low-risk groups (p < 0.0001). Gene expression patterns in the external cohort also separated clearly between the two risk groups when visualized using Cox regression-based risk plots, further confirming the biological consistency of the signature across independent datasets and platforms.
Demonstrating prognostic value in Kaplan-Meier analysis is necessary but not sufficient to establish clinical utility. A signature that correlates with prognosis only because it tracks a known clinical variable (such as metastatic status) adds no independent information beyond what clinicians already assess. The authors therefore performed both univariate and multivariate Cox proportional hazard analyses, incorporating clinicopathological characteristics alongside the 11-gene risk score to determine whether the signature retains prognostic value after adjusting for confounders.
Univariate Cox results: In the training set, metastatic status showed a borderline significant association with prognosis (HR = 1.479, p = 0.016), while gender (HR = 0.902, p = 0.724) and age (HR = 0.892, p = 0.892) did not reach significance. The 11-gene risk score was strongly significant in univariate analysis in both the training set (HR = 1.185, p < 0.001) and the external validation set (HR = 2.105, p < 0.001). The hazard ratio of 2.105 in the external cohort indicates that each unit increase in risk score was associated with more than doubling of the hazard of relapse or death.
Multivariate Cox results: After adjusting simultaneously for gender, age, and metastatic status, the 11-gene signature retained strong independent prognostic significance in both cohorts: training set HR = 1.185 (p < 0.001) and external validation set HR = 2.098 (p = 0.003). Notably, metastatic status lost significance in multivariate analysis (HR = 1.356, p = 0.291 in training; HR = 1.477, p = 0.346 in validation), suggesting that the gene signature may partially capture the biological mechanisms underlying metastatic risk. Gender and age remained non-significant in all multivariate models.
The consistency of the independent prognostic effect across two datasets, two measurement platforms, and two geographical patient populations (Chinese hospital cohort vs. ICGC international consortium samples) substantially strengthens confidence that the 11-gene signature reflects genuine tumor biology rather than dataset-specific noise. This independence from standard clinical variables is a prerequisite for the signature to add value in real-world clinical decision-making.
Beyond constructing the prognostic gene signature, the authors comprehensively characterized the immune landscape of ES using single-sample gene set enrichment analysis (ssGSEA). The ssGSEA algorithm, implemented in the R GSVA package, quantifies the degree to which a predefined set of genes (corresponding to each immune cell type) is coordinately up- or down-regulated within a single sample. This allows immune cell infiltration levels to be estimated from bulk gene expression data without requiring histological staining or flow cytometry. The analysis covered 24 distinct immune cell populations defined by the Bindea et al. (2013) immune cell marker gene sets.
PCA-based immune cell clustering: Principal component analysis of the ssGSEA-estimated immune cell infiltration scores revealed a clear separation between ES tumor samples and normal skeletal muscle controls, confirming that the tumor microenvironment is immune-distinct from normal tissue and that the ssGSEA approach captures biologically meaningful differences.
Immune cell correlation and interaction networks: Correlation analysis across all 24 immune cell types revealed significant positive correlation between B cells and cytotoxic T cells, and negative correlations involving CD8 T cells with NK cells, NK cells with NK CD56bright cells, and NK CD56bright cells with Th2 cells. An interaction network built using the igraph package (threshold: p < 0.05 and |correlation coefficient| > 0.15) identified NK cells and NK CD56bright cells as having the strongest interactions with other immune populations, while Th2 and effector memory T cells (TEM) showed the weakest interactions.
Prognostically significant immune cell types: Kaplan-Meier analysis for each of the 24 cell types identified ten with significant survival associations. Infiltration of cytotoxic T cells (HR = 0.45, 95% CI: 0.22-0.89, p = 0.004), macrophages (HR = 0.43, 95% CI: 0.24-0.76, p = 0.001), mast cells (HR = 0.49, 95% CI: 0.28-0.84, p = 0.01), central memory T cells (TCM; HR = 0.38, 95% CI: 0.22-0.65, p < 0.001), and TEM cells (HR = 0.45, 95% CI: 0.24-0.85, p = 0.003) were associated with favorable prognosis. In contrast, higher infiltration of B cells (HR = 2.27, p = 0.003), CD8 T cells (HR = 2.15, p = 0.01), NK CD56bright cells (HR = 2.02, p = 0.01), NK cells (HR = 2.33, p = 0.001), and Th2 cells (HR = 3.67, 95% CI: 2.03-6.63, p < 0.001) predicted poor prognosis.
Signature genes and immune cell correlations: Correlating each of the 11 signature genes against the 24 immune cell infiltration scores revealed that FABP4 (fatty acid binding protein 4) had the strongest positive correlation with macrophages, while NDRG1 (N-myc downstream regulated gene 1) had the strongest negative correlation with Th2 cells. FABP4 is involved in lipid metabolism and macrophage function, and high macrophage infiltration is associated with favorable ES prognosis (HR = 0.43). The negative NDRG1-Th2 relationship is notable because Th2 cell infiltration carries the worst prognostic hazard ratio in the analysis (HR = 3.67). These gene-immune cell pairings suggest mechanistic hypotheses for how the signature genes influence tumor immunity, and provide candidate entry points for future experimental validation.
Dataset size and clinical data constraints: The training cohort of 88 samples and external validation cohort of 57 samples are small by contemporary genomic standards, particularly for a multi-gene model with 11 features. While the iterative Lasso approach with 1000 iterations and frequency filtering mitigates overfitting risk, the relatively small sample sizes mean that confidence intervals around hazard ratios are wide and precise coefficient estimation remains uncertain. Clinical annotation in the available datasets was limited to gender, age, and metastatic status, precluding adjustment for treatment regimen, institutional practice differences, or molecular subtype, all of which are known to influence ES outcomes.
Technology and functional validation gaps: The study is entirely computational, built on existing public gene expression datasets. The 11-gene signature was derived from Affymetrix microarray data and validated on RNA-seq, but has not been translated into a clinical assay format such as RT-qPCR panels or NanoString CodeSets. The functional roles of the signature genes in ES biology remain largely uncharacterized at the wet-lab level. Genes such as GNRH2, PLTP, PTGDS, and ZC3HAV1 have no prior reported connection to ES, and whether their differential expression reflects a causal role in immune evasion or tumor progression versus a correlative association requires experimental follow-up.
Prospective validation priority: The authors note that large-scale prospective clinical validation with independent, multi-institutional ES samples is essential before the 11-gene signature can be applied in clinical practice. This would ideally include formalin-fixed paraffin-embedded (FFPE) tumor samples profiled on a clinically deployable platform, with outcome data extending to both overall survival and event-free survival. Integration of the signature with established clinicopathological risk factors and emerging molecular classifications of ES would also improve its translational utility.