Prognostic Model and Tumor Immune Microenvironment Analysis of Complement-Related Genes in Gastric Cancer

Introduction

One key component of the innate immune response is the complement system. In addition to maintaining tissue homeostatic balance, the complement system regulates humoral immunity and pathogenic immune surveillance.^1,2 The complement system is involved in host defense mechanisms and pathological processes of inflammatory and immunological disorders and adaptive immune responses.³ In recent decades, the complement system is thought to act as an effector arm of the immune system, helping to destroy cancerous cells.⁴ Extensive research corroborates the clinical potential of complement inhibitors as a novel therapeutic approach for a variety of malignancies.^5,6 An increasing number of studies have demonstrated the role of each complement component and its modulators in anti-tumor immunity.⁷ Activation of complement components C5 and C5a in cancer tissues and C5aR expression in breast cancer cells is associated with poor prognosis. Studies show that inhibiting the C5a/C5aR pathway helps treat breast cancer patients.⁸ At the same time, C5a is a potential new marker of prognosis in colon cancer and targeting the C5aR pathway might help treat colon cancer.⁹ Complement component C3 could also be used as a marker for the early diagnosis of pancreatic cancer.¹⁰ C3 plays a crucial role in the development and growth of cSCC as a core component of the complement system.¹¹ The complement system is also critical in treating ovarian and lung cancers.^12,13

Gastric cancer (GC) is the world’s fifth most prevalent malignancy and the fourth principal contributor to cancer-associated mortality globally. GC was responsible for almost 769,000 fatalities in 2020.¹⁴ Immunotherapies are used to treat GC in addition to radiotherapy, gastrectomy, chemotherapy, and targeted treatment.¹⁵ Despite recent advances in medicinal approaches, the death rate of GC individuals is still high.¹⁶

In both healthy gastric tissues and GC tissues, the expression of complement family genes was analyzed in this study. We also examined how these genes’ expression levels correlated with the survival outcomes of patients suffering from GC. We established a risk score model using complement family genes, constructed a prediction nomogram, and assessed the connection between the model, immune cells, and immune checkpoints. The results validated the application prospect of a novel complement family gene-based prediction model for GC.

Materials and Methods

GC Datasets and Pre-Processing

Supplementary Figure 1 displays our research flowchart. The GEO (http://www.ncbi.nlm.nih.gov/geo/) and TCGA databases were searched for common gene expression data and clinical information on the mortality and survival rate of GC patients.^17,18 The GEO dataset (GSE84437) was coupled as the study object. We modified the background and carried out quantile normalization by combining the multi-array averaging techniques of simpleaffy and Affy. A thorough study of the cytogenetic mutation data and RNA sequencing (FPKM values) from the TCGA was carried out using the TCGAbiolink tool in.^19,20 Transcripts per kilobase million (TPM) were converted from FPKM values. The R Bioconductor package was applied to analyze the data in R (version 4.1.2).

Copy Number and Mutation Analysis of Complement Genes

An analysis of complement family gene expression levels was conducted based on the TCGA database. The complement family genes’ chromosomal distribution was retrieved with the RCircos tool,²¹ whereas CNV analysis was conducted with the aid of the Perl software (5.32.1.1) and R. (4.1.2).^22,23 The mutations were subjected to waterfall mapping using the Maftools software. The impact of individual mutations on the overall gene expression patterns was then assessed via the Student’s t-test. Additionally, the ggplot2 tool was employed to display the direction and magnitude of connections between gene expression and mutations.²⁴

Model Development for Risk Scoring

By using the expression profiles and survival data of complement family genes, the full TCGA set was employed as the training dataset and the full GSE84437 set as the test dataset. The complement family genes model was developed using the training dataset and validated using the full and test datasets. The association of pyroptosis-related genes with OS rate was investigated via the Lasso and Cox regression (“glmnet” and “survival” packages).²⁵ Cross-validation was implemented during the development of the Lasso model to guarantee reliability. We identified eight genes linked to survival using the penalty parameter (λ) and utilized these genes to develop a multivariate Cox regression model. The ideal gene set was selected and utilized to predict survival based on a forward-backward Cox regression method. The training and testing sets’ survival curves were generated utilizing a Kaplan-Meier (KM) method. The formula below was used to obtain the risk score:

The median risk score was used as cutoff for high-risk group (with higher risk score) and low-risk group (with lower risk score).

Nomogram Development and Validation for GC Patients

Based on the results of univariate and multivariate Cox regression analyses, a nomogram model was formulated for the prediction of GC patients’ OS at 1, 3, and 5 years. The nomogram was evaluated utilizing the “RMS” package²⁶ in “R”. The KM method was applied to create calibration curves comparing the real values to the nomogram-predicted probabilities. The scatter points were found to follow a 45° diagonal line by a model for predicting nomograms that are properly calibrated. The Harrell concordance index (C-index), which has a range of 0.5 to 1, was utilized to evaluate the accuracy of the predictions made by the nomogram model. All two-tailed statistical analyses in the current study had a significance threshold of p < 0.05 for each test.

Gene Set Enrichment Analysis (GSEA)

To conduct GSEA, the mean value of the risk score was applied to classify the candidate genes into two categories: high- and low-risk.²⁷ The Molecular Signatures Database (MSigDB [https://www.gsea-msigdb.org/gsea/msigdb]) served as the source for the functional pre-configured gene sets.²⁸ The parameter threshold of p < 0.05 and false discovery rate (FDR) < 0.25, was adopted to pinpoint the candidate genes involved in the pathway with significant enrichment. The significantly enriched signaling pathways were identified utilizing the adjusted p-value and normalized enrichment score (NES).

Assessment of the TME’s Immunological Properties in GC

To examine the association of risk score with immune cell infiltration, the “preprocessCore” software and CIBERSORT approach were utilized to evaluate the infiltration levels of 22 different immune cells.²⁹ In addition, the association of the risk score with 22 immune checkpoint molecules was examined in the TCGA. The significance criterion was established at a p-value of < 0.05.

Statistical Analysis

R software (v4.1.2),²³ GSEA (v4.2.3),³⁰ and Perl (v5.32.1.1)³¹ were used for analyzing the data in this study. Along with a description of the key R software packages applied, the statistical techniques used in the analysis were also provided. A p-value of < 0.05 indicated the cutoff criterion for significance.

Results

Expression of Complement Family Genes and Three Classical Pathways of the Complement System

The complement system consists of more than 50 protein components, and the exact number of complement genes is still unknown.³² In this study, a total of 53 complement-related genes were obtained as candidate genes by collating past literature and reviews^32–35 (Supplementary Table 1), of which the expression of 42 genes was extracted after combining the TCGA and GSE84437 databases. The complement system is activated by the three pathways: the lectin, alternative, and classical pathways (Figure 1A). To determine if abnormal expression patterns were linked to malignancy, we examined the mRNA levels of these 42 genes in GC versus adjacent non-cancerous samples. The expression of thirty complement genes varied significantly between GC and non-cancerous samples (Figure 1B). Moreover, these results lend credence to the notion that complement gene mRNA expression is strikingly diverse in GC and that aberrant gene expression is involved in the onset and progression of GC.

Figure 1 Expression of complement family genes and three classical pathways of the complement system. (A) Three classical pathways of the complement system. (B) Expression of 42 complement genes in normal gastric tissue (dark blue) and tumor tissue (dark Orange). Interquartile ranges are shown by box plots. Boxed rows represent medians, and P-values are depicted by asterisks (*p < 0.05; **p < 0.01; ***p < 0.001).

Complement Family Genes’ Genetic Variation Patterns in GC

Forty-two complement genes had high copy number variations (CNV) frequencies, which were primarily manifested as copy number amplification. In contrast, CNV deletion frequencies were prevalent in C1S, C1R, C1RL, CFD, C3AR1, CPN1, C8A, C8B, CD209, C3, and CLU (Figure 2A). The locations of CNVs in the complement genes on the chromosome are shown in Figure 2B. Additionally, we summarized the frequency of CNV and somatic mutations of 42 complement regulators in GC. Of 433 samples screened, 144 (32.33%) exhibited genetic alterations of complement regulators. Among the complement genes, C3, CFH, CR1, CR2, CFI, VTN, C1S, C5, C6, CFHR4, CFHR5, C3AR1, C8B, C9, CFHR2, CPN1, C1R, C4BPB, C8A, CFHR1, CLU, C4BPA, MASP2, C1QA, C1RL, C5aR1, CD209, C4A, and CD55 were mutated in GC, while other regulators were not mutated (Figure 2C). By means of the univariate Cox proportional risk regression analysis, 10 prognostic genes were detected. Forest plots were created and hazard ratios were determined (Figure 2D).

Figure 2 Molecular characterization and expression variants of complement family genes in GC. (A) Analyzing the frequencies of CNV in the complement family of genes in GC. The frequencies of deletions are shown by green dots, red dots represent amplification frequencies, and alteration in frequency is depicted by the column height. (B) Location of complement family gene CNV on 23 chromosomes. (C) The frequency of mutations in genes belonging to the complement family in 433 TCGA patients with GC. TMB is represented by the top bar, and each column denotes a single patient. In addition, each regulator’s mutation rate is shown to the right in the figure. The visual representation of the distribution of different types of regulatory systems is shown on the bar graph to the right. The percentage of conversions in each group is shown as a stack of horizontal bars. (D) The genes of the complement family were examined utilizing a univariate Cox regression analysis.

Risk Score Model Constructed Using Genes in the Complement Family

The validation set consisted of 433 samples from the GSE84437 cohort, whereas the training set consisted of 371 samples from the TCGA-GC cohort. The coefficients of the chosen complement genes in the validation set were estimated by Lasso-Cox regression analysis followed by the construction of a risk score model. We arrived at the risk score equation shown below for this model: Risk score = (0.067 * C5aR1 exp.) + (0.016 * C7 exp.) + (0.099 * C9 exp.) + (-0.030 * CD46 exp.) + (0.185 * CD59 exp.) + (0.050 * CFH exp.) + (0.523 * CFHR4 exp.) + (0.081 * VTN exp.) (Figure 3A and B). Samples from both groups were classified into high- and low-risk categories as per the median risk score. Notably, the training and validation datasets included measurements of gene expression, risk score, and survival. Specifically, in the training set, the mortality rate was 51.89% in the high-risk category and 26.34% in the low-risk category; in the validation group, the mortality rate was 52.78% in the high-risk category and 43.78% in the low-risk category. Importantly, the results of the survival study illustrated that patients with low risk exhibited substantially longer OS compared to those with high risk (Figure 3C and D; p < 0.05). The distribution of scores, changes in gene expression patterns, and eight independent prognostic genes were employed to design heat maps showing patient survival rates. These heat maps were utilized for the evaluation of the risk model developed using information from the TCGA platform (Figure 3E–J). Remarkably, the prognostic model exhibited remarkable prognostic power, suggesting that it may anticipate the onset and advancement of GC.

Figure 3 Construction of a prognostic model based on complement family genes. (A) The log (lambda) sequence is shown against the Lasso coefficient profiles of the 10 genes chosen for the OS analysis. (B) The 10-fold cross-validation error (the minimal error is represented by the first vertical line, and the error within 1 standard error of the minimal is represented by the second vertical line). (C and D) Survival probabilities for patients with high-risk score in contrast with those with low-risk score, as shown by Kaplan-Meier plots. (E–G) Distribution of risk scores, associated survival data, and expression of eight core genes in the TCGA cohort. (H–J) Distribution of risk scores, associated survival data, and expression of eight core genes in the GSE84437 cohort.

The Development and Evaluation of a Nomogram

A nomogram was developed to examine the complement-based risk model’s clinical efficacy in clinical settings (Figure 4A). Patients were assigned a total score by adding the points of all prognostic parameters. When the sum of the points was greater, the patients’ clinical outcomes were less favorable. The calibration plot showed that the nomogram’s findings were somewhat similar to those of a perfect model (Figure 4B). The nomogram’s ability to make accurate predictions was further demonstrated by the ROC curve (Figure 4C). Cox regression analysis confirmed the associations of the GC patients’ prognoses with the prognostic risk score model, stage, age, and grade.

Figure 4 Developing a nomogram and analyzing prognostic factors. (A) A nomogram model was constructed to predict the 1-, 3-, and 5-year OS probabilities of patients with GC. (B) Nomogram model-related calibration curves for estimating 1-, 3-, and 5-year OS for GC patients. (C) Nomogram model-related ROC curves for estimating 1-, 3-, and 5-year OS for GC patients. (D–E) The results of univariate (D) and multivariate Cox regression (E) for the OS of patients with EC are shown in forest plots. (**p <0.01; ***p < 0.001).

We identified the model’s function in independently predicting the prognosis by comparing the risk score model with clinical prognostic variables in the two datasets using univariate and multivariate analyses. Additionally, it was established that age, risk score, and stage were linked to the GC patients’ prognosis, with the risk score being an independent prognostic predictor (Figure 4D–E).

Gene Set Enrichment Analysis (GSEA)

To probe which signaling pathways are most important in relation to the eight key prognostic genes in GC, a GSEA was conducted (Figure 5A). The findings demonstrated the enrichment of these genes in pathways linked to cell adhesion molecules cams, complement and coagulation cascades, cytokine cytokine receptor interaction, ecm receptor interaction, hematopoietic cell lineage, and neuroactive ligand receptor interaction (Figure 5B–G). The eight key prognostic genes could therefore be crucial in the inflammatory, immune, and hematologic systems, enhancing the onset and advancement of GC.

Figure 5 Gene set enrichment analysis. (A) Top 10 pathways related to risk score-associated gene enrichment. (B) Cell adhesion molecules cams. (C) Complement and coagulation cascades. (D) Cytokine cytokine receptor interaction. (E) Ecm receptor interaction. (F) Hematopoietic cell lineage. (G) Neuroactive ligand receptor interaction.

The Prognostic Model’s Risk Score Was Correlated with 46 Prevalent Immune Checkpoint Genes (ICGs) in GC

In the TCGA-GC cohort, 46 ICGs and the risk score interactions were characterized (Figure 6). As opposed to the low-risk category, the high-risk patients’ ICG expression levels were remarkably elevated (Figure 6A). Moreover, the risk score and CTLA4, LAG3, CD274, and PDCD1 all showed favorable correlations (Figure 6B–F).

Figure 6 Correlational analysis of risk score with 46 GC common immune checkpoint genes. (A) Differential expression of the 46 immune checkpoint genes between high- and low-risk score groups. (B) Correlation of risk score with 46 immune checkpoint genes in GC. (C–F) Correlation of risk score with the expression of PDCD1, CD274, LAG3, and CTLA4 in TCGA-GC dataset. (*p < 0.05; **p <0.01; ***p < 0.001).

Risk Score Correlation with Immune Cell Infiltration

By studying the connection between the eight major prognostic genes and the 22 distinct types of TIICs, the CIBERSORT method was utilized to examine the risk model correlation with immune cell infiltration (Figure 7). Patients with high and low risk were compared regarding the infiltration levels of immune cells. The proportion of regulatory T cells CD4 memory activated, and T cells follicular helper was significantly higher in the low-risk group, whereas that of T cells regulatory, dendritic cells resting, and monocytes was significantly higher in the high-risk group (Figure 7A). The risk score was negatively correlated with monocytes, dendritic cells resting, macrophages M2, and eosinophils and positively correlated with T cells CD4 memory activated and T cells follicular helper (Figure 7B–H). As opposed to the low-risk patients, the high-risk population had a considerably higher Tumor Immune Dysfunction and Exclusion (TIDE)³⁶ score, illustrating the immune escape within the GC tumor immune microenvironment (TIME) is mediated by the high-risk score and that less successful immunotherapy treatment was achieved for patients with high-risk scores (Figure 7I).

Figure 7 Link between risk score and immune cell infiltration. (A) Box plots showing the differences in the proportions of 22 different types of immune cells in GC tumor samples between high- and low-risk groups. (B) Correlation between risk score and 22 immune cell types in GC (p<0.05 in red). (C–H) Correlation of risk score with (C) monocytes, (D) dendritic cells resting, (E) macrophages M2, (F) eosinophils, (G) T cells CD4 memory activated, and (H) T cells follicular helper in GC. (I) Assessment of immune evasion efficacy in the high- and low-risk-expression groups. (*p < 0.05; **p <0.01; ***p < 0.001).

Eight Key Genes’ Expression in GC Samples and Adjoining Normal Samples, Along with Their Link to the Prognosis of Patients with GC

Analysis of mRNA expression patterns from 32 matched pairs of GC and adjoining normal samples from the TCGA-GC cohort revealed that C5aR1, C9, CFHR4, and CD46 were all expressed at considerably greater levels in the tumor cells (p < 0.05, Figure 8A, C, D, and G), while C7, CD59, CFH, and VTN were not (Figure 8B, E, F, and H). The OS correlation with eight major prognostic genes was examined using the KM plotter. Patients with lower levels of C5aR1, CD59, and VTN expression had a favorable prognosis as opposed to those with higher levels (Figure 8I, M, and P; p < 0.05), which was not the case with C7, C9, CFHR4, CFH, and CD46 (Figure 8J, K, L, N, and O). The findings demonstrated that C5aR1 was overexpressed in GC samples in comparison to adjoining non-cancerous tissues and was linked to the unfavorable 10-year OS in GC patients. In summary, C5aR1 can be used as a biomarker of poor prognosis in GC.

Figure 8 Expression of 8 core genes in GC and adjacent tissues and prognosis in GC. (A–H) C5AR1, C7, C9, CFHR4, CD59, CFH, CD46, and VTN were overexpressed in the tumor tissues as compared to the paraneoplastic tissues in TCGA-GC cohort. (I–P) Survival curve showing the impact of expressions of C5AR1, C7, C9, CFHR4, CD59, CFH, CD46, and VTN on the OS in TCGA-GC dataset. (*p < 0.05; ***p < 0.001; nsp > 0.05).

Discussion

The complement system, which regulates tissue homeostasis and enhances immune surveillance, is an important component of the human immune system.¹ The complement system primarily comprises effectors, receptors, and modulators, with each complement component playing a unique role. Other than activating complement, complement components, and inhibitors have other activities. According to recent studies, complement proteins interact with cells in the TME and metastatic target organs, helping to regulate tumor progression, anti-tumor immunity, and angiogenesis.³⁷ We also found that activating the complement system improves tumor cell escape from the immune response, produces angiogenesis, and boosts the proliferation, migration, and invasiveness of cutaneous squamous cell carcinoma (cSCC). Components of the complement system found in cancerous cells might be employed as biological markers for diagnosing and classifying cancer as well as predicting its prognosis.³⁴ Complementing innate and adaptive immunity can influence the quality and amount of T-cell activation.³⁸ Immunotherapies that block CTLA-4 and PD-1/PD-L1 immune checkpoints have been quite successful in the past ten years. Patients with Hodgkin’s lymphoma and melanoma have benefited the most.³⁹ Recent studies show that combined immunotherapy targeting CTLA-4 and PD-1 has been more effective in treating GC than targeting them individually.⁴⁰ Future research should focus on enhancing immunotherapy’s efficacy.

In previous studies, the complement system has been explored in the prognosis and immune infiltration of colon and GC,³³ but no relevant prediction model has been constructed to provide insights into the influence of complement genes in the prognosis of GC and in the mechanism of tumor microenvironment. In this work, expression data taken from the TCGA-GC cohort was used to thoroughly analyze the molecular patterns of these genes. Additionally, the expression, somatic mutations, CNVs, and prognostic significance of genes from the complement family in GC were examined. The entire TCGA dataset was used as the training group, whereas the GSE84437 dataset was used as the validation group. Eight complement genes (C5aR1, C7, C9, CFHR4, CD59, CFH, CD46, and VTN) were utilized to create a prediction model in a Lasso-Cox analysis. The model has significant advantages in analyzing gene expression data. In addition, the LASSO regression model reduces the complexity of the model and improves the accuracy of the model’s predictions, resulting in good predictive efficiency for other datasets.^41,42 The participants from both groups were subsequently classified into high- and low-risk categories predicated on the median risk score. High-risk patients recorded a worse prognosis as opposed to the low-risk population. Our results show that the constructed model can distinguish GC patients well and predict prognosis, which can help to personalize the treatment plan according to the patient’s risk. The univariate and multivariate Cox regression analyses demonstrated that the risk score independently served as a prognostic indicator. A nomogram was used to verify the complement-based risk model’s clinical efficacy. The CIBERSORT tool was utilized to assess the TIICs distribution found in low- and high-risk patients. Monocytes, resting dendritic cells, macrophages M2, and eosinophils were inversely linked to the risk score, whereas activated memory CD4 T cells and T cells follicular helper were favorably correlated. Monocytes, which are innate immune cells, are crucial to the onset and spread of cancer.⁴³ Dendritic cells are powerful antigen-presenting cells linked to immunity.⁴⁴ Macrophages, an essential constituent of the human immunological system, are often differentiated into M1-like and M2-like activated macrophages, both of which are intimately related to inflammation. Pro-inflammatory and anti-inflammatory activities are primarily mediated by M1 and M2 macrophages, respectively.^45,46 M1 phenotypic responses that are anti-tumor and pro-inflammatory may be neutralized by M2 macrophages.⁴⁷ Eosinophils infiltrate a variety of malignancies and can control the growth of such tumors either directly by interfacing with the malignant cells or shaping the TME.⁴⁸ In addition, CD4 T cells serve a crucial function in tumor prevention and immunotherapeutic response.⁴⁹ In lymphoid tissues, a subset of T cells referred to as T follicular helper cells (T(FH) cells) serves as a helper for B cells and is among the most abundant and essential groups of effector T cells.⁵⁰

Tumor immune infiltration and immune escape are correlated with cancer prognosis and treatment response.^51–53 Most tumor cells express antigens that can mediate recognition by host CD8+ T cells. Depending on TME, immune evasion could be classified into two: cancers without a T-cell-inflamed phenotype that defend themselves against immunological assault by immune system ignorance or exclusion and those with a T-cell-inflamed phenotype that defend themselves via dominating inhibition activities of the immune system-suppressive pathways.⁵⁴ The TIDE³⁶ score results showed that the immune escape potential was higher in the high-risk category and that Immunotherapy was less effective for participants in the high-risk category. In addition, among the eight core prognostic genes, C5aR1 upregulated in GC tissues compared with adjacent normal tissues, and patients with GC with high expression of C5aR1 had a lower 10-year OS.

This study may have several limitations. First, the model shown here might be overfitted for the specific sample set population. In addition, more clinical studies should be conducted to validate the predictive model established using the eight complement-associated genes. Additionally, the study only used data from publicly available sources, which might have caused selection bias.

Conclusion

Complement family genes in GC were systematically examined for prognostic significance, mutation, and expression patterns. A risk model was constructed, and patients with GC were divided into the high- and low-risk groups. Patients in the low-risk category had a remarkably more favorable prognosis than those in the high-risk category. Furthermore, a nomogram was constructed, which anticipated the recurrence of GC in an accurate manner, thus facilitating the development of individualized treatment for GC. The construction of a prognostic model allows for a comprehensive assessment of the characteristics of the corresponding TME cell infiltration. Our study provides important information for the development of more targeted treatments for GC.