Construction of Immune-related lncRNA Prognosis Models and Development Related Networks and Drugs in LGG

doi:10.21203/rs.3.rs-757554/v1

Download PDF

Primary research

Construction of Immune-related lncRNA Prognosis Models and Development Related Networks and Drugs in LGG

https://doi.org/10.21203/rs.3.rs-757554/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Glioma is the most important tumor of the nervous system. LncRNA plays an important role in the occurrence, development, and metastasis of glioma. The immune status of glioma plays an important role in its treatment and prognosis. In lower-grade glioma (LGG), based on the results of the Weighted correlation network analysis (WCGNA) of lncRNA, we used lasso regression, random forest model, and stepwise regression to establish the lncRNA prognostic model: 0.58349*CYTOR + 0. 47804*LINC01831 + 0. 24933*HOTAIRM1 + 0. 73600*AC022034 + 0. 62351*AC104407 = lncRNA model prognosis score (LMPS). Concomitantly, we explored the correlation of LMPS with the Estimation of STromal and Immune cells in MAlignant Tumours using Expression data (ESTIMATE), Tumor IMmune Estimation Resource (TIMER) infiltrating immune cells. We also established the reciprocal network of the model and verified the significant prognostic significance of the hub genes of the network and the significant impact on immunity. Finally, by connectivity Map (cMap) analysis, we obtained the drugs about the WCGNA module where the model is. Our study highlights the close role of the lncRNA model with patient prognosis, immunity, and clinical staging, and can guide immunotherapy. The network is constructed with the model as the core. Its core genes also have a guiding role in patient prognosis, immunity, and treatment. We identified loratadine and fidaxomicin as possible drugs.

Cancer Biology

WCGNA

LGG

prognosis

immune cells

ESTIMATE

clinical

network

cMap

Glioma is the most common malignant primary tumor of the central nervous system (CNS)^[1]. Gliomas contain two main subgroups: diffuse gliomas and gliomas that show a more circumscribed growth pattern ('non-diffuse gliomas'). Based on the presence/absence of IDH mutations and 1p/19q codes, diffuse gliomas are classified into astrocytomas and oligodendrogliomas. Non-diffuse gliomas are divided into ventricular distal tumors, mixed gliomas, and other gliomas. WHO grade I and II gliomas are low-grade gliomas (LGGs). WHO grade III glioma is anaplastic glioma, but at TCGA it is contained in the LGG. WHO grade IV glioma is glioblastoma^[2].

LncRNAs include a variety of RNAs larger than 200 nucleotides (nt) in size that are not protein-translatable^[3]. Their structures are similar to mRNAs, but usually have a 5'm7G cap and a 3′polyA tail. And mRNAs are more tissue-specific and dynamic, predicting that these molecules have different biological roles^[4]. In the past few years, significant advances have been not only in the identifying of lncRNAs but also in the understanding of the regulatory aspects of gene expression, protein translation, and stability^[5]. In addition, lncRNAs play a vital role in metabolism and DNA repair^[6]. The role of lncRNA in tumor development has been a hot topic this year. Different studies have previously reported the ability of lncRNAs in various aspects such as proliferation^[7], growth inhibition^[8], motility^[9], viability ^{[10, 11]}, and angiogenesis^[10] in glioma. Thus, lncRNA can be used as a prognostic factor to determine the prognosis of patients.

The communication between cells and their microenvironment is vital for maintaining normal or cancerous tissue growth. In cancer tissues, the exchange of substances and signal transduction between cancer cells and their stromal cells significantly affects the development and metastasis of cancer and even the prognosis of cancer patients^[12]. Vascular, fibroblasts, M2 macrophages, and immunosuppressive cells are the main components that support tumor cell growth in the tumor microenvironment. On the contrary, dendritic cells, NK cells, and CD8T cells are mainly anti-tumor cells^[13]. Various lncRNAs, including MALAT1^[14], LINC00152^[15], and BCYRN1^[16], have been reported in the literature to have a significant impact on patient prognosis and may serve as new therapeutic targets.

The brain is a dynamic immune organ. At the baseline level, the brain is at immune rest. The immune situation of the central nervous system is different from the periphery, with the blood-brain barrier as the first line of immune defense, allowing only a few small hydrophobic molecules to pass and excluding the majority of bacterial invasion ^{[17, 18]}. Most peripheral immune cells do not enter the brain in the quiescent state. Myeloid cell differentiated microglia are the resident immune cells in the brain. In the inflammatory status, chemokines induced by interferon direct immune cells into the brain via pathways such as cervical lymphatic^[19] to participate in the immune response. Primary tumors in the brain shape their tumor microenvironment in a different way than metastases. Their tumor immune microenvironment also predicts to some extent the prognosis of tumor immunotherapy^[20].

2.1 Data processing

We downloaded GBM data, LGG, and its clinical data from the UCSC-Xena (https://xena.ucsc.edu/) and got the CGGA data from the GDC portal (https://portal.gdc.cancer.gov/). Gene expression data are in Fragments Per Kilobase Million (FPKM) format. We used the R software (Version 4.0.3; https://www.R-project.org) to process all data. We screened 523 LGG gene expression data samples with matching clinical data for the next stage of analysis.

2.2 Weighted Correlation Network Analysis (WCGNA)

We used the WCGNA package to screen for survival-associated immune gene modules. We then performed correlation analysis with the screened survival-related immunogene modules with lncRNA and miRNA data. We obtained immune-related lncRNAs with r = 0.3 and p = 0.001 as the cut point. LncRNAs were finally subjected to WCGNA analysis to obtain survival-related immune-related lncRNAs.

2.3 Model Construction

We repeatedly constructed the regression model 200 times and counted the lncRNAs for the model (counts > 40). Then we used a random forest model to calculate the model weights for each gene (value > 0.001). Further, we used a univariate regression model to obtain the meaningful lncRNAs for constructing the model (p < 0.05). Finally, lncRNAs of the intersection of the three methods were used for model construction by stepwise regression (p < 0.05).

2.4 Differential expression of LMSP and constitutive model genes in different clinical phenotypes

To verify the clinical value of the model and to guide the clinical use of the model, we studied the expression of lncRNA model prognosis score (LMPS) and different genes in different ages, pathological types, and clinical diagnoses. The ‘limma’ package and the ’ggpubr’ package were respectively used for data processing and plotting.

2.5 Relationship between TME and LMPS

We used the R 'estimate' package to perform the Estimation of STromal and Immune cells in MAlignant Tumours using Expression data (ESTIMATE) algorithm to estimate the patient's immune as well as stromal scores. Then we used the 'stats' package to plot the correlation between the calculated LMPS and immune and stromal scores, and the 'ggplot2' and 'ggExtra' packages to plot the graphs.

2.6 Correlation between immune cells and LMPS

We obtained the estimated immune cell infiltration expression of LGG data among Tumor IMmune Estimation Resource (TIMER; https://cistrome.shinyapps.io/timer/). We then processed the data with the 'Lima' package, the 'stats' package to obtain the Spearman correlation coefficient of the function, used 'ggplot2' to plot the scatter plot and add trend lines.

2.7 Construction of lncRNA-related action network and screening of core genes

We obtained highly correlated mRNA-miRNA pairs (r > 0.3) based on immune protein genes with a vital impact on patient survival by WGCNA. We took the same approach to obtain miRNA-lncRNA pairs (r > 0.3). We constructed mRNA-lncRNA pairs (r > 0.3) using lncRNAs and immune mRNAs screened by WCGNA that were significantly associated with survival. We inputted the aforementioned RNA pairs into Cytoscape (https://cytoscape.org/; version 3.8.2), using common genes as nodes and correlations as edges to build up a correlation network. The lncRNAs that make up the prognosis model were then used as the core to construct the lncRNA-associated action network. Finally, we used the ‘cytohub’ plugin to screen for hub genes.

2.8 survival analysis of hub genes

We explored different survival curves for different hub genes. We used the 'survival' and 'ggplot2' packages to process the data and plot survival curves.

2.9 TME analysis of hub genes

We analyzed the immune microenvironment of the hub gene and explored the relationship between the expression of the hub gene and the immune and stromal scores of the tumor. The procedure was as described above.

2.10 Gene Set Enrichment Analysis (GSEA) and co-expression of enriched pathway genes with prognostic models

We classified the samples according to the level of LMPS and then enriched the expression data for analysis. We used the 'limma', 'org. Hs. Eg. DB' and 'clusterProfiler' packages to process the data and perform enrichment analysis, and 'enrichplot ' to visualize the results. For the significantly enriched immune-related pathways, we downloaded the corresponding gene sets from the GSEA website and then performed co-expression analysis between the cancer genes in them and LMPS.

2.11 Connectivity Map (cMap) analysis

We performed a differential analysis of LGG genes by the line package. We then intersected differential genes and genes with evident relationships to survival by WCGNA. Immediately after, we performed map analysis in the L1000 platform (https://clue.io) using differentially expressed immune genes, and the obtained meaningful results were presented using 3D scatter plots. We found the drug structure from Pubchem (https://pubchem.ncbi.nlm.nih.gov/).

2.12 IDH1mutation and TP53mutation

We downloaded the relevant IDH1 mutation data and TP53 mutation data from the official website of TCGA. We then processed the data using the 'string' package and presented the results using the 'ggpubr' package.

2.13 Statistical methods

If not specified, the correlations we used were Pearson correlation coefficients, and the cutoff values we used were generally 0.05.

3.1 the WCGNA analysis of immune protein-coding genes

We construct the WCGNA model of the immune protein-coding genes in LGG. As shown in Fig. 1A, we divided the immune protein-coding genes into six modules. The grey module was the collection of unclassifiable genes. The turquoise module had the most genes, but the red module contained the least. The blue, green, and yellow modules had more concentrated gene distribution than the other modules. The correlation among modules was shown in Fig. 1B. Red and yellow modules had the strongest correlation with each other, and green modules had a strong correlation with both red and brown modules. Figure 1C demonstrated the correlation between the modules and the clinical phenotype of the patients. We could easily find that the turquoise module, the green module, and the brown color were all significantly associated with survival status and survival time in LGG patients. Considering the correlation with the patient's diagnosis, age, and the correlation status between the modules, we chose the green module as our study subject. To further clarify the effect on survival, we further analyzed the gene expression of green module genes for correlation with survival time. The results obtained are presented in the form of a scatter plot (Fig. 1D).

3.2 The WCGNA analysis of the immune-related lncRNA

We calculated the correlation between the protein genes of the green module and the lncRNA of LGG and selected the significant lncRNAs with a correlation as the target for further study. We used the lncRNAs to build up the WGCNA model.

We construct the WCGNA model of the immune protein-coding genes in LGG. As shown in Fig. 2A, we divided the immune protein-coding genes into nine modules. The turquoise module had the most genes, but the magenta module contained the least. The black, red, magenta and pink modules had more concentrated gene distribution than the other modules. The correlation among modules was shown in Fig. 2B. Pink and magenta modules had a positive correlation with each other, and the green module had a positive correlation with both blue and turquoise, while the pink module had a negative correlation with the black. Figure 2C demonstrated the correlation between the modules and the clinical phenotype of the patients. We could easily find that the pink module, the magenta module, the turquoise, and the black modules were all significantly associated with survival status and survival time in LGG patients. Considering the correlation with the patient's diagnosis, age, and the correlation status between the modules, we preliminary chose pink and black modules as our study subjects. To further clarify the effect on survival, we further analyzed the gene expression of green module genes for correlation with survival time. The results obtained are presented in the form of a scatter plot (Fig. 2D, Fig. 2E). Finally, we selected the pink module as our study object for its obvious correlation between the expression of the genes and the survival time.

3.3 The prognosis model of lncRNA in the pink module

We then built the lncRNA prognosis model to predict the survival outcome of patients. We divided the patients' survival data into training and test sets according to a 2:1 ratio. Then we used a combination of random forest and lasso regression models and univariate Cox regression to model cancer prognosis. LMPS (lncRNA model prognostic score) = 0.58349*CYTOR + 0.47804*LINC01831 + 0.24933*HOTAIRM1 + 0.73600*AC022034 + 0.62351*AC104407. In the training dataset (Fig. 3A), the ROC AUC for one year, three years, and five years were 0.852, 0.862, and 0.767, respectively. The ROC AUC of the test dataset (Fig. 3B) for one year, three years, and five years were 0.935, 0.718, and 0.643. Because of the lack of corresponding lncRNA samples, we tested our model among glioblastoma samples and obtained the AUC values of the ROC curves. The ROC AUC of CGGA (Fig. 3E) for one year and three years were 0.598 and 0.631. And The ROC AUC of TCGA-GBM (Fig. 3F) for one year and three years were 0.608 and 0.704. Significant differences were found in patient survival between the high lncRNA model predict score group (LMPS) and the low LMS group in the training set (Fig. 3C; p < 0.0001), the test set (Fig. 3D; p = 0.0053), CGGA (Fig. 3G; p = 0.039) and TCGA_GBM (Fig. 3H; p = 0.011) datasets.

3.4 The Correlation of LMPS and clinical phenotypes

We next analyzed the relationship between the clinical phenotypes of the patients and LMPS and the genes which built up the model. We were able to significantly find that LMPS (P = 0.0011), AC022034 (p = 0.0073), CYTOR(p = 0.0029), HOTAIRM1 (p = 0.012), LINC01831 (p = 0.0007), and expression was higher in patients aged 65 years or older (Fig. 4A). The expression of LMPS (p = 1.7e-7), CYTOR (p = 4.7e-12), HOTAIRM1 (p = 3.5e-8), LINC01831 (p = 0.0003) was significantly higher in patients with pathology G3 than in patients with pathology G2 (Fig. 4B). There were significant differences in the expression of LMPS, CYTOR, HOTAIRM1, and LINC01831 among the patients diagnosed with anaplastic astrocytoma, astrocytoma NOS, mixed glioma, anaplastic oligodendroglioma, and oligodendroglioma NOS (Fig. 4C). Patients with IDH1 mutation have lower LMPS than non-mutation (Supfig4A, p < 2.2e-16). And patients with TP53 mutation have lower LMPS than non-mutation as well (Supfig4B, p = 5.6e-8).

3.5 Relationship between LMPS and tumor immune microenvironment (TME)

We used the estimate algorithm to calculate the scores of the immune microenvironment of patients and obtained the relationship between LMS and immune scores and stromal scores. The LMPS and genes which made up the model had a positive correlation with immune and stromal scores. The correlation coefficients of LMPS, HOTAIRM1, LINC01831, AC104407, AC022034 and CYTOR with immune score and stroma score (Fig. 5) were 0.35 and 0.16 (p < 2.2e-16, p = 0.00023), 0.49 and 0.55 (p < 2.2e-16, p < 2.2e-16), -0.029 and − 0.082(p = 0.51, p = 0.05), 0.13 and − 0.016(p = 0.0027, p = 0.72), 0.22 and 0.022 (5.8e-7, 0.62), 0.39 and 0.22(2.2e-16,2.7e-7).

3.6 The correlation of immune cells and lncRNA prognosis model.

We can clearly fund that the predicted results by the lncRNA model had a significant positive correlation with CD8 T cells, neutrophils, B cells, dendritic cells and macrophages (R = 0.43, 0.33, 0.23, 0.26, and 0.22; p = 2.2e-15, 1.9e-14, 1.2e-7, 2.5e-9, and 6.6e-7; by order; Fig. 6A). The correlation of AC022034 with CD8 T cells, neutrophils, B cells, dendritic cells and macrophages was shown in Fig. 6B (R = 0.38, 0.27, 0.23, 0.21, and 0.24; p = 2.2e-16, 8.3e-8, 2.1e-10, 2e-6, and 5.5e-8; by order). The correlation of AC104407 with CD8 T cells, neutrophils, B cells, dendritic cells and macrophages was shown in Fig. 6C (R = 0.37, 0.24, 0.21, 0.14, and 0.18; p = 2.2e-16, 2.2e-8, 1.2e-6, 0.0012, and 5e-5; by order). The correlation of HOTAIRM1 with CD4 T cells, neutrophils, B cells, dendritic cells and macrophages was shown in Fig. 6D (R = 0.47, 0.43, 0.31, 0.47, and 0.5; p = 2.2e-16, 2.2e-16, 4.3e-13, 2.2e-16, and 2.2e-16; by order). The correlation of CYTOR with CD4 T cells, neutrophils, B cells, dendritic cells and macrophages was shown in Fig. 6E (R = 0.38, 0.32, 0.19, 0.28, and 0.15; p = 2.2e-16, 4.3e-14, 1.8e-5, 1.4e-10, and 0.00085; by order). The correlation of LINC01831 with CD4 T cells and CD8 T cells was shown in Fig. 6F (R = 0.13, 0.22; p = 0.0024, 2.8e-7; by order). We also estimated the intracellular immune cell content by the CIBERSORT algorithm (Supfig1). Then immune cells with significant correlation (p < 0.001) with LMPS, AC022034, AC104477, HOTAIRM1, CYTOR, LINC01831 were selected and plotted. We found a significant positive correlation between LMPS and the content of CD8T cells (r = 0.33, p = 5.8e-19), neutrophils (r = 0.21, p = 0.00019), and macrophages (r = 0.41, p = 1.5e-13).

3.7 the GSEA analysis of the lncRNA prognosis model

We performed GSEA enrichment analysis according to LMPS in two high and low groups. KEGG analysis (Fig. 7A) revealed that calcium signaling pathway, cytokine-cytokine receptor interaction, JAK_STAT signaling pathway, natural killer cell-mediated cytotoxicity, and T cell receptor signaling pathway were up in high LMPS patients. GO analysis (Fig. 7B) pointed out that humoral immune responses, immune responses regulating cell surface receptor signaling, positive regulation of cytokine production, regulation of immune effector processes, and activation of T cells were up in high LMPS patients. Further analysis revealed that LMPS was positive correlated with the overlapping portion of the T-cell receptor signaling pathway (Fig. 7D) with oncogenic and immune signaling. Finally, we also found that the LMPS and core genes of T cell activation pathways (Fig. 7C) showed overall positive co-expression.

3.8 Construction of an immune-related lncRNA network and validation of the prognosis value of the core genes of the network

We used the Cytoscape to build the correlation network (Fig. 8A) according to the correlation analysis among lncRNA, miRNA, and protein. We also further looked up the hub genes (Fig. 8B) of this network in conjunction with the prediction model of lncRNA above. Finally, we obtained a small network structure containing ten hub genes. We validated the survival curves of hub genes in cancer patients. patients with high expression of MIR4477B (Fig. 8C; p = 0.00045), MIR4666A (Fig. 8E; p = 0.00046) and MIR6071 (Fig. 8D; p < 0.0001) had better survival outcomes. In contrast, patients with high expression of MIR4648 (Fig. 8H; p = 0.00025), AC104407 (Fig. 8F; p = 0.00026), C8G (Fig. 8I; p = 0.036), SEMA6D (Fig. 8G; p = 0.02), AC093726 (Fig. 8K; p < 0.0001), LGR4 (Fig. 8L; p = 0.026) and MIR635 (Fig. 8J; p = 0.0057) had worse survival outcomes.

3.9 Validation of TME and immune cells of hub genes inside the lncRNA network

We next explored the co-expression of immune infiltrating cells (Fig. 9) and hub genes. AC093726(r = 0.29, p = 1.2e-11), MIR635(r = 0.1, p = 0.019) and MIR4648(r = 0.13, p = 0.0032) was positively correlated with sample immunization scores, while MIR4666A (r= -0.27, p = 3.6e-10), MIR6071(r= -0.32, p = 1.2e-13) and SEMA6D(r= -0.11, p = 0.013) was negatively correlated with immunization scores. AC093726(r = 0.35, p = 0.02), AC104407(r = 0.13, p = 0.0027), MIR635(r = 0.24, p = 2e-8) and LGR4(r = 0.15, p = 0.00086) was positively correlated with sample matrix scores, while MIR4477B (r= -0.25, p = 8.2e-9), MIR4666A (r= -0.34, p = 1.1e-15) and MIR6071(r= -0.43, p < 2.2e-16) was negatively correlated with stromal score.

The expression of AC093726, LGR4, MIR635, MIR4648, and SEMA6D were significantly positively correlated with B cells, CD8T cells, neutrophils, macrophages, and dendritic cells, while the expression of MIR6071, MIR4666A, and MIR4477B were significantly negatively correlated with B cells, CD8T cells, neutrophils, macrophages, and dendritic cells (Supfig3).

3.10 The cMap analysis of the LGG

The cMap is where genes, drugs, and disease states are linked by a common gene expression signature^[21]. The L1000 platform^[22] is an expanded cMap. It can be applied to discover the mechanism of action of small molecules, to annotate the function of genetic variants in disease genes, etc. It can be employed to find possible small molecule drugs for diseases based on gene expression as well. We also performed further cMap analysis of the module's immune genes to look for potential drugs that target immunity. Potential therapeutic agents are displayed in Fig. 10A. We obtained the chemical structure (Fig. 10) formulae of loratadine, fidaxomicin, BRD-K86712468, BRD-K89010697, BRD-M52804417, BRD-K52369815, and BRD-K49675259 from PubChem (https://pubchem.ncbi.nlm.nih.gov/).

The WGCNA R package is capable of processing already normalized data for weighted correlation network analysis^[23]. The package can also be used to describe the correlation structure of high-level data structures such as gene expression profiles, proteomics data, etc. We used WCGNA to analyze normalized expression data to obtain a set of immune genes that are closely associated with LGG survival. And we obtained lncRNAs closely related to survival immune genes by correlation analysis between survival-related immune genes and lncRNAs. Finally, by WCGNA of lncRNAs associated with survival-related immune genes, we obtained lncRNA gene modules closely related to survival and immune.

The model made up by AC022034, CYTOR, HOTAIRM1, LINC01831, and AC104407 showed a good prediction of survival in LGG patients, especially at one year. It even has good results in forecasts of glioblastoma. Patients with differential LMPS had different prognoses in both the training and validation models, even in the externally validated dataset. The AUC of ROC is close to 0.6 or above in training and validation sets or TCGA and CGGA datasets. All above suggest a great predictive value of the present model for LGG (even glioma).

IDH1 mutation is a classical LGG molecular marker. IDH1 mutation patients have a better prognosis^[24]. LMPS is lower in both IDH1 mutation and TP53 mutation patients. It is consistent with the better survival of lower LMPS. It indicates the strong association of the model with common LGG molecular markers and the accuracy of predicting patient prognosis by our model.

There were significant differences in LPMS scores between the different clinical phenotypes. LMPS, AC022034, CYTOR, HOTAIRM1, LINC01831 were significantly more highly expressed in older patients, suggesting an effect of age grouping on the model. Histologically typed G3 patients have higher expression of LMPS and CYTOR, HOTAIRM1, LINC01831 compared to G2 patients. It suggests that histological type is a vital factor in the applicability of the model. At the same time, the model is a reference for histological typing. There were differential expressions of LPMS, CYTOR, and HOTAIRM1 in astrocytoma, oligodendroglioma, and mixed cell tumors. It suggests that this model may contribute to the diagnostic staging of gliomas. At the same time, It indicates that different tumor stratification combined with LMPS is further beneficial in determining the prognosis of patients.

TME is vital for revealing the immune status of patients and predicting their immune outcomes^[25]. The ESTIMATE algorithm for estimating the patient's immune microenvironment has been cited in various literature^{[26, 27]}. We explored the correlation between LMPS and patients' TME. The Result suggested that LMPS scores were positive correlated with both patients' immune and stromal scores. It is consistent with previous literature that patients with high CD86^[28] positively relating to stromal and immune scores had worse survival in LGG. It indicates that the model had a deep relationship with TME and could guide the immunotherapy.

Lei X et al. report the role of tumor immune infiltrating cells in tumors and suggest that CAR-T therapies may play a significant role in future treatments^[29]. By TIMER, our study showed that the expression of LMPS, AC0220354, AC104477, and CYTOR had significantly positive correlations with contents of CD8 T cells, neutrophils, dendritic cells, macrophages, and B cells. It suggests a close relationship between this model and immune-cell infiltration. We subsequently obtained the same results in CIBERSORT. The importance of cytotoxic CD8 T cells in mediating adaptive immune responses in various cancers^[30], including gliomas^[31], has been well documented. Liang J et al. ^[32]reported that increased recruitment of neutrophils during anti-vascular endothelial growth factor treatment promotes glioma progression and may increase resistance to treatment. Dendritic cells induce an immune response through antigen presentation, and dendritic cell immunotherapy with inducing anti-tumor immunity has been in clinical trials since the 1990s^[33]. Tumor-related macrophages are recruited into the glioma environment and release a range of growth factors and cytokines in response to those recruitment factors produced by cancer cells. Through this feedback, tumor-associated macrophages promote tumor proliferation, survival, and migration^[34]. The specific role of B lymphocytes in the development of brain tumors remains unclear. In other types of cancer, tumor-infiltrating B cells are associate with the recognition of various tumor antigens^[35]. On the one side, the induction of CD4 + T cell-dependent CD8 + memory T cells contributes to the control of tumor invasion and metastasis^[36]. On the other side, this may activate the phenotype of M2 macrophages and promote tumorigenesis^[37]. The interaction between the role of anti-cancer immune cells, represented by CD8 T cells, and pro-cancer immune cells, expressed by neutrophils and M2 macrophages, may explain the poorer prognosis of patients with high LMPS. Or this may be related to the overproliferation of dysfunctional CD8 T cells^[38]. In any case, the present immune model, and the genes that make up the model, are indisputably linked to the infiltration of tumor immune cells. In other words, this model has the potential to serve as an indicator for LGG immunotherapy. Further, it has the potential to be used as a target for LGG immunotherapy.

To further explore the relationship between this model and immunity, we performed an enrichment analysis of LGG expression data according to LMPS. KEGG analysis revealed significant positive enrichment in calcium signaling pathways, cytokine-cytokine receptor interactions, JAK_STAT signaling pathway, and T cell receptor signaling pathways in the high LMPS group. All are closely related to the immune response and the progression of cancer^[39–43]. GO analysis yielded positive enrichment of humoral immune responses, immune response regulating cell surface receptor signaling, positive regulation of cytokine production, regulation of immune effector processes, and activation of T cells in the high LMPS group. All five of these pathways are closely associated with cancer immunity and significantly inhibit cancer growth^{[44, 45]}. We found that T cell-related pathways were significantly up in KEGG as well as GO analysis, so we performed co-expression analysis of related pathway genes. We caught that LMPS showed an overall positive correlation with the core genes of both activation of T cells and T cell receptor signaling pathways. It suggests that the model is closely associated with T-cell function.

To further explore the potential of this model, through Cytoscape and a pre-correlation analysis, we were able to obtain lncRNA action networks that were highly associated with this immune model. Patients with high expression of MIR4477B, MIR4666A, and MIR6071 had better survival outcomes. In contrast, patients with higher expression of MIR4648, AC104407, C8G, SEMA6D, AC093726, LGR4, and MIR635 had worse survival outcomes. The same with it. Expression of AC093726, LGR4, MIR635, MIR4648, and SEMA6D were significantly positively associated with B cells, CD8T cells, neutrophils, macrophages, and dendritic cells, while the expression of MIR6071, MIR4666A, and MIR4477B were significantly negatively associated with B cells, CD8T cells, neutrophils, macrophages, and dendritic cells. We also analyzed the relationships between hub genes and TME, immune-infiltrating cells. The AC104407 is common in genes that composing model and hub genes. This gene performed well in terms of survival prognosis, immune score, and TME. This suggests that the AC104407 may have a potential role to be explored. It has been reported that MIR6071 inhibits the development of glioblastoma by suppressing the PI3K/AKT/mTOR pathway^[46]. It is another validation of the effectiveness of our lncRNA network. MIR4666A has been reported to have a possible role in femoral necrosis ^[47], atrial fibrillation^[48], and congestive heart failure^[49], but the present study suggests a possible role in glioma immunity and prognosis. AC093726 also performed well in prognostic and immunological analyses, suggesting a possible association with tumor development. We have also performed a drug sensitivity analysis on the immunoprotein contained in hub genes. SEMA6D was positively correlated with the expression of multiple anti-cancer drugs, suggesting that SEMA6D is associated with the anti-cancer effects of multiple cancer drugs, on the contrary, LGR4 was negatively correlated with the expression of multiple anti-cancer drugs, suggesting that LGR4 is associated with cancer drug resistance.

Finally, we further explored possible immunotherapeutic compounds. In total, seven compounds have been the subject of cMap identification. Among these, loratadine has been used for combination chemotherapy for cancers such as breast cancer^[50], melanoma^[51], and gastrointestinal tumor^[52]. However, it has not been used in glioma yet. The results of this study suggest that loratadine may also be an effective agent for adjuvant chemotherapy in LGG. Fidaxomicin acts as an RNA polymerase inhibitor, binding to the RNAP-DNA complex and thereby inhibiting transcription initiation^[53]. The drug is currently in clinical use for refractory clostridial infections^[54]. Combined with the mechanism of action of fidaxomicin and the results of this study, the therapeutic effect of this drug on LGG is expected in the future. The efficacy of BRD-K86712468, BRD-K89010697, BRD-M52804417, BRD-K52369815, BRD-K49675259 in the treatment of LGG is still waiting for more studies.

Overall, an immune-related lncRNA prognostic model was obtained in the LGG training set based on the available immune genes, and its utility was validated in the LGG test set and two glioblastoma datasets. In all datasets, the high LMPS group had worse overall survival. Its prognostic value was further validated by analysis with IDH1 and TP53. The analysis of different clinical phenotypes then illustrates that the model may be able to be used for diagnosis and guiding immunotherapy. We could find a significant positive correlation between LMPS and tumor immune score, stromal score, and the abundance of CD8T cells, neutrophils, macrophages, dendritic cells, and B cells. GSEA got pathways that are closely related to tumors as well as immunity. Despite the complexity of the immune context, the model is indisputable in predicting the immune status of patients. Moreover, in combination with the prognosis of the model, it can help predict the outcome of immunotherapy in LGG patients. Our study also established a core network of interacting lncRNAs and validated the prognostic impact of the hub genes in this network and the correlation with TME. AC093726, LGR4, SEMA6D, AC104477, and MIR6071 are new factors that we found to be significantly associated with LGG prognosis and immunity. Even the therapeutic value of LGR4, SEMA6D was explored. In the end, we also looked for drugs related to possible immune-related treatments by cMap. Although more research is still needed, the study indicates that loratadine and fidaxomicin may both be promising in the adjuvant chemotherapy of LGG.

Data Availability Statement

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare no conflict of interest.

Author Contributions

Xiaoyu Zhang is the first author of this paper, responsible for processing data and writing manuscripts, while others are responsible for finding data.

Funding

This work was supported by the National Natural Science Foundation of China (No.31971166).

Ostrom, Q.T., H. Gittleman, G. Truitt, A. Boscia, C. Kruchko, and J.S. Barnholtz-Sloan, CBTRUS Statistical Report: Primary Brain and Other Central Nervous System Tumors Diagnosed in the United States in 2011-2015. Neuro Oncol, 2018. 20(suppl_4): p. iv1-iv86.
Louis, D.N., A. Perry, G. Reifenberger, A. von Deimling, D. Figarella-Branger, W.K. Cavenee, et al., The 2016 World Health Organization Classification of Tumors of the Central Nervous System: a summary. Acta Neuropathol, 2016. 131(6): p. 803-820.
Ulitsky, I. and D.P. Bartel, lincRNAs: genomics, evolution, and mechanisms. Cell, 2013. 154(1): p. 26-46.
Deveson, I.W., S.A. Hardwick, T.R. Mercer, and J.S. Mattick, The Dimensions, Dynamics, and Relevance of the Mammalian Noncoding Transcriptome. Trends Genet, 2017. 33(7): p. 464-478.
Chen, L.L., Linking Long Noncoding RNA Localization and Function. Trends Biochem Sci, 2016. 41(9): p. 761-772.
Hu, X., A.K. Sood, C.V. Dang, and L. Zhang, The role of long noncoding RNAs in cancer: the dark matter matters. Curr Opin Genet Dev, 2018. 48: p. 8-15.
Li, Z., C. Xu, B. Ding, M. Gao, X. Wei, and N. Ji, Long non-coding RNA MALAT1 promotes proliferation and suppresses apoptosis of glioma cells through derepressing Rap1B by sponging miR-101. J Neurooncol, 2017. 134(1): p. 19-28.
Han, Y., Z. Wu, T. Wu, Y. Huang, Z. Cheng, X. Li, et al., Tumor-suppressive function of long noncoding RNA MALAT1 in glioma cells by downregulation of MMP2 and inactivation of ERK/MAPK signaling. Cell Death Dis, 2016. 7(3): p. e2123.
Vassallo, I., P. Zinn, M. Lai, P. Rajakannu, M.F. Hamou, and M.E. Hegi, WIF1 re-expression in glioblastoma inhibits migration through attenuation of non-canonical WNT signaling by downregulating the lncRNA MALAT1. Oncogene, 2016. 35(1): p. 12-21.
Ke, J., Y.L. Yao, J. Zheng, P. Wang, Y.H. Liu, J. Ma, et al., Knockdown of long non-coding RNA HOTAIR inhibits malignant biological behaviors of human glioma cells via modulation of miR-326. Oncotarget, 2015. 6(26): p. 21934-21949.
Jia, P., H. Cai, X. Liu, J. Chen, J. Ma, P. Wang, et al., Long non-coding RNA H19 regulates glioma angiogenesis and the biological behavior of glioma-associated endothelial cells by inhibiting microRNA-29a. Cancer letters, 2016. 381(2): p. 359-369.
Joyce, J.A. and J.W. Pollard, Microenvironmental regulation of metastasis. Nature reviews. Cancer, 2009. 9(4): p. 239-252.
Quail, D.F. and J.A. Joyce, Microenvironmental regulation of tumor progression and metastasis. Nature medicine, 2013. 19(11): p. 1423-1437.
Liao, K., Y. Lin, W. Gao, Z. Xiao, R. Medina, P. Dmitriev, et al., Blocking lncRNA MALAT1/miR-199a/ZHX1 Axis Inhibits Glioblastoma Proliferation and Progression. Mol Ther Nucleic Acids, 2019. 18: p. 388-399.
Reon, B.J., B. Takao Real Karia, M. Kiran, and A. Dutta, LINC00152 Promotes Invasion through a 3'-Hairpin Structure and Associates with Prognosis in Glioblastoma. Mol Cancer Res, 2018. 16(10): p. 1470-1482.
Mu, M., W. Niu, X. Zhang, S. Hu, and C. Niu, LncRNA BCYRN1 inhibits glioma tumorigenesis by competitively binding with miR-619-5p to regulate CUEDC2 expression and the PTEN/AKT/p21 pathway. Oncogene, 2020. 39(45): p. 6879-6892.
Coureuil, M., H. Lécuyer, S. Bourdoulous, and X. Nassif, A journey into the brain: insight into how bacterial pathogens cross blood-brain barriers. Nat Rev Microbiol, 2017. 15(3): p. 149-159.
Bauer, H.C., I.A. Krizbai, H. Bauer, and A. Traweger, "You Shall Not Pass"-tight junctions of the blood brain barrier. Front Neurosci, 2014. 8: p. 392.
Cserr, H.F., C.J. Harling-Berg, and P.M. Knopf, Drainage of brain extracellular fluid into blood and deep cervical lymph and its immunological significance. Brain Pathol, 1992. 2(4): p. 269-276.
Klemm, F., R.R. Maas, R.L. Bowman, M. Kornete, K. Soukup, S. Nassiri, et al., Interrogation of the Microenvironmental Landscape in Brain Tumors Reveals Disease-Specific Alterations of Immune Cells. Cell, 2020. 181(7): p. 1643-1660.e1617.
Lamb, J., E.D. Crawford, D. Peck, J.W. Modell, I.C. Blat, M.J. Wrobel, et al., The Connectivity Map: using gene-expression signatures to connect small molecules, genes, and disease. Science (New York, N.Y.), 2006. 313(5795): p. 1929-1935.
Subramanian, A., R. Narayan, S.M. Corsello, D.D. Peck, T.E. Natoli, X. Lu, et al., A Next Generation Connectivity Map: L1000 Platform and the First 1,000,000 Profiles. Cell, 2017. 171(6): p. 1437-1452.e1417.
Langfelder, P. and S. Horvath, WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics, 2008. 9: p. 559.
Verhaak, R.G., K.A. Hoadley, E. Purdom, V. Wang, Y. Qi, M.D. Wilkerson, et al., Integrated genomic analysis identifies clinically relevant subtypes of glioblastoma characterized by abnormalities in PDGFRA, IDH1, EGFR, and NF1. Cancer cell, 2010. 17(1): p. 98-110.
Zeng, D., M. Li, R. Zhou, J. Zhang, H. Sun, M. Shi, et al., Tumor Microenvironment Characterization in Gastric Cancer Identifies Prognostic and Immunotherapeutically Relevant Gene Signatures. Cancer immunology research, 2019. 7(5): p. 737-750.
Yoshihara, K., M. Shahmoradgoli, E. Martínez, R. Vegesna, H. Kim, W. Torres-Garcia, et al., Inferring tumour purity and stromal and immune cell admixture from expression data. Nature communications, 2013. 4: p. 2612.
Li, X., Y. Gao, Z. Xu, Z. Zhang, Y. Zheng, and F. Qi, Identification of prognostic genes in adrenocortical carcinoma microenvironment based on bioinformatic methods. Cancer Med, 2020. 9(3): p. 1161-1172.
Qiu, H., W. Tian, Y. He, J. Li, C. He, Y. Li, et al., Integrated Analysis Reveals Prognostic Value and Immune Correlates of CD86 Expression in Lower Grade Glioma. Front Oncol, 2021. 11: p. 654350.
Lei, X., Y. Lei, J.K. Li, W.X. Du, R.G. Li, J. Yang, et al., Immune cells within the tumor microenvironment: Biological functions and roles in cancer immunotherapy. Cancer letters, 2020. 470: p. 126-133.
Huff, W.X., J.H. Kwon, M. Henriquez, K. Fetcko, and M. Dey, The Evolving Role of CD8(+)CD28(-) Immunosenescent T Cells in Cancer Immunology. Int J Mol Sci, 2019. 20(11).
Murphy, K.A. and T.S. Griffith, CD8 T Cell-Independent Antitumor Response and Its Potential for Treatment of Malignant Gliomas. Cancers (Basel), 2016. 8(8).
Liang, J., Y. Piao, L. Holmes, G.N. Fuller, V. Henry, N. Tiao, et al., Neutrophils promote the malignant glioma phenotype through S100A4. Clinical cancer research : an official journal of the American Association for Cancer Research, 2014. 20(1): p. 187-198.
Anguille, S., E.L. Smits, E. Lion, V.F. van Tendeloo, and Z.N. Berneman, Clinical use of dendritic cells for cancer therapy. Lancet Oncol, 2014. 15(7): p. e257-267.
Hambardzumyan, D., D.H. Gutmann, and H. Kettenmann, The role of microglia and macrophages in glioma maintenance and progression. Nat Neurosci, 2016. 19(1): p. 20-27.
Domingues, P., M. González-Tablas, Á. Otero, D. Pascual, D. Miranda, L. Ruiz, et al., Tumor infiltrating immune cells in gliomas and meningiomas. Brain Behav Immun, 2016. 53: p. 1-15.
Nelson, B.H., CD20+ B cells: the other tumor-infiltrating lymphocytes. J Immunol, 2010. 185(9): p. 4977-4982.
Biswas, S.K. and A. Mantovani, Macrophage plasticity and interaction with lymphocyte subsets: cancer as a paradigm. Nat Immunol, 2010. 11(10): p. 889-896.
Li, H., A.M. van der Leun, I. Yofe, Y. Lubling, D. Gelbard-Solodkin, A.C.J. van Akkooi, et al., Dysfunctional CD8 T Cells Form a Proliferative, Dynamically Regulated Compartment within Human Melanoma. Cell, 2019. 176(4): p. 775-789.e718.
Roberts-Thomson, S.J., S.B. Chalmers, and G.R. Monteith, The Calcium-Signaling Toolkit in Cancer: Remodeling and Targeting. Cold Spring Harb Perspect Biol, 2019. 11(8).
Waldmann, T.A. and J. Chen, Disorders of the JAK/STAT Pathway in T Cell Lymphoma Pathogenesis: Implications for Immunotherapy. Annu Rev Immunol, 2017. 35: p. 533-550.
Teicher, B.A. and S.P. Fricker, CXCL12 (SDF-1)/CXCR4 pathway in cancer. Clinical cancer research : an official journal of the American Association for Cancer Research, 2010. 16(11): p. 2927-2931.
Sordo-Bahamonde, C., S. Lorenzo-Herrero, R. Payer Á, S. Gonzalez, and A. López-Soto, Mechanisms of Apoptosis Resistance to NK Cell-Mediated Cytotoxicity in Cancer. Int J Mol Sci, 2020. 21(10).
Chandran, S.S. and C.A. Klebanoff, T cell receptor-based cancer immunotherapy: Emerging efficacy and pathways of resistance. Immunol Rev, 2019. 290(1): p. 127-147.
Shishido, S.N., S. Varahan, K. Yuan, X. Li, and S.D. Fleming, Humoral innate immune response and disease. Clin Immunol, 2012. 144(2): p. 142-158.
Borst, J., T. Ahrends, N. Bąbała, C.J.M. Melief, and W. Kastenmüller, CD4(+) T cell help in cancer immunology and immunotherapy. Nature reviews. Immunology, 2018. 18(10): p. 635-647.
Zhou, Y., H. An, and G. Wu, MicroRNA-6071 Suppresses Glioblastoma Progression Through the Inhibition of PI3K/AKT/mTOR Pathway by Binding to ULBP2. Onco Targets Ther, 2020. 13: p. 9429-9441.
Liu, G.Z., C. Chen, N. Kong, R. Tian, Y.Y. Li, Z. Li, et al., Identification of potential miRNA biomarkers for traumatic osteonecrosis of femoral head. Journal of cellular physiology, 2020. 235(11): p. 8129-8140.
Xu, G., Y. Cui, Z. Jia, Y. Yue, and S. Yang, The Values of Coronary Circulating miRNAs in Patients with Atrial Fibrillation. PLoS One, 2016. 11(11): p. e0166235.
Gu, H., L. Chen, J. Xue, T. Huang, X. Wei, D. Liu, et al., Expression profile of maternal circulating microRNAs as non-invasive biomarkers for prenatal diagnosis of congenital heart defects. Biomed Pharmacother, 2019. 109: p. 823-830.
Kirshner, J.J., M.C. McDonald, 3rd, F. Kruter, A.S. Guinigundo, L. Vanni, C.L. Maxwell, et al., NOLAN: a randomized, phase 2 study to estimate the effect of prophylactic naproxen or loratadine vs no prophylactic treatment on bone pain in patients with early-stage breast cancer receiving chemotherapy and pegfilgrastim. Support Care Cancer, 2018. 26(4): p. 1323-1334.
Fritz, I., P. Wagner, M. Bottai, H. Eriksson, C. Ingvar, I. Krakowski, et al., Desloratadine and loratadine use associated with improved melanoma survival. Allergy, 2020. 75(8): p. 2096-2099.
Chen, T., Y. Hu, B. Liu, X. Huang, Q. Li, N. Gao, et al., Combining thioridazine and loratadine for the treatment of gastrointestinal tumor. Oncol Lett, 2017. 14(4): p. 4573-4580.
Venugopal, A.A. and S. Johnson, Fidaxomicin: a novel macrocyclic antibiotic approved for treatment of Clostridium difficile infection. Clin Infect Dis, 2012. 54(4): p. 568-574.
Rao, K. and P.N. Malani, Diagnosis and Treatment of Clostridioides (Clostridium) difficile Infection in Adults in 2020. Jama, 2020. 323(14): p. 1403-1404.

supfig1.tif
Supfig1: Correlation between LMPS, AC022034, AC104477, HOTAIRM1, CYTOR, LINC01831 and the immune cell predicted by the CIBERSORT algorithm.
supfig2.tif
Supfig2: Drug-sensitive type analysis with SEMA6D, C8G, and LGR4 expression
supfig3.tif
Supfig3: immune cell infiltration analysis of hub genes of the network
Supfig4.tif
Supfig4: LMPS expressions in not/IDH1 mutation (A) and not/TP53 mutation (B) groups

Download PDF

Version 1

posted

You are reading this latest preprint version

Construction of Immune-related lncRNA Prognosis Models and Development Related Networks and Drugs in LGG

Status:

Version 1

Abstract

Figures

1. Introduction

2. Materials And Methods

2.1 Data processing

2.2 Weighted Correlation Network Analysis (WCGNA)

2.3 Model Construction

2.4 Differential expression of LMSP and constitutive model genes in different clinical phenotypes

2.5 Relationship between TME and LMPS

2.6 Correlation between immune cells and LMPS

2.7 Construction of lncRNA-related action network and screening of core genes

2.8 survival analysis of hub genes

2.9 TME analysis of hub genes

2.10 Gene Set Enrichment Analysis (GSEA) and co-expression of enriched pathway genes with prognostic models

2.11 Connectivity Map (cMap) analysis

2.12 IDH1mutation and TP53mutation

2.13 Statistical methods

3. Result

3.1 the WCGNA analysis of immune protein-coding genes

3.2 The WCGNA analysis of the immune-related lncRNA

3.3 The prognosis model of lncRNA in the pink module

3.4 The Correlation of LMPS and clinical phenotypes

3.5 Relationship between LMPS and tumor immune microenvironment (TME)

3.6 The correlation of immune cells and lncRNA prognosis model.

3.7 the GSEA analysis of the lncRNA prognosis model

3.8 Construction of an immune-related lncRNA network and validation of the prognosis value of the core genes of the network

3.9 Validation of TME and immune cells of hub genes inside the lncRNA network

3.10 The cMap analysis of the LGG

4. Discussion

Declarations

References

Supplementary Files

Status:

Version 1