Abstract

   Introduction: Molecular chaperones and long non-coding RNAs (lncRNAs)
   have been confirmed to be closely related to the occurrence and
   development of tumors, especially lung cancer. Our study aimed to
   construct a kind of molecular chaperone-related long non-coding RNAs
   (MCRLncs) marker to accurately predict the prognosis of lung
   adenocarcinoma (LUAD) patients and find new immunotherapy targets.

   Methods: In this study, we acquired molecular chaperone genes from two
   databases, Genecards and molecular signatures database (MsigDB). And
   then, we downloaded transcriptome data, clinical data, and mutation
   information of LUAD patients through the Cancer Genome Atlas (TCGA).
   MCRLncs were determined by Spearman correlation analysis. We used
   univariate, least absolute shrinkage and selection operator (LASSO) and
   multivariate Cox regression analysis to construct risk models.
   Kaplan-meier (KM) analysis was used to understand the difference in
   survival between high and low-risk groups. Nomogram, calibration curve,
   concordance index (C-index) curve, and receiver operating
   characteristic (ROC) curve were used to evaluate the accuracy of the
   risk model prediction. In addition, we used gene ontology (GO)
   enrichment analysis and kyoto encyclopedia of genes and genomes (KEGG)
   enrichment analyses to explore the potential biological functions of
   MCRLncs. Immune microenvironmental landscapes were constructed by using
   single-sample gene set enrichment analysis (ssGSEA), tumor immune
   dysfunction and exclusion (TIDE) algorithm, “pRRophetic” R package, and
   “IMvigor210” dataset. The stem cell index based on mRNAsi expression
   was used to further evaluate the patient’s prognosis.

   Results: Sixteen MCRLncs were identified as independent prognostic
   indicators in patients with LUAD. Patients in the high-risk group had
   significantly worse overall survival (OS). ROC curve suggested that the
   prognostic features of MCRLncs had a good predictive ability for OS.
   Immune system activation was more pronounced in the high-risk group.
   Prognostic features of the high-risk group were strongly associated
   with exclusion and cancer-associated fibroblasts (CAF). According to
   this prognostic model, a total of 15 potential chemotherapeutic agents
   were screened for the treatment of LUAD. Immunotherapy analysis showed
   that the selected chemotherapeutic drugs had potential application
   value. Stem cell index mRNAsi correlates with prognosis in patients
   with LUAD.

   Conclusion: Our study established a kind of novel MCRLncs marker that
   can effectively predict OS in LUAD patients and provided a new model
   for the application of immunotherapy in clinical practice.

   Keywords: TCGA, LUAD, lncRNA, molecular chaperone-related lncRNA index,
   prognosis, immunotherapy

Introduction

   Lung cancer is a malignant tumor originating from the bronchi and
   alveoli. Worldwide, 1.77 million lung cancer deaths occur each year,
   and it is the leading cause of cancer death in the world ([36]Siegel et
   al., 2021). LUAD is a subtype of lung cancer and a highly heterogeneous
   malignancy, accounting for approximately half of all lung cancers
   ([37]Sivakumar et al., 2017; [38]Xu et al., 2020a). Studies have shown
   that the risk factors for LUAD mainly come from direct exposure to
   tobacco. LUAD tends to occur early in East Asian women who do not smoke
   ([39]Chen et al., 2020; [40]Devarakonda et al., 2021). This is related
   to the presence of Epidermal growth factor receptor (EGFR) mutations in
   East Asian LUAD patients ([41]Choong and Sung, 2021; [42]He et al.,
   2021). Lung cancer usually involves pleura and has a poor prognosis.
   The 5-year survival rate is less than 20% ([43]Zhang et al., 2021).
   LUAD is prone to distant metastasis, and the common sites of metastasis
   are brain, liver, bone, adrenal gland and pleura ([44]Klikovits et al.,
   2018). In the past few decades, treatment of LUAD has mainly included
   surgery, chemotherapy, and emerging immunotherapies. Although recent
   advances in LUAD have greatly improved the prognosis of LUAD patients,
   the OS of advanced LUAD patients is still very low. Therefore,
   developing new biomarkers to predict the prognosis of LUAD patients and
   find potential therapeutic targets for LUAD is crucial.

   Molecular chaperones are molecular assistants that assist in the
   folding and assembly of intracellular proteins and play an important
   role in intracellular life activities ([45]Shan et al., 2020; [46]Wei
   et al., 2022). Studies have shown that molecular chaperones play an
   important role in the occurrence and development of tumors, and it has
   also been confirmed as a prognostic marker for tumors ([47]Zhu et al.,
   2013; [48]Lu et al., 2020). Jia et al. ([49]Jia et al., 2021) found
   that heat shock protein 90 (HSP90) can promote the metastasis of LUAD
   cells by interacting with the oncogene EEF1A2. It ultimately leads to
   the poor prognosis of LUAD patients. In addition, drugs targeting
   molecular chaperone-related genes are also widely used in clinical
   practice. chaperones-related related genes with great clinical
   potential, such as HSP90 and P53, making immunotherapy become an
   important means of tumor treatment ([50]Armstrong et al., 2018;
   [51]Kaida and Iwakuma, 2021). Although molecular chaperones can be used
   as prognostic indicators in patients with LUAD, the utility of using
   only a single biomarker is limited. Therefore, establishing reliable
   biomarkers for the construction of LUAD prognostic models is an urgent
   clinical task.

   LncRNAs refer to a class of RNAs longer than 200 nucleotides that
   cannot encode complete proteins ([52]Vollmers and Carpenter, 2022).
   Studies have shown that lncRNAs play a crucial role in the occurrence
   and development of tumors, including LUAD ([53]Dong et al., 2018).
   Meanwhile, studies have reported that immune infiltrating cells play a
   crucial role in the progression and invasion of tumor cells ([54]Li et
   al., 2020a). LncRNA is a key regulatory element in the immune system,
   which has the functions of antigen presentation, antigen release,
   immune migration, immune infiltration, and immune activation ([55]Xia
   et al., 2020; [56]Zhang et al., 2020). With the deepening of research,
   the role of lncRNAs as ideal diagnostic markers for tumors has been
   gradually discovered ([57]Peng et al., 2016). These all suggest that
   lncRNAs may be used as a new biomarker to improve the prognosis and
   treatment of LUAD patients. Zhou et al. ([58]Zhou et al., 2022) found
   that lncRNAs can promote the occurrence and development of LUAD by
   binding to HSP90. However, the current research on the pathogenesis of
   MCRLncs in LUAD is still lacking. Therefore, we attempted to use
   transcriptome data and clinical data from the TCGA database to develop
   a marker of MCRLncs with guiding significance for immunotherapy. This
   study provides a new model with prognostic value for LUAD patients and
   establishes a new feature to predict LUAD patients’ response to
   immunotherapy.

Materials and methods

The data source

   The Genecards is a comprehensive bioinformatics database from 125 web
   sources such as NCBI and UCSC, which provides detailed information on
   all currently annotated and predictable genes. It covers information
   including genome, proteome, transcription, and function ([59]Safran et
   al., 2010; [60]Stelzer et al., 2011). Through the Genecards database
   portal ([61]https://www.genecards.org/), we obtained 233 genes related
   to chaperones (Relevance score ≥ 5). MSigDB is a treasure trove
   database for analyzing gene enrichment pathways ([62]Liberzon et al.,
   2011). We downloaded 17 GSEA functional pathways from MsigDB
   ([63]Supplementary Table S1). These pathways enriched 517 genes in
   total, and we removed duplicated genes to get 312 genes. Finally, we
   pooled the 233 genes obtained based on the Genecards database and the
   312 genes obtained based on the MsigDB, removed the duplicate genes,
   and finally obtained 417 molecular-chaperone genes.

   TCGA is a landmark Human Cancer Genome Project that has collected
   molecular numbers from more than 20,000 primary cancer samples,
   including LUAD ([64]Cai et al., 2021a). In this study, we from the TCGA
   data portal website ([65]https://portal.gdc.cancer.gov/) downloaded
   transcriptome data, clinical data, and mutation data in LUAD. Clinical
   data included OS, survival status, age, sex, tumor grade, and tumor
   node metastasis (TNM) stage. In this study, the samples were randomly
   divided into a training subgroup (n = 330) and a validation subgroup (n
   = 164). The training subgroup was used to construct a molecular
   chaperone-related lncRNAs risk model. The validation subgroup and the
   whole group were used to validate the risk model. In addition, the
   IMvigor210 dataset, a cohort of atezolizumab (anti-PD-L1 monoclonal
   antibody) for the treatment of bladder cancer, was extracted to
   evaluate the response of MCRLncs markers to immunotherapy efficacy
   ([66]Powles et al., 2014).

Identification of molecular chaperone-related lncRNAs in patients with LUAD

   We used the Spearman correlation analysis (| cor | > 0.4 and p < 0.001)
   to identify MCRLncs. The heatmap showed the expression correlation
   between molecular chaperone genes and lncRNAs.

Molecular chaperones-related lncRNAs risk models were constructed

   We divided all patients with lung adenocarcinoma into a training
   subgroup (n = 330) and a validation subgroup (n = 164). Data from the
   training subgroup were used to construct a prognostic model. We
   screened 301 MCRLncs by using univariate Cox regression analysis (p
   values less than 0.05 were considered significant). To narrow down the
   independent variables and avoid overfitting prognostic features, we
   performed the LASSO regression analysis on these MCRLncs ([67]Meier et
   al., 2008). Next, we performed multivariate Cox regression analysis on
   the MCRLncs obtained by LASSO regression analysis, and finally screened
   16 MCRLncs as candidate genes. The 16 candidate genes were used to
   construct prognostic models. To construct risk characteristics and
   calculate risk scores, the coefficients and expressive values of
   MCRLncs screened out from LASSO regression were used to calculate an
   individual risk score. The risk score represents a prognostic feature
   of chaperone-related lncRNAs, which helps us to distinguish high-risk
   LUAD patients from low-risk LUAD patients. Our risk score calculation
   formula is as follows:

   Coef[i] was the coefficient of lncRNA in LASSO regression. Coefficient
   xi is the expression value of selected MCRLncs ([68]Li et al., 2020b;
   [69]Li et al., 2020c; [70]Liu et al., 2020; [71]Liang et al., 2021).
   This formula is used to calculate the risk score.

   To build a more intuitive model, we divided LUAD patients into
   high-risk and low-risk groups using the median risk score as a cutoff
   point. Next, we plotted KM survival curves, risk score distribution
   maps, and heatmaps to identify differences in the expression of MCRLncs
   between high- and low-risk groups.

Validation of risk prognostic models

   We plotted the ROC curve and C-index curve to verify the prognostic
   value of our constructed prognostic model. Based on the results of
   multivariate Cox regression analysis, we constructed a nomogram that
   can predict the occurrence of 1-, 3-, and 5-year OS survival in
   patients with LUAD. Nomogram is widely used for graphical calculations
   of complex formulas with practical accuracy. We can obtain the score of
   each clinical feature from the nomogram and predict the 1-year, 3-year,
   and 5-year survival rates of LUAD patients through the total score.
   Next, we evaluated the performance of the nomogram by drawing a
   calibration curve.

Correlation analysis of risk score and clinicopathological features

   To further validate the accuracy and specificity of the prognostic
   model, we used univariate and multivariate Cox analyses to screen
   variables and further explore independent risk factors associated with
   LUAD prognosis. We mapped two forests based on independent prognostic
   analyses to determine whether the prognostic model can be used as an
   independent prognostic indicator without reference to other clinical
   characteristics, including age, sex, race, tumor grade, primary tumor
   (T), regional lymph nodes (N), distant metastasis (M), and risk score.

Principal component analysis was used to assess high-risk and low-risk
patients

   To assess whether LUAD patients were discriminative between high and
   low-risk groups, we visualized gene expression profiles using
   dimensionality reduction techniques. The expression of coding genes,
   lncRNA genes, all genes and high-risk lncRNAs genes in the risk model
   was analyzed by PCA analysis.

Go functional pathway and KEGG enrichment analysis of differential MCRLncs in
LUAD

   To understand the differential expression of MCRLncs between high and
   low-risk groups, we used the “limma” R package to perform differential
   expression analysis of lncRNAs in high and low-risk groups, and
   extracted differentially expressed genes (DEGs) for the risk model.
   Next, we performed GO functional annotation and KEGG pathway enrichment
   analysis on the significant DEGs using the “clusterProfiler” R package,
   and the false discovery rate (FDR) of 0.05 was considered statistically
   significant ([72]Mazandu et al., 2017; [73]Kanehisa et al., 2019;
   [74]Zou et al., 2021; [75]Liang et al., 2022). GO database standardized
   description of gene products from three levels of biological process
   (BP), cellular component (CC), and molecular function (MF). Through GO
   functional annotation, we can understand the biological functions and
   pathways of DEGs enrichment. KEGG enrichment analysis can know in which
   pathways DEGs are enriched.

Immune infiltration analysis was based on the single-sample gene set
enrichment analysis

   In order to study the immune infiltration of individual samples from
   the high-risk group and the low-risk group, we first downloaded the
   expression levels of specific marker genes under 13 immune function
   pathways. Next, we calculated the enrichment score for 13 immune
   function pathways in each LUAD sample by using ssGSEA ([76]Barbie et
   al., 2009). We used heatmaps to display the immune infiltration of
   prognostic lncRNAs in high-risk and low-risk groups.

Calculation of tumor mutational burden

   High tumor mutational burden (TMB) is defined as the total number of
   somatic gene coding errors, base substitutions, gene insertions or
   deletions detected per megabyte ([77]Tang et al., 2019; [78]Li et al.,
   2020d). Recent studies have shown that TMB is associated with OS after
   immunotherapy in multiple cancer types, and suggest that TMB can be
   used as a predictive biomarker of immune checkpoint inhibitor treatment
   response ([79]Xu et al., 2020b; [80]Lin et al., 2020; [81]Tan et al.,
   2020). As suggested by the reviewer, we cited references to supplement