Oncostatin M expression and TP53 mutation status regulate tumor- infiltration of immune cells and survival outcomes in cholangiocarcinoma

In this study, we used bioinformatics tools to analyze transcriptome data from cholangiocarcinoma (CCA) patients in multiple datasets (Sun Yat-sen University, TCGA and GSE32225 cohorts) to identify mechanisms that regulate tumor infiltration by immune cells and survival outcomes. We identified 96 differentially expressed genes (DEGs), including 13 upregulated and 83 downregulated genes, in CCA tissues as regulatory T cells were significantly higher and the proportions of activated natural killer cells and monocytes were significantly lower in CCA tissues than the precancerous tissues. The survival outcomes of CCA patients were associated with the TP53 gene mutation status, levels of Oncostatin M (OSM) expression, and the proportions of tumor-infiltrating immune cell types, including dendritic cells, monocytes, and T follicular helper cells. Functional enrichment analysis of the DEGs in the high OSM-expressing CCA tissues showed that pathways related to tumor progression and immune response were significantly upregulated. Our study demonstrates that OSM expression and TP53 mutation status regulate the tumor infiltration by immune cells and survival outcomes in CCA. OSM is thus a potential prognostic biomarker and therapeutic target in cholangiocarcinoma.


INTRODUCTION
Cholangiocarcinoma (CCA) is a highly malignant tumor originating from the intra-or extra-hepatic bile ducts that has shown increasing morbidity and mortality rates worldwide in the last few years [1]. The primary treatment option for early-stage CCA patients is surgical resection, but the 5-year overall survival (OS) rate is below 10% because majority of the CCA patients are diagnosed in the advanced stages and are not amenable for surgery [2,3]. Chemo-radio therapeutic outcomes are poor in CCA patients because the tumor is highly desmoplastic with fibrogenic connective tissue and immune cells such as T lymphocytes, natural killer (NK) cells and macrophages that infiltrate the tumor epithelium [4]. Hence, there is an urgent need to find effective targeted treatments to improve the clinical outcomes of CCA.
Advances in immunotherapy have shed greater focus on the role of tumor-infiltrating immune cells, which are vital components of the tumor immune microenvironment (TIME). CCA patients with a higher proportion of neutrophils and T-regulatory cells (Tregs) and lower proportion of CD8 + T cells in the tumor tissues are associated with poor prognosis [5]. However, only limited types of tumor-infiltrating immune cells have been analyzed in CCA tissues and the mechanisms regulating the tumor infiltration of immune cells are poorly understood.

AGING
Recent studies demonstrate that somatic mutations in the tumor tissues influence immunotherapeutic response in several cancers [6,7]. Somatic mutations can reduce or abolish the ability of immune cells to recognize neoantigens on the tumor cells [8]. Some studies have reported that somatic mutations can influence immunotherapeutic outcomes [6,7]. TP53 is the most frequently mutated gene in more than 50% of all human cancers [9]. TP53 mutations are also associated with the infiltration of immune cells into the tumor microenvironment [10][11][12]. In lung cancer patients, TP53 gene mutation status is associated with prognosis and therapeutic outcomes [10]. TP53 expression is also a potential diagnostic biomarker in CCA patients [13].
Oncostatin M (OSM) is a cytokine secreted by differentiated histiocytic lymphoma cells [14]. OSM regulates the secretion of cytokines such as IL-6, G-CSF and GM-CSF from the endothelial cells [15][16][17]. OSM also inhibits growth and proliferation of several types of tumor cells such as A375 melanoma [18]. Higher levels of OSM in the early stages of early gastric cancer and breast cancer regulate tumor progression [19,20]. Altered OSM expression is also associated with 22q11-q13 somatic mutations in CCA [21]. However, the mechanism by which OSM expression regulates the growth and progression of CCA has not been reported.
In this study, we used bioinformatics tools to analyze the transcriptome data from multiple CCA patient datasets (Sun Yat-sen University, TCGA and GSE32225 cohorts) and determine the relationship between OSM expression, tumor infiltration of immune cell types, TP53 gene mutational status and survival outcomes in CCA.

Identification of differentially expressed genes in three CCA patient datasets
The study workflow diagram is shown in Figure 1.

Identify the differentially expressed genes of three datasets
We first analyzed the RNA-seq data from various CCA patient datasets using the edgeR software to identify differentially expressed genes (DEGs) in CCA tissues. We identified 3625 DEGs (2063 upregulated and 1562 downregulated genes) in the TCGA CCA dataset, 1122 DEGs (455 upregulated and 667 downregulated genes) in the GSE32225 dataset, and 851 DEGs (252 upregulated and 599 downregulated genes) in the Sun Yat-Sen University CCA patient dataset (Figure 2A, Supplementary Figure 1). The volcano plots and heatmaps for all the DEGs are shown in Figure 2A.

Inflammatory response and complement pathway genes are enriched in CCA tissues
We identified 96 overlapping DEGs between the three CCA patient datasets, including 13 upregulated and 83 downregulated genes using BEG online tools ( Figure  2B). Functional and pathway enrichment analysis of the DEGS using Metascape showed that genes related to complement and coagulation pathway as well as acute inflammatory response were significantly enriched in the CCA patient tissues ( Figure 2C, 2D).

Landscape of tumor-infiltrating immune cells is altered in CCA tissues
We used the CIBERSORT algorithm to analyze the proportions of 22 types of tumor-infiltrating immune cells in the CCA and precancerous tissues. As shown in Figure 3A, 3B, the proportions of CD4 + memory T cells (p=0.045) and T regulatory cells (Tregs; p=0.001) were significantly increased and the proportions of activated NK cells (p=0.032) and monocytes (p=0.008) were significantly reduced in the CCA tissues. Furthermore, high proportions of T regulatory cells had a better prognosis (HR=0.473, p=0.024), however, high proportions of NK cells had a poorer prognosis (HR=2.63, p=0.011).
Sia et al divided the CCA patient samples in the GSE32225 dataset into inflammation (n=57, 38%) and proliferation (n=92, 62%) groups [22]. Therefore, we analyzed the differences in the proportions of tumorinfiltrating immune cells between these two groups. The results showed that proportions of 14 out of 22 types of tumor-infiltrating immune cells analyzed were significantly different between the inflammation and proliferation groups of CCA patient tissues ( Figure 3D, 3E). The proportions of CD4 + naïve T cells (p=0.014), follicular helper T cells (p=0.015), resting NK cells (p=0.001), M0 macrophages (p=0.001), M1 macrophages (p=0.002), activated myeloid dendritic cells (p=0.002), and activated mast cells (p<0.001) were significantly higher and the proportions of memory CD4 + resting cells (p<0.001), activated memory CD4 + cells (p=0.032) cells, T regulatory cells (p=0.014), activated NK cells (p<0.001), monocytes (p<0.001), M2 macrophages (p=0.005), and resting mast cells (p=0.005) were significantly lower in the inflammation group compared to the proliferation group of CCA patients. Furthermore, we generated visual plots of the relative proportions of 22 different types of immune cells in the tumor tissues from inflammation and proliferation groups of CCA patients using the 'ggplot2' R package and observed clear differences between the two groups ( Figure 3F).

TP53 gene mutations are associated with survival outcomes in CCA patients
We next explored the correlation between somatic mutations and tumor-infiltration of immune cells in the CCA tissues. We identified 20 frequently mutated genes in the ICGC cohort of 173 Japanese CCA patients ( Figure 4A) and 30 frequently mutated genes in the TCGA cohort of American CCA patients ( Figure 4B). The gene mutations were mostly missense type and were distributed on all the chromosomes ( Figure 4C). We analyzed the top30 mutant genes from all the three datasets and identified TP53, PIK3CA, BAP1, ARID1A as the most frequently mutated genes in the CCA tissues ( Figure 4E). Kaplan-Meier survival curve analysis showed that overall survival was shorter for CCA patients with TP53 mutations compared to those with wild-type TP53 ( Figure 4F). The other 3 genes, PIK3CA(p=0.385), BAP1(p=0.323), ARID1A(p=0.164) did not show prognostic significance. Hence, we focused on evaluating the correlation between TP53 mutations and the tumor-infiltration of immune cells.

Prognostic significance of OSM expression in CCA Patients with and without TP53 gene mutations
We analyzed 4 CCA patient samples with TP53 mutations from 51 TCGA CCA patient samples with somatic mutations. and identified 879 DEGs, including 728 upregulated and 151 downregulated genes (p≤0.05) in the TP53 mutant CCA patients ( Figure 5A). PPI network analysis of DEGs using the STRING database showed that 133 DEGs were enriched for the GO term, immune system process (GO:0002376, p=0.007; Figure  5B). We then calculated the nodal scores using Cytoscape and identified 10 immune-related genes out of the 133 DEGs for further analysis ( Figure 5C). Kaplan-Meier survival curve analysis and log-rank test showed that overall survival was significantly shorter for CCA patients with low OSM expression AGING (OSM low group) compared to those in the OSM high group ( Figure 5D). The remaining 9 genes did not show prognostic significance. Hence, we selected OSM for further analysis.

OSM expression is an independent predictor of prognosis in CCA patients
Western blot assays in 12 paired CCA tumor and adjacent non-tumor liver tissues showed OSM expression in CCA tumor tissues is lower than non-tumor liver tissues ( Figure 6A). Then we analyzed the correlation between clinicopathological characteristics and the levels of OSM expression in CCA patients. Towards this, we analyzed the OSM expression in the Sun Yat-Sen university cohort of 203 CCA patients by IHC and divided the patients into high-and low-OSM expression groups (n=101 and n=102, respectively) using a median IHC expression score of 4 ( Figure 6B). Kaplan-Meier survival curve analyses and log-rank test showed that the overall survival (OS) and disease-free survival (DFS) of the OSM low group CCA patients was significantly shorter compared to the CCA patients from the OSM high group (OS: p<0.001; DFS: p = 0.006; Figure 6D). These results demonstrate that high OSM expression indicates better prognosis in CCA. We then performed the Wilcox test, chi-square test and logistic regression analysis to evaluate the correlation between the levels of OSM expression and the clinicopathological characteristics of CCA patients. As shown in Figure 6C and Table 1 Multivariate Cox regression analysis showed that tumor numbers (HR=1.611, p=0.020), tumor encapsulation (HR = 0.626, p=0.018), CA199 levels (HR = 1.519, p=0.032) and the levels of OSM expression (HR = 0.541, p=0.001) were independent predictors of prognosis in CCA patients (Table 2). Furthermore, tumor numbers (HR=1.665, p=0.012), CA125 levels (HR = 1.595, p=0.014), and tumor stage (HR =1.770, P=0.028) were independent predictors of poor RFS ( Table 2).

OSM expression correlates with tumor infiltration of immune cells in CCA tissues
We then analyzed the correlation between the levels of OSM expression and the 14 types of tumor-infiltrating immune cells that are significantly different between the proliferation and inflammation groups of CCA patients. TIMER database analysis showed that OSM expression was positively associated with the proportions of memory B cells, resting CD4 + T cells, M2 macrophages, resting mast cells, monocytes, and activated myeloid dendritic cells (p<0.05; Figure 7A). We also analyzed the correlation between OSM expression and the proportion of cancer-associated fibroblasts and endothelial cells, which are known to regulate the tumor immune landscape. We observed no significant relationship between OSM expression and the proportion of cancer-associated fibroblasts (p=0.103) and endothelial cells (p=0.174) (Supplementary Figure 2).
We then analyzed the correlation between OSM expression and the immune cell markers in the CCA tissues relative to adjacent non-tumor liver tissues using the TIMER (Table 3) and GEPIA (Table 4) databases, respectively. The data were then adjusted for purity. The results showed a strong correlation between the levels of OSM expression and levels of markers related to DCs, monocytes, and Tfh cells in the CCA tissues (Tables 3, 4). Overall, these data suggest that OSM expression regulates the infiltration of immune cells such as DCs, monocytes, and Tfh cells into the tumor microenvironment in CCA tissues.

OSM expression correlates with three immune regulatory checkpoints in CCA tissues
We then investigated the correlation between four immune regulatory checkpoints PDL1, CTLA4, HAVCR2, LAG3 and the levels of OSM expression in the CCA tissues by GEPIA. The correlation analysis method is spearman. Results showed that the expression of OSM in the CCA patient tissues showed strong correlation with CTLA4 (cor=0.33, p=0.049), HAVCR2 (cor=0.45, p=0.006), PDL1 (cor=0.43, p=0.009) in

Functional enrichment analysis of OSM-related genes in CCA tissues
We then performed functional enrichment analyses of OSM-related DEGs using data from OSM high (n=18) and OSM low (n=18) patients from the TCGA dataset. KEGG pathway analysis showed that pathways related to tumor progression and immune response, such as JAK-STAT signaling pathway, B-cell receptor (BCR), T-cell receptor (TCR), and TOLL-like receptor (TLR) signaling pathways, cytokine receptor signaling pathways, and natural killer (NK) cell-mediated cytotoxicity were significantly enriched in the OSM high group, whereas, pathways related to Huntington's and  Parkinson's disease were significantly enriched in the OSM low group ( Figure 8A). Furthermore, GO terms such as inflammatory response, regulation of neutrophil migration, innate immune response activating cell surface receptor signaling pathway, cell membrane and immunological synapse, GABA receptor activity, immunoglobulin binding, chemokine activity and cytokine receptor binding were significantly enriched in the OSM high group ( Figure 8B; Table 6).

DISCUSSION
Improved understanding of the immunobiology of cancer tissues and the tumor microenvironment (TIME) has led to the emergence of highly effective anti-cancer therapies related to immune checkpoint blockers (ICB) such as PD1, PDL1, and CTLA-4 [23]. Currently, several studies are focused on identifying mechanisms that regulate tumor infiltration of immune cells because   [24]. The tumor-infiltrating immune cells show a strong association with the prognosis of lung, bladder, and pancreatic cancer patients [10,25,26]. Several types of immune cells, such as lymphocytes, macrophages, neutrophils and natural killer cells, have been reported to regulate cholangiocarcinoma genesis [27]. Elevated levels of pre-operative peripheral blood neutrophil count relative to lymphocyte counts are associated with poor prognosis of intra-hepatic and extra-hepatic CAA patients [28]. CCA is a highly desmoplastic tumor with abundant amounts of various immune cell types, including tumor-associated macrophages (TAMs) and myeloid-derived suppressor cells or MDSCs [3]. TAMs interact with the cancer stem cell niche and modulate adhesion and invasive properties of tumor cells that are associated with tumor progression [29]. Our study demonstrates that the TAM numbers are significantly higher in the CCA tissues compared to the precancerous tissues.

AGING
The efficacy of ICB therapies is associated with the proportions of tumor-infiltrating immune cells. CCA patients with higher proportions of tumor-resident CD8 + T cells respond better to ICB therapies [3].  CD8 + -tumor-infiltrating lymphocytes (TILs) and the expression of human leukocyte antigen (HLA) class I molecules [10]. Furthermore, upregulation of PD-1 and PD-L1 in the tumor cells is associated with increased tumor invasiveness, poor prognosis, worse disease and metastasis-free survival, and lower infiltration of CD3 +and CD8 + -TILs [11][12][13]18]. In contrast, low expression of PD-L1 correlates with favorable prognosis in CCA patients [19]. This suggests that lymphocytic apoptosis induced by the PD-L1/PD-1 pathway promotes CCA progression.
The chemokines and other immune cell activation factors play an important role in the recruitment of tumor-infiltrating immune cells [23]. Our study demonstrates that the overall survival of extrahepatic AGING cholangiocarcinoma (eCCA) patients significantly correlates with the proportions of DCs, neutrophils, and CD8 + T cells in the tumor tissues. Hence, the proportions of these three types of immune cells in the CCA tissues can potentially predict the prognosis of CCA patients. Moreover, enhanced infiltration of CD4 + and CD8 + T cells is associated with better overall survival, decreased lymph node metastases, and reduced venous and perineural invasion in CCA patients [5,30,31]. The neutrophils, CD8 + T cells, Tregs, and M2 macrophages also regulate inflammation during tumorigenesis [5]. However, the mechanisms regulating the differential infiltration of immune cell types in CCA are not clear. A recent study suggests that CCA cells increase the levels of TGF-β in the tumor tissues by activating Tregs, thereby inhibiting the immune response against CCA tumor cells [32].
Our study demonstrates significant differences in the proportions of activated CD4 + memory T cells, Tregs, resting NK cells and monocytes between the CCA and precancerous liver tissues. Studies have reported that CD4 + T cells inhibit tumor growth by secreting cytokines [33], and tumor-infiltrating CD4 + T cells are associated with increased overall survival of CCA patients [34]. Increased tumor-infiltration of Tregs is associated with worse DFS rates in the CCA patients [35]. Moreover, DFS positively correlates with the proportions of tumor-infiltrating neutrophils and tumorassociated macrophages in the CCA patients [5]. Preclinical data suggests that NK cells inhibit tumor growth by inducing apoptosis of CCA cells [36,37]. Infiltration of TIE2-expressing monocytes (TEMs) combined with high Ang1 expression positively correlates with overall survival of hilar cholangiocarcinoma patients [38]. Furthermore, the proportions of different types of tumor-infiltrating immune cells vary significantly in the proliferation and the inflammation groups of CCA patients from the GSE32225 dataset [22].
Somatic gene mutations in the tumor cells affect the efficacy of immunotherapy [39]. Somatic mutations in the EP300 gene activate the immune-regulatory signaling pathways and improve antitumor immunotherapeutic outcomes in bladder cancer patients [39]. We studied the correlation between somatic mutations in CAA tissues and tumor-infiltrating immune cells using the somatic mutation data for CCA tissues from the TCGA, cBioportal and ICGC databases, and observed that TP53 gene mutations correlate with overall survival of CCA patients ( Figure  4). Functional enrichment analyses showed that TP53 mutant CCA patients were significantly enriched in pathways regulating immune cell functions. In lung cancer, genes related to TP53 gene mutations are associated with tumor infiltration of immune cells [10]. We demonstrate that TP53 mutation status correlates with the levels of OSM expression and the proportions of tumor-infiltrating immune cells in CCA tissues (Figure 7).
OSM is a cell growth regulating polypeptide that was first isolated from the serum-free supernatants of U937 histiocytic lymphoma cells [18]. OSM inhibits the proliferation and differentiation of liver cancer stem cells (LCSCs) and increases the sensitivity of liver cancer cells to the chemotherapeutic agent, 5fluorouracil [40]. The role of OSM in cancers is controversial. OSM promotes epithelial to mesenchymal transition (EMT) of cancer stem cell (CSC) properties by activating Stat3/TGF-β/SMAD3 signaling, and therefore, promotes tumor metastasis and tumor recurrence [41]. OSM maintains hematopoietic progenitor cells in the bone marrow by regulating G-CSF and SDF-1 cytokine levels [14]. Bone marrow hematopoeisis is significantly reduced in the OSMdeficient mice because of elevated circulating levels of G-CSF [14]. Neutralizing anti-GM-CSF antibody reduces the levels of OSM secretion by the neutrophils, but anti-G-CSF antibody has no effects; moreover, coculturing neutrophils with GM-CSF increases OSM secretion, but, co-culturing with G-CSF has no effects on OSM secretion [42]. OSM inhibits the proliferation of breast cancer cells and promotes their detachment and motility [42,43]. IL-6 and IL-6 related cytokines promote OSM-dependent suppression of epithelialspecific genes and enhance epithelial cell death by phosphorylating STAT5 and STAT3 [42,43]. OSM reduces the growth of human melanoma xenograft tumors in a mouse model in an IL-6-dependent manner [44]. In contrast, OSM and IL-6 are expressed in breast, prostate, and lung cancer cell lines, and promote tumor development by modulating lipid metabolism, matrix degradation, angiogenesis, and cellular dedifferentiation [45,46]. In early stages of gastric cancer, higher OSM expression correlates with poor prognosis [19]. Furthermore, OSM plays an essential role in the early stages of breast cancer metastasis [20].
Our study has several limitations. Firstly, we did not have clinical data for the CCA patients from the GEO database. Furthermore, complete clinical data was available only for 30 CCA patients in the TCGA dataset, whereas, for other samples, only survival time and status were available. Therefore, we analyzed the correlation of OSM expression and prognosis of CCA with only limited number of clinicopathological characteristics. Secondly, we performed immunohistochemical analysis of the CCA patient samples from the Sun Yat-Sen University cohort to validate OSM expression and then investigated the AGING correlation between clinicopathological characteristics and OSM expression based on those findings. However, further experiments such as RT-qPCR and functional assays are necessary to verify the conclusions of our study.
In conclusion, we demonstrate that tumor infiltration of immune cells is associated with OSM expression and the status of TP53 gene mutations in CCA patients. We also demonstrate that OSM expression is an independent predictor of prognosis in CCA patients.
Our study suggests that OSM expression is a novel prognostic biomarker and therapeutic target for CCA.

CCA patient datasets
We enrolled 208 consecutive CCA patients that underwent hepatectomy between January 2007 and June 2016 at the First Affiliated Hospital of Sun Yatsen University, Guangzhou, China. The inclusion criteria were: (1) patients diagnosed with CCA based on histology; (2) patients that underwent curative R0 resection of the primary tumor; and (3) clinicopathological information and follow-up data were available [47,48]. Patients that received neoadjuvant chemotherapy or radiotherapy before hepatectomy were excluded. Therefore, we enrolled 203 out of 208 CCA patients for further analysis of clinicopathological characteristics and OSM immunohistochemistry. We also recruited 5 out of 208 patients between January 2019 and October 2019 for RNA high-throughput sequencing and 12 for western blotting. This study was approved by the Ethics Committee of the First Affiliated Hospital of Sun Yatsen University. All recruited participants volunteered to participate in this study and signed an informed consent form before enrollment [49]. The clinicopathological characteristics of 17 participants are shown in the Supplementary Table 2.
We obtained transcriptome data from the TCGA database that included 9 precancerous and 36 CCA tissue samples and the Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/) database (GSE32225 dataset with 6 precancerous and 149 CCA tissue samples). The CCA samples from the GSE32225 dataset were classified into inflammation and proliferation groups and used for analyzing immune cell infiltration. We also downloaded somatic mutation data for 51 American CCA patient samples from the TCGA portal, 173 Japanese CCA patient samples from the International Cancer Genome Consortium (ICGC) database (https://dcc.icgc.org/), and 361 CCA samples from the cBioportal database (https://www.cbioportal.org/). We obtained available clinical information of the CCA patients from the TCGA and ICGC databases for further analysis.

RNA high-throughput sequencing
The RNA high-throughput sequencing of five matched primary tumor and precancerous tissues of CCA patients was performed using the Illumina HumanRef-8 WG-DASL v3.0 expression beadchip microarray platform. The sequencing libraries were generated with the NEBNext® UltraTM RNA Library Prep Kit for Illumina® (NEB, USA). We obtained 150 bp paired-end reads by sequencing the library preparations on an Illumina Novaseq platform. We calculated fragments per kilobase of transcript per million fragments mapped (FPKM) for the whole 58736 genes.

Western blotting
Total tissue protein lysates were prepared from the snap-frozen CCA patient samples using the RIPA buffer (PC101, EpiZyme, shanghai, China). The protein concentrations were determined using the BCA protein assay. Then, equal amounts of protein samples were separated on SDS-PAGE and transferred onto PVDF membranes. Then, the membranes were blocked with 5% skimmed milk for 1 h at room temperature. Then, the membranes were incubated with primary mouse anti-human OSM antibody (NBP1-47904, Novus Biologicals, USA, 1:1000) and anti-GAPDH antibody (Cell Signaling Technology, Danvers, MA, USA, 1:1000) overnight at 4°C. Then, the membranes were incubated with the horseradish peroxidase (HRP)-conjugated secondary antibodies (Cell Signaling Technology, Danvers, MA, USA, 1:5000) at room temperature. Then, the blots were developed with the ECL system (Thermo Fisher Scientific, MA, USA).

RNA-seq data processing
We used the R software package (v.3.5.3 and v.3.6.1; http://www.rproject.org) to analyze the RNA transcriptome and somatic nucleotide variation data. We identified differentially expressing genes (DEGs) between the precancerous and CCA tissues as well as the wild-type TP53 vs. mutant TP53 groups of CCA patients using the edgeR package (3.18.1). The P values for the differences in gene expression between precancerous and CCA tissues were adjusted using the Benjamini and Hochberg method, and genes with |logFC| ≥1 and adjusted p-value <0.05 were filtered and designated as differentially expressed. Furthermore, the heat maps of the DEGs were visualized using the pheatmap R package, and the volcano plots were generated using the tools in the online sangerbox (http://sangerbox.com/Tool). The overlapping DEGs were visualized with the Venn diagram generated using the Bioinformatics and Evolutionary Genomics website (http://bioinformatics.psb.ugent.be/webtools/Venn/).

Estimation of tumor infiltrating immune cells
We used the TIMER (http://timer.cistrome.org/) database [51] and CIBERSORT algorithm [24,52] to estimate the proportions of 22 different types of tumorinfiltrating immune cells in 5 pairs of CCA samples from our CCA patient dataset and the GSE32225 dataset. The proportions of 22 immune cell subtypes in CCA tissues, correlations, and the principal component analysis (PCA) were visualized using the Vioplot, corrplot and ggplot2 R packages, respectively.

Analysis of somatic mutations
The single nucleotide variants (SNVs) in genes from the CCA tissues in the TCGA dataset were analyzed with the VarScan R packages [49]. The SNVs for the CCA tissues from the TCGA and ICGC datasets were visualized using maftools [53] and GenVisR [54] packages, respectively. The overlapping genes with somatic mutations were analyzed using the Kaplan-Meier survival curves and the log-rank test at the GEPIA website (http://gepia.cancer-pku.cn/).

Functional enrichment analysis
The altered pathways and biological functions in the CCA tissues were determined by performing functional enrichment analysis of the overlapping DEGs in the Sun Yat-Sen University CCA patient cohort, and the GSE32225 and TCGA datasets using the Metascape (https://metascape.org/) database [55]. The functional enrichment analysis (GO and KEGG pathways) of the overlapped DEGs between the wild-type TP53 and mutant TP53 containing CCA patient tissues were performed using the Broad Institute GSEA software 4.0. P <0.05 was used as the threshold to identify enriched GO terms in the cellular component (CC), molecular function (MF), and biological process (BP) categories. The results of the KEGG pathway and GO analyses were visualized using the clusterProfiler [56] and ggplot2 R packages, respectively.

Protein-protein interaction (PPI) network construction and co-expression module analyses
We used the String version 11.0 software (http://stringdb.org/) to generate PPI networks of the DEGs [57]. The co-expression modules were visualized using the Cytoscape v.3.7.2 software [58]. The important modules and the top10 nodes in the PPI network were selected by calculating the node scores.

Statistical analysis
Statistical analyses were performed with the R software (v.3.5.3 and v.3.6.1) and IBM SPSS statistical software (version.25). The Kaplan-Meier survival curves and the log-rank tests were performed using the survival R package.
The correlations between the clinicopathological characteristics and the OSM expression levels were analyzed using the chi-square test, Wilcox test, univariate and multivariate Cox regression analyses. P < 0.05 was considered statistically significant.

AUTHOR CONTRIBUTIONS
CW and LQ contributed to research design; LQ, LT and SYX performed data extraction and statistical analysis; CW, CJP and LQ drafted the manuscript. All authors consented to the final draft of the manuscript.