ABSTRACT
-
Background/Aims
Identifying patients with intrahepatic cholangiocarcinoma (ICC) likely to benefit from immunochemotherapy, the new front-line treatment, remains challenging. We aimed to unveil a novel radiotranscriptomic signature that can facilitate treatment response prediction by multi-omics integration and multiscale modelling.
-
Methods
We analyzed bulk, single-cell and spatial transcriptomic data comprising 457 ICC patients to identify an immune-related score (IRS), followed by decoding its spatial immune context. We mapped radiomics profiles onto spatial-specific IRS using machine learning to define a novel radiotranscriptomic signature, followed by multi-scale and multi-cohort validation covering 331 ICC patients. The signature was further explored for the potential therapeutic target from in vitro to in vivo.
-
Results
We revealed a novel 3-gene (PLAUR, CD40LG, and FGFR4) IRS whose down-regulation correlated with better survival and improved sensitivity to immunochemotherapy. We highlighted functional IRS-immune interactions within tumor epithelium, rather than stromal compartment, irrespective of geospatial locations. Machine learning pipeline identified the optimal 3-feature radiotranscriptomic signature that was well-validated by immunohistochemical assays in molecular cohort, exhibited favorable external prognostic validity with C-index over 0.64 in resection cohort, and predicted treatment response with an area under the curve of up to 0.84 in immunochemotherapy cohort. We also showed that anti-uPAR/PLAUR alone or in combination with anti-programmed cell death protein 1 therapy remarkably curbed tumor growth, using in vitro ICC cell lines and in vivo humanized ICC patient-derived xenograft mouse models.
-
Conclusions
This proof-of-concept study sheds light on the spatially-resolved radiotranscriptomic signature to improve patient selection for emerging immunochemotherapy and high-order immunotherapy combinations in ICC.
-
Keywords: Intrahepatic cholangiocarcinoma; Radiogenomics; Multi-omics profiling; Machine learning; Prediction model
Study Highlights
• We integrated multi-omics analyses to establish a spatial-based radiotranscriptomic signature that distinguished immune contexture, prognosis and response to immunochemotherapy for ICC. Targeting uPAR, the kernel of signature, potentiated the efficacy of anti-PD-1 therapy, and may be a promising therapeutic strategy for ICC at predicted high-risk immunochemotherapy resistance. Radiotranscriptomics enables rapid translation of precision medicine for ICC.
Graphical Abstract
INTRODUCTION
Intrahepatic cholangiocarcinoma (ICC) represents the second most common primary liver cancer after hepatocellular carcinoma, with increasing incidence and mortality worldwide [
1,
2]. Only around 20–30% of patients are eligible for curative‐intent resection with 5-year survival rate of 20– 35% [
1]. For the remaining 70–80% of patients with locally advanced or metastatic disease, systemic chemotherapy remains the standard-of-care treatment but has limited benefit, with median survival duration of <12 months. Although cancer immunotherapies have yielded unprecedented results over the last decade, early-phase studies demonstrated the limited activity of immune-checkpoint inhibitor (ICI) monotherapy in patients with advanced-stage ICC, underscoring the need for novel immune-based treatment strategies [
3,
4]. Recent evidence from phase III clinical trials and real-world data suggests programmed cell death protein 1 (PD-1) inhibitor pembrolizumab or programmed death-ligand 1 (PD-L1) inhibitor durvalumab plus chemotherapy as the new front-line standard-of-care treatment for advanced biliary tract cancer, especially ICC [
5,
6]. Nevertheless, identifying patients who are more likely to benefit from such combinations remains a key challenge.
Emerging studies have pinpointed the effect of tumor-immune interactions on tumor growth and progression, immune evasion and resistance to chemotherapy, and limited efficacy of ICI-based combination therapies [
7]. Recent highthroughput transcriptomic and single-cell RNA sequencing (RNA-seq) analyses have shed light on ICC immune profiles; however, bulk and single-cell approaches fail to provide spatial context and cellular composition across distinct regions [
8,
9]. The emergence of GeoMx™ Digital Spatial Profiling (DSP) platform addresses the above challenges based on user-defined regions of interest (ROIs) [
10]. This platform allows the delineation of distinct spatial immune micro-landscapes, which may provide additional prognostic information and opportunities for novel ICI-based combination strategies to improve patient outcome [
10,
11]. Beyond ICC genomic aspects, multi-scale imaging is an integral component of tumor staging and therapeutic response assessment. The distinct biology of the tumor and its interaction with tumor-immune microenvironment (TIME) has prompted distinct imaging phenotypes that can be mined to predict clinical outcomes and genomic hallmarks, known as radiomics and radiogenomics [
12]. However, adaptations to high-dimensional and heterogeneous data need to be addressed when interrogating either radiomics or genomics with the intent to develop prediction models and optimize risk stratification. Fortunately, machine learning that utilizes sophisticated computational modeling to learn from data shines in handling such issues in biomedicine [
13]. Here, we present a multidisciplinary study on multi-scale modelling by integrated analysis of multi-omics data to define immune-related radiotranscriptomic signatures with machine learning. We then assess whether such signatures are complementary to known clinicopathological prognostic factors and predict response to immunochemotherapy in ICC.
MATERIALS AND METHODS
Study design and data sources
Overall study design is presented in
Fig. 1. We collected bulk RNA-seq data and patient-level clinical data of ICC from publicly available datasets with a minimum of 40 tumor samples and follow-up data. A total of 401 patients (244 from Fudan-ICC [
14], 81 from GSE89749 and 76 from EMTAB-6389) were obtained. The immune-related score (IRS) was constructed on the Fudan-ICC cohort, and validated in GSE89749 and E-MTAB-6389 cohorts. Single-cell RNA-seq analysis was performed on 16 human ICC samples and 6 adjacent normal liver samples from two GEO datasets (GSE138709 and GSE189903). We profiled 120 samples from 40 patients with 3 bulk regions each, designated as tumor core (inner one-third portion), intermediate zone (transition one-third portion from the core to the margin) and invasive margin (outer one-third portion), using GeoMx™ DSP technology, to characterize the spatial heterogeneity of tumor-immune interactions and map spatial radiomics onto corresponding genomics. A total of 240 (120 epithelial and 120 stromal) ROIs annotated by a certified pathologist were profiled from all samples.
Three-level validation that included a total of 331 ICC patients from multiple high-volume institutions was purposefully designed as follows: (1) Molecular cohort (156 specimens from 52 patients); (2) Resection cohort (243 patients); (3) Immunochemotherapy cohort (36 patients). Participant inclusion and exclusion criteria as well as follow-up and immunohistochemistry (IHC) protocols are detailed in
Supplementary Methods. Baseline characteristics of multi-scale validation cohorts are summarized in
Supplementary Table 1. This study was approved by the Institutional Ethics Committee of the First Affiliated Hospital of Nanjing Medical University (No. 2024-SRFA-011), and the requirement to obtain informed written consent was waived.
We evaluated immune infiltration profiles based on single-sample gene set enrichment analysis (ssGSEA) to define the resampling-based method termed consensus clustering for immune cluster discovery. Weighted correlation network analysis (WGCNA) was then performed to identify immune-related hub genes. Overlaps between WGCNAdefined immune-related hub genes and immunologically relevant genes downloaded from the ImmPort (https://www.immport.org/shared/home) were included in ranking the importance of prognosis-related genes based on least absolute shrinkage and selection operator (LASSO) and random survival forest (RSF) algorithms, followed by backward stepwise Cox regression analysis to compile the IRS. We validated the robustness of IRS by multi-cohort bioinformatics analysis of bulk and single-cell RNA-seq data. Details are described in
Supplementary Methods.
Prepared slides were incubated with immunofluorescent antibodies for morphology staining: pan-cytokeratin for epithelial cells, CD45 for immune cells, α-SMA for fibroblast compartments, and Syto13 for nuclei. Stained slides were loaded onto Nanostring GeoMx human whole-transcriptome atlas (WTA; 18,677 genes) platform and scanned with a ×20 image of up to four-channel fluorescence. ROIs were targeted to freeform polygon-shaped regions with densest immune signal within tumor epithelium or stroma. Sequencing libraries were generated from photo-released indexing oligos and on Illumina sequencing platform by Genedenovo Biotechnology Co., Ltd. (Guangzhou, China). Details are summarized in
Supplementary Methods.
Technical specifications of contrast-enhanced computed tomography (CT) imaging and tumor segmentation method are detailed in
Supplementary Methods. Tumor sub-regions, which corresponded to tumor core, intermediate zone and invasive margin in tissue sampling, were automatically segregated based on dedicated in-house software with code deposited on the Github website (
https://github.com/DrZhenggangXu/SpatialRT_ICC). Two main protocols for high-throughput feature extraction exist: handcrafted and deep-learning radiomics (
Supplementary Fig. 1). A total of 851 handcrafted features and 2,048 deep-learning features were extracted from each three-dimensional (3D) segmentation; details are provided in
Supplementary Methods. Values of extracted features were standardized with Z-score transformation in each cohort. Stable and reproducible radiomic features, defined as features with intraclass correlation coefficient above 0.90 in both test-retest and inter-reader settings, were retained. Then, features that were highly correlated (correlation coefficient above 0.30) with any component of IRS were prioritized for further study.
The machine learning pipeline contained four steps. First, the study cohort was randomly split into a 70% training/internal validation set and a 30% external testing set. Second, pairwise features with a mutual correlation above 0.8 were removed, retaining the one with higher correlation with the task. Third, 6 machine learning classifiers (support vector machine, random forest, logistic regression, regularized logistic regression, AdaBoost and decision tree) wrapped by recursive feature elimination were employed for feature filtering and prediction modelling. Each classifier was optimized over 5-fold cross-validation to fine-tune hyper-parameters. Finally, the trained model that exhibited the best performance in internal validation set was selected. The machine learning pipeline was integrated into an open-source in-house software shared without access restrictions in a public repository (
https://github.com/salan668/FAE).
Cell proliferation was measured by CCK-8 and clone formation assays. Cell migration and invasion were analyzed by transwell and wound healing assays. We established a patient-derived xenograft (PDX) model in humanized C-NKG mice (NOD-Prkdc
scidIl2rg
em1/Cyagen; 6–8 weeks old). Mice were injected with anti-uPAR monoclonal antibody (mAb) (prepared by Lanzhou Huazhitiancheng Biotechnologies Co. Ltd.; 10 mg/kg) either alone or in combination with anti-PD-1 mAb (Camrelizumab; 10 mg/kg) intraperitoneally every 5 days, as previously described [
15]. Experimental protocols and reagents are detailed in
Supplementary Methods and
Supplementary Table 2. Animal experiments were performed according to the guidelines of the Institutional Animal Use and the Animal Experimentation Ethics Committee of The First Affiliated Hospital of Nanjing Medical University.
Between-group comparisons were analyzed by χ2, Fisher’s exact test, Student’s t-test or Mann–Whitney U-test based on data type and distribution. The Kruskal–Wallis test was performed to compare more than two groups. Correlations between two continuous variables were assessed using Pearson correlation coefficients. Area under the receiver operating characteristic curve (AUC) and its confidence interval (CI) were used to determine prediction performance. The prognostic discrimination was quantified using C-index. Survival curves were estimated by the Kaplan– Meier method and compared with the log-rank test. Tumor response assessment was performed according to Response Evaluation Criteria in Solid Tumors (RECIST) v.1.1. All statistical analyses were performed using R version 4.2.3 and associated packages. Statistical significance was established at P<0.05.
RESULTS
Discovery and multi-level validation of immunerelated score
We found an optimal 3-cluster solution (C1, C2, and C3) that demonstrated significant differences in immune infiltration, with C1 and C3 having a markedly higher immune infiltration abundance than C2, in the discovery (Fudan-ICC) cohort (
Fig. 2A). WGCNA identified four key modules (blue, green, plum2, and skyblue3) that exhibited high correlation with immune clusters, and then 706 immune-related hub genes (
Fig. 2B,
Supplementary Fig. 2A,
2B). With the intersection of ImmPort template, a total of 172 overlapping genes were retained for subsequent analysis. Next, univariate Cox analysis identified 46 prognostic genes that were subjected to LASSO-Cox and RSF algorithms for feature ranking, and top 10 candidate genes from both algorithms were extracted, which identified 5 overlapping genes (
Fig. 2C). We then pruned the multivariate Cox model by backward stepwise selection, yielding a simple 3-gene IRS: risk score=(0.44×PLAUR)+(–0.35×CD40LG)+(–0.25×FGFR4). Notably, urokinase-type plasminogen activator receptor (uPAR/PLAUR) was the only independent risk factor (
Fig. 2D).
IRS indicated favorable prediction of survival with a C-index of 0.72 (95% CI 0.67–0.77) in the Fudan-ICC cohort. Multivariate Cox regression demonstrated that IRS (hazard ratio [HR] 2.31; 95% CI 1.69–3.10;
P<0.001) remained a powerful independent predictor of prognosis after adjusting for available clinical traits (
Supplementary Fig. 3). IRS also stratified external validation cohorts (GSE89749 and EMTAB-6389) with respective C-index of 0.64 (95% CI 0.44–0.79) and 0.65 (95% CI 0.47–0.79). Stratifying patients by the median value demonstrated that score-high patients had a significantly poor prognosis across cohorts (
Fig. 2E). We further investigated its association with immune infiltration estimated by CIBERSORT, which demonstrated positive connection with immunosuppressive barriers (e.g., M0 macrophages and neutrophils) but negative connection with anti-tumor immunity (e.g., naive B cells and CD8 T cells) (
Fig. 2F). IRS was negatively correlated with ImmuCellAI score that predicts comprehensive T-cell subsets, and differed significantly between ImmuCellAIpredicted immunotherapy responders and non-responders (
Supplementary Fig. 4A,
4B). Furthermore, we observed a significant positive association between IRS and hypoxia that modulated tumor proliferation (P53 pathway), metabolism (glycolysis), invasion and metastasis (TNF-α signaling via NF-κB) (
Supplementary Fig. 4C). Drug sensitivity analysis highlighted that low-IRS correlated with remarkably increased sensitivity to current preferred chemotherapeutic (gemcitabine and oxaliplatin) and molecular-targeted (lenvatinib and regorafenib) agents (
Supplementary Fig. 4D,
4E).
Besides, we obtained single-cell transcriptomes for 51,642 cells from GEO database after quality filtering while these cells were partitioned into 10 main cell clusters using known marker genes. We found that IRS was top enriched in monocyte-macrophages, followed by malignant cells; PLAUR, CD40LG and FGFR4 were top enriched in mono-cyte-macrophages, T cells and malignant cells, respectively (
Fig. 3A). Intercellular communication networks reflected intensive cellular interactions between malignant cells, monocyte-macrophages and T cells (
Fig. 3B). We highlighted the high probability values of numerous immuneregulatory interactions between malignant cells, monocytemacrophages and T cells, such as MDK-NCL, SPP1-CD44, and LGALS9-CD44/45 (
Fig. 3C), which have been proved to inhibit T cell activation and limit anti-tumor immune responses [
16-
18]. To provide more detailed mechanistic insights of the 3-gene IRS at single-cell level, we examined the expression of selected markers in distinct sub-populations of enriched cell types. We found that malignant cells that expressed PLAUR were mutually exclusive with FGFR4 expression, suggesting their different roles in ICC progression (
Fig. 3D). Increased communication probability of the above-mentioned immunosuppressive ligand-receptor pairs was noted in PLAUR+ cells compared with FGFR4+ cells (
Supplementary Fig. 5). Moreover, CD40LG was enriched in memory CD4+ T cells that enhance antitumor immunity while PLAUR was enriched in all monocyte-macrophages clusters, characterized by M2-like phenotype and expression of inhibitory molecules (e.g., CD163, MS4A4A, and VSIG4) (
Fig. 3E, 3F) that have been proven to suppress T-cell function [
19].
We harnessed DSP to characterize the immune landscape in distinct niches associated with IRS, which demonstrated that the 3-gene IRS was almost comparable between different intra-tumor regions, but substantially differed between epithelial and stromal regions (
Fig. 4A,
4B). Representatively, FGFR4 was overexpressed in tumor epithelium versus stroma, which was consistent with single-cell RNA-seq results. We next interrogated associations between IRS and immune infiltration calculated by CIBERSORT that were almost identical to bulk RNA-seq results within tumor epithelium, rather than stromal compartment (
Fig. 4C). Thus, epithelium-based DSP data were used to create radiotranscriptomic connections.
We extracted a total of 5,798 radiomic features from each spatial sub-region of the whole tumor that corresponded to each DSP sampling region, on arterial and portal venous phase scans. Among 3,827 features with high stability and reproducibility, 288 (154 handcrafted and 134 deep learning) features that exhibited high correlation with any IRS component, were subjected to the machine learning pipeline. We designed three models defined by feature types (handcrafted versus deep-learning versus combined) for the probabilistic classification of IRS (score-low versus score-high) by using a median split. Although all three models exhibited good performance, using 3 handcrafted features with logistic regression algorithm resulted in the highest accuracy with AUC of 0.95 (95% CI 0.85–1.00) in the internal validation subset and was exported as the radiotranscriptomic signature (
Fig. 4D):
Score=4.19×AP_wavelet-LHL_glcm_MCC–3.07×AP_wavelet-LLL_firstorder_RootMeanSquared–2.97×AP_wavelet-LHH_gldm_DependenceEntropy
We confirmed strong correlations between selected radiomic features, spatial-specific gene expressions and immune infiltration in turn (
Fig. 4E). Strikingly, the 3-feature radiotranscriptomic signature revealed the highest positive correlation with PLAUR expression and macrophage infiltration within tumor epithelium.
Initially, we found a significant correlation between IRS calculated on IHC-based density quantification and the signature (Spearman’s rho 0.49; P＜0.001) in the molecular cohort, consistent with the DSP data (
Fig. 5A). We then observed that the signature stratified patients into distinct prognostic subgroups, and predicted overall and recurrence-free survival with C-index of 0.67 (95% CI 0.63–0.71) and 0.64 (95% CI 0.60–0.69), respectively, in the resection cohort (
Fig. 5B). Stepwise multivariate Cox regression analysis demonstrated that the signature was an independent prognostic factor, even when major semantic imaging features were included (
Supplementary Table 3). Additionally, median progression-free and overall survival times were 6.3 (95% CI 3.4-not applicable) and 16.2 (95% CI 10.3-not applicable) months, respectively, after a median follow-up of 16.2 (interquartile range 6.7–26.1) months for the immunochemotherapy cohort. The objective response rate (ORR) was 36.1% (13/36); the disease control rate (DCR) was 52.8% (19/36). Compared with PD-L1 expression, the signature exhibited significantly better performance metrics (all
P<0.05 by DeLong test) in predicting objective tumor response and disease control, with respective AUC of 0.84 (95% CI 0.69–0.99) and 0.81 (95% CI 0.67–0.96) (
Fig. 5C-D). We also observed statistically significant differences for progression-free survival (HR 6.73; 95% CI 2.58–17.54;
P<0.001) and overall survival (HR 14.35; 95% CI 3.98–51.77;
P<0.001) between the predicted high-risk versus low-risk groups (
Fig. 5E).
Fig. 5F shows examples of the signature to predict response to anti-PD-1/PD-L1 plus chemotherapy.
Above, we highlighted uPAR expression, the kernel of signature, in malignant cells and monocyte-macrophages, and, therefore, designed preclinical
in vitro and
in vivo studies to examine the potential of therapeutic strategies targeting uPAR. Although high uPAR expression in 3 different ICC cell lines was verified, we chose HuCCT1 cells that exhibited relatively higher expression and were interfered by short small interfering RNA (siRNA)-mediated knockdown, and RBE cells that exhibited relatively lower expression and were interfered by plasmid-mediated over-expression (OE) (
Fig. 6A,
Supplementary Fig. 6A). Using CCK8, colony formation, transwell and wound healing assays, we showed that uPAR promoted ICC progression and THP-1 migration. To elucidate its potential mechanisms, we performed RNA-seq of RBE and THP-1 cells before and after uPAR OE, followed by KEGG pathway analysis, which revealed that PI3K-Akt signaling pathway was enriched, also consistent with the results of single-cell data analysis (
Fig. 6B–
6D,
Supplementary Figs. 6,
7). Furthermore, both PI3K inhibitor (LY294002) and anti-uPAR mAb could markedly suppress PI3K-Akt signaling pathway and, subsequently, ICC progression as well as THP-1 migration (
Figs. 6E,
6F,
7A). Anti-uPAR mAb strongly induced the production of THP-1-derived CXCL9 that has been proven to regulate the recruitment of CD8+ T cells [
20].
Besides, using humanized mice bearing uPAR-positive PDX, we observed that tumor growth was markedly impeded by anti-uPAR or anti-PD-1 alone compared with placebo, and further suppressed by the combination of two mAbs. Moreover, IHC staining showed that such combination therapy induced a strong synergistic inhibition of cell proliferation and macrophage infiltration with increased tumor-infiltrating CD8+ T cells (
Fig. 7B–
7F).
DISCUSSION
ICC is a highly heterogeneous disease characterized by a desmoplastic TIME that is poorly immunogenic, and the abundance of immunosuppressive cells that facilitate tumor progression, immune evasion and resistance to chemotherapy [
7]. Here, we integrated bulk, single-cell and spatial transcriptomics to propose an IRS that distinguished immune characteristics, prognosis, and potential benefits from immunochemotherapy for ICC. Following this discovery, we established and validated a novel signature that mapped radiomic profiles onto the spatially heterogeneous IRS using a machine learning-based radiotranscriptomic approach. Using ICC cell lines and humanized PDX mouse models, we demonstrated that uPAR, the kernel of signature, and particularly in combination with PD-1, may be promising therapeutic targets for patients at predicted highrisk immunochemotherapy resistance.
Recent molecular profiling studies on bulk tumor samples have highlighted distinct ICC subtypes with theoretically therapeutic vulnerabilities; however, each core subtype bears a spectrum of marker genes, which limits its future clinical translation [
8,
21]. Through meticulously designed methodology, we identified a simple 3-gene IRS compris-ing PLAUR, CD40LG and FGFR4, and low-IRS indicated significantly improved survival and response to PD-1/PDL1 inhibitors plus chemotherapy. Note that CD40LG and FGFR4 act as protective factors in our IRS. CD40LG, a member of the tumor necrosis factor superfamily, is mainly expressed on key immune cells, confirming its involvement in anti-tumor immune response, while our single-cell data analysis showed that CD40LG expression was enriched in memory CD4+ T cells, consistent with previous reports [
22]. Our bulk and spatial RNA-seq data also highlighted the remarkable positive correlation between CD40LG and activated-memory CD4 T cells, which can induce cytotoxic T lymphocyte anti-cancer responses, suggesting that CD40LG is associated with the response to immunochemotherapy in cancer patients, in concordance with recent studies [
22,
23]. Besides, FGFR4, a member of the FGFR family, plays a pivotal role in bile acid biosynthesis, glucose metabolism as well as tumor invasiveness and angiogenesis [
24]. Note that both lenvatinib and regorafenib are smallmolecule multi-kinase inhibitors, with more potent activity against vascular endothelial growth factor receptors and the FGFR family. The combination of gemcitabine-oxaliplatin-lenvatinib-anti-PD-1 therapy may have additive or synergistic inhibitory effects for patients with advanced ICC, with a promising ORR of 80% and median overall survival of 22.5 months, in a recent report [
25]. Inspired by subtypespecific therapeutic strategies, we therefore presented a classifier that enabled identification of ICC patient subgroup (low-risk), who may be the ideal candidate for immunochemotherapy, even combined with lenvatinib or regorafenib.
On the other hand, PLAUR, also known as uPAR, acts as the key risk factor in our radiotranscriptomic signature, and triggers extracellular matrix proteolysis that has been proved to promote cell migration [
26]. Our results suggested that uPAR activated PI3K-Akt signaling pathway to promote ICC progression and macrophage infiltration, which exerts immunosuppressive effects in ICC [
7]. A recent study on diffuse-type gastric cancer has demonstrated that therapeutic strategies targeting uPAR can potentiate anti-PD-1 efficacy [
15]. Here we examined the prospects of anti-uPAR alone or in combination with PD-1 blockade against ICC using human cancer cell lines and humanized PDX mouse model, which is considered an optimal
in vivo model for immuno-oncology research [
26]. We found that uPAR knockdown or blockade significantly suppressed proliferation and migration of ICC cells
in vitro and
in vivo. Anti-uPAR therapy alone significantly reduced tumor burden in humanized mice bearing uPAR-positive ICC tumors, while anti-uPAR combined with anti-PD-1 therapy further enhanced the antitumor efficacy, indicating their additive effects on immune activation and tumor growth inhibition. Taken together, these results highlight uPAR as a potential therapeutic target for ICC.
Clinical imaging is integral to routine cancer care and exhibits enormous potential to facilitate further advances in oncology practice. Despite significant efforts in decoding biological underpinnings of radiomics, recent studies have so far relied on bulk tissue analysis that obscures the intratumoral spatial heterogeneity, and, to our knowledge, no study has yet attempted to encapsulate genomic or molecular characteristics within radiomics at geospatial level. Note that the radiotranscriptomic signature is composed entirely of arterial-phase wavelet-based features. Reassuringly, these selected radiomics features describe the IRS within tumor epithelium, while increased perfusion induced by ICC cells contributes to hyper-enhanced regions observed in the arterial phase [
27]. We have previously reported strong connections between arterial-phase enhancement patterns and spatial abundance of tumor-infiltrating immune cells [
28]. Accordingly, the application of targeted image filters to isolate spatial patterns from arterial-phase images may be better suited to detect spatial heterogeneity of biological activities for ICC patients. Besides, this study revealed that both handcrafted and deep-learning radiomics offered great potential to identify the genomic subtype of interest, but a combined approach did not exhibit synergistic effects. Moreover, previous studies have shown the potential value of radiomics-based biomarkers in ICI response prediction, yet no study has focused on those predicted non-responder. We demonstrated that targeting uPAR along with PD-1 may be a promising strategy for non-responders predicted by our radiotranscriptomic signature. Although radiotranscriptomics has all ingredients fitted to translational cancer research, it is subject to a unique set of challenges in clinical implementation, such as technical variability, confounding factors and regulatory permission [
12].
The limitations of this study include its retrospective nature that leaves it susceptible to potential bias, and the limited size of training/validation data which can give rise to potential overfitting issues. Further prospective validation with larger sample size is warranted to strengthen our findings. Furthermore, CT images obtained from different scanners and institutions may limit the quality of radiomics data, although we conducted data standardization to minimize the effect of this issue. Additionally, semi-manual tumor delineation is time-consuming and may influence radiomics feature stability, despite our efforts to address this issue by selecting stable and reproducible features.
In conclusion, our study serves as a proof-of-principle for integrating spatial transcriptome with radiomics to describe intra-tumoral TIME heterogeneity and predict treatment response in ICC. Integrated analysis of bulk, single-cell and spatial transcriptomics identifies a novel region-specific IRS that predicts prognosis and immunochemotherapy benefit, and can be decoded by non-invasive imaging using computational and machine learning methods. The final radiotranscriptomic signature has been validated in multiscale independent datasets that may motivate current precision oncology practice in treatment decision-making for ICC patients.
FOOTNOTES
-
Authors’ contribution
Conceptualization: GWJ, XFW, XCL, and XHW; Methodology: GWJ and ZGX; Formal Analysis: GWJ and ZGX; Investigation: SCL, SYC, CYJ, and ML; Resources: CYJ, BZ, YY, QX, and YXX; Writing-Original Draft: GWJ and ZGX; Writing - Review & Editing: KW, XCL, and XHW; Funding Acquisition: GWJ; Supervision: XHW.
-
Acknowledgements
This study was supported by Natural Science Foundation of Jiangsu Province (BK20210968) and National Natural Science Foundation of China (82102150).
-
Conflicts of Interest
The authors have no conflicts to disclose.
SUPPLEMENTAL MATERIAL
Supplementary material is available at Clinical and Molecular Hepatology website (
http://www.e-cmh.org).
Supplementary Figure 2.
WGCNA analysis. (A) Determination of soft-thresholding power (the left panel reveals the influence of softthreshold power on the scale-free topology fit index; the right panel reveals the influence of soft-threshold power on the mean connectivity). (B) Immune-related genes with both high gene significance and module membership in key modules.
Supplementary Figure 4.
Multi-level bioinformatics analysis of IRS. (A) Correlation between IRS and immuCellAI result-based score (n=244). (B) Comparison of IRS between ImmuCellAI-predicted immunotherapy responders and non-responders. (C) Correlation between IRS and pathway activities scored by GSVA. (D) Correlation of risk groups with sensitivity to chemotherapeutics. (E) Correlation of risk groups with sensitivity to targeted therapeutics. IRS, immune-related score; GSVA, gene set variation analysis.
Supplementary Figure 5.
Characterization of malignant cells at the single-cell level. (A) Weighted cell-cell interaction network between PLAUR/FGFR4+ malignant cells and other cell types. (B) Inferred top-ranked ligand-receptor interactions between PLAUR/FGFR4+ malignant cells, monocyte-macrophages and T cells.
Supplementary Figure 6.
uPAR promotes cell proliferation and migration in vitro. (A) The efficiency of uPAR siRNAs in HuCCT1 cells was certified by qRT-PCR. (B) CCK8 assay was performed in HuCCT1 cells transfected with uPAR siRNAs. (C) Colony formation assay of HuCCT1 cells. (D) Transwell assay was performed in HuCCT1 cells transfected with uPAR siRNAs. (E) Wound healing assay was used to evaluate the migration in HuCCT1 cells. Data are shown as mean±standard deviation. *P<0.05, **P<0.01, ***P<0.001.
Supplementary Figure 7.
KEGG pathway analysis to characterize PLAUR in malignant cells and TAMs using single-cell RNA-seq data. TAM, tumor-associated macrophage.
Figure 1.Schematic representation of the study design. (A) IRS identification. (B) Mapping spatial radiomics onto genomics. (C) Multilevel validation of the radiotranscriptomic signature. (D) Experimental investigations to determine the potential of a novel therapeutic target. IRS, immune-related score; RNA-seq, RNA sequencing.
Figure 2.Development and validation of IRS via an integrative procedure. (A) Consensus score matrix of the Fudan-ICC cohort for the optimal 3-cluster solution with the abundance of 28 infiltrating immune cell subsets. (B) Gene co-expression modules identified by WGCNA dendrogram and four key modules that exhibited high correlation with immune clusters. (C) Overlapping genes between WGCNA and ImmPort. filtered by RSF and LASSO-Cox algorithms. (D) Forest plot with hazard ratios of retained genes in the multivariate backward stepwise Cox regression analysis. (E) Kaplan–Meier plots of IRS risk groups in different cohorts. (F) Correlations between IRS and immune infiltration. IRS, immune-related score; LASSO, least absolute shrinkage and selection operator; NS, not significant; RSF, random survival forest; VIMP, variable importance; WGCNA, weighted correlation network analysis. *P<0.05, **P<0.01, ***P<0.001.
Figure 3.Analysis of IRS using single-cell RNA-seq data. (A) t-SNE plots for cell-type identification and all single cells colored by IRS as well as constitutive genes. (B) Weighted cell-cell interaction network between the identified cell types. (C) Inferred ligand-receptor interactions between malignant cells, monocyte-macrophages and T cells. (D) t-SNE plots for PLAUR and FGFR4 expressions in distinct malignant sub-clusters. (E) t-SNE plots for CD40LG expression in distinct T sub-clusters. (F) t-SNE plots for PLAUR expression, M1 and M2 phenotypes as well as violin plots of marker gene expression in distinct monocyte-macrophages sub-clusters. IRS, immune-related score; RNA-seq, RNA sequencing.
Figure 4.Training-validation-testing of spatially-resolved radiotranscriptomic signature using machine learning. (A) DSP of formalin-fixed paraffin-embedded tissue sections harvested from representative spatially separated regions. (B) Violin plots for the spatial distribution of IRS and constitutive genes. (C) Correlation analysis between defined IRS and immune infiltration estimated by CIBERSORT in DSP. (D) Machine learning pipeline with wrapped feature selection based on recursive feature elimination method, selection of the best classifier, weights of radiomics features and receiver operating characteristic curves for the exported models. (E) Correlation analysis between radiomics, IRS and immune infiltration. AUC, area under the curve; DSP, Digital Spatial Profiling; IRS, immune-related score; ROI, regions of interest. *P<0.05, **P<0.01, ***P<0.001.
Figure 5.Multi-level validation of radiotranscriptomic signature. (A) Associations between computationally derived signature and IRS determined by quantitative analysis of IHC slides. (B) Kaplan–Meier plots showing survival of patients following resection. (C) Comparison of treatment response proportions stratified by PD-L1 expression and radiotranscriptomic signature by using χ2 test. (D) Performance of PD-L1 expression and radiotranscriptomic signature in predicting treatment response. (E) Kaplan–Meier plots stratified by the signature. (F) Example implementation of the radiogenomics signature. CAPOX, capecitabine and oxaliplatin; GEMOX, gemcitabine and oxaliplatin; GC, gemcitabine and cisplatin; IHC, immunohistochemistry; IRS, immune-related score; PD-L1, programmed death-ligand 1. *P<0.05, **P<0.01, ***P<0.001.
Figure 6.uPAR promotes cell proliferation and migration by activating PI3K-Akt signaling pathway in vitro. (A) The mRNA and protein expression levels of uPAR in HiBEC and 3 ICC cell lines were measured by qRT-PCR and Western blotting while the efficiency of uPAR plasmid in RBE cells was certified by qRT-PCR. (B) CCK8, colony formation, transwell and wound healing assays were performed in RBE cells transfected with uPAR plasmid. (C) The efficiency of uPAR plasmid or siRNAs in THP-1 cells was certified by qRT-PCR while transwell assay was performed in THP-1 cells transfected with uPAR siRNAs or plasmid. (D) KEGG pathway analysis of differentially expressed genes in RBE and THP-1 cells with uPAR OE and control group. (E) Western blotting analysis showed the levels of PI3K, p- PI3K, AKT and p-AKT in RBE cells before and after PI3K inhibitor (LY294002) or anti-uPAR treatment. (F) CCK8 and Transwell assays were performed in RBE cells and TAMs before and after indicated treatments. Data are shown as mean±standard deviation. ICC, intrahepatic cholangiocarcinoma; OE, over-expression; qRT-PCR, quantitative reverse-transcription polymerase chain reaction; siRNA, small interfering RNA; TAM, tumor-associated macrophage. *P<0.05, **P<0.01, ***P<0.001.
Figure 7.Anti-uPAR alone or in combination with anti-PD-1 treatment in vitro and in vivo. (A) Colony formation and transwell assays of RBE and THP-1 cells treated with anti-uPAR mAb before and after plasmid-mediated uPAR OE with mRNA expression levels of T cellassociated chemokines in THP-1 macrophages before and after anti-uPAR treatment. (B) The uPAR IHC staining of PDX tumor. (C) Schematic diagram of experimental procedure. (D) Tumor images of humanized PDX mice in each group at the end of treatment. (E) IHC staining of Ki-67, CD68+ macrophages and CD8+ T cells. (F) The huPBMC reconstitution efficiency, tumor growth, and tumor weight as well as IHC quantification of Ki-67, CD68+ macrophages and CD8+ T cells at the end of treatment in the respective group. Data are shown as mean±standard deviation (A) or standard error of the mean (F). huPBMC, human peripheral blood mononuclear cell; IHC, immunohistochemistry; OE, over-expression; PD-1, death protein 1; PDX, patient-derived xenograft; *P<0.05, **P<0.01, ***P<0.001.
Abbreviations
intrahepatic cholangiocarcinoma
immune-checkpoint inhibitor
patient-derived xenograft
Response Evaluation Criteria in Solid Tumors
single-sample gene set enrichment analysis
tumor-immune microenvironment
weighted correlation network analysis
REFERENCES
- 1. Moris D, Palta M, Kim C, Allen PJ, Morse MA, Lidsky ME. Advances in the treatment of intrahepatic cholangiocarcinoma: An overview of the current and future therapeutic landscape for clinicians. CA Cancer J Clin 2023;73:198-222.
- 2. EASL-ILCA Clinical Practice Guidelines on the management of intrahepatic cholangiocarcinoma. J Hepatol 2023;79:181-208.
- 3. Harding JJ, Khalil DN, Fabris L, Abou-Alfa GK. Rational development of combination therapies for biliary tract cancers. J Hepatol 2023;78:217-228.
- 4. Yoo C, Hyung J, Chan SL. Recent advances in systemic therapy for advanced intrahepatic cholangiocarcinoma. Liver Cancer 2024;13:119-135.
- 5. Kelley RK, Ueno M, Yoo C, Finn RS, Furuse J, Ren Z, et al. Pembrolizumab in combination with gemcitabine and cisplatin compared with gemcitabine and cisplatin alone for patients with advanced biliary tract cancer (KEYNOTE-966): a randomised, double-blind, placebo-controlled, phase 3 trial. Lancet 2023;401:1853-1865.
- 6. Rimini M, Fornaro L, Lonardi S, Niger M, Lavacchi D, Pressiani T, et al. Durvalumab plus gemcitabine and cisplatin in advanced biliary tract cancer: An early exploratory analysis of real-world data. Liver Int 2023;43:1803-1812.
- 7. Ilyas SI, Affo S, Goyal L, Lamarca A, Sapisochin G, Yang JD, et al. Cholangiocarcinoma - novel biological insights and therapeutic strategies. Nat Rev Clin Oncol 2023;20:470-486.
- 8. Job S, Rapoud D, Dos Santos A, Gonzalez P, Desterke C, Pascal G, et al. Identification of four immune subtypes characterized by distinct composition and functions of tumor microenvironment in intrahepatic cholangiocarcinoma. Hepatology 2020;72:965-981.
- 9. Zhang M, Yang H, Wan L, Wang Z, Wang H, Ge C, et al. Single-cell transcriptomic architecture and intercellular crosstalk of human intrahepatic cholangiocarcinoma. J Hepatol 2020;73:1118-1130.
- 10. Vandereyken K, Sifrim A, Thienpont B, Voet T. Methods and applications for single-cell and spatial multi-omics. Nat Rev Genet 2023;24:494-515.
- 11. Lewis SM, Asselin-Labat ML, Nguyen Q, Berthelet J, Tan X, Wimmer VC, et al. Spatial omics and multiplexed imaging to explore cancer biology. Nat Methods 2021;18:997-1012.
- 12. Bera K, Braman N, Gupta A, Velcheti V, Madabhushi A. Predicting cancer outcomes with radiomics and artificial intelligence in radiology. Nat Rev Clin Oncol 2022;19:132-146.
- 13. Goecks J, Jalili V, Heiser LM, Gray JW. How machine learning will transform biomedicine. Cell 2020;181:92-101.
- 14. Dong L, Lu D, Chen R, Lin Y, Zhu H, Zhang Z, et al. Proteogenomic characterization identifies clinically relevant subgroups of intrahepatic cholangiocarcinoma. Cancer Cell 2022;40:70-87.e15.
- 15. Qin L, Wang L, Zhang J, Zhou H, Yang Z, Wang Y, et al. Therapeutic strategies targeting uPAR potentiate anti-PD-1 efficacy in diffuse-type gastric cancer. Sci Adv 2022;8:eabn3774.
- 16. Li K, Zhang R, Wen F, Zhao Y, Meng F, Li Q, et al. Single-cell dissection of the multicellular ecosystem and molecular features underlying microvascular invasion in HCC. Hepatology 2024;79:1293-1309.
- 17. Ma L, Wang L, Khatib SA, Chang CW, Heinrich S, Dominguez DA, et al. Single-cell atlas of tumor cell evolution in response to therapy in hepatocellular carcinoma and intrahepatic cholangiocarcinoma. J Hepatol 2021;75:1397-1408.
- 18. Peng H, Wu X, Liu S, He M, Tang C, Wen Y, et al. Cellular dynamics in tumour microenvironment along with lung cancer progression underscore spatial and evolutionary heterogeneity of neutrophil. Clin Transl Med 2023;13:e1340.
- 19. Christofides A, Strauss L, Yeo A, Cao C, Charest A, Boussiotis VA. The complex role of tumor-infiltrating macrophages. Nat Immunol 2022;23:1148-1156.
- 20. House IG, Savas P, Lai J, Chen AXY, Oliver AJ, Teo ZL, et al. Macrophage-derived CXCL9 and CXCL10 are required for antitumor immune responses following immune checkpoint blockade. Clin Cancer Res 2020;26:487-504.
- 21. Martin-Serrano MA, Kepecs B, Torres-Martin M, Bramel ER, Haber PK, Merritt E, et al. Novel microenvironment-based classification of intrahepatic cholangiocarcinoma with therapeutic implications. Gut 2023;72:736-748.
- 22. Alexandroff AB, Jackson AM, Paterson T, Haley JL, Ross JA, Longo DL, et al. Role for CD40-CD40 ligand interactions in the immune response to solid tumours. Mol Immunol 2000;37:515-526.
- 23. Lei X, de Groot DC, Welters MJP, de Wit T, Schrama E, van Eenennaam H, et al. CD4(+) T cells produce IFN-I to license cDC1s for induction of cytotoxic T-cell activity in human tumors. Cell Mol Immunol 2024;21:374-392.
- 24. Gadaleta RM, Moschetta A. Dark and bright side of targeting fibroblast growth factor receptor 4 in the liver. J Hepatol 2021;75:1440-1451.
- 25. Shi GM, Huang XY, Wu D, Sun HC, Liang F, Ji Y, et al. Toripalimab combined with lenvatinib and GEMOX is a promising regimen as first-line treatment for advanced intrahepatic cholangiocarcinoma: a single-center, single-arm, phase 2 study. Signal Transduct Target Ther 2023;8:106.
- 26. Chuprin J, Buettner H, Seedhom MO, Greiner DL, Keck JG, Ishikawa F, et al. Humanized mouse models for immunooncology research. Nat Rev Clin Oncol 2023;20:192-206.
- 27. Joo I, Lee JM, Yoon JH. Imaging diagnosis of intrahepatic and perihilar cholangiocarcinoma: Recent advances and challenges. Radiology 2018;288:7-13.
- 28. Ji GW, Xu Q, Jiao CY, Lu M, Xu ZG, Zhang B, et al. Translating imaging traits of mass-forming intrahepatic cholangiocarcinoma into the clinic: From prognostic to therapeutic insights. JHEP Rep 2023;5:100839.
