Transformative insights from transcriptome analysis of colorectal cancer patient tissues: identification of four key prognostic genes

Patient cohorts and sample preparations

A total of 49 matched tumour and normal FFPE sporadic colorectal cancer tissue samples (discovery dataset, Fig. 1: IIA1) obtained from the Ankara University, Faculty of Medicine, were used for transcriptome profiling. A second cohort of 64 tumours and matched normal CRC samples (validation dataset, Fig. 1: IIA2) was subsequently used to validate and test the reproducibility of the expression results. Table 1 presents the demographic and clinical features of participants in both the discovery and validation cohorts. Ethical approval for both the discovery (Ref: 153-4854) and validation (Ref: 15-200-19) cohorts was obtained from the Ankara University School of Medicine Clinical Research Ethics Committee. All participants provided informed consent, and the study was conducted in accordance with the Declaration of Helsinki.

The FFPE tissue samples were fixed in 10% neutral-buffered formalin immediately after collection, embedded in paraffin blocks. The blocks used in discovery and validation datasets were 2–8 and 4–10 years old, respectively. A pathologist reviewed H&E slides to select sections containing ~90% tumour cells. Four 8-μm sections per block were cut, placed on slides (2 per slide), one stained with H&E to mark tumour areas and stored at −20 °C until used for transcriptome profiling analysis.

FFPE RNA extraction and microarray transcriptome profiling

The marked tumour regions were macrodissected from serial sections and used for RNA extraction with the RNeasy FFPE kit (Qiagen, Hilden, Germany), following a modified deparaffinization step as previously described (Belder et al., 2016). Genomic DNA was removed using the gDNA Eliminator spin column included in the kit. RNA concentration and purity were assessed with a NanoDrop ND-1000 Spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA), and integrity was evaluated with an Agilent Bioanalyzer 2100 using the RNA 6000 Nano Assay (Agilent Technologies, Santa Clara, CA, USA).

According to the manufacturer’s instructions, 100 ng of input RNA was used to generate the amplified RiboSPIA product using the Ovation^® FFPE WTA System (NuGEN, San Carlos, CA, USA). The fragmented and labelled samples (3.5 μg) were then hybridised to Affymetrix GeneChip^® Human X3P Arrays (Affymetrix, Santa Clara, CA, USA) with rotation for 19 h at 45 °C. Subsequently, the arrays were washed and stained using user-prepared wash and staining reagents in a Fluidics Station 450 (Affymetrix) and scanned using a GeneChip Scanner 3000 (Affymetrix).

Microarray data pre-processing and analysis

The Affymetrix CEL files were pre-processed with Genomics Suite Version 6.6 (Partek Inc., St. Louis, MO, USA) using the robust multi-array analysis (RMA) algorithm, including background correction, quantile normalisation, logging of probes using base 2, and probe summarisation (Irizarry et al., 2003). A paired t-test was performed to determine which genes were differentially expressed between the tumour and matched normal samples, followed by Benjamini–Hochberg correction. The threshold for identifying differentially expressed genes (DEGs) was set as a false discovery rate (FDR) of <0.001 and |log₂FC| ≥1. Unsupervised hierarchical cluster analysis was performed in the Partek Genomic Suite using the Euclidean distance and average linkage algorithm. The gene expression data were deposited in the NCBI’s Gene Expression Omnibus Database (GEO: GSE271719).

Pathway and gene set enrichment analysis

The Database for Annotation, Visualization and Integrated Discovery (DAVID, https://davidbioinformatics.nih.gov/) (Huang da, Sherman & Lempicki, 2009) was used to perform KEGG pathway enrichment analysis of the DEGs. A threshold-adjusted p-value of ≤ 0.05 was considered statistically significant. Gene Set Enrichment Analysis (GSEA) software (version 4.1.0) (http://software.broadinstitute.org/gsea/index.jsp) was employed to validate the ECM-receptor interaction pathway enrichment in CRC using our discovery dataset. The number of gene set permutations was set at 1,000. NES (normalized enrichment score) >1 and FDR (false discovery rate) q-value ≤ 0.05 were defined as the significant cut-off criterias of enrichment.

Integrated mRNA microarray dataset analysis

Differential gene expression and subsequent pathway enrichment analysis were conducted using independent colorectal cancer gene expression datasets to support our microarray data and pathway analysis. Three independent datasets (GSE9348, GSE24550, and GSE21510) were downloaded from the NCBI GEO repository (https://www.ncbi.nlm.nih.gov/geo/) to identify DEGs between tumour and normal samples for the integrated microarray dataset analysis (Fig. 1: I). Detailed information about the datasets, including the platform and the number of samples, is provided in Table S1. A threshold of |log₂FC| ≥ 1 and an adjusted p-value of <0.05 were set to screen the DEGs. To analyse overlapping DEGs, we used an online Venn diagram tool (https://bioinfogp.cnb.csic.es/tools/venny/). KEGG pathway enrichment analysis was performed on the common DEGs using the DAVID bioinformatics resources.

miRNA microarray data set analysis

An integrative meta-analysis of publicly available miRNA microarray datasets was conducted to identify the miRNA-regulated genes involved in the ECM-receptor interaction pathway. Our goal was to narrow down the candidate gene list from our transcriptome analysis and identify key molecules relevant to CRC pathogenesis for further investigation.

We analysed a total of 161 samples (76 controls and 85 tumours) from the following datasets: NCBI GEO GSE68377 and GSE35982; Array Express E-MTAB-752, E-GEOD-35834, and E-MTAB-813. Detailed information on these datasets is provided in Table S2 (Fig. 1: III). Differential expression of miRNAs (DEMs) was assessed for each platform using cut-off values of adjusted p ≤ 0.05 and |FC| ≥1.5. DEMs were categorized as up-regulated or down-regulated for further analysis.

miRNA-mRNA target interactions

Potential target genes of DEMs were predicted using the miRTarBase (http://mirtarbase.mbc.nctu.edu.tw/php/index.php) and miRWalk2.0 (http://zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk2/) databases. Overlapping analysis was performed using a Venn diagram to obtain a comprehensive list of miRNA target genes. Common genes were then identified by intersecting the miRNA target genes with the DEGs from our discovery dataset for further analysis. miRNA-mRNA pairs showing inverse correlation expression trends were filtered out, and genes involved in the ECM-receptor interaction pathway were selected as key genes from miRNA-mRNA target data analysis.

Determination of candidate target genes

Functional pathway enrichment analysis of our CRC gene expression microarray data and integrated microarray analysis of independent GEO CRC gene expression datasets revealed that the DEGs were primarily implicated in the ECM-receptor interaction pathway. We, therefore, focused on this pathway using an integrated approach to identify candidate genes more accurately by combining genes in the ECM pathway from our microarray study, the GEO mRNA datasets, and critical genes mentioned in the miRNA meta-analysis. The Venn diagram tool (http://bioinfogp.cnb.csic.es/tools/venny) was used to analyse overlapping genes. Five genes (THBS2, FN1, COL1A1, COL5A1, CD44) intersecting in the ECM-receptor interaction pathway were identified as candidate genes for further analysis. Univariate Cox regression analysis was applied to examine their association and patient overall survival (OS) using the GSE17536 and TCGA datasets (Fig. 1: II-2). Four of these genes (THBS2, FN1, COL1A1, and COL5A1), which showed a statistically significant association (p < 0.05) with patient OS, were identified as the critical target genes for subsequent analysis.

Quantitative real-time PCR (qRT-PCR)

To validate expression changes in the four target genes, qRT-PCR was conducted using our discovery dataset of RNA samples (n = 49). Retrospective validation was also performed using the validation cohort, which consisted of independently collected FFPE colorectal cancer samples and matched controls (n = 64). cDNA was synthesized from 0.5 μg total RNA using the Transcriptor First Strand cDNA Synthesis Kit (Roche, Basel, Switzerland) with random primers in a total volume of 20 μl. The reverse transcription reaction was carried out at 55 °C for 30 min, followed by 85 °C for 5 min. qRT-PCR was performed in triplicate using the SYBR Green I Master Kit (Roche, Basel, Switzerland) on a LightCycler480 instrument (Roche, Basel, Switzerland), following the manufacturer’s instructions. Amplification was performed at 95 °C/180 s followed by 40 cycles at 95 °C/30 s and 60 °C/30 s, 72 °C/30 s. Gene expressions were normalized to the geometric mean of HPRT and B2M mRNA levels. Ct values were acquired using the LightCycler^® 480 software, and data were analyzed using the ΔΔCT method (Livak & Schmittgen, 2001). Primer sequences are listed in Table S3.

Verification of expression patterns and prognostic value of four candidate key genes in publicly available CRC datasets

External validation cohorts from the GEO and TCGA datasets were used to confirm the reproducibility of the overexpression of the four target genes in CRC (GSE44076, GSE24551, GSE39582, GSE32323, TCGA COAD-READ). Associations between elevated gene expression and advanced tumour progression (examined in GSE14333, GSE37892, GSE40967, GSE17538), as well as recurrence (examined in GSE40967, GSE17538), were also assessed, as shown in Fig. 1 (IIC1a). Statistical significance was set at p < 0.05. Additionally, the prognostic value of these genes in CRC (TCGA, GSE17538) was evaluated using Kaplan-Meier analysis, with p-values calculated via the log-rank test (GraphPad Prism 6.6). A log-rank p-value < 0.05 was considered statistically significant. Differences in gene expression for THBS2, FN1, COL1A1, and COL5A1 across different cancers were assessed using Oncomine (https://www.oncomine.org), with threshold parameters of p < 0.05, fold-change >2, and gene rank within the top 10%. To further investigate the potential biological roles of candidate prognostic genes in CRC, gene set enrichment analysis was performed using independent datasets (GSE17536 and GSE39582), as previously outlined. For each gene, patients were dichotomized into high- and low-expression groups based on the median expression value, and GSEA was carried out separately using KEGG pathway gene sets (c2.cp.kegg.v7.4.symbols.gmt).

Correlation analysis

The discovery dataset was subjected to Spearman correlation analysis in GraphPad Prism 6.6 to evaluate the interactions among the key genes at the transcription level. Correlation coefficients were calculated to assess the strength and direction of the correlations. Correlations among the genes were validated using external datasets, including GSE17536, GSE40967, and TCGA datasets.

Protein-protein interaction (PPI) network analysis

The online resource Search Tool for the Retrieval of Interacting Genes (STRING) database (http://stringdb.org/) was utilised to explore relationships for our crucial target proteins. After importing the candidate genes (THBS2, FN1, COL1A1 and COL5A1) into STRING, a confidence score of >0.4 was applied as the threshold for the PPI network information.

Determining the prognostic risk score for patient risk stratification

A multivariate Cox regression analysis was carried out on the discovery dataset with the log2 normalised expression values of the four candidate genes (THBS2, FN1, COL1A1, and COL5A1) and disease stage as the variables for each patient (p = 0.025). A combined gene score (∑(βi × Expi)) was calculated, where Expi represented the expression level of prognostic gene i, and βi was the regression coefficient for gene i obtained from the multivariate Cox analysis (Fig. 1: IIC2). The cohort was proportionally split into Low- and High- score values across the mean. The method was similarly applied to the GSE datasets and overall survival of the patients with high- and how- risk scores was compared by Kaplan-Meier survival and log-rank (Mantel-Cox) test to assess significance (Fig. 1: IIC2a).

Statistical analyses

All statistical analyses were conducted using GraphPad Prism 6.6 (GraphPad Software, Inc., San Diego, CA, United States), IBM SPSS^® Statistics (version 20.0; IBM Company, New York, NY, USA) software package, and R v. 3.4.3. Univariate and multivariate Cox regression analyses were conducted to evaluate the predictive value of the target genes, risk scores, and clinicopathological features in CRC. All tests were two-tailed, and p-values of <0.05 were considered statistically significant.

Gene expression profile associated with sporadic colorectal cancer

Affymetrix GeneChip Human Genome U133 X3P expression arrays were used to identify differences in gene expression between the 49 cases of sporadic CRC FFPE tissue and their paired normal tissue samples. Applying the criteria of an adjusted p-value of <0.001 and a fold change (|log₂FC|) of ≥1, a total of 845 DEGs were identified, comprising 320 up-regulated genes and 525 down-regulated genes (Table 2, Table S4). The distribution of the DEGs was visualised using the volcano plot in Fig. 2A, in which genes meeting the criteria of a p-value < 0.001 and |log₂FC| ≥ 1 are represented by green and red dots, respectively. The PCA plot (Fig. 2B) shows a distinct separation between the tumour and control samples. Additionally, hierarchical clustering of the 845 DEGs across the 49 paired samples (Fig. 2C) demonstrated that DEG expression patterns effectively distinguished tumour tissues from normal counterparts.

Table 3 lists the top 10 upregulated and downregulated DEGs. Among these, THBS2 (Thrombospondin-2), a member of the thrombospondin family, exhibited the highest upregulation (fold change = 20.20) in CRC compared to matched controls. The most downregulated gene (fold change = −39.40) was SLC26A3 (solute carrier family 26, member 3), which functions as an anion exchanger.

DEGs are primarily involved in the ECM-receptor interaction pathway

To explore the functional relevance of the 845 DEGs, KEGG pathway enrichment analysis was performed using the DAVID (v6.8) online tool and identified 13 enriched pathways (FDR < 0.05). The ECM-receptor interaction pathway was the most dysregulated in CRC, followed by pathways in cancer, focal adhesion, and PI3K-Akt signalling (Fig. 3A). Most upregulated genes were linked to focal adhesion and ECM-receptor interaction, while downregulated genes were mainly associated with metabolic processes (Fig. 3B). GSEA further confirmed significant enrichment of ECM-receptor signalling in tumour samples (Fig. 3C).

To support our findings, three independent CRC microarray datasets (GSE9348, GSE24550, GSE21510) from the NCBI GEO repository were analysed (Fig. 1: I)., resulting in 4,874, 5,634, and 7,853 DEGs (adjusted p < 0.05 and |log2FC| > 1), respectively (Table S5). The volcano plots for the DEGs in all three datasets are displayed in Figs. S1A–S1C. A Venn diagram revealed 971 common DEGs across these datasets (Fig. S2). Pathway analysis of these common DEGs highlighted four significant pathways in CRC: cell cycle, DNA replication, ECM-receptor interaction, and small-cell lung cancer (Fig. 3D). Considering both the discovery dataset and the external GEO CRC datasets for the pathway enrichment analysis, only the ECM-receptor interaction pathway emerged as a common and significantly enriched pathway. Given its crucial role in CRC, this pathway was selected as the starting point for identifying target genes and conducting further analyses.

THBS2, FN1, COL1A1, and COL5A1 were identified as candidate target genes

We conducted an integrative meta-analysis of selected publicly available miRNA microarray datasets to narrow down the enriched ECM-receptor interaction pathway gene list derived from the discovery dataset and integrated microarray analysis (Fig. 1: III). GSE68377 and GSE35982 (GEO datasets), E-MTAB-752, E-GEOD-35834, and E-MTAB-813 (Array Express datasets) (Table S2), were analysed separately. Differentially expressed miRNAs (DEMs) were identified for each platform (adj. p ≤ 0.05, |FC| ≥ 1.5) and classified as upregulated or downregulated miRNAs for further analysis. Overall, 110 miRNAs from the analysis of the Affymetrix E-MTAB-752, E-GEOD-35834, and GSE68377 datasets showed significant expression differences, with 63 upregulated and 47 downregulated. Additionally, we obtained 90 DEMs, including 75 upregulated and 15 downregulated miRNAs, from the analysis of the Agilent GSE35982 and E-MTAB-813 datasets (Table S6). Forty-two miRNAs were common across both platforms (Fig. S3).

Target genes for these forty-two miRNAs were predicted using miRTarBase and miRWalk2.0, and a Venn diagram identified overlapping genes with the DEGs from our discovery dataset (Table S7). miRNA-mRNA pairs showing an inverse correlation in expression trends were filtered out. Genes involved in the ECM-receptor interaction pathway and supported by miRNA analysis were prioritised (Fig. 4A). By intersecting ECM pathway genes from our microarray study, integrated analysis of GEO mRNA datasets, and miRNA analysis, we identified five candidate genes (THBS2, FN1, COL1A1, COL5A1, CD44) for further investigation (Fig. 1: II-1; Fig. 4B).

To assess their prognostic significance, a univariate Cox regression analysis was performed to examine the association between these five genes and patient OS using the GSE17536 dataset (Fig. 4C) and TCGA dataset (Fig. 4D). Among the five genes, four (THBS2, FN1, COL1A1, and COL5A1) showed a statistically significant association (p < 0.05) with patient OS and were therefore selected as essential target genes for subsequent analysis (Fig. 1: II-2). All four candidate genes were significantly upregulated (p < 0.05) in tumour samples compared to normal samples in both the discovery microarray dataset and the three independent GEO CRC datasets used for the integrated microarray analysis (Fig. 4E). The dysregulation and subsequent upregulation of these genes likely contribute to CRC tumorigenesis and prognosis.

THBS2, FN1, COL1A1, and COL5A1 gene expression changes in CRC were technically and externally verified

All four candidate genes (THBS2, FN1, COL1A1, COL5A1) were found to be upregulated in tumour samples compared to matched normal samples in our microarray analysis (Fig. 5A). To ensure the reliability of the differential expression of these four genes identified by our microarray profiling, we used qRT-PCR, the gold standard for confirming gene expression. The qRT-PCR results confirmed the upregulation of these genes in the same set of samples (p < 0.05) (Fig. 1: IIA1; Fig. 5B).

To assess the biomarker potential of these target genes, it is essential to validate their expression changes or upregulation in different and independent tumour samples. Therefore, we evaluated the robustness of the gene panel in four independent GEO colorectal cancer datasets (GSE44076, GSE24551, GSE39582, and GSE32323) and the TCGA COAD-READ dataset (Fig 1: IIA3). All four candidate genes were significantly upregulated in CRC patients compared to healthy controls, except for FN1 in the TCGA cohort (p < 0.05) (Fig. 5C). This validation across multiple and different cohorts highlights their potential importance in colorectal carcinogenesis.

Overexpression of THBS2, FN1, COL1A1, and COL5A1 in CRC tissues was further confirmed biologically

To validate the expression profile of the four candidate prognostic genes, a new set of independent cohort of 64 CRC tumours and matched normal samples was collected (Fig. 1 IIA2). The qRT-PCR results demonstrated that THBS2, FN1, COL1A1, and COL5A1 were significantly overexpressed in CRC tissues compared to their matched controls (p < 0.05) (Fig. 6A). Confirming the previously identified up-regulation in independent clinical samples highlights the validity and reliability of the results obtained from the microarray analysis. Furthermore, utilising the Oncomine database, we conducted a pan-cancer differential gene expression analysis for THBS2, FN1, COL1A1, and COL5A1 (Fig. 1 IIA3). As shown in Fig. 6B, our candidate genes were significantly upregulated in many cancer types compared to normal tissues. In particular, THBS2, COL1A1, and COL5A1 genes are upregulated in most cancer types, except for bladder, kidney, melanoma, and ovarian cancer. These results suggest that our essential genes may predominantly exhibit oncogenic characteristics in multiple types of cancer, particularly in CRC. To further investigate the functional roles of these genes in CRC progression, GSEA was performed using the GSE17536 and GSE39582 datasets. Patients were stratified into high- and low-expression groups based on the median expression values. Notably, the high-expression groups exhibited significant enrichment in several pathways, including focal adhesion, actin cytoskeleton regulation, and cell adhesion molecules (CAMs), in addition to the ECM-receptor interaction pathway. Furthermore, the KEGG mTOR signalling pathway was also enriched, suggesting a potential link between these genes and mTOR signalling in CRC progression. These findings offer a potential mechanistic link between the candidate genes and CRC progression (Tables S8–S15).

High expression of THBS2, FN1, COL1A1 and COL5A1 is associated with advanced stage and worse prognosis in CRC patients

We analysed the mRNA expression of THBS2, FN1, COL1A1, and COL5A1 across CRC stages using four GEO datasets (GSE14333, GSE17538, GSE40967, GSE37892) (Fig. 1 IIC1a). A significant increase in expression was observed from early to late stages (p < 0.05), suggesting their role in CRC progression (Figs. 7A–7D). We then examined whether the overexpression of the four genes was associated with the prognosis of the CRC patients in the GSE17538 and TCGA CRC datasets (Fig. 1 IIC1b). The results indicated that high expression of THBS2 (p = 0.0016, p = 0.018), FN1 (p < 0.0001, p = 0.0532), COL1A1 (p = 0.018, p = 0.0239), and COL5A1 (p = 0.0001, p = 0.067) was significantly associated with poor OS in GSE17538 and TCGA (Figs. 7E, 7F). Additionally, these genes were significantly upregulated in the recurrent CRC group compared to the non-recurrent group, further highlighting their prognostic relevance. (p < 0.05) (Figs. 7G, 7H). Spearman’s correlation analysis results revealed clear positive correlations (r > 0.6, p < 0.001) between genes in colorectal cancer in the discovery dataset (Fig. 1: IIB; Fig. 8A), which were confirmed in the independent GEO and TCGA datasets (GSE17536: r > 0.6, p < 0.001; GSE40967: r > 0.6, p < 0.001; TCGA COAD-READ: r > 0.8, p < 0.001) (Figs. 8B–8D). A heatmap further visualized these co-expressions, highlighting their collaborative role in CRC (Figs. 8E–8H). Notably, the strongest positive correlation across all datasets was observed between COL1A1 and COL5A1 (discovery dataset: r = 0.86, p < 0.0001; GSE17536: r = 0.95, p < 0.0001; GSE40967: r = 0.96, p < 0.0001; TCGA COAD-READ: r = 0.97, p < 0.0001). Finally, a PPI network analysis in the STRING confirmed a highly significant interactions, highlighting the biological interconnectivity of these genes (PPI enrichment p = 1.15e−05) (Fig. 8I).

A four-gene signature-based risk assessment model was established to predict the prognosis of colorectal cancer

Given the limitations of using a single gene to predict prognosis in a complex disease like cancer, we evaluated the combined prognostic value of the four genes. Since tumour stage is a crucial factor in patient survival (Brierley et al., 2019), we developed a prognostic risk score model that incorporates the expression levels and regression coefficients of these four genes alongside tumour stage to predict overall survival (OS) in our discovery dataset (Fig 1: IIC2). We then calculated the risk score for each patient and divided the CRC patients into high-risk and low-risk groups using the median cut-off of the risk score. Kaplan-Meier plots indicated significant differences in OS between the two groups. Patients in the high-risk group showed significantly poorer OS compared to those in the low-risk group (p = 0.0443; Fig. 9A).

To assess the robustness and effectiveness of our risk score model and verify its reproducibility, we used two independent GEO datasets (GSE17536, n = 145 and GSE40967, n = 560) with OS and gene expression information of the CRC patients (Fig 1. IIC2a). As shown by the Kaplan–Meier OS curves in Figs. 9B and 9C, patients in the high-risk group had shorter overall survival than those in the low-risk group (GSE17536: p < 0.0001; GSE40967: p = 0.0446) (Fig. 9B; Fig. 9C). Gene expression for all four prognostic genes was significantly higher in the high-risk group (p < 0.0001) (Figs. 9D–9I). Overall, the results confirmed that our risk score model, based on the four-gene scoring system and incorporating disease stage, could provide an accurate prognosis for CRC patients.