Hypoxia-related LncRNAs signature predicts prognosis and is associated with immune infiltration and progress of head and neck squamous cell carcinoma.

Background: Disclosing prognostic information is necessary to enable good treatment selection and improve patient outcomes. Previous studies suggest that hypoxia is associated with an adverse prognosis in patients with HNSCC and that long non-coding RNAs (lncRNAs) show functions in hypoxia-associated cancer biology. Nevertheless, the understanding of lncRNAs in hypoxia related HNSCC progression remains confusing.
Methods: Data were downloaded from TCGA and GEO database. Bioinformatic tools including R packages GEOquery, limma, pheatmap, ggplot2, clusterProfiler, survivalROC and survcomp and LASSO cox analysis were utilized. Si-RNA transfection, CCK8 and real-time quantified PCR were used in functional study.
Results: GEO data (GSE182734) revealed that lncRNA regulation may be important in hypoxia related response of HNSCC cell lines. Further analysis in TCGA data identified 314 HRLs via coexpression analysis between differentially expressed lncRNAs and hypoxia-related mRNAs. 23 HRLs were selected to build the prognosis predicting model using lasso Cox regression analyses. Our model showed excellent performance in predicting survival outcomes among patients with HNSCC in both the training and validation sets. We also found that the risk scores were related to tumor stage and to tumor immune infiltration. Moreover, LINC01116 were selected as a functional study target. The knockdown of LINC01116 significantly inhibited the proliferation of HNSCC cells and effected the hypoxia induced immune and the NF-κB/AKT signaling. Conclusions: Data analysis of large cohorts and functional experimental validation in our study suggest that hypoxia related lncRNAs play an important role in the progression of HNSCC, and its expression model can be used for prognostic prediction.

Background

Head and neck cancer (HNC) is one of the most common cancers worldwide with a poor five-year overall survival rate of 50% [1]. More than 90% of HNCs are histologically identified as head and neck squamous cell carcinoma (HNSCC). Although progress has been made and therapies have been enhanced, the survival rates of HNSCC have not been improved significantly [2]. It is necessary to identify novel biomarkers and predict models that can identify patients at high risk of death, stratify patients for different treatments and optimize the therapeutic effect.

Hypoxia is defined as a disparity between cellular oxygen demand and supply. Tumor hypoxia causes decreased sensitivity to radiotherapy and is thought to induce tumor progression and a more aggressive phenotype [3]. Therefore, the hypoxic status of a tumor may contribute to personalized treatment options that offer the best prognosis to an individual patient. Long non-coding RNAs (lncRNAs) are transcripts with a length greater than 200 nucleotides and contain no protein-coding sequence. LncRNAs play imperative roles in gene regulation at the transcriptional, posttranscriptional and chromosomal levels and are involved in diverse cellular processes [4,5]. Numerous studies have highlighted the irreplaceable roles of lncRNAs in tumor progression, lymph node metastasis, clinical stage and poor prognosis in HNSCC [6]. Recently, several lncRNAs have been reported to participate in hypoxia-associated cancer biology, implying a potential role of lncRNAs in maintaining cellular homeostasis and enabling adaptive survival under hypoxia [7,8]. However, the roles and mechanisms underlying hypoxia and lncRNAs in HNSCC regulation remain unclear. Thus, this study aimed to profile and evaluate the signature of hypoxia related lncRNAs in HNSCC and to develop an HRL prognostic model for predicting the prognosis of patients with HNSCC.

Material and methods

Data source

The transcriptome data of cell lines (GSE182734) were downloaded from GEO database (https://www.ncbi.nlm.nih.gov/geo/). The transcriptome data of HNSCC and normal samples were downloaded from The Cancer Genome Atlas database ((TCGA), https://cancergenome.nih.gov/). The relevant clinical materials were also downloaded from the database and were summarized in Table 1. Only patients with clinical follow-up/mortality>30 were included in this study. Hypoxia-related gene information was collected from the Molecular Signatures Database V7.2 (https://www.gsea-msigdb.org/gsea/msigdb, Hypoxia M10508, Hypoxia cancer M7547) [9]. Immune infiltration cell information was downloaded from CIBERSORTX (https://cibersort.stanford.edu/) [10].

Table 1

Baseline characteristics of HNSCC patients (n = 491).

Characteristic	n (%)
Age
≤ 65 years	321 (65.4%)
>65 years	170 (34.6%)
Unknown	0 (0%)
Sex
Female	130 (26.5%)
Male	361 (73.5%)
Unknown	0 (0%)
Tumor Grade
1&2	353 (71.9%)
3&4	119 (24.2%)
Unknown	19 (3.9%)
Tumor Stage
Ⅰ	25 (5.1%)
Ⅱ	69 (14.1%)
Ⅲ	78 (15.9%)
Ⅳ	251 (51.1%)
Unknown	68 (13.8%)
Pathologic T Stage
T0&1&2	174 (35.4%)
T3&4	262 (53.4%)
Unknown	55 (11.2%)
Pathologic N Stage
N0	167 (34.0%)
N1∼3	231 (47.1%)
Unknown	93 (18.9%)
Pathologic M Stage
M0	181 (36.9%)
M1	1 (0.2%)
Unknown	309 (62.9%)

Data processing

For all the transcriptome datasets, GSE182734 and TCGA-HNSCC, the R packages “GEOquery” and “limma” were used to identify the differentially expressed mRNAs and lncRNAs groups. The heatmap and volcano map were produced using R packages “pheatmap” and “ggplot2”. The R package “clusterProfiler” was used to perform functional enrichment analysis with Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO). The cutoff valve was set as p < 0.05. Besides, we conducted a Spearman correlation analysis between 139 hypoxia marker mRNAs and 1131 differentially expressed lncRNAs in TCGA-HNSCC dataset [11,12]. LncRNAs with a Spearman correlation coefficient >0.4 and a p value < 0.001 were defined as hypoxia related lncRNAs (HRLs).

Construction of the HRLs prognostic model

All data of the patients from TCGA-HNSCC were classified as a whole to the training dataset. Among the differentially expressed lncRNAs, we repeated the LASSO Cox regression fitting process of cross-validation evaluations 1000 times. Twenty-three genes with nonzero coefficient estimate greater than 900 times were selected for the prognostic model. We then calculated the HRL risk score of each sample according to the average coefficient and expression levels of each lncRNA. The formula is as follows:

HRLs risk score

Validation of the HRLs prognostic model

We divided the patients into the 1st and 2nd validation datasets at a ratio of 1:1 following the method described before [13]. In training and validation datasets, patients were allocated to the low- and high-risk groups based on the cutoff values evaluated by the R package “survminer”, and Kaplan–Meier survival curve analyses were then conducted. We further drew receiver operating characteristic (ROC) curves [14] and calculated the area under the ROC curves (AUC) via the R package “survivalROC” to assess the accuracy of the model. The R package “rms” was used to establish a nomogram to visualize the prognostic values of the HRL risk score. We also used the R package “survcomp” to calculate the concordance index (C-index) and to test the predictive accuracy of the nomogram [15].

Cells and reagents

SCC15 and FADU cell lines were obtained from the Chinese Academy of Sciences Cell Banke. LINC01116 siRNAs were synthesized by Sangon Biotech Co., Ltd. The sequences are as follows: si-Control sense: 5′- GAA ATT ATA GCG AAG GAC TGA -3′; si-LINC01116 #1 sense: 5′- GAG CAG TGT ATT AGA AGA CAA -3′; si-LINC01116 #2 sense: 5′- GGT AAT TGC CGC ATA GTG TAA -3′; si-LINC01116 #3 sense: 5′- GCA GTG TAT TAG AAG ACA A -3′. Transfections were performed according to the instructions of lipofectamine 2000 regent (Invitrogen). Cell counting Kit-8 (CCK8) (Dojindo Molecular Technologies) was used to evaluate cell proliferation according to the manufacturer's instructions.

Real time quantitative PCR

Trizol reagent (Invitrogen) and PrimeScript RT Reagent Kit with gDNA Eraser (Takara) were used for the isolation of RNA and to reverse transcribed RNA to cDNA. Roche A700 sequence detection machine and SYBR GREEN (Roche) were used to quantify LINC01116 levels. Has-U6 was used as an internal control for normalization. The following specific primer pairs were used: LINC01116: 5′- CGC TTT GCT GAA GAC GAG C -3′(sense), 5′- ATA TTG AAC TGA GCG GGG CT -3′ (antisense); Has-U6: 5′-GCT TCG GCA GCA CAT ATA CTA AAA T -3′(sense), 5′-CGC TTC ACG AAT TTG CGT GTC AT -3′ (antisense).

Statistical analysis

All statistical analyses were performed using R software (version 3.6.1). Spearman correlation analysis was used to identify HRLs. Kaplan–Meier survival curves and log-rank tests were used to evaluate the effects of HRLs on survival between the high- and low-risk groups. The Wilcoxon test was used to evaluate the significance of the immune cell subset distribution between the HNSCC and normal groups. The student's t-test was used to evaluate the significance of rt-PCR and CCK8 data. Two-tailed p values < 0.05 were considered statistically significant.

Results

Hypoxia treatment effects lncRNA regulation in vitro

To investigate the regulation of hypoxia on HNSCC, we first performed analysis in dataset GSE182734. In this dataset, HNSCC cells (LK0858, LK0863 and UT-SCC-14) were incubated in normal or hypoxic (1% O2) and were analyzed by Affymetrix transcriptome array 2.0. As shown in Fig. 1A, 1827 mRNAs and lncRNAs were confirmed differentially expressed after hypoxia induction (FoldChange >1.5 and p < 0.05). GO enrichment showed these genes were mainly enriched in non-coding RNA metabolic process and non-coding RNA processing (Fig. 1B, KEGG enrichment result in Fig. S1A). Therefore, we became more interested in the regulation of lncRNAs and believe that their interaction with hypoxia may be a key mechanism of HNSCC progression.

Fig. 1

A) Volcano diagram visualizing the differentially expressed mRNAs and lncRNAs between hypoxic and normal HNSCC cell line, names of top 10 were labeled (Data were obtained from GSE182734). B) GO analysis to explore enriched pathways basing on the differentially expressed genes between hypoxic and normal HNSCC cell lines. C-D) Heatmap (C) and volcano diagram (D) visualizing the differentially expressed mRNAs and lncRNAs between HNSCC and normal tissue specimens (Data were obtained from TCGA-HNSCC project, n = 546). E) Veen diagram showed the intersection of the significant lncRNAs from the two analyses above. F) The survival analysis of LINC01116 and RFPL1S in high and low lncRNA expressed groups basing on TCGA data.

Further, we explore this hypothesis in the TCGA data. We identified 2753 differentially expressed mRNAs and 1131 differentially expressed lncRNAs from the comparison between HNSCC (n = 491) and normal specimens (n = 44) (Fig. 1C and D) (|log2FoldChange| > 1.5 and p < 0.05). In addition, we found 6 lncRNAs showed significant expression changes (p < 0.05) both in the cell lines dataset and the TCGA-HNSCC (Fig. 1E). Among the 6 intersected lncRNAs, RFPL1S and LINC01116 displayed significant influence on HNSCC patient survival while other not (Fig. 1F, Fig. S1B).

We conducted a Spearman correlation analysis between 139 hypoxia marker mRNAs and differentially expressed lncRNAs and identified 314 lncRNA with p value < 0.001 and Spearman correlation coefficient >0.4 as HRLs. The correlations between hypoxia markers and HRLs were visualized by nets in Fig. 2 A and Fig. S2.

Fig. 2

A) Parts of the lncRNAs – hypoxia marker mRNAs association net. Red nodes represented the hypoxia marker mRNAs and green nodes represented the lncRNAs co-expressed with these genes (Spearman correlation coefficient >0.4 and p value < 0.001). LINC01116 was pointed out by a red arrow. B) The 23 optimal hypoxia-related lncRNAs (HRLs) identified by 1000 × LASSO Cox regression analyses. C-E) Survival-dependent receiver operating characteristics (ROCs) for predicting survival outcomes in the training dataset (C) and in the 1st (D) and 2nd (E) validation datasets. ROC curves for predicting the 2-, 3-, 4- and 5-year survival are represented by four different colors. F–H) Kaplan–Meier curves of the training dataset (F) and the 1st (G) and 2nd (H) validation datasets. I–K) The survival status distributions of samples in the training dataset (I) and the 1st (J) and 2nd (K) validation datasets. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

A total of 1000 lasso penalty analyses were further performed on 62 HRLs, which were significantly related to the OS of HNSCC (p ≤ 0.01), more than 900 times. As shown in Fig. 2B, 23 optimal prognostic HRLs were ultimately obtained to construct the risk model. The HRLs risk score for each patient was calculated by the following formula: risk score = (0.522 × LINC00922 expression) + (0.327 × DEPDC1-AS1 expression) + (0.279 × LINC00896 expression) + (0.227 × LINC01356 expression) + (0.185 × RP11-25L3.3 expression) + (0.127 × RP11-401O9.3 expression) + (0.123 × LINC02434 expression) + (0.082 × GSEC expression) + (0.021 × LINC00460 expression) + (0.018 × CTD-2510F5.4 expression) + (0.002 × MYOSLID expression) + (−0.024 × RP11-338N10.3 expression) + (−0.028 × RP11-44K6.2 expression) + (−0.066 × LA16c-390E6.4 expression) + (−0.091 × RP11-485G7.6 expression) + (−0.121 × RP11-93B14.10 expression) + (−0.317 × LINC01281 expression) + (−0.391 × LINC01352 expression) + (−0.413 × DDX11-AS1 expression) + (−0.531 × CTB-50L17.5 expression) + (−0.909 × LINC00567 expression) + (−0.991 × SIRPG-AS1 expression) + (−4.113 × LCT-AS1 expression).

All TCGA-HNSCC cases were classified as a whole to the training dataset (n = 491). In addition, we divided the patients into two validation datasets at a ratio of 1:1: the 1st validation dataset (n = 246) and the 2nd validation dataset (n = 245). In the training dataset, ROC curves were plotted and the AUCs for predicting 2-, 3-, 4- and 5-year survival were 0.72, 0.729, 0.74 and 0.737, respectively, demonstrating excellent prognostic capacity (Fig. 2C). Similar results were observed in the 1st (Fig. 2D) and 2nd (Fig. 2E) validation datasets.

The patients were then allocated to the high- and low-risk groups based on the cutoff of HRLs risk score. Kaplan–Meier curves showed poorer survival was observed in the high-risk group than in the low-risk group (p < 0.0001) (Fig. 2F) and the high-risk groups obviously showed a lower survival rate (Fig. 2I). Similar results were observed in the 1st and 2nd validation datasets (Fig. 2G and H, Fig. 2J and K). The HRLs risk score distributions in each dataset were observed (Figs. S3A–C). Meanwhile, with the increase of HRLs risk score, the expression levels of the protective lncRNAs (MYOSLID, LINC02434, LINC01356, LINC00460, RP11-25L3.3, LINC00922, CTD-2510F5.4, DEPDC1-AS1, RP11-401O9.3, LINC00896 and GSEC) decreased, while those of the risk lncRNAs (SIRPG-AS1, LINC01281, LCT-AS1, RP11-93B14.10, RP11-44K6.2, LINC00567, DDX11-AS1, CTB-50L17.5, RP11-485G7.6, LINC01352, LA16c-390E6.4and LA16c-390E6.4) increased (Figs. S3D–F).

We used univariate analysis to evaluate the relationships between clinical characteristics, HRLs risk scores and OS. The results showed that the HRLs risk score (p < 0.001) was significantly associated with prognosis in all datasets (Fig. 3A–C). The risk score was then confirmed as an independent prognostic factor via multivariate analysis (Fig. 3D–F). To predict the 3-, 4-, and 5-year survival outcomes of the patients, we then developed a nomogram based on candidate independent prognostic factors (age, gender, stage and HRLs risk score) (Fig. 3G).

Fig. 3

A-C) Univariate analysis the training dataset (A) and the 1st (B) and 2nd (C) validation datasets. D-F) Multivariate analysis of the training dataset (D) and the 1st (E) and 2nd (F) validation datasets. G) Nomogram established to predict the 3-, 4-, and 5-year survival outcomes of the patients.

Immune cells and signal pathways may include in hypoxia related lncRNAs regulation

Samples with advanced stages attained obviously higher HRLs risk scores (p = 0.002208, Fig. 4A). This finding suggests that the HRLs risk score also reflects tumor progression. We detected the differentially expressed genes between the high- and low-HRL risk groups. KEGG (Fig. S4) and GO enrichment analyses (Fig. 4B) were performed to identify the enriched pathways of these genes. It is worth mentioning that several immune-related signal pathways and genes were enriched. We further plotted the immune landscape of all the samples (Fig. 4C) and compared the difference in the number of immune cells between the low- and high-risk groups (Fig. 4D). Seventeen types of immune cells with differences in infiltration between the two groups were ultimately recognized, such as naïve B cells, CD8 T cells, CD4 naïve T cells, CD4 memory resting T cells and so on.

Fig. 4

A) Comparisons between the risk score and tumor stage in HNSCC. B) GO analysis explored enriched pathways basing on the differentially expressed genes between high- and low-risk score groups. C) Immune landscape of the patients with HNSCC. D) Correlation between the risk score and immune cell infiltration in HNSCC.

LINC01116 effect immune and NF-κB/AKT signaling and inhibits proliferation of HNSCC

To explore whether the regulation of hypoxia and lncRNAs is involved in the progression of HNSCC and whether the enriched pathways also play roles in it, we investigated the effects of LINC01116 in HNSCC progress. This choice mainly relies on: LINC01116 showed significant changes in both cell lines and TCGA datasets and its expression levels were associated to patient prognosis; LINC01116 is one of HRLs; Roles of LINC01116 in hypoxia and immune signaling has not been studied in HNSCC before. We first reviewed the expression changes of LINC01116 in cell lines cultured with hypoxia and normal conditions and in TCGA-HNSCC dataset (Fig. 5A). After that, SCC15, a tongue squamous cell carcinoma cell line, and FADU, a hypopharynx squamous cell carcinoma cell line, were used to detect the effects of LINC01116 knockdown in HNSCC. Rt-PCR confirmed that the expression of LINC01116 was remarkable down regulated by si-LINC01116 (Fig. 5B). The CCK8 assays showed that the LINC01116 knockdown inhibited the proliferation of SCC15 and FADU cells (Fig. 5C). As shown in Fig. 4D, hypoxia increased HIF1α, TNF-α and phosphorylated P65 while si-LINC01116 inhibited the induction, suggesting LINC01116 takes part in hypoxia induced immune and NF-κB signaling. Besides, MMP13, a hypoxia marker gene and the potential target of LINC01116, decreased by si-LINC01116 treatment under both normal and hypoxia conditions (Fig. 4E). We also found the activations of PI3K and ATK signaling were upregulated by hypoxia while inhibited by si-LINC01116.

Fig. 5

A) The expression fold changes of LINC01116 in TCGA data and cell lines data. B) Real-time quantitative PCR to evaluate LINC01116 levels in SCC15 and FADU cell lines transfected with si-LINC01116 or si-Control. C) The effect of LINC01116 knockdown on the proliferation of SCC15 and FADU cells was determined by CCK8 assay. D) The effects of si-LINC01116 and hypoxia on HIF1α and immune related proteins IL-1β, TNF-α and pP65 were detected by western blotting, protein samples derive from the same experiment and the blotting were processed in parallel. E) The effects of si-LINC01116 and hypoxia on MMP13, pPI3K, PI3K, pAKT and AKT were detected by western blotting, protein samples derive from the same experiment and the blotting were processed in parallel.

Discussion

Effective prognostic models and biomarkers are urgently needed to benefit HNSCC patients. In the present study, we analyzed the signature of hypoxia-related lncRNAs in HNSCC and identified an effective HRL prognostic model based on these lncRNAs for predicting HNSCC outcomes. Besides, we initially identified LINC01116 play roles in hypoxia related HNSCC progression.

HNSCC prognostic models have emerged in recent years. Transcriptional signatures were the most common models reported before [16,17], and filter conditions such as metabolism and immune checkpoints were also included [18,19]. The present study clarified that hypoxia regulation and lncRNA regulation may be two interactive factors closely related to HNSCC progression basing on cell lines data. Constructing a hypoxia related lncRNAs predicting model may be a good choice for HNSCC. Currently, models based on lncRNAs are still limited. A 4-lncRNA signature based on GEO data, an immune-related lncRNA signature prognostic model including 13 lncRNAs and a three autophagy-related lncRNA model for predicting the survival of HNSCC has been presented [[20], [21], [22]]. Consistent with the last two models, we used data from TCGA, which contains numbers of HNSCC patients and complete clinical and survival data, suggesting better statistical power than a model constructed using a small cohort. Our model took hypoxia, which plays an obviously significant role in poor patient prognosis and therapy resistance of HNSCC, as a filtrating factor and may be useful as a supplement to the tumor stage for better stratifying patients and for providing a more individualized treatment.

Tumor hypoxia can lead to resistance to chemoradiotherapy and targeted therapy, thus predisposing patients to tumor metastases, and contributes to altered metabolism and genomic instability [23]. LncRNAs have also been reported to be involved in the development of HNSCC by regulating the hypoxia pathway [7,8,24,25]. Our results that the differently expressed genes induced by hypoxia mainly enriched on ncRNA process pathways supported the ncRNA related regulations may be the main meachnisms underlying hypoxia. Besides, the coexpression analysis identified a total of 314 lncRNAs were significantly coexpresed with hypoxia marker genes, indicating the crucial relationship between hypoxia and lncRNAs in HNSCC. Some of these HRLs have been found to regulate the progression of other cancers, such as LINC00922 in liver cancer [26], breast cancer [27] and lung cancer [28] and DDX11-AS1 in esophageal carcinoma [29], liver cancer [30] and glioma [31]. However, to date, research on lncRNAs in hypoxia and HNSCC is rare, indicating a field worthy of in-depth study. We also noticed that some of these HRLs have been used in other prognostic prediction models, such as CTD-2510F5.4 in hepatocellular carcinoma [32], DEPDC1-AS1 in lung cancer [33], and LINC01352 in endometrial cancer [34]. More importantly, LINC02434 has been used in a three-lncRNA overall survival of HNSCC patients prediction model [35]. These overlapping findings with previous studies contribute to the reliability of our model.

Tumor hypoxia also changes the surrounding tumor microenvironment, leading changes in immune cell infiltration [36]. KEGG and GO functional enrichment analyses showed that several immune-related pathways were enriched, such as cytokine interactions, T cell-related actions, and lymphocyte differentiation. Seventeen types of immune cells were found to be differentially infiltrated in HNSCC and are closely associated with tumorigenesis and progression. The poor prognosis associated with HRLs may be rely on the decrease or resting of naïve B cells, CD8 T cells, CD4 memory T cells, follicular helper T cells, regulatory T cells, NK cells and macrophages and the activation of dendritic cells, mast cells and neutrophils. The mechanisms underlying hypoxia related lncRNAs and immune cells and pathways are worth to be studied in further studies. Here, combining the HRLs analysis of TCGA data and the lnRNAs analysis of hypoxic and normal cell lines, we identified LINC01116 as a valuable mechanism study target. The promotion of LINC01116 in tumor progress has been reported before in kinds of cancer including HNSCC and the underlying mechanisms had been reported as regulating the epithelial-mesenchymal transition process before [37,38]. The present study preliminarily examined the role of LINC01116 in the proliferation of HNSCC, as well as immune factors IL-1β, TNF-α and tumor-related pathways NF-κB，PI3K and ATK. The effects of LINC01116 performed on hypoxia and immune related signaling has not been reported before, so out observations were novel and supported further studies about the hypoxia induced and LINC01116-immune mediated HNSCC progress. In all, our study still has several limitations. We did not find an available independent lncRNA dataset to verify the helpfulness of our prognostic model. The validity of our prognostic model should be further evaluated in a large number of clinical samples, and the underlying mechanisms of lncRNAs influence the prognosis of HNSCC should be explored by more detailed wet experiments.

Conclusion

The present study detected the hypoxia-related lncRNA signature of HNSCC and provided an effective assessment model to improve the prognostic prediction of HNSCC patients. HRLs, for example LINC01116, will be valuable targets for investigating the tumorigenesis of HNSCC and discovering personalized treatment strategies for HNSCC.

Funding

The authors declare that no funds, grants, or other support were received during the preparation of this manuscript.

Author contributions

All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by ML, LL and SM. The first draft of the manuscript was written by MW and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.

Availability of data and materials

Data used in this study can be found here: the Cancer Genome Atlas database (TCGA), https://cancergenome.nih.gov/), the Molecular Signatures Database V7.2 (https://immport. niaid.nih.gov), CIBERSORT (https://cibersort.stanford.edu/), the STRING database (https://STRING-db.org/), the Encyclopedia of RNA Interactomes (ENCORI) (http://starbase.sysu.edu.cn/index.php), the GEO database (https://www.ncbi.nlm.nih.gov/geo/).

Ethics approval and consent to participate

Not applicable.

Consent to publish

Not applicable.

Declaration of competing interest

The authors have no relevant financial or non-financial interests to disclose.