An integrated bioinformatic investigation of focal adhesion-related genes in glioma followed by preliminary validation of COL1A2 in tumorigenesis

Focal adhesions (FAs) allow cells to contact the extracellular matrix, helping to maintain tension and enabling signal transmission in cell migration, differentiation, and apoptosis. In addition, FAs are associated with changes in the tumor microenvironment (TME) that lead to malignant progression and drug resistance in tumors. However, there are still few studies on the comprehensive analysis of focal adhesion-related genes (FARGs) in glioma. Expression data and clinical information of glioma samples were downloaded from public databases. Two distinct molecular subtypes were identified based on FARGs using an unsupervised consensus clustering algorithm. A scoring system consisting of nine FARGs was constructed using integrated LASSO regression and multivariate Cox regression. It not only has outstanding prognostic value but also can guide immunotherapy of glioma patients, which was verified in TCGA, CGGA, GSE16011, and IMvigor210 cohorts. The results of bioinformatics analysis, immunohistochemistry staining, and western blotting all revealed that the expression of COL1A2 was up-regulated in glioblastoma and related to poor prognosis outcomes in patients from public datasets. COL1A2 promotes the proliferation, migration, and invasion of glioblastoma cells. A positive correlation between COL1A2 and CD8 was determined in GBM specimens from eight patients. Moreover, the results of cell co-cultured assay showed that COL1A2 participated in the killing of GBM cells by Jurkat cells. Our study indicates that the FARGs have prominent application value in the identification of molecular subtypes and prediction of survival outcomes in glioma patients. Bioinformatics analysis and experimental verification provide a direction for further research on FARGs.


INTRODUCTION AGING
signal transmission [8,9] during cell migration, differentiation, and apoptosis. FA proteins are classified according to their functions [10] as integrins, FA kinases, paxillin, etc. Some of these have been testified linked to cancer progression [11][12][13]. However, the predictive value of FA proteins for prognosis and treatment response in cancer is still unknown. TME is closely related to tumorigenesis and malignant progression [14]. The change in TME, which reduces the adhesion of tumor cells, has an important impact on the drug resistance of tumors [15] and can cause metastasis [16,17] by promoting epithelialmesenchymal transition [18][19][20][21]. As an important component of TME, macrophages account for 30% among the immune cells, which is the highest proportion of cells in the glioma [22]. Immune checkpoint blockades (ICBs) can prevent immune escape by modulating the function of T cells [23,24]. At present, CTLA-4, PD-1, and PD-L1 are the main immunotherapeutic targets for patients with advanced tumors, although they do not always prolong the survival time of patients [23,25]. The failure of tumor immunotherapy is thought to be related to an immunosuppressive TME and low tumor mutational burden (TMB) [26,27]. Up to now, glioma patients have neither effective immunotherapeutic targets nor biomarkers that can effectively predict immunotherapeutic response.
In this paper, we not only identified two distinct molecular subtypes but also constructed a scoring system with outstanding clinical application value. In addition, we screened the gene COL1A2, which is upregulated in GBM tissue and closely related to the poor prognosis of glioma patients. In addition, we conducted intensive study on the COL1A2 gene and found that it not only participates in the regulation of biological behavior of GBM but also may be a key molecule in the TME of glioma.

Molecular classification based on FARGs in TCGA cohort
First, we show the workflow of this study in Figure 1. According to the flow chart, 192 FARGs were obtained from the intersection of TCGA, CGGA, and GSE16011 cohorts (Figure 2A). A novel molecular classification was identified by using an unsupervised clustering algorithm based on the expression profiles of the 192 FARGs extracted from the TCGA cohort. According to the area under the CDF curves, probably approximately correct algorithm and correlations between clusters, we determined the optimal number of clusters is two (k=2) ( Figure 2B-2E). Principal components analysis (PCA) can effectively distinguish glioma patients with the FARG1 subtype and FARG2 subtype ( Figure 2F).

Basic analysis between different subtypes in the TCGA cohort
The OS and progression-free survival (PFS) of glioma patients, low-grade glioma (LGG) patients, and GBM patients with the FARG1 subtype were significantly longer than those with the FARG2 subtype ( Figure 2G, 2H). Except for gender, there are dramatic differences in the distribution of the other four clinicopathological characteristics between FARG1 and FARG2 subtypes. In addition, most FARGs are differentially expressed between the two subtypes ( Figure 2I). We simultaneously performed GSVA in TCGA, CGGA, and GSE16011 cohorts and identified 167 functional pathways with different relative activity between the two subtypes ( Figure 2J and Supplementary Table 3), of which 65 functional pathways were more active in FARG1 subtype, such inositol phosphate metabolism, glycerol lipid metabolism and lysine degradation, and the other 102 functional pathways were more active in FARG2 subtypes, such as cell lung cancer, pancreatic cancer, and primary immunodeficiency.

Immune-related analysis between different subtypes in TCGA cohort
Given our above finding, we further investigated whether different subtypes have different immunological characteristics. The results revealed that there were remarkable differences in most immune signatures (except for dendritic cells and mast cells) and all four TME-related scores among the two subtypes ( Figure 3A, 3C-3F). Compared with the FARG1 subtype, the glioma tissue of the FARG2 subtype contains more immune cells and stromal cells. We also estimated the content of immune cells and found that the content of almost all immune cells was different between the two subtypes ( Figure 3B and Supplementary Table 4). Moreover, we also found that except for TNFSF9, other immune checkpoints (ICPs) were differentially expressed between the two subtypes of FARG1 and FARG2 ( Figure 3G). Considering the effect of gene mutations on tumorigenesis and progression, we found that the TMB was dramatically different between the two subtypes ( Figure 3H).

Prognostic analysis of scoring system
Since there are significant differences in the prognosis outcomes between GBM and LGG patients, IDH mutant and wild-type glioma patients, we should not only explore the prognostic value of the scoring system in pan-gliomas but also explore its prognostic value among different grades and different IDH mutation types. If the risk score of glioma patients was higher than the median value, they were classified as a high-risk group, otherwise, they were classified as a low-risk group. K-M survival curves showed that in the three cohorts, the OS of patients with pan-gliomas ( Figure 4D and Supplementary Figure Figures 3B, 4B, 4H) in the low-risk group tended to be longer than that of patients in the high-risk group. Moreover, the risk score can effectively predict the 1, 3, and 5-year OS rate of patients with pan-gliomas ( Figure 4F and  Table 6). In addition, in the TCGA cohort, patients with pangliomas, IDH mutant and IDH wildtype in the low-risk group also tended to have longer PFS than patients in the high-risk group ( Figure 4E and Supplementary Figure 3G, 3H), and the risk score can also effectively predict the 1, 3, 5-year PFS rate of patients with pangliomas, IDH mutant and IDH wildtype ( Figure 4G, and Supplementary Figure 3I, 3J and Supplementary  Table 6). Next, we explored the differential distribution of five clinicopathological characteristics and the differential expression of nine FARGs between the two risk groups in the three cohorts. The results showed that except for gender, there were significant differences in the distribution of the other four clinicopathological characteristics and the expression of nine FARGs between different risk groups ( Figure 4H and Supplementary Figure 1E, 1G). Finally, we explored whether the scoring system can be used as an independent prognostic factor for patients with pangliomas, LGG, GBM, IDH mutant, or wildtype. The results showed that the independent prognostic value of risk score was superior to that of age, gender, grade, IDH status, and 1p19q status in TCGA, CGGA, and GSE16011 cohorts ( Figure 4I, 4J, and Supplementary Figures 1F, 1H, 3E, 3F, 3K, 3L, 4E, 4F, 4K, 4L).

Genomic variation analysis between high and lowrisk groups
Because genomic variation can impact tumor immunity and immune cell infiltration patterns [28], we explored the association between genomic variation and risk score in this study. The top-16 genes with the highest mutation frequency have significant differences between the two risk groups ( Figure 5A, 5B). There are also obvious differences in the distribution of mutation and wildtype of six wellknown genes (TP53, PTEN, IDH1, EGFR, ATRX, TTN) between high and low-risk groups ( Figure 5C). Besides, we found that glioma patients with low-risk scores inclined to have a lower TMB ( Figure 5D). Moreover, patients in the low-risk group had a lower frequency of CNV, either amplification or deletion ( Figure 5E, 5F).

Immune-related analysis between high and low-risk groups
By performing GSVA in TCGA, CGGA, and GSE16011 cohorts, we found that the relative activities of 156 functional pathways were different between the two risk groups: the relative activity of 53 functional pathways increased in the low-risk group, mainly focusing on metabolic-related pathways, such as "butanoate metabolism", "propanoate metabolism" and were lower in the low-risk group than that in the high-risk score group. (G) GSVA between FARG1 and FARG2 subtypes. Red and blue represent the relative activation and inhibition of the pathways, respectively. AGING "glyoxylate and dicarboxylate metabolism"; the relative activity of other 103 functional pathways increased in the high-risk group, including immune-related pathways ( Figure 5G and Supplementary Table 7). Next, we began an in-depth analysis of the relationship between immunologic characteristics and risk scores. We found that the TME-related scores were not only highly correlated with the risk score (Supplementary Figure  5C-5E) but also significantly different between the high and low-risk groups in the three cohorts ( Figure 6A and  Figure 5A, 5B). Then, we also found that the enrichment scores of most immune signatures calculated by ssGSEA were different between the two risk groups in the three cohorts ( Figure 6A and Supplementary Figure 5A, 5B). Considering that immune cells play a key role in tumorigenesis, tumor progression, and tumor immunity, we further explored the association between the scoring system and immune cells. The content of most immune cells was not only correlated with the risk score (Supplementary Figure 6A-6C) but also significantly different among different risk groups ( Figure 6B and Supplementary Figure 6D, 6E). The results of the analysis between ICPs and risk scores ( Figure 6C and Supplementary Figure 7A-7E) prompted us to continue to explore whether the scoring system can effectively predict the immunotherapy response of glioma patients. There was a significant difference of the risk score between the responder and non-responder groups ( Figure 6D and Supplementary Figure 8A-8D).
The proportion of responders in the high-risk group is much higher than those in the low-risk group, which indicates that glioma patients in the high-risk group are more likely to benefit from immunotherapy ( Figure 6E and Supplementary Figure 8B-8E). The subclass mapping analysis indicated that glioma patients in the high-risk group are more sensitive to anti-PD-1 therapy (Bonferroni P = 0.016,0.032, 0.039 in GSE16011, TCGA, and CGGA cohorts, respectively) ( Figure 6F and Supplementary Figure 8C-8F). If the Bonferroni corrected P > 0.05, the P-value would not be marked in the subgroup map, but from the color comparison of the corresponding patches in the subgroup map, it can be seen that the P-value of anti-CTLA4 therapy in the highrisk group is lower than that in the low-risk group, indicating that glioma patients in the high-risk group are more likely to respond to anti-CTLA4 therapy ( Figure 6F and Supplementary Figure 8C-8F).

Validation of the clinical application value of the scoring system in the IMvigor210 (mUC) cohort
We selected the IMvigor210 (mUC) cohort to further validate the clinical application value of the scoring system. K-M survival curves indicated that the OS of patients with a low-risk score was longer than those with a high-risk score ( Figure 6H). Moreover, the risk score can effectively predict the 1, 3, and 5-year OS rates of patients in the IMvigor210 cohort ( Figure  6G). Then, the differential expression of PDL1 (CD274) between high and low-risk groups was also explored. The results showed that the expression level of PDL1 in the high-risk group was higher than that in the low-risk group ( Figure 6I). Furthermore, the proportion of CR/PR and PD/SD in the high-risk group is higher than that in the low-risk group ( Figure 6J, 6K).

COL1A2 is up-regulated in GBM
First, compared with the other eight FARGs, the mRNA expression of COL1A2 in GBM tissues was significantly higher than that in normal brain tissues (NBT), which was verified on the GEPIA website and three independent GBM cohorts ( Figure 7A, 7B). Then, immunohistochemical staining of eight pairs of GBM tissues and adjacent tissues showed that COL1A2 was significantly up-regulated in GBM tissues ( Figure 7C, 7D). In addition, at the tissue protein level, the expression COL1A2 in GBM tissues was also significantly higher than that in corresponding adjacent tissues ( Figure 7E, 7F). Finally, the K-M survival curves of COL1A2 in four independent GBM cohorts (Rembrandt, GSE16011, TCGA, and CGGA cohorts) showed that the prognosis outcomes of GBM patients in the high expression group were significantly shorter than that of GBM patients in the low expression group ( Figure 7G).

COL1A2 promotes the malignant progression of GBM cells in vitro
By comparing the differential expression of COL1A2 between different cell lines at the transcription level and protein level, we found that the expression level of COL1A2 was the highest in U87 cells and the lowest in U251 cells ( Figure 8A-8C). Then, we constructed knockdown and overexpression plasmids of COL1A2 to further explore the effect of COL1A2 on the biological behavior of GBM in vitro. The result of RT-qPCR and western blotting showed that the constructed knockdown and overexpression plasmids could significantly reduce or increase the expression of COL1A2 ( Figure 8D-8I). The results of colony formation assay ( Figure 8J-8M), CCK-8 assay ( Figure  8N, 8O), and EdU staining ( Figure 8P, 8Q) suggest that COL1A2 can significantly enhance the proliferation of GBM cells. In addition, the outcomes of wound healing and transwell assay support COL1A2 promoting the migration and invasion of GBM cells ( Figure 8R-8U). These results suggest that COL1A2 plays an important role in the progression of GBM cells.

Detection of the viability of GBM cells in vitro
In eight cases of GBM tissues, we found that the tissues with high expression of COL1A2 were accompanied by high expression of CD8. On the contrary, the tissues with low expression of COL1A2 were accompanied with low expression of CD8 ( Figure 9A), which to some extent suggested that the expression of COL1A2 might induce the infiltration of CD8 T cells or the infiltration COL1A2. Then, we studied whether COL1A2 played a of CD8 T cells might induce the expression of role in AGING the process of Jurkat cells acting on GBM cells in vitro. Both methods, PMA, and ionomycin or anti-CD3 and anti-CD28 antibodies, can effectively activate Jurkat cells into CD8 + Jurkat cells ( Figure 9B). CD8 T cells transiently overexpress CD69 in the early stage of activation, which is considered to be a costimulatory signal of T cell proliferation. However, we did not detect the expression of CD69 protein in Jurkat cells stimulated by anti-CD3 and anti-CD28 antibodies ( Figure 9C), so we utilized PMA and ionomycin to stimulate Jurkat cells. The results of the cell viability assay showed that activated Jurkat cells could effectively inhibit the viability of U87 ( Figure 9D) and U251 ( Figure 9E) cells. With the increase of COL1A2 expression, the lethality of activated Jurkat cells to GBM cells decreases, which may be related to COL1A2 promoting the proliferation of GBM cells.

DISCUSSION
The prognosis of glioma patients remains poor. As such, it is critical to identify new prognostic biomarkers and develop a more effective treatment. Although there have been multiple studies on immunotherapy, especially by ICB in GBM patients, the results of phase III clinical trials have not been satisfactory compared to other tumors [29,30]. Neither anti-CTLA-4 antibody alone nor in combination with anti-PD-1 antibody yielded long-term survival benefits for glioma patients [31]. Numerous factors can affect the response of glioma to immunotherapy including factors in the TME and TMB. The latter is a biomarker for evaluating the therapeutic effect of anti-PD-1 antibody [32]; and converting the immunosuppressive TME to a immunostimulatory one is an effective treatment strategy [33,34]. In this study, we developed a novel molecular classification for glioma patients based on FARGs. Not only the OS and PFS outcomes of glioma patients are significantly AGING different between different subtypes, but also the characteristics of TME are significantly different between different subtypes. However, we did not further validate the properties of the new molecular classification in multiple glioma cohorts, making it difficult to effectively validate its clinical application value. In addition, a scoring system was constructed based on nine FARGs (SPP1, THBS4, ERBB2, ELK1, COL1A2, PTEN, CDC42, FLNC, and PDGFA). The expression levels of these FARGs differed between different WHO grades. PTEN regulates signaling pathways related to cell growth and survival [35,36] and cell metabolism [37]. COL1A2 participates in collagen synthesis [38] but has been implicated in the immune response [39]. ELK1 is a transcription factor that activates target genes via some protein kinase/ regulatory kinase pathways [40,41]. ERBB2 is a marker gene in breast cancer. SPP1 is known to be overexpressed in many malignant tumors including glioma [42][43][44][45][46] and regulates cell growth, proliferation, apoptosis, and migration [47]. THBS4 is a member of the thrombospondin protein family and plays important roles in wound healing and tissue repair [48-51], intracellular migration, adhesion, and proliferation [52][53][54]. FLNC is an action-binding filamin protein that regulates actin reorganization-dependent processes such as differentiation, migration, and proliferation of cells [55]. As a Rho family GTPase, CDC42 plays a key role in the activation of signaling cascades that regulate cell adhesion, cytoskeletal composition, proliferation, and migration and is therefore important for the malignant transformation of tumors [56]. In this study, COL1A2 was screened from nine FARGs used to construct the scoring system using bioinformatics analysis methods for further research, but no corresponding studies were conducted on the other eight FARGs, which will be the focus of our future research.
The clinical application value of the scoring system has been verified in TCGA, CGGA, and GSE16011 cohorts.
The scoring system can not only effectively predict the prognosis outcome of glioma patients, but also can be used as an independent prognostic factor for glioma patients. TMB can predict the immunotherapeutic response of some types of tumors and has a close relationship with the scoring system. There are also significant differences in the distribution of mutation types of some classic genes between high and low-risk groups. Moreover, the TME between high and low-risk groups is also significantly different, which will directly affect the immunotherapeutic response of glioma patients. In addition to verifying the clinical application value of the scoring system in different glioma cohorts, we also found that the scoring system has an important clinical application value in the imvigor210 (mUC) cohort. Although we have conducted a comprehensive bioinformatics analysis of the scoring system, we have not conducted a comparative study on its clinical application value with other existing scoring systems, resulting in a lack of horizontal comparison.
To further study the role of FARGs in the progression of GBM, COL1A2 was screened and further studied. The screening process of COL1A2 is simple, and if it. can be screened from multiple perspectives, it would make the screening process more convincing. COL1A2 is upregulated in GBM tissue and promotes malignant progression of GBM cells, but its specific mechanism has not been thoroughly studied, which will be our future research focus. It has been suggested that the poor response to immunotherapy in glioma patients is attributable to immunosuppression in the brain [27], however, the specific mechanism is still unclear at present. COL1A2 was still taken as the research object. It was found that it may promote the infiltration of CD8 T cells, but at the same time may inhibit the lethality of CD8 T cells on tumor cells. As for the specific role of COL1A2 in the tumor microenvironment of glioma and in the immunotherapy process, if further research can be made in this study, it would make this study more valuable in clinical application.
In summary, we identified two distinct subtypes based on FARGs. More importantly, we established a scoring system with great clinical application value, which can not only effectively predict the prognosis outcome of glioma patients, but also predict the immunotherapy response for glioma patients. In addition, we screened COL1A2 and verified its involvement in the progression of GBM cells in vitro. It is expected that this study can contribute to the diagnosis and treatment of glioma patients.

Data sources
The flow diagram of our study was shown in Figure 1.  [57]. The genomic identification algorithm of important targets in tumors is used to process CNV data [58]. We searched 199 FARGs from the Molecular Signatures Database and the details of these FARGs can be found in Supplementary Table 2. In addition, the IMvigor210 cohort treated with the PD-L1 inhibitor was also included in this study to confirm our findings.

Acquisition of GBM tissue samples
The current research was authorized by the Ethics Committee of The First People's Hospital of Fuzhou City and informed consent was obtained from the patients. The surgical specimens of GBM patients were removed to liquid nitrogen immediately after surgical excision.

Unsupervised consensus clustering based on FARGs
Different glioma cohorts have different sequencing methods and gene annotation information. To facilitate follow-up analysis, we took the intersection of 199 FARGs in three independent glioma cohorts and finally obtained 192 FARGs. A novel molecular classification was identified based on the 192 FARGs extracted from the TCGA cohort by unsupervised consensus clustering with 1000 iterations using the "ConsensusClusterPlus" package [59]. The optimal value of molecular classification should consider not only the rate of increase of the area under the cumulative distribution function (CDF) curves but also consider the correlation between different molecular classifications.

Gene set variation analysis (GSVA) between different subtypes
After downloading "c2.cp.kegg.v7.2.symbols" from the Molecular Signatures Database, GSVA between different subtypes was performed using the R package "GSVA" [60]. Differences in the relative activity of functional pathways between different subtypes were explored by using the R package "limma" [61].

Immune-related analysis between different subtypes
Different types of cell components in tumor tissues can be estimated by calculating TME-related scores using the R package "ESTIMATE" [62]. The content of different types of immune cells was reckoned using MCPCOUNTER, CIBERSORT, and QUANTISEQ algorithms [63,64]. Enrichment scores of 29 immune signatures were determined by using a single sample gene set enrichment analysis (ssGSEA), which to some extent represents the immune activity within the tumor. [65]. The analyzed immune checkpoint proteins (ICPs) were selected from a previous study [66]. The above immune-related characteristics were analyzed in different subtypes.

Construction of a scoring system and verification of its prognostic value
We screened 175 FARGs with prognostic significance from 192 FARGs using uni-Cox regression analysis. AGING 175 FARGs were further screened by lasso algorithm and multi-Cox regression analysis, and nine FARGs were finally obtained to construct a scoring system: risk score = sum (gene expression × coefficient). High and low-risk groups were classified according to the median value of the risk score. The prognostic value of the scoring system was evaluated by using Kaplan-Meier (K-M) survival curves (R package "survminer"), receiver operating characteristic (ROC) curves (R package "survivalROC"), uniand multi-Cox regression analysis (R package "survival" and "forestplot"). The different distribution of different clinicopathological characteristics between high and low-risk groups was analyzed by R package "limma".

Genomic variation analysis between high and lowrisk groups
The mutation type and frequency of genes between high and low-risk groups were analyzed by the R packages "maftool" [57] and "GenVisR" [67]. The relationship between TMB, as a marker for the therapeutic efficacy of anti-PD-1 antibodies in other cancers [32,68], and risk score was also analyzed by R package "ggpubr". As genomic alterations can impact tumor immunity and immune infiltration patterns [69,70], we compared amplifications and deletions between high and low-risk groups and visualized the results as circle graphs using the R package "RCircos" [71].

Immune-related analysis of risk score
At present, only a sub-fraction of patients achieved longlasting clinical benefits from ICBs treatment. To effectively predict the response of tumor patients to ICBs, the TIDE algorithm was developed mainly to model two primary mechanisms of tumor immune evasion: inducing the dysfunction of cytotoxic T lymphocytes (CTLs) in tumors and preventing the infiltration of CTLs in tumor tissues [72]. T cell dysfunction was identified by measuring the interaction between each gene and the infiltration level of CTL to influence patient survival. The TIDE algorithm explored the association between gene expression data and markers of T cell dysfunction. Tumor samples that were highly positively correlated with markers of T cell dysfunction were identified as non-responders and otherwise as responders. Finally, the TIDE algorithm and an unsupervised subclass mapping method [73] were used together to forecast the response of glioma patients to anti-PD-1 and anti-CTLA-4 immunotherapy.

RT-qPCR and Western blotting
The The extraction and detection methods of protein and total RNA in different cell lines were completely consistent with those described in our previous study [74]. Here, we will describe the extraction steps of tissue protein in detail. First, the tissues were weighed and sheared into 2ml EP tubes. RIPA and PMSF were mixed in a ratio of 100:1 to prepare tissue lysates. 1ml of lysate was added to every 50mg of sheared tissue and homogenized with a homogenizer, and then lysed on ice for 30 minutes. Finally, the mixture of lysate and tissue was centrifuged with a low-temperature high-speed centrifuge at 4° C and 12000 rpm for 10 minutes. The subsequent steps are the same as the extraction of cellular proteins. Both RT-qPCR and western blotting assays were repeated three times.

Immunohistochemical and immunofluorescence staining
First, eight pairs of GBM and adjacent tissues were fixed with 10% formalin for seven days, and then these tissues were embedded with paraffin and sectioned. After dewaxing and dehydration, the tissue sections were treated with 3% hydrogen peroxide for 10 minutes. Subsequently, the tissue sections were blocked with 5% BSA and then incubated overnight with a primary antibody against COL1A2 (Servicebio, GB13022-2, 1:1000 dilution) at 4° C. Then the tissue sections were treated with the corresponding secondary antibodies at room temperature for one hour. Finally, DAB staining, target molecule detection, and hematoxylin re-staining were carried out in turn. For immunofluorescence staining, tissue sections were immune-stained overnight with primary antibodies against COL1A2 (Servicebio, GB13022-2, 1:1000 dilution) and CD8 (Proteintech, 1:400) at 4° C, and then incubated with fluorochrome-conjugated antibodies. DAPI was added as a nuclear counterstain.

Cell proliferation assay
Cells were seeded into 6-well plates at 1000 cells/well and cultured for about 15 days. It was predicted that the cells would proliferate for about 5-7 generations. The culture medium was changed every three days, and the colony formation was closely observed.  D-35578). In addition, we also performed the migration assay in the same way as the invasion assay, except that the upper chamber was not covered with Matrigel. Both invasion and migration assays were repeated three times.

Tumor cell viability assay
In light of the instructions offered by the manufacturer, the lentivirus expressing luciferase designed by Sheweisi Biotechnology Company (Tianjin, China) AGING was transfected into U251 and U87 cells and screened with neomycin. Then COL1A2 knockdown and overexpression plasmids were transfected into cells stably expressing luciferase, respectively. When these cells grew to about 90%, they have inoculated with Jurkat cells in a ratio of 1:1 to 96 well plates with a white and clear bottom (total 1.0×10 4 cells/ml). Cocultured cells were taken at different time points (8,16,24,36,48 h), the supernatant was discarded, washed with PBS, cell lysate was added, and then D-luciferin potassium salt (D1009, UElandy, China) was added. After the above operations are completed, the fluorescence intensity is detected using a multifunctional microplate reader. The experiment was repeated three times.

Statistical analyses
Spearman and Pearson analysis methods were used to determine correlations between different grouping variables. The unpaired Student's t-test and Mann-Whitney U-test were used to analyzing normally and non-normally distributed data, respectively. The Wilcoxon rank-sum test was used to compare two groups or categories; for more than two groups or categories, the Kruskal-Wallis test was used. Survival was assessed by K-M analysis with the log-rank test. Cox regression analyses were carried out to evaluate the prognostic value and stability of the risk prediction model. ROC curves were used to assess the prognostic value of 1-, 3-, and 5year OS. TIIC infiltration level was analyzed using three algorithms (CIBERSORT, QUANTISEQ, and MCPCOUNTER). All statistical tests were two-sided and P<0.05 was considered statistically significant.

AUTHOR CONTRIBUTIONS
Yufan Zhou and Xiang Zhou designed this study and directed the writing of the article. Guojun Yao did all the experiments of this study and participated in writing the article. Ling Deng performed the data analysis and the figures were plotted. Xinquan Long performed the data acquisition and critical reading of the manuscript. All authors have carefully read and consented to the final manuscript.

CONFLICTS OF INTEREST
The authors state that the study was conducted without any direct or indirect conflict of interest.

ETHICAL STATEMENT AND CONSENT
This study was carried out in accordance with the principles of the Declaration of Helsinki. Approval was granted by the Ethics Committee of The First People's Hospital of Fuzhou City. Informed consent was obtained from the patients included in the study. Written informed consent for publication was also obtained from the patients enrolled in the study.

FUNDING
There is no fund to support this study, but it has received support from the First People's Hospital of Fuzhou City. AGING Supplementary Figure 8. Differential response to immunotherapy between high and low-risk groups in the CGGA and GSE16011 cohorts. (A) The difference in risk score between two immunotherapy response groups in the CGGA cohort. (B) The proportion of responders and non-responders to immunotherapy between high and low-risk groups in the CGGA cohort. (C) Subgroup map of predicted response to ICB therapy in the two risk groups in the CGGA cohort. (D) The difference in risk score between two immunotherapy response groups in the GSE16011 cohort. (E) The proportion of responders and non-responders to immunotherapy between high and low-risk groups in the GSE16011 cohort. (F) Subgroup map of predicted response to ICB therapy in the two risk groups in the GSE16011 cohort. In (A, D), the upper and lower lines of the boxes indicate the interquartile range of values, and the lines in the boxes represent the median value, and the black dots show outliers.