Single-cell RNA sequencing (scRNA-seq) data processing was performed using Seurat (v.5.2.1) (14) in R (v.4.4.0). The dataset comprised 19 samples (12 post-operative, 7 pre-operative biopsies). Low-quality cells were filtered based on the following criteria: nFeature_RNA > 300, nCount_RNA > 800, and percentage of mitochondrial genes (percent.mt) < 15. Gene expression counts were normalized using a scale factor of 10,000. Highly variable features were identified using the FindVariableFeatures function, selecting the top 5000 features. Data scaling was performed using default parameters. In the batch effect correction process, considering that the sequencing data were derived from two different platforms, we used the sequencing platform (BD/XGY) as the primary covariate for Harmony correction to eliminate systematic differences introduced by the varying technical platforms. We also recognized that other factors such as the source of patient samples and sample preparation dates might introduce variations. In this study, sample information was incorporated as additional correction covariates into the Harmony algorithm. This ensured that the corrected data not only removed technical batch effects but also minimized the potential confounding effects from differences in sample collection sites. Dimensionality reduction was carried out using Principal Component Analysis (PCA), with the first 30 principal components used for downstream analysis. Clustering was performed based on these Harmony-corrected embeddings at a resolution of 0.4. Major cell populations were annotated based on the expression of canonical marker genes: Epithelial cells (EPCAM, KRT19), Endothelial cells (VWF, CDH5), Mast cells (GATA2, KIT), Fibroblasts (PDGFRA, FAP), Pericytes (ACTA2, RGS5), Myeloid cells (CD14, FCGR3A), T cells (CD3D, CD3E), B cells (MS4A1, CD79B), and Plasma cells (MZB1, JCHAIN).

VISIUM Spatial transcriptomics (ST) data were also processed using Seurat (v.5.2.1). Spots with total counts less than 200 and genes detected in fewer than 10 spots were filtered out. Data were read using Seurat.v5, specifying high-resolution images and their corresponding scale factors. Data layers were merged using JoinLayers, and normalization was performed using the SCTransform method. Spatial feature plots were generated to visualize gene expression and cell type distributions. For cell type annotation of spatial spots, deconvolution was performed using spacexr [Robust Cell Type Decomposition (RCTD)] (v.2.2.1) package (15). A reference was built using the raw count matrix from the annotated scRNA-seq data, focusing on the 17 most abundant and representative cell types. The total UMI count per cell (nUMI) was calculated for the reference. For each spatial sample, the raw count matrix was extracted, and nUMI per spot was calculated. Spatial coordinates for each spot were obtained using GetTissueCoordinates. A SpatialRNA object was created for each spatial sample by combining the count matrix, nUMI vector, and coordinate data frame. RCTD was run in “full” mode with default parameters. The resulting cell type weight matrices for each spot were added as an assay to the Seurat object. Cell type localization was visualized using SpatialFeaturePlot. Additionally, the vizAllTopics function from the STdeconvolve(1.3.2) package (16) was used to visualize RCTD results. Spatial co-localization analysis was performed using SpatialFeaturePlot with default parameters.

For subcluster analysis of major lineages, cells from each population were subset, re-normalized, and scaled. Harmony was re-run on each subset to mitigate residual batch effects, using a clustering resolution of 0.4. Myeloid cells were subclassified into Monocytes (FCN1, VCAN), Macrophages (CD68, CD163), Neutrophils (FCGR3B, CXCR2), and Dendritic cells (CD1C, CLEC10A) based on their respective canonical markers. Macrophages were further subdivided into five subsets: CCL4L2⁺TAMs, CYP27A1⁺TAMs, FGL2⁺TAMs, LPL⁺TAMs, and SPP1⁺TAMs. T cells were subclassified into CD8⁺T cells (CD8A, CD8B), CD4⁺ T cells (IL7R, CCR7), T regulatory cells (Treg; FOXP3, IL2RA), Cycling T cells (MKI67, TOP2A), and NK cells (FCGR3A, GNLY). Fibroblasts were subclassified into inflammatory Cancer-Associated Fibroblasts (iCAFs; IL6, CXCL12) and myofibroblastic Cancer-Associated Fibroblasts (myCAFs; ACTA2, POSTN). Clusters exhibiting mixed expression of markers from distinct lineages, potentially representing doublets despite prior doublet removal, were identified and excluded from downstream analysis due to the inability to precisely estimate the doublet rate post-filtering.

To reliably distinguish malignant epithelial cells from normal epithelial cells, copy number variation (CNV) was inferred using the inferCNV package (v.1.22.0). To account for platform-specific batch effects, inferCNV was run individually for each sample. Endothelial cells were used as reference cells. Analysis was performed with default parameters. The output files infercnv.observations and infercnv.references were read, and CNV scores were calculated per cell. The top 5% of cells with the highest CNV aberration scores were selected. A similarity matrix was computed for all cells against this malignant reference set. A threshold of 0.20 was applied, and cells within the top 20% of similarity scores were classified as malignant for each sample.

Spatial niches, termed Composition Clusters (CCs), were identified using Identifying Spatial Co­occurrence in Healthy and InflAmed tissues (ISCHIA) (V1.0.0.0) (17), which detects spatial co-occurrence patterns of cell types. The cell type weight matrix generated by RCTD and the integrated multi-sample Seurat object were used as input. The optimal number of niches was determined by evaluating the Calinski-Harabasz Index, Gap Statistic, and Elbow Method, collectively indicating 9 as the optimal number of clusters for subsequent analysis. For analyzing co-occurrence within specific niches (e.g., CC4 and CC6), a probability threshold (prob.th) of 0.05 was used. The distribution of niches across sample groups was visualized using the dittoSeq package (v.1.18.0).

Pathway activity scores for different spatial niches were inferred using the PROGENy(v.1.28.0) (18) method on the integrated spatial transcriptomics data. Activity scores for hallmark pathways were calculated using the SCT assay, without additional scaling (scale = FALSE), using the top 1000 variable genes for model building. The resulting PROGENy score matrix was centered and scaled using the ScaleData function and visualized via heatmap.

Cell type proportions within each niche were analyzed using two methods. Direct Proportion Calculation: RCTD-derived cell type annotations were aggregated by niche, and the proportion of each cell type relative to the total number of cells in each niche was calculated. Results were visualized using boxplots.Gene Set Enrichment-based Scoring. To mitigate issues with very low cell type proportions (e.g., B cells in CC4/CC6), cell type signature scores were calculated. The top 25 significant marker genes (by logFC) for each cell type, identified using FindAllMarkers, were compiled into gene sets. Enrichment scores for these gene sets per spot were calculated using the escape.matrix function from the ESCAPE package (v.2.2.4) (19). The distribution of these enrichment scores across niches was then visualized using boxplots.

To investigate macrophage state transitions and pseudotemporal ordering, macrophages from pathological complete response (pCR) patients were subset from the Seurat object. Pseudotime analysis was performed using Monocle (v.2.34.0). A CellDataSet object was constructed following the package vignette. Highly variable genes were identified with min_expr = 0.3. Genes were retained for dimensionality reduction if their mean expression was >= 0.5 and their empirical dispersion was at least 1.5 times that expected under a null model. Dimensionality reduction was performed using DDRTree, and cell ordering was established using orderCells with default parameters. For branch analysis, branch1 was selected as the branch point. Gene Ontology Biological Process (GO-BP) enrichment analysis was performed on genes associated with identified states/branches.

Spatial trajectory analysis was performed using SPATA2 (v.3.1.4) (20). The Seurat object was converted to a SPATA2 object using asSPATA2. A spatial trajectory was constructed spanning CC2, CC6, and CC4 using createSpatialTrajectories, with a width of 9.78 units and a length of 137.96 units. Gene expression data and trajectory coordinates for spots along the trajectory were extracted. Hallmark gene sets (obtained from msigdbr v.7.5.1) were used for Gene Set Variation Analysis (GSVA) along the trajectory. Trajectory coordinates were transformed into a one-dimensional pseudotime variable ordered along the trajectory direction. A trajectory heatmap was generated, adapting the plot_trajectory_heatmap function from the referenced methodology [Multimodal decoding of human liver regeneration (21)] and customizing aesthetics due to SPATA2 version differences. Spatially variable genes were identified using runSPARKX and filtered (threshold_pval = 0.05, fdr < 0.5) for visualization.

To characterize functional states of macrophages subpopulations, GSVA(v.1.52.3) (22) was performed. Gene sets related to macrophages function were selected from the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway collection. The average expression matrix for each macrophage subtype was calculated using AverageExpression, and the RNA assay matrix was extracted. Pathway enrichment scores were computed using the gsvaParam method and the gsva function with default parameters. To compare functional states before and after immunotherapy, the add.ident parameter was used to group cells by timepoint before calculating the average expression matrix, using the same gene sets and GSVA parameters.

Cell-cell communication analysis was performed using CellChat(v.2.1.2) (23), which incorporates spatial distance constraints. Analyses were run separately for different niches and overall cell types. RCTD-derived cell type weights and the SCT normalized expression matrix were used. Spatial coordinates for each sample were provided separately. A spatial factors object was created for each sample by calculating the micron-per-pixel ratio (based on each spot has a diameter of 55 micrometers, with a spacing of 10 micrometers between adjacent spots). The interaction.range and contact.range parameters were set to 250 and 100, respectively. The databases “Secreted Signaling”, “ECM-Receptor”, and “Cell-Cell Contact” were used. Communication probabilities were computed, filtering out interactions with fewer than 10 potential links (filterCommunication). Pathway-level communication (computeCommunProbPathway) and aggregated networks (aggregateNet) were calculated. Co-expression of ligand-receptor pairs (e.g., CXCL9-CXCR3) was visualized using SpatialFeaturePlot with a minimum cutoff min.cutoff = 0.05.

A quantitative assessment of the spatial association between CYP27A1⁺tumor-associated macrophages (TAMs) and CD8⁺T cells was conducted. First, based on the deconvolution analysis results, positive spots were defined using a thresholding method (spots with CYP27A1⁺TAMs/CD8⁺T cell expression exceeding the 25th percentile of the corresponding expression in all spots were considered positive spots). The spatial coordinates of positive spots were extracted, and samples with fewer than 40 spots for any cell type were filtered out. Multi-type point pattern data were constructed using the spatstat package (v.3.3-1). Cross-K-function (Kcross) and L-function (Lcross) analyses were employed to examine the association patterns between the two cell types at different spatial distances. To assess statistical significance, Monte Carlo simulations were used to calculate 95% confidence intervals, with the null hypothesis being spatial independence between the two cell types. Effect sizes were quantified by calculating the standardized mean deviation (Cohen’s d) of L(r)-r from 0 and the aggregation index (Aggregation Index, proportion of distances where L(r)-r > 0). Spatial distances were converted to actual micrometer units for analysis based on the original image scale.

Spatial co-localization between cell types was analyzed using MISTY (v.1.14.0) (24). Based on the image scale factor, the view distances were set as follows: juxta view = 200 µm, para view = 3000 µm. The zoi parameter for the para view was set to 200 µm to mask the areas already covered by the intra and juxta views. Prior to analysis, an offset of 10,000 units was added to the X and Y coordinates of each sample to prevent coordinate overlap. Analyses were run with default parameters otherwise. Results were summarized across samples by calculating the median standardized importance for each interaction. Results were visualized via heatmap and community interaction plots (plot_interaction_communities), using a cutoff of 1.5 for importance values in all views for community detection.

RNA-seq data (IlluminaHiSeq, n=576 samples) and corresponding clinical/survival data for the TCGA-LUAD cohort were downloaded from UCSC Xena (http://xena.ucsc.edu/).For CYP27A1 expression analysis, samples were dichotomized into high and low expression groups using an optimal cutoff determined by the survminer package (v.0.5.0). Overall survival (OS) was compared between groups using Kaplan-Meier curves generated with the survival package (v.3.8-3), and the log-rank test was used to calculate p-values. For CYP27A1⁺TAMs signature analysis, the marker genes (by logFC) for the CYP27A1⁺TAMs subset were defined as a gene signature. A single-sample Gene Set Enrichment Analysis (ssGSEA) score was calculated for each sample in the TCGA cohort based on this signature. Samples were split into high and low signature score groups based on the median score. Kaplan-Meier survival analysis was performed as described above.

Signature scores for TAMs-LXRhi, TAMs-LXRlow, and CD8⁺T cells (defined by the top marker genes from scRNA-seq, ranked by logFC) were calculated in the TCGA-LUAD cohort using ssGSEA. Pairwise Pearson correlations and associated p-values between these signature scores were computed using cor.test. scRNA-seq Cohort: The proportions of TAMs-LXRhi and TAMs-LXRlow cells among all macrophages, and the proportion of CD8⁺T cells among all T cells, were calculated for the scRNA-seq cohort. Pairwise Pearson correlations between these proportions were computed using cor.test.

External validation datasets were processed with standardized pipelines: scRNA-seq Cohort GSE131907: Cells were filtered (nFeature_RNA > 300, nCount_RNA > 1000, percent.mt < 15), resulting in 201,003 cells. scRNA-seq Cohort GSE223203: Cells were filtered (nFeature_RNA > 750, nCount_RNA > 1500, percent.mt < 15), resulting in 33,443 cells. Spatial Transcriptomics Cohorts (E-MTAB-13530, GSE267960): Spots were filtered (genes with < 10 counts removed, spots with < 200 total counts removed). E-MTAB-13530 yielded 57,238 spots; GSE267960 yielded 9,403 spots. Cell type annotation for the spatial validation cohorts was performed using RCTD (“full” mode), with the annotated GSE131907 scRNA-seq data serving as the reference. MISTY analysis was performed using the same parameters as for the primary spatial cohorts.

During the development of the prognostic model, differential gene expression analysis was first performed. Differentially expressed genes (DEGs) between tumor and normal tissues in the TCGA-LUAD dataset were identified using the criteria of LogFC > 1.5 and Padj < 0.05. Concurrently, marker genes for CYP27A1⁺TAMs were selected based on thresholds of LogFC > 0.25 and Min.pct = 0.25. The intersection of these two gene sets yielded 108 candidate genes.

Univariate Cox regression analysis was then applied to these 108 genes within the TCGA-LUAD cohort to identify those with a significant impact on prognosis. We selected the 24 genes demonstrating the strongest prognostic associations. However, due to the absence of CD302 and C4orf48 expression data in several validation cohorts, a final set of 22 genes was utilized for model construction.

For the model building phase, we systematically evaluated 101 distinct combinations of machine learning algorithms, including Random Forest (25), CoxBoost (26), Elastic Net (27), Gradient Boosting Machine (GBM) (28), Lasso (29), plsRcox (30), Ridge (31), StepCox, SuperPC (32), and Survival-SVM (33). The optimal algorithm combination was selected based on the highest average C-index across cohorts, forming the final CMRS (CYP27A1⁺Macrophage Risk Score). To assess the robustness of the CMRS model and address potential overfitting concerns, we performed comprehensive stability and calibration analyses. Bootstrap resampling with 1000 iterations was conducted on the TCGA-LUAD training cohort to evaluate the consistency of C-index values. Additionally, we examined C-index distributions across 500 different random seeds (80:20 training-validation splits) to assess model performance variability. Calibration curves were generated for 1-, 3-, and 5-year survival predictions using the bootstrap method across all validation cohorts. The Hosmer-Lemeshow test was applied to statistically evaluate the agreement between predicted and observed survival probabilities at each time point. The robustness and superior performance of CMRS were subsequently validated through multiple approaches: comparative analysis of C-indices against previously published models, assessment of AUC values across different validation cohorts, and demonstration of significant prognostic differences between high-risk and low-risk patient groups stratified by the median CMRS. The coefficients (COEF) of the genes in the CMRS have been provided in Supplementary Material 1.

We achieved overexpression of CYP27A1 in RAW by lentiviral transduction. The treatment group was infected with CYP27A1-overexpressing lentivirus at an MOI of 10 (RAW-OE), the negative control group was treated with NC lentivirus at an MOI of 10(RAW-NC), and the blank control group received no additional treatment (RAW). All cells were maintained in a humidified incubator at 37 °C with 5% CO₂.

All animal procedures in this study were approved by the Animal Ethics Committee of Tianjin Chest Hospital. The C57BL/6 female mice were obtained from Si Bei Fu Biotechnology (Beijing, China) at 6–8 weeks of age were subcutaneously injected in the right inguinal region with Lewis lung carcinoma (LLC) cells resuspended in ice-cold PBS, at a density of 1×10⁶ cells per mouse. Tumor growth was monitored weekly by measuring the longest diameter with a vernier caliper. Starting from the second week after LLC injection, 1×10⁶ RAW were separately administered to mice in the experimental group, vector control group, and blank control group. On day 28 after LLC cell injection, mice were euthanized by inhalation of excess CO2, with the CO2 flow rate maintained at 30–50% of the chamber volume per minute. Death was confirmed by cessation of respiration, dilated pupils, and no response to toe pinch. Subsequently, tumors were immediately dissected and removed in full.

Total RNA was extracted from cells using TRIzol reagent (Thermo, #15596018) following the manufacturer’s instructions. cDNA was synthesized using the PrimeScript RT Master Mix kit (Vazyme, #R232-01). Quantitative real-time PCR was performed with SYBR Green Master Mix (Vazyme, #Q111-02). GAPDH mRNA was used as an internal control for normalization. Relative expression levels were calculated using the 2−ΔΔCt method. The primer used in the paper have been provided in Supplementary Material 2.

Tumor tissues were fixed, paraffin-embedded, and sectioned. Deparaffinization was performed using xylene, followed by rehydration through a graded ethanol series. Antigen retrieval was carried out by incubating the sections in citrate buffer at 65 °C for 2 hours. Endogenous peroxidase activity was quenched with 3% hydrogen peroxide. The sections were then blocked with 5% bovine serum albumin (BSA) and incubated with specific primary antibodies (1:500). The remaining procedures were performed using an immunohistochemistry kit (ZSBIO, PV-9000).

Tumor tissues were fixed, paraffin-embedded, and sectioned. Deparaffinization was performed using xylene, followed by rehydration through a graded ethanol series. The sections were then blocked with 5% BSA. Immunofluorescence staining was carried out using an Immunofluorescence System Kit (Yeasen,60410), with primary antibodies diluted at 1:200. Three antibodies were incubated using the TSA method. Cell nuclei were stained with DAPI. After staining, the sections were imaged using a 3DHISTECH fluorescence slide scanner and analyzed with SlideViewer software.

Western blotting was performed according to standard protocols. Proteins were extracted using RIPA lysis buffer, and concentrations were determined via the BCA assay. Primary antibodies used included those against CYP27A1 (1:2000) and GAPDH (1:5000). A goat anti-rabbit IgG HRP-conjugated antibody (1:5000) was used as the secondary antibody.

Bioinformatics analyses were performed using R (v4.4.0). Basic experimental data were processed with GraphPad Prism and ImageJ. For comparisons between groups, Student’s t-test or one-way ANOVA was applied for normally distributed data, while the Wilcoxon rank-sum test or Kruskal–Wallis test was used for non-normally distributed data. Survival analysis was conducted using the Kaplan–Meier method, and differences were assessed with the log-rank test. Correlations were evaluated using Spearman’s rank correlation coefficient. A p-value < 0.05 was considered statistically significant, with levels of significance denoted as *P < 0.05, **P < 0.01, and ***P < 0.001.

Single-cell and spatial transcriptomics data were obtained from Yan et al. and downloaded from Zenodo (9). The single-cell dataset comprised 19 samples, all from tumor tissues of NSCLC patients, including 7 from pre-operative biopsies and 12 from surgical resections. All patients had received chemotherapy combined with ICB therapy. Based on pathological response, they were categorized into 6 pathological responders (5 with pathological complete response, pCR, and 1 with major pathological response, MPR, collectively referred to as the pathological responders, PR group) and 13 non-major pathological responders (NMPR group). We performed re-quality control on the data, resulting in 148,812 high-quality cells. Batch effects were removed, followed by re-dimensionality reduction, clustering, and cell population identification (Supplementary Figure 1A). Based on classical markers, cells from the 19 samples were categorized into 10 clusters (Figures 1A, B). Epithelial cells were further classified into normal epithelial and malignant epithelial cells using inferCNV based on copy number variations (Supplementary Figures 1B, C). T cells and myeloid cells were subjected to separate re-dimensionality reduction, clustering, and subpopulation analysis (Supplementary Figures 1D, E). Macrophages were specifically named according to their corresponding markers as: CCL4L2⁺TAMs, CYP27A1⁺TAMs, FGL2⁺TAMs, LPL⁺TAMs, SPP1⁺TAMs (Figure 1C).

Fifteen corresponding VISIUM spatial transcriptomics samples were all obtained post-operatively, including 6 tumor samples and 9 stromal samples. Among the 6 tumor samples, 2 were from the tumor core (T_C), 3 were isolated tumor masses, and 1 was invasive tumor (all collectively referred to as T_I). The 9 stromal samples comprised 5 from MPR and 4 from NMPR patients. After filtering out low-quality spots and genes, 24,810 genes and 59,387 spots remained. To characterize the TIME by integrating the gene expression matrices and spatial distribution information from both datasets, we applied RCTD. Using the single-cell data corresponding to the spatial samples, we selected 16 representative clusters to annotate the spatial transcriptomics samples. The annotation results were satisfactory: tumor samples were predominantly populated by tumor cells, with clear distinction between normal and malignant epithelium, while stromal samples were highly enriched with fibroblasts (Figures 1D–G). Notably, the deconvolution results were equally robust in identifying tertiary lymphoid structures (TLS), accurately pinpointing their immune cell composition, which was primarily comprised of B cells and T cells, consistent with literature reports (34–36).

By integrating and deconvoluting matched single-cell and spatial transcriptomics data, we reconstructed a detailed spatial cellular atlas. This comprehensive approach enables us to explore the tumor immune microenvironment holistically, from the complex genetic sequences to the intricate spatial architecture.

To better understand the immune microenvironment of spatial samples at VISIUM resolution, we employed ISCHIA niche analysis to investigate the TIME from the perspective of cell communities. Using a co-occurrence analysis method to identify spatial neighborhoods, this clustering approach can recapitulate spatial aggregation and interaction tendencies, clustering a total of 59,387 spots into 9 distinct niches. Niches are local microenvironments within the tumor formed by specific cell populations (e.g., cancer cells, immune cells, stromal cells) and their interactions, termed Composition Cluster 1-9 (CC1-9) (Figures 2A–C; Supplementary Figures 2A, B). Based on the predominant cell types and spatial distribution tendencies, we categorized the 9 CCs into 4 classes. PROGENy pathway enrichment analysis and Uniform Manifold Approximation and Projection (UMAP) plots helped us quickly and intuitively understand the different niches (Figures 2D, E).

CC2, CC4, and CC6 belong to the first class, designated as tumor-enriched niches. Among these, CC2 has the highest proportion of malignant cells and is located at the furthest edge in the UMAP distribution. Tumor-enriched niches highly express pathways such as MAPK, EGFR, PI3K, and estrogen signaling, reflecting the active growth, proliferation, and invasive state of tumor cells (37–40). Furthermore, CC2 exhibits the lowest expression of tumor suppressor pathways like P53 and Trail, and the highest expression of VEGF (41–43), providing further evidence that it possesses the highest degree of tissue instability and metastatic potential.

CC7 and CC8 constitute the second class, termed tumor-stroma niches. They contain abundant iCAFs, plasma cells, and B cells (Figure 2C; Supplementary Figures 2C, D) and are primarily located in the stroma of T_I samples. Compared to tumor niches, they upregulate relatively more pro-apoptotic pathways, as well as inflammatory pathways like NFκB and TNFα (44–46). They foster complex inflammatory components within the stroma that regulate the tumor microenvironment.

CC5 is the third class, identified as the normal epithelium-enriched niche. It contains relatively few immune cells, resembling a normal lung epithelial environment. Its functional state also exhibits complexity: on one hand, it highly expresses the P53 pathway to maintain genomic stability; on the other hand, it highly expresses JAK-STAT signaling to regulate the immune system and cell proliferation. Additionally, it shows the highest expression of androgen signaling, which is opposite to the tumor-enriched niches. This reflects the ability of the CC algorithm to distinguish the tendencies of different cell types at the functional level.

CC1, CC3, and CC9 form the fourth class, termed immune-stromal niches. They are heavily clustered in stromal samples and enriched with various immune cells, particularly macrophages. It can be said that the functional differences of macrophages primarily constitute their main distinctions (Figure 2D). Notably, these niches showed minimal expression of invasiveness markers. Conversely, they also highly express P53. Their points of difference lie in the fact that CC9 has the highest expression of the Trail pathway, suggesting that immune cells within CC9 may tend to exert anti-tumor effects by promoting apoptosis. In contrast, CC1 and CC3 show highly activated inflammatory pathways, including NFκB and TNFα, indicating they might perform anti-tumor functions by activating more inflammatory responses. It is important to note that while CCs of the same type share functional similarities, different CCs possess their own unique heterogeneity, particularly in terms of cell functional states, which warrants further exploration.

Although both CC1 and CC9 are distributed within the tumor stroma, their distribution tendencies differ across samples with different pathological responses (Figure 2B). This quantifiable, clinically correlated distributional difference demonstrates that spatial niches can not only help us quickly understand the composition of the tumor microenvironment from a holistic perspective but can also serve as a comparable basic unit for deconstructing the heterogeneity of the tumor microenvironment. In summary, using 15 spatially transcriptomic samples annotated via deconvolution, we applied the ISCHIA algorithm to construct 9 spatial niches. By integrating transcriptomic functional changes, cellular composition, and spatial tendencies, we reconstructed the spatial architecture of the TIME from a new perspective and preliminarily compared their differences.

To decipher the tumor ecosystem, we first focused on analyzing the tumor-enriched spatial niches. In the tumor ecosystem, CC2 is the most accumulated with tumor cells, predominantly composed of tumor cells and myCAFs, with very few immune cells, exhibiting a typical immune desert landscape. CC2 strictly exists in the tumor core and a small portion of infiltrating tumors. It is not present in stromal samples, consistent with pathological annotations, reflecting the reliability of annotations and ecological niche analysis. Similarly, CC6 is only enriched in the tumor core and infiltrating tumors, but compared to CC2, it is more distributed within infiltrating tumors and less in the tumor core. CC4 is more complex; its distribution in the tumor is similar to CC6, mainly present in infiltrating tumors. Interestingly, CC4 is the only one among the tumor-enriched ecosystems that is relatively distributed in the stroma, especially in MPR samples (Figure 2B). Since CC4 and CC6 are located at the tumor margins, often at the interface between tumor and normal tissue, understanding their roles is crucial for deepening our exploration of the tumor microenvironment.

In terms of cellular composition, both are lacking T/B cells and are enriched with immunosuppressive cells such as SPP1⁺TAMs, Tregs, and myCAFs (Figures 3A, Supplementary Figure 3A). Notably, myCAFs tends to associate with tumor infiltration, whether in the tumor core or infiltrating tumor, consistent with existing studies (47). Comparing the spatial distributions, CC6 and CC2 are more directly related, often encapsulating the tumor core, whereas CC4 is relatively more marginal and often located at the periphery, surrounding CC4 and CC2. Analyzing ecosystem co-occurrence internally, we found that Tregs and other cell types, especially macrophages, are closely associated. Additionally, in CC4, co-occurrence of SPP1⁺TAMs and tumor cells was observed, These spatial patterns align with the known pro-tumor functions of these immunosuppressive cells. Specifically, SPP1⁺TAMs are reported to collaborate with myCAFs in remodeling the extracellular matrix and establishing an immunosuppressive microenvironment (48). Tregs, by suppressing effector T cell activity, create an immune-permissive environment for tumor growth and infiltration (49). Thus, the spatial co-occurrence of SPP1⁺TAMs, Tregs, and myCAFs likely creates a synergistic immunosuppressive circuit that enhances tumor cell survival and invasive potential within the TME. (Supplementary Figure 3B). To further confirm cell interactions within ecosystems, we used CellChat for analysis of nine different CCs. To avoid errors and overinterpretation, we set interaction.range = 250 and shielded signals at longer distances, focusing only on intra-ecosystem interactions (Supplementary Figure 3C). Results for CC4, which has fewer tumor cells, indicated that cell interactions mainly promote tumor migration and Extracellular matrix (ECM) remodeling (Figure 3B). Tumor invasion is mainly mediated by ligands such as VEGF, COMP, ANGPTL2, and CEACAM6, which promote angiogenesis, inhibit immune responses, and reduce tumor cell adhesion, facilitating infiltration (50, 51). Specifically, VEGF activating FLT1 (VEGFR1) reduces T cell infiltration and cytotoxic immune responses (52); COMP expression aids the infiltration of M2 macrophages and fibroblasts; CEACAM6 and ANGPTL2 promote tumor infiltration by decreasing cell adhesion and downregulating MHC-I expression. Some studies suggest that stromal fibrosis promotes tumor infiltration (53), and ECM remodeling in CC4 mainly operates through collagen cross-linking and integrin pathways, which facilitate fibrosis (54, 55). Notably, COL4A1 promotes tumor migration via FAK-Src signaling, and its expression is associated with poor prognosis (56, 57).

Another infiltrating ecosystem, CC6, also primarily involves immunosuppression and ECM remodeling, but with a different regulatory network (Figure 3C). Besides collagen-integrin signaling, CD44 plays an important role in matrix reprogramming; activation of CD44 via Nanog signaling maintains stem cell properties (58). SDC1 serves as a central node in CC6 metabolic network; as a heparan sulfate proteoglycan, its high expression is associated with poor tumor prognosis (59). Its activation increases tumor migration, directly aiding infiltration (60, 61); it can also be activated by COL1A1, exacerbating fibrosis (62). Many previous studies have linked increased extratumoral matrix rigidity with worse prognosis (63). From the perspective of ecological niche integration with cellular interactions, we re-examine this process, especially emphasizing the crucial functions of SDC1 at the edges of tumor infiltration. It may act as a bridge in tumor invasion and ECM formation. Overall, the internal interactions of CC4 and CC6 provide new insights into mechanisms of immune evasion and metastasis in infiltrating tumors, as well as the unique pro-fibrotic role of myCAFs in ECM.

To investigate the potential progression from the tumor core to the invasive margin, we constructed a spatial trajectory from CC2 through CC6 to CC4 based on SPATA2, originating from CC2, passing through CC6, and ending at CC4, to compare spatial changes in tumor ecosystem development (Figure 3D). Genes along this trajectory were subjected to HALLMARK pathway enrichment; results showed that pathway enrichment perfectly followed the trajectory changes, with CC6 exhibiting overexpression of gene sets characteristic of CC2 and CC4 (Figure 3E). Compared to CC6, CC2 upregulated pathways mainly involved MYC_TARGETS and oxidative phosphorylation, indicating vigorous proliferative capacity. Conversely, anti-tumor pathways such as inflammation, interferon, antigen presentation, P53, complement, and apoptosis were almost uniformly under-expressed, consistent with the immune desert landscape characterized by sparse immune cells. CC4 specifically upregulated pathways related to angiogenesis, epithelial-mesenchymal transition (EMT), and myofiber formation, aligning with typical infiltrating tumor features (64), consistent with cell interaction results. We further examined the spatial expression patterns of ligands and receptors critical in cell communication within CellChat (Figures 3F, G). The results were intriguing: SDC1 and ITGA2 are highly expressed specifically in the tumor microenvironment; COMP is specifically distributed at the tumor margins. COMP is considered related to immune suppression, and some studies suggest that activation of COMP promotes EMT (65). Interestingly, in CC4, COMP acts as a ligand interacting with SDC1, indicating it may play a crucial role in tumor infiltration and progression (66). Given SDC1’s central role in communication within CC4, inhibiting their interaction might hinder tumor infiltration.

In summary, our analysis delineates a spatial continuum from the proliferative, immune-desert tumor core (CC2) to the invasive margin (CC4), which is characterized by EMT, angiogenesis, and specific ligand-receptor interactions such as COMP-SDC1. However, immune cells are impeded by fibroblast-formed fibrous barriers and a series of immunosuppressive mechanisms at the tumor margin. These immune cells cannot enter the tumor center to exert anti-tumor effects, while infiltrating tumors can generate new blood vessels and continue expanding via mechanisms like EMT. Therefore, finding therapeutic strategies targeting the tumor margin microenvironment might greatly enhance anti-tumor efficacy. We propose that inhibiting SDC1 and COMP could be potential therapeutic approaches. In conclusion, by dissecting different tumor components and integrating spatial interactions and trajectories, we offer a new perspective on the infiltration characteristics of the tumor microenvironment.

The stroma in the TIME also plays a critical role in determining therapeutic efficacy (67). To investigate the key cellular communities (CCs) influencing immunotherapy outcomes, we compared the composition of CCs between different treatment response groups (Figure 2B). To evaluate the independent predictive value of CC9, we performed multivariate logistic regression analysis in the IMvigor210 cohort. The CC9 signature score for each sample, quantified using ssGSEA, was incorporated into the model along with tumor stage and Immune Phenotype. The analysis confirmed that a higher CC9 signature score was significantly associated with increased likelihood of treatment response, identifying it as an independent favorable factor(Supplementary Figure 3D). These results suggest that deconstructing CC9 may help uncover the mechanisms underlying differential immunotherapy responses.

In the UMAP visualization, CC9 and CC1 were closely distributed, and both were predominantly located in the stromal region. Additionally, the proportions of cell types within these two communities were generally similar. However, CC9 exhibited a significantly higher proportion of CD8⁺T cells than CC1 (Figure 3A). To dissect the cellular basis for this differential T cell abundance and explore the potential multi-cellular cooperation within CC9, we performed GSVA to score the activity of cytotoxic T lymphocytes (CTLs), as well as the presence of TLS and other key cell types within this community (Supplementary Figure 4). Subsequent correlation analysis revealed a strong positive association between CYP27A1⁺TAMs and CTL activity. In contrast, SPP1⁺TAMs and myCAFs showed the expected negative correlations. Notably, the TLS signature displayed the strongest positive correlation with CTL activity. Significant positive correlations were also observed between CTL activity and signatures for B cells and iCAFs. These findings collectively suggest that within CC9, multiple stromal and immune components—including TLS, B cells, iCAFs, and CYP27A1⁺TAMs—likely act in concert to foster an anti-tumor microenvironment, constituting a multi-cellular cooperative immune niche. Given its strongest correlation with CTLs among individual cells (Supplementary Figure 4) and its marked enrichment in CC9 versus CC1 (Figure 2C), we focused onCYP27A1⁺TAMs.

Strikingly, the Post-MPR group showed a significantly higher number of CYP27A1⁺TAMs than the Post-NMPR group (Figures 4A), implying that CYP27A1⁺TAMs may play a crucial role in immunotherapy. To further investigate their functional significance, we performed GSVA on macrophages using key KEGG gene sets (Figure 4B). As expected, FGL2⁺TAMs exhibited M1-like properties, demonstrating strong pro-inflammatory and cytotoxic capabilities, along with antigen presentation. In contrast, SPP1⁺TAMs showed minimal pro-inflammatory and antigen-presenting functions, consistent with previous findings (68). However, the tumor microenvironment is more complex than previously thought, and understanding TAMs solely from a pro-/anti-inflammatory perspective is insufficient (7, 69). Interestingly, although CYP27A1⁺TAMs did not express traditional M1 macrophages markers, they exhibited the highest expression of MHC pathway genes. Moreover, in the TCGA-LUAD cohort, patients with high expression of CYP27A1 and CYP27A1⁺TAMs had significantly better prognosis than those with low expression (Figures 4C, D).

We then utilized Monocle2 to analyze macrophages functional dynamics from a developmental perspective, focusing on patients with pCR (Figure 4E). The results showed that CYP27A1⁺TAMs and LPL⁺TAMs had similar state distributions and were almost mutually exclusive with SPP1⁺TAMs, suggesting that CYP27A1⁺TAMs may represent a macrophage subset with dual roles in lipid metabolism and tumor suppression within the TIME. Next, we examined macrophages developmental trajectories by analyzing gene expression differences at branching point 1, which revealed two distinct clusters (Figures 4F, G). Functional enrichment analysis of these clusters showed that BP terms were primarily associated with two major categories: cluster 1 was enriched in lipid metabolism-related processes, while cluster 2 was associated with adaptive immune responses (Figure 4H). Based on the state distributions of CYP27A1⁺TAMs and FGL2⁺TAMs, we propose that CYP27A1⁺TAMs are a macrophage subset with dual functions in immune response and lipid metabolism regulation. Here, we first rigorously assessed the specificity of CYP27A1 to ensure it serves as a unique marker for CYP27A1⁺TAMs (Supplementary Figure 5A). Since traditional alveolar macrophages are known to possess both immune and lipid metabolic functions, we further compared the expression of established alveolar macrophage markers across several macrophage subtypes. This analysis confirmed that CYP27A1⁺TAMs are distinct from alveolar macrophages (Supplementary Figure 5B). Having established this distinction, we proceeded to analyze the characteristics of CYP27A1⁺TAMs. Given that CYP27A1⁺TAMs exhibited the strongest antigen-presenting capacity, we hypothesized that they may primarily exert their immune functions through antigen presentation and T cell activation. To validate whether CYP27A1⁺TAMs interact with T cells, we performed cell communication analysis. Since spatial-level CellChat accounts for actual spatial distances and provides more realistic interaction insights, we applied it to spatial transcriptomics samples. The results were consistent with single-cell data: CYP27A1⁺TAMs played a critical role in antigen presentation to T cells, particularly CD8⁺T cells, via both MHC-I and MHC-II pathways (Figures 5A, C).Previous studies have suggested that after PD-1 immunotherapy, macrophages interact with CD8⁺T cells via CXCL9/CXCR3 to enhance T cell recruitment and improve immunotherapy efficacy (70). We therefore examined whether CXCL9/CXCR3 were highly expressed in CYP27A1⁺TAMs and validated their co-localization in spatial samples. The results confirmed that CYP27A1⁺TAMs specifically highly expressed CXCL9/CXCR3 (Figure 5B). More importantly, CYP27A1⁺TAMs were spatially co-localized with CD8⁺T cells (Supplementary Figures 5C, D), and the spatial patterns of CXCL9/CXCR3 significantly overlapped with this co-localization (Figure 5E). Given the important role of CD86 in T cell activation and functional differentiation, we analyzed CD86 interaction patterns. At a significance threshold of p = 0.01, CYP27A1⁺TAMs dominated CD86 interactions, extensively activating CD8⁺T and CD4⁺T cells but not Tregs (Figure 5D; Supplementary Figures 5E, F). These CellChat results indicate that CYP27A1⁺TAMs recruit and activate T cells, especially CD8⁺T cells, through multiple pathways, which may be the primary mechanism underlying their anti-tumor effects.

To investigate whether the interaction between CYP27A1⁺TAMs and CD8⁺T cells is widespread, we used the spatstat package and calculated their co-localization in each sample based on Ripley’s K function method. The results indicated that in every sample, the tendency for co-localization between the two far exceeded random distribution, and both the Aggregation Index and Cohen’s d were highly significant(Supplementary Figure 6). Then we performed multi-view spatial co-localization analysis using MISTY, a machine learning-based framework (24). MISTY employs interpretable machine learning algorithms to quantify the co-localization patterns of cells or genes across multiple spatial transcriptomics samples within defined spatial contexts. The definitions of the three distance views were adapted from the method described by Ateeq M. Khaliq et al. (47): the intra-view (55 µm, within a spot), the juxta-view (200 µm, capturing direct and short−range interactions), and the para-view (3000 µm, assessing long−range or tissue−level organization). The para−view excluded co−localization within 200 µm to avoid overlap with the juxta−view (Figure 5F). Using a cutoff value of 1.5, tumor cells, iCAFs, and CYP27A1⁺TAMs emerged as core players in cellular co-localization, playing crucial roles (Figures 5G, H). These results further support the co-localization between CYP27A1⁺TAMs and CD8⁺T cells, and their joint co-localization with tumor cells highlights the indispensable role of CYP27A1⁺TAMs in the anti-tumor microenvironment. Notably, although the number and pattern of spatial co-localizations changed under juxta-view and para-view, CYP27A1⁺TAMs consistently maintained a relationship with CD8⁺T cells, demonstrating robust spatial dependency (Supplementary Figures 5G–J). Additionally, cells of the same type tended to co-localize, reflecting harmonious spatial organization.

Integrating these findings, we compared the reasons for differences in immunotherapy efficacy and identified CYP27A1⁺TAMs as a key factor. These TAMs exhibit characteristics distinct from traditional M1 macrophages and play a critical anti-tumor role in the TIME by recruiting CD8⁺T cells.

To further investigate the impact of immunotherapy on the regulation of CYP27A1 expression, we performed GSVA on CYP27A1⁺TAMs in different states. Surprisingly, the antigen presentation capacity did not increase after immunotherapy (Figure 6A). To explore this discrepancy, we compared the differentially expressed genes in CYP27A1⁺TAMs before and after immunotherapy. The results showed significant upregulation of CYP27A1, LXR, and its target genes (71) post-treatment (Figure 6B). Additionally, the expression of SFTPA1 and SFTPA2 was markedly increased. Previous studies suggest that high expression of SFTPA1 and SFTPA2 in macrophages enhances their survival and resistance to apoptosis (72). In line with these findings, following combined chemoimmunotherapy, various apoptosis pathways were downregulated in CYP27A1⁺TAMs, while pathways such as CXCR GNB G PI3K AKT SIGNALING PATHWAY and KEAP1 NRF2 SIGNALING PATHWAY (73, 74) were significantly upregulated. These pathways promote cell survival and enhance resistance to oxidative stress, which may be key to the sustained anti-tumor role of CYP27A1⁺TAMs in patients receiving combination therapy. LXR (Liver X Receptor α and β), nuclear hormone receptors encoded by NR1H2 and NR1H3, are ligand-activated transcription factors specifically expressed in the liver and macrophages (75). Recently, the team of Srustidhar Das demonstrated that LXR has anti-tumor effects mediated through the recruitment of B cells, T cells, and the formation of tertiary lymphoid structures. They also identified CYP27A1 as the most important agonist of LXR (76). Interestingly, in our study, both ligands and receptors of LXR were significantly upregulated in CYP27A1⁺TAMs after immunotherapy. This led us to investigate whether high LXR expression enhances the ability of CYP27A1⁺TAMs to recruit T cells.

We divided CYP27A1⁺TAMs into two groups based on LXR expression levels: TAMs-LXRhi and TAMs-LXRlow. Using RCTD for re-deconvolution and annotation of spatial transcriptomics samples, we found that the distribution of TAMs-LXRhi and TAMs-LXRlow in the stroma showed no significant difference (Figures 6C, D). However, within tumor samples, TAMs-LXRhi was significantly more abundant than TAMs-LXRlow (Figures 6E, F). Moreover, we observed that TAMs-LXRhi exhibited stronger co-localization with CD8⁺T cells (Figure 6G). Compared to TAMs-LXRlow, TAMs-LXRhi displayed a more activated phenotype, suggesting that LXR activation may enhance the function of CYP27A1⁺TAMs. To broadly validate this observation across all samples, we used MISTY to compare the strength of co-localization between CYP27A1⁺TAMs and T cells under different LXR expression levels (Figure 6J). Furthermore, we validated the correlation between LXR expression in CYP27A1⁺TAMs and CD8⁺T cells in both single-cell and TCGA cohorts (Figures 6H, I). The results confirmed that TAMs-LXRhi had significantly stronger co-localization with T cells than TAMs-LXRlow, indicating that high LXR expression is associated with enhanced T cell recruitment ability. To identify the downstream mediators responsible for LXR-driven T cell recruitment, we analyzed the correlation between canonical LXR target genes (including APOE, ABCG1, and ABCA1) (77) and CD8⁺T cell markers (CD8A and CD8B) in the TCGA-LUAD dataset. This analysis revealed significant positive correlations between these LXR targets and CD8⁺T cell abundance (Supplementary Figure 7), with APOE showing the most pronounced association. This suggests APOE may be a pivotal LXR-induced factor enhancing the T cell recruitment capacity of these macrophages. Therefore, our data indicate that APOE, among other LXR-induced factors, may play a pivotal role in enhancing the CD8⁺T cell recruitment capacity of CYP27A1⁺TAMs following LXR activation.

In summary, by comparing the genetic and pathway changes in CYP27A1⁺TAMs before and after treatment, we found that anti-apoptotic capacity and LXR expression were significantly increased post-treatment. TAMs-LXRhi represents an activated state of CYP27A1⁺TAMs with stronger survival ability and enhanced CD8⁺T cell recruitment capacity.

To assess the generalizability of our findings, we performed validation across multiple independent external datasets. The first validation dataset was GSE131907. After quality control (Supplementary Figure 5A), dimensionality reduction, clustering, and batch effect removal (Supplementary Figure 5B), we identified 10 major clusters (Figure 7A; Supplementary Figure 5C). Consistent with our own data, CYP27A1 was specifically highly expressed in macrophages, and NR1H2 and NR1H3 showed strong correlation with CYP27A1 expression levels (Supplementary Figure 5D). To assess the impact of CYP27A1 expression levels on macrophage functional states, we divided mononuclear phagocytes into quartiles based on CYP27A1 expression (Supplementary Figures 5E, F). Notably, macrophages in lung tissue exhibited higher overall CYP27A1 expression compared to those in pleural effusion, brain tissue, and lymph nodes. The groups were named, from highest to lowest expression: CYP27A1-high, CYP27A1-intermediate, CYP27A1-low, and CYP27A1-negative. We then scored their functional states using KEGG pathways (Figure 7B). Mirroring our internal results, the CYP27A1-high group demonstrated the strongest antigen presentation capacity for both MHC-I and MHC-II pathways and a greater reliance on β-oxidation rather than glycolysis at the metabolic level. Subsequently, using the annotations from GSE131907, we deconvoluted and annotated 57,238 spots from a spatial transcriptomics dataset (E-MTAB-13530), comprising lung adenocarcinoma and normal tissue. The results revealed a high spatial correlation between CYP27A1 and CD8A/CD8B (Figure 7C; Supplementary Figures 5G, H). We further employed MISTY to quantify co-localization differences between cell types across all samples, which again confirmed that macrophages with high CYP27A1 expression are more likely to co-localize with CD8⁺T cells in space (Figure 7D).

The second validation cohort’s single-cell data came from GSE223203. After quality control and batch effect removal (Supplementary Figures 6A, B), we applied the same clustering method as the original publication (Figure 7E; Supplementary Figure 6C). CYP27A1 expression was remarkably consistent across datasets, specifically appearing in myeloid cells and even serving as a potential marker for this lineage. Cluster 8, lacking specific markers, was ultimately defined as ‘Unknown’. For macrophages within the myeloid cell population, we applied the same quartile-based grouping (Supplementary Figure 6D). Notably, non-responder (NR) patients had a relatively higher proportion of macrophages in the CYP27A1-negative group, suggesting that low CYP27A1 expression in macrophages may be associated with resistance to immunotherapy. Subsequent KEGG pathway analysis of the quartile-grouped macrophages (Figure 7F) showed that the CYP27A1-high group remained the primary contributor to antigen presentation and preferred β-oxidation for energy production. Considering the presence of ‘Unknown’ cells in GSE223203, we used GSE131907 to deconvolute and annotate another spatial dataset, GSE267960 (two samples, 9401 spots). Both at the gene and cell type level, CYP27A1 and CD8⁺T cells were in close spatial proximity (Figure 7G; Supplementary Figure 6E).

Most importantly, our findings were strongly corroborated by the Human Protein Atlas. First, we examined the organ and cellular distribution of CYP27A1 (Figure 7H), which confirmed prior reports of its high expression in hepatocytes (78). In the lung, CYP27A1 was specifically highly expressed in macrophages. Furthermore, only macrophages in the lung and pancreas highly expressed CYP27A1; this dual specificity makes CYP27A1⁺TAMs easier to identify and the results more reliable. In the Human Protein Atlas single-cell clustering, CYP27A1 is part of cluster 58 (Macrophages - Innate immune response) with a confidence score of 1 (Supplementary Figure 6F). (Confidence is the fraction of times a gene was assigned to the cluster in repeated clustering, and therefore reflects how strongly associated it is to the cluster. A confidence of 1 indicates that the gene was assigned to this cluster in all repeated clusters.) Cluster 58, Macrophages - Innate immune response, specifically highly expresses 420 genes (Supplementary Materials 3). Using these 420 genes, GSVA and Jaccard similarity scores were calculated for each cluster in the single-cell dataset. Both analyses showed that the highest scores were observed in the CYP27A1⁺macrophage cluster. (Supplementary Figures 6G, H), suggesting very similar functional states despite different microenvironments.

Examining the GO enrichment results for Cluster 58 (Figure 7I), its biological processes are primarily focused on innate immune response—namely, non-specific immune functions like phagocytosis and antigen presentation. Cellular component analysis reveals exceptionally well-developed lysosomes and endosomes, consistent with active phagocytic activity. Given that lysosomes are crucial for antigen processing and presentation in antigen-presenting cells, their proper function is essential for activating anti-tumor T cell immunity. Furthermore, its molecular functions are enriched in pattern recognition receptors and lipid metabolism, mirroring the characteristics of CYP27A1⁺TAMs in our study. This underscores the stability of CYP27A1 as a marker for a specific macrophage subset within the TIME that exerts anti-tumor effects by recruiting T cells.

In summary, through two external single-cell datasets, two spatial transcriptomics datasets, and integration with the extensive single-cell atlas from the Human Protein Atlas, we firmly establish that CYP27A1 is specifically highly expressed in lung macrophages and demonstrates a robust correlation with strong antigen presentation capacity.

To translate our findings into a clinically applicable tool, we developed a prognostic model based on the CYP27A1⁺TAMs gene signature. Given that high expression of genes associated with CYP27A1⁺TAMs correlates with improved patient survival (Figures 4C, D), we sought to develop a gene expression-weighted prognostic model. This model, designated the CYP27A1⁺Macrophage Risk Score (CMRS), stratifies patients into distinct risk groups.

Our analysis identified 2,285 and 1,289 differentially expressed genes from the TCGA-LUAD dataset and CYP27A1⁺TAMs, respectively. The intersection of these gene sets yielded 108 candidate genes (Figure 8A). Subsequent univariate Cox regression analysis of these genes within the TCGA-LUAD cohort identified the top 24 genes most significantly associated with prognosis (Figure 8B), which were utilized for model construction. Notably, among these 24 prognostically significant genes, only one was identified as a hazard factor, underscoring the predominant protective role and robust prognostic predictive value of CYP27A1⁺TAMs.

For model development, the TCGA-LUAD cohort served as the training set, while the GSE13213, GSE31210, GSE42127, and GSE50081 datasets were employed as independent validation cohorts. By systematically evaluating 101 distinct machine learning algorithm combinations, we developed a robust predictive model. The predictive performance of each model was assessed by calculating the average C-index across all cohorts (Figure 8C). Our analysis identified the Elastic Net algorithm with an alpha of 0.9 as the optimal model, demonstrating superior integrated performance across the TCGA-LUAD, GSE13213, GSE31210, GSE42127, and GSE50081 cohorts, with a mean C-index of 0.648. Consequently, this algorithm was selected to construct the final CMRS scoring system, which incorporates 7 protective genes. The relative feature importance within the model is depicted in Figure 8D.

We first compared the performance of CMRS against previously published models by evaluating their C-indices across different cohorts. The results demonstrated that CMRS consistently achieved superior C-index values in both the training and all validation sets (Figure 9E). To further evaluate the prognostic significance of the CMRS, Kaplan-Meier survival analysis was performed. Patients in each cohort (including a meta-cohort combining GSE13213, GSE31210, GSE42127, and GSE50081) were stratified into high-risk and low-risk groups based on the median CMRS score. The analysis revealed that the low-risk group exhibited significantly more favorable overall survival compared to the high-risk group in every cohort. Furthermore, time-dependent receiver operating characteristic (ROC) analysis for 1-, 3-, and 5-year survival confirmed the model’s predictive accuracy against actual observed outcomes (Figures 9A–D). These consistent findings across multiple independent datasets support the robustness, generalizability, and reproducibility of the CMRS across diverse clinical populations.

To further evaluate the robustness of the CMRS and address potential concerns regarding overfitting, we conducted comprehensive stability and calibration analyses. First, bootstrap resampling (n=1000) within the TCGA-LUAD training cohort was performed to assess the consistency of the C-index, which yielded a stable and tightly distributed C-index (Supplementary Figure 8A), indicating good internal stability of the model. Second, we examined the distribution of C-indices obtained from 500 repeated iterations with different random seeds for splitting the training and validation sets (Supplementary Figure 8B); the results showed minimal variance, confirming the robustness of our model selection against random data partitioning. Finally, calibration curves for 1-, 3-, and 5-year survival predictions were constructed. The Hosmer-Lemeshow test was applied to statistically evaluate the agreement between predicted and observed survival probabilities at each time point across cohorts. The curves demonstrated that the survival probabilities predicted by CMRS were consistently in close alignment with the actual observed outcomes (Supplementary Figures 8C, D), indicating excellent predictive calibration.

In summary, we constructed and validated the CMRS using the TCGA cohort and multiple independent validation sets. This model not only demonstrates performance superior to existing prognostic tools but also provides a robust instrument for prognostic stratification in LUAD, potentially offering a transformative approach for patient risk assessment and treatment planning.

To functionally validate the anti-tumor role of CYP27A1⁺TAMs, we performed in vivo experiments. To investigate the impact of macrophages with high CYP27A1 expression on tumor growth, we established a C57 mouse model. Briefly, We selected 6–8-week-old C57BL/6 mice and subcutaneously injected them with LLC cells to establish subcutaneous tumors. Two weeks later, RAW from different groups were injected.

First, we measured the transcriptional level of CYP27A1 via PCR. The lentiviral overexpression was successful, resulting in a significant increase in CYP27A1 expression in the RAW-OE group (Figure 10C). We then compared differences at the protein level using Western blot (WB), which confirmed that the CYP27A1 protein level in the RAW-OE group was also significantly higher than in the RAW-NC group (Figure 10A).

Subsequently, we intravenously administered these three types of macrophages via the tail vein into C57 mice that had borne subcutaneous tumors for two weeks. We then monitored tumor size over four weeks and measured tumor weight at the fourth week. The results demonstrated that RAW-OE significantly inhibited tumor cell growth, indicating that macrophages with high CYP27A1 expression within the TIME indeed exert an anti-tumor effect (Figures 10B, D, E).

To investigate whether this effect is mediated by CD8⁺T cells, we examined the differences from a morphological perspective. We performed serial sectioning of paraffin-embedded tumor tissues for subsequent Hematoxylin and Eosin (H&E) staining and immunohistochemistry (IHC) with an anti-CD8 antibody. In H&E-stained sections, potentially due to the nature of the subcutaneous tumor model, we did not observe significant differences in tertiary lymphoid structure formation. Furthermore, although the number of CD8⁺T cells in the tumor core of the RAW-OE group was somewhat higher than in the RAW-NC group, this difference was not statistically significant. The most notable distinction was observed in the inflammatory reaction zone at the tumor margin. This area was infiltrated by numerous neutrophils. In the RAW-OE group, a substantial number of CD8⁺T cells were clustered around the looser regions of this inflammatory zone. In contrast, very few T cells were visible inside or outside the inflammatory area in the RAW-NC group (Figures 11A–D).

To further explore whether CYP27A1⁺TAMs play a key role in recruiting CD8⁺T cells, we performed multiplex immunofluorescence (IF) on tumors from the RAW-OE and RAW-NC groups. The results revealed that the total number of CD68⁺macrophages was similar between the two groups. However, the RAW-OE group exhibited a significantly greater number of CYP27A1⁺, CD68⁺cells, alongside a markedly increased population of CD8⁺T cells. Spatial analysis confirmed the co-localization of these CYP27A1⁺TAMs and CD8⁺T cells (Figures 11E, F). This indicates that CYP27A1⁺TAMs play a crucial role in recruiting CD8⁺T cells within the TIME.

In conclusion, through our animal model, we have functionally and morphologically validated the discoveries from our single-cell and spatial transcriptomic analyses: macrophages with high CYP27A1 expression exert their anti-tumor function, at least in part, by recruiting more CD8⁺T cells to kill tumor cells.