Abstract
Introduction:
The complexity of tumor cell subclonal structure has been extensively investigated in hepatocellular carcinoma. However, the role of subclonal complexity in reshaping the tumor microenvironment (TME) remains poorly understood.
Methods:
We integrated single-cell transcriptome sequencing data from four independent HCC cohorts, involving 30 samples, to decode the associations between tumor subclonal complexity and the TME. We proposed a robust metric to accurately quantify the degree of subclonal complexity for each sample based on discrete copy number variations (CNVs) profiles.
Results:
We found that tumor cells in the high-complexity group originated from the cell lineage with FGB overexpression and exhibited high levels of transcription factors associated with poor survival. In contrast, tumor cells in low-complexity patients showed activation of more hallmark signaling pathways, more active cell-cell communications within the TME and a higher immune activation status. Additionally, cytokines signaling activity analysis suggested a link between HMGB1 expressed by a specific endothelial subtype and T cell proliferation.
Discussion:
Our study sheds light on the intricate relationship between the complexity of subclonal structure and the TME, offering novel insights into potential therapeutic targets for HCC.
Introduction
Hepatocellular carcinoma (HCC) is the third most lethal cancer in the world, with patients often experiencing post-surgical recurrence and metastasis, leading to poor survival outcomes (). Tumor heterogeneity is the main cause of drug resistance and treatment failure (; ; ), which encompasses the cellular diversity and dynamics both among tumors in different patients and within individual tumor. The tumor heterogeneity is manifested through genomic mutations, aberrations in epigenetic modifications, transcriptional alterations, and changes at the protein level (). Additionally, extrinsic factors such as hypoxia, pH levels, and the interactions between tumor cells and other stromal components within the tumor microenvironment (TME) contribute to the diversity and dynamics in tumor genotypes and phenotypes (). Previous studies, such as those conducted by TCGA and ICGC, have integrated multi-omics profiles to delineate the heterogeneous molecular features of cancers and explored the subtyping of cancers, guiding personalized therapeutic strategies (). Tumor heterogeneity also underpins tumor evolution, as tumor cells continuously accumulate genetic and epigenetic variations, forming heterogeneous subclones which possess different survival fitness and undergo dynamic processes of subclonal selection and elimination under selective pressures such as nutrient availability and metabolism (; ). Therefore, understanding the complexity of subclonal structure within tumor is crucial for elucidating tumor biology, disease progression, and response to therapy.
The recent advancements in single-cell sequencing technology (; ) have provided unprecedented insights into the composition, function, and spatial organization of immune and stromal cells within the tumor microenvironment (TME), leading to a deeper understanding of its complexity and diversity. In the process of tumor progression, tumor heterogeneity continuously remodels the tumor microenvironment by altering the transcription levels of target genes in non-malignant cells within the TME, which leads to divergences in developmental trajectories, immune landscapes, and intercellular networks (; ; ; ). Specifically, heterogeneous secretion of cytokines by tumor cells can also modulate immune-tumor interactions, as demonstrated by studies in mouse models (; ). Tumor cells can induce profound phenotypic changes in non-immune stromal components within the TME, impacting the immune components and their activities. For instance, oncogenic BRAFV600E signaling in human melanoma cells has been shown to interfere with T cell-mediated anti-tumor responses by modulating the phenotype of cancer-associated fibroblasts (; ). Furthermore, immune cells within the microenvironment can eliminate clones with high mutational burdens (high immunogenicity) through the process of immune editing, maintaining a balance that ultimately favors the selection of tumor variants capable of evading immune surveillance (; ; ; ; ). Thus, the genotype and phenotype of tumor cells and the TME are inextricably linked. However, the intricacies of their interactions and how the subclonal complexity reshapes the TME remain largely unexplored.
In this study, we performed an integrated single-cell analysis of HCC patients with untreated primary tumors. We developed a metric to quantify the tumor subclonal complexity based on Shannon entropy theory. Our investigation focused on systematically elucidating the relationship between the subclonal complexity and cancer cell state, the functional activity of the immune cells. We found that cell-cell interaction landscape and TME polarization varied with subclonal complexity. Our findings underscore the significance of subclonal complexity in modulating the tumor microenvironment and its implications for immune responses.
Materials and methods
scRNA-seq datasets
The datasets of HCC were acquired from GEO, including 4 scRNA-seq cohorts (GSE156625 (24), GSE149614 (25), GSE151530 (26), GSE189903 (27)). Specifically, GSE156625 (https://www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA658541&o=acc_s%3Aa) was the raw sequencing data and the other 3 datasets are processed read count matrices. To accurately analyze the impact of tumor subclonal complexity on the remodeling of the tumor microenvironment, we established the following criteria for sample selection: 1) all samples are from the 10X Genomics platform; 2) availability of raw sequencing data or read count expression matrices; 3) samples are untreated; 4) samples are primary tumors without invasion, recurrence and metastasis; 5) for multi-region sequencing data, we selected samples from the tumor core region whenever possible. Finally, only 30 samples were retained for further analysis (Supplementary Table S1).
Pre-processing of scRNA-seq data
Raw reads from GSE156625 were processed to generate gene expression matrices by using the standard internal pipeline based on the Cell Ranger toolkit (3.1.1). Then the expression matrices of each individual sample were converted to Seurat object by the “CreateSeuratObject” function from Seurat package (4.3.0) (). We firstly filtered out cells that had either lower than 200 or higher than 8,000 expressed genes. Furthermore, the cells with the percent of mitochondrial genes over 40% of total expressed genes were discarded. Additionally, we initially removed doublet cells by DoubletFinder (2.0.3) () to avoid the potential influence on analytical results. A total of 67,125 cells were retained for further analysis.
Dimensionality reduction, clustering, DEGs identification, and major cell type annotation
To better reflect the biological features from the data, we used the “NormalizeData” function to normalize the expression matrices for the Seurat object. And the top 2,000 highly variable genes (HVGs) identified by “FindVariableFeatures” function were used to scale data through “ScaleData” function. The principal component analysis (PCA) was carried out by “RunPCA” function to generate 50 PCs, followed by “RunHarmony” function implanted in Harmony package (0.1.1) () with the top 30 PCs to eliminate batch effect. The clustering analysis was implemented depended on the integrated joint embedding produced by Harmony with the Louvain algorithm after computing a shared nearest-neighbor (SNN) graph. The UMAP technique was used to create a 2D map on which the discovered clusters were displayed. Then we used “FindAllMarkers” function to identify over-expressed genes in specific cluster when compared with the other clusters (adjusted P-value <0.05, only.pos = T and logfc.threshold = 0.25). The well-known cell markers (; ) were used to annotate the clusters: Endothelial (ENG, VWF, PECAM1), Fibroblasts (ACTA2, COL1A2, PDGFRB), T cells (CD2, CD3D, CD3E, CD3G), B cells (CD79A, CD79B, MS4A1, BANK1), Myeloid (CD14, CD68), Plasma (MZB1, DERL3, SDC1), Hepatocytes (APOA2, ALB, APOA1, AMBP), proliferative T (TOP2A, MKI67, TUBB).
Cell subtype annotation
We further clustered T cells, endothelial cells and myeloid cells following a similar pipeline as described above, including normalization, highly variable genes identification, dimensionality reduction, batch effect correction with Harmony and clustering. Within each cell subtype, an iterative procedure that a cell cluster displaying at least two canonical marker expression of 7 major cell type was determined as doublet cells to be removed, if any, and reclustered the remaining cells, was adopted to ensure the reliability of the study. Specifically, we used the following well-known markers (; ) for subtype identification. NK (GNLY, KLRF1, NKG7), CD4+_CCR7 (CD4, CCR7), CD4+_IFNG+ (CD4, IFNG), CD4+_FOXP3_Treg (CD4, FOXP3, IL2RA), CD8+_GZMK (CD8A, GZMK), CD8+_TNFSF9 (CD8A, TNFSF9), macrophage (C1QC, APOE, CD68), monocyte (FCN1, LYZ), cDC1_CLEC9A (CLEC9A, FLT3), cDC2_CD1C (CD1C, FCER1A), cDC3_LAMP3 (LAMP3, CCR7), mast (KIT). Similarly, endothelial cell were also confirmed by canonical markers (; ), such as activated_PCV (ACKR1, SELP, VCAM1, POSTN), Capillaries_EDNRB (EDNRB, CA4, HPGD, IL1RL1), Capillaries_FCN3 (FCN3, BTNL9, NOSTRIN, EDN), arteries_GJA5 (GJA5, EFNB2, SOX17, DKK2), tip_CXCR4 (CXCR4, ADM, ANGPT2, APLN), RBP4+_EC (RBP4, IGFBP7), CXCL10+_EC (CXCL10, GBP1, SOD2). To accurately annotate proliferative T cells, we calculated the averaged Pearson correlation coefficient and averaged p-value for each proliferative T cell against other T cell subtypes and NK subtypes based on single-cell expression profiles. Proliferative T cells were excluded if the mean correlation coefficients were all less than 0.4 or the mean p-values all exceeded 0.05. Otherwise, they were assigned to the subtype with the highest average correlation.
CNVs estimation and malignant cells identification
We inferred copy number variations (CNVs) for each cell through the InferCNV (1.10.1) package of R. The non-malignant cells were applied as the reference to estimate the CNVs values of the hepatocytes (observation). The genes were sorted by their genomic locations on each chromosome. Specifically, The InferCNV analysis was implemented with parameters “denoise = T, HMM = T, analysis_mode = “subclusters”, cutoff = 0.1, cluster_by_groups = T”.
To identify the malignant cells from the hepatocytes, we designed the CNVs score quantifying the degree of CNVs fluctuation for each cell as follows:where represented the CNVs value of gene i in the cell, was the mean CNVs value of gene i in the cell, and n was the number of the genes in the cell. If the CNVs score of a hepatocyte exceeded the threshold determined by the intersection of the distributions of CNVs score in hepatocytes and non-malignant cells, then the hepatocyte was considered malignant (Supplementary Figure S1C).
Identification of intratumor NMF programs
The non-negative Matrix Factorization (NMF) algorithm was applied to identify the underlying expression programs from the malignant cells based on the NMF R package. We employed NMF (rank = 2:6, nrun = 10) to the relative expression matrices of the top 2000 HVGs in each sample with all negative values converted to 0. We chose the optimal rank value at which the cophenetic coefficient started producing the maximum descent (Supplementary Figure S1E). In total, we identified 98 programs from 30 samples (Supplementary Table S2). The 30 genes of each expression program with the top NMF scores were input into “AddModuleScore” function to generate the program score for malignant cells in each sample. We calculated the correlations between the 98 program scores in each sample individually. After excluding clustering groups containing less than half of samples, we finally extracted 3 meta-programs from the 98 programs according to hierarchical clustering of averaged correlations of pairs of programs across all samples (Figure 1F; Supplementary Table S3). Each meta-program, retaining 30 genes with the highest average NMF score, was used to perform pathway enrichment analysis.
FIGURE 1
CNV complexity score based on single cell CNVs profiles
Shannon entropy was used to quantify the CNV complexity for each sample based on discrete CNVs profiles inferred from InferCNV. Firstly, we calculated the frequency of copy number gains and losses for each gene within the tumor as follows:where n represented the total number of malignant cells, was the number of cells in which the copy number of the gene was lost, and was the number of cells in which the copy number of the gene was gained. Specifically, Discrete CNVs value of the gene greater than (less than) 3 was considered to be gained (lost). The frequency of copy number gain and loss were divided into 10 equal-length frequency intervals separately (i.e., Loss:0–0.1, 0.1–0.2, …, 0.9–1; Gain:0–0.1, 0.1–0.2, …, 0.9–1). After mapping all the genes to equal-length frequency intervals according to the and of the genes, we then counted the number of genes belonging to each equal-length frequency interval and calculated the frequency of the number of genes in each interval. Finally, the CNV complexity score could be calculated for each sample as follow:where i was the ith interval, was the number of genes falling into the ith interval, and was the frequency of the number of genes in the ith interval.
Simulation data and method evaluation
To evaluate whether the CNV complexity score could accurately distinguish between polyclonal and monoclonal samples, we used computer simulations to generate single-cell CNV profiles for 100 polyclonal samples and 100 monoclonal samples as the ground truth data, respectively. Briefly, the size of cell population N ∼ U [200, 2000] and the total number of clones C ∼ U [5, 20] were generated for each sample. For each clone, the size of clone was generated through binomial distribution B (N, p). For monoclonal samples, we assumed that p = 0.9 for one dominant clone and p = 0.1/(c-1) for other rare clones. We considered int (N/2) subclones in the tumor cells, where int takes the nearest integer in the polyclonal samples. We used p = 2/N to generate the size of subclone. For detailed steps of the entire simulation data generation process, please refer to . We then applied the three methods (CNV complexity score, Guo et al.’s and Ma et al.’s) to the simulated single-cell CNV profiles which had been labeled as either monoclonal or polyclonal to calculate scores. A higher score for this sample means that the algorithm is more likely to predict the sample as a polyclonal sample. Finally, we used R package pROC () to calculate precision, recall, and the area under the receiver operating characteristic curve (AUC).
Gene set enrichment analysis, pathway enrichment analysis and ssGSEA
Gene set enrichment analysis (GSEA) and pathway enrichment analysis were performed using R package clusterProfiler (). The ssGSEA algorithm was used to calculate the activity score of the gene sets for each cell based on GSVA package (). We obtained 13 cancer functional states from CancerSEA (), as well as a collection of 50 cancer Hallmark gene sets downloaded from MsigDB (). The signature gene lists of macrophages (e.g. Angiogenesis, Phagocytosis, MHC I molecules, MHC II molecules; Supplementary Table S5) and T cells (e.g. Cytotoxic, Exhaustion, Terminally exhaustion, Proliferation, TCR signaling pathway, Type II interferon response; Supplementary Table S6) were collected from previous studies (; ; ; ; ).
CytoTRACE and pseudotime analysis
The R package CytoTRACE v.0.3.3 (47) was applied to predict the differentiation state of cells from the single-cell RNA-seq (scRNA-seq) profiles. The CytoTRACE score range from 0 to 1, with higher scores being associated with greater stemness (less differentiation) and vice versa, which help us determine the starting point of differentiation. We employed Monocle2 () to discover the differentiation trajectory for T cells and endothelial cells. The top 2,000 DEGs between the cell subgroups were used to create the DDRTree. Utilizing the root state determined by the CytoTRACE, each cell was given a pseudotime value using the “order_cells” function. The cells were ordered and visualized with the “plot_cell_trajectory” function.
TF regulatory network construction and critical TFs identification
The SCENIC () was used to construct TF regulatory network for malignant cells in different subclonal complexity. Briefly, GRNBoost2 estimated the co-expression network and RcisTarget was used to identify the regulons. The regulon activity was then measured by AUCell for every cell. The differently activated TFs regulons between different subclonal complexity were identified by wilcox.test(). The TF-gene regulatory network was visualized by Cytoscape ().
Evaluation of associations between critical TFs and survival outcomes
To explore the prognostic power of TFs, we downloaded two gene expression datasets from cBioPortal () and GEO (GSE76427) for HCC patients, which contained clinical information. HCC patients from TCGA were treated as a training dataset, and HCC patients from GSE76427 were used as an independent validation dataset. A univariate Cox proportional hazards regression analysis was performed to evaluate the association between the critical TFs and patients’ OS. Only TFs with p-value <0.05 were selected to further conduct variable selection according to stepwise Cox proportional hazards regression analysis. Finally, we created a risk-score formula based on TFs expression weighted by Cox regression coefficients.where N was the number of selected TFs, was the expression value of the ith gene, and was the Cox regression coefficient of the ith gene in the univariate Cox regression analysis. Patients were divided into high-risk and low-risk group based on the median risk score. KM curves were performed to compare the OS and DFS between two groups. To evaluate whether risk score could be an independent prognostic factor, Multivariate Cox proportional hazards regression model was conducted with risk score, age, sex, height, weight, race as covariates.
Cell-cell communications analysis
We applied CellPhoneDB () to infer cell-cell communication between different cell types. The specific ligand-receptor interactions between cell types were identified based on permutation test. The interaction score represented the total mean of the individual ligand-receptor partner average expression value in the corresponding interacting pairs of cell types. To predict ligands driving the transcriptomic changes of target cells, we further performed NicheNet analysis () between RBP4+_EC and T cells with the “nichenet_seuratobj_aggregate” function. Specifically, the differentially expressed genes (DEGs) among T cells between different subclonal complexity with high complexity group as control were used as gene sets of interest and all expressed genes in T cells were used as background of genes.
Cytokines signaling activity calculation
We used CytoSig algorithm () to quantify the signaling activity of cytokines for T cells. The CytoSig was a linear model predicting cytokines signaling activity score for each cell and covered 51 cytokines provided by the curated training data. We implemented CytoSig for single-cell transcriptional profiles of T cells in Python 3.
Statistical analysis
Continuous variables were compared by the Wilcoxon rank-sum test and categorical variables were compared by the chi-squared test. Pearson correlation was calculated by cor () and cor. test (). The log-rank test was used to compare Kaplan–Meier curves. When p < 0.05, the differences were considered to be significant. All statistical analyses were conducted using R software version 4.1.2 (http://www.rproject.org).
Results
The transcriptional diversity analysis reveals transcriptional programs in primary HCC
We downloaded 30 untreated primary HCC samples from four cohorts of single-cell RNA-seq data (; ; ; ) (Supplementary Table S1). After quality control, we obtained a total of 66,609 high-quality cells with a median of 2,220 cells per sample. Cell types were annotated by canonical cell markers, identifying 31,832 hepatocytes and 34,777 non-malignant cells (Figure 1A; Supplementary Figure S1A). We inferred copy number variations (CNVs) based on scRNA-seq data using InferCNV(()) (Supplementary Figure S1B). Malignant and non-malignant hepatocytes were distinguished by their accumulated CNV scores (see Section Methods, Supplementary Figure S1C), ultimately identifying 30,955 malignant cells.
We performed clustering and UMAP visualization separately for malignant and non-malignant cells. Consistent with previous findings, non-malignant cells exhibited distinct separation by cell type but showed a mixture of different samples (Supplementary Figure S1D; Figure 1B), whereas malignant cells exhibited high inter-tumor heterogeneity (Figure 1C) (). The transcriptomic similarity among inter-tumor cells was lower than that among intra-tumor cells (Figure 1D). However, the similarity of intra-tumor malignant cells varied across patients; for instance, HCC04T showed relatively high similarity, while H72 exhibited lower similarity among malignant cells (Figure 1E).
To identify underlying transcriptional programs, we applied non-negative matrix factorization (NMF) () to each sample. The identified gene programs tend to be mutually exclusive among malignant cells within sample, and were negatively correlated (Supplementary Figure S1E–G). A total of 98 gene programs were identified across the 30 tumor samples (Supplementary Table S2). We then performed hierarchical clustering on these 98 expression programs, extracting three meta-programs (Figure 1F; Supplementary Table S3) that encompassed more than half of the samples. These three meta-programs were associated with distinct signaling pathways: Meta-program 1 was enriched in cell cycle-related pathways, Meta-program 2 was enriched in metabolic pathways, and Meta-program 3 was enriched in stress response-related pathways (Figure 1G).
Quantifying intra-tumor subclonal complexity using CNV profiles
To explore intra-tumor subclonal complexity, we proposed a metric to quantify CNV complexity based on the discrete CNVs profiles inferred from InferCNV (Figure 2A, method). We quantified CNV complexity for each HCC sample and classified samples into high (H) and low (L) groups based on the median CNV complexity score. We found that CNV complexity scores were significantly correlated with two other heterogeneity scores proposed by () and () respectively (Figures 2B, C). Moreover, the CNV complexity score was more robust than the other two heterogeneity scores, both of which were negatively correlated with the number of the cells (Figures 2D, E). Simulation data further demonstrated the robustness of the CNV complexity score (Supplementary Figure S2A).
FIGURE 2
Finally, we evaluated the three methods on simulated datasets generated in silico (Supplementary Figure S2B, C). The CNV complexity score showed excellent performance in reflecting the degree of subclonal complexity, as demonstrated by a receiver operating characteristic (ROC) curve with the highest AUC value, which was superior to the other methods (Figure 2F). Overall, our constructed CNV complexity score demonstrates excellent stability and sensitivity, and more importantly, it is not influenced by the number of cells, providing a more accurate quantification of subclonal complexity in samples.
Phenotypic regulatory plasticity in HCC samples with varying subclonal complexity
To explore the phenotypic plasticity underlying different subclonal complexity, we used cancer functional gene sets (Supplementary Table S4) along with ssGSEA and identified 37 gene sets with significant differences between the H and L groups, with 34 functional gene sets being significantly upregulated in the L group (p < 0.05, Wilcoxon rank-sum test, Figure 3A; Supplementary Figure S3). For instance, the L group was significantly enriched for hypoxia, invasiveness, metastasis, emphasizing varying signaling activity between the H group and L group. Malignant cells are susceptible to epithelial-mesenchymal transition under hypoxic conditions, gaining invasive and metastatic abilities, accompanied by a transition from low to high stemness. Consistent with hypoxia, invasion and metastasis activities mentioned above, it was discovered that the stemness of malignant cells in L group was also significantly higher than that in the H group (p = 0.0049, Wilcoxon rank-sum test, Figure 3B) (), which highlights that stemness of malignant cells is closely related to hypoxia, invasiveness and metastasis.
FIGURE 3
We further investigated the cell lineage origin of malignant cells using scHCL(), and found significantly different distributions of cell lineages between the H and L group (p < 0.05, chi-squared test, Figure 3C), with a significantly higher proportion of cell lineage overexpressing FGB in the H group (p < 0.05, Wilcoxon rank-sum test, Figure 3D). FGB was significantly upregulated in malignant cells from the H group (p < 2.22e-16, Wilcoxon rank-sum test, Supplementary Figure S4A). Altered expression of this gene leads to several disorders, including afibrinogenemia, dysfibrinogenemia, and thrombotic tendency. In summary, malignant cells from samples with varying CNV complexity might undergo diverse transcriptional reprogramming, resulting in different cancer functional states and cell lineage origins.
To investigate the regulatory mechanisms contributing to the different cell lineages and functional phenotypes of tumor cells between the H and L group, we employed SCENIC to identify transcription factors (TFs) regulatory modules (Figure 3E). We then constructed TFs regulatory networks corresponding to the H and L groups using differentially expressed and differential CNV genes targeted by identified TFs (Supplementary Figure S4B; Supplementary Figures S5A, B). Functional enrichment analysis of the TF regulatory network revealed that the upregulated genes in the L group’s TFs regulatory network were enriched in pathways such as hypoxia, apoptosis, TNFA signaling via NFkB and regulation of Wnt signaling. In contrast, the upregulated genes in the H group were involved in pathways regulating epithelial cell proliferation and metabolism-associated pathways (Supplementary Figure S4C).
Among these TFs, we identified four transcription factors MAFG, YY1, EZH2, AR, associated with overall survival (OS) which were upregulated in the H group. Samples were stratified into high-risk and low-risk groups based on the median risk score quantified by these four TFs. Kaplan-Meier curve analysis demonstrated poorer survival in the high-risk group compared to the low-risk group in terms of OS and disease-free survival (DFS) (OS: p < 0.0001, DFS: p < 0.0001, log-rank test, Figure 3F). Further multivariate Cox proportional hazards regression analysis confirmed the independent prognostic value of the risk score (Figure 3G; Supplementary Figure S4D). Finally, we validated the prognostic power of the 4-TF signatures in an external validation dataset (GSE76427). The KM curves of the two groups were significantly different (p = 0.0016, log-rank test, Supplementary Figure S4E) and the risk score showed the independent prognostic value (Supplementary Figure S4F).
Characterizing cell-cell communications in the tumor microenvironment across different subclonal complexity
To investigate the characteristics of cell-cell communications within the tumor microenvironment under different complexity of subclonal structure, we performed cell-cell interaction analysis using CellPhoneDB () between malignant cells and non-malignant cells in the H and L groups, respectively. We found that fibroblast and endothelial cells strongly interacted with tumor cells in the L group rather than the H group (Figure 4A). Additionally, myeloid cells, including macrophage, monocytes and dendritic cells (Supplementary Figure S6A, B), showed strong interactions with tumor cells in both the H and L groups, with increased communications observed in the L group (Figure 4B).
FIGURE 4
We further explored the ligand-receptor interactions between tumor cells and myeloid cells in the L group and H group. We found that many chemokine-mediated interactions were specifically present in the L group, such as the CCL4L2-VSIR, FCGR2A-CXCL9, CCL3L3-CCR1, CCL3L1-CCR11, IDE-CCL23, and CXCL11-DPP4, CCL2-CCR10 pathways which suggests that tumor cells in the L group might recruit more myeloid cells. In the H group, there were also a few specific interactions, such as the interactions between the dendritic cell-expressed ligands CCR6 and CXCR3 with the tumor cells-expressed receptor CCL20 (Figure 4C). Additionally, we confirmed that coinhibitory interactions mediated by pathways like CSF1-SIRPA, LAG3-FGL1, CCL4L2-VSIR, and LGALS9-HAVCR2 were specifically activated in the L group (Figure 4D). Similar results were also observed in the costimulatory-mediated interactions between malignant cells and myeloid cells (Supplementary Figure S6C). Interestingly, the L group was enriched for both coinhibitor- and costimulator-mediated ligand-receptor interaction pairs. In summary, these results revealed that varying subclonal complexity could alter the cell-cell interactions landscape.
Given the differential interaction strength of myeloid cells and tumor cells between the H and L groups, we investigated the functional characteristics of myeloid cells. Myeloid cells from the L group showed upregulation of MHC I and MHC II molecules and higher angiogenesis and phagocytosis signature scores (p < 0.05, Wilcoxon rank-sum test, Supplementary Table S5; Figures 4E, F). Based on gene set enrichment analysis (GSEA) of the differentially expressed genes identified in myeloid cells between the L and H groups (Supplementary Figure S6D), more hallmark signaling pathways were enriched in myeloid cells from the L group, including hypoxia, EMT, TNFA signaling via NFkB, DNA repair, apoptosis, and angiogenesis (Figure 4G; Supplementary Figure S6E, F). These results were consistent with what we observed in tumor cells from samples with different complexity of subclonal structure.
Differential activation and functional states of T Cells with varying subclonal complexity
Although there was no significant difference in tumor-T cells interaction strength between the H and L groups, T cells exhibited significantly different activation statuses based on ssGSEA analysis of functional gene sets (Supplementary Table S6). Both CD4+ and CD8+ T cells from the L group had higher scores for TCR signaling pathway, exhaustion and proliferation but lower type II interferon response score (p < 0.05, Wilcoxon rank-sum test, Figure 5A). Moreover, cytotoxicity-related genes (e.g., GZMA, GZMB, GZMH, NKG7, PRF1) and immune checkpoint genes (e.g., LAG3, PDCD1, LCK) were significantly upregulated in T cells from the L group (p < 0.05, Wilcoxon rank-sum test, Supplementary Figure S7A, B). The scores of these functional gene sets were significantly correlated with CNV complexity score, implying the potential of subclonal complexity in predicting the activation status of T cells (Figure 5B).
FIGURE 5
We performed subclustering on T/NK cells, resulting in seven cell subpopulations: CD4_CCR7, CD4_IFNG+, CD4_FOXP3+_Treg, CD8_GZMK+, CD8+_TNSF9+, proliferative T cells, and NK cells (Figure 5C; Supplementary Figure S7C). The proliferative T cells were further annotated according to their similarity with other T subtypes and NK subtype (see methods, Figure 5D). Based on trajectory analysis of CD4+ and CD8+ T cells, we found that both CD4+ and CD8+ T cells from the L group tended to locate in the intermediate of the trajectory whereas those from the H group were primarily located at the terminal of trajectory (Figure 5E, F). Additionally, T cells from L group showed increased CytoTRACE score, indicating a heightened activation status of T cells in the L group (CD8+: p = 0.043, CD4+: p = 0.075, Wilcoxon rank-sum test, Figures 5E, F; Supplementary Figure S7D, E).
Immune activation mediated by endothelial cell subtypes with varying subclonal complexity
To investigate the regulatory mechanisms underlying the different activation states of T cells, we estimated the signaling activities of 51 cytokines on T cells and their correlation with T cell proliferation using the CytoSig model (). The activity of the cytokine HMGB1 showed a significant positive correlation with the proliferative score of T cells (p = 0.014, R = 0.45, Figure 6A). Moreover, HMGB1 signaling activity was significantly higher in the L group compared to the H group (p = 0.011, Wilcoxon rank-sum test, Figure 6B). HMGB1 was primarily expressed in endothelial and fibroblast cells (Figure 6C). However, the expression levels of HMGB1 in endothelial cells, rather than fibroblast cells, were upregulated in the L group (p = 0.0023, Wilcoxon rank-sum test) and significantly correlated with T cell proliferation (Figures 6D, E; Supplementary Figure S8A). Additionally, the expression levels of HMGB1 in endothelial cells was positively correlated with the activity scores of T cell exhaustion and TCR signaling pathway, and negatively correlated with the activity score of type II interferon response score (Supplementary Figure S8B).
FIGURE 6
We further dissected the subtypes of endothelial cells, identifying seven distinct endothelial cell subtypes (Figure 6F; Supplementary Figure S8C). The proportion of RBP4 positive endothelial cells (RBP4+EC) was significantly higher in the L group compared to the H group (p < 0.05, Wilcoxon rank-sum test, Figure 6G). Moreover, HMGB1 expression which was present only in RBP4+ ECs was significantly increased in the L group compared to the H group (p = 0.009, Wilcoxon rank-sum test, Figure 6H; Supplementary Figure S8D). Trajectory analysis demonstrated that RBP4+ECs persisted throughout the differentiation trajectory and showed highly expressed HMGB1 in the L group along with the pseudotime (Figures 6I–K). Moreover, CytoTRACE analysis suggested that RBP4+ECs in the L group showed a less differentiated status (p = 0.0079, Wilcoxon rank-sum test, Supplementary Figure S8E–G).
To further investigate the regulatory mechanisms of RBP4+ECs on T cell activation, we identified regulators from RBP4+ECs for T cell activation through Nichenet analysis () (Figure 7A). We found that RBP4+ECs regulated the expression of 29 target genes in T cells by upregulating HMGB1 in the L group, including IFNG, which enhances tumor immune function and promotes T cell exhaustion (), and CCL3, which facilitates dendritic cells (DC) recruitment for antigen presentation and T cell activation (). Interaction analysis between RBP4+ECs and T cells showed more and stronger chemokine-mediated interactions between RBP4+ECs and CD8+T in the L group including CXCR6-CXCL16, NR3C1-CCL2, CXCL12-CXCR3, CXCL12-CXCR4, and PGRMC2-CCL4L2 (Figure 7B). Similar results were observed in the interactions between RPB4+ECs and CD4+ T cells (Figure 7C).
FIGURE 7
Discussion
Tumor heterogeneity pervades across various cancers, profoundly influencing the survival of patients (). In this study, we observed sample-specific clustering patterns among malignant cells, whereas non-malignant cells exhibited clustering patterns predominantly defined by cell types. The correlation analysis of malignant cells revealed that intra-tumor heterogeneity was variable among samples and intra-tumor cells similarity was generally greater than that of inter-tumor cells. Previous research has identified complex transcriptional programs across diverse malignancies, including cell cycle, interferon response, partial epithelial-mesenchymal transition (pEMT), and complete EMT (cEMT) (). Here, we identified three meta-programs that were associated with distinct hallmark functions using NMF. Several strategies have been proposed to quantify the degree of tumor heterogeneity. For instance, employed PCA to project all malignant cells to the eigenvector space to derive PCs and calculated the averaged distance of malignant cells from the centroid (arithmetic mean of PCs of all malignant cells within the corresponding tumor) in principal component space. This averaged distance was used as a measure of heterogeneity for that particular sample. Guo et al. first calculated inter cell distances based on the Pearson correlation between any two malignant cells within a sample. The distance between cells was defined as D = 1- R, where R is the Pearson correlation coefficient between the CNV profiles of any two cells. Then, the metric of heterogeneity was defined as the median distance of distribution of pairwise cell distances (). However, these methods implicitly involved the cell number in the calculation of heterogeneity, thus making them potentially influenced by the number of cells. Therefore, we proposed a new metric that avoided the involvement and interference of cell number and accurately quantified the complexity of subclonal structure without bias from cell numbers.
Data originating from glioblastoma (; ), oligodendroglioma (), head and neck cancer (), and melanoma () consistently demonstrated a spectrum of differentiation within tumors, ranging from stem-like or progenitor states to fully differentiated cells, alongside a plethora of heterogeneous functional states in malignant cells, such as stress responses, interferon reactions, and hypoxia (; ; ; ; ). We identified a total of 37 significantly different cancer functional states between low and high-complexity groups. The low-complexity group was significantly enriched in functional states associated with hypoxia, invasiveness, metastasis, stemness, apoptosis, and DNA repair, suggesting their potential roles in tumor progression and therapeutic resistance, driven by highly competitive dominant clones (; ). Furthermore, an investigation into lineage origins revealed that distinct compositions of cellular lineages in samples with varying subclonal complexity. Additionally, the TFs regulatory networks associated with the low and high-complexity groups reveal a correlation between subclonal complexity and survival outcomes in HCC.
Previous research has established that tumor heterogeneity is a critical cornerstone for reconstructing tumor evolution and driving the dynamics and diversity of the tumor microenvironment (). Genetic and epigenetic alterations collectively define the transcriptomic and phenotypic heterogeneity of malignant cells, leading directly or indirectly to the reprogramming of the TME (; ; ). In this study, we observed markedly distinct interaction strength between malignant cells and other cell types within the microenvironment in high- versus low-complexity group. Specifically, malignant cells exhibited a substantially higher number of ligand-receptor interactions with myeloid cells, fibroblasts, and endothelial cells in low-complexity group compared to the high-complexity group. Notably, endothelial cells displayed the largest differences between the two groups, suggesting a potential role for endothelial cells in shaping the TME according to the complexity of subclonal structure. Our subsequent analysis revealed a remodeled immunological landscape in the low-complexity group, characterized by enhanced antigen-presenting capabilities of macrophages and heightened T cells immunoreactivity.
Cytokines play an important role in remodeling the tumor microenvironment (). In this study, we identified HMGB1, a critical cytokine expressed by RBP4+ECs, which exerts immunostimulatory effects on T cells. Our results demonstrated that both the expression and signaling activity of HMGB1 were significantly higher in the low complexity group compared to the high complexity group. A notable correlation was observed between the expression levels of HMGB1 in endothelial cells and immune activity score of T cells. We further confirmed that RBP4+ECs had more frequent chemokine-mediated ligand-receptor interaction pairs with T-cells in the low complexity group. These findings collectively highlighted the intricate interplay between the subclonal complexity and the TME, particularly emphasizing the role of subclonal complexity in orchestrating immunological landscape. The increased interaction strength between RBP4+ECs and T cells and activation of T cells in the low complexity group, mediated by HMGB1, suggests a potential mechanism through which the cytokines can either promote or suppress antitumor immunity. Future studies should focus on elucidating the molecular mechanisms underlying the varying expression of HMGB1 and its impact on T cells, potentially leading to the identification of novel biomarkers and therapeutic targets for cancer treatment.
Statements
Data availability statement
Publicly available datasets were analyzed in this study. This data can be found here: GSE156625 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE156625 GSE149614 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE149614 GSE151530 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE151530 GSE189903 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE189903.
Ethics statement
Ethical approval was not required for the study involving humans in accordance with the local legislation and institutional requirements. Written informed consent to participate in this study was not required from the participants or the participants’ legal guardians/next of kin in accordance with the national legislation and the institutional requirements.
Author contributions
SJ: Writing–original draft, Conceptualization, Data curation, Formal Analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing–review and editing. YZ: Data curation, Formal Analysis, Software, Visualization, Writing–review and editing. LX: Conceptualization, Data curation, Formal analysis, Visualization, Writing–review and editing. FW: Writing–original draft, Writing–review and editing, Conceptualization.
Funding
The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. This work was supported by the National Natural Science Foundation of China (Grant Nos. 32300543).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2024.1467682/full#supplementary-material
References
1
AibarS.Gonzalez-BlasC. B.MoermanT.Huynh-ThuV. A.ImrichovaH.HulselmansG.et al (2017). SCENIC: single-cell regulatory network inference and clustering. Nat. Methods14 (11), 1083–1086. 10.1038/nmeth.4463
2
AlizadehA. A.ArandaV.BardelliA.BlanpainC.BockC.BorowskiC.et al (2015). Toward understanding and exploiting tumor heterogeneity. Nat. Med.21 (8), 846–853. 10.1038/nm.3915
3
BarkleyD.MoncadaR.PourM.LibermanD. A.DrygI.WerbaG.et al (2022). Cancer cell states recur across tumor types and form specific interactions with the tumor microenvironment. Nat. Genet.54 (8), 1192–1201. 10.1038/s41588-022-01141-9
4
BenciJ. L.JohnsonL. R.ChoaR.XuY.QiuJ.ZhouZ.et al (2019). Opposing functions of interferon coordinate adaptive and innate immune responses to cancer immune checkpoint blockade. Cell178 (4), 933–948. 10.1016/j.cell.2019.07.019
5
BinnewiesM.RobertsE. W.KerstenK.ChanV.FearonD. F.MeradM.et al (2018). Understanding the tumor immune microenvironment (TIME) for effective therapy. Nat. Med.24 (5), 541–550. 10.1038/s41591-018-0014-x
6
BrowaeysR.SaelensW.SaeysY. (2020). NicheNet: modeling intercellular communication by linking ligands to target genes. Nat. Methods17 (2), 159–162. 10.1038/s41592-019-0667-5
7
Cancer Genome Atlas ResearchN. (2013). Comprehensive molecular characterization of clear cell renal cell carcinoma. Nature499 (7456), 43–49. 10.1038/nature12222
8
CastellinoF.HuangA. Y.Altan-BonnetG.StollS.ScheineckerC.GermainR. N. (2006). Chemokines enhance immunity by guiding naive CD8+ T cells to sites of CD4+ T cell-dendritic cell interaction. Nature440 (7086), 890–895. 10.1038/nature04651
9
ChengS.LiZ.GaoR.XingB.GaoY.YangY.et al (2021). A pan-cancer single-cell transcriptional atlas of tumor infiltrating myeloid cells. Cell184 (3), 792–809 e23. 10.1016/j.cell.2021.01.010
10
Dagogo-JackI.ShawA. T. (2018). Tumour heterogeneity and resistance to cancer therapies. Nat. Rev. Clin. Oncol.15 (2), 81–94. 10.1038/nrclinonc.2017.166
11
DavisA.GaoR.NavinN. (2017). Tumor evolution: linear, branching, neutral or punctuated?Biochim. Biophys. Acta Rev. Cancer1867 (2), 151–161. 10.1016/j.bbcan.2017.01.003
12
DentroS. C.LeshchinerI.HaaseK.TarabichiM.WintersingerJ.DeshwarA. G.et al (2021). Characterizing genetic intra-tumor heterogeneity across 2,658 human cancer genomes. Cell184 (8), 2239–2254.e39. 10.1016/j.cell.2021.03.009
13
DunnG. P.BruceA. T.IkedaH.OldL. J.SchreiberR. D. (2002). Cancer immunoediting: from immunosurveillance to tumor escape. Nat. Immunol.3 (11), 991–998. 10.1038/ni1102-991
14
EfremovaM.Vento-TormoM.TeichmannS. A.Vento-TormoR. (2020). CellPhoneDB: inferring cell-cell communication from combined expression of multi-subunit ligand-receptor complexes. Nat. Protoc.15 (4), 1484–1506. 10.1038/s41596-020-0292-x
15
FlavahanW. A.GaskellE.BernsteinB. E. (2017). Epigenetic plasticity and the hallmarks of cancer. Science357 (6348), eaal2380. 10.1126/science.aal2380
16
GaoJ.AksoyB. A.DogrusozU.DresdnerG.GrossB.SumerS. O.et al (2013). Integrative analysis of complex cancer genomics and clinical profiles using the cBioPortal. Sci. Signal6 (269), pl1. 10.1126/scisignal.2004088
17
GaujouxR.SeoigheC. (2010). A flexible R package for nonnegative matrix factorization. BMC Bioinforma.11, 367. 10.1186/1471-2105-11-367
18
GoveiaJ.RohlenovaK.TavernaF.TrepsL.ConradiL. C.PircherA.et al (2020). An integrated gene expression landscape profiling approach to identify lung tumor endothelial cell heterogeneity and angiogenic candidates. Cancer Cell37 (1), 421–436 e13. 10.1016/j.ccell.2020.03.002
19
GreavesM.MaleyC. C. (2012). Clonal evolution in cancer. Nature481 (7381), 306–313. 10.1038/nature10762
20
GulatiG. S.SikandarS. S.WescheD. J.ManjunathA.BharadwajA.BergerM. J.et al (2020). Single-cell transcriptional diversity is a hallmark of developmental potential. Science367 (6476), 405–411. 10.1126/science.aax0249
21
GuoL.YiX.ChenL.ZhangT.GuoH.ChenZ.et al (2022). Single-cell DNA sequencing reveals punctuated and gradual clonal evolution in hepatocellular carcinoma. Gastroenterology162 (1), 238–252. 10.1053/j.gastro.2021.08.052
22
HanX.ZhouZ.FeiL.SunH.WangR.ChenY.et al (2020). Construction of a human cell landscape at single-cell level. Nature581 (7808), 303–309. 10.1038/s41586-020-2157-4
23
HanzelmannS.CasteloR.GuinneyJ. (2013). GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinforma.14, 7. 10.1186/1471-2105-14-7
24
HaoY.HaoS.Andersen-NissenE.MauckW. M.3rdZhengS.ButlerA.et al (2021). Integrated analysis of multimodal single-cell data. Cell184 (13), 3573–3587.e29. 10.1016/j.cell.2021.04.048
25
JiangP.ZhangY.RuB.YangY.VuT.PaulR.et al (2021). Systematic investigation of cytokine signaling activity at the tissue and single-cell levels. Nat. Methods18 (10), 1181–1191. 10.1038/s41592-021-01274-5
26
KashimaY.SakamotoY.KanekoK.SekiM.SuzukiY.SuzukiA. (2020). Single-cell sequencing techniques from individual to multiomics analyses. Exp. Mol. Med.52 (9), 1419–1427. 10.1038/s12276-020-00499-2
27
KhaliliJ. S.LiuS.Rodriguez-CruzT. G.WhittingtonM.WardellS.LiuC.et al (2012). Oncogenic BRAF(V600E) promotes stromal cell-mediated immunosuppression via induction of interleukin-1 in melanoma. Clin. Cancer Res.18 (19), 5329–5340. 10.1158/1078-0432.CCR-12-1632
28
KnocheS. M.LarsonA. C.SlikerB. H.PoelaertB. J.SolheimJ. C. (2021). The role of tumor heterogeneity in immune-tumor interactions. Cancer Metastasis Rev.40 (2), 377–389. 10.1007/s10555-021-09957-3
29
KorsunskyI.MillardN.FanJ.SlowikowskiK.ZhangF.WeiK.et al (2019). Fast, sensitive and accurate integration of single-cell data with Harmony. Nat. Methods16 (12), 1289–1296. 10.1038/s41592-019-0619-0
30
LeaderA. M.GroutJ. A.MaierB. B.NabetB. Y.ParkM. D.TabachnikovaA.et al (2021). Single-cell analysis of human non-small cell lung cancer lesions refines tumor classification and patient stratification. Cancer Cell39 (12), 1594–1609.e12. 10.1016/j.ccell.2021.10.009
31
LiY.JiangT.ZhouW.LiJ.LiX.WangQ.et al (2020). Pan-cancer characterization of immune-related lncRNAs identifies potential oncogenic biomarkers. Nat. Commun.11 (1), 1000. 10.1038/s41467-020-14802-2
32
LiberzonA.BirgerC.ThorvaldsdottirH.GhandiM.MesirovJ. P.TamayoP. (2015). The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst.1 (6), 417–425. 10.1016/j.cels.2015.12.004
33
LiuR. B.EngelsB.SchreiberK.CiszewskiC.SchietingerA.SchreiberH.et al (2013). IL-15 in tumor microenvironment causes rejection of large established tumors by T cells in a noncognate T cell receptor-dependent manner. Proc. Natl. Acad. Sci. U. S. A.110 (20), 8158–8163. 10.1073/pnas.1301022110
34
LiuY.XunZ.MaK.LiangS.LiX.ZhouS.et al (2023). Identification of a tumour immune barrier in the HCC microenvironment that determines the efficacy of immunotherapy. J. Hepatol.78 (4), 770–782. 10.1016/j.jhep.2023.01.011
35
LongZ.SunC.TangM.WangY.MaJ.YuJ.et al (2022). Single-cell multiomics analysis reveals regulatory programs in clear cell renal cell carcinoma. Cell Discov.8 (1), 68. 10.1038/s41421-022-00415-0
36
LuY.YangA.QuanC.PanY.ZhangH.LiY.et al (2022). A single-cell atlas of the multicellular ecosystem of primary and metastatic hepatocellular carcinoma. Nat. Commun.13 (1), 4594. 10.1038/s41467-022-32283-3
37
MaL.HeinrichS.WangL.KeggenhoffF. L.KhatibS.ForguesM.et al (2022). Multiregional single-cell dissection of tumor and immune cells reveals stable lock-and-key features in liver cancer. Nat. Commun.13 (1), 7533. 10.1038/s41467-022-35291-5
38
MaL.HernandezM. O.ZhaoY.MehtaM.TranB.KellyM.et al (2019). Tumor cell biodiversity drives microenvironmental reprogramming in liver cancer. Cancer Cell36 (4), 418–430. 10.1016/j.ccell.2019.08.007
39
MaL.WangL.KhatibS. A.ChangC. W.HeinrichS.DominguezD. A.et al (2021). Single-cell atlas of tumor cell evolution in response to therapy in hepatocellular carcinoma and intrahepatic cholangiocarcinoma. J. Hepatol.75 (6), 1397–1408. 10.1016/j.jhep.2021.06.028
40
MarusykA.JaniszewskaM.PolyakK. (2020). Intratumor heterogeneity: the rosetta stone of therapy resistance. Cancer Cell37 (4), 471–484. 10.1016/j.ccell.2020.03.007
41
McGinnisC. S.MurrowL. M.GartnerZ. J. (2019). DoubletFinder: doublet detection in single-cell RNA sequencing data using artificial nearest neighbors. Cell Syst.8 (4), 329–337. 10.1016/j.cels.2019.03.003
42
MoncadaR.BarkleyD.WagnerF.ChiodinM.DevlinJ. C.BaronM.et al (2020). Integrating microarray-based spatial transcriptomics and single-cell RNA-seq reveals tissue architecture in pancreatic ductal adenocarcinomas. Nat. Biotechnol.38 (3), 333–342. 10.1038/s41587-019-0392-8
43
NamA. S.ChaligneR.LandauD. A. (2021). Integrating genetic and non-genetic determinants of cancer evolution by single-cell multi-omics. Nat. Rev. Genet.22 (1), 3–18. 10.1038/s41576-020-0265-5
44
NeftelC.LaffyJ.FilbinM. G.HaraT.ShoreM. E.RahmeG. J.et al (2019). An integrative model of cellular states, plasticity, and genetics for glioblastoma. Cell178 (4), 835–849. 10.1016/j.cell.2019.06.024
45
O’DonnellJ. S.TengM. W. L.SmythM. J. (2019). Cancer immunoediting and resistance to T cell-based immunotherapy. Nat. Rev. Clin. Oncol.16 (3), 151–167. 10.1038/s41571-018-0142-8
46
PatelA. P.TiroshI.TrombettaJ. J.ShalekA. K.GillespieS. M.WakimotoH.et al (2014). Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma. Science344 (6190), 1396–1401. 10.1126/science.1254257
47
PolyakK.HavivI.CampbellI. G. (2009). Co-evolution of tumor cells and their microenvironment. Trends Genet.25 (1), 30–38. 10.1016/j.tig.2008.10.012
48
Project. iotTC (2019). Available at: https://github.com/broadinstitute/inferCNV.
49
PuramS. V.TiroshI.ParikhA. S.PatelA. P.YizhakK.GillespieS.et al (2017). Single-cell transcriptomic analysis of primary and metastatic tumor ecosystems in head and neck cancer. Cell171 (7), 1611–1624. 10.1016/j.cell.2017.10.044
50
QiuX.MaoQ.TangY.WangL.ChawlaR.PlinerH. A.et al (2017). Reversed graph embedding resolves complex single-cell trajectories. Nat. Methods14 (10), 979–982. 10.1038/nmeth.4402
51
RambowF.RogiersA.Marin-BejarO.AibarS.FemelJ.DewaeleM.et al (2018). Toward minimal residual disease-directed therapy in melanoma. Cell.174 (4), 843–855. 10.1016/j.cell.2018.06.025
52
RobinX.TurckN.HainardA.TibertiN.LisacekF.SanchezJ. C.et al (2011). pROC: an open-source package for R and S+ to analyze and compare ROC curves. BMC Bioinforma.12, 77. 10.1186/1471-2105-12-77
53
RooneyM. S.ShuklaS. A.WuC. J.GetzG.HacohenN. (2015). Molecular and genetic properties of tumors associated with local immune cytolytic activity. Cell160 (1-2), 48–61. 10.1016/j.cell.2014.12.033
54
SchuppJ. C.AdamsT. S.CosmeC.Jr.RaredonM. S. B.YuanY.OmoteN.et al (2021). Integrated single-cell atlas of endothelial cells of the human lung. Circulation144 (4), 286–302. 10.1161/CIRCULATIONAHA.120.052318
55
ShannonP.MarkielA.OzierO.BaligaN. S.WangJ. T.RamageD.et al (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res.13 (11), 2498–2504. 10.1101/gr.1239303
56
SharmaA.SeowJ. J. W.DutertreC. A.PaiR.BleriotC.MishraA.et al (2020). Onco-fetal reprogramming of endothelial cells drives immunosuppressive macrophages in hepatocellular carcinoma. Cell183 (2), 377–394. 10.1016/j.cell.2020.08.040
57
SunX. X.YuQ. (2015). Intra-tumor heterogeneity of cancer cells and its implications for cancer treatment. Acta Pharmacol. Sin.36 (10), 1219–1227. 10.1038/aps.2015.92
58
TiroshI.VenteicherA. S.HebertC.EscalanteL. E.PatelA. P.YizhakK.et al (2016). Single-cell RNA-seq supports a developmental hierarchy in human oligodendroglioma. Nature539 (7628), 309–313. 10.1038/nature20123
59
TsaiC. H.ChuangY. M.LiX.YuY. R.TzengS. F.TeohS. T.et al (2023). Immunoediting instructs tumor metabolic reprogramming to support immune evasion. Cell Metab.35 (1), 118–133.e7. 10.1016/j.cmet.2022.12.003
60
VasanN.BaselgaJ.HymanD. M. (2019). A view on drug resistance in cancer. Nature575 (7782), 299–309. 10.1038/s41586-019-1730-1
61
VitaleI.ShemaE.LoiS.GalluzziL. (2021). Intratumoral heterogeneity in cancer progression and response to immunotherapy. Nat. Med.27 (2), 212–224. 10.1038/s41591-021-01233-9
62
VoitR. A.TaoL.YuF.CatoL. D.CohenB.FlemingT. J.et al (2023). A genetic disorder reveals a hematopoietic stem cell regulatory network co-opted in leukemia. Nat. Immunol.24 (1), 69–83. 10.1038/s41590-022-01370-4
63
WangW.GuoY.ZhongJ.WangQ.WangX.WeiH.et al (2021). The clinical significance of microvascular invasion in the surgical planning and postoperative sequential treatment in hepatocellular carcinoma. Sci. Rep.11 (1), 2415. 10.1038/s41598-021-82058-x
64
WangY.ZhangM.ShiJ.ZhuY.WangX.ZhangS.et al (2023). Cracking the pattern of tumor evolution based on single-cell copy number alterations. Brief. Bioinform24 (6), bbad341. 10.1093/bib/bbad341
65
WuF.FanJ.HeY.XiongA.YuJ.LiY.et al (2021). Single-cell profiling of tumor heterogeneity and the microenvironment in advanced non-small cell lung cancer. Nat. Commun.12 (1), 2540. 10.1038/s41467-021-22801-0
66
YuF.CatoL. D.WengC.LiggettL. A.JeonS.XuK.et al (2022). Variant to function mapping at single-cell resolution through network propagation. Nat. Biotechnol.40 (11), 1644–1653. 10.1038/s41587-022-01341-y
67
YuF.SankaranV. G.YuanG. C. (2021). CUT&RUNTools 2.0: a pipeline for single-cell and bulk-level CUT&RUN and CUT&Tag data analysis. Bioinformatics38 (1), 252–254. 10.1093/bioinformatics/btab507
68
YuG.WangL. G.HanY.HeQ. Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS16 (5), 284–287. 10.1089/omi.2011.0118
69
YuanH.YanM.ZhangG.LiuW.DengC.LiaoG.et al (2019). CancerSEA: a cancer single-cell state atlas. Nucleic Acids Res.47 (D1), D900-D908–D8. 10.1093/nar/gky939
70
ZhangA.MiaoK.SunH.DengC. X. (2022a). Tumor heterogeneity reshapes the tumor microenvironment to influence drug resistance. Int. J. Biol. Sci.18 (7), 3019–3033. 10.7150/ijbs.72534
71
ZhangY.VuT.PalmerD. C.KishtonR. J.GongL.HuangJ.et al (2022b). A T cell resilience model associated with response to immunotherapy in multiple tumor types. Nat. Med.28 (7), 1421–1431. 10.1038/s41591-022-01799-y
72
ZhengL.QinS.SiW.WangA.XingB.GaoR.et al (2021). Pan-cancer single-cell landscape of tumor-infiltrating T cells. Science374 (6574), abe6474. 10.1126/science.abe6474
Summary
Keywords
intra-tumoral heterogeneity, tumor microenvironment, subclonal structure, immune activation, cell-cell communications
Citation
Shi J, Zhang Y, Xu L and Wang F (2024) Single-cell transcriptomics reveals tumor microenvironment remodeling in hepatocellular carcinoma with varying tumor subclonal complexity. Front. Genet. 15:1467682. doi: 10.3389/fgene.2024.1467682
Received
20 July 2024
Accepted
19 August 2024
Published
29 August 2024
Volume
15 - 2024
Edited by
Zhe Lei, The First Affiliated Hospital of Soochow University, China
Reviewed by
Yanshuo Chu, University of Texas MD Anderson Cancer Center, United States
Chuan-Xing Li, Karolinska Institutet, Sweden
Updates
Copyright
© 2024 Shi, Zhang, Xu and Wang.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Jian Shi, shijian27@mail.sysu.edu.cn; Fang Wang, wangf359@mail.sysu.edu.cn
† These authors share first authorship
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.