$$\rightleftharpoonup{xx}$$
$$\longleftharp{xx}$$,
$$\longrightharp{xx}$$,
All procedures involving human tissues complied with institutional guidelines and the Declaration of Helsinki and were approved by the Institutional Review Board of Fujian Medical University (Approval No. 2021KYB089). Written informed consent was obtained from all participants prior to tissue procurement.
Gene expression and survival analysis
RNA sequencing data and corresponding clinical information were obtained from multiple public databases. 1) TCGA cohort: RNA-seq (FPKM) data for 175 glioblastoma multiforme (GBM) and 534 low-grade glioma (LGG) samples were downloaded from The Cancer Genome Atlas (https://portal.gdc.cancer.gov/); 2) Normal controls: Expression profiles of 211 normal brain tissues and 662 glioma tissues were downloaded from the UCSC Xena database (https://xenabrowser.net/datapages/); 3) External validation: Data from CGGA693 and CGGA325 cohorts were obtained from the Chinese Glioma Genome Atlas (http://www.cgga.org.cn); 4) GEO dataset: The GSE43378 dataset, containing expression and clinical data for 50 glioma samples, was downloaded from the Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/). All raw count data were converted to transcripts per million (TPM) and log2-transformed. For datasets already normalized, expression matrices were examined to ensure comparable distributions. Genes with TPM values < 1 in more than 80% of samples were excluded. Missing clinical information (age, IDH status, 1p/19q codeletion, MGMT methylation) was removed using complete-case filtering. Batch effects among datasets were adjusted using the ComBat algorithm implemented in the R package sva. Expression values were standardized by z-score transformation within each dataset. Survival analyses were performed using the R packages survival and survminer. Patients were dichotomized into high- and low-expression groups according to the median expression level of IRAIN. Kaplan-Meier survival curves were generated, and statistical significance was evaluated by the log-rank test. Hazard ratios (HRs) and 95% confidence intervals (CIs) were estimated using Cox proportional hazards regression models.
Definition of immune and metabolic gene sets
Immune-related genes (IRGs, n = 2,483) were obtained from the ImmPort database (https://www.immport.org/shared/), and metabolic-related genes (MRGs, n = 948) were obtained from the Molecular Signatures Database (MSigDB, https://www.gsea-msigdb.org/). The combined set of these genes was defined as immunometabolic-related genes (IMRGs). These gene lists served as references for subsequent differential expression and network analyses.
Differential expression and weighted gene co-expression network analysis
Differentially expressed genes (DEGs) between normal brain and glioma tissues were identified using the R package limma. Expression data were fitted with a linear model followed by empirical Bayes moderation. Genes with |log₂ fold change| > 1.5 and false discovery rate (FDR) < 0.05 were considered significantly differentially expressed. Weighted gene co-expression network analysis (WGCNA) was conducted using the R package WGCNA. Outlier samples were excluded through hierarchical clustering. The soft-thresholding power was set to β = 8 to achieve a scale-free topology fit index (R2 ≥ 0.85) while maintaining adequate mean connectivity. Topological overlap matrices (TOM) were constructed, and genes were grouped into modules with a minimum size of 50 using the dynamic tree-cut algorithm. Module eigengenes were correlated with clinical traits, and the module most strongly associated with glioma (Pearson's r > 0.7, P < 1×10-10) was selected for hub gene identification.
Machine learning-based prognostic model construction
A comprehensive leave-one-out cross-validation (LOOCV) framework integrating ten machine-learning algorithms was applied to construct and evaluate prognostic models. In total, 101 combinatorial workflows were implemented using the TCGA cohort as the training dataset. Prognosis-associated immunometabolic-related genes (IMRGs) were first identified by univariate Cox regression (P < 0.05). The optimal model was determined by maximizing the mean Harrell's concordance index (C-index) across three validation datasets (CGGA693, CGGA325, and GSE43378). The resulting RSF-Enet model (α = 0.3) demonstrated the highest predictive performance and maintained robust generalizability across independent cohorts22.
TME and immune infiltration
To comprehensively characterize the immunogenomic landscape, we employed a multi-tiered analytical approach. First, immune and stromal infiltration levels were quantified using the ESTIMATE algorithm23. Differential expression of key immune checkpoint molecules, including PDCD1, CTLA4, and LAG3, was then assessed through limma-based analysis, and the correlations among checkpoint genes were visualized using correlation matrices. Somatic mutation profiles from 903 glioma samples in the TCGA cohort were used to calculate tumor mutational burden (TMB), microsatellite instability (MSI), and tumor immune dysfunction and exclusion (TIDE) scores to predict potential responses to immunotherapy. Patients were subsequently stratified into four prognostic groups according to combined TMB status (high/low) and risk scores (high/low), and survival outcomes were compared using Kaplan-Meier analysis.
Functional enrichment analysis
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were conducted using the R package clusterProfiler. Enrichment results with adjusted P values < 0.05 were considered statistically significant. Overrepresented biological processes, cellular components, and molecular functions were visualized using dot plots and bar plots. Protein-protein interaction (PPI) networks were constructed using the STRING database (≥ 0.4) and visualized in Cytoscape. Functional modules within the PPI network were identified using the MCODE algorithm. Gene-gene interaction and co-expression networks were further analyzed using GeneMANIA (https://string-db.org; confidence score ≥ 0.4) and visualized in Cytoscape. Functional modules within the PPI network were identified using the MCODE algorithm. Gene–gene interaction and co-expression networks were further analyzed using GeneMANIA (https://genemania.org), which integrates information on physical and genetic interactions, shared pathways, and co-expression patterns to infer potential functional associations.
Clinical specimens
Fresh glioma tissues (n = 6) and paired adjacent non-tumorous brain tissues (n = 6; located at least 3 cm from the tumor margin and histologically confirmed as tumor-free) were collected from patients undergoing primary glioma resection at the Zhangzhou Affiliated Hospital of Fujian Medical University. None of the patients had received chemotherapy or radiotherapy prior to surgery. All pathological diagnoses were independently verified by two neuropathologists according to the 2021 World Health Organization (WHO) classification of central nervous system tumors. Immediately after surgical excision, tissue specimens were rinsed with ice-cold phosphate-buffered saline (PBS) to remove residual blood, snap-frozen in liquid nitrogen (-196 °C), and stored at -80 °C until RNA extraction.
Cell lines and cell culture
Human glioblastoma cell lines SHG44, U251, A172, and T98G, as well as normal human glial cells (HEB), were obtained from authenticated repositories and confirmed to be free of mycoplasma contamination prior to use. Cells were maintained in Dulbecco's Modified Eagle's Medium (DMEM, high glucose) supplemented with 10% fetal bovine serum (FBS), 2 mM L-glutamine, and 1% penicillin-streptomycin, at 37 °C in a humidified incubator with 5% CO₂. Cells were passaged every 4-5 days upon reaching 80-90% confluence. To establish IRAIN-overexpressing and control cell lines, cells were transduced with lentiviral vectors carrying the full-length IRAIN transcript or an empty vector as control. Stable clones were selected using puromycin (2 µg/mL) for 14 days. Overexpression efficiency was confirmed by quantitative reverse transcription PCR (qRT-PCR) prior to downstream assays.
3- (4,5-dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide (MTT) cell proliferation assay
Cells were seeded in 96-well plates at a density of 1 × 104 cells per well in 100 µL of complete culture medium. At 24, 48, and 72 h after seeding, 20 µL of MTT solution (5 mg/mL in phosphate-buffered saline) was added to each well and incubated for 4 h at 37 °C. The supernatant was then removed, and 150 µL of dimethyl sulfoxide (DMSO) was added to dissolve the formazan crystals. The plate was gently agitated for 10 min to ensure complete solubilization. Absorbance was measured at 490 nm using a microplate spectrophotometer. Background readings from blank wells were subtracted. Cell viability was calculated relative to the 24-h or control group (set as 1.0). All experiments were performed with six technical replicates and three independent biological replicates. Data are expressed as mean ± standard deviation (SD), and statistical significance was determined using a two-tailed t-test.
Flow cytometry for apoptosis (Annexin V - FITC/PI staining)
Cells were seeded at 60-70% confluence and treated for 24 h under the indicated conditions. Floating and adherent cells were collected using EDTA-free trypsin, combined, and washed twice with ice-cold PBS. Cell pellets were resuspended in Annexin V binding buffer (10 mM HEPES pH 7.4, 140 mM NaCl, 2.5 mM CaCl2) at 1 × 106 cells/mL. For each sample, 100 µL of suspension was incubated with 5 µL of Annexin V-FITC and 5 µL of propidium iodide (PI; 50 µg/mL stock) in the dark for 15 min at room temperature. Following the addition of 400 µL of binding buffer, samples were kept on ice and analyzed within 1 h on a flow cytometer (488 nm excitation; 530/30 nm for FITC and >585 nm for PI). Appropriate single-stain and fluorescence-minus-one controls were included for compensation. At least 10,000 events were recorded per sample. Data were analyzed by quadrant gating: live (Annexin V⁻/PI⁻), early apoptotic (Annexin V⁺/PI⁻), late apoptotic (Annexin V⁺/PI⁺), and necrotic (Annexin V⁻/PI⁺) populations. Percentages of early + late apoptotic cells were reported (mean ± SD, n = 3).
Quantitative real-time PCR (qRT-PCR)
Total RNA was isolated using an acid phenol-guanidinium reagent according to the manufacturer's protocol. RNA purity was verified by spectrophotometry (A₂₆₀/A₂₈₀ = 1.8-2.1), and integrity was confirmed by gel electrophoresis (RNA integrity number ≥ 7). One microgram of total RNA was treated with DNase I and reverse-transcribed in a 20 µL reaction using random hexamers and oligo(dT) primers. The reaction was carried out at 25 °C for 10 min, 50 °C for 30 min, and 85 °C for 5 min. Quantitative PCR was performed in a 10 µL system containing 5 µL of 2× SYBR Green Master Mix, 0.3 µM each primer, and 1 µL of cDNA (≈ 20 ng RNA equivalent). Thermal cycling conditions were 95 °C for 5 min, followed by 40 cycles of 95 °C for 15 s and 60 °C for 30 s, then a melt-curve analysis from 65 °C to 95 °C in 0.3 °C increments. All reactions were performed in triplicate, together with no-template and minus-RT controls. Ct values > 35 or technical replicate SD > 0.5 were excluded. Relative expression was calculated using the 2⁻ΔΔCt method, with GAPDH as the internal control. Mean ± SD values from three independent biological replicates were reported, and group differences were analyzed using a two-tailed t-test.
Western blot analysis
Cells were lysed on ice in RIPA buffer (50 mM Tris-HCl, pH 7.4, 150 mM NaCl, 1% NP-40, 0.5% sodium deoxycholate, 0.1% SDS) supplemented with protease and phosphatase inhibitors. Lysates were incubated for 30 min on ice with intermittent vortexing and cleared by centrifugation at 12,000 × g for 15 min at 4 °C. Protein concentrations were measured by BCA assay, adjusted to 1-2 µg/µL, and mixed 1:3 with 4× Laemmli buffer (final 1× buffer containing 100 mM DTT). Samples were denatured at 95 °C for 5 min. Equal amounts of protein (50 µg) were resolved by 12 % SDS-PAGE at 100 V for 90 min and electro-transferred to PVDF membranes at 250 mA for 90 min. Membranes were blocked with 5 % non-fat milk in TBST (0.1 % Tween-20) for 1 h at room temperature (or 5 % BSA for phosphoproteins) and incubated overnight at 4 °C with primary antibodies against IGF1, IGF1R, JAK2, p-JAK2 (Y1007/1008), STAT3, p-STAT3 (Y705), BIRC5, and β-actin (typical dilution 1:1000, β-actin 1:5000). After three 10-min washes in TBST, membranes were incubated with HRP-conjugated secondary antibody (1:5000) for 1 h at room temperature, washed again, and developed using chemiluminescent substrate. Band intensities were quantified with ImageJ, normalized to β-actin or total protein, and expressed as mean ± SD from three independent experiments.
Immunocytochemistry
Cells grown on sterile glass coverslips were rinsed twice with PBS and fixed in 4% paraformaldehyde for 15 min at room temperature. After three PBS washes, cells were permeabilized with 0.2% Triton X-100 for 10 min, blocked with 5% bovine serum albumin (BSA) for 1 h, and incubated overnight at 4 °C with primary anti-CD31 antibody (1:200 dilution in 1% BSA). Following three PBS washes, cells were incubated with Alexa Fluor-conjugated secondary antibody (1:500 dilution) for 1 h in the dark, counterstained with DAPI (1 µg/mL, 5 min), and mounted in antifade medium. Images were captured using a fluorescence microscope under identical exposure and gain settings. The percentage of CD31-positive area was quantified in five randomly selected non-overlapping fields per sample using ImageJ software. This assay was performed in cell models rather than tissue sections.
Statistical analysis
Statistical analyses were performed utilizing R version 4.3.0 along with its associated packages. To compare categorical variables, the chi-squared test was employed, while continuous variables were assessed using either the Wilcoxon rank-sum test or the T test. The evaluation of continuous variables was accomplished through Pearson's correlation coefficient. Survival analyses were carried out using the survival package, which included Cox proportional hazards modeling and the generation of Kaplan-Meier curves, with optimal stratification thresholds established by the survminer package and the formula Riskscore = 
. The CompareC package was used to assess the C-indices of various variables. The receiver operating characteristic curve (ROC), aimed at predicting binary categorical variables, was generated using the pROC package. Additionally, the time-dependent area under the ROC curve (AUC) for survival metrics was analyzed using the timeROC package. All statistical tests were conducted with a two-sided approach. A significance level of P < 0.05 was considered statistically significant.