Exploration of bacterial lipopolysaccharide-related genes signature based on T cells for predicting prognosis in colorectal cancer

Acquisition and processing of data

The scRNA-seq dataset GSE200997 was obtained from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov) database, which included 16 CRC samples and 7 adjacent samples of normal colon tissue. There was a total of 49,859 cells in all samples, including 31,586 cells from tumor samples and 18,273 cells from normal samples. A total of 5286 LPS-related genes were downloaded from CTD database (http://ctdbase.org/).

The Cancer Genome Atlas Colon Adenocarcinoma (TCGA-COAD) cohort was downloaded from UCSC Xena platform (https://xenabrowser.net/datapages/), which is used to construct the prognosis model, including the mRNA expression data, clinical information (gender, age, microsatellite status, overall survival (OS), pathological stage, TNM stage), and mutation profiling data. The detailed information was shown in Supplementary Table 1. The validation dataset (GSE39582) was downloaded from the GEO database, and the detailed information was shown in Supplementary Table 2.

According to the published report [15], we gained an immunotherapeutic cohort (IMvigor210), which was treated with atezolizumab (anti-PD-L1 agent). Based on the guideline on http://research-pub.gene.com/IMvigor210CoreBiologies, the IMvigor210CoreBiologies R package was used to acquire Expression sets and clinical information of this cohort [16]. The samples were divided into platinum-treated cohort and non-platinum-treated cohort based on whether or not they received platinum-based chemotherapy.

scRNA-seq dataset analysis and TALRGs identification

The cell filtration criteria in the scRNA-seq dataset (GSE200997) were consistent with previous study [17]. According to the previous report [17], the filtered data was used to reduce the dimensionality of the features. The cells were then clustered using the FindNeighbors and FindClusters functions. The t-distributed stochastic neighbor embedding, and uniform manifold approximation and projection (UMAP) algorithms were used to reduce the dimensionality and visualize clustering classification based on the first 20 principal components selected. Then, we used SingleR package of R software to annotate the cell types of clusters, with “HumanPrimaryCellAtlasData” as reference [18]. Subsequently, we identified DEGs between tumor and normal cells in T cells using the FindMarkers function with logfc.threshold = 0.585. Moreover, the Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of DEGs in T cells were conducted using the Database for Annotation, Visualization and Integrated Discovery (DAVID) database (https://david.ncifcrf.gov/conversion.jsp). Finally, we selected the intersection of DEGs in T cells and bacterial LPS-related genes to obtain TALRGs.

Development and validation of the TALRGs prognostic model

TCGA-COAD cohort was used to recognize the prognostic TALRGs and established a prognostic risk model. Univariate Cox analysis was used to screen the prognostic signature using Survival R package. Kaplan–Meier analysis was performed to evaluate the relationship between mRNA expression of prognostic TALRGs and overall survival (OS). According to the each TALRG’ coefficients (which obtained from the multivariate Cox results), the risk score was calculated by the following formula [17]:

Kaplan–Meier analysis and time-dependent receiver operating characteristic (ROC) curve were applied to validate the predictive ability of the TALRGs prognostic model using the survminer R package and survivalROC R package, respectively. To further confirm the prediction model, we introduced a validation dataset: GSE39582. The prediction ability of the model was verified by the above method.

Assessment of clinicopathological characteristics

TCGA-COAD cohort and GSE39582 dataset were performed to validate the prognostic capability of TALRGs prognostic model. According to the constructed risk scoring model, the dataset was divided into high-risk group and low-risk group. Kruskal-Wallis test and Wilcoxon test were used to investigate the relationship between TALRGs prognostic model and clinicopathological characteristics, including pathological stage and TNM stage.

Establishment of nomogram

Nomogram is a multivariate regression analysis based on clinicopathological information and risk score. Based on the risk score and clinical factors such as age, sex, microsatellite status, and pathological stage, the nomogram was constructed using the TCGA-COAD cohort using rms R package. The calibration curve was employed to evaluate the suitability of clinical use.

Estimation of immunological properties of TALRGs signature

Similar to the previous approach [16], we employed CIBERSORT algorithm to evaluate the infiltration of 22 types of tumor-infiltrating immune cells [19]. The Wilcoxon test was applied to compare to the proportions of immune cells between low- and high-risk groups (P-value < 0.05). Then, we employed unpaired t-test and Wilcoxon test to assess the difference of the mRNA levels of immune checkpoints and IC₅₀ value of chemotherapy drugs between the two risk groups, respectively.

Analysis of mutation characteristics of TALRGs signature

We employed the maftools R package to analyze and visualize the somatic mutation profiles. Subsequently, the correlation between the mutation status of gene and OS was evaluated. The value of tumor mutation burden (TMB) was calculated and visualized by the maftools R package. The Kruskal-Wallis test was performed to compare the TMB value between the two risk groups. The surv_cutpoint algorithm of survival R package was used to determine the optimal cutoff value of TMB, and then all samples were divided into high-TMB group and low-TMB group. Kaplan–Meier analysis was used to evaluate the relationship between TMB values and OS. P < 0.05 was considered statistically significant.

The role of the TALRGs signature in immunotherapy

Similar to the previous study [16], we validated the constructed TALRGs signature using platinum-treated cohort and non-platinum-treated cohort, respectively. Secondly, patients with complete response (CR) or partial response (PR) were classified as responders and compared with non-responders (stable disease (SD) or progressive disease (PD)), and the risk score for each patient was calculated according to the built risk score model. Kaplan–Meier analysis and ROC curve were applied to assess predictive ability of the risk model. Finally, we employed the Wilcoxon test to compare the tumor mutation load and neoantigen load between the two risk groups.

scRNA-seq analysis of CRC and identification of TALRGs

In all, 47,958 high-quality cells were obtained after strict quality control (Figure 1A), including 30,341 tumor cells and 17,617 normal cells. The dimensionality reduction by principal component analysis (PCA) showed that tumor cells and normal cells were significantly different (Figure 1B). Subsequently, 22 distinct cell subgroups were identified and visualized by using UMAP method (Figure 1C). Through the annotation of 22 cell subgroups, a total of 6 major cell types were identified, mainly including T cells, B cells, NK cells, endothelial cells, hepatocytes, and monocytes (Figure 1D). Subsequently, we compared and analyzed the genes between tumor and normal cells in individual T cell clusters. T cells revealed 82 DEGs between tumor and normal cells (Supplementary Table 3). KEGG enrichment analysis showed that DEGs in T cells are mainly involved in “Antigen processing and presentation”, “Intestinal immune network for IgA production”, “Staphylococcus aureus infection” and “Cell adhesion molecules” (Figure 1E), which are correlative signaling pathways of bacterial LPS-related genes. A total of 5,286 bacterial LPS-related genes were obtained from the CTD database. Then LPS-related genes were overlapped with DEGs of T cells, and 44 TALRGs were obtained (Figure 1F).

Figure 1. TALRGs identification from scRNA-seq data and CTD. (A) The quality control of scRNA-seq data. (B) The dimensionality of tumor cells and normal cells was reduced by PCA. (C) UMAP distribution of the 22 cell clusters. (D) The cell clusters were annotated into 6 main cell types. (E) KEGG analysis of DEGs in T cells. (F) The Venn diagram of DEGs in T cells and LPS-related genes.

Establishment and validation of TALRGs signature

We further studied the association of TALRGs with prognosis, and 5 optimal prognosis-related TALRGs with P < 0.05 were identified using univariate Cox analysis (Table 1), namely HSPA1A, HSPB1, SERPINA1, SLC2A3, and TIMP1. Compared with normal samples, the expressions of HSPA1A, HSPB1 and SERPINA1 were down-regulated in tumor samples, while TIMP1 was up-regulated in tumor samples. SLC2A3 expression had no significant difference between normal and tumor samples (Figure 2A). Furthermore, Kaplan–Meier survival analyses demonstrated that patients with increased HSPA1A (P = 0.018, Supplementary Figure 1A), SLC2A3 (P = 0.0073, Supplementary Figure 1D), and TIMP1 (P = 0.0031, Supplementary Figure 1E) showed poor survival, while patients with increased SERPINA1 indicated favorable survival (P = 0.02, Supplementary Figure 1C). There was no significant relationship between HSPB1 expression and OS (P = 0.28, Supplementary Figure 1B). Subsequently, multivariate Cox analysis established a TALRGs signature based on the TCGA-COAD cohort, and the hazard ratio of each gene was shown in Figure 2B. All patients were then divided into high- and low-risk groups based on the median risk score. The formula for the risk score model was as follows: risk score = (0.15872 × expression value of HSPA1A) + (0.04431 × expression value of HSPB1) + (−0.19961 × expression value of SERPINA1) + (0.0404 × expression value of SLC2A3) + (0.22075 × expression value of TIMP1). The Kaplan–Meier analysis manifested low-risk patients had a better survival advantage (P < 0.0001, Figure 2C). In addition, the area under the curve (AUC) of the TCGA-COAD cohort was 0.624 (at 1-years), 0.639 (at 3-years), and 0.648 (at 5-years), respectively (Figure 2D). The model was validated in the GSE39582 cohort. The OS rate of the low-risk group was higher than that of the high-risk group (P = 0.00019, Figure 2E). The AUC of anti-PD-L1 cohort was 0.59 with P = 0.0366463 (Figure 2F). These results suggested that TALRGs prognostic model is effective in predicting the prognosis of CRC patients.

Figure 2. Construction and validation of the TALRGs signature. (A) Comparison of the TALRGs expression between tumor and normal samples. (B) The multivariate Cox analysis. (C) Kaplan–Meier analysis of the TALRGs signature in TCGA-COAD cohort. (D) ROC curve analysis of the TALRGs signature in TCGA-COAD cohort. (E) Kaplan–Meier analysis of the TALRGs signature in GSE39582 cohort. (F) ROC curve analysis of the TALRGs signature in anti-PD-L1 cohort.

Evaluation of clinicopathological characteristics

To explore the relationship between prognostic values and clinicopathological characteristics of CRC patients, clinical traits in TCGA-COAD and GSE39582 cohorts were analyzed. The CRC patients with more advanced clinical stages in TCGA-COAD cohort, such as advanced pathological stage (P = 0.00022, Figure 3A), later T stage (P = 0.0032, Figure 3B), later M stage (P = 0.0042, Figure 3C), and later N stage (P = 0.0017, Figure 3D) had higher risk scores, indicating patients with low-risk scores have a better prognosis. Analogously, there were differences in the risk scores of different pathological stages (P = 0.018, Figure 3E) and different T stages (P = 0.0051, Figure 3F) in GSE39582 cohort, suggesting that risk scores were correlated with pathological stages and T stages. Besides, the risk score in M1 was not significantly different from M0 (P = 0.26, Figure 3G), while the risk score in N2/N3 was significantly higher than N0/N1 (P = 0.00045, Figure 3H). These results indicated that TALRGs signature has good clinical performance.

Figure 3. Analysis of clinicopathological characteristics. The Kruskal-Wallis test was used to analyze the relationship between risk score and (A) advanced pathological stages, (B) T stages, (C) M stages, and (D) N stages in TCGA-COAD cohort. The Kruskal-Wallis test was used to analyze the relationship between risk score and (E) advanced pathological stages and (F) T stages in GSE39582 cohort. The Wilcoxon test was performed to analyze the relationship between risk score and (G) M stages and (H) N stages in GSE39582 cohort.

Construction of a nomogram

As shown in Figure 4A, age (≤ 60), stage III and IV, and risk score were observably associated with OS (P < 0.05), while gender, microsatellite status, stage I and II were not (P > 0.05). Then, we constructed a nomogram plot to ameliorate the accuracy of TALRGs signature and provided a quantitative and intuitive way to predict 1-, 3-, and 5-years survival probabilities for CRC patients (Figure 4B). We can determine the total score based on the individual score and thus predict the probability of survival of the patient at 1-, 3-, and 5-years. In addition, the calibration curves of the model showed that the predicted curves are close to the ideal curves, indicating that the nomogram has a good prediction performance (Figure 4C–4E).

Figure 4. Construction and validation of a nomogram. (A) The forest map demonstrated the association between risk score, clinical characteristics, and OS. (B) Nomogram based on risk score and clinical characteristics in predicting 1-, 3-, and 5-years OS. Calibration curves were used to predict the probability of 1-year (C), 3-year (D), and 5-year (E) OS.

Immune infiltrate analysis and efficacy prediction of TALRGs signature

Given that the characteristics of tumor immune infiltrating cells are associated with tumor development and progression and may affect the prognosis of CRC, we used the CIBERSORT method to compare the immune infiltration of 22 immune cells in different risk groups. The results showed B cells, plasma cells, resting CD4+ T cells, activated CD4+ T cells, and resting dendritic cell (DC) cells were more abundant in the low-risk group, while CD8+ T cells, T regulatory cells (Tregs), and M0/M1/M2 macrophages were more abundant in the high-risk group (Figure 5A). In addition, we analyzed correlations between the risk model and CD274 (PDL1), PDCD1 (PD1), CTLA4, HAVCR2 (TIM3), LAG3, and TIGIT. We found observably higher expressions of all immune checkpoint molecules in the high-risk group (Figure 5B), indicating that immune function is suppressed in high-risk patients. Moreover, we discovered the IC₅₀ of dactinomycin, docetaxel, mitoxantrone, vinorelbine, camptothecin, and teniposide was higher in the high-risk group, while the IC₅₀ of rapamycin was higher in the low-risk group (Figure 5C–5J), suggesting that patients in the high-risk group are more resistant to most chemotherapy drugs.

Figure 5. Immune infiltrate analysis and efficacy prediction of TALRGs signature. (A) The Wilcoxon test was applied to compare the proportions of immune cells between low- and high-risk groups. (B) Comparing the difference of the mRNA levels of immune checkpoints unpaired between the two risk groups using unpaired t-test. (C–J) Comparing the difference of IC50 value of chemotherapy drugs between the two risk groups using Wilcoxon test.

Mutation analysis of TALRGs signature

Genetic alteration analysis showed that APC, TP53, TTN, KRAS, SYNE1, MUC16, PIK3CA, FAT4, RYR2 and OBSCN were the mutation rates of the top 10 genes in CRC patients (Figure 6A). The mutant profiles of high-risk and low-risk patients were shown in Supplementary Figure 2A, 2B, respectively. Besides, patients in the high-risk group had a higher median number of mutations than those in the low-risk group (118 vs. 110). In the high-risk group, single nucleotide polymorphism (SNP) was the most variant type, missense mutation was the most common SNP classification, and C > T was the main single nucleotide variants (SNV) class, and TTN, APC, MUC16, SYNE1, TP53, RYR2, DNAH11, ABCA13, KRAS, PIK32CA were the top 10 mutated genes in patients in the high-risk group (Figure 6B). Except for the top 10 mutated genes (TTN, APC, MUC16, SYNE1, TP53, FAT4, KRAS, OBSCN, PIK3CA, ZFHX4), the mutation profiles of the low-risk group were similar to that of the high-risk group (Figure 6C). Moreover, the mutation rate of TP53 in the high-risk group was significantly higher than in low-risk group (P = 0.0257, Supplementary Figure 2C). The mutation status of TP53 didn’t affect patients’ OS (P = 0.183, Supplementary Figure 2D), while TP53 combined with mutations in other genes (e.g. LRP18, SYNE1, TTN) significantly reduced survival time (Figure 6D). In addition, the median TMB value was 2.22/Mb in CRC patients (Supplementary Figure 2E), and there was no significant difference in TMB value between low-risk group and high-risk group (P = 0.8, Supplementary Figure 2F). However, patients with lower TMB had better prognosis (P = 0.028, Figure 6E).

Figure 6. Analysis of mutation characteristics of TALRGs signature. (A) The mutation profiles of CRC patients. (B) The summary information of high-risk patients’ mutation profiles. (C) The summary information of low-risk patients’ mutation profiles. (D) Kaplan–Meier analysis of the two genes mutation status. (E) CRC patients with higher TMB values had poor OS.

Effect of immunotherapy on prognosis

In order to evaluate the predictive effect of TALRGs signature on immunotherapy, we further validated constructed risk score model using the platinum-treated cohort and the non-platinum-treated cohort. The results exhibited patients with low-risk scores had significantly longer survival times (platinum-treated cohort: P = 0.0032, Figure 7A; non-platinum-treated cohort: P = 0.00017, Figure 7B). The AUCs of constructed model in anti-PD-L1 cohorts were 0.639 (platinum-treated cohort, Figure 7C) and 0.627 (non-platinum-treated cohort, Figure 7D), respectively. The clinical benefit of immunotherapy was not significant in patients receiving platinum therapy (P = 0.27, Figure 7E), while patients receiving non-platinum therapy marked clinical benefits from immunotherapy in the low-risk score group compared to those with a high-risk score group (P = 0.023, Figure 7F). In addition, higher tumor mutation load and higher neoantigen burden in platinum-treated patients was not significantly associated with risk score (P = 0.29, Figure 7G; P = 0.94, Figure 7H), while non-platinum-treated patients were associated with a low-risk score (P = 5.1e-06, Figure 7I; P = 1.7e-05, Figure 7J). These findings demonstrated that TALRGs signature is effective in predicting the prognosis of immunotherapy patients.

Figure 7. Effect of immunotherapy on prognosis. Kaplan–Meier analysis of the TALRGs signature in the platinum-treated cohort (A) and non-platinum-treated cohort (B). ROC curve analysis of anti-PD-L1 cohort in the platinum-treated cohort (C) and non-platinum-treated cohort (D). Comparing the difference of the risk score between responders group and non-responders group in the platinum-treated cohort (E) and non-platinum-treated cohort (F). Comparing the difference of the tumor mutation load and neoantigen load between low- and high-risk groups in the platinum-treated cohort (G, H) and non-platinum-treated cohort (I, J).