Identification of Crucial Cancer Stem Cell Genes Linked to Immune Cell Infiltration and Survival in Hepatocellular Carcinoma

Abstract

Hepatocellular carcinoma is characterized by high recurrence rates and poor prognosis. Cancer stem cells contribute to tumor heterogeneity, treatment resistance, and recurrence. This study aims to identify key genes associated with stemness and immune cell infiltration in HCC. We analyzed RNA sequencing data from The Cancer Genome Atlas to calculate mRNA expression-based stemness index in HCC. A weighted gene co-expression network analysis was performed to identify stemness-related gene modules. A single-sample gene set enrichment analysis was used to evaluate immune cell infiltration. Key genes were validated using RT-qPCR. The mRNAsi was significantly higher in HCC tissues compared to adjacent normal tissues and correlated with poor overall survival. WGCNA and subsequent analyses identified 10 key genes, including minichromosome maintenance complex component 2, cell division cycle 6, forkhead box M1, NIMA-related kinase 2, Holliday junction recognition protein, DNA topoisomerase II alpha, denticleless E3 ubiquitin protein ligase homolog, maternal embryonic leucine zipper kinase, protein regulator of cytokinesis 1, and kinesin family member C1, associated with stemness and low immune cell infiltration. These genes were significantly upregulated in HCC tissues. A functional enrichment analysis revealed their involvement in cell cycle regulation. This study identified 10 key genes related to stemness and immune cell infiltration in HCC. These genes, primarily involved in cell cycle regulation, may serve as potential targets for developing more effective treatments to reduce HCC recurrence and improve patient outcomes.

Keywords: cancer stem cells; cell cycle genes; hepatocellular carcinoma; mRNA expression-based stemness index; tumor microenvironment.

Figure 1

Relationships between the mRNAsi and the clinicopathological characteristics and prognosis of HCC patients. (A) mRNAsi score and (B) EREG-mRNAsi score in HCC group and adjacent normal group. (C) Relationship between mRNAsi and clinical stage. (D–F) Relationship between mRNAsi and tumor status. (G) Kaplan–Meier analysis of the relationship between mRNAsi and overall survival of HCC. mRNAsi: mRNA expression-based stemness index; EREG-mRNAsi: epigenetic regulation based-index.

Figure 2

Summary of the TMB information. (A) Missense mutation was the most common variant classification and had the highest frequency. (B) SNP occurred most frequently in variant types. (C) C > T accounted for the most fraction in SNV. (D) The number of tumor mutation burdens in specific samples. (E) The top 10 mutated genes in HCC. (F) Relationship between mRNAsi and TMB. (G) Landscape of mutation profiles in HCC. SNP: single-nucleotide polymorphism; INS: insertion; DEL: deletion; SNV: single-nucleotide variant.

Figure 3

The mRNAsi index associated with WGCNA of HCC. (A,B) A volcano plot and heat map of DEGs in HCC; green indicates downregulated genes and red indicates upregulated genes. |Log2FoldChange| >3, p < 0.05. (C,D) Determination of soft threshold for the similarity matrix. The scale-free correlation coefficient and the mean connectivity for soft threshold powers were analyzed. The number represents the power value and the horizontal axis represents the soft threshold power = 6. (E,F) Gene clustering and gene module partition results. The different branches of the cluster dendrogram correspond to different gene modules that are represented by different colors. (G,H) The correlation between the gene modules and mRNAsi. The Pearson correlation coefficient of the gene module and the traits was plotted as a heat map. (I) The intersection of the DEGs and mRNAsi-related WGCNA-derived genes. DEG: differential expression genes; WGCNA: weighted gene co-expression network analysis.

Figure 4

The immune cell infiltration level was inferred using ssGSEA. (A) Derived immune cell infiltration scores using 28 immune cell types. (B) Classifying samples into two groups (cluster 1 and cluster 2) using unsupervised learning. (C–E) Differences in immune score, stromal score, and ESTIMATE score between clusters 1 and 2.

Figure 5

Differences in 25 key genes expression between different clusters. A total of 10 out of 25 key genes were significantly upregulated in the low group.

Figure 6

Functional enrichment analysis of 10 key genes in HCC. (A) The clustered heat map in HCC and adjacent normal tissue. (B) PPI networks were performed for the key genes by STRING. (C) RT-qPCR was performed to measure the expression of 10 key genes. The enrichment analysis of the 10 key genes of (D) KEGG pathway, (E) biological process, (F) cellular component, and (G) molecular function. * indicated a significance of p < 0.05.