Identification of molecular pattern and prognostic risk model based on ligand-receptor pairs in liver cancer

Introduction: The tumor microenvironment of hepatocellular carcinoma is composed of multiple cells, and the interactive communication between cells drives tumor progression and characterizes the tumor. Communication between cells is mainly achieved through signal transduction between receptor ligands, and the rise of single-cell technology has made it possible to analyze the communication network between cells.

Methods: We applied a train of bioinformatic techniques and in vitro experiments. We analyzed the composition of the microenvironment of liver cancer by combining single-cell sequencing data and transcriptome sequencing data from liver cancer to construct molecular typing and risk models for LRs. Then, we analyzed association of it with prognosis, mutation, KEGG, tumor microenvironment (TME), immune infiltration, tumor mutational burden (TMB) and drug sensitivity in liver cancer. qPCR and was used to identify SLC1A5 expression in LIHC cell lines and CCK8, transwell and cell colony formation were performed to validate the function of SLC1A5. Meanwhile, we also performed polarization of macrophages.

Results: In this experiment, we found that liver cancer tissues are rich in immune and mesenchymal cells, and there is extensive signaling between individual cells, so we constructed molecular typing and risk models for LRs. Combining clinical data revealed significant differences in clinical characteristics, prognosis and mutated genes between the molecular typing of receptor-ligand pairs, as well as in sensitivity to drugs; similarly, there were significant prognostic differences between the risk models. There were also notable differences in activated signaling pathways, infiltrating immune cells and immune subtypes. Subsequently, we used siRNA to knock down SLC1A5 in hepatocellular carcinoma cells and found that cell proliferation, migration and invasion were diminished.

Conclusions: In conclusion, our LRs model may become a marker to guide clinical treatment and prognosis.

Introduction

HCC ranks as the sixth most prevalent new tumor in the world, but the second most common new death (1). Due to the insidious character of liver cancer, the detection rate of early-stage liver cancer is quite poor, while the highly malignant characteristics of liver cancer frequently result in the identification of liver cancer at an advanced stage. The prognosis for advanced liver cancer is terrible, with a lack of surgical opportunities and a high risk of metastasis within and outside the liver, resulting in worse quality of living and a shorter survival time for patients (2).

The tumor microenvironment, a hothouse of research in recent years, is composed of tumor cells, immune cells and stromal cells that together modulate tumor growth, metastasis, drug resistance and other properties (3). In the last few years, the rapid development of single-cell sequencing has made it possible for scientists to peek into the details of the tumor microenvironment (4). Through the use of gene expression profiles of individual cells, we can identify the specific roles played by different cell types in the tumor microenvironment, refining the function of cells to investigate the microenvironment in a specific target (5).

Cell-to-cell signaling moderates the function and state of the cell. With the continuous enhancement of precision in the research of tumors, receptor-ligand interactions to deliver signals between tumor tissues have attracted the attention of scientists (6). In particular, in the tumor microenvironment, a majority of the interactions between tumor cells, immune cells and stromal cells are mediated by receptors and ligands, of which investigations on PD1 and PD-L1 have made a significant contribution to clinical diagnosis and treatment (7).

Currently, there is no effective method to forecast the prognosis of patients with hepatocellular carcinoma, but receptor-ligand interactions could, up to some extent, anticipate the malignancy of the tumor and thus predict the prognosis of patients. Consequently, we developed a receptor-ligand pairs (LRs) model based on hepatocellular carcinoma to predict the risk of patients and ultimately improve their survival.

Materials and methods

Datasets

The single cell sequencing data (GSE146115) was obtained from the GEO database, the gene expression profiling data, clinical data and mutation data from TCGA and ICGC. The ICGC-LIHC sample (231) set was regarded as the external validation set, and the TCGA-LIHC sample (365) set as the training set.

scRNA-seq data analysis and cell type definition

“Seurat” R package was applied to analyze the expression matrix of single cells, and we screened for cells with the optimal number of genes expressed (50~20000) (8). The mitochondrial genes have been removed from the expression matrix. After controlling for average expression and dispersion relationships, all highly variable genes in single cells were identified. Subsequently, the highly variable genes were used to perform principal component analysis to identify significant principal components as a means of eliminating batch effects based on the “jackstraw” function. Cells were cohorted into 12 different cell types at a resolution of 0.5 using the “FindClusters” function. The “FindAllMarkers” function was applied to discern differentially expressed genes (DEGs). Additionally, a few traditional markers for defining cell subsets were gathered from earlier studies (Table 1) and manually annotated in accordance with marker expression.

TABLE 1

Table 1 The marker genes for the cells.

Cell to cell communication

The various cells of the tumor microenvironment interact with each other to exert tumor promoting or inhibiting effects through activation between various ligands and receptors. Cellular communication was accomplished through “cellphonedb”, a public database containing ligands, receptors, and their interactions, and by annotating the membrane, secreted, and peripheral proteins of each cell subgroup at various time points(9). We conducted the “cellphonedb” to unpack the matrix of cellular communication and we filtered the receptors and ligands that appeared in the expression matrix in TCGA.

Selection of receptor-ligand pairs

Cellular interactions depend on the simultaneous expression of receptors and ligands, and receptors can only be stimulated to mediate intercellular communication when the number of receptors and ligands is at a parallel elevation. We screened for LRs with receptor-ligand co-expression correlations greater than 0.3 (p < 0.05) and used these LRs for clustering to determine molecular types (10).

Molecular subtyping calculation

We took the sum of the gene expression of the receptor and ligand as the expression of the LR. We combined the receptor-ligand expression data with the prognostic data to filter out LRs with prognostic significance and then used the R package “ConensusClusterPlus” to generate consensus matrix, “euclidean” was selected as the distance metric for PAM algorithms. A random subset of the TCGA data was selected from the TCGA data, the size of the subset was 80% of the original data set and 500 replications were performed (11). The amount of clusters was varied ranging from 2 to 10 and the most appropriate number of clusters was determined through the calculation of the consensus matrix and the cumulative distribution function (CDF).

Presumption of drug sensitivity

The R package “pRRophetic” was applied to calculate the sensitivity of different drugs based on expression matrix (12). We determined the appropriate drugs for different classifications of patients based on the sensitivity of the drug, as well as the classification of the different patients.

Risk model

The R package “glmnet” was conducted to screen for LRs using lasso regression, combined with prognostic data, to build a risk model to predict the risk of patients by classifying them into different risk groups using the median risk score. First, we chose “cv.glmnet” fuction to filter the λ that minimised the discrepancy. Then, based on the value of λ taken at this point, we got to filter the best LRs to build the prognostic model. The risk score was calculated using the formula LR. Score = ∑β(i) ×Exp(i), where i refers to the LR pair, Exp represents the level of LR pair expression and beta is the coefficient of the LR pair in the model. The median value of the training set (TCGA) was selected as the truncation value (13). The R package “survival” was employed to depict survival curves and to compare survival differences between high and low risk groups. The R package “timeROC” to portray the ROC curves of risk scores and traditional prognostic indicators, and to calculate the AUC values to assess the accuracy of the prognosis prediction (14).

Function enrichment and analysis of mutations and immunity

The “maftools” package was used to visualize the mutation data. We presented the twenty genes with the most significant mutations and compared mutations in patients of different subtypes. We applied “clusterProfiler” to analysis KEGG pathway (15). The ‘hallmark’ gene set collection from the molecular signature database was used for pathway enrichment analysis. Immune cell infiltration assessment was carried out with the “ssGSEA” package (16).

In vitro experimental validation

The cells used in this experiment were obtained from the laboratory of the Department of Oncology, Tongji Hospital, Huazhong University of Science and Technology. LO2, SNU398, Huh7 and HLF were cultured in DMEM medium with 10% fetal bovine serum added. We extracted RNA from the cells, reverse transcribed them into DNA and qPCR detected the expression of SLC1A5 in the cells. Subsequently, SLC1A5 was knocked down in Huh7 and SNU398 cells using siRNA, and CCK8, cell colony formation and transwell were used to detect the proliferation, migration and invasion ability of the cells (17). Meanwhile, we also performed polarization of macrophages. Detailed experimental steps are in the Supplementary Data Sheet 1, and all experiments were repeated three times.

Statistical analysis

Dichotomous variables were tested using the chi-square test, survival analysis was performed using the log-rank test, and comparison between the two groups was performed using the Wilcoxon test. p<0.05 was considered to be statistically different.

Results

The single-cell transcriptome landscape of hepatocellular carcinoma

Figure 1 showed the overall design and flow chart of this study. Since the gene expression data from single cells excluded mitochondrial genes, we calculated the correlation between the number of unique molecular identifiers and mRNA, which suggests that a significant positive correlation was shown between the number of unique molecular identifiers and mRNA (Figure 2A). Gene numbers for the vast majority of cells are between 0 and 8000(Figure 2B). After filtering the cells, a total of 3200 cells were included in the subsequent analysis. The differential genes in the various cell types were calculated for a total of 12 cell types after normalizing the expression data and filtering the first 2000 highly variable genes for the subsequent principal component analysis (Supplementary Figure 1A, Figure 2C). Based on the marker genes in the different cells, we identified the cell types that needed to be labelled. The marker genes for the cells were derived from databases and earlier studies. As a result, six cells were identified, such as B cells, fibroblast cells, hepatocytes, myeloid cells, NK cells, T cells (Figure 2D).

FIGURE 1

Figure 1 The flowchart of the study design.

FIGURE 2

Figure 2 (A) Correlation between the number of unique molecular identifiers and mRNA (B) Violin plot of features (C) tSNE of cell profiles; different color blocks represent related cell clusters. (D) tSNE of cell profiles; different color blocks represent related sample sources. (E) LR interactions between different cell subsets. (F) Network overview for the interaction between different cell subsets.

Intercellular communication networks in hepatocellular carcinoma

The tumor microenvironment in solid tumors is composed primarily of stromal cells and immune cells, in addition to the tumor cells themselves. Within the microenvironment, various cells communicate with each other to transmit information to influence tumor progression. When speculating on intercellular communication, we used “cellphonedb” and found that tumor cells contact mainly with fibroblasts and fibroblasts with myeloid cells in liver cancer (Figure 2E). In Figure 2F, thicker lines indicate more interactions, and the numbers on the lines represent the number of interacting nodes.

Molecular typing based on ligand-receptor pairs

In cellular communication, receptor-ligand interactions play an essential role, and through “cellphonedb” simulations we obtained the corresponding LRs. Based on the TCGA expression profile and taking into account the need for synergistic expression of receptors and ligands, we chose LRs with co-expression R-values above 0.3 and P-values less than 0.05. 81 LRs were selected in total (Supplementary Table 1). The sum of the expression values of the receptor and the ligand took the place of the LRs’ expression values. Molecular typing analysis was executed on the sample set of TCGA using the “ConsensusClusterPlus”. Based on the CDF value (Figure 3A), splitting into two clusters was the preferred candidate when the k value was taken as 2 (Figures 3B, C). The ICGC validation dataset was subjected to the same data processing procedure as the TCGA training dataset, and the outcomes were identical. Similar distinctions between the two groups of patients were made based on the LR expression pattern (Supplementary Figures 2A–C).

FIGURE 3

Figure 3 (A) CDF curve of samples from TCGA cohort. (B) Delta area curve of samples from TCGA cohort. (C) TCGA clustering heatmap of samples from TCGA cohort when consensus k = 2. (D) Overall survival curves of molecular subtypes based on LR pairs.

Comparison of clinical information in different molecular subtypes

According to the results of the clustering, combined with the prognosis analysis of the patients, we analyzed the survival of the patients in both clusters and the survival curves depicted are shown in the figure. The log-rank test showed that there is a significant difference in survival between the two clusters, with patients in the A subtype having a better prognosis than those in the B subtype in TCGA (Figure 3D). The outcomes of the prognostic analyses in the ICGC validation set were agreed with those in the TCGA training set. However, the number of patients in B subtype was too small, perhaps as a result of the limited sample size of patients in ICGC, leading to a P value (0.052) of more than a little over 0.05 (Supplementary Figure 2D). Meanwhile, we went on to analyze the gene expression and clinical characteristics of patients in different clusters and found that the expression levels of risk genes were significantly higher in B subtype patients than in A subtype (Figure 4A). Regarding the clinical traits, we discovered that patients in the B subtype had greater tumor grade and TNM stage as well as more fatalities than those in the A subtype (Figures 4B–D).

FIGURE 4

Figure 4 (A) Heatmap of expression of LRs in molecular subtypes (B) Grade in molecular subtypes (C) Fustat in molecular subtypes (D) Stage in molecular subtypes (E) TMB in molecular subtypes (F) Somatic mutation variation analysis in A cluster in the TCGA-LIHC cohort. (G) Somatic mutation variation analysis in B cluster in the TCGA-LIHC cohort.

Mutational characteristics of different molecular subtypes

Genetic alterations in tumor cells, mutational inactivation of anti-oncogenes, amplified overexpression of oncogenes, etc., all of which contribute to tumor development, are intimately associated to tumor growth. When we analyzed the mutation data of patients in A and B subtype patients, we found that there were more mutations in the A subtype patients (Figures 4F, G). The mutations in the B subtype patients were mainly in the antioncogenes and the types of mutations were mainly missense mutations, such as TP53 (Figure 4G).

We then analyzed the tumor mutational burden of patients to predict the efficacy of immunotherapy in patients, and we found that the tumor mutational burden (TMB) was higher in A subtype patients than in B subtype patients (Figure 4E), which means that the effectiveness of immunotherapy may be better in A subtype patients than in B subtype patients.

Drug sensitivity of different molecular subtypes

We also attempted to see whether there were changes in drug sensitivity amongst the various patient clusters, in addition to examining differences in clinical features and gene expression. Using the R package “pRRophetic” to predict patient sensitivity to chemotherapeutic drugs, we analyzed prominent liver cancer drugs and found that sorafenib, a first-line drug for liver cancer, was more efficacious in B subtype patients. At the same time, we observed that some drugs were more effective in A subtype patients, such as “Bleomycin” (Figure 5A),” Doxorubicin” (Figure 5B), “Gemcitabine” (Figure 5C), “Mitomycin” (Figure 5D), “Paclitaxel” (Figure 5F), and conversely, some drugs were more targeted in B subtype, such as “Methotrexate” (Figure 5E), “Rapamycin” (Figure 5G), “Sorafenib” (Figure 5H), “Temsirolimus” (Figure 5I).

FIGURE 5

Figure 5 (A) Bleomycin. (B) Doxorubicin. (C) Gemcitabine. (D) Mitomycin. (E) Methotrexate. (F) Paclitaxel. (G) Rapamycin. (H) Sorafenib. (I) Temsirolimus.

Establishment of the risk model based on the ligand–receptor pair score

Molecular subtypes based on LRs had different clinicopathological characteristics as well as drug sensitivity. We sought to construct a prognostic model using LRs to assess patient risk. A total of 18 LRs associated with prognosis were considered (Supplementary Figure 1B), and we used LASSO cox regression to shrink the number of LRs, with optimal results occurring when λ was 0.0159 Figures 6A, B). We then used stepwise multivariate regression for optimization, and a total of nine LRs were screened out for model construction (Supplementary Table 2). In the training set (TCGA), the median risk score was selected as the cut-off value to classify patients into high and low risk groups, while patients in the validation set (ICGC) were similarly classified into high and low risk groups using the cut-off value in the training set as the boundary.

FIGURE 6

Figure 6 (A) Plotting of multinomial deviance versus log(λ). (B) LASSO coefficient profiles of the LRs. (C) Survival benefit of LR.score in the high and low LR.score groups in the TCGA-LIHC cohort. (D) Survival benefit of LR.score in the high and low LR.score groups in the ICGC-LIHC cohort. (E) The predictive value of LR.score in patients among the TCGA-LIHC cohort. (F) The predictive value of LR.score in patients among the ICGC-LIHC cohort. (G) Univariate cox regression analysis of LR.score, age, TNM stage and grade for overall survival (OS) in the TCGA-LIHC cohort. (H) Multivariate cox regression analysis of LR.score, and TNM stage for OS in the TCGA-LIHC cohort. (I) Univariate cox regression analysis of LR.score, age and TNM stage for overall survival (OS) in the ICGC-LIHC cohort. (J) Multivariate cox regression analysis of LR.score, and TNM stage for OS in the ICGC-LIHC cohort. (K) Fustat in the high and low LR.score groups in the TCGA-LIHC cohort. (L) Fustat in the high and low LR.score groups in the ICGC-LIHC cohort.

Correlation between the risk model and clinical features

Based on the risk groupings, we plotted the survival curves for the training and validation sets separately and used log-rank test to analyze whether there was a difference in survival, and we found that the high-risk group had a worse prognosis in both datasets (Figures 6C, D). We also plotted ROC curves based on risk scores and clinicopathological characteristics, and calculated AUC values for different factors separately, and observed that risk scores performed better in predicting patient risk in the training and validation sets (Figures 6E, F). In the TCGA training set, the risk model’s sensitivity was 0.784 and its specificity was 0.652; in the ICGC validation set, these values were 0.823 and 0.744, respectively. Both univariate and multivariate analyses showed that the risk score was an independent risk factor to assess patients’ prognosis in TCGA database (Figures 6G, H). In the validation set, we similarly found that risk score was also an independent prognostic factor in both univariate and multivariate analyses (Figures 6I, J).

The results of the principal components analysis showed that the grouping of high and low risk was able to separate the characteristics of the patients in TCGA (Supplementary Figure 3A), and the same phenomenon was observed in the validation set (Supplementary Figure 3B). Similarly, we found that the risk scores of surviving patients were significantly lower than those of the deceased in the training set (Figure 6K), and the same was witnessed in ICGC (Figure 6L). This suggested that the risk score in the prognostic model was highly effective in predicting patient survival, and that patients could be separated according to different risk groupings. In the training set, the proportion of patients who died increased as the risk score increased (Supplementary Figures 3C, D). Similarly, we observed the same phenomenon in the training set (Supplementary Figures 3E, F).

Pathway analysis and immune characterization of risk model

Additionally, we conducted an enrichment analysis of patients in the high-risk and low-risk groups and found that compared to low-risk patients, high-risk patients mainly activate a number of immune-related, metabolism-related and hypoxia-related pathways (Figures 7A, B).

FIGURE 7

Figure 7 (A) The results of the KEGG in TCGA-LIHC cohort. (B) The results of the KEGG in ICGC-LIHC cohort. (C) Analysis of immune cell scores in the TCGA-LIHC cohort calculated using the ssGSEA algorithms. (D) Comparison of immune subtypes in different risk groups (E) Correlation of LR.score and TIDE score (F) The expression level of CD274 (G) Alluvial diagram of the distribution of different gene cluster, risk score and survival outcome subtypes. *p< 0.05; **p< 0.01; ***p< 0.001.

We conducted the R package “ssGSEA” to analyze the immune infiltration of patients and we found a decrease in the infiltration of B cells, CD8+ T cells, mast cells, NK cells and, at the same time, an increase in the infiltration of activated DC cells and macrophages in high-risk patients (Figure 7C). The reduction of tumor-killing immune cells and addition of tumor-promoting immune cells may be one of the reasons for the worse prognosis in high-risk patients. Combined with the immune phenotyping of the patients, we found that the immune subtypes in the high-risk group were focused on C1 & C2 and in the low-risk group were mostly in C3 (Figure 7D).

We compared the levels of PD-L1 expression in the high and low risk groups and discovered that patients in the high-risk group had higher levels of PD-L1 expression (Figure 7F). This finding is also consistent with what we found of the ssGSEA above, which revealed that the high-risk group was primarily in an immunosuppressed state. As a result, the high-risk group may be more responsive to immunotherapeutic treatment targeting PD-1/PD-L1.

The TIDE scores are one of the basis for evaluating the current immunotherapy, so we also looked at the TIDE scores of the high and low risk groups. As a result, we discovered that there was a difference in the TIDE scores of the high and low risk groups (Figure 7E). Lower TIDE scores in the high-risk group imply that patients in the high-risk group will have a better immunotherapy outcome, which is consistent with the results of the PD-L1 expression level. Patients in subtype A were dominantly in the low-risk group and the state of survival, while the opposite was observed in subtype B patients (Figure 7G).

Functional experiments in vitro

In the LR model, the coefficient of LGALS9_SLC1A5 was the largest. Considering that SLC1A5 is a receptor, we took a series of experiments to test whether SLC1A5 affects the function of tumor cells. To clarify the expression of SLC1A5, we selected three hepatocellular carcinoma cells and normal hepatic epithelial cells to perform qPCR experiments. We found that the expression level of SLC1A5 was higher in all hepatocellular carcinoma cells than in normal hepatic epithelial cells, with the higher expression level in Huh7 and SNU398 cells (Supplementary Figure 1C). We performed the cell proliferation assay (CCK8 and cell colony formation), cell migration and invasion assay (transwell) using siRNA to knock down the expression level of SLC1A5 in Huh7 and SNU398 cells, and observed that the knockdown levels of si1 and si3 were the best among the three siRNAs (Figures 8A, B). The results indicated that the proliferation (Figures 8C–H), migration (Figures 9A–D) and invasion (Figures 9E–H) of hepatocellular carcinoma cells were significantly reduced after knockdown of SLC1A5.

FIGURE 8

Figure 8 (A) Validation of knockdown efficiency by qPCR in HuH7 cells. Results represent mean ± SD. n = 3. ****p < 0.0001; two-tailed t-test. (B) Validation of knockdown efficiency by qPCR in SNU398 cells. Results represent mean ± SD. n = 3. ****p < 0.0001; two-tailed t-test. (C) SLC1A5 knockdown inhibited colony formation of HuH7 cells. (D) SLC1A5 knockdown inhibited colony formation of SNU398 cells. (E) The colony number of HuH7 cells. Results represent mean ± SD; n = 3; ***p < 0.001; two-tailed t-test. (F) The colony number of SNU398 cells. Results represent mean ± SD; n = 3; ***p < 0.001; two-tailed t-test. (G) SLC1A5 siRNA displayed reduced proliferation of HuH7 cells. Results represent mean ± SD; n = 3; ***p < 0.001; two-tailed t-test. (H) SLC1A5 siRNA displayed reduced proliferation of SNU398 cells. Results represent mean ± SD; n = 3; ****p < 0.0001; two-tailed t-test.

FIGURE 9

Figure 9 (A) SLC1A5 knockdown inhibited migration of HuH7 cells. (B) SLC1A5 knockdown inhibited migration of SNU398 cells. (C) The cell number of migration of HuH7 cells. Results represent mean ± SD; n = 3; ****p < 0.0001; two-tailed t-test. (D) The cell number of migration of SNU398 cells. Results represent mean ± SD; n = 3; ****p < 0.0001; two-tailed t-test. (E) SLC1A5 knockdown inhibited invasion of HuH7 cells. (F) SLC1A5 knockdown inhibited invasion of SNU398 cells. (G) The cell number of invasion of HuH7 cells. Results represent mean ± SD; n = 3; ***p < 0.001; **p < 0.01; two-tailed t-test. (H) The cell number of invasion of SNU398 cells. Results represent mean ± SD; n = 3; ****p < 0.0001; two-tailed t-test. (I) Relative expression of gene markers of M2 (CD206, ARG1) macrophages by stimulation of different CM in THP1 cell lines.

After culturing macrophages using conditioned medium from tumor cells, we examined the expression of marker molecules in macrophages. Subsequently, we found that after culturing macrophages using conditioned medium from hepatocellular carcinoma cells with knockdown of SLC1A5, the expression of M2-type macrophage marker molecules was decreased, and the expression levels of CD206 and ARG1 were markedly down-regulated in macrophages (Figure 9I). This experimental phenomenon suggests that the expression of SLC1A5 by hepatocellular carcinoma cells may have an effect on the immune microenvironment.

Discussion

The tumor microenvironment is a burning research topic in the field of oncology. Immune cells are drawn to tumor-associated inflammatory changes, and these immune cells combine with stromal cells in the tumor tissue to form the tumor microenvironment (18). The tumor microenvironment is a double-edged sword in the process of tumor development. Tumor cells will be killed by immune cells to suppress tumor progression. Simultaneously, tumor cells will evade the surveillance of immune cells and gradually remodel the tumor microenvironment so that it promotes tumor metastasis and drug resistance (19).

According to the findings of single cell sequencing, fibroblasts, myeloid cells, B cells, T cells, and NK cells have been found infiltrating liver cancer tissue. The infiltration of T cell kills tumor cells, but recent research has revealed that T cell killing depends on non-depleted T cells and that tumor cells evade the immune system by causing T cell depletion (20, 21). Myeloid cells are predominantly macrophages and tumor cells can induce the conversion of macrophages into tumor-associated macrophages, which are able to interact with other cells either through direct contact or by secreting various effector molecules (22, 23). Tumor-associated fibroblasts contribute to tumor extracellular matrix remodeling, stemness characteristics, angiogenesis and drug resistance (24, 25).

Due to the sophisticated exploration of tumors, the treatment of tumors has evolved from targeting the tumor cells themselves in the early days, directly killing them through radiotherapy or chemotherapy (26), to currently targeting the tumor microenvironment and reversing the cancer-promoting microenvironment to eliminate the tumor cells. Immunotherapy targeting PD-1 on T cells interacting with PD-L1 on tumor cells has been widely used in clinical practice for a considerable period of time (27). Not all patients benefit from them, despite the fact that their efficacy is promising in some tumors (28). Immunotherapies targeting other destinations have also been explored in advance, and we were trying to identify other cellular interactions in the tumor microenvironment that could improve the prognosis of tumor patients. We therefore employed current single cell sequencing to analyze and identify cell types in hepatocellular carcinoma, utilizing software to predict cell-to-cell communication and sought out LRs that were prognostically relevant in hepatocellular carcinoma. Constructing molecular subtypes with different expression patterns of LRs, we observed significant differences in prognosis as well as clinicopathological features of patients in different subtypes, suggesting that we targeted these LRs as possible targets for future clinical therapy. At the same time, we observed that the sensitivity of patients to chemotherapeutic drugs differed between subtypes, which could be the basis for more targeted drug delivery to different patients in the clinical field.

The function of the 18 pairs of LRs for which we constructed molecular typing has been studied in a number of ways. For example, the IL15_IL15RA interaction produces a two-sided effect, with IL-15 promoting the proliferation and maintaining the survival of certain T cells as well as consistently promoting anti-tumor responses and being crucial for controlling tumor growth and metastasis in vivo. However, chronic inflammatory stimulation of IL-15 increases tumor growth and metastasis (29). The connection between SELP and CD34 shows that patients gain from the stimulation of the innate immune response to improve anti-tumor immunity, eliminate tumor cells, and hinder the growth of tumors (30). The interaction of the remaining 16 pairs of LRs, in contrast, is what is causing the progression of the tumor, and this interaction allows the tumor microenvironment to change in a way that is pro-cancer, causing M2 type macrophages to transform and inhibiting the activity of NK cells. The main mechanisms of interaction between tumor cells and tumour-associated macrophages are LGALS9-SLC1A5, CCR1_CCL23, CSF1R _CSF1, GRN_SORT1, CSF1_SIRPA, and SIRPA_CD47, which primarily induce macrophage differentiation to M2 type while inhibiting macrophage phagocytosis of tumor cells, cause immune escape, and promote tumor cell proliferation, invasion, and metastasis (31–35). Specifically, LAIR1_LILRB4, EPHB6_EFNB1, and LGALS9_CD44 have immune-suppressive effects. The interactions of LGALS9_CD44 promote the differentiation and maintenance of Treg cells (36). The protein encoded by the EFNB1 gene in EPHB6_EFNB1 is a ligand for type I membrane proteins and Eph-associated receptor tyrosine kinase, and its binding to the ligand primarily exerts suppressive immune effects(37). Tyrosine phosphatases SHP-1 and SHP-2 and/or inositol phosphatase SHIP, which are detrimental to immune activation and promote tumor growth, are recruited by LAIR1_LILRB4 activation (38). WNT7B_FZD1, ICAM1_SPN and EPHA2_EFNA5 are mainly involved in facilitating tumorigenesis and invasive metastasis. Among them, ICAM1_SPN plays an important role in cell-cell interactions, and circulating tumor cells with stem cell properties may be able to use the adhesion protein ICAM1 to promote the formation of circulating tumor clusters that migrate from the body’s primary tumor sites to other organs (39, 40); EPHA2 in EPHA2_EFNA5 induces inhibition of the focal adhesion kinase (FAK), extracellular regulatory protein kinases (ERK) and Akt phosphorylation, thereby regulating motility, viability and proliferation of a variety of malignant cell lines (41).

SLC1A5 is a mitochondrial glutamine transporter and glutamine regulates energy metabolism, signal transduction and redox status in cells (42). Previous studies have indicated that SLC1A5 might affect how immune cells behave and infiltrate tumor microenvironment, which can promote cancer. In breast cancer, SLC1A5 can promote tumor progression (43). Additionally, SLC1A5 accelerates the growth of lung and colorectal malignancies by forcing tumor cells to undergo metabolic reprogramming.

In hepatocellular carcinoma, we observed that the expression level of SLC1A5 was dramatically up-regulated, while the ability of hepatocellular carcinoma cells to proliferate was significantly down-regulated after SLC1A5 expression was silenced using siRNA. Additionally, hepatocellular carcinoma cells’ capacity for migration was noticeably suppressed. The results suggested that the malignancy of hepatocellular carcinoma cells increased when the ligand agonized SLC1A5, and targeting SLC1A5 could alleviate tumor progression in hepatocellular carcinoma patients in the future (44, 45).

In the prognostic model we constructed, there were disparities in immune cell infiltration between the high and low risk groups, with fewer anti-cancerous immune cells and more pro-cancerous immune cells in the high-risk group, which contributed to the worse prognosis in the high-risk group. Meanwhile, patients in the low-risk group for hepatocellular carcinoma were primarily clustered in the C3 immune subtype, and the C3 type was inflammatory, in line with the results of the previous analysis (46). Most of the patients with the molecular subtype A belonged to the low-risk category. Patients in subtype A had a higher tumor mutational burden, and in relation to the previous analysis, we hypothesized that immunotherapy was more effective in subtype A patients. Moreover, both in the validation and training sets, the high-risk group focused on activating tumor metabolism (47), cell cycle, hypoxia, and immune-related pathways compared to the low-risk group. Activation of these pathways tended to make tumor worse, implying a poorer prognosis.

The functions of the 9 pairs of LRs used to construct the prognostic models also vary, among which LGALS9_HAVCR acting primarily as a promoter of lymphocyte activation. As a co-stimulatory molecule during lymphocyte activation, HAVCR enhances the anti-tumor effects of lymphocyte and induces changes in the tumor microenvironment, leading to efficient anti-tumor immunity (48). KLRB1 in KLRB1_CLEC2D is the gene encoding human CD161. KLRB1 gene inactivation or antibody-mediated KLRB1 blockade enhances T cell-mediated glioma cell killing in vitro, and the CD161_CLEC2D pathway defines a potential target for immunotherapy of glioma and other human cancers (49). Therefore, increased expression of KLRB1_CLEC2D and LGALS9_HAVCR suggests that patients may have better efficacy of immunotherapy.

Considering the limitations piece, there were three flaws in this paper. Firstly, the studies in the article were retrospective studies that underwent analysis after data collection, the actual clinical value of which had not yet been ascertained in genuine clinical practice. Additionally, the article contained just a limited amount of data, subsequently a larger sample size is supposed to be incorporated into the study. Moreover, the data in the article were partially biased.

Conclusion

In conclusion, we have analyzed the interactions between cells in hepatocellular carcinoma, thereby establishing molecular subtypes of cellular communication as well as a prognostic model. The expression pattern of LRs may be able to predict the effect of chemotherapy and immunotherapy in patients with hepatocellular carcinoma and to forecast the prognosis of patients. Our findings highlight the clinical implications of LRs and provide a basis for subsequent clinical translation.

Data availability statement

Publicly available datasets were analyzed in this study. The datasets can be found here: GEO database (HYPERLINK "https://www.ncbi.nlm.nih.gov/geo/"Home - GEO - NCBI (nih.gov)), ICGC database (HYPERLINK "https://dcc.icgc.org/"Welcome | ICGC Data Portal) and TCGA database (HYPERLINK "https://portal.gdc.cancer.gov/"GDC (cancer.gov)). Other datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author contributions

All authors contributed to data analysis, drafting or revising the article, have agreed on the journal to which the article will be submitted, gave final approval of the version to be published, and agree to be accountable for all aspects of the work. All authors contributed to the article and approved the submitted version.

Funding

This study was supported by the National Natural Science Foundation of China (81372644).

Acknowledgments

We gratefully acknowledge contributions from Ze Zhang for giving us LIHC cell lines and Duo Xu for modification of images. Thanks for the members’ help in the process of data analysis.

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/fimmu.2023.1187108/full#supplementary-material

Supplementary Figure 1 | (A) Heatmap showing marker genes for each cluster. (B) The HR values of LRs. (C) Quantitative polymerase chain reaction (qPCR) displayed upregulation of SLC1A5 in liver cancer cells compared to normal cell line (LO2). Results represent mean ± SD; n = 3. ****p < 0.0001; ***p < 0.001; two-tailed t-test.

Supplementary Figure 2 | (A) CDF curve of samples from ICGC cohort. (B) Delta area curve of samples from ICGC cohort. (C) ICGC clustering heatmap of samples from ICGC cohort when consensus k = 2. (D) Overall survival curves of molecular subtypes based on LR pairs.

Supplementary Figure 3 | (A) The results of the PCA in TCGA-LIHC cohort. (B) The results of the PCA in ICGC-LIHC cohort. (C) Distribution of risk score (high or low) in TCGA. (D) Distribution of status (dead or alive) in TCGA. (E) Distribution of risk score (high or low) in ICGC. (F) Distribution of status (dead or alive) in ICGC.

Supplementary Table 1 | The selected LRs with R-values above 0.3 and P-values less than 0.05.

Supplementary Table 2 | LRs associated with prognosis.

Supplementary Data Sheet 1 | Detailed procedures of cell transfection, qPCR, CCK8, transwell, colony formation, western blotting, polarization of macrophages, collection of conditioned media and LR.score.

References

1. Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, et al. Global Cancer Statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin (2021) 71(3):209–49. doi: 10.3322/caac.21660

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Xu M, Yang L, Lin Y, Lu Y, Bi X, Jiang T, et al. Emerging nanobiotechnology for precise theranostics of hepatocellular carcinoma. J Nanobiotechnol (2022) 20(1):427. doi: 10.1186/s12951-022-01615-2

CrossRef Full Text | Google Scholar

3. Hinshaw DC, Shevde LA. The tumor microenvironment innately modulates cancer progression. Cancer Res (2019) 79(18):4557–66. doi: 10.1158/0008-5472.Can-18-3962

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Yan X, Xie Y, Yang F, Hua Y, Zeng T, Sun C, et al. Comprehensive description of the current breast cancer microenvironment advancements via single-cell analysis. J Exp Clin Cancer Res (2021) 40(1):142. doi: 10.1186/s13046-021-01949-z

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Cao H, Huang T, Dai M, Kong X, Liu H, Zheng Z, et al. Tumor microenvironment and its implications for antitumor immunity in cholangiocarcinoma: future perspectives for novel therapies. Int J Biol Sci (2022) 18(14):5369–90. doi: 10.7150/ijbs.73949

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Ghoshdastider U, Rohatgi N, Mojtabavi Naeini M, Baruah P, Revkov E, Guo YA, et al. Pan-cancer analysis of ligand-receptor cross-talk in the tumor microenvironment. Cancer Res (2021) 81(7):1802–12. doi: 10.1158/0008-5472.Can-20-2352

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Wang DR, Wu XL, Sun YL. Therapeutic targets and biomarkers of tumor immunotherapy: response versus non-response. Signal Transduct Target Ther (2022) 7(1):331. doi: 10.1038/s41392-022-01136-2

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Maynard A, McCoach CE, Rotow JK, Harris L, Haderk F, Kerr DL, et al. Therapy-induced evolution of human lung cancer revealed by single-cell RNA sequencing. Cell (2020) 182(5):1232–1251.e1222. doi: 10.1016/j.cell.2020.07.017

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Liu F, Wang P, Sun W, Jiang Y, Gong Q. Identification of ligand-receptor pairs associated with tumour characteristics in clear cell renal cell carcinoma. Front Immunol (2022) 13:874056. doi: 10.3389/fimmu.2022.874056

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Cao ZX, Weng X, Huang JS, Long X. Receptor-ligand pair typing and prognostic risk model for papillary thyroid carcinoma based on single-cell sequencing. Front Immunol (2022) 13:902550. doi: 10.3389/fimmu.2022.902550

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics (2010) 26(12):1572–3. doi: 10.1093/bioinformatics/btq170

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Maeser D, Gruener RF, Huang RS. oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data. Brief Bioinform (2021) 22(6). doi: 10.1093/bib/bbab260

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Chu F, Liang T, Chen CLP, Wang X, Ma X. Compact broad learning system based on fused lasso and smooth lasso. IEEE Trans Cybern (2023). doi: 10.1109/tcyb.2023.3267947

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Iasonos A, Schrag D, Raj GV, Panageas KS. How to build and interpret a nomogram for cancer prognosis. J Clin Oncol (2008) 26(8):1364–70. doi: 10.1200/jco.2007.12.9791

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics (2012) 16(5):284–7. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Barbie DA, Tamayo P, Boehm JS, Kim SY, Moody SE, Dunn IF, et al. Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1. Nature (2009) 462(7269):108–12. doi: 10.1038/nature08460

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Li P, Huang T, Zou Q, Liu D, Wang Y, Tan X, et al. FGFR2 promotes expression of PD-L1 in colorectal cancer via the JAK/STAT3 signaling pathway. J Immunol (2019) 202(10):3065–75. doi: 10.4049/jimmunol.1801199

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Yao M, Ventura PB, Jiang Y, Rodriguez FJ, Wang L, Perry JSA, et al. Astrocytic trans-differentiation completes a multicellular paracrine feedback loop required for medulloblastoma tumor growth. Cell (2020) 180(3):502–520.e519. doi: 10.1016/j.cell.2019.12.024

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Sas Z, Cendrowicz E, Weinhäuser I, Rygiel TP. Tumor microenvironment of hepatocellular carcinoma: challenges and opportunities for new treatment options. Int J Mol Sci (2022) 23(7). doi: 10.3390/ijms23073778

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Chow A, Perica K, Klebanoff CA, Wolchok JD. Clinical implications of T cell exhaustion for cancer immunotherapy. Nat Rev Clin Oncol (2022). doi: 10.1038/s41571-022-00689-z

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Ramadori P, Kam S, Heikenwalder M. T cells: Friends and foes in NASH pathogenesis and hepatocarcinogenesis. Hepatology (2022) 75(4):1038–49. doi: 10.1002/hep.32336

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Cheng K, Cai N, Zhu J, Yang X, Liang H, Zhang W. Tumor-associated macrophages in liver cancer: From mechanisms to therapy. Cancer Commun (Lond) (2022) 42(11):1112–40. doi: 10.1002/cac2.12345

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Xu W, Cheng Y, Guo Y, Yao W, Qian H. Targeting tumor associated macrophages in hepatocellular carcinoma. Biochem Pharmacol (2022) 199:114990. doi: 10.1016/j.bcp.2022.114990

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Zhang J, Gu C, Song Q, Zhu M, Xu Y, Xiao M, et al. Identifying cancer-associated fibroblasts as emerging targets for hepatocellular carcinoma. Cell Biosci (2020) 10(1):127. doi: 10.1186/s13578-020-00488-y

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Peng H, Zhu E, Zhang Y. Advances of cancer-associated fibroblasts in liver cancer. biomark Res (2022) 10(1):59. doi: 10.1186/s40364-022-00406-z

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Wang H, Xu Y, Zuo F, Liu J, Yang J. Immune-based combination therapy for esophageal cancer. Front Immunol (2022) 13:1020290. doi: 10.3389/fimmu.2022.1020290

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Liu X, Yang L, Tan X. PD-1/PD-L1 pathway: A double-edged sword in periodontitis. BioMed Pharmacother (2023) 159:114215. doi: 10.1016/j.biopha.2023.114215

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Yuan Y, Adam A, Zhao C, Chen H. Recent advancements in the mechanisms underlying resistance to PD-1/PD-L1 blockade immunotherapy. Cancers (Basel) (2021) 13(4). doi: 10.3390/cancers13040663

CrossRef Full Text | Google Scholar

29. Ma S, Caligiuri MA, Yu J. Harnessing IL-15 signaling to potentiate NK cell-mediated cancer immunotherapy. Trends Immunol (2022) 43(10):833–47. doi: 10.1016/j.it.2022.08.004

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Liu X, Qin X, Qin H, Jia C, Yuan Y, Sun T, et al. Characterization of the heterogeneity of endothelial cells in bleomycin-induced lung fibrosis using single-cell RNA sequencing. Angiogenesis (2021) 24(4):809–21. doi: 10.1007/s10456-021-09795-5

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Chao MP, Weissman IL, Majeti R. The CD47-SIRPα pathway in cancer immune evasion and potential therapeutic implications. Curr Opin Immunol (2012) 24(2):225–32. doi: 10.1016/j.coi.2012.01.010

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Quaranta V, Rainer C, Nielsen SR, Raymant ML, Ahmed MS, Engle DD, et al. Macrophage-derived granulin drives resistance to immune checkpoint inhibition in metastatic pancreatic cancer. Cancer Res (2018) 78(15):4253–69. doi: 10.1158/0008-5472.Can-17-3876

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Czako B, Marszalek JR, Burke JP, Mandal P, Leonard PG, Cross JB, et al. Discovery of IACS-9439, a potent, exquisitely selective, and orally bioavailable inhibitor of CSF1R. J Med Chem (2020) 63(17):9888–911. doi: 10.1021/acs.jmedchem.0c00936

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Krishnan V, Tallapragada S, Schaar B, Kamat K, Chanana AM, Zhang Y, et al. Omental macrophages secrete chemokine ligands that promote ovarian cancer colonization of the omentum via CCR1. Commun Biol (2020) 3(1):524. doi: 10.1038/s42003-020-01246-z

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Ma L, Heinrich S, Wang L, Keggenhoff FL, Khatib S, Forgues M, et al. Multiregional single-cell dissection of tumor and immune cells reveals stable lock-and-key features in liver cancer. Nat Commun (2022) 13(1):7533. doi: 10.1038/s41467-022-35291-5

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Wu C, Thalhamer T, Franca RF, Xiao S, Wang C, Hotta C, et al. Galectin-9-CD44 interaction enhances stability and function of adaptive regulatory T cells. Immunity (2014) 41(2):270–82. doi: 10.1016/j.immuni.2014.06.011

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Lu P, Shih C, Qi H. Ephrin B1-mediated repulsion and signaling control germinal center T cell territoriality and function. Science (2017) 356(6339). doi: 10.1126/science.aai9264

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Li Z, Deng M, Huang F, Jin C, Sun S, Chen H, et al. LILRB4 ITIMs mediate the T cell suppression and infiltration of acute myeloid leukemia cells. Cell Mol Immunol (2020) 17(3):272–82. doi: 10.1038/s41423-019-0321-2

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Ma P, Zhang P, Chen S, Shi W, Ye J, Chen S, et al. Immune cell landscape of patients with diabetic macular edema by single-cell RNA analysis. Front Pharmacol (2021) 12:754933. doi: 10.3389/fphar.2021.754933

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Taftaf R, Liu X, Singh S, Jia Y, Dashzeveg NK, Hoffmann AD, et al. ICAM1 initiates CTC cluster formation and trans-endothelial migration in lung metastasis of breast cancer. Nat Commun (2021) 12(1):4867. doi: 10.1038/s41467-021-25189-z

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Wilson K, Shiuan E, Brantley-Sieders DM. Oncogenic functions and therapeutic targeting of EphA2 in cancer. Oncogene (2021) 40(14):2483–95. doi: 10.1038/s41388-021-01714-8

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Yoo HC, Park SJ, Nam M, Kang J, Kim K, Yeo JH, et al. A variant of SLC1A5 is a mitochondrial glutamine transporter for metabolic reprogramming in cancer cells. Cell Metab (2020) 31(2):267–283.e212. doi: 10.1016/j.cmet.2019.11.020

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Ma H, Qu S, Zhai Y, Yang X. circ_0025033 promotes ovarian cancer development via regulating the hsa_miR-370-3p/SLC1A5 axis. Cell Mol Biol Lett (2022) 27(1):94. doi: 10.1186/s11658-022-00364-2

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Hou L, Hou Y, Liang Y, Chen B, Zhang X, Wang Y, et al. Anti-tumor effects of P-LPK-CPT, a peptide-camptothecin conjugate, in colorectal cancer. Commun Biol (2022) 5(1):1248. doi: 10.1038/s42003-022-04191-1

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Liu X, Qin H, Li Z, Lv Y, Feng S, Zhuang W, et al. Inspiratory hyperoxia suppresses lung cancer metastasis through a MYC/SLC1A5-dependent metabolic pathway. Eur Respir J (2022) 60(6). doi: 10.1183/13993003.00062-2022

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Thorsson V, Gibbs DL, Brown SD, Wolf D, Bortone DS, Ou Yang TH, et al. The immune landscape of cancer. Immunity (2018) 48(4):812–830.e814. doi: 10.1016/j.immuni.2018.03.023

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Zhang J, Zou S, Fang L. Metabolic reprogramming in colorectal cancer: regulatory networks and therapy. Cell Biosci (2023) 13(1):25. doi: 10.1186/s13578-023-00977-w

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Evans JP, Liu SL. Multifaceted roles of TIM-family proteins in virus-host interactions. Trends Microbiol (2020) 28(3):224–35. doi: 10.1016/j.tim.2019.10.004

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Mathewson ND, Ashenberg O, Tirosh I, Gritsch S, Perez EM, Marx S, et al. Inhibitory CD161 receptor identified in glioma-infiltrating T cells by single-cell analysis. Cell (2021) 184(5):1281–1298.e1226. doi: 10.1016/j.cell.2021.01.022

PubMed Abstract | CrossRef Full Text | Google Scholar