Malignant cell receptor-ligand subtypes guide the prediction of prognosis and personalized immunotherapy of liver cancer

Single-cell RNA clustering

After rigorous quality control, we retained 24,329 cells, comprising eight cell subpopulations (Figure 1A). The number of UMIs and mRNAs showed a significant correlation, while the number of UMI/mRNAs was not significantly correlated with the content of mitochondrial genes (Supplementary Figure 1A). Supplementary Figure 1B, 1C depict the violin plots before and after quality control, and Supplementary Figure 1D demonstrates significant elimination of batch effects between samples.

Figure 1. Single-cell RNA clustering. (A) UMAP plot of the distribution of 10 samples, UMAP plot of the distribution of 20 subpopulations, and UMAP plots of the subpopulations after cell annotation. (B) After CNV comment, the proportion and cell number of the subpopulation in each sample, and the UMAP map of malignant cells. (C) Dot plots of the expression of the first 5 marker genes of the subpopulations after CNV annotation. (D) The differential gene enrichment analysis of each subpopulation of cells.

UMAP dimensionality reduction analysis was performed on 24,329 cells using RunUMAP function, revealing the presence of eight cell types: B cell, Endothelial cells, Hepatocytes, Macrophages, Monocytes, NK cell, Fibroblasts, and T cells. Aneuploid cells were included in cell cluster marker identification for ten cell types. We further analyzed the 10 samples based on the proportion of the ten subpopulations and identified the location of aneuploid malignant cells on the UMAP (Figure 1B). The top-five marker genes with the most outstanding contribution in each subgroup are shown in Figure 1C, and the results of marker genes are provided in Supplementary Table 1. The differential gene enrichment analysis of each subpopulation of cells is shown in Figure 1D, and the aneuploid cells differential genes were mainly enriched in the complement and coagulation cascades, cholesterol metabolism, PPAR signaling pathway, and pentose and glucuronate interconversion.

Cellchat analysis and construction of cellular communication network

After rigorous quality control, we retained 49,701 cells, comprising 22,920 cancer cells and 26,781 normal liver tissue cells. Uniform manifold approximation and projection (UMI) and the number of mRNAs showed significant correlation, while there was no significant correlation between number of UMI/mRNAs and the content of mitochondrial genes, as seen in Supplementary Figure 1A. The violin plot before and after quality control can be seen in Supplementary Figure 1B–1D demonstrates significant elimination of batch effects between samples.

We used copykat to identify malignant cells and analyzed the intersubpopulation communication network using the CellChat (version 1.5.0) package. Results showed that the communication between malignant cells and macrophages had the highest number and intensity (Figure 2A, 2B). Additionally, visual analysis revealed that the first six signaling pathways (MIF, MHC-II, APP, MHC-I) of subgroup communication had signal outputs of malignant cells, with the MIF signaling pathway having the strongest output and Monocyte being the main signal receiver cell (Figure 2C, 2D).

Figure 2. CellChat analysis and construction of cellular communication network. (A) Circle diagrams of that numb and strength of subpopulation interaction. (B) Number of communication interactions of each subgroup. (C) Communication and interaction strength of each subgroup. (D) Aneuploidy cell subpopulations and ligand receptor profiles among cell subpopulations.

Construction of molecular subtypes based on malignant cell receptor ligand genes

To analyze the role of malignant cell-associated receptors in Bulk RNA-seq data, we screened 72 malignant cell-associated receptor genes using the results of the single-cell dataset CellChat. We found that 7 of the 72 related genes were associated with prognosis (p < 0.05), and protective gene 3 and risk gene 4 are shown in Figure 3A. We observed a relatively stable clustering result for Cluster 2 in the cumulative distribution function (CDF) Delta area curve (Figure 3B, 3C) and selected K = 2 to obtain two molecular subtypes (Figure 3D). Furthermore, analysis of the prognostic characteristics of these two molecular subtypes revealed a significant difference in prognosis (Figure 3E), with intersecting clust2 and clust1 subtypes having worse prognosis. The tcga. Subtype. CLI. Txt contains data for the TCGA dataset subtype. We performed the same clustering method on the independent ICGC-JP dataset and obtained similar results (Figure 3F). The icgc. Subtype. CLI table provides data on subtypes in the ICGC-JP dataset.

Figure 3. Construction of molecular subtypes based on malignant cell receptor ligand genes. (A) Single factor Cox analysis HR distribution map of malignant cell associated receptor gene set. (B) CDF curve of TCGA cohort sample. (C) CDF Delta area curve of TCGA cohort sample, the horizontal axis represents the number of clusters K, and the vertical axis represents the relative change of the area under the CDF curve. (D) Sample cluster heat map of consensus K = 2. (E) KM curve of prognosis relationship between two subtypes of TCGA. (F) KM curve of prognosis relationship between two subtypes of ICGC-JP.

Mutation characteristics of molecular subtypes

We obtained SNV mutation data from TCGA using Mutect2, as shown in Figure 4A, and presented the top 15 genes with the most significant mutations in each subtype. Moreover, we surveyed “Homologous Recombination Defects”, “Fraction Altered”, “Number of Segments”, “Nonsilent Mutation Rate”, and “Aneuploidy Score,” as well as the distribution of the Silent Mutation Rate. We noticed that the Homologous Recombination Defects and Aneuploidy scores varied between subtypes, as demonstrated in Figure 4B (PMC5982584).

Figure 4. Genomic alterations in molecular subtypes of the TCGA cohort. (A) Somatic mutation analysis of different molecular subtypes in the TCGA cohort. (B) Comparison of “Homologous Recombination Defects”, “Fraction Altered”, “Number of Segments” “in different molecular subtypes of the TCGA cohort”. “Nonsilent Mutation Rate”, “Aneuploidy Score”, “difference in Silent Mutation Rate”.

Pathway analysis of molecular subtypes

We involved performing differential analysis for TCGA gene expression profiles grouped by cluster subtypes. This analysis revealed that when comparing clust1 to clust2, the expression of 188 genes increased while 1083 genes exhibited decreased expression (Figure 5A). Using the Hallmark gene sets in the Misgbd database as the background set, we performed a difference analysis GSEA (Figure 5B). The analysis results revealed that the Peroxisome and Histidine metabolism were activated in the clust1 subtype, whereas Salmonella infection, Shigellosis, and Epithelial cell signaling in Helicobacter pylori infection were inhibited, as demonstrated in Figure 5C. Our next step utilized GSVA pathway scoring on KEGG pathway enrichment obtained from the Misgbd database. Our analysis showed that 33 pathways were up-regulated and 23 were down-regulated, as demonstrated in Figure 5D.

Figure 5. Pathway analysis of molecular subtypes. (A) Volcano map of differential analysis of genes between the two subtypes. (B) GSEA-GO analysis of differential genes between the two groups. (C) Bubble map of related pathways activated/inhibited in the comparison of the two subtypes. (D) Heat map of KEGG-related pathways with differences between the two subtypes.

Immunological characteristics of molecular subtypes and differences in immunotherapy/chemotherapy

To determine the differences in immune microenvironment among molecular subtypes, we evaluated the infiltration level of immune cells using the expression levels of genes in immune cells in the TCGA cohort. We calculated the scores of immune cells using CIBERSORT, EPIC, Estlmate, MCPcounter, quanTlseg, TIMER, and xCell. Results indicated that the majority of immune cells had noteworthy differences between the two subtypes, with higher scores observed in clust2 (Figure 6A), which indicates that clust2 (with poor prognosis) exhibited a high immune score. We also looked at the expression levels of immune checkpoint genes in the two subtypes and discovered that there were significant variations between the two subtypes, with clust1 having a higher expression level than clust2 (Figure 6B). Patients with high TIDE prediction scores are unlikely to benefit from immunotherapy, indicating a higher probability of immune escape. A Wilcox Test on Figure 6B showed that the TIDE score was highest in the TCGA cohort’s clust2 subtype, suggesting that this subtype has an elevated likelihood of immune escape and is less likely to benefit from immunotherapy.

Figure 6. Immunological characteristics of molecular subtypes and differences in immunotherapy/chemotherapy. (A) The difference of 28 immune cell scores between different molecular subtypes in the TCGA cohort. (B) The difference of ESTIMATE immune infiltration between different molecular subtypes in the TCGA cohort. The immune checkpoint of differential expression between different subgroups in the TCGA cohort. Difference of TIDE analysis results between different groups in TCGA queue.

Construction of risk model

We conducted Cox regression analysis on the differentially expressed genes among the cluster subtypes and identified 285 prognostically relevant genes present in both the TCGA and ICGC datasets. Using Lasso regression, we compressed the 285 genes further in the TCGA dataset by reducing the number of genes present in the risk model. The change trajectory of each independent variable was analyzed to determine the optimal model based on the lambda parameter; the figure (Figure 7A) indicates that with the gradual increase of lambda, the number of independent variable coefficients tending to 0 also increases gradually. We used a 10-fold cross-validation technique to build the model and analyzed the confidence interval for each lambda parameter, as shown in Figure 7B. The optimal lambda value was found to be 0.04711966, and 16 genes were selected as target genes based on this value. Based on these 16 genes related to prognosis, nine genes (LPCAT1, SLC2A1, NEIL3, SRXN1, TRNP1, SLC7A11, STC2, ZNF239, CBX2, CDCA8, EPO, PON1, PBK, PFN2, ACOT12, ADH4) were identified using Lasso regression. We developed a final 16-gene signature with a Cox proportional hazards regression model with the RiskScore calculated as follows:

RiskScore = −0.205 × LPCAT1 −0.068 × SLC2A1 −0.194 × NEIL3 + 0.508 × SRXN1 + 0.117 × TRNP 1 + 0.29 × SLC7A11 + 0.209 × STC2 + 0.087 × ZNF239 + 0.25 × CBX2 + 0.111CDCA8 + 0.332EPO + 0.141PON1 + 0.377PBK + 0.186PFN2 + −0.175ACOT12 + −0.052 × ADH4.

Figure 7. Construction of risk model. (A) Trajectory of each independent variable with lambda; (B) Confidence interval under lambda. (C) KM curve of high and low risk of risk model constructed by 16 genes in TCGA data set. (D) ROC curve of risk model constructed by 16 genes in TCGA data set. (E) KM curve of high and low risk of risk model built by 16 gene in ICGC data set. (F) ROC curve of risk models built by 16 gene in ICGC data set.

We used the TCGA dataset as our training dataset and calculated the risk score of each sample through 16 gene expression levels to analyze the prognostic prediction efficiency classification at 1, 3 and 5 years. The AUC for each year was 0.7. The Riskscore was standardized by zscore, and samples with Riskscore greater than zero after zscore were divided into a high-risk group and a low-risk group. The KM curve was drawn, and we found a significant difference between them (p < 0.05) (Figure 7C, 7D). To validate the accuracy of our model, we verified the ICGC dataset using the same method and found similar results (Figure 7E, 7F).

RiskScore in different clinicopathological features

We analyzed the relationship between RiskScore scores and tumor clinical characteristics by studying the differences in scores between various clinical phenotypes in the TCGA dataset (Figure 8A). The findings revealed that the risk score increased with clinicopathological characteristics such as clinical grade, as indicated in Figure 8A, 8B.

Figure 8. RiskScore in different clinicopathological features. (A) Comparison of clinical phenotypes between RiskScore groups in the TCGA cohort. (B) Differences between RiskScores of different phenotypes in the TCGA cohort (Wilcox. Test, ^*P < 0.05; ^**P < 0.01; ^***P < 0.001; and ^****P < 0.0001).

Relationship between RiskScore and Immunity

To further investigate the relationship between RiskScore and immunity, we used ESTIMATE to predict the immunity score and, employing Spearman’s method, found a significant positive correlation between the immunity score and RiskScore, as depicted in Figure 9A. Then we analyzed the scores of 28 immune cells to determine their association with high and low risk groups, discerning significant differences in the scores of 12 immune cells, as shown in Figure 9B. The results of TIDE prediction score demonstrated that the TIDE score was higher among those with higher risk scores, as indicated in Figure 9C. We also evaluated the expression of checkpoints CTLA4, PD-1, and PD-L1 in the high-risk and low-risk groups, finding that the expression of these immune checkpoint genes was higher in the high-risk group (Wilcox. Test, Figure 9D). We calculated the TRS score using the ssGSEA method based on TRS-associated characteristic genes, discovering that the Risk-High group had a higher TRS score (Wilcox. Test, Figure 9E). Furthermore, we computed the interferon (IFN) score using ssGSEA method and identified the gamma score (Th1/IFNγ score) to be significantly higher in the Risk-high group than in the Risk-low group, as presented in Figure 9F.

Figure 9. Relationship between RiskScore and immunity. (A) Correlation analysis of RiskScore and immune score. (B) Comparison of 28 immune cell scores in high and low risk groups. (C) Correlation analysis of RiskScore for TIDE. (D) Comparison of immune checkpoint expression in high and low risk groups. (E) Comparison of tumor reactivity scores in high and low risk groups. (F) Th1/IFNγ score comparison between high and low risk group.

RiskScore combined with clinicopathological features further improves the prognostic model and survival prediction

We demonstrated RiskScore as the most significant prognostic factor by performing univariate and multivariate Cox regression analysis on RiskScore and clinical characteristics, as shown in Figure 10A, 10B. To estimate the risk assessment and survival probability of patients, we integrated RiskScore with other clinicopathological characteristics to design a nomogram presented in Figure 10C. The results of the model showed that RiskScore had the most substantial impact on survival prediction. Moreover, we employed the Decision curve analysis (DCA) to evaluate the model’s reliability and observed that both RiskScore and Nomogram benefits were significantly higher than the extreme curve compared to other clinicopathological features. Both the nomogram and RiskScore demonstrated significant survival predictive power, as presented in Figure 10D, 10F. Additionally, we assessed the prediction accuracy of the model using the Calibration curve, illustrated in Figure 10E, and observed that the calibration curve of the three calibration points at 1, 3, and 5 years was close to the standard curve. This indicates that the nomogram provided good prediction performance.

Figure 10. RiskScore combined with clinicopathological features further improves the prognostic model and survival prediction. (A) Univariate Cox analysis of RiskScore and clinical characteristics. (B) Multivariate Cox analysis of RiskScore and clinical characteristics. (C) Nomogram model. (D) Calibration curve of nomogram at 1, 3 and 5 years. (E) Decision curve of nomogram. (F) AUC line chart, the horizontal axis is the time unit year, and the vertical axis is the model AUC value.

Validation of the risk score in HCC

To validate the risk-scoring model in HCC, we first conducted a survival analysis to identify SRXN1 with prognostic significance in the TCGA cohort (Figure 11A). Next, we examined one of the previously identified SRXN1 mRNA that constructs predictive model in several tumor cells. The result of the Quantitative reverse transcriptase-polymerase chain reaction (qRT-PCR) reveals SRXN1 mRNA levels were increased significantly in all tumor cell lines (Figure 11B). Finally, we examined the protein level of SRXN1 in cells of THLE-2 and HEP3B. Our Western blot results show that the SRXN1 protein level was increased significantly in the tumor cell line (Figure 11C).

Figure 11. Validation of the risk score in HCC. (A) Survival analysis of TCGA-LIHC based on SRXN1 levels. (B) qRT-PCR images of tumor cell line (HEP3B, HCCLM3, HEPG2 and HUH7). SRXN1 levels were analyzed by qRT-PCR; ^**P < 0.01 versus liver cell line. (C) Representative images of Western blots of the liver cell line (THLE-2) and tumor cell line (HEP3B). SRXN1 levels were determined by Western blot (mean ± SEM; n = 5; ^**P < 0.001 versus THLE-2).

Single cell data downloading and preprocessing

We utilized the Seurat (Version 4.1.1) R package to process single cell data obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) for a reliable liver cancer single cell expression profile dataset (GSE149614) [24]. The sequencing platform used was GPL24676 Illumina NovaSeq 6000 for Homo sapiens. We filtered the single-cell data by setting the criteria that each gene should be expressed in at least three cells and each cell should express at least 250 genes. We further ensured that the gene expressed by each cell was more than 500 and less than 4000, the unique molecular identifier (UMI) of each cell was less than 15000, mitochondrial gene content was less than 10%, ribosome expression was more than 3%, and red blood cell expression was less than 1%, using the Percentage Feature Set function. We obtained a total of 18 single-cell sequencing samples from the GSE149614 dataset, including 10 hepatocellular carcinoma tissues and 8 normal liver tissues.

TCGA database data and ICGC database data download

We retrieved clinical phenotype data of liver cancer from The Cancer Genome Atlas (TCGA) database and removed samples without survival time and survival status. We ensured that all patients had a survival time of more than 0 days. Additionally, we downloaded TCGA expression profile data and obtained 368 tumor samples. Moreover, we downloaded the liver cancer data, expression profile data, and survival data of the ICGC-JP cohort from the International Cancer Genome Consortium (ICGC) database. We removed samples lacking survival time and survival status and ensured that all patients had a survival time of more than 0 days. Eventually, we obtained a total of 232 samples. We calculated snv mutation information of the gene of TCGA-LIHC using Mutect2 software through TCGAbiolinks (version 2.24.3) download [25]. Homologous Recombination Defects, Fraction Altered, Number of Segments, Nonsilent Mutation Rate, Aneuploidy Score, Silent Mutation Rate Data were obtained from a previous study [26].

Single-cell dimensionality reduction and annotation

To normalize the data of the 18 samples, we applied log-normalization separately. Then, we found hypervariable genes using the FindVariable Features function, which identified variable features based on variance stabilization transformation (VST). We scaled all genes using the ScaleData function followed by PCA dimensionality reduction using RunPCA to find anchor points, choosing dim = 30. The cells were clustered through the FindNeighbors and FindClusters functions with a resolution of 0.45 to avoid batch effect of the samples, a Harmony package (version 0.1.0) was used to batch correct the samples. For cell type annotation, we used the HumanPrimary CellAtlasData dataset in the SingleR (version 1.10.0) R package. We also used the tools Copycat [27] (version 1.0.8). Subsequently, marker genes of each sub-population were screened using the FindAllMarkers function with logfc = 0.5 (fold difference) and Minpct = 0.5 (minimum expression ratio of differential genes), with a screening threshold of a corrected p-value less than 0.05.

Cell communication analysis

We used the R package CellChat (version 1.5.0) for the analysis of cell-to-cell communication [28]. The results of CellChat were used to screen out the receptor ligand genes associated with malignant cells.

Consistent clustering for receptor-ligand genes

We performed univariate Cox regression analysis to identify genes associated with prognosis (p < 0.05), from the genes identified in CellChat. Based on the expression profile of these prognostically relevant genes, we classified the patients using the consensus clustering of tumor tissues in the TCGA dataset through ConsensusClusterPlus (version 1.60.0) package in R programming language [29]. We used the HC algorithm with “Pearson” as the metric distance and performed 500 bootstraps with each bootstrap procedure including 80% of the training set of patients. The number of clusters was set from 2 to 10, and the optimal number of clusters was determined by calculating the consistency matrix and the consistency cumulative distribution function. The CDF was used to determine the best classification.

Path analysis and variance analysis

We obtained Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway genes from the msigdb (http://www.gsea-msigdb.org/) database. To calculate pathway scores for samples of different molecular subtypes, we utilized the GSVA package [30]. For differential analysis of the pathway scores between the two subtypes, we used the limma package (version 3.52.2) [31]. Here, we selected pathway visualization with |logFC| >0.15 and P-value < 0.05 criteria. Furthermore, we conducted differential analysis of gene expression by subtype using the limma package and identified differential genes by screening for |logFC| >log2 (1.2) and P-value < 0.05. We then employed gene set enrichment analysis (GSEA) using the Hallmark gene sets in the msigdb database as the background set, and we merged these results with the previous differential analysis to reveal further insights into the pathways and genes involved.

LASSO regression analysis and risk model construction

Genes significantly associated with survival outcomes were identified through univariate Cox regression models to construct a prognostic model. In addition, reliable predictors were selected using LASSO analysis [32], a compression estimator that constructs a penalty function to compress coefficients and set others to zero. This methodology retains the advantages of subset shrinkage and acts as a biased estimator that reduces multicollinearity problems encountered in regression analysis. We conducted LASSO Cox regression using the R software package glmnet (version 4.1-4). Risk scores for each patient in the TCGA and ICGC databases were calculated using the formula: Risk score = Σ coefficient mRNA × expression level mRNA. Finally, we analyzed the correlation between patient risk scores and prognosis.

Immune infiltration and prediction of immunotherapy

To assess the distribution of specific cellular components in the immune microenvironment, we calculated the scores of 28 immune cells using ssGSEA based on characteristic genes obtained in a previous study [33]. Differences among subtypes were assessed using Kruskal test. In addition, we assessed immune cell infiltration using ESTIMATE [34] (version 1.0.13). We used TIDE to evaluate the potential clinical effects of immunotherapy across groups [35]. A previous study [36] showed that the TRS score could predict a patient’s response to immunotherapy. CD8+ T cells in the tumor microenvironment can produce interferon-γ (IFNγ), which can up-regulate PD-1/PD-L1 and IDO1 gene expression [37, 38]. Up-regulated IDO1 expression is positively associated with tumor progression and poor prognosis [39]. We extracted Th1/IFNγ gene signatures from a previous study [32] and calculated IFNγ scores using ssGSEA for each patient.

Cell culture

Here, we purchased normal liver cell line (THLE-2, Shanghai Academy of Life Science), and tumor cell lines, including HEP3B (catalog No. ZQ0024), HEPG2 (catalog No. ZQ0022), HCCLM3 (catalog No. ZQ0023) and HUH7 (catalog No. ZQ0025). Culture dishes with Dulbecco’s Modified Eagle’s Medium (DMEM) containing 10% fetal bovine serum (FBS) were used to inoculate all cell lines to maintain growth at 37°C with5% CO₂. The cell culture medium was also supplemented with 1% penicillin/streptomycin. The culture dishes were purchased from Guangzhou Jet Biofiltration (Guangzhou, China). FBS, and penicillin/streptomycin were purchased from BI (BI, Ridgefield, CT, USA).

RNA extraction and quantitative reverse transcription PCR

Using the total RNA isolation kit (Foregene, Chengdu, China) and the PrimeScript RT kit (Takara Biomedical Technology, Beijing, China), RNA was extracted from the cells and reverse transcribed into cDNA. TB Green Premix Ex Taq II (Takara Biomedical Technology, Beijing, China) and was tested by real-time PCR with 2^−ΔΔCt. GAPDH was the reference gene in the process. The qPCR primers are listed as follows: SRXN1 Forward Primer CAGGGAGGTGACTACTTCTACTC; SRXN1 Reverse Primer CAGGTACACCCTTAGGTCTGA.

Western blot

RIPA lysis buffer was used to extract the protein, and the BCA Protein Assay Kit was used to calculate the protein concentration. 10% SDS-PAGE was used to separate equal amounts of protein, and the resultant membrane, polyvinylidene difluoride, was used. The membrane was incubated with the primary antibody overnight at 4°C after being blocked with 5% BSA in PBST for 1 hour at room temperature. It was then subjected to three PBST washes before being exposed to an appropriate horseradish peroxidase (HRP)-conjugated secondary antibody for one hour at room temperature. With the use of an ECL chemiluminescent reagent, the reaction was observed. Using Image Lab, the blot’s intensity was measured.