Mining and exploration of rehabilitation nursing targets for colorectal cancer

Download the colorectal cancer data sets

Colorectal cancer data sets GSE32323 and GSE113513 were extracted respectively from GEO database, specifically generated from GPL6102 and GPL14951. GSE32323 consisted of 17 bowel cancer cases and 17 normal cases, while GSE113513 consisted of 14 bowel cancer cases and 14 normal cases. The primary objective was to investigate differentially expressed genes (DEGs) in relation to colorectal cancer.

To batch processing

An orchestrated approach was used for the integration and disentangling of the multiple data sets. Initially, the R software package facilitated the amalgamation of the GSE32323 and GSE113513 datasets. Using the R package In SilicoMerging, a consolidated matrix was obtained by merging these data sets. Batch effects were then removed with the help of the remove patch effect routine in limma R (version 3.42.2). The resulting batch-independent matrix was used as the basis for subsequent analysis.

Screening of DEGs

The fusion matrices of GSE32323 and GSE113513 were subjected to a fusion clustering and alignment using the R package limma. The Benjamini-Hochberg approach was employed to refit p-values to increase statistical robustness. False discovery rate (FDR) was applied to calculate fold change (FC) using a stringent set of criteria with a threshold of P < 0.05 and FC >1.5. The identification of differentially expressed genes (DEGs) was facilitated by generating a visually informative volcano map.

Weighted gene co-expression network analysis (WGCNA)

First, a batch correction process was performed to use the genome description matrix by merging GSE32323 and GSE113513. Then, excluding the upper 50% of the smallest genes, the Median Absolute Deviation (MAD) was calculated for each gene. Outlying genes and outlying samples were identified and excluded using the goodSamplesGenes method from the R package WGCNA. A coexpression scaled network was then assembled using WGCNA. For further refinement of the module analysis, the genetic differences within the modules have been calculated and cut lines have been selected on the module tree. In particular, modules with a merging distance of less than 0.25 were consolidated. Importantly, the grey module was considered unassignable to a gene module collection.

Functional enrichment analysis

The goal of the study was to identify functional differences within the gene lists through extensive computational analysis using GO and KEGG methods. By submitting the gene list to the following KEGG Web API (https://www.kegg.jp/kegg/rest/keggapi.html), the current KEGG Annotation for this gene was retrieved. Using this information as a background, genes were assigned and subjected to clustering analysis using the R package clusterProfiler (version 3.14.3). Gene set enrichment results were considered statistically significant if they met the criteria of a lowest gene set of 5, a highest gene set of 5000, and a P-value < 0.05 with a FDR <0.25.

In addition to this, the Metascape database, which is a valuable resource for the comprehensive annotation and analysis of genes, was also used. The gene enrichment analysis results were further explored and visualised using the Metascape database (http://metascape.org/gp/index.html). This process provided a thorough understanding of the functional implications of the identified gene differences and allowed the export of relevant gene lists for more detailed study.

GSEA

Gene Set Enrichment Analysis (GSEA) was purchased from the official GSEA homepage. Two different groups of colorectal cancer and normal tissue samples were stratified. Subsets were then downloaded from the Molecular Signatures Database, which allowed the evaluation of the corresponding pathways and underlying mechanistic pathways. Using the genetic expression profile and group of phenotype, GSEA was performed with a lowest gene set size of 5 and a highest gene set size of 5000, using 1000 resampling iterations. Results were considered statistically significant if they met the criteria of a P-value of < 0.05 and FDR of < 0.25. In addition, GSEA was extended to include Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) annotations to provide a comprehensive exploration of the colorectal cancer genetic landscape.

Immunoinfiltration analysis

The computational tool CIBERSORT (http://CIBERSORT.stanford.edu/) was used as a commonly accepted technique to assess infiltration of immune-system cells in colorectal cancer. An integrated bioinformatics approach using the CIBERSORT software package was adopted to comprehensively analyse the combined gene expression matrix of the colorectal cancer datasets GSE32323 and GSE113513. Using the principles of linear regression, this method facilitated the unfolding of the immune cell component expressivity matrices and provided estimates of immune cell abundance. A confidence level (P < 0.05) was used to ensure that data were included with sufficient reliability and accuracy.

Construction and analysis of protein-protein interaction (PPI) networks

The Search Tool for the Retrieval of Interacting Genes (STRING, http://string-db.org/) database is a comprehensive platform for aggregating, evaluating and integrating publicly available data on protein-protein interactions, complemented with computational predictions. In this investigation, the index of differentially expressed genes was entered into the STRING database to establish a PPI network to predict key genes (confidence >0.4). Cytoscape software, renowned for its capabilities in biological network analysis and two-dimensional (2D) visualisation, was used to visualise and predict these key genes within the PPI networks.

The process involved importing the PPI mesh into Cytoscape and identifying the clusters with the strongest correlations using the MCODE algorithm. Five different tools (MCC, MNC, EPC, Eccentricity, Betweenness) were applied to predict the most robustly correlated genes and their intersections. A list of core genes contributing to a broader understanding of the underlying molecular interactions was identified through visualisation.

Survival analysis

When examining clinical survival data for colorectal cancer from The Cancer Genome Atlas (TCGA), we employed the R software package maxstat (version: 0.7–25) to determine the ideal RiskScore truncation value associated with core genes. By defining a lowest sample size of 25% and a highest sample size of 75%, we derived the ideal truncation values. Patients were then stratified into a hyper and hypo group based on these values. The surpass function of the R software package survival was then employed to examine the prognographic differences between these groups, and the difference was assessed for significance using the test for significance (log-rank test).

In order to assess the impact of each of the independent core genes on the prognosis of colorectal cancer, we observed any notable changes in effect. Receiver operating characteristic (ROC) analyses were also performed by means of the R package pROC (version 1.17.0.1). Patient follow-up at specific time points (365, 1095 and 1825 days) was included in the ROC analyses. The pROC ROC function and ci function assessed area under the curve (AUC) alongside confidence intervals, providing a comprehensive assessment of core gene prognostic score diagnostic value.

Heat map of gene expression

Visual representations in the form of heat maps illustrating the level of regulation of the core genes clustered within the PPI network were generated using the R package heatmap. This approach allowed us to visually explore the nuanced differences in expressing colorectal cancer and regular tissue samples within the GSE32323 and GSE113513 datasets.

CTD analysis

The Comparative Toxicogenomics Database (CTD) is an invaluable resource that seamlessly integrates a wide range of data, including chemical, genetic and functional information. This comprehensive database facilitates the exploration of interactions between phenotypes and disease-related environmental exposures, providing considerable convenience for investigating potential disease and drug mechanisms.

For further analysis, the identified ‘key genes’ were entered into the CTD website and the primary diseases associated with these ‘key genes’ were identified. The complex relationships between each gene and their respective diseases were further clarified by generating a differential expression radar map for each gene using Excel. This process has added a layer of visual clarity to the complex interplay between the core genes and their associated diseases and has improved our understanding of the potential mechanisms and relationships involved.

The miRNA

In our investigation, we utilised the capabilities of TargetScan (http://www.targetscan.org). This online database is highly regarded for its ability to predict and analyse miRNAs and their target genes. Our use of TargetScan focused on the meticulous screening of miRNAs that are intimately involved in the regulation of key differentially expressed genes (DEGs). The goal was to reveal the complex regulatory relationships between miRNAs and the key genes we identified.

Data availability

The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

Exploration of genes with altered expression

In this investigation, we used pre-defined cut-off criteria to identify genes with altered expression levels within the merged dataset of the colorectal cancer datasets GSE32323 and GSE113513. As shown in Figure 1, the result of this analysis was a comprehensive set of 1117 DEGs.

Figure 1. Exploration of differentially expressed genes (DEGs). A comprehensive set of 1117 DEGs was identified and analyzed, shedding light on the molecular distinctions underlying colorectal cancer.

Functional enrichment analysis

DEGs

As elegantly illustrated in Figure 2A–2C, based on Gene Ontology (GO) analysis, the enriched terms predominantly revolved around key biological functions including cell cycle, cell reproduction and nuclear content containing protein complexes.

Figure 2. Analysis of the functional enrichments. (A–C) GO. (D) Analysis of the KEGG. (E–H) GSEA.

The KEGG analysis highlighted significant pathways, with particular emphasis on calcium signalling, P53 signalling, colorectal cancer and the intestinal immune network of IgA, as shown in Figure 2D.

Metascape enrichment analysis

In Figure 3A, the GO enhancement project unfolds graphically in the Metascape enrichment exploration. Simultaneously, we generated an output visualising the enrichment network, where the colouring reflects the enrichment terms and their corresponding p-values (Figures 3B, 3C and 4). This visual representation artfully conveys the interconnectedness and confidence of each enrichment term, enhancing clarity and depth of understanding.

Figure 3. Enrichment analysis in the metascape. (A) In the metascape enrichment project, the GO enrichment project is visible. (B) The enrichment network coloured by enrichment terms. (C) The enrichment network coloured by p-values.

Figure 4. Analysis of the metascape enrichment. (A) The different modules of MCODE. (B) PPI network in the MCODE.

Immunoinfiltration analysis

Insights into immune cell composition were obtained by comprehensive analysis of the combined matrix from the colorectal cancer datasets GSE32323 and GSE113513 by means of the CIBERSORT software package. The proportions of immunological cells within the total genomic pattern were determined at a 95% confidence level. In particular, M2 and activated mast cells were found to be dominant components in these samples (Figure 5A).

Figure 5. Analysis of immune infiltration. (A) Activated M2 macrophages and mast cells comprised a high proportion of the sample. (B) A Heatmap of immune cell expression in the dataset. (C) The correlative analysis was conducted on immune cell infiltration, the common expression patterns between components of immune cells, macrophages M2 high expression, activated mast cells can express is low.

In a hematogram depicting the expression levels of immune cells in the dataset, a particularly increased representation of activated mast cells was observed in colorectal cancer samples (Figure 5B). To unravel potential relationships between immune cell infiltrations, correlation analysis was performed, which revealed common expression patterns. Notably, an inverse correlation was observed between high M2 macrophage expression and low activated mast cell expression. As shown in Figure 5C, this complex interplay may have implications for colorectal cancer progression.

WGCNA

A critical step in WGCNA is the soft threshold power selection. Using network structure analysis, we set the performance soft-threshold to 9, which is the bottom performance that achieves a scale-freeness topological fit index of 0.9 (Figure 6A). A clustered hierarchical tree was then established, revealing 11 distinct assemblies of genes (Figure 6B). In summary, in addition to generating a correlation heatmap between modules and phenotypes (Figure 7A), in-depth analysis focused on interactions within an essential module (highlighted in Figure 6C). A scatterplot (Figure 7B) was also generated to explore the correlation between gene significance (GS) and module membership (MM).

Figure 6. WGCNA. (A) β = 16, 0.87. β = 16, 5.23. (B) All the genes involved in the production of the 11 modules. (C) The interaction between the most important module.

Figure 7. WGCNA. (A) A module and a phenotypic correlation heatmap. (B) The corresponding hub gene of GS and MM correlation scatter plot. (C) Venn diagram of WGCNA and the differentially expressed genes screened by DEGs, and their intersections.

We performed correlations between module signature vectors and gene abundance, emphasising a strict cutting standard (|MM| >0.8), to identify critical hub genes. Ten strongly associated genes in clinically meaningful modules were clustered as hub genes. To further elucidate the molecular landscape, a Venn diagram was constructed to illustrate the overlap between genes identified by WGCNA and those identified as DEGs (Figure 7C).

Protein - protein interaction (PPI) network building and analysis

A complex PPI network was established by integrating WGCNA with gene lists from DEGs. This network was meticulously generated with the STRING online database and analysed using Cytoscape software (Figure 8A), yielded a recognisable core set of genes as highlighted in Figure 8B.

Figure 8. Protein - protein interaction (PPI) network building and analysis. (A) PPI net. (B) The nuclear cluster. (C–G) Five algorithms (MCC, MNC, EPC, EcCentricity, Betweenness) were used for the identification of the core genes. (H) To obtain the combination as core genes, the Venn diagram was made.

In an effort to identify central genes within this network, we applied five different algorithms (MCC, MNC, EPC, EcCentricity, Betweenness), which resulted in the identification of potential central genes, each of which is shown in Figure 8C–8G. An elaborate Venn diagram analysis merging the results of these algorithms accurately pinpointed the core genes, as depicted in Figure 8H. At the end of this analysis, six genes (MYC, MAD2L1, CENPF, UBE2C, NUF2 and NCAPG2) emerged as the focal points of interest within the complex molecular landscape.

Survival analysis

We obtained colorectal cancer survival data from TCGA and explored the complex association of colorectal cancer risk scores with gene expression patterns. The research produced a comprehensive map showing the relationship between prognostic scores and differences in the methylation of key genes between colorectal cancer and regular test tissue samples. Notably, there was a significant decrease in patient survival as the risk score increased, with the low-risk group having significantly longer periods and rates of survival than the high-risk group (Figure 9A).

Figure 9. Survival analysis. (A) Prognosis score relationship map and expression difference heat map of core genes between colorectal cancer and normal tissue specimen. (B) The core genes were protective factors, showing decreasing expression trends with increasing risk score. (C) Risk score ROC curve.

In addition, the identified core genes functioned as protective factors, showing a consistent downregulation as risk score increased (Figure 9B). To estimate the predictive ability of the risk scores, we constructed an ROC curve, which showed promising results with AUC values indicating the forecasting potential of the risk scores (Figure 9C).

Distinct expression patterns were observed in the box plot analysis of the core genes in colorectal cancer. Notably, the core genes (MYC, MAD2L1, CENPF, UBE2C, NUF2, NCAPG2) had differential expression in colorectal cancer samples in comparison to normal samples. As elegantly illustrated in Figure 10, these core genes showed an increased expression in colorectal cancer samples, contrasting with a decreased expression in normal samples.

Figure 10. Survival analysis. Box results of colorectal cancer core genes.

Heat map of core gene expression

By visualising the expression levels in the combined matrix of GSE32323 and GSE113513, a striking pattern was revealed. The core genes, namely MYC, MAD2L1, CENPF, UBE2C, NUF2 and NCAPG2, showed elevated expression in colorectal cancer samples, in stark contrast to their comparatively lower expression in normal samples (shown in Figure 11A). This intriguing observation leads to speculation about the potential regulatory role of these master genes in the context of colorectal cancer.

Figure 11. (A) Heat map of the expression of core genes. (B) Analysis of the CTDs.

CTD analysis

As part of this research, we submitted the cluster gene list to the CTD website, which provided valuable insights into the diseases closely associated with the cluster genes. Specifically, the core genes (MYC, MAD2L1, CENPF, UBE2C, NUF2, NCAPG2) showed associations to colorectal neoplasia, tumour cell transformation, inflammation and necrosis, as elegantly illustrated in Figure 11B. This nuanced exploration has deepened our understanding of the intricate interplay between these core genes and a variety of disease processes.

Hub gene associated miRNA prediction and functional annotation

Here, the index of target genes was incorporated into TargetScan, which revealed relevant miRNAs to improve the understanding of gene expression regulation (Table 1). Interestingly, MAD2L1 was found to be associated with miRNAs such as hsa-miR-6088, hsa-miR-143-3p and hsa-miR-4770. CENPF genes showed associations with miRNAs such as hsa-miR-302c-3p and hsa-miR-520f-3p. In addition, the NUF2 gene was associated with hsa-miR-599. The intricate regulatory networks involving these genes and their associated miRNAs are revealed by this nuanced investigation.