Abstract Background Intrahepatic cholangiocarcinoma (iCCA) is a highly heterogeneous and lethal hepatobiliary tumor with few therapeutic strategies. The metabolic reprogramming of tumor cells plays an essential role in the development of tumors, while the metabolic molecular classification of iCCA is largely unknown. Here, we performed an integrated multiomics analysis and metabolic classification to depict differences in metabolic characteristics of iCCA patients, hoping to provide a novel perspective to understand and treat iCCA. Methods We performed integrated multiomics analysis in 116 iCCA samples, including whole‐exome sequencing, bulk RNA‐sequencing and proteome analysis. Based on the non‐negative matrix factorization method and the protein abundance of metabolic genes in human genome‐scale metabolic models, the metabolic subtype of iCCA was determined. Survival and prognostic gene analyses were used to compare overall survival (OS) differences between metabolic subtypes. Cell proliferation analysis, 5‐ethynyl‐2'‐deoxyuridine (EdU) assay, colony formation assay, RNA‐sequencing and Western blotting were performed to investigate the molecular mechanisms of diacylglycerol kinase α (DGKA) in iCCA cells. Results Three metabolic subtypes (S1‐S3) with subtype‐specific biomarkers of iCCA were identified. These metabolic subtypes presented with distinct prognoses, metabolic features, immune microenvironments, and genetic alterations. The S2 subtype with the worst survival showed the activation of some special metabolic processes, immune‐suppressed microenvironment and Kirsten rat sarcoma viral oncogene homolog (KRAS)/AT‐rich interactive domain 1A (ARID1A) mutations. Among the S2 subtype‐specific upregulated proteins, DGKA was further identified as a potential drug target for iCCA, which promoted cell proliferation by enhancing phosphatidic acid (PA) metabolism and activating mitogen‐activated protein kinase (MAPK) signaling. Conclusion Via multiomics analyses, we identified three metabolic subtypes of iCCA, revealing that the S2 subtype exhibited the poorest survival outcomes. We further identified DGKA as a potential target for the S2 subtype. Keywords: diacylglycerol kinase α, intrahepatic cholangiocarcinoma, MAPK signaling, metabolic classification, multiomics analysis, phosphatidic acid metabolism __________________________________________________________________ Abbreviations DGKA diacylglycerol kinase α iCCA intrahepatic cholangiocarcinoma WES whole‐exome sequencing Human‐GEMs human genome‐scale metabolic models EdU 5‐Ethynyl‐2’‐deoxyuridine KRAS Kirsten rat sarcoma viral oncogene homolog ARID1A AT‐rich interactive domain 1A PA phosphatidic acid MAPK mitogen‐activated protein kinase TP53 tumor protein p53 FGFR2 fibroblast growth factor receptor 2 IDH1/2 isocitrate dehydrogenase 1/2 BAP1 BRCA1‐associated protein‐1 OS overall survival QC quality control BWT Burrows‐Wheeler transform GATK Genome Analysis Toolkit SNV single nucleotide variant Vcf Variant call format ANNOVAR annotate variation VAF variant allele frequency MS mass spectrum ITMS Ion trap mass spectrometry iBAQ intensity‐based absolute quantification CV coefficient of variation NMF non‐negative matrix factorization KEGG Kyoto Encyclopedia of Genes and Genomes GSVA gene set variation analysis DFS disease‐free survival HR hazard ratio BH Benjamini‐Hochberg PermFIT permutation‐based feature importance test SVM support vector machines DNN deep neural networks RF random forest AUC area under the curve TCA tricarboxylic acid cycle tINIT task‐driven integrative network inference for tissues ACHR artificially centered hit‐and‐run t‐SNE t‐distributed stochastic neighbor embedding scRNA‐seq single‐cell RNA sequencing scCATCH single‐cell cluster‐based automatic annotation toolkit for cellular heterogeneity TPM transcripts per million CNV copy number variation shRNAs short hairpin RNAs CCK8 cell counting kit‐8 DMSO dimethyl sulfoxide TMA Tissue microarray IHC immunocytochemistry DAB 3,3 N‐diaminobenzidine tetrahydrochloride RIPA radioimmunoprecipitation assay buffer PVDF polyvinylidene fluoride BSA bovine serum albumin ERK1/2 extracellular regulated protein kinases 1/2 MEK1/2 mitogen‐activated protein kinase kinase1/2 GAPDH glyceraldehyde‐3‐phosphate dehydrogenase HRP horseradish peroxidase ECL enhanced chemiluminescence GSEA gene set enrichment analysis HK1 Hexokinase 1 PFKM muscle phosphofructokinase NT5E ecto‐5’‐nucleotidase ENTPD1 ectonucleoside triphosphate diphosphohydrolase 1 ATP adenosine triphosphate AMP adenosine monophosphate PRODH proline dehydrogenase GLS2 glutaminase 2 HEPH hephaestin ALOX5 arachidonate 5‐Lipoxygenase PLA2G4A phospholipase A2, group IVA CNDP1 cytosolic non‐specific dipeptidase 1 GFPT2 glutamine–fructose‐6‐phosphate aminotransferase 2 HIF‐1 hypoxia‐inducible factor‐1 FcεRI Fc epsilon Receptor I LCAT lecithin‐cholesterol acyltransferase SLC3A2 solute carrier family 3 member 2 GPD2 glycerol‐3‐phosphate dehydrogenase 2 ALDH1L2 aldehyde dehydrogenase 1 family member L2 DC dendritic cell Treg regulatory T PLOD1 2‐oxoglutarate 5‐dioxygenase 1 ASPH aspartate beta‐hydroxylase TTN titin MUC16 mucin 16 EPHA2 ephrin type‐A receptor 2 DG phosphorylate diacylglycerol DEGs differentially expressed genes TNBCs triple‐negative breast cancer MPSs metabolic pathway‐based subtypes PD‐1 programmed cell death protein 1 TGF‐β‐Smad4 transforming growth factor beta‐SMAD family member 4 1. BACKGROUND Intrahepatic cholangiocarcinoma (iCCA), a highly lethal cancer, is the second most prevalent primary liver malignancy, comprising approximately 20% of all hepatic tumors [[68]1]. The incidence of iCCA has increased by more than 140% over the past decades [[69]2] and was estimated to rise by up to 10 folds worldwide in the following decades [[70]3]. Due to silent clinical symptoms at an early stage, most patients are diagnosed at advanced stages. With limited treatment strategies, the 5‐year survival rate of iCCA patients remains at approximately 9% [[71]4]. However, the mechanism underlying the pathogenesis of iCCA remains unclear. Recent studies using high‐throughput sequencing technology have revealed the genetic landscape of iCCA, indicating that tumor protein p53 (TP53), KRAS, fibroblast growth factor receptor 2 (FGFR2), isocitrate dehydrogenase 1 (IDH1)/IDH2, and BRCA1‐associated protein 1 (BAP1) mutations (fusion) are the primary driver gene variations [[72]5]. Based on these findings, precise treatments targeting FGFR2 and IDH1/IDH2 have been developed and have shown preliminary clinical efficacy [[73]6, [74]7]. Therefore, multiomics technology may help build more profound insight into the mechanisms underlying iCCA and provide a novel therapeutic strategy. Through integrated multiomics analysis, the results of molecular studies have led to the development of diverse classifications of iCCA. Dong et al. [[75]5] identified four distinct proteomic subgroups with different biological and clinical features, including inflammatory (S1), mesenchymal (S2), metabolic (S3) and differentiated (S4) features. Among these subgroups, patients with S4 showed the best overall survival (OS). Deng et al. [[76]8] classified cholangiocarcinoma, including iCCA and extrahepatic cholangiocarcinoma, into the following three proteomic subtypes: S‐I (metabolism), S‐II (proliferation), and S‐III (stromal), and patients with S‐I showed the best OS. Lin et al. [[77]9] described the immunogenomic traits of iCCA, and Song et al. [[78]10] depicted the molecularly distinct features of two histological subtypes of iCCA. In addition, Cho et al. [[79]11] identified three clinically supported iCCA subtypes (stem‐like, poorly immunogenic, and metabolism) and found that NCT‐501 (an IDH1 family member inhibitor) exhibited synergism with nanoparticle albumin‐bound paclitaxel for the stem‐like subtype in an iCCA organoid model. These classifications provide us with molecular phenotypes to understand prognostic differences among iCCA patients. However, these studies provide little information on metabolic dysregulation and heterogeneity in iCCA patients. Metabolic reprogramming is recognized as a hallmark of malignancy [[80]12]. On the one hand, metabolic reprogramming promotes tumorigenesis by promoting the rapid proliferation, survival, invasion, metastasis, and therapeutic resistance of cancer cells. On the other hand, with the development of tumors, cancer cells acquire more mutations and changes, further enhancing metabolic reprogramming, which in turn accelerates tumor growth and development [[81]13, [82]14]. In addition, accumulating evidence demonstrates that metabolic reprogramming also has a crucial role in regulating the proliferation, differentiation, and function of immune cells, which eventually leads to tumor immune escape [[83]15] and is related to patient prognosis [[84]16]. Nevertheless, little is known about the metabolic differences among iCCA patients with different prognoses. Thus, we performed an integrated multiomics analysis of 116 iCCA and paired adjacent normal tissues and performed metabolic classifications to depict differences in iCCA patient survival, hoping to provide a novel perspective to understand and treat iCCA. 2. MATERIALS AND METHODS 2.1. Clinical summary We collected iCCA samples and the paired adjacent normal tissues from 116 patients who received liver resection without any anticancer treatments prior to surgery and were pathologically diagnosed with iCCA at Zhongshan Hospital (Shanghai, China), and the follow‐up ended in December 2022. We obtained tissue samples and promptly stored them in liquid nitrogen, within 30 min of resection. All collected samples were intratumoral specimens and confirmed by pathologists from Zhongshan Hospital. Additionally, there were two interruptions in perfusion, each lasting for 15 min during the surgical procedure. The sample collection procedure was consistent across all cases. Patients enrolled in this study all signed the informed consent form to allow the use of their data and samples. The medical ethics committee of Zhongshan Hospital reviewed and approved the use of all the human tissues involved in this research (B2021‐611). The baseline clinicopathological characteristics of the 116 patients with iCCA are shown in Supplementary Tables [85]S1‐S2. 2.2. Sequencing Each sample was weighed to approximately 0.2 g, crushed and mixed in liquid nitrogen, and the samples were taken separately for the following analyses: 1. Genomic analysis (for 116 iCCA samples and the paired adjacent normal tissues), wherein whole‐exome sequencing (WES) was performed using the NovaSeq 6000 platform (Illumina, San Diego, CA, USA) with a single sample size of more than 10 G and a sequencing depth >100×. 2. Transcriptomic analysis (for 116 iCCA samples), wherein RNA was extracted for sample quality control, and quality control‐compliant samples were subjected to bulk RNA‐sequencing using the NovaSeq 6000 platform (Illumina), yielding 10 G raw data per sample. 3. Proteomic analysis (for 116 iCCA samples), wherein approximately 0.1 g issue was lysed in urea solution (8 mol/L urea in 0.1 mol/L Tris‐HCl, pH 8.5, containing 1× protease and phosphatase inhibitor), and the whole proteins were extracted by 1 min of sonication (3 s on and 3 s off, amplitude 25%). The peptide segments were isolated by filter‐aided sample preparation enzymes and fractionated by C18 reversed‐phase liquid chromatography (each sample was divided into 3 fractions). Five hundred nanograms of peptides from each fraction were subjected to mass spectrometry (MS) identification using an Orbitrap Fusion Lumos Tribrid mass spectrum spectrometer (Thermo Fisher Scientific, Carlsbad, CA, USA) coupled online to a nanoflow LC system (EASY‐nLC 1200, Thermo Fisher Scientific). The polypeptide samples were dissolved in 10 μL of solvent A (0.1% formic acid in high‐performance liquid chromatography grade water). Then, 5 μL of dissolved polypeptide samples were loaded onto the precolumn at a flow rate of 3 μL/min on the EASY‐nano‐LC chromatography system, and separated on the column at a flow rate of 600 nL/min. The gradient was as follows: 0‐10 min, 5%–11% solvent B; 10‐110 min,11%–25% solvent B; 110‐140 min, 25%–40% solvent B; 140‐141 min, 40%–95% solvent B; and 141‐150 min, 95% solvent B. 2.3. Quality control and quantification of sequencing data Genomic and transcriptomic data were sequenced for fragment quality using FastQC software ([86]https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) to filter out low‐quality reads using the FASTX‐Toolkit (v 0.0.13, [87]http://hannonlab.cshl.edu/fastx_toolkit/). WES data were analyzed as follows: the Burrows‐Wheeler Alignment Tool [[88]17] (v0.7.17, [89]https://bio‐bwa.sourceforge.net/) was used to map all sequenced read segments to the hg38 reference genome, with the main algorithm implemented using either spaced‐seed indexing or Burrows‐Wheeler transform (BWT) techniques. Data pre‐processing, including the removal of repetitive sequences, was performed according to the Best Practices Workflows recommended by the Genome Analysis Toolkit (GATK) [[90]18]. We used GATK (v4.1.5.0, [91]https://gatk.broadinstitute.org/hc/en‐us) to identify single nucleotide variants (SNVs) in tumor tissues and adjacent normal tissues, using the GATK‐recommended af‐only gnomad.hg38.vcf with somatic‐hg38_1000g_pon.hg38.vcf as a reference. Variant call format (Vcf) was further performed to screen for normal population germline variants and to annotate them with annotate variation (ANNOVAR) (v2019‐10‐24, [92]http://www.openbioinformatics.org/annovar). Finally, 18,593 SNVs were obtained by screening variant allele frequency (VAF) > 5% and supporting mutations of reads > 10 in tumor tissues, VAF < 1% and variant reads < 5 in adjacent normal tissues. Tumor cellularity and tumor ploidy were estimated by sequenza‐utils and sequenza, as well as copy number detection (allele specificity and total copy number) and quantification using the adjacent normal tissues as references.