COADREADx: A comprehensive algorithmic dissection of colorectal cancer unravels salient biomarkers and actionable insights into its discrete progression

Data preprocessing

Normalized and log₂-transformed Illumina HiSeq RNA-Seq gene expression data for Colorectal Adenocarcinoma (COADREAD) processed by the RSEM pipeline (Li & Dewey, 2011) were obtained from TCGA via the http://firebrowse.org/ portal (accessed 06-01-2019) (Broad Institute TCGA Genome Data Analysis Center, 2016). The patient barcode (uuid) of each sample encoded in the variable called ‘Hybridization REF’ was parsed and used to annotate the controls and cancer samples. To annotate the stage information of the cancer samples, we obtained the corresponding clinical dataset from http://firebrowse.org/ and merged the clinical data with the expression data by matching the “Hybridization REF” in the expression data with the aliquot barcode identifier in the clinical data. The cancer staging is encoded in the attribute “pathologic_stage” of the clinical data. The sub-stages (A,B,C) were collapsed into the parent stage, resulting in four stages of interest (I, II, III, IV). We retained a handful of clinical variables related to demographic features, namely age, sex, height, weight, and vital status. Using this merged dataset, we filtered out genes that showed little change in expression across all samples (defined as σ <1). We also removed cancer samples that were missing stage annotation (value ‘NA’ in the “pathologic stage”) from our analysis. Data pre-processing was done with R v4.2.3 (www.r-project.org) and the final dataset was processed through voom function in R limma v3.54.2 to prepare for linear modeling (Law et al., 2014).

Linear modelling

Linear modeling of expression across cancer stages relative to the baseline expression (i.e., in normal tissue controls) was performed for each gene using R limma v3.54.2 (Ritchie et al., 2015). The following linear model was fit for each gene’s expression based on the design matrix shown in Fig. 2A: (1)${\displaystyle y=\alpha +{\beta }_1{x}_1+{\beta }_2{x}_2+{\beta }_3{x}_3+{\beta }_4{x}_4}$ where the independent variables are indicator variables of the sample’s stage, the intercept α is the baseline expression estimated from the controls, and β_i are the estimated stagewise log fold-change (lfc) coefficients relative to controls. The linear model was subjected to empirical Bayes adjustment to obtain moderated t-statistics (McCarthy & Smyth, 2009). To account for multiple hypothesis testing and the false discovery rate, the p-values of the F-statistic of the linear fit were adjusted using the method of Hochberg & Benjamini (1990). The linear trends across cancer stages for the top significant genes were visualized using boxplots to ascertain the regulation status of the gene relative to the control.

Pairwise contrasts

To perform contrasts, a slightly modified design matrix shown in Fig. 2B was used, which would give rise to the following linear model of expression for each gene: (2)${\displaystyle y={\beta }_0{x}_0+{\beta }_1{x}_1+{\beta }_2{x}_2+{\beta }_3{x}_3+{\beta }_4{x}_4}$ where the controls themselves constitute one of the indicator variables, and the β_i are all coefficients estimated only from the corresponding samples. Our first contrast of interest, between each stage and the control, was achieved using the contrast matrix shown in Table 1. Four contrasts were obtained, one for each stage vs control. A threshold of —lfc—>2 was applied to each contrast to identify genes differentially expressed with respect to the control. Genes could be differentially expressed in any combination of the stages. In the first pass, we identified genes whose —lfc—>2 for any stage. For the genes that passed, we identified the stage that showed the highest —lfc— for each gene and assigned the gene as specific to that stage for the rest of our analysis.

Significance analysis

We applied four-pronged criteria to establish the salience of the stage-specific differentially expressed genes. (i) Adj. p-value of the contrast with respect to the control < 0.001 (ii)–(iv) P-value of the contrast with respect to other stages < 0.05. Six such contrasts are possible.

To obtain the above p-values (ii)–(iv), we used the contrast matrix shown in Table 2, which was supplied as an argument to the contrastsFit function in limma.

To deal with any sparsity of progression-significant genes salient to any stage,we defined the “pval_pdt” of a given gene in a certain stage as the product of the p_values of its expression contrast in that stage vs each of the other stages (e.g., pval_pdt of gene x in stage 1 is (pval(gene x in st1 vs st2))*(pval(gene x in st1 vs st3))*(pval(gene x in st1 vs st4)).

Monotonic expression

The linear model in Eq. (1) would not be sufficient to identify genes with an monotonic, trend of expression in sync with disease progression, which could uncover stage-agnostic expression of progression-significant driver genes. Towards this end, we used a model of gene expression where the cancer stage was treated as a numeric variable: (3)${\displaystyle y=aX+b}$ where X takes a value in (0,1,2,3,4) corresponding to the sample stage: (control, I, II, III, IV), respectively. It was noted the mean gene expression could show the following patterns of monotonic expression across cancer stages:

(i) monotonic upregulation, where mean expression follows:

control <I <II <III <IV.

(ii) monotonic downregulation, where mean expression follows:

control >I >II >III >IV.

The sets of genes conforming to either (i) or (ii) were identified to yield monotonically upregulated and monotonically downregulated genes. These two sets were merged, and the final set of genes was evaluated using the adj. p-values from the model given by Eq. (3) to yield genes with significant monotonic patterns of expression.

Models for cancer screening and prognosis

Benchmarking

Principal component analysis (PCA) was performed using prcomp in R. We used the rand function to choose 100 random genes. In order to visualize significant outlier genes with a large effect size, volcano plots were obtained by plotting the (−log₁₀)-transformed p-value vs. the log fold-change of gene expression. Heat maps of significant stage-salient differentially expressed genes were visualized using R pheatmap v1.0.12 and clustered using R hclust function. Novelty of the identified stage-salient genes was ascertained by screening against curated databases, including the Cancer Gene Census (CGC at https://cancer.sanger.ac.uk/cosmic; accessed 01-12-23) (Futreal et al., 2004), Network of Cancer Genes NCG7.0 (accessed 01-12-23) (Repana et al., 2019), and the Clinical Trials Registry (http://www.clinicaltrials.gov; accessed 01-12-23). STRINGdb was used to translate the findings into network-level insights (Szklarczyk et al., 2021). To perform immuno-cyte infiltration analysis, we used Cibersort and estimated the proportion of tumor-infiltrating immune cells in TCGA COADREAD samples based on gene expression signatures (Amin et al., 2017; Newman et al., 2019). Cibersort’s inbuilt LM22 signature estimated the proportion of 22 standard immune cell types; setting the number of permutations to 100 allowed the calculation of sample-wise statistical significance with respect to the estimated values. The immuno-cyte patterns of significant samples were analyzed to provide a snapshot of immune ecotypes at play in significant tumor and normal samples, which would increase our basic understanding of colorectal cancer pathologies and advance rational therapies. The cell-type correlation matrix computed from the proportions of cell-types across significant samples was used to identify substantial co-occurrence patterns. The relative abundance of immunocytes between tumor and normal samples was compared to pinpoint significant differentially elevated or depressed tumor-infiltrating immune cells.

Normals-augmented validation

To examine any negative results with the inclusion of more controls in teasing out stage-specific markers, we augmented the dataset using RNAseqDB, which added 339 normal colorectal samples. We noted that the RNAseqDB preprocessing protocol eliminated non-coding transcripts from consideration, ignoring possible expression salience of non-coding RNA biomarkers like HOTAIR. Application of our whole protocol to this controls-augmented dataset yielded a linear model, 1925 stage-specific DEGs (755 stage-I, 418 stage-II, 163 stage-III and 589 stage-IV), and 105 stage-salient markers (40 stage-I, 6 stage-II, 2 stage-III and 57 stage-IV). These are presented in File S12. We found a substantial consensus of stage-salient genes between the two datasets, with 70 biomarkers in common (Table 8; highlighted in File S12). Notably six of the top stage-I salient genes and nine of the top stage-IV salient genes were identified as salient to the respective stages with the normals-augmented dataset as well, providing robust validation for these biomarkers.

In addition, we identified a colonic cancer dataset with stage-annotation from the Gene Expression Omnibus (GEO) database (Barret et al., 2013), namely GSE39582, provided by the Carte d’identité des tumeurs, Ligue Nationale contre le Cancer, France (accessed 12-08-2022) (Marisa et al., 2013). The dataset had a large number of stage-II (271) and stage-III samples (210), relative to stage-I (38) and stage-IV (60) samples. However, only two normal samples were available, so the dataset was augmented with 308 normal colonic tissue samples from the GTEx (accessed 12-10-2022). The augmented dataset was subjected to batch correction using ComBat (Leek et al., 2012), and antilog₂ was taken to obtain the necessary counts for input to voom and the protocol described in the Methods was applied. The results are presented in File S13. Five stage-IV salient genes, namely CYP24A1, FGF19, NKD1, COL9A3, and EDNRA are common to both the analyses. In addition, six stage-I salient genes, namely CPXM2, NPR3, PALM, PRDM6, TAGLN, and TPM2 are identified as stage-IV salient here. However the concordance between the markers from the reference TCGA dataset and GSE39582 is not extensive, and merits discussion. Foremost, GSE39582 is limited to colon cancer samples, which might differ in some features from rectal cancers, thereby missing some variation that is captured in the TCGA COADREAD dataset. Second, we would like to note that out-of-domain cohorts might be sensitive to distribution shifts in gene expression, which require measurement calibration with an adequate number of normals from the same (new) cohort. Since there were few normal samples in the original GSE39582 dataset, this might significantly skew the extension of gene signatures established with the reference TCGA cohort. The addition of 308 normal colonic samples available in the GTEx does not mitigate this issue, since (i) these are from an entirely different cohort, and (ii) normal rectal tissue samples remain unaccounted for. In addition, the applicability of candidate biomarker signatures to new cohorts might be bounded by bioinformatic problems pertaining to data curation and processing. The contrarian findings prompted us to seek robust validation of the models developed below.

Development of a diagnostic aid for colorectal cancer screening

We combined the 157 stage-salient genes, top ten genes from linear modeling, and the 18 genes that were both linear and monotonically expressed into a single expression feature-space of 185 genes. The TCGA dataset was randomly split into a train dataset of 287 cancer and 41 normal samples, and a holdout testset of 71 cancer and 10 normal samples. Application of the feature selection techniques yielded a consensus feature space of just seven essential features, viz. four of the top ten linear modelling genes (ESM1, DHRS7C, OTOP3, AADACL2), two stage-salient genes (stage-2 salient LPHN3 and stage-4 salient GABRD) and one linearly monotonic gene (LPAR1). Using these features, four different ML models were trained, and hyperparameters optimized. The models were ranked on their performance on the training and holdout test sets (Table 9), and the Random Forest and 2-layer neural network models were identified for blind external validation.

Two external datasets were chosen for blind validation: (i) Rectal_cancer_MSK (Chatila et al., 2022) with 113 mRNA-Seq expression samples, obtained from 738 primary rectal tumors (https://www.cbioportal.org/; accessed 02-09-2023; accessed 02-09-2023) and (ii) 308 normal colon samples from the GTEx (accessed 12-10-2022). It is noted that the microarray-based GEO datasets benchmarked in our study, namely GSE25071, GSE21510, and GSE39582 were limited in the coverage of the gene-space, lacking expression values for some of the seven features used in the models, and not further considered. The hyperparameter-optimized Random Forest and 2-layer neural network models were re-built on the full TCGA dataset and evaluated on the external datasets (Table 10). All the cancer samples were correctly predicted by the Random Forest model, yielding ‘perfect’ recall. There were just eleven misclassified instances out of the 421 samples in the combined external dataset, and all such instances were normal colon tissue samples, leading to a balanced accuracy of 98.27%. The Random Forest model outperformed the 2-layer Neural network model on all the metrics considered, including sensitivity, specificity, F1-score, and Matthews correlation coefficient (MCC).

Development of a prognostic model for colorectal cancer

All the 157 stage-salient genes were subjected to univariate Cox regression analysis, and the significant biomarkers (P < 0.05) are presented in File S14. Of the top stage-salient genes, five emerged significant, namely JPH3, HOTAIR, CST6, GABRD, and DKK1 (all P < 0.03). HOTAIR, CST6, GABRD, and DKK1 are stage-IV salient, while JPH3 is stage-I salient (Fig. 8). Multivariate Cox regression analysis with feature selection yielded an optimal panel of three genes, namely HOTAIR, GABRD, and DKK1, with a model p-value ∼5e−04, and individual significances ∼0.0086, 0.0053, and 0.0238, respectively (i.e., all p-values < 0.05). The multivariate risk model was given by:

Risk-score = 0.14872 * HOTAIR + 0.4423 * GABRD + 0.10877 * DKK1

The hazard rate for all the prognostic factors significantly exceeded 1.0, indicating that the constituents of the biomarker panel elevated the prognostic risk, suggesting possible oncogenic roles in line with their overexpression. The distribution of risk scores yielded a median maxstat value of 2.74 for patient risk stratification. Further, the Kaplan–Meier curve of the multivariate model suggested that the high-risk group was significantly associated (p-value <0.0014) with a poorer overall survival than the low-risk group (Fig. 12D). The model yielded an acceptable Concordance index (C-index) ∼0.71 ± 0.05, suggesting further application as a novel prognostic panel (Svoboda et al., 2014; Niu et al., 2020; Sui et al., 2019). It is significant (and perhaps not surprising) that the identified prognostic panel is entirely composed of stage-IV salient biomarkers, suggesting that the distance to metastasis is the single dominant factor in the stratification and determination of prognosis of colorectal cancer.

Stage-1 salient DEGs

The genes CALB2 and TMEM59L cluster together in Fig. 7 with the least stage-I expression, suggesting the hypothesis that they function as tumor suppressor genes. This is supported in the literature, specifically that cells in which CALB2 is silenced do not respond to 5-flourouracil, a favored treatment for CRC, indicating that CALB2 expression is necessary for 5-flourouracil induced apoptosis (Stevenson et al., 2011). Another study found that heterozygosity in SNP513 of Intron 9 of the gene CALB2 might be a predictive marker for CRC (Vonlanthen et al., 2007). It has also been noted that increased TMEM59L expression was a pro-apoptotic indicator of cell death during oxidative stress in neuronal cells (Zheng et al., 2017). Regarding SOX2 and SOX10, it is noteworthy that the Cancer Genome Atlas Network observed SOX9 as a novel gene with significant recurrent mutations in COADREAD (The Cancer Genome Atlas Network, 2012).

Stage-2 salient DEGs

KLHL34 was found to be hypermethylated in Locally Advanced Rectal Cancer, and knockdown of KLHL34 lowered colony formation, increased cytotoxicity, and increased radiation induced caspase 3 activity in LoVo cells (Ha et al., 2015). CCBP2, encoding the Chemokine decoy receptor D6, has an inhibitory effect on breast cancer malignancies due to its action to sequester pro-malignant chemokines (Yang et al., 2013). The lncRNA PLAC2 induces cell cycle arrest in glioma by binding to Ribosomal Protein RP L36 in a mechanism involving STAT1 (Hu et al., 2018). GPC5 was found to be overexpressed in the lung cancer phenotype (Li & Yang, 2011), in lymphoma, and in gastric cancer.

Stage-3 salient DEGs

Copy number polymorphisms of TRY6 gene have been found in breast cancer (Wagner et al., 2007). HABP2 gene overexpression has been observed in lung adenocarcinoma and has been proposed as a novel biomarker for the same (Wang et al., 2002).

Stage-4 salient DEGs

The lncRNA HOTAIR was found to be significantly overexpressed in HCC, and a potential biomarker for lymph node metastasis in HCC (Geng et al., 2011), and later implicated in different cancers (Hajjari & Salavaty, 2015). Another widely-cited study (Gupta et al., 2010) showed that enforced HOTAIR gene expression in epithelial cancer cells induces chromatin reprogramming and an increased metastatic state, while inhibition of HOTAIR inhibits cancer invasiveness. These accounts of the role of HOTAIR in metastasis accord with our findings that HOTAIR is a stage-4 salient significantly monotonically expressed biomarker. GWAS analysis identified a strong association of C6orf15 with occurrence of follicular lymphoma (Skibola et al., 2009). Promoter methylation of cell free DNA of the CST6 gene was found to be a potential plasma biomarker for Breast Cancer (Chimonidou et al., 2013). Expression of VGLL1 and its intronic miRNA miR-934 are associated with sporadic and BRCA1-associated triple negative basal-like breast carcinomas (Castilla et al., 2014). Expression of DKK1, an inhibitor of osteoblast differentiation, was found to be associated with the presence of bone lesions in patients with multiple myeloma (Tian et al., 2003). TMEM40 has been found to be a potential biomarker in patients with Bladder cancer, serving as an oncogene and a possible therapeutic target (Zhang et al., 2018). The emergence of the C,E, and F members of the Lymphocyte Antigen 6 (LY6) family (Loughner et al., 2016; Upadhyay, 2019) as monotonically expressed proto-oncogenes holds promise for immunotherapy. There is a substantial evidence base for GABRD (Gross, Kreisberg & Ideker, 2015), which is a key component of both the screening and prognostic models developed here. Consistent expression trends in GABRD and other stage-salient MEG DEGs provide unmistakable evidence for the existence of molecular signatures in CRC progression.

Benchmarking with curated databases

We found 13 of the top 200 genes from the linear model documented in the CGC v84 as known cancer genes (Table 11). Two genes, MACC1 and SALL4, were specifically documented for colorectal cancer. HSP90AB1 had been earlier identified as a top MEG in HCC (Sarathi & Palaniappan, 2019). Screening the 157 stage-salient genes against the NCG7.0, which is a curated database of cancer drivers and healthy drivers, yielded 28 genes, of which eight were in the top 40 stage-salient genes (File S17). All the hits were documented to carry mutations in their coding region (vs noncoding region). Three were canonical oncogene drivers, namely HOXC11, SOX2, and KCNJ5, while the rest 25 are putative oncogenes and putative tumor suppressors in almost equal measure. Two stage-salient genes, namely CNTN1 and BAI3 (ADGRB3) were documented as putative tumor suppressor genes involved in gastric adenocarcinoma, providing specific support for our findings. PIGR is identified as an essential healthy driver (Olafsson et al., 2020), signifying that mutations in this gene confer an exceptional protective effect, and its down-regulation could drive tumorigenic processes. Intriguingly, the stage-salient genes C5ORF23 (NPR3), SOX2, and KCNJ5 are the only instances where the documentation is dissonant with our primary findings; these three were marked as putative oncogenes, though they are identified as down-regulated here. Further investigations in this direction are warranted to set the literature straight. Documentary evidence for drugs targeting any of these genes is absent, emphasizing the value of the present study in pinpointing novel candidates for diagnosis, therapy and prognosis. To perform a systematic analysis of therapeutic interventions based on these targets, we consulted ClinicalTrials.gov for clinical trials targeting stage-salient genes. Ten genes from the top stage-salient genes are being pursued in clinical trials, either as target or endpoint, colorectal or other cancers. Details of clinical trials along with the current status/phase of each trial are provided in File S18. DKK1 and HOTAIR are the only stage-4 salient genes implicated as targets/endpoints in clinical trials. DKK1 is involved in three clinical trials for colorectal and gastric cancers. HOTAIR is the target of a clinical trial for thyroid cancer (ctgov:NCT03469544) (Abudoureyimu et al., 2016). HOTAIR is documented in the NONCODE database (http://www.noncode.org/) as disease-associated, specifically with colorectal cancer (ID: NONHSAG011264.3), validating its role in oncogenic processes. It is notable that GABRD is not a target in any of the registered clinical trials, flagging a prime potential interest for future efforts. LPHN3, a stage-2 salient gene, is targeted in four clinical trials aimed against metastatic colorectal cancer, to explore possible therapeutic efficacy in thwarting cancer progression prior to irreversible outcomes. FADS6 (a stage-II salient gene) is an endpoint in a clinical trial to treat colorectal adenomatous polyps, which is a precursor to malignant lesions. CALB2 and C5orf23 (NPR3) are each involved in one clinical trial related to colorectal cancer. Some stage-salient genes are being pursued in treatment of cancers in other cell types/tissues, underlining the role played by certain genes in contributing to general cancer hallmark processes (Hanahan & Weinberg, 2011). Specifically NAT2 is a target in nine different clinical trials against diverse cancers, significantly highlighting its essential role in driving hallmark processes in unrelated cancers.

Insights from network analysis

Stage-wise network analysis of colorectal cancer progression has shed light on certain genes potentially underlying progression (Rahiminejad et al., 2022). The strength of the computational evidence for the candidate biomarkers identified herein urged a network analysis to examine the findings in a larger context. The intersection between the sets of all stage-salient biomarkers and the significant MEGs might highlight monotonically enriched pathways essential to the pathophysiology of colorectal cancers. Hence the 31 stage-salient MEGs were chosen to reconstruct the STRING network, with 50 interactors in the first shell and 10 interactors in the second shell. This yielded a PPI with 235 edges with an extremely significant enrichment p-value < 1.0e−16 (Fig. 10). A Gene Ontology (Ashburner et al., 2000) analysis of this reconstructed network showed enrichment for the Wnt-Frizzled-LRP5/6 complex component at p-value < 1E-04. An analysis with KEGG (Kanehisa et al., 2016) showed enrichment for 2-oxocarboxylic acid metabolism at p-value ∼0.001, indicating a Warburg-shift in metabolism. An analysis with Reactome (Sidiropoulos et al., 2017) showed significant enrichment of SMAD2/3 and SMAD4 MH2 Domain Mutants in Cancer (p-value < 0.01). These observations in toto provide striking evidence for the involvement of these biomarkers in driving CRC progression.

Isolated nodes in the network included GABRD, DLX3, ISM2, LY6G6C & E, SPYDC, UPK2, C2orf48, PIGR, KRTAP3-1, C7orf52 (NAT16), SPERT, and PLAC1. All the isolated nodes are proto-oncogenic (see Table 7), hence could provide targets for inhibition in personalized cancer medicine. An outlier component (not in the giant connected component) was made of the CXCL chemokine family, stemming from CXCL13—a recently recognized immune checkpoint with a key role in tumor progression (Yang et al., 2021; Ren et al., 2022). This component could constitute a novel target for upregulation in CRC immunotherapy. A drug-repurposing search with the DrugGeneBadger (Wang et al., 2019) for each of the 31 stage-salient MEGs yielded drugs (small molecules with q-values < 0.05) to pharmacologically alter the expression of these identified biomarkers. The search revealed that curcumin is effective against at least 13 of these targets, and piperlongumine against eight of these targets. Six biomarkers (HOTAIR, ISM2, KRTAP3-1, SPDYC, LY6G6F, and NKD1) found no drug available in LINCS1000 (Subramanian et al., 2017) to modulate their expression, and these constitute potential novel targets for drug discovery against metastatic transition in CRC.

A network specific to colon cancer could be obtained using the results for GSE39582, and is presented in File S19. Enrichment analysis of this network with Gene Ontology indicated significance for Arp2/3 complex-mediated actin nucleation (p-value ∼1e−4), which is known to contribute to invasive colorectal cancer (Zheng et al., 2023). A KEGG analysis showed enrichment for oxidative phosphorylation (p-value ∼1e−20), with a prominent clustering of NDUF and COX gene families. A Reactome analysis showed a minor enrichment of enzymatic protein conjugation processes (UBE2I, UBA2, SAE1) that monitor intracellular proteins and cell states (p-value ∼0.02). These findings indicate an enrichment of proliferation-independent metabolism-rewiring pathways necessary for colorectal cancer progression, and could be contrasted with the analyses in Marisa et al. (2013).

Immune-cell infiltration analysis

Deconvolution of the TCGA samples based on the LM22 immno-cyte signature with 100 permutations yielded 107 samples with significance (P < 0.05), including eleven controls. These samples, with their TCGA identifiers, are presented in File S20. The significant samples were analyzed for the relative abundance of the 22 immune cell types. A heatmap of the sample-wise immune cell-type proportions was generated (File S21A), and the clustering patterns of the cell-types across samples was visualized using a dendrogram. We observed the following clusters: mast cells resting and plasma cells; mast cells activated and neutrophils; T cells CD8, T cells follicular helper, and macrophages M1; T cells CD4 memory resting and B cells naive. The macrophages M0 and M2 were clear outgroups in the dendrogram. A normalized stacked bar chart of the sample-wise immnuo-cyte fractions revealed substantial variations in immune cell-type composition between normal and cancer samples (File S21B). To investigate further, we analyzed the differences in distribution of cell proportions between normal and tumor samples for each immune cell-type (File S21C; data presented in File S20). Eight of the 22 immunocyte types showed significant distribution differences (adj. P < 0.05). Specifically, we found the infiltration of four immune cell-types preferentially enriched in tumor samples, namely macrophages M0, T cells CD4 memory activated, mast cells activated, and neutrophils, while four other immune cell-types were preferentially depleted in tumor samples, namely macrophages M2, T cells CD4 memory resting, mast cells resting, and plasma cells. In particular, macrophages M0 exhibited both the largest effect size (>2.0) and the greatest significance (<1E-07) of infiltration in tumor samples. The preferential enrichment of mast cells activated and T cells CD4 memory activated versus the preferential depletion of mast cells resting and T cells CD4 resting suggested that tumorigenesis activates resting immune cell-types, potentiating their infiltration of the tumor microenvironment. To integrate these observations, we computed the correlation matrix of the immune cell-types based on their sample-wise proportions over both normal and tumor samples (File S21D). The largest positive correlations were exhibited by T cells follicular helper with T cells CD8 (Pearson’s ρ ∼0.52), and with macrophages M1 (Pearson’s ρ ∼0.45), reinforcing their clustering in the dendrogram. Intriguingly, the largest negative correlation (in magnitude) was exhibited by macrophages M0 and T cells CD4 memory resting (Pearson’s ρ ∼−0.51) (File S21D). Given that macrophages M0 are preferentially enriched in tumor samples whereas T cells CD4 resting and mast cells resting (Pearson’s ρ ∼−0.47 with macrophages M0) are both preferentially depleted, these observations cohere and could hold preliminary significance for immunotherapy. Discovery of multicellular immunocyte community structures could pave the way for personalized immunotherapy in CRC treatment (Sarathi & Palaniappan, 2019; Ge et al., 2019).

COADREADx

Based on the external validation, the Random Forest model was identified as the best model for screening early-stage cancer. Coupled with the prognostic model, these could aid the risk stratification of patient samples. With this application in mind, we have deployed COADREADx, an experimental web service for the screening of patient samples as ‘cancer’ or ‘normal’, and subsequent prognostication in the case of ‘cancer’. COADREADx has been implemented using R-Shiny (https://shiny.rstudio.com/), and is available for academic use at: https://apalanialab.shinyapps.io/coadreadx/. A help document with sample input files for different use-cases, and a companion how-to video have been made available on the landing page. To aid the effective interpretation of COADREADx predictions, the prediction probability for the predicted diagnostic class is provided, yielding a level of confidence in the prediction. Similarly the risk stratification of ‘cancer’ samples is accompanied by the quantile of the estimated risk-score as well as its fold-change from the median value of the risk score distribution. These values suggest the strength of evidence for the predicted risk class.

In summary, we have performed a novel de novo analysis of the TCGA COADREAD gene expression dataset, and identified multiple interesting classes of biomarkers. The biomarkers have been validated with alternative datasets, network analysis and immune cell infiltration analysis. Some of the biomarkers could suggest novel hypotheses for targeted therapy and immunotherapy. Using purifying techniques, we have carved feature spaces from these biomarkers to build screening and prognostic models of colorectal cancer. The screening model has been externally validated, while the prognostic model has been bootstrapped for confidence. Both the models have been deployed via COADREADx, a web-server designed to return confidence estimates for all its predictions. Phenomena of distribution drift and shift in new samples and out-of-domain cohorts challenge the applicability of COADREADx, which might need refinement in the light of such evidence. Further validation of the models on colon cancer samples might also be warranted. Enabling risk stratification is vital to treatment strategy and clinical management of the cancer. Thus experimental validation and further improvement of COADREADx is necessary to demonstrate its clinical utility for screening and prognosis purposes. It is reckoned that the availability of such software-as-medical-devices could ease the accessibility to effective surveillance technologies for early detection of colorectal cancer (Muthamilselvan, Ramasami Sundhar Baabu & Palaniappan, 2023).