System analysis based on the ER stress-related genes identifies WFS1 as a novel therapy target for colon cancer

Data collection

A total of 419 ER stress-associated genes are derived from the Gene Set Enrichment Analysis (GSEA) website. Gene expression profiles of GSE40967 (training set) and its corresponding clinical characteristics are downloaded from the gene expression omnibus, containing 19 normal samples and 566 tumor samples. Furthermore, gene expression profiles of 41 samples from normal colon tissue and 473 samples from tumor specimens are derived from The Cancer Genome Atlas (TCGA) database (Downloaded 2022.01.01), which is used as validation set.

Removing batch effect

We adjusted the test set to the training set using R “limma” package and R “sva” package.

Identifying differentially expressed genes (DEGs)

We detected DEGs between the tumor and normal tissues by R “limma” package. False discovery rate (FDR) <0.05 and |log₂foldchange| >1 were set as the cutoff criteria. Then, heatmap was plotted with R heatmap package. Protein-protein interaction (PPI) network was built via STRING database. Cyto-hubba plug-in from Cytoscape software was performed to detect hub genes.

Developing and validating a prognostic model relied on ERSRGs

To determine which genes were associated with death, univariate Cox regression analyses were conducted using the training set. Candidates for ERSRGs associated with overall survival were those with p-value <0.05. By performing lasso regression, we were able to reduce the model complexity and multicollinearity by shrinking the coefficients. For the purpose of developing a survival model, we conducted multivariate Cox regression analyses. Applying the formula below, we calculated the risk scores for each patient:

$$\text{risk score}={\displaystyle \sum _{j=1}^n(Coe{f}_j\times Xj)}$$

Xj measures gene expression value of genes, and Coefj represents the regression coefficient. High-risk and low-risk patients were grouped based on median scores calculated by the risk score formula. We compared the overall survival gaps between high-risk and low-risk groups using the Kaplan-Meier survival curves produced by the R package “Survminer”.

Developing a nomogram and its corresponding calibration curve

A nomogram is established based on clinical characteristics combined with risk score by R “rms” package. Calibration curves are created to determine the deviation between predicted and actual survival status.

Gene set enrichment analysis (GSEA)

An analysis of GSEA is conducted in order to detect enriched biological pathways in low- and high-risk groups, respectively. (GSEA 4.2.1.).

Immune cell infiltration analysis

“Estimation of STromal and Immune cells in MAlignant Tumor tissues (ESTIMATE) algorithm” is used to evaluate proportion of immune and stromal cells of each sample. Survival curves were plotted via the Kaplan-Meier method to predict the prognosis of COAD patients in the low- and high-stromal/immune score groups. Moreover, “Cell-type Identification by Estimating Relative Subsets of RNA Transcripts (CIBERSORT)” are performed to calculate the relative proportion of 22 immune cell infiltration of COAD samples from the training set. Then, the results calculated by CIBERSORT were screened according to p-value <0.05. The relationship between risk scores and different types of immune cells of each sample was evaluated via R “corrplot” package.

Immune checkpoints analysis

Spearman correlation test is performed to determine the relationship between immune checkpoint molecules and risk scores.

Prediction of drug therapy response

The drug-response prediction was evaluated with the R pRRophetic package and 2D conformations of drugs were visualized by PubChem website.

Cell culture

HCT116 and DLD-1 cells are derived from ATCC. In this experiment, cells are cultured with RPMI-1640 medium supplemented with L-glutamine and 10% fetal bovine serum.

Western blotting

WFS1 antibody (11558-1-AP) is purchased from Proteintech. The protein level of WFS1 in cells was measured through typical western blotting process. Beta-actin served as a loading control for western samples. By using ECL reagent, the protein concentration of WFS1 and beta-actin is determined. Two independent siRNAs are employed to decrease the WFS1 expression. The sequences of these two siRNAs used in this study are indicated as follows: si-WFS1-1: 5′-GCA GCG AGU CCA AGA ACU ACA-3′, si-WFS1-2: 5′-GCG UGA CUG ACA UCG ACA ACA-3′.

Cell viability assay

HCT116 and DLD-1 cells are transfected with siRNAs targeting WFS1 or scrambled in 6-well plate. After 24 hours transfection, we seeded 4000 cells per well in 96-well cell culture plates. The OD value at 450 nm is measured at indicated time points via a CCK8 kit.

Clone formation tests

HCT116 and DLD-1 cells are transfected with siRNAs targeting WFS1 or scrambled in 6-well plate. After 48 hours transfection, we seeded 500 cells per well in 6-well cell culture plate. After 2 weeks culture, the cells are stained with crystal violet.

EdU assay

HCT116 and DLD-1 cells were seeded in 24-well plates at about 50% in confluency. In cells with 70% cell density, cells were transfected with scrambled or two independent siRNA targeting WFS1. After 48 hours culture, cells were incubated with EdU for 2 hours. Next, a solution of 4% paraformaldehyde was used to fix the cells. According to the manufacturer's instructions, the staining process was performed. Images were obtained through Nikon microscope and the ratio of EdU positive cells was measured via the ImageJ software.

Consent for publication

All authors consent to the publication of this study.

Availability of data and material

All data and R script in this study are available from the corresponding author upon reasonable request. Publicly available datasets were analyzed in this study, these can be found in The Cancer Genome Atlas (https://portal.gdc.cancer.gov/) and Gene Expression Omnibus (GSE40967). All authors read and approved the final manuscript.

Identify differentially expressed ERSRGs signature

The overall workflow of this study is depicted in Figure 1. According to the screening criteria of | log2FC | >1 and adjust P < 0.05, we identified 31 DEGs between COAD samples and colon normal tissues in the training set. Within COAD samples, 21 ERSRGs were dramatically down-regulated and 10 up-regulated compared to normal colon tissue (Supplementary Figure 1A–1C). A PPI network was obtained through uploading these 31 differentially expressed ERSRGs to STRING website (Supplementary Figure 1D), and CXCL8, CCND1, PSAT1, SLC7A5, EIF4EBP1, DDIT4, TRIB3, RRP9, NOLC1, and CEBPB were determined as hub genes with top ten scores (Supplementary Figure 1E).

Figure 1. Flow chart of experimental design in this study.

Construction of a prognostic model relied on ERSRGs

First, 26 ERSRGs were screened out via univariate Cox regression analysis univariate using the cutoff of P < 0.05, which were associated with overall survival of patients with COAD in the training set, including 9 risky genes and 17 protective genes (Figure 2A). In order to further shrink the variables, LASSO regression analysis was conducted based on these 26 ERSRGs (Figure 2B, 2C). Finally, 9 ERSRGs were determined to construct the prognostic model according to the multivariate Cox regression, including 3 risky genes and 6 protective genes (Figure 2D). The formula to calculate the risk score is indicated as follows: Risk score = (−0.34508 × AQP11) + (−0.5975 × BCL2) + (−0.46549 × EDEM2) + (−0.48618 × EXOSC7) + (0.373054 × FLOT1) + (−0.15286 × PPP1R1B) + (−0.39254 × PPP1R8) + (0.134627 × TXNIP) + (0.223501 × WFS1). COAD patients in training set were classified into high-risk and low-risk groups according to the median risk scores (Figure 2E). COAD patients in the high-risk group had worse prognosis compared to those in the low-risk group (Figure 2G). The heatmap (Figure 2H) presented expression levels of 9 ERSRGs in different groups. Kaplan-Meier survival curves suggested that COAD patients in the low-risk group had significantly better prognosis than those in the high-risk group (Figure 2F). The area under the receiver operating characteristic (AUC) curve was 0.669, showing that prognostic model we built had a good predictive performance (Figure 2I). Next, we detected whether the expression levels of these 9 ERSRGs were associated with prognosis of patients with COAD. We found that the expressions of APQ11, EDEM2, EXOSC7, PPP1R8, TXNIP, and WFS1 were strongly correlated with the OS of COAD patients (p < 0.05) (Supplementary Figure 2A–2I).

Figure 2. Construction of risk model in the GSE40967 cohort. (A) Univariate Cox regression analysis of ERSRGs associated with OS of COAD patients. (B) Lasso regression analysis of the ERSRGs based on univariate Cox regression analysis. (C) Cross-validation for tuning the parameter selection in the LASSO regression. (D) Multivariate cox regression analysis of the ERSRGs based on LASSO regression analysis. (E) Distribution of patients based on the risk score. (F) Kaplan-Meier curve of survival probability of patients in the high-risk group and low-risk group. Statistical tests were performed using the Chi-square test with statistical significance at P < 0.05. (G) Survival time and survival status of patients with different risk scores. (H) The heatmap of the expression of prognostic ERSRGs between high-risk group and low-risk group. (I) ROC curve of risk score and other indicators.

Validation of the prognostic model

To validate the performance of the prognostic model in independent cohort, we downloaded the gene expression files of patients with COAD from TCGA as the test set. The risk scores of patients from test set were calculated using the same formula. COAD patients from the test set were separated into high-risk (n = 202) and low-risk (n = 215) groups according to the median risk score of the training set as the cutoff point. (Figure 3A). Similar with patients in training set, those with high-risk scores had shorter OS than those with low-risk scores (Figure 3B). The heatmap (Figure 3C) showed expression levels of 9 ERSRGs of each patient in test set. Kaplan-Meier survival curves suggested that COAD patients in the low-risk group had dramatically better prognosis compared with those in the high-risk group (Figure 3D). Our prognostic model had a good predictive performance, with an area under the receiver operating characteristic curve of 0.633 (Figure 3E).

Figure 3. Verification of risk signature in the TCGA cohort. (A) The distribution of patients based on the risk score. (B) Survival time and survival status of patients with different risk scores. (C) The heatmap of the expression of prognostic ERSRGs in the high-group and low-risk group (D) Kaplan-Meier curve of survival probability of patients in the high-risk group and low-risk group. Statistical tests were performed using the Chi-square test with statistical significance at P < 0.05. (E) ROC curve of risk score and other indicators.

Determination of independent prognostic factors

To confirm whether the risk score of the prognostic risk model has independent prediction ability, we evaluate the prognostic value of risk score and clinical features in the training set. Univariate Cox regression analysis revealed that risk score, age, stage, T, M, and N could be used to predict outcomes of COAD patients (Supplementary Figure 3A). Furthermore, multivariate Cox regression verified that risk score, age, T, M, and N could still be regarded as independent prognostic predictors (Supplementary Figure 3B). These data demonstrated that prognostic risk model, age, T, M, and N were independently correlated with the outcomes of COAD patients.

Construction of the nomogram

To better predict outcomes of COAD patients, we built a nomogram relied on risk scores combined with clinical characteristics to assess the 1-, 3- and 5-year survival probability (Supplementary Figure 4A). The calibration curves verified that the predicted prognosis from the nomogram have good consistency with the actual outcomes, which indicated that nomogram has good prediction performance (Supplementary Figure 4B–4D).

Gene set enrichment analysis (GSEA)

For the purpose of uncovering the biological pathways that differentiate patients of high-risk from low-risk group, we conducted GSEA on the expression data from GSE40967 and TCGA cohorts, respectively. The result of KEGG enrichment analysis revealed that Arrhythmogenic right ventricular cardiomyopathy, Citrate cycle, and Dilated cardiomyopathy were significantly enriched in high-risk group in GSE40967 cohort. Meanwhile, Huntington’s disease, Hypertrophic cardiomyopathy, Oxidative phosphorylation, Parkinson’s disease, and Pyruvate metabolism pathway were dramatically enriched in low-risk group in GSM40967 cohort (Supplementary Figure 5A). Moreover, genes from high-risk group in TCGA were enriched in Arrhythmogenic right ventricular cardiomyopathy, Cell cycle, Complement and coagulation cascade. However, Homologous recombination, Hypertrophic cardiomyopathy, Purine metabolism, Pyrimidine metabolism, and RNA degradation pathway were abundant among low-risk group in TCGA (Supplementary Figure 5B).

Identifying the differences of infiltration of immune cells

Next, we performed tumor microenvironment analysis using ESTIMATE algorithm and evaluated the proportion of immune cell infiltrated in tumor tissue using CIBERSORT. Stromal scores of patients in the high-risk group were significantly higher than those in the low-risk group in the training set (Supplementary Figure 6A). Accordingly, COAD patients in the high-stromal score group had worse prognosis compared to those in low-stromal score group (Supplementary Figure 6B). In addition, the infiltration of B cells memory, Plasma cells, T cells CD4 memory activated, and Neutrophils had statistical difference between high-risk and low-risk groups, respectively (Supplementary Figure 6C). Scatter plots (Supplementary Figure 6D–6H) showed that risk scores were positively associated with infiltration levels of B cells memory, Neutrophils, and Mast cells activated (R > 0, P < 0.05), and negatively associated with infiltration levels of Plasma cells, and T cells CD4 memory activated in the immune microenvironment (R < 0, P < 0.05).

Identifying the different expression level of immune checkpoint genes

The expressions level of CD27, CD48, HHLA2, ICOSLG, IDO1, NRP1, TNFRSF14, TNFRSF18, TNFRSF9, TNFSF14, TNFSF4, and TNFSF9 had statistical difference between high-risk and low-risk groups, especially for NRP1, TNFRSF14, and TNFSF4 (P < 0.001) (Supplementary Figure 7A). The Spearman correlation analysis indicated CD27 was markedly positively associated with CD48, IDO1, and TNFRSF9. Furthermore, CD48 was dramatically positively associated with DO1 and TNFRSF9. In addition, TNFRSF9 was positively correlated with IDO1. TNFSF4 was positively associated with NRP1 (Supplementary Figure 7B).

Response of patients with COAD to candidate drugs

We estimated the response of COAD patients in these two groups to candidate drugs. IC50 values of 60 compounds were significantly different in high-risk and the low-risk patients (Supplementary Table 1). We visualized the top eight drugs with the greatest difference in P values. AZD.0530, Bryostatin.1, CHIR.99021, Imatinib, LFM.A13, and CCT007093 were more sensitive to high-risk compared with low-risk patients (Supplementary Figure 8A–8F). However, PF.4708671 and EHT.1864 showed the opposite trend (Supplementary Figure 8G, 8H).

Knockdown of WFS1 molecule inhibits HCT116 and DLD-1 cells growth

We investigated the WFS1 expression in COAD and normal colon tissues in TCGA. Expression levels of WFS1 in COAD were dramatically higher than those in normal colon tissues (Figure 4A). Accordingly, patients with lower WFS1 expression had a longer OS (Figure 4B). These findings indicate that WFS1 has the potential to be a promising target for the therapy of colon cancer.

Figure 4. Knockdown of WFS1 inhibits COAD cell proliferation. (A) Box plots showed that WFS1 expression was significantly higher in COAD than in normal samples. (B) Survival curves showed that patients in the low WFS1 expression group had a better prognosis than those in the high WFS1 expression group. (C) Western blot analysis of WFS1 expression in HCT116 and DLD-1 cell lines. actin is internal control. (D) Cell viability of HCT116 or DLD-1 cells transfected with control or WFS1siRNAs was measured by CCK8 assay. (E) Clone formation assay of HCT116 or DLD-1 cells transfected with control or WFS1siRNAs. (F) EdU staining for cell proliferation transfected with control or WFS1siRNAs.

In order to evaluate the potential role of WFS1 in COAD, we inhibited the expression of WFS1 with two independent siRNAs in HCT116 and DLD-1 cell lines. These two siRNAs significantly decrease the protein expression of WFS1 (Figure 4C). WFS1 depletion markedly inhibited cell viability and clone formation capability of COAD cells (P < 0.05, Figure 4D–4F). In summary, these results indicate that WFS1 has a potential carcinogenic effect in COAD.