Abstract
Background and purpose: Radioresistance remains a major reason of
radiotherapeutic failure in esophageal squamous cell carcinoma (ESCC).
Our study is to screen the immune-related long non-coding RNA
(ir-lncRNAs) of radiation-resistant ESCC (rr-ESCC) via Gene Expression
Omnibus (GEO) database and to construct a prognostic risk model.
Methods: Microarray data ([33]GSE45670) related to radioresistance of
ESCC was downloaded from GEO. Based on pathologic responses after
chemoradiotherapy, patients were divided into a non-responder (17
samples) and responder group (11 samples), and the difference in
expression profiles of ir-lncRNAs were compared therein. Ir-lncRNA
pairs were constructed for the differentially expressed lncRNAs as
prognostic variables, and the microarray dataset ([34]GSE53625) was
downloaded from GEO to verify the effect of ir-lncRNA pairs on the
long-term survival of ESCC. After modelling, patients are divided into
high- and low-risk groups according to prognostic risk scores, and the
outcomes were compared within groups based on the COX proportional
hazards model. The different expression of ir-lncRNAs were validated
using ECA 109 and ECA 109R cell lines via RT-qPCR.
Results: 26 ir-lncRNA genes were screened in the [35]GSE45670 dataset
with differential expression, and 180 ir-lncRNA pairs were constructed.
After matching with ir-lncRNA pairs constructed by [36]GSE53625, six
ir-lncRNA pairs had a significant impact on the prognosis of ESCC from
univariate analysis model, of which three ir-lncRNA pairs were
significantly associated with prognosis in multivariate COX analysis.
These three lncRNA pairs were used as prognostic indicators to
construct a prognostic risk model, and the predicted risk scores were
calculated. With a median value of 2.371, the patients were divided
into two groups. The overall survival (OS) in the high-risk group was
significantly worse than that in the low-risk group (p < 0.001). The
1-, 2-, and 3-year prediction performance of this risk-model was 0.666,
0.702, and 0.686, respectively. In the validation setting, three
ir-lncRNAs were significantly up-regulated, while two ir-lncRNAs were
obviouly down-regulated in the responder group.
Conclusion: Ir-lncRNAs may be involved in the biological regulation of
radioresistance in patients with ESCC; and the prognostic risk-model,
established by three ir-lncRNAs pairs has important clinical value in
predicting the prognosis of patients with rr-ESCC.
Keywords: radioresistance, esophageal squamous cell carcinoma, lncRNA,
prognostic model, bioinformatics
Introduction
Esophageal cancer (EPC) is one of the most lethal tumors in China and
worldwide ([37]Sung et al., 2021). Esophageal squamous cell carcinoma
(ESCC) is the major pathologic type in Chinese population, accounting
for more than 90% cases, and more than 45%–60% of them were diagnosed
at an advanced stage ([38]Yang et al., 2018). In general, the prognosis
of EPC is very poor, the cancer-specific mortality rate ranks fourth,
and the 5-year overall survival rate is less than 30% ([39]Chen et al.,
2016; [40]Zeng et al., 2018). Radiotherapy is one of the main
treatments for EPC, especially in patients at an advanced stage.
However, after radiotherapy with/without chemotherapy, the total
objective response rate (ORR) is only 60%–76%, and radiotherapy
resistance is an important reason for the local failure of radiotherapy
for EPC ([41]Chen et al., 2021). Although multiple factors may affect
the radiosensitivity, the specific mechanisms of long non-coding RNA
(lncRNAs) in radioresistance are still worth exploring. The outcomes of
radiotherapy are heterogeneous, and no clinical or pathological method
could predict tumor response of radiotherapy.
LncRNA referred to a type of RNAs whose length is greater than 200 bp
and does not encode or translate proteins after transcription. LncRNAs
are involved in various physiological and pathological processes, such
as the cellular replication, transcription, translation and so on. In
the process of gene expression and transcription, about 70% of human
genes are regulated by lncRNAs ([42]Batista and Chang, 2013; [43]Shi et
al., 2013). It’s reported that lncRNAs are widely involved in the
biological processes of cancer, including tumor radioresistance ([44]Li
and Chen, 2013). Many studies believed that the prognosis of tumors is
closely related to the tumor microenvironment (TME), and lncRNAs may
participate in the regulation of TME through molecular biological
functions such as chromatin modification, transcription, and
post-transcriptional processing ([45]Sounni and Noel, 2013; [46]Singh
et al., 2016). In addition, some lncRNAs are also involved in the
regulation of immune function in the TME, where may involve various
types of cells and cytokines. Such lncRNAs are often referred to as
immune-related lncRNAs (ir-lncRNAs) ([47]Bremnes et al., 2016; [48]Chen
M. M. et al., 2020). In order to explore whether ir-lncRNAs are
involved in the biological process of radio-resistance in patients with
ESCC, our study analyzed the differences in the expression of
ir-lncRNAs between ESCC-patients with radioresistance and complete
tumor remission based on gene expression profiling data in the Gene
Expression Omnibus (GEO) database. Ir-lncRNA pairs, which had high
correlation with immune genes, were selected to establish a prognostic
risk-model to explore the value of ir-lncRNAs in predicting the
prognosis of ESCC patients with radioresistance, and further to find
potential biomarkers of radioresistance.
Materials and methods
Dataset collection and preparation
Using “esophageal cancer, esophagus cancer, and radioresistance” as
keywords, mRNA and lncRNA expression profiling data related to
radioresistance of esophageal cancer was searched in the GEO database
([49]https://www.ncbi.nlm.nih.gov/). Two public esophageal cancer
microarray profiling datasets ([50]GSE45670 and [51]GSE53625) were
downloaded from GEO and were finally selected for data mining. The
[52]GSE45670 microarray profiling datasets were provided by Wen et al.,
where gene expression analyses were performed on pretreatment cancer
biopsies from 28 ESCCs who received neoadjuvant chemoradiotherapy (CRT)
and surgery and 10 normal esophageal epithelia using Affymetrix U133
Plus 2.0 arrays ([53]Wen et al., 2014). After preoperative
chemoradiotherapy among 28 ESCCs, complete remission of tumor occurred
in 11 patients, who were divided into responder group, while no obvious
tumor regression occurred in the other 17 patients, who were considered
to be radioresistant and divided into non-responder group. The average
age of the patients in the non-responder group was (55.65 ± 5.53) years
old, of whom 16 patients were male (accounting for 94.11%); four
patients were in T2N1M0 stage, and 13 patients were in T3N1M0 stage;
while average age of the patients in the responder group was (57.45 ±
6.80) years old, of whom nine patients were male (accounted for
81.82%); four patients were in T2N1M0 stage and seven patients were in
T3N1M0 stage.
The [54]GSE53625 microarray profiling datasets were provided by Li et
al.([55]Shi et al., 2016; [56]Li et al., 2017; [57]Liu et al., 2020),
including 179 patients with esophageal cancer. The profiling datasets
were derived from cancer tissue and adjacent normal tissue. The dataset
contained detailed clinical information of 179 ESCC patients that can
be used for prognostic validation analysis. Both datasets were
annotated with the platform files provided by GEO to obtain Ensemble
ID, and then were annotated with the “org.Hs.eg.db” package to obtain
gene symbols, thus ensuring that the two datasets have similar
annotation conditions.
Screening, differential expression analysis and matched pairs of ir-lncRNAs
The immune gene list was download from the ImmuPORT database
([58]https://www.immport.org/). Subsequently, the subsets of immune
genes or lncRNAs were extracted from the [59]GSE45670 dataset and
[60]GSE53625, respectively. The correlation test was performed on the
immune gene expression matrix and the lncRNA expression matrix, and the
related lncRNAs were screened as ir-lncRNAs. The screening criteria
were: correlation coefficient ρ ≥ 0.4 or ρ ≤ −0.4 and p ≤ 0.0001.
Differential expression analysis (DEA) was performed on eligible
ir-lncRNAs via the “limma” and “SVA” packages, and the expression
difference was defined as: the absolute value of the log2 value (fold
change) of the difference fold was greater than (mean ± 2 times the
standard deviation of expression); p value ≤ 0.05. The differentially
expressed lncRNAs were paired for each other. Taking the
LINC01121|FAM167A-AS1 gene as an example, if the expression of
LINC01121 gene was greater than that of FAM167A-AS1 gene, the matched
index was recorded as 1, otherwise as 0. The ir-lncRNA pairs were
introduced as independent variables into the COX prediction model to
establish a prognostic risk-model.
Construction of the predictive model
The clinical survival data of ESCC samples were extracted from the
[61]GSE53625 dataset, and the ir-lncRNA pairs at the common
intersection between the [62]GSE45670 dataset and the [63]GSE53625
dataset were considered as components of the survival data.
Univariate and multivariate Cox regression survival analyses were
performed using the “survival” package in R 4.1.2 software to screen
ir-lncRNAs with significant impact on prognosis. The least absolute
shrinkage and selection operator (LASSO) regression analysis was
carried out to narrow down the prognostically significant ir-lncRNA
pairs. An ir-lncRNAs pairs-based risk prediction model was established
with ir-lncRNAs of statistically significant differences in both
univariate and multivariate Cox regression survival analysis.
The model building formula is ([64]Huang et al., 2021):
[MATH:
Riskscore=∑i=1<
/mn>n(ln
cRNApai<
/mi>rsi×coe
mi>fi).
:MATH]
Where n is the counts of ir-lncRNA pairs, and
[MATH:
lncRNAp
airsi :MATH]
and
[MATH: coefi
:MATH]
represent the related matched results (values of 1 or 0) and
coefficients of modeled ir-lncRNA pairs, respectively.
After the riskscores of all the included samples in [65]GSE53625
dataset were calculated, the samples were divided into low-risk group
and high-risk group according to the median value of riskscores.
Furthermore, the Kaplan-Meier survival curve was applied to analyze the
difference of survival prognosis between the high- and low-risk groups
with Survival package. The ROC curves of 1-, 2-, and 3-year of overall
survival were plotted with ROC package. Univariate and multivariate Cox
regression survival analyses were performed with ESCC baseline data and
riskscores via the “survival” package in R 4.1.2 software,
respectively, to explore the independent prognostic factors of ESCC.
Differential analysis of riskscore for different clinical characteristics
The clinical characteristics such as gender, age, tumor TNM stage, etc.
of ESCC samples were extracted from [66]GSE53625 dataset and combined
with the corresponding risk scores. Wilcoxon rank sum test was
performed via the “ggpubr” package and the boxplots was drawn to show
the correlations between risk scores and clinical characters.
Differential expression analysis of immune-related genes in different
radiosensitivity groups
Based on those ir-lncRNAs screened from the [67]GSE53625 dataset with
prediction model, the immune-related genes in the [68]GSE45670 dataset
were reversely extracted, and the differential expression analysis of
immune-related genes were compared between the non-responder group and
the responder group. A heat map and a volcano plot were drawn.
Gene ontology, Kyoto encyclopedia of genes and genomes functional enrichment
analysis
For the purpose of exploring the molecular mechanism of ir-lncRNAs and
ir-genes related to radio-resistance in ESCC, we implemented gene
ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG)
functional enrichment analysis via the “clusterProfiler” R package. In
these analyses, a p value < 0.05 was considered statistically
significant. GO enrichment analysis was composed of cellular components
(CC), molecular functions (MF), and biological processes (BP), which
showed the biological functions of genes at different levels,
respectively. KEGG pathway enrichment analysis was used to evaluate the
enrichment degree of genes in different pathways. The function of
selected ir-genes whose expression were significantly different between
the non-responder group and the responder group could indirectly
speculate on the biological roles and mechanisms of immune-related
lncRNAs.
lncRNAs isolation, cDNA synthesis, and RT-qPCR
Cellular total RNAs were isolated via TRIzol reagent (Thermo, United
States) from ECA-109 cell lines and ECA-109R cell lines, where ECA-109R
were identified as non-responder esophageal cancer and ECA-109 was
identified as responder esophageal cancer. To increase the specificity
of the real-time quantitative PCR (RT-qPCR), first-strand cDNA was
synthesized from 1 mg total RNA via RevertAid First Strand cDNA
Synthesis Kit (Invitrogen, China). The relative expression level of the
target lncRNAs was detected by three-step RT-qPCR via the Fast SYBR
Green Master Mix (Applied Biosystems Inc., CA, United States). The
cycling conditions were 5 min of pre-degeneration at 95°C and 30 s of
denaturation (polymerase activation) at 95°C followed by 40 cycles of
annealing of primers at 95°C for 5 s and extension of primers at 60°C
for 30 s. GAPDH was applied as an internal reference control. The
relative expression level was calculated by the relative quantification
2^−ΔΔCT method. The experiment was repeated three times and the primer
sequences were listed in [69]Supplementary Table S1.
Statistical analysis
All statistical analyses were performed via R software (version 4.1.2).
The R packages involved the “Biobase” package, the “GEOquery” package,
the “survival” package, the “ggpubr” package, and so on. The one-way
ANOVA, Student’s t-test (2-tailed) methods or Wilcoxon rank sum test
were applied to assess the statistical significance. Quantitative data
were shown as the mean ± SE. A p value < 0.05 was considered
statistically significant.
Results
Screening results of ir-lncRNAs in [70]GSE45670 and [71]GSE53625 datasets
Flow chart of data collection and analysis is shown in [72]Figure 1.
According to the screening criteria, a total of 681 ir-lncRNAs were
screened in the [73]GSE45670 dataset. Further DEA showed that 12
ir-lncRNAs were up-regulated and 14 ir-lncRNAs were down-regulated in
the non-responder group (p < 0.05). The results are indicated in
[74]Table 1 and [75]Figure 2.
FIGURE 1.
[76]FIGURE 1
[77]Open in a new tab
Flow chart of data collection and analysis.
TABLE 1.
Differential expression results of immune-related lncRNAs in
[78]GSE45670.
Genes log FC Ave Expr t p value adj p value Change
LINC01121 −1.196 3.898 −3.679 0.001 0.318 Down
LINC00592 −1.306 7.669 −2.93 0.006 0.589 Down
CTD-3080P12.3 −1.076 3.852 −2.595 0.015 0.589 Down
WAKMAR2 −1.041 7.665 −2.569 0.016 0.589 Down
H19 −1.902 8.538 −2.536 0.017 0.589 Down
DLGAP4-AS1 −0.982 5.319 −2.497 0.018 0.589 Down
SCAT1 −1.459 5.612 −2.487 0.019 0.589 Down
LOC101928557 −1.139 4.483 −2.417 0.022 0.589 Down
PURPL −1.619 3.829 −2.403 0.023 0.589 Down
IQCF5-AS1 −1.061 3.662 −2.302 0.029 0.641 Down
ELFN2 −1.126 6.26 −2.291 0.029 0.641 Down
MIR124-2HG −0.991 4.119 −2.217 0.034 0.69 Down
LINC01102 −1.018 2.774 −2.088 0.046 0.825 Down
LOC101928389 −1.09 3.768 −2.048 0.05 0.825 Down
ADAMTS9-AS2 1.714 5.225 4.129 0 0.187 Up
SOX2-OT 1.731 5.417 2.935 0.006 0.589 Up
GRK3-AS1 1.601 3.371 2.864 0.008 0.589 Up
DELEC1 1.11 3.089 2.691 0.012 0.589 Up
FAM167A-AS1 1.391 3.184 2.536 0.017 0.589 Up
ZNF503-AS1 1.112 7.336 2.532 0.017 0.589 Up
RNF217-AS1 0.984 5.366 2.486 0.019 0.589 Up
MGC12916 1.079 4.97 2.446 0.021 0.589 Up
LOC101927798 1.036 5.289 2.409 0.022 0.589 Up
LINC00551 1.432 6.052 2.391 0.023 0.589 Up
LINC00942 1.545 6.705 2.109 0.044 0.824 Up
FSIP2-AS2 1.335 5.226 2.081 0.046 0.825 Up
[79]Open in a new tab
FIGURE 2.
[80]FIGURE 2
[81]Open in a new tab
Differential expression results of immune-related lncRNAs in
[82]GSE45670 cohort. (A) Heatmap, (B) Volcano map.
The 26 differentially expressed ir-lncRNAs genes constituted 325
ir-lncRNA pairs, of which 180 ir-lncRNA pairs were extracted according
to the matching rate [pairRaito in the range of (0.2–0.8)], and the
data are shown in [83]Supplementary Table S2. Subsequently, we used the
same method to construct the ir-lncRNA pairs in the [84]GSE53625
dataset. Finally, 74 ir-lncRNA pairs successfully matched between the
[85]GSE45670 and the [86]GSE53625 dataset after intersection. The
detailed information of the 74 ir-lncRNA pairs is shown in
[87]Supplementary Table S3.
Clinical characteristics of esophageal squamous cell carcinoma samples in the
[88]GSE53625 dataset
The detailed clinical characteristics of [89]GSE53625 ESCC are shown in
[90]Table 2.
TABLE 2.
Clinical characteristics of esophageal squamous cell carcinoma samples
in [91]GSE53625.
Items Overall (n = 179, %) Items Overall (n = 179, %)
Gender Age
Female 33 (18.44) Mean (SD) 53.9 (9.03)
Male 146 (81.56) Median (min, max) 59.7 (36.0,82.0)
Tobacco Alcohol
Yes 114 (63.69) yes 106 (59.22)
T stage N stage
T1 12 (6.7) N0 83 (46.37)
T2 27 (15.08) N1 62 (34.64)
T3 110 (61.45) N2 22 (12.29)
T4 30 (16.76) N3 12 (6.7)
TNM stage Tumor grade
Stage I 10 (5.59) Poorly 49 (27.37)
Stage II 77 (43.02) Moderately 98 (54.75)
Stage III 92 (51.4) Well 32 (17.88)
[92]Open in a new tab
Univariate and multivariate COX regression analysis of differential
ir-lncRNAs
The results of univariate COX analysis showed that there were six
ir-lncRNA pairs of statistical significance associated with survival
outcomes (p < 0.05). In order to verify the positive variables obtained
by univariate COX regression, LASSO regression (Least absolute
shrinkage and selection operator) was further used to screen
significant ir-lncRNAs as independent variables for multivariate COX
analysis. As a result, five ir-lncRNA pairs, that is,
LINC01121|FAM167A-AS1, ADAMTS9-AS2|MGC12916, MIR124-2HG|FAM167A-AS1,
LINC00942|ADAMTS9-AS2, and PURPL|FAM167A-AS1, were the most powerful
explanatory set of ir-lncRNAs as potential independent variables in
multivariate COX model. A further multivariate COX regression was
performed, and the results showed that three ir-lncRNA pairs were
statistically significant (p < 0.05), namely LINC01121|FAM167A-AS1,
ADAMTS9-AS2|MGC12916, and MIR124-2HG|FAM167A-AS1. The above analysis
revealed that the three ir-lncRNA pairs could serve as independent
prognostic factors for ESCC. Since the hazard ratios of the three
ir-lncRNAs were all greater then 1, they were considered to be risk
factors for the prognosis of ESCC (as shown in [93]Table 3 and
[94]Figure 3). Lasso regression road map and selection map were shown
in [95]Supplementary Figures S1, S2. ROC curves were used to assess the
predictive ability and accuracy of the predictive models based on the
three ir-lncRNA pairs. The results showed that the best cut-off point
was 2.371 and the area under the curve (AUC) at 1-, 2-, and 3-year were
0.666, 0.702, and 0.686, respectively, as shown in [96]Figures 4A,B.
TABLE 3.
Univariate and multivariate COX regression analysis of differential
ir-lncRNAs.
Ir-lncRNA pairs Model β se HR HR.95L HR.95H p
LINC01121|FAM167A-AS1 Univariate 0.789 0.237 2.201 1.385 3.5 0.001
ADAMTS9-AS2|MGC12916 Univariate 0.664 0.288 1.943 1.105 3.418 0.021
MIR124-2HG|FAM167A-AS1 Univariate 0.567 0.196 1.764 1.202 2.588 0.004
LINC00942|ADAMTS9-AS2 Univariate −0.84 0.424 0.432 0.188 0.991 0.048
PURPL|FAM167A-AS1 Univariate 0.598 0.249 1.819 1.117 2.961 0.016
ZNF503-AS1|H19 Univariate 2.662 1.041 14.327 1.863 110.187 0.011
LINC01121|FAM167A-AS1 Multivariate 0.687 0.252 1.989 1.213 3.259 0.006
ADAMTS9-AS2|MGC12916 Multivariate 0.863 0.293 2.371 1.335 4.211 0.003
MIR124-2HG|FAM167A-AS1 Multivariate 0.419 0.209 1.52 1.009 2.291 0.045
[97]Open in a new tab
Note: HR.95L: lower limit of 95% confidence interval; HR.95H: upper
limit of 95% confidence interval.
FIGURE 3.
[98]FIGURE 3
[99]Open in a new tab
Forest plots of Cox regression analysis for ir-lncRNA pairs: (A)
Univariate Cox regression analysis, (B) Multivariate Cox regression
analysis.
FIGURE 4.
[100]FIGURE 4
[101]Open in a new tab
Time ROC curves on overall survival prediction in the [102]GSE53625
cohort (A) Time ROC curves with cutoff value, (B) Time ROC curves at
one year, two year and three year.
Additionally, [103]Supplementary Figure S3 provides differential
distribution of riskscore with different clinical characters.
Construction of ir-lncRNAs-based risk model and survival analysis of high-
and low-risk groups
The model formula obtained from multivariate COX regression analysis
was:
[MATH:
Riskscore=0.687×
LINC01121|FAM167A−AS1+
mo>0.863×ADAMTS9−AS2|MGC129
16+0.419×MIR124
−2HG|FAM
167A−AS1
:MATH]
The median riskscore was 1.989, and according to the median riskscore,
179 cases were divided into the high- (N = 76) and low-risk (N = 103)
groups in the [104]GSE53625 cohort. The Kaplan-Meier survival curves
shows that patients in the high-risk group had a worse OS than patients
in the low-risk group with more death (p < 0.001, [105]Figure 5). The
1-, 3-, and 5-year survival rates in the low-risk group were higher
than those in the high-risk group.
FIGURE 5.
FIGURE 5
[106]Open in a new tab
Kaplan-Meier curves for the overall survival of patients in the high-
and low-risk group.
Differential expression analysis of ir-lncRNAs in high- and low-risk groups
Based on the aforementioned results, independent sample t-test were
further applied to assess the expression of ir-lncRNAs in different
riskscore groups. The distribution of risk scores is presented in
[107]Figure 6. Among the five significant lncRNAs, only FAM167A-AS1 had
a decreased expression in the high-risk group, suggesting that
FAM167A-AS1 is a protective factor.
FIGURE 6.
[108]FIGURE 6
[109]Open in a new tab
Different expression of five ir-lncRNAs between the low- and the
high-risk group.
Univariate and multivariate COX regression analysis of clinical
characteristics
The clinical characteristics parameters in [110]Table 2 combined with
riskscore were selected to perform univariate and multivariate COX
regression analysis. Univariate COX regression analysis demonstrated
that age (HR = 1.681, 95% CI = 1.147–2.463, and p = 0.008), N-stage
(Lymph node staging, N2 vs. N0, HR = 2.051, 95% CI = 1.137–3.702, and p
= 0.017; N3 vs. N0, HR = 2.973, 95% CI = 1.426–6.200, and p = 0.004),
TNM-stage (Stage III vs. Stage I, HR = 3.626, 95% CI = 1.138–11.548,
and p = 0.029) and risk score (HR = 1.394, 95% CI = 1.225–1.587, and p
< 0.001), which were all negative prognostic factors of OS in the
[111]GSE53625 cohort ([112]Supplementary Table S4; [113]Figure 7A).
After adjusting for age, tumor grade, N-stage and TNM-stage,
multivariate Cox analysis demonstrated that only risk score was a
negative prognostic factor of OS (HR = 1.305, 95% CI = 1.139–1.495, and
p < 0.001) ([114]Supplementary Table S5; [115]Figure 7B).
FIGURE 7.
[116]FIGURE 7
[117]Open in a new tab
Forest plots of clinical characteristics and risk score for univariate
and multivariate COX regression analysis (A) Univariate Cox regression
analysis, (B) Multivariate Cox regression analysis.
CIBERSORT immune infiltration analysis
To explore the difference of tumor immunity landscape between patients
in the non-responder group and patients in the responder group, the
CIBERSORT algorithm was utilized to evaluate immunity infiltration in
the GSE 45670 dataset. The main results are shown in [118]Figures 8A,B.
According to the 22-classification method ([119]Newman et al., 2019),
the proportions of infiltrating activated mast cells were significantly
higher in the responder group (p < 0.05). In addition, the proportions
of infiltrating macrophages. M0 were also increased in the responder
group, although the difference was not statistically significant (p >
0.05). The 4-classification method showed that the infiltration level
of macrophages in the responder group was significantly higher than
that in the non-responder group (p < 0.05) ([120]Li B. et al., 2019).
FIGURE 8.
[121]FIGURE 8
[122]Open in a new tab
Differences of infiltrating immune cell types between the non-responder
group and the responder group of CIBERSORT in [123]GSE45670 cohort. (A)
22-classification method, (B) 4-classification method.
Differential expression analysis of immune genes
The correlation coefficient ρ = ±0.4 and p < 0.05 were used to screen
the immune genes related to the above five ir-lncRNAs, and a total of
137 immune genes were screened, of which 11 genes had significant
differences in expression levels between the non-responder group and
the responder group (p < 0.05). Six genes (IL12RB2, IL32, MMP9, NGF,
OASL, and TNFRSF12A) were down-regulated and 5 genes (BMP4, CHP2,
OSGIN1, PAK5, and POMC) were up-regulated. The distribution of immune
genes expression was presented in [124]Figures 9 and [125]10. Further,
we analyzed the expression differences of 5 ir-lncRNAs from
[126]GSE45670 cohort. The distribution of ir-lncRNAs expression was
presented in [127]Figure 11. Among the above five ir-lncRNAs,
ADAMTS9-AS2, and FAM167A-AS1 were up-regulated in the non-responder
group, while LINC01121 and MIR124-2HG were down-regulated (p < 0.05).
FIGURE 9.
[128]FIGURE 9
[129]Open in a new tab
Differential expression results of immune genes in [130]GSE45670
cohort. (A) Heatmap, (B) Volcano map.
FIGURE 10.
[131]FIGURE 10
[132]Open in a new tab
Differential expression analysis of immune genes in [133]GSE45670
cohort.
FIGURE 11.
[134]FIGURE 11
[135]Open in a new tab
Differential expression analysis of five ir-lncRNAs in [136]GSE45670
cohort.
Functional gene ontology and Kyoto encyclopedia of genes and genomes
enrichment analysis of immune genes
The GO and KEGG enrichment analysis results are shown in [137]Table 4
and [138]Figure 12. In the GO molecular function enrichment analysis,
the differentially expressed immune genes were mainly enriched in
receptor ligand activity, signaling receptor activator activity, growth
factor activity, and other pathways. In biological process enrichment
analysis, differential immune genes were mainly enriched in extrinsic
apoptotic signaling pathway, regulation of apoptotic signaling pathway,
regulation of extrinsic apoptotic signaling pathway, etc. In the
cellular component enrichment analysis, the differential immune genes
were mainly enriched in endosome lumen, tertiary granule lumen, and
Golgi lumen, but the corrected p value was not statistically
significant, suggesting that these pathways were not significant. In
the KEGG enrichment analysis, the differential immune genes were mainly
enriched in Cytokine-cytokine receptor interaction, Estrogen signaling
pathway, Fluid shear stress and atherosclerosis. The immune genes
mainly involved in enrichment analysis included NGF, IL32, BMP4,
TNFRSF12A, and IL12RB2.
TABLE 4.
Main results of GO and KEGG enrichment analysis for significantly
expressed immune genes.
Type Description p p adjust Q value gene ID Count
GO_BP extrinsic apoptotic signaling pathway 0.0000 0.004322 0.0021
NGF/BMP4/TNFRSF12A/PAK5 4
GO_BP regulation of apoptotic signaling pathway 0.0000 0.014634 0.0071
MMP9/BMP4/TNFRSF12A/PAK5 4
GO_BP regulation of extrinsic apoptotic signaling pathway 0.0001
0.020688 0.0100 BMP4/TNFRSF12A/PAK5 3
GO_MF receptor ligand activity 0.0000 0.000115 0.0001
NGF/IL32/BMP4/OSGIN1/POMC 5
GO_MF signaling receptor activator activity 0.0000 0.000115 0.0001
NGF/IL32/BMP4/OSGIN1/POMC 5
GO_MF growth factor activity 0.0001 0.001442 0.0009 NGF/BMP4/OSGIN1 3
GO_CC endosome lumen 0.0195 0.211736 0.1932 NGF 1
GO_CC tertiary granule lumen 0.0305 0.211736 0.1932 MMP9 1
GO_CC Golgi lumen 0.0570 0.211736 0.1932 NGF 1
KEGG Cytokine-cytokine receptor interaction 0.0000 0.000317 0.0003
NGF/IL32/BMP4/TNFRSF12A/IL12RB2 5
KEGG Estrogen signaling pathway 0.0095 0.150833 0.1419 MMP9/POMC 2
KEGG Fluid shear stress and atherosclerosis 0.0096 0.150833 0.1419
MMP9/BMP4 2
[139]Open in a new tab
FIGURE 12.
[140]FIGURE 12
[141]Open in a new tab
Functional GO and KEGG Enrichment Analysis of Immune Genes. (A) GO
biological process enrichment, (B) GO cellular component enrichment,
(C) GO molecular function enrichment, (D) KEGG pathway enrichment.
lncRNAs isolation, cDNA synthesis, and RT-qPCR
The expression levels of ADAMTS9-AS2, FAM167A-AS1, LINC01121,
MIR124-2HG, MGC12916 and the internal reference RNA were detected by
RT-qPCR using three of the above-mentioned cell lines. Compared with
those lncRNAs in the non-responder group, expression levels of ADAMTS9-
AS2, FAM167A-AS1, and MGC12916 lncRNA in the responder group were
significantly up-regulated, while expression levels of LINC01121 and
MIR124-2HG were down-regulated. The results are shown in
[142]Supplementary Figure S4 and [143]Supplementary Table S6.
Discussion
Many recent studies have focused on establishing the signatures of
coding genes with or without non-coding RNAs to assess prognosis in
patients with malignancies ([144]Zhu et al., 2016; [145]Qu et al.,
2018; [146]Hong et al., 2020). This research strategy divides
prognostic groupings based on the absolute expression level of certain
genes of interest, of which the simplicity of operation is an important
advantage, as the prognostic model was set up on quantifying the
expression levels of transcripts ([147]Hong et al., 2020). However, the
shortcomings are also obvious, because the expression of most genes
varies greatly, especially for non-coding RNAs. This means that some
basic experiments (such as rt-qPCR) should be used for practical
verification to ensure that the genes of interest had sufficient
expression ([148]Zhu et al., 2021). To overcome the shortcomings of the
above model, new prognostic models based on the strategy of the
relative expression of immune-related gene pairing were proposed, which
became a mainstream bioinformatics research method ([149]Sun et al.,
2020; [150]Wu et al., 2020; [151]Wang Y. et al., 2021). In this study,
we were inspired by the new strategy and attempted to construct a
reasonable model with ir-lncRNA pairs assess the prognostic value of
ir-lncRNAs. Therefore, we did not apply rt-qPCR to verify their exact
expression levels in the signatures in our study. To the best of our
knowledge, there are no similar studies investigating the
radioresistance of esophageal cancer.
First, we explored the expression level of ir-lncRNAs between
non-responder ESCC patients and responder ESCC patients to investigate
the potential mechanism of radioresistance in ESCC after
chemoradiotherapy. We first screened 26 differentially expressed
ir-lncRNAs from the [152]GEO45670 dataset. The differential expression
of these lncRNAs indicates that ir-lncRNAs had played an important role
in the radioresistance of ESCC. Second, previous studies have shown
that lncRNAs achieve biological functions through a variety of target
genes, which may involve immune-related genes ([153]Wang Y. et al.,
2021). Therefore, we further explored the differential expression of
immune genes, which were not only differentially expressed in
esophageal cancer tissues with different response characteristics of
radiotherapy, but also correlated with the expression of immune
lncRNAs, and thus could be considered as target genes of ir-lncRNAs.
Third, differential co-expression analyses were performed to classify
ir-lncRNAs and immune genes, and the prognostic performance of
ir-lncRNAs were validated via single-pairing approach along with a 0 or
1 matrix by cyclical calculation. Fourth, we performed univariate
analysis combined with Lasso penalized regression to identify the
significant irlncRNA pairs. Next, multivariate regression was used to
identify ir-lncRNAs with independent effects. Fifth, the best model was
set up by calculating AUC value with ROC, where the best cut-off point
was applied to distinguish high- or low-risk groups of EPC patients.
Sixth, we further evaluated this new model in a variety of clinical
settings, including survival, clinical or pathological characters. At
last, tumor-infiltrating immune cells between different radioresistance
of esophageal cancer.
Radiotherapy is one of the main and effective treatments for advanced
EPC, but the therapeutic outcomes of which are still unsatisfactory
([154]Li et al., 2016; [155]Yan et al., 2022), because the complete
response rate of radiotherapy is less than 35%–40% ([156]Mori et al.,
2021). Radioresistance is the biggest problem and obstacle faced by
radiotherapy in esophageal cancer, as poor response to radiotherapy
could result in local failure after radiotherapy ([157]Chen H. et al.,
2020). Therefore, it is particularly important to search for
radioresistance-related molecular markers and new radiosensitizers to
enhance radio-sensitivity in ESCC cells and improve the survival of
ESCC patients with radioresistance. Radioresistance is the focus and
difficulty in cancer radiobiology research, which is also a clinically
urgent issue ([158]Wang et al., 2019; [159]Yu et al., 2020; [160]Sun et
al., 2021; [161]Liu et al., 2022). Genes are an intrinsic determinant
of tumor radioresistance, and studies have shown that multiple genes
can affect the radioresistance of esophageal cancer ([162]Huang et al.,
2019; [163]Wang et al., 2019; [164]Hua et al., 2020; [165]Yu et al.,
2020; [166]Han et al., 2021; [167]Sun et al., 2021; [168]Liu et al.,
2022). Practice has shown that it is common in clinical practice that
even patients with esophageal cancer of the same stage and the same
pathological type have great differences in the effect of radiotherapy
after receiving the same radiotherapy regimen. The fundamental reason
is the genetic differences between different individuals ([169]Huang et
al., 2019; [170]Wang et al., 2019; [171]Hua et al., 2020; [172]Yu et
al., 2020; [173]Han et al., 2021; [174]Sun et al., 2021; [175]Liu et
al., 2022). If the decisive genes of radioresistance can be screened
out, it is of great significance for the study of radiosensitization,
targeted therapy, and prediction of radiotherapy effect to guide
individualized therapy. Although more molecular studies have been
reported on the radioresistance or radiosensitivity of ESCC, the
clinical significance of most molecular markers remains unclear and
inconsistent due to the complexity and variability of detection
methods.
Recent studies had found that a variety of lncRNAs can affect the
radiosensitivity of EPC by regulating gene expression and key signal
transduction pathways ([176]Liu et al., 2021). Many lncRNAs are
involved in the regulation of the tumor immune microenvironment, which
are often referred to as immune-related lncRNAs ([177]Huang et al.,
2018). However, there is few reports on the relationship between
ir-lncRNAs and radioresistance and prognosis of ESCC. Furthermore, the
mechanism of ir-lncRNAs involved in the radioresistance is still
unclear. In recent years, the use of bioinformatics methods for data
mining at the molecular level provides new ideas for the study of
molecular pathogenesis of various diseases including tumors ([178]Zheng
et al., 2020). Our present study used bioinformatics methods to
re-analyze the radioresistance-related microarray data of esophageal
cancer from GEO, and screened for differentially expressed ir-lncRNA
genes. Through biological process annotation and signal pathway
enrichment analysis, we have excavated some ir-lncRNA genes and
immune-genes and potential signal pathways related to radioresistance
to investigate the mechanism of radioresistance at the molecular level
in ESCC. [179]Zhu et al. (2021) had identified that a total of 111
immune-related lncRNAs were different expressed in ESCC, in which,14
lncRNAs markedly related to prognosis of ESCC were identified via
univariate analysis. Finally, a model based on the 8-lncRNA signature
was identified in the multiple regression model, which demonstrated
that the 8-lncRNA signature has certain power in predicting the
prognosis of ESCC patients ([180]Zhu et al., 2021). Since the
expression abundance of ir-lncRNAs is extremely low, which is different
from the expression of protein genes, we used the relative expression
status between different ir-lncRNAs to explore the biological functions
of lncRNAs. Based on 26 differentially expressed ir-lncRNAs, we
constructed 325 ir-lncRNA pairs in the [181]GSE45670 dataset, and 180
of which were extracted to further study. As seen in our study, we
found that there were differences in gene expression between
non-responder ESCC and responder ESCC in terms of LINC01121,
FAM167A-AS1, ADAMTS9-AS2, MGC12916, MIR124-2HG. ADAMTS9-AS2 and
FAM167A-AS1 were up-regulated in the non-responder group, while
LINC01121 and MIR124-2HG were down-regulated (p < 0.05). Among the five
lncRNAs we identified; some studies have demonstrated their prognostic
effects on malignant tumors. [182]Shen et al. (2020) had found that a
worse 5-year overall survival was detected in ESCC-patients with
low-expressed ADAMTS9-AS2. Similar results were seen in clear cell
renal cell carcinoma and bladder cancer ([183]Song et al., 2019;
[184]Zhang et al., 2020), which indicates that ADAMTS9-AS2 is a tumor
suppressor gene. However, the opposite results were found in patients
with tongue squamous cell carcinoma, where high-expression of LncRNA
ADAMTS9-AS2 promotes proliferation, migration and
epithelial-mesenchymal transition (EMT) with poor prognosis, and
low-expression was detected in patient with lymph node metastasis
([185]Li Y. et al., 2019). In our study, the expression of ADAMTS9-AS2
was increased in non-responsive EPC patients, suggesting that
ADAMTS9-AS2 was involved in the radioresistance of EPC, but the
specific mechanism has not been reported. Some scholars have found that
ADAMTS9-AS2 have been shown to play essential roles in temozolomide
(TMZ) resistance in glioblastoma (GBM) ([186]Yan et al., 2019). It
should be mentioned that, although multiple lncRNAs, such as HOTAIR
([187]Lv et al., 2013), CCAT2 ([188]Zhang et al., 2015) and MALAT1
([189]Deng et al., 2016), have shown potential prognostic value in
ESCC, the role of immune-related lncRNA signatures in prognosis has not
been elucidated in the literature.
A few studies had addressed the role of immune-related lncRNAs in
survival and prognosis of ESCC. A previous study had established an
immune gene-based prognostic model for ESCC and esophageal
adenocarcinoma (EAC) ([190]Fei et al., 2021). Prognosis-related
immune-gene-based model based on BMP1, EGFR, S100A12, HLA-B, TNFSF18,
IL1B, and MAPT had proved to be useful for prognosis in ESCC ([191]Fei
et al., 2021). To verify the relationship between ir-lncRNA pairs and
survival from ESCC, we adopted the [192]GES53625 dataset as validation
cohort. We constructed a novel prognostic prediction model consisting
of three ir-lncRNA pairs, including LINC01121|FAM167A-AS1,
ADAMTS9-AS2|MGC12916, and MIR124-2HG|FAM167A-AS1, which involved five
ir-lncRNAs. We further confirmed by Wilcoxon rank sum test that there
were significant differences in the expression levels of the above
ir-lncRNAs in different risk score groups in the [193]GES53625 cohort.
In addition, the prognostic value of the novel model was confirmed by
the multivariate Cox analysis. The overall survival in the high-risk
group was significantly worse than that in the low-risk group (p <
0.001). The 1-year, 2-year, and 3-year prediction performance of this
risk-model was 0.666, 0.702, and 0.686, respectively. Univariate and
multivariate Cox analysis confirmed that prognostic riskscores based on
three ir-lncRNA pairs were important prognostic factors for ESCC.
Univariate and multivariate COX regression analysis indicated that the
three ir-lncRNA pairs were associated with the prognosis of ESCC, and
in the high-risk group ADAMTS9-AS2, LINC01121, MGC12916, and MIR124-2HG
was up-regulated, FAM167A-AS1 was up-regulated. All the three ir-lncRNA
pairs were markers of poor prognosis with worse overall survival for
ESCC. A previous research showed that ADAMTS9-AS2 is a prognostic
biomarker correlated with immune infiltrates and predicted a poorer
overall survival when it was low expressed in lung adenocarcinoma
([194]Lin et al., 2021). [195]Li W. et al. (2020) had identified the
expression levels of eight lncRNAs to establish a signature for
predicting the survival of patients with ESCC. [196]Li Z. Y. et al.
(2020) found that the expression of lncRNA Rpph1 in patients with EPC
was significantly higher than that in healthy participants (p < 0.05),
and was positively correlated with cancer tissues (r = 0.681, p <
0.05). In vitro experiments confirmed that silencing lncRNARpph1 could
up-regulate radio-induced pro-apoptotic-related proteins such as Bax,
down-regulate anti-apoptotic-related proteins such as Bcl-2, thereby
increasing radiotherapy-induced apoptosis of EPC cells. In addition,
silencing lncRNA Rpph1 can also improve the radiosensitivity of EPC
cells by reducing radiation-induced G2/M phase arrest and
epithelial-mesenchymal transition (EMT) ([197]Li Z. Y. et al., 2020).
Wang et al. found that lncRNA CCAT2 was highly expressed in EPC cells,
which can negatively regulate the expression of miR-145 and inhibit the
phosphorylation of Akt, ERK, and p70 s6K1 to increase radioresistance.
In vitro experiments confirmed that knockout of lncRNA CCAT2 can
significantly increase radiation-induced apoptosis, thereby increasing
radiosensitivity in EPC ([198]Wang et al., 2020). In addition to
mediating radioresistance, some lncRNAs also have radiosensitizing
effects. Lin et al. reported that compared with the radioresistant ESCC
cell line TE-1-R, the expression of lncRNA GAS5, and RECK was higher in
the radiosensitive cell line TE-1, while the expression of miR-21 was
lower in cell line TE-1 ([199]Lin et al., 2020). Further research found
that up-regulation of lncRNA GAS5 can increase the expression of RECK
by inhibiting miR-21, reduce the viability and colony formation ability
of EPC cells under radiation exposure, and increase radiation-induced
apoptosis of cancer cells ([200]Lin et al., 2020). All these evidences
had suggested that lncRNAs play an important role in radioresistance or
radiosensitivity of ESCC. In general, high-abundance lncRNAs have
significant biological functions, especially lncRNAs with significant
differences in expression ([201]Yan et al., 2021). Different from the
above studies, we focused on the ir-lncRNAs related to the immune
environment or immune regulation of radioresistance. As expected, we
found that ir-lncRNAs are involved in the bidirectional regulation of
radiosensitivity, that is, some lncRNAs could promote radioresistance,
while others may increase radiosensitivity.
In recent years, the study of the tumor immune microenvironment has
taken a leading role in field of cancer research ([202]Wang et al.,
2022). The differential expression of immune genes is an important
molecular mechanism leading to changes in the tumor immune
microenvironment ([203]Anderson, 2020; [204]Peña-Romero and
Orenes-Piñero, 2022). Several previous studies reported the prognostic
value of a single immune-related gene in EPC or lung cancer, such as
FGFR1, TNFRSF10B, and IL1B ([205]Schabath et al., 2013; [206]Takase et
al., 2016; [207]Zhang Y. et al., 2021). However, the study focused on
the prognostic role of the ir-lncRNAs with microenvironment and immune
cells in ESCC is lacking, especially for ESCC with radioresistance.
Identification of prognostic value of the tumor microenvironment in
esophageal cancer is necessary ([208]Tan et al., 2021). In our study,
we identified 11 immune genes with differential expression in the GSE
45670 cohort, and our findings further confirmed the involvement of
immune alterations in the biological process of radiation resistance in
ESCC. To further explore the changes and biological pathways in the
immune microenvironment between patients in the non-responder group and
patients in the responder group, CIBERSORT was carried out in this
study as well, which was an analytical tool developed by [209]Newman et
al. (2015) to provide an estimation of the abundances of member cell
types in a mixed cell population via gene expression data, and provided
the possibility of identifying immune biomarkers for diagnosis and
prognosis. The GSE 45670 dataset was utilized to evaluate immunity
infiltration via the CIBERSORT algorithm. Growing evidence pointed to
the underlying mechanisms by which the local immune microenvironment
and immune cells drive tumorigenesis in many cancers ([210]Greten and
Grivennikov, 2019; [211]Wang Q. et al., 2021; [212]Kumar et al., 2022).
Previous studies have shown that DNA damage repair-related genes,
apoptosis-related genes, cellular hypoxia-related genes, cell
cycle-related genes, and autophagy genes play important roles in
radiosensitivity by changing the microenvironment ([213]Chen et al.,
2017; [214]Tang et al., 2018; [215]Zhang H. et al., 2021).
Tumor-infiltrating immune cells (TIICs) in esophageal cancer tissue may
be an important determinant of prognosis and therapy response ([216]Lu
et al., 2020). In our study, we found that there were some differences
in infiltrating immune cells between the non-responder group and the
responder group, which mainly manifested in differences in infiltration
of activated mast cell and macrophage. Mast cells are an important
member of innate immune cells. Circulating mast cells contribute to the
growth and metastasis of many tumors, while mast cell infiltration in
tumor tissue is closely related to tumor survival in some cancer
([217]Welsh et al., 2005). In a study of lung cancer, it has been shown
that more human mast cells infiltrated in cancer tissue improved the
survival of cancer, suggesting that mast cells can participate in the
antitumor immune process ([218]Welsh et al., 2005). In ESCC cases with
tumor complete remission, the number of mast cells infiltrated was
higher, indicating that mast cells were involved in the therapeutic
effect of radiotherapy.
There are a large number of tumor-associated macrophages (TAMs) in the
tumor microenvironment, which have a high degree of interaction with
tumor cells, tumor stem cells, epidermal cells, fibroblasts, T/B cells,
and NK cells ([219]Welsh et al., 2005). Although macrophages
theoretically have the ability to destroy tumors, there is growing
experimental evidence that TAMs promote tumor progression
([220]Mantovani and Allavena, 2015). Resilience and diversity are two
characteristics of macrophages, which means that macrophages have dual
effects in tumor development ([221]Cassetta and Pollard, 2020).
According to the activation type of macrophages and their different
roles in the tumor microenvironment, TAMs are generally classified into
two functionally opposite subtypes, classically activated M1
macrophages and alternately activated M2 macrophages, both of which
represent one of the main tumor-infiltrating immune cell types
([222]Cassetta and Pollard, 2020). Basic research has found that a
large number of TAMs proliferate after radiation, and at the same time
release a large number of inflammatory signals (IL-1) and
immunosuppressive signals (TGF-b). Unfortunately, the massive
accumulation of macrophages can lead to tumor recurrence, which is very
similar to TAMs-guided tissue damage repair after chemotherapy. M1
macrophages have antitumor effects, while M2 macrophages mainly play a
role in promoting tumor growth, invasion and metastasis. Our study
found that macrophages infiltrated more in esophageal cancer tissues in
responder group, among which M0 macrophages were the main infiltrating
cells. The infiltrating number of M2 macrophages was more often in the
non-responder group, although there was no statistical difference
compared with that in responder group. It is worth noting that both M1
and M2 macrophages have high degree of plasticity, which makes it
possible to design appropriate methods to re-induce and re-educate
them, thereby becoming an effective weapon against tumors. Different
types of macrophages can be converted into each other upon tumor
microenvironment changes or therapeutic interventions. In view of the
important role of macrophages in tumor radioresistance, it needs to do
in-depth research on the realization of their functions and their
regulatory mechanisms, in order to find new anti-tumor targets.
Some limitations should be mentioned in our study. First, the
prognostic model based on ir-lncRNA pairs was established through
bioinformatics analyses from data available in the GEO databases.
Hence, some further prospective trials or experimental data should be
performed to validate the findings of this study. Second, our study
found that most of clinical characters were not correlated with the
riskscore, and tumor grade had a weak correlation with the riskscore.
We speculate that preoperative adjuvant therapy altered gene expression
status, resulting in clinical factors not correlated with gene
expression. In addition, insufficient sample size may also be an
important reason. Third, our study only initially explored the
potential relationship between ir-lncRNAs risk signatures and immune
cell infiltration, so further studies are needed to reveal the
underlying mechanisms. Fourth, the effect of ir-lncRNAs on the immune
microenvironment was indirectly speculated through immune genes, and
there is currently a lack of direct evidence. More experiments are
needed to confirm the impact of ir-lncRNAs on the
radioresistance-related microenvironment of ESCC. Finally, as a
preliminary exploratory study, it only has a qualitative role in the
identification of ir-lncRNAs and the prognosis risk of ESCC patients.
The relationship between ir-lncRNAs and the prognosis of ESCC patients
has not been accurately quantified, which still needs to be further
verified by multicenter studies with large samples.
In conclusion, ir-lncRNAs may be involved in the biological regulation
of radioresistance in patients with ESCC. Changes in macrophage
infiltration and immune gene expression are potential mechanisms of
radiotherapy resistance, which are worthy of further study. This study
had successfully established a prognostic risk model based on three
ir-lncRNAs pairs, which is an important attempt to identify and predict
the prognosis of ESCC. More importantly, it is a useful supplementary
method to predict the prognosis of patients with esophageal squamous
cell carcinoma based on TNM staging.
Data availability statement
The original contributions presented in the study are included in the
article/[223]Supplementary Material, further inquiries can be directed
to the corresponding author.
Author contributions
Conception and design: JZ and XC; Provision of study materials or
patients: JZ and BH; Collection and assembly of data: JZ and BH; Data
analysis and interpretation: JZ and JL; Manuscript writing and editing:
JZ, XC, BH, and JL. All authors contributed to the article and approved
the submitted version.
Funding
This study was supported in part by the National Clinical Key Specialty
Construction Program (Grant No. 2021), and the Fujian Provincial
Clinical Research Center for Cancer Radiotherapy and Immunotherapy
(Grant No. 2020Y2012), Training and Nurturing Project for Young and
Middle-aged Leading Talents in Healthcare in Fujian Province (Fujian
Healthcare Personnel Document 2022-954 to XC), Joint Funds for the
Innovation of Science and Technology, Fujian Province (Grant No:
2020Y9036 to XC), Fujian provincial health technology project (Youth
Scientific Research Project, 2019-1-50 to JZ) and the Nursery Fund
Project of the Second Affiliated Hospital of Fujian Medical University
(Grant No: 2021MP05 to JZ).
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:
[224]https://www.frontiersin.org/articles/10.3389/fgene.2022.921902/ful
l#supplementary-material
Supplementary Figure S1
Selection of the optimal candidate genes in the LASSO model.
Supplementary Figure S2
LASSO coefficients of prognosis-associated ir-lncRNA pairs, each curve
represents a gene.
Supplementary Figure S3
Differential distribution of risk score with different clinical
characters.
Supplementary Figure S4
The expression levels of ir-lncRNAs detected by RT-qPCR, *p < 0.05, **p
< 0.01, ***p < 0.001.
[225]Click here for additional data file.^ (1.3MB, JPEG)
[226]Click here for additional data file.^ (104.4KB, JPEG)
[227]Click here for additional data file.^ (144.3KB, JPEG)
[228]Click here for additional data file.^ (988.3KB, TIFF)
[229]Click here for additional data file.^ (81.8KB, xlsx)
References