Integrating single-cell and bulk transcriptomes identifies B cell features associated with neoadjuvant chemoradiotherapy sensitivity in rectal cancer.

Introduction
Rectal cancer ranks among the four most prevalent cancers with the highest mortality rates globally, accounting for 3.2% of all cancer-related deaths1.For patients with Locally Advanced Rectal Cancer (LARC), Neoadjuvant chemoradiotherapy (nCRT) followed by radical surgery—specifically total mesorectal excision —has been established as the standard treatment modality according to the National Comprehensive Cancer Network (NCCN) clinical practice guidelines2. This treatment strategy has demonstrated efficacy in reducing local recurrence, with 8%–20% of complete responders achieving organ preservation3–5. However, not all LARC patients benefit from nCRT, as only 15% to 27% of patients receiving nCRT achieve pathological complete response (pCR)3. approximately 20% of patients with LARC demonstrate a suboptimal response to chemoradiotherapy. Consequently, extensive efforts have been undertaken to identify sensitive, specific, and reliable biomarkers for predicting nCRT response, particularly pCR status6. to date, no precise and universally accepted biomarkers for response prediction have been identified. Consequently, there remains an urgent clinical need to discover novel biomarkers for patient stratification and disease subtyping.
In recent years, a growing body of research has indicated that tumor resistance is not solely attributed to the tumor itself; the tumor immune microenvironment also plays a crucial role in this process. nCRT not only directly kills tumor cells but also enhances antitumor efficacy by promoting the release of tumor-associated antigens, thereby boosting the immune response7,8.B cells, as professional antigen-presenting cells (APCs), can uptake and process tumor antigens while providing co-stimulatory signals to CD4⁺ T cells. This interaction enhances T cell-mediated antitumor immunity, thereby contributing to an effective immune response against tumors9. Transcriptomic sequencing of patients who are sensitive or refractory to neoadjuvant treatment for LARC has revealed that B cell- and interferon signaling-related features are highly expressed in sensitive patients. Immunohistochemical analysis has validated that B-cell infiltration levels are associated with the treatment efficacy of nCRT in patients with LARC10. Although B cells can promote antitumor immunity in many cases, studies have also identified specific B cell subsets—particularly regulatory B cells (Bregs)—that may contribute to tumor progression. Bregs can produce immunosuppressive factors such as IL-10 and TGF-β, which inhibit the cytotoxic function of CD8⁺ T cells and promote T cell exhaustion11. Therefore, the role of B cells in the nCRT for LARC requires further investigation.
The latest single-cell sequencing technologies can comprehensively characterize cells within complex tissues, unveiling critical subpopulations that impact tumor metastasis, treatment response, and survival outcomes. Nevertheless, the clinical application of single-cell sequencing is hindered by its high cost and labor-intensive nature, particularly in large cohorts. Leveraging phenotypic information from large-scale RNA-seq studies to guide the identification of subpopulations and biomarkers via single-cell sequencing has yielded remarkable insights into detecting highly disease-relevant cell subpopulations12. In this study, we employed the “Scissor” algorithm to integrate bulk data and single-cell transcriptome data, leading to the identification of B cells subpopulations associated with nCRT response. and further elucidated their interactions within the ME and their clinical significance.

Materials and methods

Data acquisition
The scRNA-seq dataset was downloaded from the NCBI Gene Expression Omnibus (GEO) database under the accession number GSE13246513.
bulk RNA expression dataset was downloaded from the GEO database under the accession number GSE40492, and GSE45404. In the GSE45404 dataset, patients were graded according to Tumor Regression Grade (TRG), with TRG 3, 4, and 5 considered resistant and TRG 1 and 2 considered sensitive. In the GSE40492 dataset, TRG information was provided by the original authors; patients with TRG 0 were classified as PCR, while all others were classified as non-PCR, TRG 0, 1, and 2 were considered the sensitive group, and TRG 3, 4, and 5 the resistant group.

Patients and samples
Patients with LARC who underwent nCRT followed by surgery in Fujian Medical University Union Hospital (Fujian, China) between 2015-01-01 and 2018-12-31 were consecutively enrolled in this study based on the following inclusion criteria: histologically proven rectal adenocarcinoma, no distant metastasis; Eastern Cooperative Oncology Group performance status of 0–1; fresh frozen biopsy tumor tissue samples available; neoadjuvant concurrent chemoradiotherapy with total mesorectal excision surgery. To further limit potential confounding factors along these lines, included patients were given a similar chemoradiotherapy strategy. A total of 64 patients with LARC met the criteria for inclusion. All the patients treated with preoperative radiotherapy plus Capecitabine standardized neoadjuvant chemoradiotherapy regime consisting of Capecitabine (Xeloda®, 825 mg/m2, Welwyn Garden City, United Kingdom), twice daily and one dose of radiotherapy (50.4 Gy in 28 fractions of 180 cGy/fr). The baseline characteristics of the patients are shown in Table 1. TRG was used to evaluate the pathological response to nCRT. Specifically, TRG 0 indicated no residual tumor cells; TRG 1, near-complete regression with tumor cells present individually or in small clusters; TRG 2, residual tumor cells accompanied by a desmoplastic response; and TRG 3, minimal or no regression. Patients with TRG 0 were classified into the pCR group, whereas all others (TRG 1–3) were categorized as non-pCR. The median follow-up was 27.4 months for PFS and 87.2 months for OS. Both the median PFS and median OS were not reached. The 3-year PFS rate was 68.9%, and the 5-year OS rate was 82.7%.

Single-cell RNA-seq data processing
The Seurat R package was used to perform unsupervised clustering of the single cells using the read count matrix as input. The batch effects among patients were eliminated with the harmony package14. Finally, 9146 B cells from the tumor and normal tissue were obtained. Using the “Scissor” algorithm, bulk sequencing from GSE45404 patients was correlated with single-cell data. The parameter alpha was set at 0.2 to filter out the most relevant nCRT Response. “Scissor” is specifically designed to integrate bulk clinical phenotypes with single-cell transcriptomic data, enabling the identification of cell subpopulations that are most strongly associated with treatment response. Traditional single-cell clustering methods cannot directly link cellular heterogeneity to patient-level outcomes, whereas Scissor provides a statistically robust framework to bridge this gap12. We used Monocle 2 to explore the pseudotime of each cell; the dynamic changes in gene expression were constructed using the plot pseudo time heatmap function15, and we used CellChat to investigate the ligand-receptor crosstalk between B cells and other cell clusters16.We used scMetabolism to analyze differences in metabolic pathway activity in single-cell datasets17.

Functional enrichment analysis
After findmarker identified subgroups of differentially expressed genes, ClusterProfiler was used to analyze between-group differences. Functional analysis of subgroup marker genes was performed using Gene Ontology (GO).

Cluster analysis
Subpopulation trait Genes were screened to identify molecular subtypes by COX-LASSO. The Package “ConsensusClusterPlus” was used to perform cluster analysis to identify nCRT Response subtypes18.

Gene set scores
The IOBR R package was used to assess the immune features and immune cell infiltration in subtypes. The IOBR package in R integrates 6 commonly used algorithms (MCPcounter, TIMER, xCell, CIBERSORT, EPIC, and quanTiseq) to separately analyze tumor-infiltrating immune cells (TILs) in the TME. To assess pathway activity at the sample level, we applied Gene Set Variation Analysis (GSVA) to the normalized gene expression matrix. hallmark pathways, were used to calculate enrichment scores for each sample. Differential pathway activity between groups was evaluated using linear modeling, and multiple comparisons were controlled using the Benjamini–Hochberg false discovery rate (FDR) correction. Pathways with an adjusted P value < 0.05 were considered statistically significant.

Machine learning model development
This study used the mlr3 package in R to train and compare multiple machine learning models to predict the pCR to nCRT in LARC patients. First, several classification models were constructed, including k-nearest neighbors (kNN), linear discriminant analysis (LDA), logistic regression, naive Bayes, random forest (ranger), decision tree (rpart), and support vector machine (SVM).

Visualization
All the visualization were performed using ggplot2 (v3.3.5) (https://github.com/tidyverse/ggplot2), ggpubr(https://github.com/kassambara/ggpubr), scRNAtoolVis (https://github.com/junjunlab/scRNAtoolVis).

Results

Identifying NCRT response associated B cell
We retrieved data from the public database GSE132465 and reanalyzed 23 tumor samples from primary rectal cancer patients, along with 10 matched normal mucosal samples. Based on cell markers, we identified six distinct cell clusters, including T cells (CD3D, CD3E), B cells (MS4A1, CD79A), myeloid cells (LYZ, CD68), mast cells (KIT, CPA3), stromal cells (COL1A1, COL1A2, DCN), and epithelial cells (EPCAM, KRT19) (Supplementary Fig. 1). We extracted 9,146 B cells for further analysis17. All B cells were subjected to further sub-clustering analysis, which identified two main B cell populations: CD20⁺ B cells and CD138⁺ plasma cells (Figs. 1 A, B). CD20⁺ B cells comprise naïve B cells (CD20, CD27, CD38, IGHD), memory B cells (CD20, TNFRSF13B), and germinal center (GC) B cells (CD20, CD27, CD38, MKI67) (Fig. 1C).

By applying the “Scissors” algorithm, we identified 625 B cells as Scissors⁺ cells, which were associated with a favorable nCRT response, and 50 B cells as Scissors⁻ cells, which were associated with an unfavorable nCRT response (Fig. 1D). Most germinal center B cells were identified as Scissors⁺ cells (Fig. 1E). Additionally, previous studies have shown that germinal center B cells are associated with longer patient survival and improved therapeutic responses to nCRT19. Furthermore, we compared the gene expression profiles of Scissors⁺ cells with those of Scissors⁻ cells and identified the differentially expressed genes (DEGs) associated with an effective nCRT response (defined as nCRT-responsive genes) (Supplementary Table 1). Our gene set variation analysis (GSVA) revealed that the gene signature scores of nCRT-responsive genes were significantly higher in nCRT responders than in non-responders (Fig. 1F, G; GSE40492, p = 0.036; GSE45404, p = 0.072), demonstrating that the Scissors algorithm successfully identified B cells associated with the efficacy of nCRT.

nCRT-associated B cells subpopulation trait and interactions within microenvironment
Scissors⁺ B cells showed high expression of genes including ACTG1, CFL1, HLA-DQB1, CD74, HLA-DQA2, and HLA-DRB5. In contrast, Scissors⁻ B cells exhibited high expression of genes such as LMNA, HIST1H2BG, KLF2, IER3, FOS, and BTG2. Next, we investigated the functions of the subgroup marker genes. Marker genes of Scissors⁺ B cells were enriched in pathways including antigen processing and presentation of exogenous antigens, and MHC protein complex assembly. Marker genes of Scissors⁻ B cells were enriched in pathways related to cellular processes, including the response to tumor necrosis factor (Fig. 2A). This indicates that the presence of a higher number of Scissors⁺ B cells in the TME prior to nCRT may enhance the immune response by recognizing additional tumor antigens released by nCRT.

Emerging evidence suggests that cellular crosstalk between tumor cells and non-malignant cells within the TME may be associated with nCRT treatment efficacy. Therefore, the “CellChat” algorithm was used to analyze and construct the potential cellular communication network between various cell populations in the LARC TME. We observed that Scissors⁺ B cells barely emit signals affecting T cells or epithelial cells, whereas Scissors⁻ B cells emit signals targeting both T cell and epithelial cell populations. Further analysis revealed that Scissors⁻ B cells primarily activate the MIF and IGF pathways, which regulate T cells and epithelial cells, respectively (Fig. 2B, C, D, E). Accumulating evidence has indicated that targeting the IGF and MIF signaling pathways can improve the efficacy of cancer treatment20.This may be one of the reasons why Scissors⁻ B cells are associated with resistance to neoadjuvant therapy.

Characterization of B cells developmental trajectories and dynamic changes revealed by pseudotime
To explore the developmental trajectory characteristics of B cells, we used the R package Monocle 2 for pseudotime analysis. We found that Scissors⁺ B cells are localized at the start of the Monocle trajectories, whereas Scissors⁻ B cells are localized at the end (Fig. 3A, B, C). We speculate that within the TME, Scissors⁺ B cells may differentiate into Scissors⁻ B cells. Next, we analyzed the single-cell transcriptomes along the trajectory and observed increased expression of genes such as ISG20, PSME2, RPS5, C1QBP, TMSB10, and HLA-C. In contrast, we observed decreased expression of genes such as SEC61G, IGLC7, SDF2L1, HSP90B1, and PPIB. We extracted the top 200 differentially expressed genes with the most significant changes, which could be clustered into two distinct expression patterns. Group 1 consists of genes with decreasing expression levels along the trajectory. Pathway enrichment analysis revealed that these genes were associated with the aerobic electron transport chain and ATP synthesis-coupled electron transport. Group 2 includes genes with increasing expression levels along the trajectory, and these genes are involved in B cell activation (Fig. 3D).

Metabolism features of B cell across the subsets
Given the significant impact of metabolic reprogramming on immune cell functionality, we investigated the metabolic landscape using the “scMetabolism” package to analyze the metabolic characteristics of the two B cell subtypes. We found that the activities of glycolysis, the TCA cycle, and oxidative phosphorylation were significantly higher in Scissors⁺ B cells compared with Scissors⁻ B cells (Fig. 4A). The expression of genes related to glycolysis and the TCA cycle is upregulated in Scissors⁺ B cells (Fig. 4B, C). Further analysis revealed that the scores for glycolysis, the TCA cycle, and oxidative phosphorylation decreased along the trajectory (Fig. 4D, E, F).

The clinical significance of B-cell-related characteristic genes in the treatment of rectal cancer with nCRT
To clarify the clinical significance of nCRT-responsive genes, 19 genes (Supplementary Table 2) associated with nCRT responsiveness in LARC were further identified from these nCRT-responsive genes via univariate COX-Lasso regression analysis based on the GSE40492 dataset.
Consensus clustering was performed samples from the GSE40492 and GSE45404 datasets, with k = 3 selected based on the CDF, resulting in three stable clusters (Clusters A, B, and C). A higher proportion of patients in Cluster C responded positively to nCRT (76% responders). In contrast, 52% of patients in Cluster A exhibited a response to nCRT, while only 46% of patients in Cluster B benefited from nCRT. Cluster C showed a significantly better nCRT response compared to Cluster B (Chi-square test, p < 0.05; Fig. 5A). Kaplan-Meier analysis showed that patients in Group A had better overall survival (OS) compared to patients in Groups C and B. The five-year survival rate for Group C was 88%, while Group A had a rate of 0.6710 and Group B had a rate of 0.5815. The risk in Group B was significantly higher than that in Group C (HR = 2.93(95Cl% 1.1203–7.648), p = 0.028), while the risk in Group A was higher than in Group C (HR = 2.40(95Cl% 0.8562–6.741), p = 0.096) (Fig. 5B).

These data suggest that the molecular subtyping characteristics of B cell subsets may serve as potential biomarkers for nCRT. The immune infiltration landscape of the three clusters was further assessed. Cluster A exhibited higher fibroblast infiltration. Cluster B showed reduced infiltration of CD8⁺ T cells and cytotoxic lymphocytes. Cluster C demonstrated increased CD8⁺ T cell infiltration (Fig. 5C). PD-L1 expression is well known to correlate with responsiveness to immunotherapy. Higher PD-L1 expression was observed in Cluster A, suggesting that this cluster may benefit from immunotherapy (Fig. 5D). GSVA enrichment analyses were further performed to explore the features of different subtypes. Cluster A exhibited higher scores for hypoxia, epithelial-mesenchymal transition, and angiogenesis. Cluster C had the lowest hypoxia score and the highest oxidative phosphorylation score (Fig. 5E).

Predicting pathological complete response in patients receiving nCRT based on various machine learning models
A predictive model for evaluating nCRT responsiveness in LARC patients was constructed based on nCRT response-related genes. The GSE40492 dataset was used as the training set, and stratified repeated cross-validation was performed using the mlr3 package. Following 3 repetitions of 5-fold cross-validation, the performance of different models was evaluated. The logistic regression method exhibited the best performance in predicting nCRT efficacy in LARC (Fig. 6A). We ultimately chose logistic regression to construct the model. The final regression model showed an AUC of 0.911 in the training set GSE40492 and an AUC of 0.651 in the GSE45404 cohort. We collected Rectal Cancer patients from our hospital who underwent surgery after nCRT as an external validation cohort, achieving an AUC of 0.819 (Fig. 6B, C, D). In our internal cohort of 64 samples, Kaplan–Meier analyses OS and PFS did not reach statistical significance. Nevertheless, the survival curves showed a clear separation, suggesting a potential trend that may not have achieved significance due to the relatively short follow-up period (Supplementary Fig. 2 A, B). In the GSE40492 dataset, Kaplan–Meier analyses were performed to compare survival outcomes between the predicted PCR and non-PCR groups. The predicted PCR group showed markedly better OS and disease-free survival DFS than the non-PCR group. Specifically, the 5-year OS rate was 65% in the non-PCR group and 100% in the PCR group. Similarly, the 5-year DFS rate was 69.2% in the non-PCR group and 100% in the PCR group (Supplementary Fig. 2 C, D).

Discussion
In recent years, a growing body of evidence has supported an association between TME and responsiveness to nCRT21–23. A comprehensive understanding of how the TME influences the mechanism of action of nCRT is crucial for improving the prognosis of rectal cancer. Therefore, bioinformatics approaches were employed to investigate the relationship between B cell subsets and nCRT responsiveness.
In this study, the “Scissor” algorithm was employed to integrate single-cell data with bulk data, identifying two critical B cell subsets associated with nCRT responsiveness. Specifically, the Scissors⁺ subset—characterized by enhanced tumor responsiveness to nCRT—exhibited high expression of antigen presentation genes. Additionally, a significant portion of germinal center B cells was identified as related to nCRT responsiveness. Within the TME, B cells primarily reside within intratumoral tertiary lymphoid structures (TLS). Studies have shown that in mature TLS lacking germinal centers, B cells tend toward immunosuppressive pathways, marked by high expression of genes encoding immunosuppressive cytokines such as IL-10 or TGFβ. In contrast, B cells within mature TLS rich in germinal centers primarily function in antigen recognition and generate disease-related antibodies. These antibodies can mark antigen-expressing cells for clearance through opsonization, complement-mediated lysis, or antibody-dependent cellular cytotoxicity27,28. Based on this, it is speculated that a higher presence of Scissors⁺ subsets before nCRT may facilitate the recognition of tumor antigen exposure induced by radiation, thereby enhancing the body’s anti-tumor immune response.
Over the past decade, numerous studies have demonstrated that immune cell function is closely associated with cellular metabolic patterns26. From a bioenergetic perspective, immune responses against pathogens appear to be quite costly. Changes in B cells—such as proliferation, differentiation, subsequent population contraction, and the eventual formation of memory—resulting from antigen recognition, require a sustained energy supply to match the varying cellular demands across different stages27. Pseudotime inference revealed a trajectory consistent with a possible transition from Scissors⁺ to Scissors⁻ B cells. This transition is accompanied by changes in adenosine triphosphate (ATP) energy metabolism, glycolysis, and the tricarboxylic acid cycle, as well as a reduction in oxidative phosphorylation levels. Such metabolic shifts may be key factors driving the functional transformation of these subsets. these metabolic features require further experimental validation to determine whether they directly drive the proposed phenotypic transition. Currently, it is widely recognized in the field of oncology that cancer cells attempt to suppress immune cell metabolism by competing for nutrients in the crowded and nutrient-restricted TME28,29. Highly invasive tumors often outcompete immune cells in this environment, ultimately modulating immune cell function in a manner that promotes tumor progression. Therefore, targeting these metabolic programs and alterations that govern these changes could uncover new targets for cancer immunotherapy.
Very little is known about how B cells interact with other components in the rectal cancer TME. In the present study, a subset of B cells unresponsive to nCRT was identified, which displays high expression of the IGF signaling pathway and MIF signaling pathway. The IGF pathway primarily affects tumor cells, while the MIF pathway influences T cells30. In a murine rectal cancer xenograft model, the combination of an IGF1R small molecule inhibitor with 5-fluorouracil (5-FU) chemoradiotherapy significantly reduced tumor volume and weight when compared to treatment with nCRT alone, accompanied by an increase in DNA double-strand breaks (DSBs)31MIF can promote the formation of an immunosuppressive microenvironment by inducing the expansion of regulatory T cells (Tregs) or inhibiting the anti-tumor activity of CD8⁺ T cells32,33, Thereby weakening the body’s anti-tumor immune response. These data suggest that the MIF or IGF pathway may be a potential mechanism underlying B cells-mediated resistance to nCRT.
To elucidate the clinical relevance of B cells in rectal cancer treatment, we stratified rectal cancer into three molecular subtypes using signature genes obtained from single-cell data. The subtype exhibiting a favorable response to nCRT showed higher oxidative phosphorylation characteristics and lower hypoxia scores. Conversely, the nCRT-resistant group displayed elevated hypoxia scores. These findings suggest that the hypoxic microenvironment may trigger metabolic reprogramming in B cells, leading to functional alterations. Interestingly, Subtype A exhibited higher levels of angiogenesis, epithelial-mesenchymal transition (EMT), and PD-L1 expression, indicating that Subtype A may be more likely to benefit from immunotherapy.
This study has several limitations. First, the findings are primarily based on computational analyses, including pseudotime inference and transcriptomic profiling, which may introduce model-dependent bias. Second, the study lacks mechanistic investigations to directly confirm the proposed functional transitions between B-cell subsets and the associated metabolic changes. Despite these limitations, the preliminary findings still offer valuable and constructive insights. The identification of a B cell subpopulation associated with response to nCRT provides a potential avenue for personalized treatment in rectal cancer. Patients exhibiting this signature pathways may benefit more from nCRT or immunotherapy. Moreover, integrating these B cell-related molecular signatures into predictive models could aid in patient stratification, treatment selection, and follow-up planning.

Supplementary Information
Below is the link to the electronic supplementary material.