RNA-seq data from READ patients at the fragments per kilobase million gene level with clinical information and produced by Illumina HiSeq 2000 RNA sequencing platform were downloaded as training datasets from TCGA websites using the genomic data commons data transfer tool (https://gdc-portal.nci.nih.gov/) before 19 October 2021. DEGs were detected using 158 READ tissues and nine normal controls. Supplementary File 1 contains the sample name obtained from TCGA. Simultaneously, the validation data set for GSE56699 was downloaded from the GEO database (http://www.ncbi.nlm.nih.gov/geo/) using the Illumina HumanHT-12 WG-DASL V4.0 R2 expression bead chip. We included the samples with genome-wide expression profile data and clinical prognostic information. Overall, a total of 61 samples were included.

The flowchart of this study is depicted in Supplementary Figure S1. The DEGs in the TCGA training data set were identified using the limma package version (v) 3.34.7 of R v3.6.1 (https://bioconductor.org/packages/release/bioc/html/limma.html) with a false discovery rate (FDR) of 0.05 and |log2 (fold change) | > 1. Then, a two-way hierarchical clustering analysis was performed using heatmap v1.0.8 in R v3.6.1 (https://cran.r-project.org/web/packages/pheatmap/index.html) on the DEG expression levels obtained in the training data set.

Immune-related DEGs were downloaded for further analysis from the AmiGO 2 (http://amigo.geneontology.org/amigo) and KEGG databases (https://www.kegg.jp/). The function of these DEGs was then determined using the GO biology process and KEGG signal pathway enrichment using DAVID v6.8 (https://david.ncifcrf.gov/), with an FDR threshold of <0.05.

Univariate and multivariate Cox regression analyses were used to identify independent DEGs associated with overall survival (OS) using the survival package v2.41-1 of R v3.6.1 (http://bioconductor.org/packages/survivalr/). Significant DEGs were identified using a log-rank p-value threshold of <0.05.

The optimal DEGs combination was then determined using the LASSO Cox regression model in R v3.6.1 of penalized package v0.9.50 (https://cran.r-project.org/web/packages/penalized/index.html). The screening model’s optimal parameter “lambda” is obtained through 1,000 cross-validation likelihood algorithm calculation cycles. The following prognostic risk model was constructed using the prognostic coefficients from the LASSO Cox regression model and the DEG expression level:

Prognostic risk score = ∑β_DEGs × Exp_DEGs

Here, β_DEGs denote the DEGs coefficient derived from the LASSO Cox regression model, whereas ExpDEGs denote the target DEGs expression level in the training dataset.

Each sample in the TCGA training and GSE56699 validation datasets was analyzed and assigned a prognostic risk score. The median value was used to classify the samples as high or low risk. In the TCGA training and GSE56699 validation datasets, the correlation between actual survival prognosis and that predicted by the prognostic risk model was evaluated using the Kaplan–Meier curve method in R v3.6.1 with the survival package v2.41-1GSE56699.

The gene modules of the Tumor Immune Estimation Resource (TIMER; cistrome. shinyapps.io/timer/) were used to investigate the correlation between the expression of prognostic DEGs, and the abundance of six immune infiltrates (B cell, CD4⁺ T cell, CD8⁺ T cell, neutrophil, macrophage, and dendritic cell). A heatmap was generated using the partial correlation index for each DEG in each immune infiltrate and several scatterplots.

Independent clinical factors such as age (years), gender, pathologic M (M0/M1/-), pathologic N (N0/N1/N2/-), pathologic T (T1/T2/T3/T4/-), pathologic stage (I/II/III/IV/-), history of colon polyps, lymphatic invasion, radiotherapy, prognostic model, death, and OS time (months) screened patients with READ in the TCGA training data set using univariate and multivariate Cox regression analyses with log-rank p-value <0.05 as the threshold. The analysis was stratified by age (>65 and ≤65 years of age) and pathologic stage (N0, N1, and N2).

In order to evaluate the prognostic risk prediction model, stratified analysis was performed on the samples that were divided into different sample comparing groups. The nomogram displaying 3-year and 5-year OS was constructed to further reveal the correlation between independent factors and actual prognosis using the rms package v5.1-2 (https://cran.r-project.org/web/packages/rms/index.html) in R v3.6.1. survcomp version 1.34.0 (http://www.bioconductor.org/packages/release/bioc/html/survcomp.html) was used to calculate the C-index in R v3.6.1.

There were 1,772 DEGs identified, with 768 upregulated and 1,004 downregulated across all genes (Figure 1A). From the heatmap, we observed the clustering of tumor and control samples clustered separately, ensuring the reliability of the original data (Figure 1B).

Simultaneously, we downloaded 3,020 and 817 unique immune-related genes from the AmiGO two and KEGG databases, respectively, leaving 3,255 union immune-related genes. When TCGA DEGs were compared to immune-related genes, 326 immune-related DEGs were retained for further investigation (Figure 1C). We provided detailed information on 326 DEGs (log2 FC, p-value, and FDR) in Supplementary File 2.

Additionally, we examined the functions of 326 DEGs using GO and KEGG analyses, identifying 36 BP and 22 KEGG under a <0.05 p-value (Figure 2; Table 1). The DEGs were mainly enriched in the GO term of the immune response, inflammatory response, chemokine-mediated signaling pathway, chemotaxis, innate immune response, positive regulation of the ERK1 and ERK2 cascades, response to lipopolysaccharide, adaptive immune response, positive regulation of transcription from the RNA polymerase II promoter, and positive regulation of cell proliferation (Figure 2A). Additionally, KEGG pathways also involved cytokine–cytokine receptor interactions; the chemokine signaling, Ras signaling, complement, coagulation cascades, cancer, Fc epsilon RI signaling pathways; natural killer cell-mediated cytotoxicity; leukocyte transendothelial migration; B cell receptor signaling pathway, and serotonergic synapse (Figure 2B).

The overall 326 DEGs were subjected to univariate Cox regression analysis, and 41 were identified as prognostic-related DEGs. Following the multivariate Cox regression analysis, a total of 22 DEGs remained. After that, nine optimal DEG combinations related to immunity were identified, including galectin 9C [LGALS9C; hazard ratio (HR) = 0.930, 95% confidence interval (CI) = 0.850–0.973], coagulation factor XIII A chain (F13A1; HR = 1.012, 95% CI = 1.004–1.554), ADAM-like decysin 1 (ADAMDEC1; HR = 0.987, 95% CI = 0.854–0.992), macrophage receptor with collagenous structure (MARCO; HR = 1.008, 95% CI = 1.002–1.504), L3MBTL histone methyl-lysine binding protein 1 (L3MBTL1; HR = 0.883, 95% CI = 0.803–0.970), solute carrier family 7 member 11 (SLC7A11; HR = 0.881, 95% CI = 0.814–0.954), UL16 binding protein 3 (ULBP3; HR = 0.884, 95% CI = 0.800–0.976), complement component 4 binding protein alpha (C4BPA; HR = 1.019, 95% CI = 1.001–1.458), and Cbp/p300 interacting transactivator with Glu/Asp rich carboxy-terminal domain 1 (CITED1; HR = 1.065, 95% CI = 1.008–1.125) (Table 2). The prognostic risk model was constructed based on LGALS9C, F13A1, ADAMDEC1, MARCO, L3MBTL1, SLC7A11, ULBP3, C4BPA, and CITED1 expression levels. Then, using the constructed model, the prognostic risk score of each sample was calculated as: prognostic risk score = (−0.2332) × ExpLGALS9C + 0.1177 × ExpF13A1 + (−0.1821) × Exp _ADAMDEC1 + 0.1287 × Exp_MARCO+ (−0.4846) × Exp_L3MBTL1+ (−0.6165) × Exp_SLC7A11+ (−0.5122) × Exp_ULBP3+ 0.0929 × Exp_C4BPA+ 0.3338 × Exp_CITED1.

The highly expressed genes in tumor cells typically have positive associations with tumor purity. As such, we examined the association between prognosis DEG expression and six immune infiltrates (Supplementary Figure S2). ADAMDEC1 correlated positively with six immune infiltrates, particularly B cell and dendritic cell infiltrates (partial correlation = 0.421, p = 2.52e-07; Figures 3A,B). The results indicated that increased ADAMDEC1 expression was associated with a higher purity of READ tumor cells in B and dendritic cells. F13A1 also had positive correlation with six immune infiltrates, most notably macrophages (partial correlation = 0.423, p = 2.06e-07; Figures 3A,C) and dendritic cells (partial correlation = 0.598, p = 7.35e-15; Figures 3A,C). Apart from the B cells, LGALS9C had a negative correlation with immune infiltrates (Figure 3A).

The prognostic risk model was evaluated by classifying samples into high- and low-risk groups in the TCGA training and GSE56699 validation datasets. Then, in two data sets, the high- and low-risk groups were compared according to their prognostic risk models to those classified according to their actual status. Both TCGA training [HR = 9.989 (3.382–29.50), p = 3.373e-07; Figure 4A] and GSE56699 validation datasets [HR = 8.428 (1.074–66.12), p = 8.077e-03; Figure 4B] revealed significant differences between high- and low-risk groups. Additionally, the prognostic risk model demonstrated a strong correlation with the actual situation in both the TCGA training [AUC = 0.906 (0.908, 0.893); Figure 4C] and GSE56699 validation datasets [AUC = 0.836 (0.860, 0.727); Figure 4D].

Univariable Cox regression analysis was used to eliminate independent clinical prognosis factors such as age (p = 7.97e-05), pathologic M (p = 2.30e-03), pathologic N (p = 4.61e-03), pathologic T (p = 3.07e-02), pathologic stage (p = 1.03e-03), and prognostic model (p = 3.37e-07) (Table 3). After multivariable Cox regression analysis, only three factors were retained for further investigation: age (p = 1.27e-03), pathologic stage (p = 4.91e-02), and prognostic risk model (p = 4.88e-03; Table 3).

The stratified analysis based on age (>65 and ≤65 years old) revealed that patients aged 65 years and older had a significantly lower survival rate [HR = 3.812 (1.537–9.449), p = 1.409e-03; Figure 5A]. According to the risk score proposed above, the samples in each subgroup were divided into low- and high-risk groups. In patients ≤65 years of age, those with high-risk scores had a significantly shorter OS time than those with low-risk scores [HR = 5.522 (1.926–22.43), p = 2.846e-04; Figure 5B]. In patients over the age of 65, those with a high-risk score had a significantly shorter OS time than those with a low-risk score [HR = 6.190 (2.073–18.49), p = 1.595e-04; Figure 5C].

The stratified analysis based on pathologic stage (N0, N1, and N2) revealed that patients with a higher pathologic stage [N1, N2; HR = 1.886 (1.195–2.976), p = 4.609e-03; Figure 5D] had a significantly lower survival rate than those who were (N0). The patients with low-risk scores at pathologic N0 have a significantly shorter OS time than patients with high-risk scores [HR = 4.870 (1.932–25.44), p = 3.854e-02; Figure 5E]. There were no significant differences in OS times between pathologic N1 and N2 patients classified as high- or low-risk (Figures 5F,G). The results indicated that a lower age (≤65 years) and a more advanced pathologic stage are associated with a better prognosis for READ patients, consistent with their actual status.

The survival nomogram model analysis of TCGA training dataset samples revealed that age and mRNA prognostic factors were the most significant predictors of 3- and 5-year survival (Figure 6A). The 3-year (C-index = 0.759) and 5-year (C-index = 0.724) survival probabilities predicted by the model were generally consistent with actual survival rates (Figure 6B).