Construction of a competing endogenous RNA network in head and neck squamous cell carcinoma by pan-cancer analysis.

Background: Head and neck squamous cell carcinoma (HNSC) is the sixth most common cancer worldwide, and new cases are anticipated to reach 1.08 million in 2030. Our study aimed to identify the competing endogenous RNAs (ceRNAs) involved in HNSC tumorigenesis.
Methods: First, a pan-cancer correlation analysis was conducted on the expression and survival conditions of sideroflexin (SFXN3) based on data downloaded from the Xena database. Second, the upstream regulatory microRNAs (miRNAs) and long non-coding RNAs (lncRNAs) of SFXN3 were predicted using the Encyclopedia of RNA Interactomes (ENCORI) database. Expression and survival analyses were subsequently used to construct lncRNA-miRNA-mRNA ceRNA network that correlated with HNSC. Third, the proportion of various types of immune cells in HNSC was calculated using the CIBERSORT algorithm. Finally, a correlation analysis was performed on SFXN3, including immune cell infiltration (ICI), clinical stage, and immune checkpoints.
Results: The pan-cancer analysis suggested that SFXN3 was up-regulated in HNSC, and it correlated with poor prognosis. The ceRNA regulatory network MIR193BHG-miR-29c-3p-SFXN3 was identified as one of the potential biological regulatory pathways of HNSC. The upstream lncRNA MIR193BHG was associated with a poor prognosis in HNSC, and its target gene SFXN3 was correlated with tumor ICI, immune cell biomarkers, and immune checkpoints. Conclusions: By performing ceRNA analysis, our study demonstrated that MIR193HG-miR-29c-3p-SFXN3 is significantly involved in HNSC, and this action axis markedly affect the therapeutic effect and prognosis. 2022 Translational Cancer Research. All rights reserved.

Introduction

Head and neck squamous cell carcinoma (HNSC) is a heterogeneous group of neoplasms, comprising tongue carcinoma (TSCC), oral carcinoma (OSCC), laryngeal carcinoma (LSCC), and nasopharyngeal carcinoma (NPC) forms of the disease, mainly caused by exposure to certain pathogenic factors, such as smoking. HNSC is the sixth most common cancer worldwide, with 890,000 new cases and 450,000 deaths in 2018, and its incidence is anticipated to reach 1.08 million new cases annually in 2030 (1).

Long non-coding RNAs (lncRNAs) can function as competing endogenous RNAs (ceRNAs) to protect messenger RNAs (mRNAs) by acting as molecular sponges for microRNAs (miRNAs) that specifically inhibit the target mRNAs. The ceRNA hypothesis posits that mRNAs, lncRNAs, and circular RNAs (circRNAs) competitively connect with miRNAs, forming a ceRNA network and regulating gene expression. A previous study has reported that ceRNAs were related to the occurrence and development of cancer and could function as a diagnostic marker or therapeutic target (2). For example, infarction-associated transcript lncRNA (MIAT) acts as a ceRNA and promotes papillary thyroid carcinoma (PTC) cell invasion (3). A previous study compared the gene expression changes between NPC and noncancerous human immortalized nasopharyngeal epithelial cell lines, Epstein-Barr virus+ (EBV+) compared with Epstein-Barr virus− (EBV−) NPC cell lines, and metastatic vs. nonmetastatic patient samples, and constructed a circRNA-miRNA-mRNA ceRNA network. By ceRNA network analysis, they successfully identified the non-coding RNA axis that was closely biologically linked with NPC tumorigenesis (4). However, the role of ceRNA network in HNSC tumorigenesis is still unclear. It was reported by clinical trials that tumor immune cell infiltration (ICI) is correlated with immunotherapy and prognosis (5). The Cancer Genome Atlas (TCGA) is one of the most practical and successful cancer genomics programs. Biomolecular data of 33 types of cancers from more than 11,000 individuals are available in TCGA, including genomic sequence, expression, and survival condition, which can be used to analyze the deeper relationships between genes and diseases (6).

Sideroflexin (SFXN3) is a member of the Sfxn homologs of sideroflexin 1 (Sfxn1) and a coding protein involved in oral squamous cell carcinoma. Tissue-specific expression of SFXN3 causes the occurrence and development of retinal degeneration and squamous cell carcinoma (7-9). However, no investigation was conducted toward examining whether the effect of SFXN3 is associated with HNSC. The goal of our study was to identify the ceRNA network involved in HNSC tumorigenesis and find out the potential therapeutic target. According to the ceRNA hypothesis, we predicted a lncRNA–miRNA–SFXN3 axis, which is the potential pathogenic mechanism with poor prognosis in HNSC.

In this study, we performed a pan-cancer analysis of data on 39 kinds of tumors using TIMER 2.0. Through a co-expression analysis and a survival analysis, we found that SFXN3 exhibited tissue-specific up-regulated expression only in HNSC, along with poor outcomes. Therefore, we speculated that SFXN3 might be closely linked to HNSC tumorigenesis. The upstream miRNAs and lncRNAs of SFXN3 were predicted, then used to construct a ceRNA action network of lncRNA-miRNA-SFXN3. The HNSC ceRNA network consisted of a MIR193BHG–hsa-miR-29c-3p–SFXN3 axis. The relationship between SFXN3 and ICI, immune checkpoints, immune cell markers, and clinical stage, were assessed. SFXN3 proved to be an important biological factor related to poor prognosis and tumor immune infiltration of HNSC. These findings demonstrated that the MIR193BHG–hsa-miR-29c-3p–SFXN3 axis is a key therapeutic target of HNSC. We present the following article in accordance with the REMARK reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-22-632/rc).

Methods

Data collection, preprocessing, and differential expression analysis

The University of California Santa Cruz (UCSC) Xena (https://xena.ucsc.edu/) is a tool widely used for visual integration and exploration for multi-omics data and associated clinical and phenotypic annotations (10). The raw data were downloaded from UCSC Xena, including the gene expression RNA sequencing data and the survival data of 33 kinds of cancers, which were originally from Genomic Data Commons (GDC) in TCGA, and then organized by the perl script. The data from all the samples in the database were downloaded and used. As for the gene expression RNA sequencing data, the same samples with different vials/portions/analytes/aliquotes were averaged, while different samples were combined into genomic matrix. Regarding the survival data, clinical phenotype data related to prognosis in TCGA were collected, including overall survival and disease-free survival. The differential expression analysis was conducted using R (R Foundation for Statistical Computing, Vienna, Austria) with packages (BiocManager and limma). We used the Wilcoxon signed-rank test to identify differential expression genes (DEGs) with an adjusted P value of less than 0.05 and |log2FC| of greater than 1. The survival analysis was performed with the survival package of R using DEGs and the survival data of HNSC, and P<0.05 was the statistical significance criterion.

TIMER platform analysis

TIMER 2.0 (http://timer.cistrome.org/) is a biological tool for evaluating associations between immune infiltrates and genetic or clinical regulations, which can provide comprehensive analysis and visualization functions of tumor-infiltrating immune cells (11). It was used to perform a pan-cancer expression analysis of SFXN3 in the total 39 kinds of cancers. The results of differential expression analysis were output as a boxplot following the Wilcoxon signed-rank test, which was annotated by the number of stars (*, P value <0.05; **, P value <0.01; ***, P value <0.001).

In addition, the immune checkpoint genes, including programmed cell death 1 (PDCD1), programmed cell death ligand 1 (CD274), and cytotoxic T lymphocyte antigen 4 (CTLA4), are widely used to check on the immune situation, especially in individuals with cancer (12,13). We searched for the association between SFXN3 and the three immune checkpoints with purity adjustment, where P<0.05 was considered to be statistically significant.

Construction of the ceRNA network

The upstream miRNAs and lncRNAs were identified by conducting a differential expression analysis, co-expression analysis, and survival analysis. P<0.001 was considered statistically significant. Then, the miRNAs and lncRNAs as well as SFXN3 would be used to build up a network using Cytoscape.

Gene Expression Profiling Interactive Analysis (GEPIA) database analysis

GEPIA (http://gepia.cancer-pku.cn/) is an online tool based on TCGA, and the Genotype-Tissue Expression (GTEx) data that offer fast and customizable functionalities, including differential expression analysis, patient survival analysis, and analysis of immune checkpoints and targeted genes (14). |R| >0.1 and P<0.05 were set to identify the selection criteria of statistical significance.

Statistical analysis

In this study, statistical analysis was performed with R packages, Perl, or online databases, depending on the analysis needs. |log2FC| >1 and P<0.05 were considered to be statistically significant.

Ethical statement

The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Results

Pan-cancer expression analysis of SFXN3

The expression condition of SFXN3 in 39 types of cancers with 11,057 samples is shown in . Only 15 types of cancers in the total 39 kinds of cancers were considered statistically significant (P<0.001) using the Wilcoxon signed-rank test (). SFXN3 was significantly differently expressed in 14 cancers, with upregulation in HNSC, cholangiocarcinoma (CHOL), colon adenocarcinoma (COAD), kidney renal clear cell carcinoma (KIRC), kidney renal papillary cell carcinoma (KIRP), liver hepatocellular carcinoma (LIHC), rectum adenocarcinoma (READ), stomach adenocarcinoma (STAD), and thyroid carcinoma (THCA), while it was down-regulated in glioblastoma multiforme (GBM), kidney chromophobe (KICH), lung squamous cell carcinoma (LUSC), prostate adenocarcinoma (PRAD), and uterine corpus endometrial carcinoma (UCEC). On the contrary, SFXN3 was not significantly differently expressed in other cancers, including adrenocortical carcinoma (ACC), bladder urothelial carcinoma (BLCA), breast invasive carcinoma (BRCA), cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC), lymphoid neoplasm diffuse large B-cell lymphoma (DLBC), esophageal carcinoma (ESCA), acute myeloid leukemia (LAML), lung adenocarcinoma (LUAD), mesothelioma (MESO), ovarian serous cystadenocarcinoma (OV), pancreatic adenocarcinoma (PAAD), pheochromocytoma and paraganglioma (PCPG), sarcoma (SARC), skin cutaneous melanoma (SKCM), testicular germ cell tumors (TGCT), thymoma (THYM), uterine carcinosarcoma (UCS), uveal melanoma (UVM). Then, these 14 kinds of cancers were analyzed using the GEPIA database. As shown in , SFXN3 was significantly up-regulated in HNSC, CHOL, and KIRC, but markedly down-regulated in KICH, LUSC, PRAD, and UCEC. These data indicate that SFXN3 might play a key role in seven kinds of cancer tumorigenesis, including HNSC, CHOL, KIRC, KICH, LUSC, PRAD, and UCEC.

Figure 1

Differential expression analysis of the SFXN3 gene in 39 types of cancer. (A) The expression of SFXN3 in 39 types of cancer, including HNSC, CHOL, COAD, KIRC, KIRP, LIHC, READ, STAD, and THCA, while it was down-regulated in GBM, KICH, LUSC, PRAD, UCEC, ACC, BLCA, BRCA, CESC, DLBC, ESCA, LAML, LUAD, MESO, OV, PAAD, PCPG, SARC, SKCM, TGCT, THYM, UCS and UVM. (B) The expression of SFXN3 in up-regulated cancers based on both TCGA and the GTEx dataset. (C) The expression of SFXN3 in down-regulated cancers based on both TCGA and GTEx. *, P<0.05; **, P<0.01; ***, P<0.001. SFXN3, sideroflexin 3; HNSC, head and neck squamous cell carcinoma; CHOL, cholangiocarcinoma; COAD, colon adenocarcinoma; KIRC, kidney renal clear cell carcinoma; KIRP, kidney renal papillary cell carcinoma; LIHC, liver hepatocellular carcinoma; READ, rectum adenocarcinoma; STAD, stomach adenocarcinoma; THCA, thyroid carcinoma; GBM, glioblastoma multiforme; KICH, kidney chromophobe; LUSC, lung squamous cell carcinoma; PRAD, prostate adenocarcinoma; UCEC, uterine corpus endometrial carcinoma; ACC, adrenocortical carcinoma; BLCA, bladder urothelial carcinoma; BRCA, breast invasive carcinoma; CESC, cervical squamous cell carcinoma and endocervical adenocarcinoma; DLBC, lymphoid neoplasm diffuse large B-cell lymphoma; ESCA, esophageal carcinoma; LAML, acute myeloid leukemia; LUAD, lung adenocarcinoma; MESO, mesothelioma; OV. ovarian serous cystadenocarcinoma; PAAD, pancreatic adenocarcinoma; PCPG, pheochromocytoma and paraganglioma; SARC, sarcoma; SKCM, skin cutaneous melanoma; TGCT, testicular germ cell tumors; THYM, thymoma; UCS, uterine carcinosarcoma; UVM, uveal melanoma; TCGA, The Cancer Genome Atlas; GTEx, Genotype-Tissue Expression.

Pan-cancer survival analysis of SFXN3

To further confirm the role of SFXN3 in tumorigenesis, the overall survival analysis and the disease-free survival analysis of seven kinds of tumors were performed with GEPIA (). We performed Kaplan-Meier analysis and log-rank test on these seven kinds of tumors. Importantly, we only found that the expression level of SFXN3 was negatively correlated with HNSC overall survival rate (P=0.056) (). Interestingly, the expression level of SFXN3 was negatively correlated with HNSC (P=0.037) and KICH (P=0.025) in the disease-free survival analysis (). These data indicated that SFXN3 might be a biomarker of poor prognosis in patients with HNSC.

Figure 2

The overall survival analysis and disease-free survival analysis for SFXN3 in seven cancers. (A-G) The overall survival analysis for SFXN3 in head and neck squamous cell carcinoma, kidney chromophobe, lung squamous cell carcinoma, kidney renal clear cell carcinoma, prostate adenocarcinoma, uterine corpus endometrial carcinoma, and cholangiocarcinoma. (H-N) The disease-free survival analysis for SFXN3 in head and neck squamous cell carcinoma, kidney chromophobe, lung squamous cell carcinoma, kidney renal clear cell carcinoma, prostate adenocarcinoma, uterine corpus endometrial carcinoma, and cholangiocarcinoma. SFXN3, sideroflexin 3; HNSC, head and neck squamous cell carcinoma; KICH, kidney chromophobe; LUSC, lung squamous cell carcinoma; KIRC, kidney renal clear cell carcinoma; PRAD, prostate adenocarcinoma; UCEC, uterine corpus endometrial carcinoma; CHOL, cholangiocarcinoma.

Predication analysis of upstream miRNA of SFXN3

The Encyclopedia of RNA Interactomes (ENCORI) (https://starbase.sysu.edu.cn/) is an open database available for biological investigations, such as of non-coding RNA interactions (15). ENCORI contains seven miRNA target prediction software packages, including PITA, RNA22, miRmap, microT, miRanda, PicTar, and TargetScan. In this study, we identified the upstream miRNAs that could be targeted while binding to SFXN3 by ENCORI. To improve accuracy, the predicted miRNAs were selected as the upstream miRNAs of SFXN3 when they were identified by no less than two prediction programs. A total of 32 miRNAs were identified as the upstream miRNAs of SFXN3, and they were subsequently used to construct a miRNA-mRNA action network using Cytoscape (16) ().

Figure 3

The network between SFXN3 and microRNAs. SFXN3, sideroflexin 3.

The miRNA expression data in HNSC were downloaded from GDC (https://portal.gdc.cancer.gov/) (17). A total of 2,242 mature miRNAs from 44 normal samples and 525 tumor samples were downloaded from GDC. Only 32 miRNAs that predicted binding to SFXN3 were ultimately selected for co-expression analysis. As a result, three miRNAs, including hsa-miR-423-5p (P=3.60E-06), hsa-miR-29c-3p (P=5.68E-06), and hsa-miR-361-5p (P=7.05E-06) (), were co-expressed with SFXN3, which was identified using the Wilcoxon signed-rank test. The co-expression condition between SFXN3 and miRNAs is shown in . Moreover, the differential expression analysis of each miRNA between tumor and normal samples in HNSC was performed (). As there should be a negative correlation between miRNA and SFXN3, the results of expression correlation analysis identified that only hsa-miR-29c-3p was markedly negatively correlated with SFXN3 in HNSC (P<0.01, logFC <0) (). Hsa-miR-423-5p, hsa-miR-29c-3p, and hsa-miR-361-5p were further assessed for survival analysis using log-rank test. As shown in , hsa-miR-29c-3p (P=0.023) and hsa-miR-361-5p (P=0.032) were significantly correlated with the overall survival rate. Importantly, the higher expression level of hsa-miR-29c-3p was accompanied by a longer overall survival rate.

Table 1

Co-expression analysis between SFXN3 and predicted miRNAs

Gene	miRNA	Cor	P value	logFC	Diff Pval
SFXN3	hsa-miR-423-5p	−0.206572705	3.60E-06	0.57694	1.81E-07
SFXN3	hsa-miR-29c-3p*	−0.202388935*	5.68E-06*	−2.12259*	1.04E-19*
SFXN3	hsa-miR-361-5p	−0.200370216	7.05E-06	0.04434	0.831042

*, represents the targeted miRNA with logFC <0.

Figure 4

The correlation analysis of miRNAs in head and neck squamous cell carcinoma. (A-C) The co-expression analysis between SFXN3 and hsa-miR-29c-3p, hsa-miR-423-5p, and hsa-miR-361-5p. (D-F) The differential expression analysis of hsa-miR-29c-3p, hsa-miR-423-5p, and hsa-miR-361-5p. (G-I) The overall survival analysis of hsa-miR-29c-3p, hsa-miR-423-5p, and hsa-miR-361-5p. miRNAs, microRNAs.

Therefore, has-miR-29c-3p was identified to be the upstream miRNA of SFXN3, and the correlated pair of has-miR-29c-3p and SFXN3 was involved in HNSC prognosis.

Predication analysis of upstream lncRNA of has-miR-29c-3p

The ENCORI database was employed to predict the upstream lncRNA of has-miR-29c-3p, and a total of 312 potential lncRNAs were identified. The action network between hsa-miR-29c-3p and the 312 lncRNAs was constructed using Cytoscape ().

Figure 5

The network between hsa-miR-29c-3p and long non-coding RNAs.

The co-expression condition between lncRNAs and the miRNA hsa-miR-29c-3p analyzed using the Wilcoxon signed-rank test is shown in . According to the co-expression analysis (), NOP14-AS1 (P=2.9E-10), MIR193BHG (P=7.8E-16), and LINC02323 (P=1.4E-07) were negatively correlated with hsa-miR-29c-3p. Referring to the expression correlation analysis and the Wilcoxon signed-rank test between the gene SFXN3 and lncRNAs, NOP14-AS1 (P=5.1E-09), MIR193BHG (P=1.2E-14), and LINC02323 (P=0.5) showed a positive correlation with SFXN3, but only the P values of MIR193BHG and NOP14-AS1 were less than 0.001 (). In addition, these lncRNAs were found to be up-regulated in tumor samples (P<0.01), which is shown using boxplots (). The survival analysis showed that the high-level expressions of NOP14-AS1 (P=0.247), LINC02323 (P=0.103), and MIR193BHG (P=9.7E-06) tended to bring about a decrease in the overall survival rate, but only the P value of MIR193BHG was less than 0.001 (). This suggested that the high expression of MIR193BHG was associated with the poor prognosis of HNSC.

Figure 6

The correlation analysis of lncRNAs in head and neck squamous cell carcinoma. (A-C) The co-expression analysis between hsa-miR-29c-3p and NOP14-AS1, LINC02323, and MIR193BHG. (D-F) The co-expression analysis between SFXN3 and NOP14-AS1, LINC02323, and MIR193BHG. (G-I) The differential expression analysis of NOP14-AS1, LINC02323, and MIR193BHG. (J-L) The overall survival analysis of NOP14-AS1, LINC02323, and MIR193BHG. (M) The competing endogenous RNA network of MIR193BHG–has-miR-29c-3p–SFXN3. lncRNAs, long non-coding RNAs.

According to the regulation mechanism of ceRNAs, MIR193BHG might be the upstream potential lncRNA of the hsa-miR-29c-3p–SFXN3 axis in HNSC. Therefore, MIR193BHG, hsa-miR-29c-3p, and SFXN3 were the nodes of the ceRNA network ().

SFXN3 and ICI

WilcoxTest was used to perform the ICI of SFXN3. The results of ICI are shown in . There were significant differences in naïve B-cells, plasma cells, gamma delta T-cells, resting natural killer (NK) cells, activated NK cells, monocytes, M0 macrophages, activated dendritic cells, and resting mast cells between the high-expressed and low-expressed SFXN3 groups in HNSC. To be more specific, SFXN3 was negatively correlated with naïve B-cells, plasma cells, activated NK cells, and monocytes, whereas it was positively correlated with resting NK cells, M0 macrophages, gamma delta T-cells, and activated dendritic cells.

Figure 7

The relationship of immune cell infiltration with SFXN3 level in head and neck squamous cell carcinoma. SFXN3, sideroflexin 3.

Correlation analysis between GPRIN1 and immune cell biomarkers

A total of 21 marker genes from 7 immune cells were assessed for the co-expression condition with SFXN3 using Wilcoxon signed-rank test, as shown in . M1 macrophages, neutrophils, B-cells, and CD8+ T-cells had a positive correlation with SFXN3, while dendritic cells, CD4+ T-cells, and M2 macrophages had a negative correlation. In particular, it was hard to define the correlation of neutrophils given their different co-expression condition in two marker genes. It is more likely these are the correlative immune cells when all the marker genes of each immune cell meet the criteria (P<0.01). Therefore, M2 macrophages, B-cell, and CD4+ T-cell were thought to be the correlative immune cells.

Table 2

Correlation analysis between SFXN3 and immune cell biomarkers

Immune cell	Gene	Cor	P value
B cell	CD19	−0.209299517	2.24E-06
B cell	CD79A	−0.192772558	1.42E-05
CD4+ T cell	CD4	0.193381259	1.33E-05
CD8+ T cell	CD8B	−0.039715079	0.374560864
CD8+ T cell	CD8A	−0.009689175	0.828497416
Dendritic cell	NRP1	0.493238342	<2.2E-16
Dendritic cell	HLA-DPA1	0.172679551	0.000103211
Dendritic cell	HLA-DRA	0.164365099	0.00022081
Dendritic cell	ITGAX	0.155430183	0.000480874
Dendritic cell	HLA-DQB1	0.146906381	0.000973472
Dendritic cell	HLA-DPB1	0.138821483	0.001838313
Dendritic cell	CD1C	0.095248441	0.032874875
M1 macrophage	NOS2	−0.279798149	1.76E-10
M1 macrophage	IRF5	−0.148080906	0.000885216
M1 macrophage	PTGS2	−0.071335946	0.110390528
M2 macrophage	MS4A4A	0.286088749	7.98E-11
M2 macrophage	CD163	0.270906232	7.99E-10
M2 macrophage	VSIG4	0.251295861	1.28E-08
Neutrophil	CEACAM8	−0.223416122	4.25E-07
Neutrophil	ITGAM	0.080407836	0.071870988
Neutrophil	CCR7	−0.016427949	0.713390365

Immune checkpoint analysis

PDCD1, CD274, and CTLA4 were analyzed to reveal the regulation between a single gene and an immune checkpoint () because of their tumor immune escape. They are important immune checkpoints responsible for tumor immune escape. TIMER and GEPIA were employed to analyze the relationship between SFXN3 and PDCD1, CD274, and CTLA4. As shown in , only CTLA4 was positively correlated with the gene SFXN3 (P=9.2E-03 in TIMER and P=2.2E-05 in GEPIA), while PDCD1 (P=9.14E-01 in TIMER and P=0.4 in GEPIA) and CD247 (P=6.64E-01 in TIMER and P=1.7E-07 in GEPIA) did not match the necessary criteria.

Figure 8

Correlation analysis between SFXN3 and immune checkpoints. (A,B) Correlation analysis between SFXN3 and cytotoxic T lymphocyte antigen 4 using TIMER 2.0 and GEPIA. (C,D) Correlation analysis between SFXN3 and programmed cell death 1 using TIMER 2.0 and GEPIA. (E,F) Correlation analysis between SFXN3 and programmed cell death ligand 1 (CD274) using TIMER 2.0 and GEPIA. SFXN3, sideroflexin 3; GEPIA, Gene Expression Profiling Interactive Analysis.

Discussion

HNSC is the sixth most common cancer worldwide, comprising TSCC, OSCC, LSCC, and NPC. The lack of effective treatment results in high mortality rates and a low quality of life of patients. To understand the molecular mechanism of HNSC tumorigenesis, it is necessary to attain important and practical clues in seeking more efficient therapeutic targets or specific prognostic biomarkers. Multi-omics data bioinformatics analysis is useful for researchers to better understand the molecular mechanism of HNSC and develop new drugs for treatment.

In this study, we downloaded the expression data and survival data of 33 types of cancer and conducted different expression and survival analyses. We only found that the gene SFXN3 was significantly highly expressed in HNSC and correlated with a poor prognosis. A previous study has proven that SFXN3 might act as a key role in cancerization. SFXN3, coding a mitochondrial membrane protein, was found to be a tumor marker of oral cancer patients by Japanese scientists (18).

Thirty-two upstream miRNAs were predicted to bind to SFXN3 by using the ENCORI database. After co-expression analysis, it was suggested that three miRNAs (hsa-miR-423-5p, hsa-miR-29c-3p, and hsa-miR-361-5p) were involved in HNSC. Differential expression and survival analysis suggested that only one miRNA, hsa-miR-29c-3p, was relatable to a better prognosis in HNSC. Previous studies have reported that hsa-miR-29c-3p could function as a tumor suppressor during tumorigenesis (19,20). MiR-29c-3p has different expression levels at different stages of LSCC tumor progression, and downregulation of miR-29c-3p is associated with a poor prognosis in patients with LSCC (19). Importantly, lower expression of has-miR-29c-3p in HNSC tumor tissue was associated with a higher tumor grade and worse relapse-free survival (20).

LncRNA has been reported to regulate gene expression by acting as a ceRNA molecule (21). According to the ceRNA hypothesis, the upstream regulatory lncRNA of the hsa-miR-29c-3p–SFXN3 axis is the potential pathogenic factor in HNSC, which is normally positively correlated with mRNA but negatively with hsa-miR-29c-3p. So, the ENCORI database was used to predict the upstream lncRNA of the hsa-miR-29c-3p–SFXN3 axis. After the differential expression analysis, survival analysis, and co-expression analysis, it was proposed that only MIR193BHG might be the regulatory lncRNA of the hsa-miR-29c-3p–SFXN3 axis in HNSC. A previous study demonstrated that MIR193BHG was a hypoxia-induced lncRNA and involved in the fine-tuning of cholesterol metabolism (22). It was also reported that MIR193BHG was related to the tumorigenesis of ovarian cancer, lung adenocarcinoma, and KIRC (23,24). As MIR193BHG acts as a common oncogene, we speculate that it might promote malignant cell proliferation and tumor growth in HNSC. Therefore, the MIR193BHG–hsa-miR-29c-3p–SFXN3 axis might serve as the potential regulatory pathway of HNSC.

The tumor microenvironment primarily consists of tumor cells, tumor-infiltrating immune cells, and the stromal component. An increasing number of clinical trials have revealed that the ICI is associated with immunotherapy and the prognosis of HNSC, but the full nature is still not elucidated (25). After comprehensive analyses, it was suggested that SFXN3 is positively correlated with resting NK cells, M0 macrophages, gamma delta T-cells, and activated dendritic cells. Moreover, SFXN3 is also positively correlated with some of the biomarkers of infiltrating immune cells. This data suggests that an ICI mechanism also exists in the process of HNSC. Correlation analysis between SFXN3 and immune checkpoints revealed that CTLA4 was positively correlated with SFXN3. Therefore, SFXN3 can be an immune checkpoint for HNSC, and it is expected to be a biomolecule that will improve the curative effect of immunotherapy. However, the results may need more confirmatory experiments to augment its stringency.

Conclusions

In summary, this study constructed a ceRNA network for HNSC based on the ceRNA hypothesis. Our study has identified that SFXN3, with a poor prognosis, exhibits abnormal higher expression in HNSC. The miRNA hsa-miR-29c-3p was negatively correlated with SFXN3, while the lncRNA MIR193BHG was positively correlated with SFXN3. The upstream regulatory mechanism of SFXN3 was identified and the action network of the MIR193BHG–hsa-miR-29c-3p–SFXN3 axis was constructed. SFXN3 was also associated with tumor ICI and immune checkpoint expression, which was suggested as a potential therapeutic target of HNSC.

The construction of a ceRNA network via pan-cancer analysis might provide some novel thoughts for HNSC treatment. The key nodes of the ceRNA network can reveal a new target to treat this tumor. Additionally, the analysis of the ICI, immune checkpoints, and immune cell biomarkers can also help the examination and therapy of HNSC more precisely.

The article’s supplementary files as