Jinsong Wang
Xue Cui1
Yiming Weng
Jiayan Wei
Xinyi Chen
Peiwei Wang
Jian Qin
Min Peng- 1Department of Oncology, Renmin Hospital of Wuhan University, Wuhan, China
- 2Central Laboratory, Renmin Hospital, Wuhan University, Wuhan, China
Lung adenocarcinoma (LUAD) is an essential pathological subtype of non-small cell lung cancer and offers a severe problem for worldwide public health. There is mounting proof that angiogenesis is a crucial player in LUAD progression. Consequently, the purpose of this research was to construct a novel LUAD risk assessment model based on genetic markers related to angiogenesis. We accessed The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases for LUAD mRNA sequencing data and clinical information. Based on machine algorithms and bioinformatics, angiogenic gene-related risk scores (RS) were calculated. Patients in the high-risk category had a worse prognosis (p < 0.001) in the discovery TCGA cohort, and the results were confirmed by these three cohorts (validation TCGA cohort, total TCGA cohort, and GSE68465 cohort). Moreover, risk scores for genes involved in angiogenesis were independent risk factors for lung cancer in all four cohorts. The low-risk group was associated with better immune status and lower tumor mutational load. In addition, the somatic mutation study revealed that the low-risk group had a lower mutation frequency than the high-risk group. According to an analysis of tumor stem cell infiltration, HLA expression, and TIDE scores, the low-risk group had higher TIDE scores and HLA expression levels than the high-risk group, and the amount of tumor stem cell infiltration correlated with the risk score. In addition, high-risk groups may benefit from immune checkpoint inhibitors and targeted therapies. In conclusion, we developed an angiogenesis-related gene risk model to predict the prognosis of LUAD patients, which may aid in the classification of patients with LUAD and select medications for LUAD patients.
1 Introduction
Lung cancer is a global health issue and one of the top causes of morbidity and death among cancer patients, posing a severe danger to public health (Torre et al., 2012). There are several histological subtypes of lung cancer, with lung adenocarcinoma (LUAD) being the most common, accounting for around 40% of all lung malignancies (Ruiz-Cordero and Devine, 2020). The LUAD survival rate is just 4%–17% (Hirsch et al., 2017). In China, lung cancer prevention is complicated by high smoking rates and exposure to secondhand smoke (Qiu et al., 2021). Significant progress has been achieved in recent years as research into possible treatment targets for LUAD has continued to increase (Hirsch et al., 2017). Several immunotherapies targeting PD-1 and PD-L1 have been used clinically with significant survival benefits for LUAD patients (Steven et al., 2016). However, medication resistance and high recurrence rates continue to be the primary causes of treatment failure, resulting in unsatisfactory 5-year survival rates. To benefit more LUAD patients, there is an urgent need to identify new therapeutic targets and prognostic indicators to predict survival and guide clinical treatment in LUAD patients.
The presence of pathological angiogenesis is essential for tumor development and progression. On the one hand, these new blood vessels provide oxygen and nutrients to sustain rapid tumor growth and proliferation while helping tumor cells excrete metabolic waste. On the other hand, they provide pathways for tumor cells to enter the bloodstream and undergo distant metastasis (Viallard and Larrivée, 2017). Therefore, it is crucial for tumor cells’ survival, invasion, and metastasis (Jiang et al., 2020). Numerous activating and inhibitory factors govern tumor angiogenesis, control pathological angiogenesis, and impact patient prognosis. The factors vascular endothelial growth factor (Liu et al., 2021; Abdulazeem et al., 2022), platelet-derived growth factor (Zhang et al., 2019), and fibroblast-derived growth factor all play essential roles in angiogenesis (Eguchi and Wakabayashi, 2020). Serum VEGF promotes proliferation and migration, inhibits apoptosis, and regulates endothelial permeability (Ferrara et al., 2003). Anti-tumor angiogenic drugs such as bevacizumab (Syrigos et al., 2021) are currently used to treat lung adenocarcinoma. However, current research has focused on the impact of individual angiogenic genes on lung adenocarcinoma development and prognosis as a potential target for drug development. Few studies have integrated multiple angiogenesis-related genes by high-throughput biomarker sequencing and synthesized the relationship between these genes and lung adenocarcinoma prognosis and survival.
This work used TCGA’s comprehensive genome-wide gene expression profile to identify angiogenesis-related genes (AGRs) strongly correlated with lung cancer prognosis. We constructed and validated a diagnostic, prognostic, and recurrence model of lung adenocarcinoma and the corresponding nomogram. Corresponding data from GEO further validate these results. To summarize the flow chart of this study is shown in Figure 1.
FIGURE 1. Data gathering and processing flowchart.
2 Materials and methods
2.1 Selection of angiogenesis-related genes
AGRs were found on the GeneCards website (https://www.genecards.org/) by searching for “angiogenesis.” Correlation scores were used to indicate the strength of the correlation between genes and angiogenic activity, ranging from 0 to 100. High scores indicate strong correlations. For further analysis, we chose ARGs with correlation coefficients more prominent than 5. We then downloaded 36 additional tumor angiogenesis genes (hallmark-angiogenesis (Subramanian et al., 2005)) from Gene Enrichment Analysis (GSEA). By eliminating overlapping genes, a total of 137 angiogenesis-related genes were found.
2.2 Acquisition of lung adenocarcinoma datasets
The TCGA database (https://portal.gdc.cancer.gov/) was used to get the mRNA, mutation, and clinical data. There are 594 samples in the mRNA files (535 LUAD and 59 non-tumor samples). GSE68465, GSE50081 and GSEGSE31210 are validation datasets from the Gene Expression Omnibus Database (GEO, https://www.ncbi.nlm.nih.gov/geo/).
2.3 Differentially expressed gene screening and gene mutation analysis
To find angiogenic genes associated with LUAD, RNA sequencing data from the TCGA database was matched to a total of 137 ARGs. We then used the R package “limma” to identify differentially expressed genes (DEGs) with an absolute fold change (|logFC|) > 1 and an adjusted p-value <0.05. Using the “pheatmap” package, a heatmap was made to visualize differentially expressed genes. The protein-protein interaction (PPI) network is available on the STRING website (http://string.embl.de/) (Szklarczyk et al., 2017).
2.4 Functional enrichment analysis
Using the “clusterProfiler” package, KEGG (Kyoto Encyclopedia of Genes and Genomes) and GO (Gene Ontology) analyses were made to determine which pathways and functions angiogenesis-related genes are enriched in. DO (disease ontology) studies were also conducted using the R package “DOSE” to determine if angiogenesis-related genes are implicated in lung disorders, particularly LUAD (Yu et al., 2015).
2.5 Construction and validation of a prognostic signature based on angiogenesis-related genes
On a 7:3 ratio, 503 samples were separated into training and validation sets, with no samples designated as controls. Combining ARGs mRNA expression levels and clinical data, Univariate Cox regression analysis was used to find differentially expressed ARGs with significant prognostic significance. A predictive polygenic model was created using identified genes related to overall survival (OS). A minimum absolute shrinkage and selection operator (LASSO) Cox regression approach was used to construct multivariate models of ARG using the “glmnet” package of R software (Sauerbrei et al., 2007). Only genes with non-zero coefficients were chosen in LASSO regression to construct risk scores (Kidd et al., 2018). The best model was determined by maximizing the performance and using the least number of genes.
Subsequently, a predictive risk score formula was established based on a linear combination of expression levels and weighted regression coefficients obtained by LASSO Cox regression analysis.
2.6 Tumor mutation burden and gene mutation analysis
Somatic mutations in the TCGA were studied using the R package “maftools” to identify variations in somatic mutations between the HR and LR groups (Mayakonda et al., 2018). Subsequently, we estimated the two groups’ tumor mutation burden (TMB) per patient. Mutations in signature genes were searched on the cBioPortal website (https://www.cbioportal.org).
2.7 Immune microenvironment analysis
Immune cell infiltration was identified using timer 2.0 (cistrome.shinyapps.io/timer/) via the MCPCOUNTER, CIBERSORT, QUANTISEQ, Timer, CIBERSORT-ABS, EPIC, and XCELL algorithms. Infiltration levels of stromal and immune cells can be calculated with the ESTIMATE algorithm (Yoshihara et al., 2013). The enrichment score of genes in a particular gene set can be calculated by ssGSEA (single-sample gene set enrichment analysis). The process of ssGSEA includes ranking genes according to the absolute expression of genes in the sample and then calculating the enrichment score by integrating the differences between the empirical cumulative distribution functions of gene ranking (Bindea et al., 2013; Finotello and Trajanoski, 2018). Concentration scores of 16 immune cells were calculated using “GSEABase” and “GSVA” packages. The TIMER database studied six immune cells (B cells, macrophages, neutrophils, dendritic cells, CD8+ T cells, and CD4+ T cells) infiltration for its association with gene expression. Gene copy number variation in the TIMER database was studied for its potential impact on immune cell infiltration.
The expression of multiple immune checkpoint molecules was compared to determine whether there were differences in immune checkpoint blockade (ICB) therapy between the HR and LR groups. Immune checkpoints with differential expression between the two groups were visualized. Additionally, TIDE (Tumor Immune Dysfunction and Exclusion) score was calculated online following the instructions (https://tide.dfci.harvard.edu/). An inverse correlation was found between the TIDE score and ICB treatment success (Jiang et al., 2018).
2.8 Antigen presentation analysis
Human leukocyte antigen (HLA), found on numerous immune cells’ surfaces, is crucial for triggering cellular and humoral immunity (Halpert et al., 2020). To determine whether or not there were distinctions in antigen expression between the two groups, the “limma” package was used to compare the HLA expression levels of the two groups.
2.9 Cancer stem cell infiltration analysis
The UCSC Xena browser (http://xena.ucsc.edu/) was used to extract the DNA methylation-based stemness scores (DNAss) and RNA-based stemness scores (RNAss) of TCGA-LUAD patients. We conducted a comparative study at the DNA and RNA levels to examine the variations in stem cell infiltration between the two groups.
2.10 Gene set enrichment analysis (GSEA)
The c2. cp.kegg.v7.5.1. symbols and c5. go.v7.5.1. symbols collection was used to explore the function annotation in HR and LR groups using the R package “org. Hs. eg. db”. Gene sets with FDR <0.05 were considered statistically significant.
2.11 Predicting drug therapeutic response
Using the Cancer Immunome Atlas (https://tcia.at/), immunophenoscore (IPS) was derived to predict the sensitivity to immunotherapy. The IC50 of common chemotherapeutic agents in the total TCGA cohort was calculated by the administered “pRRophetic” software package to assess the predictive ability of AGRs for drug treatment response. The IC50 differences between the HR and LR groups were then compared using the Wilcoxon ranking test. Finally, we applied the “ggplot” R package to draw the results into bar charts.
2.12 Verification of signature genes in databases
The Cancer Cell Line Encyclopedia (CCLE) Line (http://www.broadinstitute.org/ccle) provides downloadable information on the mRNA expression levels of several LUAD cell lines. The Human Protein Atlas (HPA) database was used to verify protein expression levels.
2.13 Quantitative real-time PCR (qPCR)
The expression of model genes in humans was confirmed using a quantitative real-time PCR technique. Total RNA was isolated from human tissues initially using the Trizol reagent. Then, reverse transcription was used to transform the isolated RNA into cDNAs. Finally, the quantitative real-time PCR technique was used to assess the expression levels of model genes in human tissues. The primers are shown in Supplementary Table S3. Relative expression values of model genes were calculated using the 2−ΔΔCT method and normalized with beta-actin. The experiments were reviewed and approved by the authority (Ethics No: WDRY2022-K041) and executed in accordance with relevant guidelines.
2.14 Statistical analysis
R was used for all statistical analysis and graphical creation (version 4.2.1). Volcano plots were drawn with the “ggplot2″ package. Violin diagrams are drawn with “ggpubr” packets. The Mann-Whitney test was utilized for differentially expressed gene analysis, tumor mutation burden analysis, ssGSEA score, immunological checkpoint analysis, and HLA analysis. Cancer stem cell infiltration and drug sensitivity tests were conducted using a correlation test. The log-rank test and Kaplan-Meier analysis compared overall survival (OS) between groups.
3 Results
3.1 Identification and exploration of angiogenesis-related differentially expressed genes
Differential expression analysis found fifty-six genes to be expressed differently between 539 cancers and 59 normal samples (Figure 2A). Twenty-five genes were elevated, and 31 were downregulated in tumor samples relative to normal samples (Figure 2B). PPI revealed that the majority of genes are linked (Figure 2C). According GO enrichment analyses, different ARGs are crucial for LUAD angiogenesis and vascular growth (Figure 3A). KEGG pathway enrichment analysis revealed that differentially expressed ARGs were primarily implicated in tumor angiogenesis pathways such as PI3K-Akt, MAPK, and Rap1 (Figure 3B). DO research revealed that they were linked to lung illness, lung adenocarcinoma, non-small cell lung cancer, and other conditions (Figure 3C, Supplementary Table S1).
FIGURE 2. A sum of 56 genes showed significant differential expression. (A) Heatmaps of differentially expressed angiogenesis-related genes. (B) Volcano plots of differentially expressed angiogenesis-related genes. (C) PPI of 56 genes.
FIGURE 3. Functional enrichment analysis of differential genes. (A) GO. (B) KEGG. (C)DO.
3.2 Developing an RS prediction model for angiogenesis-related genes
First, we took the intersection of the total TCGA cohort with GSE68465 to obtain the common genes. Subsequently, the two cohorts intersected with differentially expressed genes, and 51 genes were obtained for subsequent analysis. Thirteen of the 51 differentially expressed genes correlate with the overall survival of LUAD patients in the TCGA cohort (Supplementary Figures S1A–N). Then ten genes in these 13 genes were finally screened by Lasso regression (Figures 4A, B). Finally, the multivariate cox regression revealed six of these ten angiogenic-associated genes (MET, PDGFB, TIMP1, PECAM1, CCND2, POSTN) as independent risk variables to build risk models, the angiogenesis-related genes (ARGs) risk model. And the ARGs risk score formula was as follows:
FIGURE 4. Establishment of an ARGs prognostic model about LUAD’s prognosis using a lasso regression model. (A) A LASSO coefficient profile of the 13 ARGs. (B) A coefficient profile plot was generated against the log(lambda) sequence. (C) Univariate COX regression analysis for RS in TCGA training cohort patients with LUAD. (D) Multivariate Cox regression analysis for RS in TCGA training cohort patients with LUAD.
FIGURE 5. Development of RS based on the 6 ARGs signature of patients with LUAD in the TCGA training and TCGA validation cohorts. (A, B) The heatmap of the 6 ARGs expression profiles between the high-risk and low-risk groups in the training or validation cohort, together with the RS distribution and patient vital status. (C, D) Kaplan-Meier analysis of the prognostic model in training or validation cohort. (E, F) The ideal AUC of the gene signature in the two cohorts was determined using a time-dependent ROC analysis.
FIGURE 6. Construction of a predictive nomogram. (A, B) The nomogram for OS prediction in LUAD patients at 1, 3, and 5 years. (C, D) Calibration curves of the nomogram for OS prediction at 1, 3, and 5 years.
FIGURE 7. Subgroup study of various clinical characteristics using the RS formula. (A) Age>65 years and age≤65 years. (B) FEMALE and MALE. (C) M0 and M1. (D) NO and N1-3. (E) Stage I-II and stage III-IV. (F) T1-2 and T 3-4.
3.3 Tumor mutation burden and gene mutation analysis
Given that TMB is linked to immunotherapy effectiveness, we calculated the difference in TMB value between the two groups. TMB was considerably more remarkable in the HR group than in the LR group (Figure 8B). However, those with elevated TMB had a favorable survival advantage (Figure 8C). The utility of integrating risk scores with TMB to predict patient outcomes was then investigated. The Kaplan-Meier analysis suggested that LR scores and high TMB were associated with more prolonged survival (Figure 8D). We analyzed mutation rates between the HR and LR groups. The findings revealed a higher frequency of mutational events within the HR group. The most prevalent alterations in the HR group, TP53 and TTN were also considerably more frequent than those in the LR group. The other six LUAD mutated genes (MUC16, CSMD3, RYR2, LRP1B, ZFHX4, and USH2A) also showed increasing trends to different degrees (Figures 8E, F). Using the LUAD cohort from the cBioportal database, we analyzed the value of six angiogenic genes (MET, PDGFB, TIMP1, PECAM1, CCND2, POSTN) in the development of diagnostic, prognostic, and recurrence models. The study found that 19.6% (566 LUAD samples) showed genetic alterations, of which 7% occurred in MET, whose primary alteration was amplification, and 6% of mutations occurred in POSTN (Figure 8A). These genes may be used as therapeutic targets, offering fresh avenues for treating LUAD.
FIGURE 8. Tumor mutation burden and gene mutation analysis. (A) Genetic alteration analysis of the ARGs in the TCGA LUAD cohort. (B) TMB between the low- and high-risk subgroups based on RS. (C) Survival study of the various TMB stratification groups. (D) Analysis of the survival of various groups stratified by TMB and RS. (E, F) Assessment of the differences in the mutational landscape between high- and low-risk.
3.4 Immune features analysis
We investigated the relationship between RS and the immune status of patients in the TCGA cohort and found a significant change in immune cells. From the results, the LR group had a better immune status (Figure 9I). While there was a negative correlation between the ESTIMATE score of immune cells and the risk score, there was a positive correlation between the ESTIMATE score of stromal cells and the risk score (Figures 9A, B). Next, we compared the immune cell scores of the two patient groups. The total TCGA cohort showed that the LR group’s infiltration scores of most immune cells were higher than in HR groups, such as aDC, B cells, and T-helper cells (Figure 9G). Considering the importance of checkpoint inhibitors in clinical treatment, we further analyzed the differences in ICBs expression and found substantial differences in CTLA4, CD28, ID02, and CD27 between the two groups (Figure 9H). We next evaluated the potential therapeutic effectiveness of immunotherapy in both patient groups using TIDE. A higher TIDE prediction score, the higher the likelihood of immune evasion, suggesting that patients are less likely to benefit from ICI therapy. We found that patients in the HR group had a lower TIDE score than those in the LR group, suggesting that the HR patients might respond better to ICI therapy (Figure 9C). In addition, we discovered that CAF and T-cell exclusion ratings were more excellent in the HR group, whereas T-cell dysfunction scores were higher in the LR group (Figures 9D–F).
FIGURE 9. Immune features analysis. (A) Correlation analysis of the risk score and immune score. (B) Correlation analysis of the risk score and stromal score. (C) Comparison of TIDE and risk score. (D) Comparison of CAF and risk score. (E) Comparison of T-cell exclusion and risk score. (F) Comparison of T-cell dysfunction and risk score. (G) Infiltration score of immune cells. (H) Expression of immune checkpoints. (I) Heatmap of immune response in high and low-risk groups based on multiple algorithms. *p < 0.05, **p < 0.01, ***p < 0.001.
Six signature genes were revealed to be linked with immune cell infiltration. Positive correlations were seen between MET and CD4+ T cells, macrophages, dendritic cells, CD8+ T cells, and neutrophils, but not B cells (Supplementary Figure S4B). PDGFB was associated favorably with the other four immune cells and negatively with B cells, but not with CD8+ T cells (Supplementary Figure S4C). POSTN was negatively connected to B cells and favorably related to the other four immune cells, but not to CD4+ T cells (Supplementary Figure S4E). Positive correlations were seen between TIMP1 and CD4+ T cells, macrophages, dendritic cells, and neutrophils, but not CD8+ T cells and B cells (Supplementary Figure S4F). Positive correlations were seen between CCND2 and PECAM1 and all six immune cells (Supplementary Figures S4A, D). Copy number variation in genes may affect immune cell infiltration compared to the diploid/normal group. As PECAM1 was not available in the database, we analyzed the copy number variation of the remaining 5 cell types about immune cells. MET copy number variants were significantly associated with decreased levels of three immune cell types (B cells, neutrophils, and macrophages), POSTN copy number variants were significantly associated with decreased levels of three immune cell types (B cells, CD4 + T cells, and macrophages), TIMP1 copy number variants were significantly associated with decreased levels of six immune cell types (B cells, CD4 + T cells, CD8 + T cells, neutrophils, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and dendritic cells), PDGFB copy number variation was significantly associated with decreased levels of three immune cells (B cells, CD4+ T cells, and macrophages), and CCND2 copy number variation was significantly associated with decreased levels of five immune cells (B cells, CD4+ T cells, neutrophils, macrophages, and dendritic cells) (Figures 10A–E). These results suggest that changes in model gene copy number can significantly modulate immune cell infiltration.
FIGURE 10. Association between CNV and immunity in 5 model genes. (A) CCND2. (B) MET. (C) PDGFB. (D) TIMP1. (E) POSTN. *p < 0.05, **p < 0.01, ***p < 0.001.
3.5 Cancer stem cell infiltration analysis
We used DNAss and RNAss to compare tumor stemness amongst distinct risk patterns. The results demonstrated a positive relationship between risk scores and DNAss (r = 0.11, p = 0.019) and RNAss (r = 0.15, p = 0.0017) values (Figures 11C, D). This shows that the HR group’s LUAD cells exhibited a more marked stem cell profile and a lower degree of cell differentiation.
FIGURE 11. Exploration of possible reasons for differences between high and low-risk groups. (A, B) KEGG Enrichment Analysis. (C, D) Tumor stem cell analysis. (E, F) GO Enrichment Analysis. (G) Expression of HLA. *p < 0.05, **p < 0.01, ***p < 0.001.
3.6 Gene set enrichment analyses
According to KEGG, pathways enriched in the HR group included the cell cycle, DNA replication, proteasome, and focal adhesion, whereas pathways enriched in the LR group included asthma, linoleic acid metabolism, systemic lupus erythematosus, and alpha-linolenic acid metabolism (Figures 11A, B). Additionally, GO revealed that genes in the LR group were engaged in complement activation, antigen binding, immunoglobulin complexes, and other processes, whereas genes in the HR group were enriched in chromosome segregation, condensed chromosomes, and nuclear chromosome segregation (Figures 11E, F).
3.7 Antigen presentation analysis
There was a substantial variation in HLA expression associated with antigen presentation between the HR and LR groups. The expression of numerous HLA classes I and II was more significant in the LR group than in the HR group in the total TCGA cohort (Figure 11G).
3.8 Drug sensitivity analysis
Through TCIA, we specifically looked at the impact of risk scores on the efficacy of immunotherapy. The results demonstrated that the HR group was more likely to respond to CTLA4-positive/PD-L marker-positive treatment than the LR group (Figures 12C, D). This shows that patients in the HR group may respond better to CTLA4-positive/PD-L marker-positive immunotherapy, leading to a better clinical outcome. Patients in the HR group were more responsive to docetaxel, erlotinib, etoposide, gemcitabine, and cisplatin, as evidenced by their lower IC50s compared to those in the LR group (Figures 12H–L). The IC50 of ABT.888 (also known as veliparib), axitinib, and ATRA was higher in the HR group than in the LR group (Figures 12E–G). We examined the association between six signature genes and medication effectiveness in NCI-60 cell lines, including the LUAD cell line, to guide gene-targeted treatment. The findings of the first sixteen studies were shown based on the order of the p-values, from smallest to most lavish. PECAM1 was sensitive to bendamustine and methylprednisolone, and zalcitabine, MET was sensitive to auranofin, lomustine, and arsenic trioxide, and CCND2 was sensitive to amiodarone hydrochloride (Figure 12M). POSTN was sensitive to zoledronic acid and caffeic acid (Figure 12M). Supplementary Table S2 also lists the names of the medicines and genes along with their connection coefficients.
FIGURE 12. Immunotherapeutic responses of the ARGs prognostic signature and drug sensitivity analysis. (A–D) Comparison of the immunophenoscore (IPS) between the low- and high-risk groups stratified by both CTLA4 and PD-1. (E–L) Drug sensitivity analysis. (M) Multiple drugs for six signature genes.
3.9 Verification of signature genes
The IHC data for six proteins in the HPA database were evaluated. Despite the inability to precisely assess the differences between normal lung tissue and malignancies, the preliminary findings suggest increased expression of CCND2 and PDGFB in tumor tissues (Figure 13G). Meanwhile, PECAM1 and POSTN levels may be lowered in tumorous tissues (Figure 13G). There was no discernible variation in the levels of MET and TIMP1 expression (Figure 13G). In several LUAD cell lines, the same gene is expressed differently (Figures 13A–F). We employed a quantitative real-time PCR technique to measure the expression levels of model genes in human tissues to confirm our prognostic model’s validity further. The findings revealed that MET, PDGFB, PECAM1, POSTN, and TIMP1 expression differed in human tissues and were consistent with the TCGA database (Supplementary Figures S6A–E). In human tissues, however, there was no variation in CCND2 expression (Supplementary Figure S6F).
FIGURE 13. HPA database and compare gene expression in LUAD cell lines. (A) CCND2. (B) MET. (C) PDGFB. (D)PECAM1. (E) POSTN. (F) TIMP1. (G) Immunohistochemical results of 6 model genes in the HPA database.
4 Discussion
The vast majority of cases of NSCLC are LUAD (Bade and Dela Cruz, 2020). LUAD has a poor prognosis due to a lack of knowledge of its pathogenesis (Peng et al., 2021). Targeted treatment and immunotherapy for LUAD have made significant strides in recent years due to researchers examining the etiology of LUAD from various routes and biological processes. Despite this progress, it remains difficult for doctors to assess the risk of LUAD and make prognostic predictions for their patients. Since the numerous scoring methods now used in clinical practice, such as the TNM staging system, have more or fewer limitations, there is an urgent need for novel molecular biomarkers to predict the survival of LUAD patients. Angiogenesis plays a crucial role in tumor formation; hence, we developed and validated a risk model to predict the prognosis of LUAD patients based on angiogenesis-related genes in this study. This study is the first to our knowledge to develop predictive indicators using a collection of angiogenesis-related genes in patients with LUAD.
Angiogenesis is essential for tissue repair, reproduction, embryonic development, and wound healing, among other physiological functions (Cui et al., 2014). In the normal resting state, endothelial cells can sense angiogenic signals and remain highly plastic under controlled conditions to participate in angiogenesis. However, in many disease states, such as cancer, rheumatoid arthritis, and atherosclerosis, uncontrolled angiogenesis will further promote disease progression and become a hallmark of these disease states (Carmeliet and Jain, 2011). Metastasis and growth of LUAD are similarly dependent on the process of angiogenesis. In hypoxic and energy-deficient conditions, tumor cells and their surroundings produce significant quantities of angiogenic factors, promoting pathological angiogenesis. These new blood vessels will offer oxygen and energy to the tumor and carry metabolic waste generated by the tumor cells, boosting tumor development.
In this study, we collected transcript expression data and corresponding clinical data from TCGA and GEO databases. We determined the differential expression of ARGs in tumor and normal tissues using the TCGA database. GO annotation and KEGG enrichment analysis revealed that these differential genes are primarily engaged in angiogenesis and vascular development and are enriched in the P13K-Akt, Rap1, and MAPK pathways. PI3K/AKT is an essential intracellular signal transduction molecule involved in the control of cell proliferation, apoptosis, and differentiation, and it can affect the production of VEGF and hypoxia-inducible factor (HIF-1) through activation of p70S6K1 and Hdm2 (Skinner et al., 2004). A vascular endothelial growth factor is an essential regulator in tumor angiogenesis, and hypoxia increases the release of angiogenic factors. Furthermore, the PI3K-Akt pathway is critical in K-ras signaling pathway-mediated hematopoiesis and angiogenesis (Liu et al., 2008). The MAPK signaling pathway is one of the essential signaling systems in living creatures, and it is crucial for cell survival, proliferation, and angiogenesis (Teng et al., 2018). Numerous studies have demonstrated that RAP1 is activated in numerous malignancies, including leukemia and solid tumors (Minato, 2013; Shah et al., 2019). aRAP1 plays a role in the invasion and metastasis of various tumor cells by regulating adhesive junctions and cytoskeletal remodeling. DO analysis has shown that these differential genes are associated with lung adenocarcinoma, non-small cell lung cancer, and lung disease.
Subsequently, we used 6 ARGs (MET, CCND2, PDGFB, POSTN, PECAM1, TIMP1) to establish model equations for risk assessment. TIMP1 is a natural inhibitor of matrix metalloproteinases, and the imbalance between TIMP1 and matrix metalloproteinases in gastrointestinal cancers is a crucial element in colon cancer (Song et al., 2016). The interaction of TIMP1 protein’s C-terminal structural domain with tetraspanin CD63 stimulates conformational activation of integrin b1 and activates MAPK signaling, resulting in cancer (Jung et al., 2006). Through the TIMP-1/CD63 signal, TIMP1 can activate fibroblast-like hepatic stellate cells (HSCs) and release SDF-1 to attract neutrophils that promote metastasis. In this process, TIMP1 significantly increases the sensitivity of the liver to circulating tumor cells and creates a tumor microenvironment that promotes tumor liver metastasis (Seubert et al., 2015; Grünwald et al., 2016). A transmembrane receptor tyrosine kinase, MET is the hepatocyte growth factor receptor. When it binds to ligands, it activates several downstream channels that govern cell growth, survival, and migration, among other things (Trusolino et al., 2010). It has been demonstrated that MET interacts with vascular endothelial growth factors to promote neovascularization (Eder et al., 2009). POSTN promotes angiogenesis in colorectal malignancies by preventing stress-induced apoptosis and increasing endothelial cell survival, thereby promoting colon cancer metastasis, leading to poor prognosis (Bao et al., 2004). In some tumors, such as glioma (Hermanson et al., 1992; Guo et al., 2003) and ovarian cancer (Czekierdowska et al., 2017), PDGFB plays a crucial role in promoting angiogenesis and stimulating malignant cell proliferation. PECAM-1 (also referred to as a cluster of differentiation 31, CD31) is primarily regarded as an adhesion molecule that participates in cell proliferation, apoptosis, migration, and cellular immunity. PECAM-1 is expressed in some tumor cells and has been linked to tumor invasion (Ranamukhaarachchi et al., 2019; Gong et al., 2020). CCND2 encodes cell cycle protein D2, which regulates cell cycle protein-dependent kinases 4 and 6 (CDK) 4/6 in G1-S (Ortega et al., 2002). Although CCND2 dysregulation is a significant source of medication resistance in breast cancer endocrine treatment (Kwapisz, 2017), its significance in LUAD remains unknown. The role of these genes in LUAD needs to be further investigated. Subsequently, through a series of validations and exploration of possible related mechanisms, we demonstrated that the risk model we constructed could predict the prognosis of LUAD patients.
Our study certainly has some limitations. To some extent, intra-tumor or intra-patient tumor heterogeneity is unavoidable, given that our study cohort was compiled from a wide variety of high-throughput sequencing platforms and public databases. Several studies demonstrate that the heterogeneity of tumors impacts the effectiveness of immunotherapy and chemotherapy. We must disregard the considerable heterogeneity of lung adenocarcinoma due to data constraints. Second, although it has been discovered that angiogenesis-related genes influence immunological interactions and survival in LUAD patients, the scientific or medical processes behind these occurrences are not fully understood. To validate and explain the significance of angiogenesis-related genes in LUAD, therefore, large-scale prospective research and functional and mechanistic tests are required. Third, while the median RS cut-off point was utilized to categorize LUAD samples into HR and HR groups, the ideal RS cut-off point may be a preferred technique for stratifying LUAD patients. Due to the lack of comprehensive clinicopathological data, we finally compiled and adjusted certain clinical data for survival analysis and Cox regression. However, this may introduce bias and uncertainty regarding whether RS is an independent predictive factor.
Non-etheless, our research identifies genes and pathways implicated in LUAD angiogenesis and develops a valid predictive model that correlates with immune infiltration features and treatment. The RS model contributes to a greater understanding of the prognosis of LUAD patients and opens up new perspectives for targeted therapy.
Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.
Author contributions
MP and JQ were responsible for the design of the entire study and the revision of the manuscript. JiW were responsible for the analysis of the data and preparation of the manuscript. XC, YmW, and JaW are responsible for data collection. XC, PW, and TW were responsible for some of the pictures and tables. MP, JQ, and JiW revised the manuscript. All authors had read and revised the final version.
Funding
This study was supported by the national natural science foundation of China (82103628), Wu Jieping Medical Foundation Project (320.6750.2021-02-127), Beijing Science and Innovation Medical Development Foundation (KC2021-JX-0186-18), National Science Foundation of China (No. 81770169), and China primary healthcare foundation (cphcf-2022-183).
Acknowledgments
The authors gratefully acknowledge the data generated by GEO and TCGA databases used in this study.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2023.1092968/full#supplementary-material
References
Abdulazeem, L., Tariq, A., and Abdalkareem Jasim, S. (2022). An investigation of vascular endothelial growth factor (VEGFR-1 and VEGFR-2) in burn wound healing. Archives Razi Inst. 77 (2), 747–751. doi:10.22092/ARI.2022.356981.1954
Bade, B. C., and Dela Cruz, C. S. (2020). Lung cancer 2020: Epidemiology, etiology, and prevention. Clin. chest Med. 41 (1), 1–24. doi:10.1016/j.ccm.2019.10.001
Bao, S., Ouyang, G., Bai, X., Huang, Z., Ma, C., Liu, M., et al. (2004). Periostin potently promotes metastatic growth of colon cancer by augmenting cell survival via the Akt/PKB pathway. Cancer Cell 5 (4), 329–339. doi:10.1016/s1535-6108(04)00081-9
Bindea, G., Mlecnik, B., Tosolini, M., Kirilovsky, A., Waldner, M., Obenauf, A. C., et al. (2013). Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer. Immunity 39 (4), 782–795. doi:10.1016/j.immuni.2013.10.003
Carmeliet, P., and Jain, R. K. (2011). Molecular mechanisms and clinical applications of angiogenesis. Nature 473 (7347), 298–307. doi:10.1038/nature10144
Cui, H., Wang, Y., Huang, H., Yu, W., Bai, M., Zhang, L., et al. (2014). GPR126 protein regulates developmental and pathological angiogenesis through modulation of VEGFR2 receptor signaling. J. Biol. Chem. 289 (50), 34871–34885. doi:10.1074/jbc.M114.571000
Czekierdowska, S., Stachowicz, N., Chróściel, M., and Czekierdowski, A. (2017). Proliferation and maturation of intratumoral blood vessels in women with malignant ovarian tumors assessed with cancer stem cells marker nestin and platelet derived growth factor PDGF-B. Ginekol. Pol. 88 (3), 120–128. doi:10.5603/GP.a2017.0023
Duan, J., Lei, Y., Lv, G., Liu, Y., Zhao, W., Yang, Q., et al. (2021). Identification of a novel autophagy signature for predicting survival in patients with lung adenocarcinoma. PeerJ 9, e11074. doi:10.7717/peerj.11074
Eder, J. P., Vande Woude, G. F., Boerner, S. A., and LoRusso, P. M. (2009). Novel therapeutic inhibitors of the c-Met signaling pathway in cancer. Clin. cancer Res. 15 (7), 2207–2214. doi:10.1158/1078-0432.CCR-08-1306
Eguchi, R., and Wakabayashi, I. (2020). HDGF enhances VEGF-dependent angiogenesis and FGF-2 is a VEGF-independent angiogenic factor in non-small cell lung cancer. Oncol. Rep. 44 (1), 14–28. doi:10.3892/or.2020.7580
Ferrara, N., Gerber, H. P., and LeCouter, J. (2003). The biology of VEGF and its receptors. Nat. Med. 9 (6), 669–676. doi:10.1038/nm0603-669
Finotello, F., and Trajanoski, Z. (2018). Quantifying tumor-infiltrating immune cells from transcriptomics data. CII 67 (7), 1031–1040. doi:10.1007/s00262-018-2150-z
Gong, G., Lin, T., and Yuan, Y. (2020). Integrated analysis of gene expression and DNA methylation profiles in ovarian cancer. J. ovarian Res. 13 (1), 30. doi:10.1186/s13048-020-00632-9
Gong, Z., Li, Q., Yang, J., Zhang, P., Sun, W., Ren, Q., et al. (2022). Identification of a pyroptosis-related gene signature for predicting the immune status and prognosis in lung adenocarcinoma. Front. Bioeng. Biotechnol. 10, 852734. doi:10.3389/fbioe.2022.852734
Grünwald, B., Harant, V., Schaten, S., Frühschütz, M., Spallek, R., Höchst, B., et al. (2016). Pancreatic premalignant lesions secrete tissue inhibitor of metalloproteinases-1, which activates hepatic stellate cells via CD63 signaling to create a premetastatic niche in the liver. Gastroenterology 151 (5), 1011–1024. e7. doi:10.1053/j.gastro.2016.07.043
Guo, P., Hu, B., Gu, W., Xu, L., Wang, D., Huang, H. J., et al. (2003). Platelet-derived growth factor-B enhances glioma angiogenesis by stimulating vascular endothelial growth factor expression in tumor endothelia and by promoting pericyte recruitment. Am. J. pathology 162 (4), 1083–1093. doi:10.1016/S0002-9440(10)63905-3
Halpert, M. M., Konduri, V., Liang, D., Vazquez-Perez, J., Hofferek, C. J., Weldon, S. A., et al. (2020). MHC class I and II peptide homology regulates the cellular immune response. FASEB J. 34 (6), 8082–8101. doi:10.1096/fj.201903002R
Hermanson, M., Funa, K., Hartman, M., Claesson-Welsh, L., Heldin, C. H., Westermark, B., et al. (1992). Platelet-derived growth factor and its receptors in human glioma tissue: Expression of messenger RNA and protein suggests the presence of autocrine and paracrine loops. Cancer Res. 52 (11), 3213–3219.
Hirsch, F. R., Scagliotti, G. V., Mulshine, J. L., Kwon, R., Curran, W. J., Wu, Y. L., et al. (2017). Lung cancer: Current therapies and new targeted treatments. Lancet (London, Engl. 389 (10066), 299–311. doi:10.1016/S0140-6736(16)30958-8
Jiang, P., Gu, S., Pan, D., Fu, J., Sahu, A., Hu, X., et al. (2018). Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat. Med. 24 (10), 1550–1558. doi:10.1038/s41591-018-0136-1
Jiang, X., Wang, J., Deng, X., Xiong, F., Zhang, S., Gong, Z., et al. (2020). The role of microenvironment in tumor angiogenesis. J. Exp. Clin. cancer Res. CR 39 (1), 204. doi:10.1186/s13046-020-01709-5
Jung, K. K., Liu, X. W., Chirco, R., Fridman, R., and Kim, H. R. (2006). Identification of CD63 as a tissue inhibitor of metalloproteinase-1 interacting cell surface protein. EMBO J. 25 (17), 3934–3942. doi:10.1038/sj.emboj.7601281
Kidd, A. C., McGettrick, M., Tsim, S., Halligan, D. L., Bylesjo, M., and Blyth, K. G. (2018). Survival prediction in mesothelioma using a scalable lasso regression model: Instructions for use and initial performance using clinical predictors. BMJ open Respir. Res. 5 (1), e000240. doi:10.1136/bmjresp-2017-000240
Kwapisz, D. (2017). Cyclin-dependent kinase 4/6 inhibitors in breast cancer: Palbociclib, ribociclib, and abemaciclib. Breast cancer Res. Treat. 166 (1), 41–54. doi:10.1007/s10549-017-4385-3
Liu, L., Zhu, S., Gong, Z., and Low, B. C. (2008). K-ras/PI3K-Akt signaling is essential for zebrafish hematopoiesis and angiogenesis. PloS one 3 (8), e2850. doi:10.1371/journal.pone.0002850
Liu, Y., Long, L., Zhang, F., Hu, X., Zhang, J., Hu, C., et al. (2021). Microneedle-mediated vascular endothelial growth factor delivery promotes angiogenesis and functional recovery after stroke. J. Control. release official J. Control. Release Soc. 338, 610–622. doi:10.1016/j.jconrel.2021.08.057
Mayakonda, A., Lin, D. C., Assenov, Y., Plass, C., and Koeffler, H. P. (2018). Maftools: Efficient and comprehensive analysis of somatic variants in cancer. Genome Res. 28 (11), 1747–1756. doi:10.1101/gr.239244.118
Minato, N. (2013). Rap G protein signal in normal and disordered lymphohematopoiesis. Exp. Cell Res. 319 (15), 2323–2328. doi:10.1016/j.yexcr.2013.04.009
Ortega, S., Malumbres, M., and Barbacid, M. (2002). Cyclin D-dependent kinases, INK4 inhibitors and cancer. Biochimica biophysica acta 1602 (1), 73–87. doi:10.1016/s0304-419x(02)00037-9
Peng, B., Li, H., Na, R., Lu, T., Li, Y., Zhao, J., et al. (2021). Identification of a novel prognostic signature of genome instability-related LncRNAs in early stage lung adenocarcinoma. Front. Cell Dev. Biol. 9, 706454. doi:10.3389/fcell.2021.706454
Qiu, H., Cao, S., and Xu, R. (2021). Cancer incidence, mortality, and burden in China: A time-trend analysis and comparison with the United States and United Kingdom based on the global epidemiological data released in 2020. Cancer Commun. Lond. Engl. 41 (10), 1037–1048. doi:10.1002/cac2.12197
Ranamukhaarachchi, S. K., Modi, R. N., Han, A., Velez, D. O., Kumar, A., Engler, A. J., et al. (2019). Macromolecular crowding tunes 3D collagen architecture and cell morphogenesis. Biomaterials Sci. 7 (2), 618–633. doi:10.1039/c8bm01188e
Ruiz-Cordero, R., and Devine, W. P. (2020). Targeted therapy and checkpoint immunotherapy in lung cancer. Surg. Pathol. Clin. 13 (1), 17–33. doi:10.1016/j.path.2019.11.002
Sauerbrei, W., Royston, P., and Binder, H. (2007). Selection of important variables and determination of functional form for continuous predictors in multivariable model building. Statistics Med. 26 (30), 5512–5528. doi:10.1002/sim.3148
Seubert, B., Grünwald, B., Kobuch, J., Cui, H., Schelter, F., Schaten, S., et al. (2015). Tissue inhibitor of metalloproteinases (TIMP)-1 creates a premetastatic niche in the liver through SDF-1/CXCR4-dependent neutrophil recruitment in mice. Hepatol. Baltim. Md) 61 (1), 238–248. doi:10.1002/hep.27378
Shah, S., Brock, E. J., Ji, K., and Mattingly, R. R. (2019). Ras and Rap1: A tale of two GTPases. Seminars cancer Biol. 54, 29–39. doi:10.1016/j.semcancer.2018.03.005
Skinner, H. D., Zheng, J. Z., Fang, J., Agani, F., and Jiang, B. H. (2004). Vascular endothelial growth factor transcriptional activation is mediated by hypoxia-inducible factor 1alpha, HDM2, and p70S6K1 in response to phosphatidylinositol 3-kinase/AKT signaling. J. Biol. Chem. 279 (44), 45643–45651. doi:10.1074/jbc.M404097200
Song, G., Xu, S., Zhang, H., Wang, Y., Xiao, C., Jiang, T., et al. (2016). TIMP1 is a prognostic marker for the progression and metastasis of colon cancer through FAK-PI3K/AKT and MAPK pathway. J. Exp. Clin. cancer Res. CR 35 (1), 148. doi:10.1186/s13046-016-0427-7
Steven, A., Fisher, S. A., and Robinson, B. W. (2016). Immunotherapy for lung cancer. Respirol. Carlt. Vic. 21 (5), 821–833. doi:10.1111/resp.12789
Subramanian, A., Tamayo, P., Mootha, V. K., Mukherjee, S., Ebert, B. L., Gillette, M. A., et al. (2005). Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. Sci. U. S. A. 102 (43), 15545–15550. doi:10.1073/pnas.0506580102
Syrigos, K., Abert, I., Andric, Z., Bondarenko, I. N., Dvorkin, M., Galic, K., et al. (2021). Efficacy and safety of bevacizumab biosimilar FKB238 versus originator bevacizumab: Results from AVANA, a phase III trial in patients with non-squamous non-small-cell lung cancer (non-sq-NSCLC). BioDrugs. 35 (4), 417–428. doi:10.1007/s40259-021-00489-4
Szklarczyk, D., Morris, J. H., Cook, H., Kuhn, M., Wyder, S., Simonovic, M., et al. (2017). The STRING database in 2017: Quality-controlled protein-protein association networks, made broadly accessible. Nucleic acids Res. 45 (D1), D362–d8. doi:10.1093/nar/gkw937
Teng, F., Xu, Z., Chen, J., Zheng, G., Zheng, G., Lv, H., et al. (2018). DUSP1 induces apatinib resistance by activating the MAPK pathway in gastric cancer. Oncol. Rep. 40 (3), 1203–1222. doi:10.3892/or.2018.6520
Torre, L. A., Bray, F., Siegel, R. L., Ferlay, J., Lortet-Tieulent, J., and Jemal, A. (2012). Global cancer statistics. CA a cancer J. Clin. 65 (2), 87–108. doi:10.3322/caac.21262
Trusolino, L., Bertotti, A., and Comoglio, P. M. (2010). MET signalling: Principles and functions in development, organ regeneration and cancer. Nat. Rev. Mol. Cell Biol. 11 (12), 834–848. doi:10.1038/nrm3012
Viallard, C., and Larrivée, B. (2017). Tumor angiogenesis and vascular normalization: Alternative therapeutic targets. Angiogenesis 20 (4), 409–426. doi:10.1007/s10456-017-9562-9
Xu, Q., and Chen, Y. (2021). An aging-related gene signature-based model for risk stratification and prognosis prediction in lung adenocarcinoma. Front. Cell Dev. Biol. 9, 685379. doi:10.3389/fcell.2021.685379
Yang, J., Chen, Z., Gong, Z., Li, Q., Ding, H., Cui, Y., et al. (2022). Immune landscape and classification in lung adenocarcinoma based on a novel cell cycle checkpoints related signature for predicting prognosis and therapeutic response. Front. Genet. 13, 908104. doi:10.3389/fgene.2022.908104
Yoshihara, K., Shahmoradgoli, M., Martínez, E., Vegesna, R., Kim, H., Torres-Garcia, W., et al. (2013). Inferring tumour purity and stromal and immune cell admixture from expression data. Nat. Commun. 4, 2612. doi:10.1038/ncomms3612
Yu, G., Wang, L. G., Yan, G. R., and He, Q. Y. (2015). Dose: An R/bioconductor package for disease ontology semantic and enrichment analysis. Bioinforma. Oxf. Engl. 31 (4), 608–609. doi:10.1093/bioinformatics/btu684
Zhang, J., Zhang, H., Chen, Y., Fu, J., Lei, Y., Sun, J., et al. (2019). Platelet-derived growth factor D promotes the angiogenic capacity of endothelial progenitor cells. Mol. Med. Rep. 19 (1), 125–132. doi:10.3892/mmr.2018.9692
Zhang, W., Wan, S., Qu, Z., Ge, J., Zhang, C., Li, C., et al. (2022). Establishment of a prognostic signature for lung adenocarcinoma by integration of 7 pyroptosis-related genes and cross-validation between the TCGA and GEO cohorts: A comprehensive bioinformatics analysis. Medicine 101 (29), e29710. doi:10.1097/MD.0000000000029710
Keywords: angiogenesis-related genes, risk model, lung adenocarcinoma, prognosis, immunotherapy
Citation: Wang J, Cui X, Weng Y, Wei J, Chen X, Wang P, Wang T, Qin J and Peng M (2023) Application of an angiogenesis-related genes risk model in lung adenocarcinoma prognosis and immunotherapy. Front. Genet. 14:1092968. doi: 10.3389/fgene.2023.1092968
Received: 08 November 2022; Accepted: 09 January 2023;
Published: 01 February 2023.
Edited by:
Apostolos Zaravinos, European University Cyprus, CyprusReviewed by:
Hua Huang, Tianjin Medical University, ChinaGuangsheng Zhu, Tianjin Medical University General Hospital, China
Copyright © 2023 Wang, Cui, Weng, Wei, Chen, Wang, Wang, Qin and Peng. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Jian Qin, qinjian621@163.com; Min Peng, mpeng320@whu.edu.cn
†These authors have contributed equally to this work