mRNA expression- and DNA methylation-based stemness indices in HCC
Firstly, to explore the relationship between the index distributions of stem cells and clinical characteristics. We ranked the HCC samples according to their mRNA gene expression-based stemness index (mRNAsi) and DNA methylation-based stemness index (mDNAsi) values and tested whether any clinical feature was correlated with low or high stemness index. The mRNAsi and mDNAsi exponential distribution of TCGA HCC samples were shown in Fig. 1a and Fig. 1b, respectively. The mRNAsi exponential distribution of ICGC HCC samples was shown in Fig. 1c.
Relationship Between Stem Cell Index And Clinical Characteristics
To further explore the relationship between mRNAsi and mDNAsi indexes and clinical features of tumor samples, such as: gender, age, pathological stage, whether HBV or HCV infection, immune cell score,TNM stage and iCluster typing, we performed a linear regression test on the above characteristics. As shown in Table 1, the mRNAsi of TCGA samples were significantly related to pathological stage, HBV infection, immune cell score, and iCluster typing (P < 0.05). And the mDNAsi was remarkably associated with immune cell score and iCluster typing (P < 0.05). Furthermore, we analyzed the data from ICGC cohort, which indicated that the mRNAsi was significantly associated with TNM stage (Table 2, P < 0.05).
Table 1
Correlation between mRNAsi and mDNAsi indices of TCGA liver cancer samples and their characteristic information and molecular typing
| Percent/Average | p(mRNAsi) | p(mDNAsi) |
Sex | Male:65.6% | 0.80 | 0.12 |
Female:34.4% |
Age | 60.13 | 0.47 | 0.23 |
Grade | G1:17.8% | 3.71E-05 | 0.67 |
G2:48.3% |
G3:33.9% |
HBV | 23.00% | 3.04E-04 | 0.47 |
HCV | 18.60% | 0.76 | 0.15 |
Immune score | 166.12 | 1.2E-03 | 0.018 |
iCluster | iCluster1:35.5% | 1.29E-07 | 5.64E-31 |
iCluster2:30.1% |
iCluster3:34.4% |
Table 2
Correlation between the mRNAsi index of ICGC liver cancer samples and the characteristic information and molecular typing of the samples
| Percent/Average | P (mRNAsi) |
Sex | Male:75.4% | 0.12 |
Female:24.6% |
Age | 67 | 0.2 |
TNMStage | I:16.3% | 9.09E-03 |
II:47.3% |
III:29.1% |
IV:7.4% |
HBV | 23.00% | 0.81 |
HCV | 18.60% | 0.64 |
Immunescore | 166.12 | 0.97 |
Table 3
HR analysis and optimal threshold for TCGA stem cell index grouping
iCluster | mRNASi |
hazard.ratio | Cut P |
Hazard ratio | lower | upper | P | Cut | P (Cut) |
C1 | 1.80 | 0.96 | 3.40 | 6.77E-02 | 0.37 | 0.23 |
C2 | 3.61 | 1.20 | 10.84 | 2.21E-02 | 0.35 | 0.08 |
C3 | 1.60 | 0.79 | 3.23 | 0.19 | 0.43 | 0.77 |
All | 1.68 | 1.11 | 2.53 | 0.0145 | 0.38 | 0.08 |
iCluster | mDNAsi |
Hazard ratio | Cut P |
Hazard ratio | lower | upper | P | Cut | P (Cut) |
C1 | 0.69 | 0.37 | 1.31 | 0.26 | 0.12 | 0.87 |
C2 | 0.52 | 0.20 | 1.34 | 1.75E-01 | 0.20 | 0.67 |
C3 | 1.85 | 0.91 | 3.76 | 0.09 | 0.30 | 0.48 |
All | 0.62 | 0.39 | 0.98 | 4.04E-02 | 0.12 | 0.51 |
In TCGA, the mRNAsi level was significantly higher in HBV-infected samples compared with uninfected samples (Fig. 2a, P < 0.05), and the mRNAsi score of iCluster 3 samples was significantly higher than that of iCluster 1 and iCluster 2 (Fig. 2b, c, P < 0.05). In addition, the patients with late pathological stage have the higher mRNAsi expression. (Fig. 2d, e, P < 0.05). The iCluster typing was significantly correlated with mDNAsi. The mDNAsi score was lowest in iCluster1 and highest in iCluster3.While the HBV-infected and pathological stage had no obviously correlation with mDNAsi. (Supplementary Fig. 1a-e, P < 0.05). In ICGC, TNM was significantly related to the mRNAsi, while HBV-infected subgroup exhibited no significantly relation with the mRNAsi. The samples with higher mRNAsi index probably in higher TNM stage (Fig. 2f, g, P < 0.05). The immune score was negatively related to mRNAsi and mDNAsi in TCGA (Fig. 2h, Supplementary Fig. 1f, P = 0.0012). Furthermore, the correlation between the mRNAsi and mDNAsi indexes of tumor samples were validated and indicating that these two subtypes of indexes were significantly positively correlated (cor = 0.30, Fig. 2i, P = 0.00012). Taken together, the stemness indices were negatively related with clinical characteristics, samples with worse clinical characteristics tended to have higher stemness indices, demonstrating that the stem cell indexes could play essential roles in tumor progression.
Correlations Of Stemness Indices With OS In HCC
To better understand the relationship between stemness indices and prognosis, we compared differences of OS with clinical features such as gender, age, pathological stage, iCluster typing, immune cell scores, and HBV or HCV infection among various samples. The results illustrated that clinical factors such as gender, age had no significant correlation with OS (Fig. 3a, b, P = 0.26, P = 0.6). And the patients over survival presented no obviously relation with HCV infection, pathological stage (Fig. 3c, d, P = 0.42, P = 0.25) and iCluster typing (Fig. 3e, P = 0.22).There were significant differences in OS of samples grouped based on HBV infection and immune cell scoring. The high-immune score group demonstrated higher overall survival rate compared with low-immune score group (Fig. 3f, P = 0.012). And the survival probability of group of HBV infected was lower than the group of not infected (Fig. 3g, P = 0.036). The OS of mRNAsi-low group was higher than the mRNAsi-high group, and the mDNAsi-low group was higher than the mDNAsi-high group (Fig. 3h, Supplementary Fig. 2a, P = 0.014). Then the samples grouped based on the stemness indices mRNAsi and mDNAsi were identified according to iCluster typing. Only in the iCluster2 sample group, there was a significant difference of OS in samples grouped by the mRNAsi (Fig. 3i, j, k, Supplementary Fig. 2b, c, d, P = 0.064, P = 0.015, P = 0.18). Shorter and fewer deaths made it difficult to distinguish.
Then we performed further verification in the ICGC database indicating that the OS had no correlation with gender, age (Fig. 4a, b, P = 0.077, P = 0.2), and also had no relevant with immune cell scores, HBV infection (Fig. 4c, d, P = 0.15, P = 0.88). The overall survival rate of HCV infection subgroup was lower and the higher TNM stage had lower OS (Fig. 4e, f, P = 0.033, P < 0.05). And the OS of mRNAsi-low group was higher than mRNAsi-high group (Fig. 4g, P < 0.001). All these hinted, the stemness indices negatively related to OS and affected the prognosis of patients, which may provide a new objective criterion for evaluating the prognosis of HCC patients.
Screening The Stem Cell-related Genes In HCC
Based on the above results, we further analyzed the differentially expressed genes (DEGs). In TCGA cohort, a total of 136 differential genes were filtered according to the threshold, among which the mRNAsi-high group was down-regulated by 110 genes and up-regulated by 26 genes compared to the mRNAsi-low group. The DEGs were mainly down-regulated differentially expression (Fig. 5a). A total of 569 differential genes were selected after filtering, of which mDNAsi-high group were down-regulated by 546 genes and up-regulated by 23 genes compared to mDNAsi-low group. The DEGs were mainly down-regulated differentially expression (Fig. 5b). Among the differential genes in the two sets, mRNAsi and mDNAsi were both down-regulated genes with intersect into four genes, respectively APOBEC3C, C1orf116, GABBR2 and TFCP2L1. In ICGC cohort, the volcano map of DEGs showed there were 190 differential genes. Among them, mRNAsi-high group have 16 genes down-regulated and 174 genes up-regulated compared to mRNAsi-low group, which were major up-regulated differential expression (Fig. 5c). The intersection of the up-regulated genes with the TCGA samples based on the mDNAsi grouping was ALDH1A3. These up- or down-regulated differential genes reflect different gene expressions in HCC samples with high and low stem cell indices, which provides new ideas for the further exploration of HCC pathological mechanisms.
Functional Analysis Of Differentially Expressed Genes
To further explore the function of differentially expressed genes and their roles in tumor genesis, metastasis, recurrence, and drug resistance, a functional analysis of them were performed. GO function annotations of 136 differential genes of TCGA tumor samples grouped by mRNAsi showed without significant enrichment (FDR < 0.05). GO function annotations of 569 differential genes of TCGA tumor samples identified by mDNAsi indicating that the pathway “metal ion trans-membrane transporter activity” and “monovalent inorganic action trans-membrane transporter activity” were significantly enriched (Fig. 5d, FDR < 0.05). KEGG pathway enrichment analysis of mRNAsi differential genes found no significant enriched pathway (FDR < 0.05). For the KEGG pathway enrichment analysis, our results demonstrated the pathway “Camp signaling pathway” was significantly enriched (Fig. 5e, FDR < 0.05). For the GO function annotations of 569 differential genes obtained by ICGC tumor samples based on mRNAsi grouping, the pathway “extracellular matrix structural constituent” and “glycosaminoglycan binding” were significantly enriched (Fig. 5f, FDR < 0.05). The KEGG pathway enrichment analysis of mRNAsi differential genes found three significantly enriched pathways. “PI3K-Akt signaling pathway” was significantly enriched (Fig. 5g, FDR < 0.05). Taken together, we found the DEGs of high and low stemness indices always enriched on pathway correlated with tumor formation, metastasis and recurrence, which help to understand the tumor genesis mechanism and guide the targeted therapy of tumors.
PPI analysis and hub genes analysis of differentially expressed genes
To further find out hub genes related to stemness indices, PPI analysis of differentially expressed genes were performed. 10 genes (SNAP25, EPCAM, CALB2, SOX2, KRT19, GABBR1, MUC1, AFP, GRIN1, and GAD1) were selected as hub genes (Fig. 6a). The correlation between these hub genes and the mDNAsi index were performed. The results suggested that SNAP25, KPT19, GABBR1 and EPCAM showed a significant negative correlation with mDNAsi (Fig. 6b-e, P < 0.0001, P = 0.0023, P = 0.039, P = 0.0006), while the other six hub genes had no prominent connection (Fig. 6f-k, P = 0.062, P < 0.039, P = 0.7, P = 0.21, P = 0.12, P = 0.94). In general, 10 genes were found as hub genes related to stemness indices through PPI analysis. And 4 of 10 hub genes were significantly negatively correlated with mDNAsi. It suggested that the tumor stemness indices could be used as a new objective indicator to predict the prognosis of patients with HCC. SNAP25, KPT19, GABBR1 and EPCAM might key genes for HCC occurrence, metastasis and recurrence, which provided new ideas for targeted treatment of HCC.