IL-8, MSPa, MIF, FGF-9, ANG-2 and AgRP collection were identified for the diagnosis of colorectal cancer based on the support vector machine model

ABSTRACT

Colorectal cancer (CRC) is one of the most common cancer, and the early detection of CRC is essential to improve the survival rate of patients. To identify diagnostic markers for colorectal cancer (CRC) by screening differentially expressed proteins (DEPs) in CRC. The DEPs were initially obtained from 12 CRC samples and 12 healthy control samples, and verification analysis was performed in another 34 CRC samples and 34 normal controls. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment with DEPs was analyzed by the R package clusterProfiler (Version 3.2.11), and the DEP-associated protein-protein interaction (PPI) network was created from the STRING database. Additionally, Support Vector Machine (SVM) model prediction and survival analyses were conducted on the key DEPs. Preliminary screening and functional analysis showed that the DEPs mainly overrepresented in pathways such as cytokine-cytokine receptor interaction, chemokine signaling pathway, Rap1, Ras, and MAPK signaling pathways. The key DEPs, including AgRP, ANG-2, Dtk, EOT3, FGF-4, FGF-9, HCC-4, IL-16, IL-8, MIF, MSPa, TECK, TPO, TRAIL R3, and VEGF-D, were used to construct a custom chip. The drug-gene interaction network suggested that TPO was a key drug target. ROC curve showed the SVM diagnostic model with the DEPs IL-8, MSPa, MIF, FGF-9, ANG-2, and AgRP had better diagnostic performance with an AUC of 0.933. Survival analysis showed the expression of FGF9, TPO, TRAIL R3, Dtk, TECK and FGF4 were associated with prognosis. This study revealed the important serum proteins in the pathogenesis of CRC, which might serve as useful and noninvasive predictors for the diagnosis of CRC.

KEYWORDS:

• Colorectal cancer
• differentially expressed proteins
• diagnosis
• svm model prediction

Introduction

Colorectal cancer (CRC), one of the common gastrointestinal cancers, is the fourth most frequent cause of cancer mortality with more than 1.2 million annual new cases globally [1–3]. Although there have been great advancements in surgical and adjuvant therapies, the prognosis of CRC still remains poor [4,5]. And due to the lack of symptoms in the early stages of CRC, more than half of CRC cases are diagnosed at an advanced stage or present with distant metastasis [6,7]. Furthermore, because CRC develops slowly from removable precancerous lesions, early detection of CRC is essential to improve the survival rate of patients. Hence, there is an urgent need for credible and noninvasive biomarkers for the early detection of CRC.

High throughput techniques have played critical roles in exploring effective and nonaggressive biomarkers for cancer detection [8]. Proteomic analysis, which is efficient, powerful, and cost-effective, is especially commonly used in clinical research [9]. Currently, several tumor biomarkers, such as CA125, carcinoembryonic antigen (CEA) and CA199 are regarded as the conventional predictors for CRC [10,11]. However, their reliability remains controversial because they are insufficiently sensitivity. In this study, we screened key differentially expressed proteins (DEPs) between CRC serum samples and normal controls based on Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway and Gene ontology (GO) functional analyses, analysis of the Protein-Protein Interaction (PPI) network, and prediction of drug-gene interactions. Subsequently, 15 key DEPs were chosen to construct a custom chip, and the candidate DEPs associated with the diagnosis and prognosis of CRC patients were determined using Support Vector Machine (SVM) model prediction and survival analysis.

Materials and methods

Ethical statement

The study was approved by the Ethics Committee of China-Japan Union Hospital. The informed consents were obtained from all participants.

Samples collection

In total, 24 serum samples, including 12 CRC samples (age: 61.4 ± 12.0; 6 men and 6 women) and 12 healthy control samples (age: 60.9 ± 12.8; 6 men and 6 women) were included in this study, and none of the patients had received surgery or chemotherapy. All patients were pathologically diagnosed with CRC. All samples were processed immediately, and the processing steps were as follows: after extracting the blood, it was immediately put into the centrifuge tube and placed vertically on the rack of the centrifuge tube. Then the blood was placed in the centrifuge for 1 hour and centrifuged at 3000P/m for 10 minutes. The protein was obtained directly from serum as follows: the serum samples were centrifuged at 10,000 RPM and 4 °C for 5 min, and the supernatant was taken to obtain serum protein. 120 proteins were collected, and then divided into two chips (C6 and C7, which were two corresponding membranes used for covering 120 proteins), with 60 proteins per chip. The proteins were then sorted by the indexes of cumulative strength, average cumulative strength, and average net cumulative strength through detecting by the chemiluminescence imaging system (Tanon, Shanghai, China).

Data preprocessing

The raw datasets with the three indexes were preprocessed to adjust for the background using the limma package (Version 3.26.9; http://bioconductor.org/packages/release/bioc/html/limma.html) [12]. After the data was normalized with the quantile method, the log (X + 1, 2) transformation (X presents the expressed value) was applied to the data.

Identification of DEPs and hierarchical clustering

The limma package was used to screen for DEPs between the CRC and normal control samples. The threshold for DEP screening was set at P < 0.05 and |log2FC| > 0.263. Hierarchical cluster analysis of the group samples was conducted with the pheatmap package (Version 1.0.8; https://cran.r-project.org/web/packages/pheatmap).

Functional enrichment analysis of DEPs

Functional enrichment analysis was performed for all the DEPs obtained from the three indexes to investigate their involved functions. KEGG pathway [13] and GO [14] functional analyses of these DEPs were performed using the R package of clusterProfiler (Version 3.2.11, http://www.bioconductor.org/packages/release/bioc/html/clusterProfiler.html). P < 0.05 was selected as the threshold value.

The construction of the Protein-Protein interaction (PPI) network

The STRING database (version: 10.0, http://www.string-db.org/) [15] was used to identify interactions between proteins. A required confidence (combined score) > 0.4 was set as the threshold for interaction, and the network was constructed with the Cytoscape [16] software.

Prediction of drug-gene interactions

In order to investigate whether the genes encoding DEPs could be targeted in CRC treatment or be used for drug development, we mined the existing information in the drug-gene interaction database (DGIdb2.0; http://www.dgidb.org/). With the preset filters of “FDA Approved”, drug-gene interaction pairs were predicted. Cytoscape (version: 3.2.0, http://www.cytoscape.org/) was used to construct the network.

Verification analysis

To further determine the clinical value of the key DEPs, we screened another 68 serum samples including 34 CRC samples (age; 61.5 ± 9.3; 21 men and 13 women) and 34 normal controls (age; 61.0 ± 11.2; 21 men and 13 women) with a custom commercial chip. The chip was added with chemiluminescence solution and tested directly by was detected by laser scanner and analyzed by Maxpix, and the samples was in duplicate. The limma package was used to screen for DEPs between the CRC and normal samples, and the threshold for DEP screening was set at P < 0.05 and |log2FC| > 0.585. Specifically, 15 proteins: AgRP, ANG-2, Dtk, EOT3, FGF-4, FGF-9, HCC-4, IL-16, IL-8, MIF, MSPa, TECK, TPO, TRAIL R3 and VEGF-D were used to construct the custom chip. Bioinformatic analysis with the chip used the same data preprocessing, identification of DEPs and hierarchical clustering as described above.

SVM model prediction

The expression levels of the significant DEPs in the two cohorts (24 samples in initial analysis and the 68 samples in verification analysis) were used as the eigenvalues for the e1071 SVM package [17]. Datasets (n = 92) were randomly divided into training set (n = 65) and independent testing set (n = 27) according to the probability ratio of 7:3. The ROC curve was then used to evaluate the classification and prediction efficiency of the data. The default parameters for the SVM model were: Kernel function: sigmoid kernel; Cross: 10-fold cross validation.

Survival analysis

Clinical information and mRNA-seq data for CRC were downloaded from the UCSC database (http://xena.ucsc.edu/), which contains the TCGA data. A total of 435 samples with survival information were obtained. The expression levels of the candidate genes corresponding to the DEPs were screened with the TCGA data, and then divided into high expression and low expression groups based on the optimal cutoff obtained by surv_cutpoint function. Survival analysis was determined using the K-M survival curve with a log-rank test, and P < 0.05 was used as the threshold.

Results

Data preprocessing

After the data was normalized with the quantile method, the internal structure of the data changed dramatically. The cumulative strength, average cumulative strength, and average net cumulative strength before and after data preprocessing are shown in Figure 1.

Figure 1. Results before and after data preprocessing for cumulative strength (a), average cumulative strength (b) and average net cumulative strength (c). The upper part of each picture is a boxplot and the lower part is the data density distribution curve. The intensity represents the protein expression value; Density represents the percentage of the total number

DEG screening and hierarchical clustering

Following data preprocessing, 6 DEPs (2 up-regulated and 4 down-regulated) were screened based on the cumulative intensity index (Table 1). Similarly, 6 DEPs (4 up-regulated and 2 down-regulated), and 4 DEPs (all down-regulated) were identified based on the mean cumulative strength and net cumulative strength, respectively. A volcano plot of the DEPs is shown in Figure 2.

Table 1. The DEPs under the cumulative intensity index, the mean cumulative strength and mean net cumulative strength

Indexes	Protein	logFC	AveExpr	t	P.Value
Cumulative intensity	dtk	1.2750	9.2921	2.5346	0.0175
	Eotaxin-3	1.2455	11.9783	3.7907	0.0008
	TRAIL-R3	−0.5561	11.1347	−2.2323	0.0342
	MIF	−0.6047	10.0441	−2.0709	0.0482
	MSP-a	−0.6571	13.3801	−2.1852	0.0379
	IL-16	−0.6700	11.1803	−2.1982	0.0368
Mean cumulative strength	Angiopoietin-2 (ANG-2)	1.1791	3.0044	3.0366	0.0054
	AgRP	0.8498	4.9343	2.7008	0.0121
	FGF-4	0.8018	4.9546	2.6135	0.0148
	FGF-9	0.6787	3.8830	2.1302	0.0429
	IL8	−0.6444	5.5619	−2.1175	0.0441
	TPO	−0.8790	3.6620	−2.3411	0.0273
Mean net cumulative strength	VEGF-D	−0.2809	4.4690	−2.1153	0.0421
	TECK	−0.3269	4.2303	−2.3720	0.0238
	HCC-4	−0.3536	6.0521	−2.5265	0.0166
	IL8	−0.5219	5.7839	−2.7077	0.0107

Figure 2. The volcano plot of differentially expressed proteins (DEPs) based on the cumulative strength (a), average cumulative strength (b) and average net cumulative strength (c). The blue dots represent down-regulated DEPs, and the red dots represent up-regulated DEPs. The x-coordinate is log₂FC, and the y-coordinate is log₁₀ (P value)

Functional enrichment analysis of the DEPs

All the DEPs obtained from the three indexes were chosen for subsequent analysis. As illustrated in Figure 3(a), there were 15 DEPs. Only IL-8 was differentially expressed by both the average cumulative strength and average net cumulative strength. GO analysis results showed that the positive regulation of the ERK1 and ERK2 cascades, cell chemotaxis, and leukocyte migration were significantly enriched with DEPs (Figure 3b). Furthermore, KEGG pathway analysis suggested cytokines-cytokine receptor interaction, and the chemokine signaling, Rap1, Ras, and MAPK signaling pathways were the main pathways enriched in DEPs (Figure 3c).

Figure 3. Functional enrichment of the differentially expressed proteins (DEPs). (a) Venn diagram of the DEPs; (b) Bubble diagram of Gene ontology (GO) terms enriched with DEPs. (c) Bubble diagram of pathways enriched in DEPs. The horizontal axis is the GeneRatio: the higher the ratio, the greater the fraction of proteins which are DEPs. The size of the dot is the number of DEPs. As the color changes from blue to red, the smaller the p value and the higher the corresponding significance

PPI network construction and drug-gene interactions

To further predict the potential interactions of these DEPs, the PPI network was constructed based on the STRING database. As shown in Figure 4(a), the network consisted of 12 nodes and 12 gene interaction pairs, of which, IL-8, CCL16, and TYRO3 had the highest degree of connectivity.

Figure 4. Protein-Protein Interaction (PPI) and drug-gene interaction networks. (a) Down-regulated DEPs are green and up-regulated DEPs are red. A blue edge shows that the protein was identified by a difference in average net cumulative strength, an orange edge shows that the protein was identified by a difference in mean cumulative strength, and a purple edge shows that the protein was identified by a difference in cumulative strength. The thickness of the line corresponds to the score of the interaction relationship with thicker lines corresponding to higher interaction scores. (b) The red oval is the key gene, the yellow rectangle is the drug, and the arrow points to the inhibitory effect

Based on the DGIdb, there were 10 drug-gene interaction pairs, including 5 genes and 10 drugs (Figure 4b). Of note, four drug-gene interactions had inhibitory effects.

Verification analysis of key DEPs with a custom chip

To further validate our findings, a custom chip was constructed with 15 DEPs and used to screen 34 CRC and 34 normal control samples. As illustrated in Figure 5, the heatmap suggested that the CRC samples and the normal controls could be distinguished based on the DEPs. Of the DEPs on the custom chip, the up-regulated DEPs: IL-8, MSPa and MIF, as well as the down-regulated DEPs: FGF-9, ANG-2 and AgRP were the key DEPs with statistical difference. Moreover, these key DEPs would be used as potential diagnostic markers or drug targets in SVM model prediction.

Figure 5. Bidirectional clustering heatmap. The color in the top column represents the grouping. Red is the healthy group and purple is the disease group. In the heat map, green represents low expression and red represents high expression

SVM model prediction and survival analysis

As illustrated in Figure 6, the ROC curve of the SVM model indicated that the 6 key DEPs could effectively distinguish the CRC samples from the normal controls, with the AUC of 0.933.

Figure 6. Receiver operating characteristic (ROC) curve for SVM model. The horizontal axis represents specificity, the vertical axis sensitivity, and the point is the cutoff threshold

Furthermore, survival analysis was performed for all DEPs in our study to investigate their prognostic value. Among the 15 DEPs, six DEPs were found to be associated with prognosis, including FGF-9, TPO, TNFRSF10C (also named TRAIL R3), TYRO3 (also named Dtk), CCL25 (also named TECK) and FGF-4 (Figure 7). It could be seen that patients with high expression of FGF-4 and TNFRSF10C had favorable prognosis than that of low expression. While patients with high expression of FGF-9, CCL25, TPO and TYRO3 had worse outcomes than that of low expression. The results suggested that these DEPs might be predictor for CRC prognosis. Of which, FGF-9 was a key DEPs in SVM diagnostic model, suggesting both diagnostic and prognostic value of FGF-9 in CRC.

Figure 7. Results for survival analysis. Survival curves shows the prognostic value of FGF9, TPO, TNFRSF10C, TYRO3, CCL25 and FGF4. Red represents the high-risk group and black the low-risk group

Discussion

Over the last decade, an increasing number of studies have been reported which investigate potential factors involved in the diagnosis and prognosis of CRC [18–20], but none is completely reliable. A protein profile can provide useful information on the pathology and molecular mechanism of a disease or disorder. In the present study, we investigated a series of DEPs by calculating the indexes of cumulative strength, average cumulative strength, and average net cumulative strength, and pathways associated with CRC development. A custom chip with the 6 key DEPs, IL-8, MSPa, MIF, FGF-9, ANG-2 and AgRP could effectively distinguished the CRC samples from normal controls. Survival analysis showed that high expression of FGF-9 was associated with worse prognosis than that of low expression. In addition, survival analysis also showed TPO, TRAIL R3, Dtk, TECK and FGF-4 were also showed prognostic value in CRC.

Tumorigenesis is a complex process that involves the interaction of various genes, proteins and signaling pathways. In the case of CRC, an increasing number of signaling pathways have been reported to play important roles in the progression of the disease [21,22]. Therefore, an analysis of the pathways related to CRC could provide a mechanistic basis for the progression of the disease. In the present study, we found that the DEPs were mainly found in cytokines-cytokine receptor interaction, and the chemokine signaling, Rap1, Ras and MAPK signaling pathways. Qi et al [21] also found that the differentially expressed genes in GSE113513 (obtained from GEO database) mainly found in the chemokine signaling pathway and cytokine-cytokine receptor interaction. A previous study suggested that the CXC family of chemokines could serve as prognostic or predictive biomarkers and possible drug targets for CRC [23]. Our data also revealed candidate proteins associated with the chemokine signaling pathway, which might provide a basis for research into the role of the DEPs in the formation of CRC. The RAS family has been reported to be among the most commonly mutated genes in all human malignancies including CRC. Coelho et al [24] suggested that oncogenic RAS signaling could promotes tumor immunoresistance by stabilizing PD-L1 mRNA. In CRC, the Ras signaling pathway is activated through different genes [25,26].

Additionally, we constructed the PPI and drug-gene interaction networks. Of the PPI network nodes, IL-8 was as the key protein with the highest degree of connectivity. It has been reported that polymorphisms in IL-8 are related to the risk of colorectal adenoma [27], and that IL-8 regulates cancer stem cell like properties in colon cells through induction of O-GlcNAc modification [28]. Our study revealed the interaction between IL-8 and other proteins such as CCL16, CCL25, TYRO3, CCL26, MIF, ANGPT2, and IL-16, which has not previously been reported. TPO was the key drug target for the drug-gene interactions. TPO expression has previously been found to be significantly different in CRC and healthy controls [29]. However, research on whether it could be a drug target is still limited, and our study provides novel potential therapeutic targets for CRC treatment.

Through the prognostic analysis of the DEPs, we found that FGF-4 could be a predictor of prognosis for patients with CRC. Currently, FGF-4 is widely studied in hepatocellular carcinoma, breast cancer, oral squamous cell carcinoma and esophageal squamous cell carcinoma. Huang [30] found that FGF-4 expression in combination with clinical stage could be used as a relatively accurate predictor of the 5-year probability of death and recurrence in ESCC patients. However, the role of FGF-4 in CRC development is still unclear. Moreover, in our study, the DEPs IL-8, MSPa, MIF, FGF-9, ANG-2 and AgRP could be used as potential diagnostic markers with an AUC of 0.933 based on the SVM model prediction. Of these proteins, the diagnostic values of MSPa, FGF-9 and AgRP in CRC have not previously been reported. Taken together, our data revealed a set of proteins (DEPs IL-8, MSPa, MIF, FGF-9, ANG-2 and AgRP) used in the panel in serum for distinguishing the CRC patients from normal controls. However, external validation with clinical samples is still required to confirm the diagnostic value of the identified proteins.

In recent year, great efforts had been made to identify biomarkers for the diagnosis of CRC. Circulating tumor DNA (ctDNA) shares primary tumor-specific genetic changes because that it is mainly released by apoptotic or necrotic tumor cells into circulation or actively secreted by circulating tumor cells [31]. Therefore, analysis of ctDNA is considered to have great potential to improve CRC diagnosis. However, there remains certain difficulties. Firstly, ctDNA accounts for only a minority of total circulating DNA which contains abundant cell-free DNA secreted from healthy cells. Secondly, the low concentration and high fragmentation rate of ctDNA in circulation present a great challenge to analysis [32,33]. Machine learning-based diagnosis had shown well classification performance in different diseases [34]. SVM diagnostic model with tumor marker had been reviewed to had potential clinical application value in the diagnosis of cancer, including CRC [35]. Study had reported that the genes screened by SVM model had the ability to accurately distinguish metastatic and non-metastatic CRC samples [36]. Therefore, we speculated that our SVM diagnostic model with the six key DEPs had potential clinical application value in CRC diagnosis.In conclusion, we conducted a bioinformatics analysis of the DEPs between CRC serum and normal control samples. The DEPs: AgRP, ANG-2, Dtk, EOT3, FGF-4, FGF-9, HCC-4, IL-16, IL-8, MIF, MSPa, TECK, TPO, TRAIL R3 and VEGF-D, as well as cytokines-cytokine receptor interaction, and the chemokine signaling, Rap1, Ras, and MAPK signaling pathways, which are enriched in them, may be involved in CRC development. Moreover, the DEPs: IL-8, MSPa, MIF, FGF-9, ANG-2 and AgRP could be potential diagnostic markers. Further survival analysis showed FGF9, TPO, TNFRSF10C, TYRO3, CCL25 and FGF4 could be biomarkers to predict the prognosis of CRC patients. Taken together, these findings have the important clinical significance because they could improve our understanding of the pathogenesis and molecular mechanisms of CRC; moreover, they provide effective tools for the prognosis and diagnosis of CRC. However, further clinical study is still needed to confirm the function of the identified proteins.

Highlights

1. A series of DEPs was found in CRC.

2. Pathways associated with DEPs were enriched in CRC.

3. A series of DEPs related to CRC diagnosis was identified.

Acknowledgments

None.

Disclosure statement

The authors declare that there is no conflict of interest.