Identify cell types and malignant cell sub-populations in HCC samples
To assess the progression of malignant cells in primary liver cancer patients, we analyzed a single-cell sequencing dataset from hepatocellular carcinoma (HCC) patients (n = 12)(17). We used 10,434 quality-controlled cells for subsequent analysis. After data quality control, we examined nCount_RNA, nFeature_RNA, and percent.mt (Supplementary Fig. 1A). Clustering analysis was performed using the top 20 principal components (Supplementary Fig. 1B) and a resolution of 0.5 based on PCA dimensionality reduction with Seurat. We used t-distributed stochastic neighbor embedding (tSNE) for non-linear dimensionality reduction and marker gene analysis, resulting in the identification of 19 major cell clusters (Supplementary Fig. 1C). Cell types were determined based on marker gene expression and upregulated genes in each cluster. Identified cell types include Myeloid cells (AIF1, RNASE1, C1QB, HLA-DRA), T cells (IL7R, CD3G, CD2, ITM2A, CD3D), NK cells (GNLY, GZMB, KLRD1, KLRF1, B3GNT7), B cells (MS4A1, BANK1, CD79A, TNFRSF13G, BCL11A), Malignant cells (APOA2, ALB, APOA1, AMBP, APOH, TTR), and HSC (RGS5, COL1A1, ACTA2, PDGFRB) (Fig. 1A and 1B). Our cell clusters were identified based on the annotated information in the source article(17) (Fig. 1C and Supplementary table 1). Malignant cells were further divided into 6 clusters based on tSNE and sub-cluster marker gene expression. (Fig. 1D, 1E and Supplementary Fig. 1F).
Heterogeneity of malignant cell sub-populations in hepatocellular carcinoma
To further analyze the impact of malignant cell subtypes on the survival of liver cancer patients. Firstly, we display the cellular proportion of malignant cell subtypes (Fig. 2A). Then, using ssGSEA, we calculate the infiltration level of each subtype in the LIHC sample of TCGA database (Fig. 2B). Subsequently, survival analysis is performed on each malignant cell subtype (Fig. 2C), among which subgroups 3, 5, and 6 exhibit significant differences in survival (P<0.05). High infiltration scores of cells in subgroups 3 and 5 are associated with lower survival rates. Interestingly, subgroup 6, as a subtype of malignant cells, exhibits a different survival trend compared to other malignant cell subtypes. At the seventh year as a survival node, two different trends in survival appear. In the stage before seven years, a higher degree of cellular infiltration in subgroup 6 is associated with better survival. However, in the stage after seven years, survival is worse. This contradicts the conventional understanding of malignant cell subtypes, as not all subtypes have negative effects in the malignant cell subtype. This warrants further exploration and analysis.
Identification of DEGs in subgroup 6 and analysis of interactions between DEGs and malignant markers
To investigate survival analysis contradictions in subgroup 6, we re-clustered its samples using the "ConsensusClusterPlus" R package (Fig. 3A), and selected samples from "Subtype 1" and "Subtype 2" for subsequent analysis. The survival of "Subtype 1" and "Subtype 2" showed opposite trends. In "Subtype 1", a higher cell infiltration score was associated with poor survival, whereas in "Subtype 2", a higher cell infiltration score was associated with better survival. Survival analysis was then performed from the first year, with a half-year interval as the node. The overall survival results in "Subtype 1" were meaningless, but the survival results were meaningful in the 6.5th and 7th years (P < 0.05). In contrast, the survival results were meaningful overall in "Subtype 2" (Fig. 3B). Differential analysis was conducted on samples with significant survival results. 92 samples in "Subtype 1" were analyzed using "limma" with logFC > ± 1 and P.Value < 0.05 as screening conditions to identify DEGs. A total of 320 DEGs were identified (Supplementary table 5), and the expression levels of all DEGs were displayed in the heatmap (Fig. 3C), which were well clustered between high and low expression groups. We wanted to know the interaction between DEGs and malignant cells. We used cytoscape to visualize the PPI analysis results of 320 DEGs and malignant cell marker genes (ALB, APOA1, AMBP, TTR, APOA2, APOH)(17), and used the "cytoHubba" plugin with the MCC method to predict the features in the network nodes and select the TOP30 genes, as shown in Fig. 3D. The same analysis was performed on all samples (77 samples) in "Subtype 2," and 275 DEGs were obtained for subsequent analysis (Supplementary table 6), with results shown in Fig. 3C and 3D. Our results indicate that there are two different categories of survival outcomes in subgroup 6, and the DEGs selected from the two groups are related to malignant markers. We further predicted drugs that may be effective against HCC based on the TOP30 genes.
Five effective drugs for HCC inhibition and Pyridoxine dosing for HP expression in HCC
To predict potential drugs for HCC, we initially screened 30 genes from "Subtype 1" and "Subtype 2" using the "survival" and "survminer" R packages for survival analysis. Kaplan-Meier plots (Supplementary Fig. 2A and 2B). revealed meaningful survival results, and we identified 7 genes with high expression associated with poorer survival in "Subtype 1" and 12 genes with high expression associated with longer survival in "Subtype 2" (Fig. 4A and 4B). Next, we predicted drug-gene interactions using the DGIdb database, with "FDA approved" and "inhibitor" as the primary filters, followed by "FDA approved" or "FDA Not Approved" but as "inhibitor" for subsequent analysis (Supplemetary Table 2 and Supplemetary table 3). We predicted drugs for three genes (MMP9, CXCL8, and CRP) in "Subtype 1" and six genes (CCL2, CYP2C9, CYP3A4, DCN, HP, and SAA1) in "Subtype 2" (Table 2). Finally, we conducted MTT experiments with these drugs (Pyridoxine, Verapamil, Benzbromarone, Streptozocin, Aspirin, Ketoconazole, Rosuvastatin Calcium, Risperidone, Naproxen, Rapamycin, Ilomastat, Marimastat, Pamidronic acid, Ritonavir) on the two cell lines at 24h, 48h, and 72h. The results showed that Rosuvastatin Calcium, Rapamycin, Pyridoxine, Naproxen, and Pamidronic acid effectively inhibited cell viability at all three time points in both cell lines (Fig. 4C).
Table 1
Name of antibody | Type | Company | Cat No | Dilution |
Anti-CCR4 | Rabbit | Bioss | Bs-1168R | 1:1000 |
Anti-LAG3 | Rabbit | Bioss | Bs-2646R | 1:1000 |
Anti-PCDC1 | Rabbit | Immunoway | YN3309 | 1:1000 |
Anti-IFN-γ | Rabbit | Immunoway | YT2279 | 1:1000 |
Anti-CXCL13 | Rabbit | Fine Test | FNab02094 | 1:1000 |
Anti-CXCR3 | Rabbit | Fine Test | FNab09859 | 1:1000 |
Anti-HP | Rabbit | Bioss | Bs-1808R | 1:1000 |
Anti-GAPDH | Rabbit | Beyotime | AF1186 | 1:3000 |
Table 2
Drug-gene interactions information.
Group | Gene | Drug | Interaction_types | FDA approved or not |
Subtype 1 | MMP9 | MARIMASTAT | inhibitor | FDA Not Approved |
| | ILOMASTAT | inhibitor | FDA Not Approved |
| CXCL8 | PAMIDRONIC ACID | / | FDA approved |
| | ASPIRIN | / | FDA approved |
| CRP | ROSUVASTATIN (CALCUIM) | / | FDA approved |
Subtype 2 | CYP3A4 | RITONAVIR | inhibitor | FDA approved |
| | KETOCONAZOLE | inhibitor | FDA approved |
| | VERAPAMIL | inhibitor | FDA approved |
| DCN | SIROLIMUS/RAPAMYCIN | / | FDA approved |
| CYP2C9 | BENZBROMARONE | inhibitor | FDA approved |
| CCL2 | RISPERIDONE | / | FDA approved |
| SAA1 | NAPROXEN | / | FDA approved |
| HP | PYRIDOXINE | / | FDA approved |
| | STREPTOZOCIN | / | FDA approved |
To evaluate the effects of these five drugs on liver cancer cells, we treated Hep3B and Huh-7 cells with varying drug concentrations. LO2 cells served as the normal control. Rosuvastatin Calcium, Rapamycin, and Pyridoxine inhibited cell proliferation at 24h, 48h, and 72h, while Naproxen affected cell viability at 24h and 72h, and Pamidronic acid at 48h and 72h (Fig. 5A). observed decreased tumor weight and volume in mice treated with the drugs (Fig. 5D-5F), indicating the effectiveness of Rosuvastatin Calcium, Rapamycin, Pyridoxine, Naproxen, and Pamidronic acid in inhibiting HCC.
Furthermore, we found that Pyridoxine (5mM) increased HP protein expression in Huh-7 and Hep-3B liver cancer cells but had no significant effect on the normal liver cell line LO2 (Fig. 5G and 5H). HP, as a potential diagnostic biomarker in HCC, warrants further investigation due to limited research in this area. Therefore, we conducted a bioinformatics analysis to explore the role of HP in HCC.
Down-regulation of HP expression in HCC is associated with effects on differentiation of Th2 cells to Th1 cells and IFNG
To assess HP's performance in HCC, we analyzed its expression levels in pan-cancer. HP was found to be highly expressed in LIHC, with lower expression in liver cancer compared to normal liver tissue (Fig. 6A). GEO data supported this finding, showing higher HP expression in normal liver tissue than in HCC (Fig. 6B). We investigated the immune cell infiltration status of HP in HCC using three methods: MCPcounter, Timer, and GSEA. MCPcounter results indicated a negative correlation between T cell and Myeloid dendritic cell infiltration and HP expression (Fig. 6C) (Supplementary table 7). Timer results showed a negative correlation between HP expression in LIHC and the infiltration levels of B cells, CD8 + cells, CD4 + cells, macrophages, neutrophils, and dendritic cells (p < 0.05) (Fig. 5D). Additionally, we explored HP's association with immune-related biological processes (BP) and KEGG pathways. Using GSEA, we assessed the functional enrichment of HP in gene sets "c2.cp.kegg.v7.5.1.entrez.gmt" and "c5.go.bp.v7.5.1.entrez.gmt." Results revealed a concentration of functional genes at the end of the sequence in pathways "KEGG_T_CELL_RECEPTOR_SIGNALING_PATHWAY" and "GOBP_B_CELL_ACTIVATION," with HP showing a negative correlation with these pathways. (Fig. 6E).
To identify relevant pathways associated with HP up-regulation, we divided TCGA HCC samples into two groups based on HP expression (median value) and selected DEGs (LogFC > ± 1 & P.Value < 0.05). Subsequent analysis revealed 361 DEGs (Supplementary table 8). Pathways associated with HP up-regulation included Metabolic pathways, Complement and coagulation cascades, Metabolism of xenobiotics by cytochrome P450, Retinol metabolism, Drug metabolism - cytochrome P450, Drug metabolism - other enzymes, Chemical carcinogenesis - DNA adducts, Bile secretion, and Biosynthesis of cofactors (Fig. 6F) (Supplementary table 4). Our results suggest that HP may act as a potential tumor suppressor gene in HCC but negatively correlates with immune cell infiltration. Previous studies have indicated HP's involvement in regulating Th1 and Th2 cell differentiation(18, 19). However, the relationship between HP and Th1/Th2 cell differentiation in HCC remains unclear, warranting further investigation into the underlying mechanisms.
Differentiation of Th2 cells to Th1 cells in the T helper cell subpopulation in HCC
To investigate the relationship between Th1 and Th2 cells in HCC, we first classified T cell subsets in HCC into 10 clusters (Fig. 7A and 7B). Activated naïve CD4 + T cells can differentiate into different types of helper T cell subsets(20). Among them, CD4 gene expression was up-regulated in T_Sub 2 and T_Sub 5 (Fig. 7C). Clustering analysis on the top 20 principal components of T_Sub 2 and T_Sub 5 successfully divided T_Sub 5 into three subgroups. Th1 cells (CD4, CXCR3, IFNG, CXCL13, LAG3, PDCD1) and Th2 cells (CD4, CCR4) were identified based on signature gene up-regulation (Fig. 7F). However, re-clustering of T_Sub 2 was unsuccessful. Considering the expression patterns of signature genes and other up-regulated genes (TNFRSF18, CCR8, LAYN, IL1R2, BEX3(21), and FOXP3(22))(Supplementary Fig. 3E), T_Sub 2 was speculated to be regulatory T cells. Using the "Monocle2" R package, we constructed the differentiation trajectory within the T_Sub 5 subset. Results indicated early appearance of Th2 cells in the trajectory, while Th1 cells were mainly observed in branches 1 and 2 (Fig. 7E), suggesting differentiation of Th2 cells into Th1 cells within the HCC T helper cell environment.
To investigate the mechanism underlying HP's impact on Th1 and Th2 cell differentiation in HCC, we used the "Seurat" R package's "FindMarkers" tool to identify DEGs in the Th1 and Th2 sub-populations. We obtained 192 up-regulated and 760 down-regulated DEGs (Fig. 8A and Supplementary table 9). Further analysis of KEGG pathways enriched by these DEGs revealed that 25 DEGs were enriched in the "Th1 and Th2 cell differentiation" pathway (P.adjust<0.05) (Fig. 8B) (Supplementary table 10), including CD3G, HLA-DRA, CD3D, IL2RG, IFNGR2, HLA-DPB1, HLA-DRB1, IFNGR1, CD3E, HLA-DPA1, LCK, HLA-DQB1, HLA-DRB5, HLA-DMA, LAT, HLA-DQA1, HLA-DQA2, IFNG, IL2RB, HLA-DOA, JAK3, ZAP70, FOS, and PRKCQ. Protein-protein interaction analysis between HP and the genes enriched in the "Th1 and Th2 cell differentiation" pathway revealed an association between HP and IFNG (Fig. 8C). Using the "corrplot" R package and "Pearson" analysis, we evaluated the correlation between HP and the 25 DEGs. HP showed a negative correlation with 17 genes (Fig. 8D), including HLA-DQA1, HLA-DOA, HLA-DMB, HLA-DPB1, LAT, JAK3, CD3D, IL2RG, IFNG, PRKCQ, ZAP70, CD3E, LCK, CD3G, IL2RB, IFNGR2, and HLA-DQA2. Our results suggest that HP's impact on Th1 and Th2 differentiation in HCC involves IFNG. We further validated our findings through western blot analysis, which showed up-regulation of Th1 markers (IFNG, LAG3, PDCD1, CXCR3) and down-regulation of the Th2 marker (CCR4) in the HP knockdown group compared to the control group (Fig. 8E). These results indicate that down-regulation of HP expression stimulates the expression of Th1 cell-associated markers and may indirectly influence the differentiation of Th2 cells into Th1 cells.