Characterization of cancer-related fibroblasts (CAF) in hepatocellular carcinoma and construction of CAF-based risk signature based on single-cell RNA-seq and bulk RNA-seq data

ScRNA-seq data of GSE149614↗ was downloaded from the Gene Expression Omnibus (GEO) database, including 10 samples of primary tumors, 2 samples of portal vein tumor thrombi, 1 sample of metastatic lymph node, and 8 samples of non-tumor liver. For scRNA-seq data, single cells were firstly screened with each gene expressing in at least 3 cells and each cell expressing at least 250 genes. Then PercentageFeatureSet function in Seurat R package was conducted to evaluate the proportion of mitochondria and rRNA. The single cells were further screened by setting each cell expressing at least 6000 genes with UMI > 100. Finally, a total of 69145 cells were remained. The transcriptome data, single-nucleotide variant (SNV) and copy number variants (CNV) data of Masked Copy Number Segment, and corresponding clinical information of HCC were obtained from The Cancer Genome Atlas (TCGA) database. For transcriptome data, the samples without survival data and outcome status were removed, and eventually, 360 tumor samples and 50 para-cancerous samples were included. GSE76427↗ cohort with 115 HCC samples was downloaded from the GEO database as a validation cohort after the removal of normal tissue samples and tumor samples without follow-up and outcome status information. Ten cancer-related pathways (Cell Cycle, HIPPO, MYC, NOTCH, NRF1, PI3K, TGF-Beta, RAS, TP53, and WNT) were retrieved from the literature (20).

We re-analyzed the scRNA-seq data of HCC using the Seurat package (21) to comprehensively characterize the CAF signature. Firstly, we removed the cells with over 6000 or below 250 expressed genes, followed by log normalization of expressed genes. The batch effects for 21 samples were eliminated using the FindIntegrationAnchors function. The non-linear dimensional reduction was conducted using the uniform manifold approximation and projection method, with 15 principal components and a resolution at 0.2. Single cells were clustered into different subgroups by using the functions of FindNeighbors and FindClusters (dim = 40 and resolution = 0.2).Then t-distributed stochastic neighbor embedding (TSNE) dimensional reduction was conducted using the RunTSNE function. Fibroblasts were annotated with 4 marker genes, including ACTA2, FAP, PDGFRB, and NOTCH3. The fibroblasts were re-clustered with the same algorithm of FindNeighbors and FindClusters functions. TSNE dimensionality reduction was further performed on fibroblasts clusters. Marker genes of each CAF cluster were identified using FindAllMarkers function by comparing one cluster with other clusters with logFC = 0.5, minpct = 0.35, and adjusted p-value<0.05. Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis on the marker genes of CAFs clusters using the clusterProfiler package (22), and the CNV characteristics among the CAFs clusters were analyzed using the CopyKAT R package to differentiate between tumor cells and normal cells in each sample (23).

Firstly, the differentially expressed genes (DEGs) between the tumor and normal tissue were screened out via limma package with a false discovery rate (FDR)<0.05 and |log2(Fold Change)|>1 (24). Then, we assessed the correlations between the DEGs and CAF clusters, and identify the key CAF-related genes with p<0.001 and cor>0.4. The prognosis-related genes were further identified using univariate Cox regression analysis in survival package with p<0.05 (https://rdocumentation.org/packages/survival/versions/2.42-3↗). To compress the gene number, we performed the least absolute shrinkage and selection operator (lasso) cox regression analysis, followed by multivariate Cox regression analysis with a stepwise regression method. According to the results of the multivariate Cox model, we constructed a risk signature with the following formula: risk score=Σβi*Expi. Where i is the gene in risk signature, expi represents the expression of the gene i, and βi represents the coefficients of gene i in multivariate Cox model. The patients were divided into the high- and low-risk groups after zero-mean normalization. The receiver operating characteristic curve (ROC) analysis was performed using the timeROC package (https://cran.r-project.org/web/packages/timeROC/index.html↗) to analyze the predictive performance of the risk signature. Similar analyses were conducted in the validation cohort.

The proportions of 22 immune cell subtypes in the TCGA cohort were evaluated by the CIBERSORT algorithm (25), a tool for assessing immune cell infiltration, and the immune and stromal scores were calculated using the ESTIMATE algorithm (https://sourceforge.net/projects/estimateproject/↗) to further explore the TME.

To construct a nomogram model for clinical use, we first perform the univariate and multivariate Cox regression analysis on clinicopathological and risk signature. The variables with p<0.05 in the multivariate Cox model were used to construct a nomogram for the prediction of HCC prognosis using the rms package (26). The calibration curve was generated to evaluate the predictive accuracy of the model. The reliability of the model was evaluated using decision curve analysis (DCA).

We downloaded the transcriptomic, and matched clinical data of patients with HCC treated with an anti-PD-L1 agent (atezolizumab) (27) from IMvigor210 cohort (http://research-pub.gene.com/IMvigor210CoreBiologies↗). Meanwhile, GSE78220↗ cohort comprised of transcriptomic data from pre-treatment melanomas receiving anti-PD-1 checkpoint inhibition therapy (28), and also download for the determination of the potential value of the risk signature score in the prediction of responsiveness to immune checkpoint blocks (ICB).

All statistical analyses were performed using the R software (v3.6.3). The correlation matrices were conducted using Pearson or Spearman correlation. Wilcoxon test was conducted for the comparisons between the two groups. Survival differences were compared using K–M curves with a Log-rank test. P-value < 0.05 was considered statistical significance.

The flow chart of this study was shown in Figure 1. A total of 69145 cells were obtained from the scRNA-seq data after initial screening (Table 1). The detailed results of data preprocessing were shown in Figure S1. After log-normalization and dimensionality reduction, 15 subpopulations were obtained, and 9 CAF populations were identified based on four marker genes, including ACTA2, FAP, PDGFRB, and NOTCH3 (Figures S2A, B). The cells of 9 CAF populations were extracted for further clustering and dimensionality reduction. The CAF populations were further clustered by using the same clustering algorithm and four CAF clusters were identified (Figures S2C, D). The epithelial cell specific gene was not expressed in all four CAF clusters, supporting the accuracy of CAF identification (Figure S3). Figure 2A showed the TSNE plot of 21 sample distributions. As a result, four CAF clusters were finally generated and used for subsequent analysis (Figure 2B). A total of 211 DEGs among the 4 CAF clusters were identified and the expression of the top 5 DEGs (determined as the marker genes of CAF clusters) in the 4 clusters was shown in Figure 2C. The proportion of the 4 clusters in each cohort were illustrated in Figure 2D. As shown in Figure 2E, the results of KEGG analysis demonstrated that these DEGs were enriched in multiple pathways, including vascular smooth muscle contraction, focal adhesion, oxytcosin signaling pathway, PPARG signaling pathway, etc. In addition, the 4 CAF clusters consist of 1533 tumor cells and normal cells according to the CNV characteristics (Figure 2F).

To elucidate the associations between the CAF clusters and tumor progression, we investigated the characteristics of ten tumor-related pathways in the four CAF clusters. The GSVA scores of the ten tumor-related pathways in different CAF clusters were shown in Figure 3A. The ratio of malignant cells in CAF_0 cluster was significantly higher than that in the other three clusters (Figure 3B). However, there were no significant differences among the CAF_1, CAF_2, and CAF_3. Furthermore, we analyzed the GSVA scores of the ten tumor-related pathways between malignant and non-malignant cells in each CAF cluster, with slight differences observed (Figures 3C-F).

To determine the associations between the CAF clusters and prognosis, we first calculated the ssGSEA score of the marker genes (the top 5 DEGs of CAF clusters defined in Figure 2C) of each CAF cluster based on the TCGA cohort. The results demonstrated that the CAF_2 cluster had a significantly higher score in tumor samples than in normal samples, whereas the other CAF clusters had an opposite trend, with a higher score in normal samples than in tumor samples (Figure 4A). The HCC samples of TCGA dataset were separated into the high- and low-CAF score groups according to the optimal cut-off value analyzed by survminer R package. The samples in the high-CAF score group had a better prognosis in the CAF_0, CAF_1, and CAF_2 clusters than those in the low-CAF score group, whereas the CAF_3 was not associated with the prognosis of HCC (Figures 4B-E). The above results suggested that CAF_3 cluster may contribute little in the HCC progression although CAF_3 enrichment was differential in HCC and normal samples.

To construct a risk signature, we firstly screened out DEGs between the tumor and normal tissues. As shown in Figure 5A, a total of 2349 DEGs were obtained, with 462 up-regulated DEGs and 1887 down-regulated DEGs. Among them, there are 423 genes that showed significant correlations with those prognosis-related CAF clusters. Furthermore, the prognostic value of each gene was assessed via univariate Cox regression analysis, with 234 genes exhibiting prognostic values (Figures 5A, B). Lasso Cox regression analysis was performed to narrow down the number of genes, with 11 genes left as lambda=0.047 (Figures 5C, D). Finally, we included 6 genes, including HMG-box containing 3 (HMGXB3), GCN1 activator of EIF2AK4 (GCN1), LUC7 like 3 pre-mRNA splicing factor (LUC7L3), ADAMTS like 2 (ADAMTSL2), solute carrier organic anion transporter family member 2A1 (SLCO2A1), and CD4 molecule (CD4), in the risk signature after multivariate Cox regression analysis with stepwise regression method (Figure 5E). The final 6-gene signature formula is as follows: RiskScore = -0.088*ADAMTSL2 - 0.121*SLCO2A1 - 0.217*CD4 + 0.249*GCN1 + 0.345*HMGXB3 + 0.271*LUC7L3. We calculated the risk score for each sample and divided them into the high- and low-risk groups after z-mean normalization. The AUC values of the model for 1- to 5-year survival range from 0.68 to 0.76 in the TCGA cohort and range from 0.65 to 0.7 in the GEO cohort (Figures 5F, G). Kaplan-Meier survival analyses revealed that high-risk patients had significantly poorer survival outcomes compared with low-risk patients in the TCGA cohort, as well as in the GEO cohort (Figures 5H, I).

Next, we checked out the SNV mutations of the six genes of the risk signature. It showed that ADAMTSL2, SLCO2A1, HMGXB3, LUC7L3, and CD4 have SNV mutations in more samples, while no SNV mutation was observed in GCN1 (Figure S4A). We analyzed the co-occurrence probability of these key genes and the 10 most mutated genes. As revealed in Figure S4B, there was no significant probability of co-occurrence of the mutations in these 5 genes, but LUC7L3 presented a significant probability of co-occurrence with ABCA13 mutation. In the 6 genes, it was found that only a very small number of samples had gain/loss of CNV (Figure S4C). To further elucidate the associations between the risk genes and HCC, we analyzed the correlations between these genes and several molecular signatures of HCC. The results demonstrated that SLCO2A1 had significantly negative correlations with Aneuploidy Score, Homologous Recombination Defects, Fraction Altered, Number of Segments, and Nonsilent Mutation Rate, whereas HMGXB3, LUC7L3, and GCN1 showed significantly positive correlations with Homologous Recombination Defects and Fraction Altered (Figure S4D). In addition, we analyzed the potential pathways associated with each risk gene. As shown in Figures 6A, B, a total of 39 pathways were significantly correlated with these six genes, including angiogenesis, apical junction, apoptosis, etc.

Our data demonstrated that ADAMTSL2, SLCO2A1, and CD4 presented significantly positive correlations with the stromal score, immune score, and estimate score, while LUC7L3 showed significantly negative correlations with the stromal, immune, and estimate scores. However, there was no significant correlations observed between the three scores and the other genes (GCN1 and HMGXB3) (). After grouping according to the median value of expression of each gene, we compared the three scores in different expressed groups. The results showed that, with regard to the ADAMTSL2, SLCO2A1, and CD4 genes, the three scores of the high expression group were significantly higher than those of the low expression group (). Correlation analysis revealed that ADAMTSL2, SLCO2A1, and CD4 presented a significantly negative correlation with the majority of T cells. Additionally, LUC7L3, GCN1, and HMGXB3 significantly positively correlated with M0 macrophages and neutrophils (). Moreover, we also observed significant differences between the high and low expression groups of risk genes in several immune cells (). 5 5 5 5

T-cell immunotherapy has emerged as an anticancer treatment with synergistic survival benefits (29). Therefore, we assessed the prognostic value of risk signature for immune-checkpoint therapy in the IMvigor210 and GSE78220↗ cohorts. The 348 patients in the IMvigor210 cohort showed varying degrees of response to anti-PD-L1 receptor blockers, including complete response (CR), partial response (PR), stable disease (SD), and progressive disease (PD). SD/PD patients presented higher risk scores than CR/PR patients (Figure 7A). In the high-risk group, the percentage of SD/PD was higher than that in the low-risk group (Figure 7B). We observed that in the IMvigor210 cohort, patients in the low-risk group showed significant clinical benefits and a significantly longer overall survival as compared with those in the high-risk group (Figure 7C, p=0.0053). Specifically, there were significant survival differences in Stage I+II patients between the different risk groups (Figure 7D, p=0.0017), but not in Stage III+IV patients (Figure 7E, p=0.5). It suggested that the risk score is more sensitive in early-stage patients. In the GSE78220↗ cohort, we also found a significantly longer overall survival of patients in low-risk group than in high-risk group (Figure 7F, p=0.036). Meanwhile, the percentage of SD/PD in the high-risk group was higher than that in the low risk group (Figure 7G).

To optimize the predictive performance of the risk signature, we integrated the clinicopathological characteristics and risk score via univariate and multivariate Cox regression analysis. Multivariate analysis demonstrated that risk signature was the most significant independent prognostic factor of osteosarcoma [hazard ratio (HR) = 1.77, 95% confidence interval (CI): 1.42 - 2.13, P < 0.001], followed by metastatic status (HR = 1.74, 95%CI: 1.22 - 2.46, P = 0.002) (Figures 8A, B). Therefore, a nomogram combining stage and risk score was constructed, as shown in Figure 8C. The calibration plot demonstrated that the nomogram can effectively forecast the actual survival outcomes (Figure 8D). Moreover, DCA revealed a better discriminative ability of the nomogram in recognizing patients at high risk than the risk score and stage, as shown in Figure 8E. TimeROC analysis showed that the AUC of the risk score and nomogram was higher than that of other indicators in the TCGA cohort (Figure 8F).

The original contributions presented in the study are included in the article/. Further inquiries can be directed to the corresponding authors. 1