The patient’s liver tumor and adjacent tissues were obtained from a patient with clinically confirmed HCC. The patient was a male, 52 years old. Para-tumor liver tissue was collected less than 1 cm from the tumor, and distant normal liver tissue was collected more than 1 cm from the tumor. Samples were collected with the approval from the Ethics Committee of the Affiliated Cancer Hospital and Institute of Guangzhou Medical University. Informed consent was signed by the patient.

The obtained tissue samples were washed two to three times in 4°C precooled 1 × PBS on ice, transferred to a sterile RNase-free culture dish, and cut into small pieces (around 0.5 mm²) with 10-cm surgical scissors. The tissues were washed with 1 × PBS. Non-purpose tissues such as blood stains, fatty layers, and connective tissues were removed as much as possible. The tissue samples were incubated in a constant temperature water bath or shaker at 37°C with the self-prepared hydrolysate or the hydrolysate kit (MACS Human Tumor Dissociation Kit DS_130-095-929). Incubation was terminated when the tissue mass disappeared. A 40-μm cell sieve was used to filter. The cell suspension was then centrifuged at 300 rpm at 4°C for 5 min. The supernatant was discarded, and 1 ml of resuspension buffer was added to the centrifuge tube to resuspend the cells. A measure of 3 ml of precooled erythrocyte lysis liquid was added to the tube (Solarbio), blowing and performing suction evenly. It is then incubated at 4°C for 5–10 min and centrifuged at 300 rpm at 4°C for 5 min immediately after completion. The supernatant is removed, and 1 ml of resuspension buffer (precooled 1 × PBS) was added to fully resuspend the cells. Testing was conducted with a cell counter (Countstar Rigel S2), and the concentration was adjusted according to the test result. The target concentration is 700–1200 cell/μL. The cell activity should be greater than 90%, and the agglomeration should be less than 15%. Once the final cell concentration and activity were reached, the cells were placed on ice, and the 10× genomics single-cell transcriptome chip-on-board experiment was carried out within 30 min.

We prepared single-cell RNA-seq libraries with Chromium Next GEM Single Cell 3′ Reagent Kits v.3.1 on the Chromium Controller (10× Genomics). We prepared single-cell suspensions from cultured cell lines, and single cells were suspended in PBS containing 0.04% BSA. The cell suspension was loaded onto the Chromium Next GEM Chip G and ran the Chromium Controller to generate single-cell gel beads in the emulsion (GEMs), according to the manufacturer’s recommendation. The captured cells were lysed, and the released RNA was barcoded through reverse transcription in individual GEMs. Barcoded, full-length cDNA was generated, and libraries were constructed, according to the performer’s protocol. The quality of libraries was assessed by Qubit 4.0 and the Agilent 2100 systems. Sequencing was performed on the Illumina NovaSeq 6000 system with a sequencing depth of at least 50,000 reads per cell and 150 bp (PE150) paired-end reads.

We performed alignment to this amended reference using 10× cellranger 6.0.0, which employs the STAR sequence aligner. The reference genome was the human genome GRCh38. Overall, 12,784 cells passed the quality control. To remove the cells with low quality, cells with gene numbers over 500 and the ratio of mitochondria lower than 0.20 were maintained, and genes with at least one feature count in more than three cells were used for the following analysis. We determined the gene expression contusion unique molecular identifiers (UMIs) for each cell barcode–gene combination. Following alignment, we filtered cell barcodes to identify those which contain cells using the approach implemented in cellranger 3.0.2, and only these barcodes were considered for downstream analysis including clustering and cell-type identification and differential expression analysis by Seurat (version 4.0.1).

To enable unsupervised clustering and cell-type identification, we perform dimensionality reduction with principal component analysis (PCA) on the combined set of samples for each tissue. To visualize the data, we further reduced the dimensionality of all 26,231 cells using Seurat and also used t-SNE to project the cells into 2D space. The steps include 1. using the LogNormalize method of the “Normalization” function of Seurat software to calculate the expression value of genes; 2. PCA (principal component analysis) was performed using the normalized expression value, within all the PCs, and the top 10 PCs were used to perform clustering and t-SNE analysis; 3. to find clusters, selecting the weighted shared nearest neighbor (SNN) graph-based clustering method. Marker genes for each cluster were identified by the “bimod” function (likelihood-ratio test) with default parameters via the FindAllMarkers function in Seurat (4.0.1). Fold Change > 1.5 and FDR < 0.1 are filtered from the calculation results of the FindAllMarkers function, and then they are sorted to the top 10 genes as the marker gene.

We applied Louvain community detection to the nearest neighbor graph constructed in PCA space to define a cluster partition. To infer cell types, we trained a neural network classifier to predict cell ontology classes, given single-cell RNA-seq mRNA abundance profiles. First, we used SingleR (4.0.4) containing seven data sets to automate the identification of cell types. In addition to these annotations, we manually added cell state annotations to SingleR to provide a level of granularity below cell ontology classes. After completing the annotations, we further corrected the cell types using the CellMarker database (http://bio-bigdata.hrbmu.edu.cn/CellMarker/).

The gene subcellular location data were downloaded from the Human Protein Atlas (https://www.proteinatlas.org/about/download) database. TCGA RNA-seq mRNA expression data on LIHC were obtained from the LIHC project of TCGA (https://tcgadata.nci.nih.gov/tcga/) database. The data were downloaded from UCSC Xena (https://xena.ucsc.edu/). The LIHC project contains 377 primary liver cancer tissue samples and 50 normal liver tissue samples. The LIHC samples were divided into high- and low-expression groups based on a gene median expression level. Further mRNA expression data and clinical pathological data were obtained from the ICGC portal (https://dcc.icgc.org/projects/LIRI-JP). These samples belong to a Japanese population primarily infected with HBV/HCV. We used the normalized read count values given in the gene expression file. The GEO datasets GSE25097 from the Sangerbox (http://sangerbox.com/) were downloaded and analyzed. GSE25097, based on Rosetta/Merck Human RSTA Affymetrix 1.0 microarray, comprises 268 HCC tumor, 243 adjacent non-tumor, 40 cirrhotic, and 6 healthy liver samples. The GTEx datasets of bulk tissue gene expression were downloaded from the Genotype-Tissue Expression Project (GTEx: https://www.gtexportal.org/home/).

The differential gene expression analysis among HCC cells in tumor tissue and normal hepatocytes in distal liver tissue was determined by single-cell sequencing. Genes were chosen with differential expression setting; the average of fold change (tumor/normal) ≥ 2 and p value ≤0.05, which encode proteins located in the cell membrane for subsequent analysis. The public databases’ differential expression profiles between tumor tissues and the normal liver tissues were generated based on the normalized expression value of RNA-seq data. Independent student’s t-test was used to compare the median expression level of two different groups. A one-way ANOVA test was used to compare the means between three subgroups. The test was performed in GraphPad Prism (La Jolla, CA, United States). Kaplan–Meier survival curves of the two risk groups were plotted, and the log-rank p value of the survival difference was calculated between them.

Clinically confirmed HCC tissue and paired para-tumor and distant normal liver tissues were collected from a 52-year-old male patient and dissociated into single cells for high-throughput sequencing. A total of 12,784 cells from HCC patients with liver cancer (4,285 cells), as well as from proximal (4,586 cells) and distal normal tissue (3,913 cells), were obtained. In total, 28 cell clusters were identified and shown with Uniform Manifold Approximation and Projection (UMAP) (Figure 1A). The distribution of cells in different samples is shown in Figure 1B.

SingleR (4.0.4) containing seven datasets was used for automatic identification of cell types. In addition to these annotations, we manually added cell-state annotations to SingleR to provide a level of granularity below cell ontology classes. As shown in Figure 2A, types of cells were identified, and most of them are immune cells, such as B cells, CD8⁺ T cells, dendritic cells, macrophages, monocytes, neutrophils, and natural killer cells (Supplementary Table S1). This finding indicated the existence of intensive immune infiltration in the tumor and para-tumor regions. In addition, fibroblasts and endothelial cells were also found in the tumor microenvironment. This indicated that the tumor mass is composed of multiple cell types with high heterogeneity. These cells may interact with hepatic cells and form a niche for tumor development. In particular, cells with hepatic origin including normal hepatocytes and HCC tumor cells were highlighted for further analysis (Figure 2B).

Cell markers for normal hepatocytes and hepatocellular carcinoma cells, as reported in the current available studies, as well as in the CellMarker database (http://biobigdata.hrbmu.edu.cn/CellMarker/), were used to extract cell populations with hepatic origins from the total cells. The UMAP plots showed the distribution of hepatic cells in different tissue samples (Figure 3A). The annotated hepatocytes and HCC cells are shown in Figure 3B. The relative expression and distribution of HCC tumor-specific marker genes (EPCAM, SOX9, AFP, KRT7, S100A6, and S100A11) and hepatic lineage marker genes (ALB, PCK1, FGG, FGA, TTR, CEBPB, APOB, CYP2E1, and APOE) are also shown in the UMAP plots (Figure 3C and Figure 3D). To identify the potential tumor-specific markers, gene expression was compared between normal hepatocytes within the distant normal liver tissue and tumor cells within the tumor tissue. A total of 769 genes were identified by differential expression analysis with Seurat (version 4.0.1).

By matching these genes to a list of genes encoding plasma membrane proteins obtained from the Human Protein Atlas (HPA) website, 110 genes were selected. A total of 42 genes with differential expression settings of average fold change (tumor/normal) ≥ 2 and p value ≤ 0.05 were retained for further analysis. The relative expression and clinical significance of 42 genes selected were investigated in TCGA database. The LIHC dataset from TCGA database with differential expression profiles between tumor tissues and normal liver tissues were generated based on the normalized expression value of RNA-seq data. Three genes including integrin subunit Beta 1 (ITGB1), melanoma cell adhesion molecule (MCAM), and protein phosphatase 1 regulatory subunit 16A (PPP1R16A) highly expressed in HCC and negatively correlated with patient survival were further identified (Figure 4A and Figure 4B). To dissect the gene expression pattern at single-cell resolution, both the average expression and the percentage expression of the selected signature genes were further analyzed in whole-cell populations. The average expression reflects the total gene expression level, while the percentage expression reflects the percentage of positive cells in the selected cell populations. As shown in Figure 4C, most of the signature genes were not only upregulated in individual cells but also expanded to a higher percentage of the whole-cell population in the tumor tissue. Conversely, the hepatic terminal differentiation markers were downregulated and have restricted expression only in a small percentage of cells in the tumor. These findings further confirmed that tumor cells harboring the signature genes are dominant clones during cancer evolution. The distribution of the three genes is shown in the UMAP plots (Figure 5A). The expression of the three genes was found to be upregulated both in tumor tissues and para-tumor liver tissues, compared with the distant normal liver tissues (Figure 5B).

To further validate the findings, the relative expression of the 3-gene signature was investigated in two other independent HCC cohorts including the LIRI-JP datasets and GEO datasets (GSE25097). The results showed that the three genes were all significantly upregulated in the two cohorts of HCC patients (Figures 6A,B). These findings further confirmed that the three genes encoding cell membrane proteins are consistently upregulated in HCC tumor cells. We noticed that the absolute expression of the three genes in distant normal liver tissues from our scRNA-seq data was very low. We further checked the expression of the three genes in normal tissues or organs from the Genotype-Tissue Expression (GTEx) database. The results showed that the expression of the three genes remains at a low level in most normal tissues (Supplementary Figure S1). Further design of antibodies or antibody-conjugated drugs recognizing the 3-gene signature might help accurately target and eliminate the tumor cells.