Prognostic Implications and Immune Infiltration in the Nod-like receptor signaling pathway: A Comprehensive Analysis across Pan-cancer

Methods

Gene expression data from TCGA related to the NOD signaling pathway were integrated with clinical data. Prognostically relevant NOD pathway genes were analyzed using univariate and multivariate Cox regression and Kaplan-Meier survival analysis. The accuracy of our prediction model was validated through receiver operating characteristic (ROC) curve analysis.Single-cell analysis of genes associated with reduced survival in patients, and single-sample immunoinfiltration analysis revealed cell-level differences between different groups.

Results

Univariate Cox regression analysis, multivariate cox regression analysis and Kaplan-Meier analysis were used to identify prognostic genes in NOD pathway. TRAF5 is an important prognostic gene in multiple cancer types, and mutation analysis showed that patients with TRAF5 mutations had reduced survival. Immune infiltration analysis revealed differences in effector memory CD8 T cells and immature B cells between high- and low-risk groups, suggesting potential druggable targets. Single-cell analysis highlighted that reduced survival was associated with overexpression of TXN in both primary and metastatic tissues.

Conclusion

NOD signaling pathway, specifically TRAF5, plays a critical role in cancer prognosis across various cancer types. Immune infiltration disparities offer therapeutic opportunities, and TXN represents a promising target for novel anticancer treatments.

Data collection and processing

Nod-like receptor signaling pathway genes were obtained from the KEGG database (https://www.genome.jp/kegg/), totaling 184 genes (Table S1). After excluding cancers with limited sample cases or those lacking normal sample controls, gene expression matrices for both normal and diseased samples for 10 cancer types were acquired from the Cancer Genome Atlas (TCGA). Additionally, clinical data, including patient gender, age, survival status, pathologic stage, and survival period, were collected. The 10 cancer types encompassed Bladder Urothelial Carcinoma (BLCA), Breast Invasive Carcinoma (BRCA), Colon Adenocarcinoma (COAD), Esophageal Carcinoma (ESCA), Glioblastoma Multiforme (GBM), Head and Neck Squamous Cell Carcinoma (HNSC), Lung Adenocarcinoma (LUAD), Lung Squamous Cell Carcinoma (LUSC), Rectum Adenocarcinoma (READ), and Stomach Adenocarcinoma (STAD). Table S2 lists the full names, tumor samples, normal samples and total samples of the ten cancers. Subsequently, the data underwent thorough cleaning and preprocessing to eliminate samples with missing critical information, ensuring the availability of a sample collection for each cancer type for subsequent analysis.

Identification of Differentially Expressed Genes

For the analysis of Nod-like receptor pathway genes, we employed the "edgeR" R package to identify differentially expressed genes between normal and tumor samples. The criteria for screening were set as |log2 fold change (FC)| > 1 and a false discovery rate (FDR) of < 0.05

Univariate and Multivariate Cox Regression Analysis

We performed univariate Cox regression analyses using the "survival" software package to assess the relationship between each gene and patient survival

^16,

¹⁷

. The upper and lower limits of the corresponding hazard ratio (HR) and 95% confidence intervals (HR95H, HR95L) are calculated. Risk score HR> 1 indicates that the expression level of this gene is negatively correlated with the prognosis, that is, the higher the expression level, the worse the prognosis; HR<1 indicates that the expression level of this gene is positively correlated with the prognosis, that is, the higher the expression level, the better the prognosis; HR=1 indicates that the gene is not associated with prognosis. Subsequently, the "glmnet", "survival" and "survminer" packages were used to conduct multivariate cox regression analysis

¹⁸

, establish the prognosis model, and calculate the sample risk score.

Mapping Survival and Risk Curves

The risk scores for patients across the 10 cancer types were computed using the risk score formula derived from the prognostic model. Subsequently, patients were stratified into high-risk and low-risk groups. To investigate the connection between risk scores and overall survival, a log-rank test was conducted for each cancer type employing the "survival" software package. Statistical significance was defined as p-values less than 0.05.

Moreover, risk curves and heat maps were generated using the "pheatmap" software package to further explore the relationship between risk scores and survival outcomes

^19,

²⁰

.

Independent Prognostic Analysis

In the independent prognostic analysis, the risk score derived from the prognostic prediction model was employed as a predictor variable, while other clinical factors were included as covariates. This analysis was conducted to evaluate the independent prognostic predictive capability of the risk score. The predictive performance of the risk score for patient survival, in consideration of other factors, was assessed by calculating the hazard ratio (HR) and its corresponding 95% confidence interval (HR95H, HR95L).

Construction of ROC curves and nomograms

Receiver Operating Characteristic (ROC) curves were constructed using the "survival," "survminer," and "timeROC" software packages to evaluate the predictive performance over time. The Area Under the Curve (AUC) was utilized to illustrate the accuracy of predictions. A higher AUC value indicates that the corresponding ROC curve is closer to the upper-left corner, signifying a higher true positive rate and a lower false positive rate. In essence, a higher AUC value corresponds to greater prediction accuracy

^21,

²²

.

The nomogram is a clinical prediction tool that leverages the outcomes of multivariate Cox regression analysis

²³

. It amalgamates the associations of multiple variables to visually represent the impact of each variable on the prognosis. The nomogram assigns scores to the range of values for each of the selected clinical factors in the model, based on their respective influence on the final outcome variable. When predicting an individual sample, the scores for each influencing factor can be summed to yield a total score. The relationship between the total score and the probability of the outcome event is then utilized to diagnose or predict the onset and progression of the disease, thus furnishing the predictive value of the individual outcome event.

Single-Sample Immunoinfiltration Analysis

To evaluate the correlation between gene expression and the abundance of eight distinct types of immune-infiltrating cells (activated CD4 T cells, activated dendritic cells, effector memory CD8 T cells, γ δ T cells, immature B cells, natural killer cells, neutrophils, and plasma cell-like dendritic cells), we employed the R package "GSVA" in conjunction with the ssGSEA algorithm. Specifically, we categorized samples into high and low expression groups based on risk scores derived from multivariate Cox regression analysis. Subsequently, the ssGSEA algorithm was employed to assess the correlation between the expression of nod-like receptor-related genes and the abundance of immune-infiltrating cells. The data were analyzed on the online tool of HiOmics Cloud Platform (https://henbio.com/en/tools)

²⁴

.

Mutation Analysis

Further analysis of mutations in nod-like receptor pathway-related genes was conducted using the online platform cBioPortal (https://www.cbioportal.org/). This analysis aimed to determine the mutation frequency and mutation types of the selected genes and to explore the association between mutations and prognosis. The findings of the mutation analysis offer valuable insights into the potential functions and repercussions of gene mutations in various cancer species.

Single-Cell Analysis

Datasets for head and neck squamous carcinoma (GSE227156/GSE173468), gastric adenocarcinoma (GSE163558), breast carcinoma (GSE161529), and thyroid carcinoma (GSE184362) were sourced from the GEO database of the National Center for Biotechnology Information (NCBI). These datasets encompassed data from carcinoma in situ (T), metastatic tissue (LN), and paraneoplastic tissue (N). Multi-sample merge analysis on these diverse cancer datasets was executed utilizing the Seurat R package (version 4.3.0).

The data underwent an initial quality control process followed by NormalizeData normalization. Subsequently, highly variable genes were identified, and data normalization was conducted using the ScaleData function. Principal component PCA analysis was employed to reduce the dimensionality of the data. The Harmony function was utilized to mitigate batch effects. Cell projection and clustering analyses were performed using nonlinear dimensionality reduction and the FindClusters function, and cell types were visualized through UMAP analysis. Cell clusters were characterized using the FindAllMarkers function to identify marker genes for each cell cluster. Literature references cited in the dataset were used for cell cluster identification. Bubble plots were generated to illustrate the expression of genes (MAPK9, MAPK10, TXN, TXN2, IFNAR2, CCL2, IL6, IL1B, IKBKE) in each sample.

Identification of Differentially Expressed Genes

To identify differentially expressed genes between cancerous and normal tissue samples, we used log2FC values obtained from differential analysis to generate heat map (Fig. 1A) and box patterns (Fig. 1B-G). Our analysis revealed differential expression of ANTXR2, CXCL2, ITPR1, GPRC6A, PYDC1, and TXNIP between cancerous lesions and normal tissue samples in multiple cancers.

Prognostic Gene Value in Pan-Cancer

Prognostic genes with relevance to pan-cancer were initially identified through Univariate Cox regression analysis (Table S3). Subsequently, these genes underwent Multivariate Cox regression analysis to pinpoint the most suitable genes for constructing the prognostic risk model (Table S4). Utilizing the risk scoring formula, the risk score for each sample was computed, and samples were ranked from smallest to largest. Using the median of the risk scores as the threshold value, the samples were categorized into high-expression and low-expression groups. The distribution of differences in survival status and survival time between the high- and low-risk groups was visualized (Fig. 2).

The results revealed a significantly higher number of deceased cases in the high-risk region for cancers such as BRCA, COAD, ESCA, HNSC, LUAD, READ, and STAD compared to the low-risk group. Furthermore, the heatmap illustrating the differential expression of prognostic genes between high and low-risk groups demonstrated that all genes associated with a favorable prognosis (HR < 1, Table S4) were significantly up-regulated in the low-risk group, while all genes unfavorable to prognosis (HR > 1, Table S4) were significantly up-regulated in the high-risk group.

Kaplan-Meier curves (Fig. 3A-J) clearly demonstrated a substantial increase in survival chances for patients in the low-risk group compared to those in the high-risk group across various cancer types, including BLCA (p < 0.001), BRCA (p < 0.001), COAD (p < 0.001), ESCA (p = 0.025), GBM (p < 0.001), HNSC (p < 0.001), LUAD (p < 0.001), LUSC (p < 0.001), READ (p < 0.001), and STAD (p < 0.001). This observation underscored the validity of the risk score as a robust prognostic indicator, with statistical significance defined by p < 0.05. Furthermore, time-dependent ROC curves (Fig. 3K-T) were employed to assess the accuracy of the 10 cancer prognosis-related genes in predicting the 1-, 3-, and 5-year overall survival of cancer patients. The area under the curve (AUC) values indicated that the prognostic model exhibited strong predictive capability, thereby identifying the aforementioned prognostic-related genes as valuable prognostic markers for their respective cancer types. In addition, within the pan-cancer prognostic genes, it was noted that TRAF5 was featured in the prognostic genes of BLCA, BRCA, COAD, and READ, while GABARAPL2 was present in the prognostic genes of ESCA and STAD within the digestive system cancers.

Assessment of Independent Prognostic Factors in Patients with 10 Types of Cancer

To gauge whether the risk scores derived from the prognostic prediction model were influenced by other clinical factors, including patient age, gender, and clinical stage, both one-way and multifactor Cox regression analyses were conducted on the risk scores and clinical factors. The results of the univariate Cox regression analysis (Figures S1) demonstrated that the risk scores (P < 0.001) could function as independent prognostic factors, independent of other clinical factors, for all cancers except GBM (P = 0.998). Importantly, multifactorial Cox regression analysis (Fig. 4A-J) further confirmed that the risk score (P < 0.001) was an independent risk factor for these cancers, except for GBM.

Subsequently, time-dependent ROC curve analysis was employed to evaluate the accuracy of the predictive model. The area under the risk score curve (AUC) is as follows: 0.745, 0.703, 0.737, 0.816, 0.704, 0.765,0.674, 0.703, 0.810 and 0.768 (Fig. 4K-T). These AUC values further underscored the predictive accuracy of the model.

Alignment Diagram (Nomogram)

The alignment diagram, commonly referred to as a nomogram, is constructed based on multifactor regression analysis. It integrates multiple predictors and utilizes line segments with scales plotted on the same plane according to a predefined scale. The length of each line segment represents the range of values that a variable can assume. This graphical representation serves to express the predictive model, incorporating variables such as age, stage, gender, and risk score, which contribute to the magnitude of the outcome event (Fig. 5).

In clinical practice, healthcare professionals can employ the column line diagram. They input the actual clinical information of the patient into the diagram, determining the corresponding scores for each variable based on the score scale. These scores are then summed to obtain a total score. This total score can be placed on the total score scale, and a vertical line can be drawn on the survival rate scale at the corresponding position of the total score scale to estimate the patient's survival rate at 1 year, 2 years, and 3 years.

Analysis of Differences in Immune Cell Infiltration Between Different Samples

Single-sample gene set enrichment analysis (ssGSEA) was utilized to compute the infiltration scores of individual immune cells. The immune infiltration heatmap revealed variations in the levels of specific immune cells. For instance, activated CD4 T cells exhibited higher immune infiltration levels in ESCA and LUSC, with infiltration scores of 0.48 and 0.42, respectively. Conversely, neutrophils demonstrated higher immune infiltration scores in COAD and READ, with scores of 0.65 and 0.37, respectively (Fig. 6A).

The analysis of differences in immune cell infiltration between high-risk and low-risk groups uncovered significant disparities. Effector memory CD8 T cells exhibited elevated infiltration in COAD (p = 8e-05), ESCA (p = 0.0025), HNSC (p = 0.033), LUAD (p = 0.012), READ (p = 0.0047), and STAD (p = 0.00014) when comparing high-risk and low-risk groups (Fig. 6B-K). In these cases, effector memory CD8 T cell infiltration levels were higher in the high-risk groups than in the low-risk groups(Supplemental Fig. 2).

Similarly, immature B cells showed variations in infiltration between high-risk and low-risk groups. Specifically, in BRCA (p = 0.03), COAD (p = 0.00016), GBM (p = 0.0058), LUSC (p = 3.7e-09), READ (p = 0.047), and STAD (p = 1.3e-07), statistically significant differences were observed (Fig. 6B-K). Immature B cell infiltration levels were notably higher in the high-risk groups than in the low-risk groups for BRCA, COAD, GBM, LUSC, READ, and STAD (Figures S2).

Single-cell analysis and mutational analysis of TRAF5 in cancer

The frequency of TRAF5 mutations in the TCGA database was explored using the CBioPortal website, covering data from ten studies encompassing 5667 samples. The findings revealed that TRAF5 mutations were present in 3% of patients (Fig. 7A). Notably, copy number variations (CNV) were predominant in BRCA, COAD, ESCA, LUAD, LUSC, and STAD, with CNV occurring in more than 8% of cases in BRCA, although no CNV was observed in READ (Fig. 7B).

A comprehensive analysis of mutation sites uncovered a total of 40 mutation sites, comprising 30 missense , 7 truncating, and 3 splice mutations (Fig. 7C). Furthermore, the correlation between TRAF5 gene mutations and clinical prognosis was examined. Samples with at least one gene change in mRNA expression were categorized as the altered group, while those without alterations were designated the unaltered group. The results demonstrated significant differences between the two groups in terms of disease-specific survival (DSS) (p = 1.475e-3), overall survival (OS) (p = 3.326e-3), and progression-free survival (PFS) (p = 1.383e-3) (Fig. 7D-F). Patients in the altered group exhibited significantly shorter survival than those in the unaltered group, indicating that TRAF5 gene mutations were associated with a poorer prognosis regarding DSS, OS, and PFS.

To elucidate the expression patterns of prognostic genes (MAPK9, MAPK10, TXN, TXN2, IFNAR2, CCL2, IL6, IL1B, and IKBKE) in different cell types of the tumor microenvironment across various cancers, we analyzed datasets of primary tumors, lymph node metastases, and adjacent non-cancerous tissues from BRCA, HNSC, STAD, and THCA, and visualized the data using bubble charts (Figure 8). Our findings revealed that these nine prognostic genes were primarily expressed in cancer cells. Notably, the TXN gene exhibited high expression levels in multiple cancer cell types across BRCA, HNSC, STAD, and THCA. Specifically, in HNSC and STAD, CD8+ T cells, epithelial cells, fibroblasts, and macrophages from primary tumor and lymph node metastasis tissues showed elevated expression of TXN, whereas its expression was lower in adjacent non-cancerous tissues. These results provide valuable insights into the complexity of cancer and the underlying biological mechanisms, paving the way for the development of novel therapeutic strategies and potential targets for future cancer treatments.