ncRNAs-mediated high expression of SEMA3F correlates with poor prognosis and tumor immune infiltration of hepatocellular carcinoma.

Hepatocellular carcinoma (HCC) is notorious for its poor prognosis. Increasing evidence has demonstrated that semaphorin 3F (SEMA3F) plays key roles in initiation and progression of several types of human cancer. However, the specific role and mechanism of SEMA3F in HCC remains not fully determined. In this study, we first performed pan-cancer analysis for SEMA3F's expression and prognosis using The Cancer Genome Atlas (TCGA) and The Genotype-Tissue Expression (GTEx) data and found that SEMA3F might be a potential oncogene in HCC. Subsequently, noncoding RNAs (ncRNAs) contributing to SEMA3F overexpression were identified by a combination of a series of in silico analyses, including expression analysis, correlation analysis, and survival analysis. Finally, the TMPO-AS1/SNHG16-let-7c-5p axis was identified as the most potential upstream ncRNA-related pathway of SEMA3F in HCC. Moreover, SEMA3F level was significantly positively associated with tumor immune cell infiltration, biomarkers of immune cells, and immune checkpoint expression. Collectively, our findings elucidated that ncRNAs-mediated upregulation of SEMA3F correlated with poor prognosis and tumor immune infiltration in HCC.

Introduction

Hepatocellular carcinoma (HCC) is the most common type of primary liver cancer and also ranks as the third leading cause of cancer-related deaths all over the world., Lots of risk factors linked to initiation and progression of HCC have been reported, such as virus infection, cirrhosis, alcohol abuse, and immune system. Despite the huge improvements that have been achieved in aspects of diagnosis, therapy, and prognosis, the outcome of patients with HCC remains unsatisfactory, with more than 700,000 deaths every year. Therefore, it is an urgent need to develop effective therapeutic targets or seek promising prognostic biomarkers in HCC.

In vertebrates, semaphores (SEMAs), first identified as chemo-repulsive molecules for axonal growth cones and followingly found to be implicated in modulating cell motility in the context of vascular growth and tumor metastasis, are a class of proteins that can be generally divided into two subclasses, including transmembrane proteins (classes 1, 4, 6, and 7) and secretory proteins (classes 2 and 3)., It has been well documented that SEMA3F is closely associated with initiation and progression of multiple types of human cancer, including colorectal cancer, breast cancer,, endometrial cancer,, oral squamous cell carcinoma, and head and neck squamous carcinoma. Also, in HCC, SEMA3F could facilitate cancer cell metastasis by activating focal adhesion pathway. SEMAs have also been reported to be involved in immune signaling and immune synapse formation. For example, Casazza et al. suggested that SEMA3A regulated localization and retention of tumor-associated macrophages (TAMs) by functioning as an attractant for TAMs in the aeras of hypoxic tumors. A recent report published by the team of Tracie Plant has also validated that SEMA3F signaling could actively retain neutrophils at sites of inflammation. However, a comprehensive study regarding the expression, prognosis, and mechanism of SMEA3F in HCC is still absent. Moreover, the association of SEMA3F with tumor immune infiltration in HCC is still not determined.

In this study, we first performed expression analysis and survival analysis for SEMA3F in multiple types of human cancer. Next, the noncoding RNA (ncRNA)-associated regulation of SEMA3F, involving microRNAs (miRNAs) and long noncoding RNAs (lncRNAs), was also explored in HCC. Finally, we determined the relationship of SMEA3F expression with immune cell infiltration, biomarkers of immune cells, or immune checkpoints in HCC. Taken together, our findings suggest that ncRNAs-mediated upregulation of SEMA3F correlates with poor prognosis and tumor immune infiltration of patients in HCC.

Results

Pan-cancer analysis of SEMA3F expression

To explore possible roles of SEMA3F in carcinogenesis, we first analyzed its expression in 18 types of human cancer. As shown in Figure 1A, compared with normal samples, SEMA3F was significantly upregulated in 12 cancer types, including BLCA, BRCA, CHOL, COAD, GBM, KICH, KIRC, HCC, LUSC, READ, THCA, and UCEC, and was markedly downregulated in 3 cancer types, involving KIRP, LUAD, and PRAD. However, no significant difference of SEMA3F in ESCA, HNSC, or STAD was observed. Next, we further validated the expression of SEMA3F in these 18 cancer types using GEPIA database. As presented in Figures 1B–1I, SEMA3F expression in BRCA, CHOL, COAD, GBM, KICH, KIRC, HCC, or READ was statistically increased when compared with corresponding normal controls. And in KIRP, LUAD, PRAD, or ESCA, SEMA3F was obviously decreased (Figures 1J–1M). Taken together, SEMA3F was upregulated in BRCA, CHOL, COAD, GBM, KICH, KIRC, HCC, and READ, and downregulated in KIRP, LUAD, and PRAD, indicating that SEMA3F may function as crucial regulator in carcinogenesis of the 11 types of cancer.

Figure 1

Expression analysis for SEMA3F in multiple cancers

(A) The expression of SEMA3F in 18 types of human cancer based on TCGA cancer and normal data. (B–M) SEMA3F expression in TCGA BRCA (B), CHOL (C), COAD (D), GBM (E), KICH (F), KIRC (G), LIHC (H), READ (I), KIRP (J), LUAD (K), PRAD (L), and ESCA (M) tissues compared with corresponding TCGA and GTEx normal tissues. ∗p value < 0.05; ∗∗p value < 0.01; ∗∗∗p value < 0.001.

The prognostic values of SEMA3F in human cancer

Next, survival analysis for SEMA3F in BRCA, CHOL, COAD, GBM, KICH, KIRC, HCC, READ, KIRP, LUAD, and PRAD was conducted. Two prognostic indices, consisting of overall survival (OS) and disease-free survival (RFS), were included. For OS, high expression of SEMA3F in GBM and HCC had unfavorable prognosis but KIRC patients with higher expression of SEMA3F indicated better prognosis (Figure 2). For RFS, among all cancer types, only increased expression of SEMA3F indicated poor prognosis in HCC (Figure 3). No statistical significance of SEMA3F for predicting prognosis of patients in other cancer types was observed. By combination of OS and RFS, SEMA3F may be utilized as an unfavorable prognostic biomarker in patients with HCC.

Figure 2

The overall survival (OS) analysis for SEMA3F in various human cancers determined by GEPIA database

(A–K) The OS plot of SEMA3F in BRCA (A), CHOL (B), COAD (C), GBM (D), KICH (E), KIRC (F), LIHC (G), READ (H), KIRP (I), LUAD (J), and PRAD (K).

Figure 3

The disease-free survival (RFS) analysis for SEMA3F in various human cancers determined by GEPIA database

(A–K) The RFS plot of SEMA3F in BRCA (A), CHOL (B), COAD (C), GBM (D), KICH (E), KIRC (F), LIHC (G), READ (H), KIRP (I), LUAD (J), and PRAD (K).

Prediction and analysis of upstream miRNAs of SEMA3F

It has been widely acknowledged that ncRNAs are responsible for the regulation of gene expression. To ascertain whether SEMA3F was modulated by some ncRNAs, we first predicted upstream miRNAs that could potentially bind to SEMA3F and finally found 13 miRNAs. To improve visualization, a miRNA-SEMA3F regulatory network was established using cytoscape software (Figure 4A). Based on the action mechanism of miRNA in regulation of target gene expression, there should be negative correlation between miRNA and SEMA3F. Thus, the expression correlation analysis was performed. As listed in Figure 4B, SEMA3F was significantly negatively correlated with let-7c-5p and positively correlated with let-7e-5p in HCC. There were no statistical expression relationships between SEMA3F and the other 11 predicted miRNAs. Finally, the expression and prognostic value of let-7c-5p in HCC were determined. As presented in Figures 4C and 4D, let-7c-5p was markedly downregulated in HCC and its upregulation was positively linked to patients’ prognosis. All these findings suggest that let-7c-5p might be the most potential regulatory miRNA of SEMA3F in HCC.

Figure 4

Identification of let-7c-5p as a potential upstream miRNA of SEMA3F in HCC

(A) The miRNA-SEMA3F regulatory network established by cytoscape software. (B) The expression correlation between predicted miRNAs and SEMA3F in HCC analyzed by starBase database. (C) The expression of let-7c-5p in HCC and control normal samples determined by starBase database. (D) The prognostic value of let-7c-5p in HCC assessed by Kaplan-Meier plotter.

Prediction and analysis of upstream lncRNAs of let-7c-5p

Next, the upstream lncRNAs of let-7c-5p were predicted using starBase database. A total of 53 possible lncRNAs were forecasted. Identically, to improve visualization, a lncRNA-let-7c-5p regulatory network was constructed by cytoscape software (Figure S1). Then, the expression levels of these lncRNAs in HCC were determined using GEPIA. As shown in Figures 5A–5D, among all the 53 lncRNAs, only HEIH, TMPO-AS1, SNHG16, and LINC00665 were significantly upregulated in HCC compared with normal controls. Subsequently, the prognostic values of the four lncRNAs in HCC were assessed. As suggested in Figures 5E–5L, only HCC patients with higher expression of SNHG16 possessed both poorer OS and RFS. Besides, overexpressed TMPO-AS1 indicated poor RFS of patients with HCC. According to the competing endogenous RNA (ceRNA) hypothesis, lncRNA could increase mRNA expression by competitively binding to shared miRNAs. Therefore, there should be negative correlation between lncRNA and miRNA or positive correlation between lncRNA and mRNA. As listed in Table 1, the expression correlation between the four lncRNAs and let-7c-5p or SEMA3F in HCC was also detected using starBase database. Taking expression analysis, survival analysis, and correlation analysis into consideration, SNHG16 and TMPO-AS1 might be the two most potential upstream lncRNAs of let-7c-5p/SEMA3F axis in HCC.

Figure 5

Expression analysis and survival analysis for upstream lncRNAs of let-7c-5p in HCC

(A–D) The expression of HEIH (A), TMPO-AS1 (B), SNHG16 (C), and LINC00665 (D) in TCGA HCC compared with “TCGA normal” or “TCGA and GTEx normal” data. (E–L) The OS analysis for HEIH (E), TMPO-AS1 (F), SNHG16 (G), and LINC00665 (H) in HCC. The RFS for HEIH (I), TMPO-AS1 (J), SNHG16 (K), and LINC00665 (L) in HCC. ∗p value < 0.05.

Table 1

Correlation analysis between lncRNA and let-7c-5p or lncRNA and SEMA3F in HCC determined by starBase database.

lncRNA	miRNA	R value	p value
HEIH	let-7c-5p	–0.156a	2.65E-03∗∗a
TMPO-AS1	let-7c-5p	–0.167a	1.25E-03∗∗a
SNHG16	let-7c-5p	–0.169a	1.10E-03∗∗a
LINC00665	let-7c-5p	–0.250a	1.07E-06∗∗∗a
lncRNA	mRNA	R value	p value
HEIH	SEMA3F	0.041	4.25E-01
TMPO-AS1	SEMA3F	0.255a	5.62E-07∗∗∗a
SNHG16	SEMA3F	0.162a	1.72E-03∗∗a
LINC00665	SEMA3F	0.185a	3.29E-04∗∗∗a

aThese results are statistically significant.

p value < 0.01; ∗∗∗p value < 0.001.

SEMA3F positively correlates with immune cell infiltration in HCC

SEMA3F is a member of SEMAs, which have an immunoglobulin domain and are reported to play a critical role in the immune system. As shown in Figure 6A, no significant change of immune cell infiltration level under various copy numbers of SEMA3F in HCC was observed. Correlation analysis could provide key clues for studying the function and mechanism of SEMA3F. Thus, the correlation of SEMA3F expression level with immune cell infiltration level was evaluated. As presented in Figures 6B–6G, SEMA3F expression was significantly positively associated with all analyzed immune cells, including B cell, CD8+ T cell, CD4+ T cell, macrophage, neutrophil, and dendritic cell in HCC.

Figure 6

The relationship of immune cell infiltration with SEMA3F level in HCC

(A) The infiltration level of various immune cells under different copy numbers of SEMA3F in HCC. (B–G) The correlation of SEMA3F expression level with B cell (B), CD8+ T cell (C), CD4+ T cell (D), macrophage (E), neutrophil (F), or dendritic cell (G) infiltration level in HCC.

Expression correlation of SEMA3F and biomarkers of immune cells in HCC

To further explore the role of SEMA3F in tumor immune, we determined the expression correlation of SEMA3F with biomarkers of immune cells in HCC using GEPIA database. As listed in Table 2, SEMA3F was significantly positively correlated with B cell’s biomarkers (CD19 and CD79A), CD8+ T cell’s biomarkers (CD8A and CD8B), CD8+ T cell’s biomarker (CD4), M1 macrophage’s biomarkers (NOS2, IRF5, and PTGS2), M2 macrophage’s biomarkers (CD163, VSIG4, and MS4A4A), neutrophil’s biomarkers (ITGAM and CCR7), and dendritic cell’s biomarkers (HLA-DPB1, HLA-DRA, HLA-DPA1, CD1C, NRP1, and ITGAX) in HCC. These findings partially support that SEMA3F is positively linked to immune cell infiltration.

Table 2

Correlation analysis between SEMA3F and biomarkers of immune cells in HCC determined by GEPIA database.

Immune cell	Biomarker	R value	p value
B cell	CD19	0.17a	8.9E-04∗∗∗a
	CD79A	0.19a	1.8E-04∗∗∗a
CD8+ T cell	CD8A	0.18a	4.5E-04∗∗∗a
	CD8B	0.13a	9.6E-03∗∗a
CD4+ T cell	CD4	0.29a	2.0E-08∗∗∗a
M1 macrophage	NOS2	0.37a	4.2E-13∗∗∗a
	IRF5	0.30a	5.4E-09∗∗∗a
	PTGS2	0.33a	4.8E-11∗∗∗a
M2 macrophage	CD163	0.11a	3.0E-02∗a
	VSIG4	0.16a	2.2E-03∗∗a
	MS4A4A	0.24a	4.1E-06∗∗∗a
Neutrophil	CEACAM8	−0.03	6.1E-01
	ITGAM	0.27a	1.2E-07∗∗∗a
	CCR7	0.23a	1.1E-05∗∗∗a
Dendritic cell	HLA-DPB1	0.25a	1.3E-06∗∗∗a
	HLA-DQB1	0.04	4.5E-01
	HLA-DRA	0.26a	4.3E-07∗∗∗a
	HLA-DPA1	0.27a	2.4E-07∗∗∗a
	CD1C	0.21a	4.6E-05∗∗∗a
	NRP1	0.48a	7.3E-23∗∗∗a
	ITGAX	0.35a	5.6E-12∗∗∗a

aThese results are statistically significant.

p value < 0.05; ∗∗p value < 0.01; ∗∗∗p value < 0.001.

Relationship between SEMA3F and immune checkpoints in HCC

PD1/PD-L1 and CTLA-4 are important immune checkpoints that are responsible for tumor immune escape. Considering the potential oncogenic role of SEMA3F in HCC, the relationship of SEMA3F with PD1, PD-L1, or CTLA-4 was assessed. As suggested in Figures 7A–7C, SEMA3F expression was significantly positively correlated with PD1, PD-L1, and CTLA-4 in HCC, which was adjusted by purity. Similar to TIMER data analysis, we also found that there was significant positive correlation of SEMA3F with PD1, PD-L1, or CTLA-4 in HCC (Figures 7D–7F). These results demonstrate that tumor immune escape might be involved in SEMA3F-mediated carcinogenesis of HCC.

Figure 7

Correlation of SEMA3F expression with PD-1, PD-L1, and CTLA-4 expression in HCC

(A) Spearman correlation of SEMA3F with expression of PD-1 in HCC adjusted by purity using TIMER. (B) Spearman correlation of SEMA3F with expression of PD-L1 in HCC adjusted by purity using TIMER. (C) Spearman correlation of SEMA3F with expression of CTLA-4 in HCC adjusted by purity using TIMER. (D) The expression correlation of SEMA3F with PD1 in HCC determined by GEPIA database. (E) The expression correlation of SEMA3F with PD-L1 in HCC determined by GEPIA database. (F) The expression correlation of SEMA3F with CTLA-4 in HCC determined by GEPIA database.

Discussion

To date, HCC is still notorious for its poor prognosis. Elucidating the molecular mechanism of HCC carcinogenesis may provide key clues for developing effective therapeutic targets or seeking promising prognostic biomarkers. Increasing evidence has demonstrated that SEMA3F plays key roles in initiation and progression of multiple human cancers, including HCC. However, the knowledge of SEMA3F in HCC remains inadequate and needs to be further investigated.

In this study, we first conducted pan-cancer analysis of SEMA3F’s expression using The Cancer Genome Atlas (TCGA) data, after which GEPIA database was further employed to validate SEMA3F’s expression. Survival analysis for SEMA3F in those cancer types of interest indicated that HCC patients with high expression of SEMA3F had poor prognosis. Ye et al. suggested that SEMA3F could activate focal adhesion pathway, thus leading to metastasis of HCC. This report together with our analytic results showed the oncogenic role of SEMA3F in HCC.

It has been well documented that ncRNAs, including miRNAs, lncRNAs, and circular RNAs (circRNAs), participated in regulation of gene expression by talking with each other through the ceRNA mechanism.18, 19, 20, 21, 22 To explore the upstream regulatory miRNAs of SEMA3F, we introduced seven prediction programs, involving PITA, RNA22, miRmap, microT, miRanda, PicTar, and TargetScan, to predict possible miRNAs that could potentially bind to SEMA3F. At the end, 13 miRNAs were finally obtained. Most of these miRNAs have been found to act as tumor-suppressive miRNAs in HCC. For example, let-7a-5p enhanced the sensitivity of HCC cells to cetuximab and inhibited self-renewal of HCC stem-like cells,, miR-98-5p suppressed HCC cell proliferation by targeting EZH2 and IGF2BP1,, and miR-328-3p inhibited malignant progression of HCC by decreasing MMP9 expression level. After performing correlation analysis, expression analysis, and survival analysis, let-7c-5p was selected as the most potential upstream tumor suppressive miRNA of SEMA3F. Previous studies also showed that let-7c-5p played inhibitory roles in modulating proliferation and migration of HCC.

Based on the ceRNA hypothesis, the potential lncRNAs of let-7c-5p/SEMA3F axis should be oncogenic lncRNAs in HCC. Next, upstream lncRNAs of let-7c-5p/SEMA3F axis were also predicted and 53 possible lncRNAs were found. By conducting expression analysis, survival analysis, and correlation analysis, two of the most potential upregulated lncRNAs, including TMPO-AS1 and SNHG16, were identified. The two lncRNAs have been reported to function as oncogenes in multiple malignancies, including HCC. For instance, TMPO-AS1 enhanced development of HCC, and SNHG16 promoted proliferation, migration, invasion, and sorafenib resistance of HCC.32, 33, 34 Taken together, TMPO-AS1 and SNHG16/let-7c-5p/SEMA3F axis were identified as potential regulatory pathways in HCC.

Numerous studies have confirmed that tumor immune cell infiltration could influence the efficacies of chemotherapy, radiotherapy, or immunotherapy and prognosis of cancer patients.35, 36, 37 Our work suggested that SEMA3F was significantly positively correlated with various immune cells, including B cell, CD8+ T cell, CD4+ T cell, macrophage, neutrophil, and dendritic cell in HCC. Moreover, SEMA3F was also markedly positively associated with biomarkers of these infiltrated immune cells. These findings indicated that tumor immune infiltration might partially account for SEMA3F-mediated oncogenic roles in HCC.

Also, the efficacy of immunotherapy not only needs adequate immune cells infiltrating in tumor microenvironment but also depends on the sufficient expression of immune checkpoints. Thus, we also assessed the relationship between SEMA3F and immune checkpoints. The results demonstrated that high expression of SEMA3F was strongly linked to PD1, PD-L1, or CTLA-4 in HCC, indicating that targeting SEMA3F might increase the efficacy of immunotherapy in HCC.

In summary, we elucidated that SEMA3F was highly expressed in multiple types of human cancer (including HCC) and positively correlated with unfavorable prognosis in HCC. We identified an upstream regulatory mechanism of SEMA3F in HCC, namely TMPO-AS1/SNHG16-let-7c-5p axis (Figure 8). Furthermore, our current findings also indicated that SEMA3F might exert its oncogenic roles through increasing tumor immune cell infiltration and immune checkpoint expression. However, these results should be validated by much more basic experiments and large clinical trials in the future.

Figure 8

The model of TMPO-AS1/SNHG16-let-7c-5p-SEMA3F axis in carcinogenesis of HCC

Materials and methods

TCGA data download, process, and analysis

The mRNA expression data of 18 cancer types (BLCA, BRCA, CHOL, COAD, ESCA, GBM, HNSC, KICH, KIRC, KIRP, LIHC, LUAD, LUSC, PRAD, READ, STAD, THCA, and UCEC) were downloaded from TCGA database (https://genome-cancer.ucsc.edu/), after which these data were normalized and then differential expression analysis was performed for SEMA3F using R package limma. p value <0.05 was considered as statistically significant.

GEPIA database analysis

GEPIA (http://gepia.cancer-pku.cn/) is a web tool for cancer and normal gene-expression profiling and interactive analyses based on TCGA and The Genotype-Tissue Expression (GTEx) data. GEPIA was used to determine SEMA3F and lncRNA expression in various types of human cancer. p value <0.05 was considered as statistically significant. GEPIA was employed to conduct survival analysis for SEMA3F in 11 various cancer types, including OS and RFS. GEPIA was also utilized to assess the prognostic values of candidate lncRNAs in HCC. Log rank p value <0.05 was considered as statistically significant. In addition, expression correlation of SEMA3F with immune checkpoints in HCC was also evaluated using GEPIA database. |R| >0.1 and p value <0.05 were set as selection criteria for identifying as statistically significant.

Candidate miRNA prediction

Upstream binding miRNAs of SEMA3F were predicted by several target gene prediction programs, consisting of PITA, RNA22, miRmap, microT, miRanda, PicTar, and TargetScan. Only the predicted miRNAs that commonly appeared in more than two programs as mentioned above were included for subsequent analyses. These predicted miRNAs were regarded as candidate miRNAs of SEMA3F.

starBase database analysis

starBase (http://starbase.sysu.edu.cn/) is a database for exploring miRNA-related studies. starBase was introduced to perform expression correlation analysis for miRNA-SEMA3F, lncRNA-let-7c-5p, or lncRNA-SEMA3F in HCC. The expression level of let-7c-5p in HCC and normal controls was also analyzed by starBase. Besides, starBase was used to predict candidate lncRNAs that could potentially bind to let-7c-5p.

Kaplan-Meier plotter analysis

Kaplan-Meier plotter (http://kmplot.com/analysis/), an online database capable of accessing the effects of genes or miRNAs on survival in more than 20 cancer types including HCC, was employed to conduct survival analysis for let-7c-5p in HCC as we previously described. Log rank p value <0.05 was considered as statistically significant.

TIMER database analysis

TIMER (https://cistrome.shinyapps.io/timer/) is a web server for comprehensive analysis of tumor-infiltrating immune cells. TIMER was used to analyze the correlation of SEMA3F expression level with immune cell infiltration level or immune checkpoint expression level in HCC. p value <0.05 was considered as statistically significant.

Statistical analysis

The statistical analysis in this study was automatically calculated by the online database mentioned above. p value <0.05 or log rank p value <0.05 was considered as statistically significant.