Identification of colorectal cancer associated biomarkers: an integrated analysis of miRNA expression

Abstract

Colorectal cancer is one of the leading causes of cancer-related deaths worldwide. This complex disease still holds severe problems concerning diagnosis due to the high invasiveness nature of colonoscopy and the low accuracy of the alternative diagnostic methods. Additionally, patient heterogeneity even within the same stage is not properly reflected in the current stratification system. This scenario highlights the need for new biomarkers to improve non-invasive screenings and clinical management of patients. MicroRNAs (miRNAs) have emerged as good candidate biomarkers in cancer as they are stable molecules, easily measurable and detected in body fluids thus allowing for non-invasive diagnosis and/or prognosis. In this study, we performed an integrated analysis first using 4 different datasets (discovery cohorts) to identify miRNAs associated with colorectal cancer development, unveil their role in this disease by identifying putative targets and regulatory networks and investigate their ability to serve as biomarkers. We have identified 26 differentially expressed miRNAs which interact with frequently deregulated genes known to participate in commonly altered pathways in colorectal cancer. Most of these miRNAs have high diagnostic power, and their prognostic potential is evidenced by panels of 5 miRNAs able to predict the outcome of stage II and III colorectal cancer patients. Notably, 8 miRNAs were validated in three additional independent cohorts (validation cohorts) including a plasma cohort thus reinforcing the value of miRNAs as non-invasive biomarkers.

Keywords: biomarker; colorectal cancer; diagnosis; microRNAs; prognosis.

Figure 1

Number of miRNAs throughout data processing and statistical analysis. The number of miRNAs after each major event (Listwise case deletion and Statistical analysis) is shown below each dataset. Different colours represent: yellow TCGA, green GSE33215, orange GSE18392 and reddish-purple GSE30454. From the top, the numbers of miRNAs represent: the initially number of miRNAs present in each dataset before any processing step; the number of miRNAs after removing the ones with more than 50% information missing (modified listwise case deletion); and the number of miRNAs after performing all the statistics and used to perform the identification of differently expressed miRNAs across datasets.

Figure 2

Venn diagram of the differently expressed miRNAs between datasets. Allocation of the 1043 differently expressed miRNAs found between the 4 datasets used in this work. Each dataset is represented by a colour, TCGA (yellow), GSE30545 (reddish purple), GSE33215 (green) and GSE18392 (orange). The number in each overlay of datasets represents the common miRNAs between those datasets.

Figure 3

Differentially expressed miRNAs between colorectal cancer and normal tissue samples. The log₂(FC) values calculated for each dataset are reported with red scale boxes for upregulated miRNAs and blue scale boxes for the downregulated miRNAs. White boxes represent the inexistence of the miRNA on the dataset. Only the miRNAs differentially expressed in at least 2 datasets while simultaneously being present in the TCGA dataset are displayed.

Figure 4

Chord-dendrogram of the interactions between differently expressed miRNAs and altered genes in CRC. The interactions between the 25 differently expressed miRNAs identified and the 173 experimentally validated target genes were obtained through MirTarBase. Target genes’ role in CRC is according to CoReCG database. All targets are represented in grey, while each miRNA is represented by one colour. MiRNA:target interactions are represented by a line of the same colour of the respective miRNA. The size of rectangle next to the name of the miRNAs and target genes is proportional to the number of interactions they perform.

Figure 5

Circular barplot evidencing the expression of the 172 experimentally validated miRNA target genes. The expression values of the 172 genes were obtained from the TCGA colon and rectal cohorts and the log₂(FC) values between the primary tumour and normal tissue samples were calculated. The log₂(FC) value of each gene is given by the length and tone of each coloured bar in accordance with each gene regulation status (red colours – upregulated; blue colours- downregulated).

Figure 6

KEGG pathway analysis – interaction between the 25 differently expressed miRNAs and the frequently altered pathways in colorectal cancer. The 173 miRNA-target genes were used to perform the enrichment pathway analysis using KEEG at DAVID database. The miRNA-target genes found to be involved in CRC were crossed with the genes intervening in each signalling pathway in order to identify the miRNAs affecting each pathway. The miRNAs interaction with each KEGG Pathway is reported in pink.

Figure 7

Areas under the ROC curve (AUC) of the 25 differentially expressed miRNAs between colorectal cancer and normal tissues samples. The miRNAs AUC values in each of the datasets TCGA (Yellow), GSE18392 (orange), GSE33125 (Green) and GSE30454 (reddish purple) are reported as blue scale boxes. MiRNAs with AUC = 1 were considered perfect diagnostic biomarkers, 0.9 < AUC < 1 highly accurate, 0.7 < AUC ≤ 0.9 moderately accurate and 0.5 < AUC ≤ 0.7 less accurate [76].

Figure 8

Individual miRNA overall survival (OS) prognosis for stage III and stage IV patients. (A) Stage III OS Kaplan-Meier curve for miR-133b (p-value = 0.047, Log rank; HR = 2.28) (B) Stage IV OS Kaplan-Meier curve for miR-129-5p (p-value = 0.021, Log rank; HR = 4.23). Time is represented in years. Higher (in red) and Lower (in blue) expression groups represent the patients with miRNA expression above and below miRNAs’ median expression, respectively. Censored data is represented by small plus signs in each group. The number of patients at risk for each group and per time point is shown in the table below each graph. HR, hazard ratio.

Figure 9

Panels of 5 miRNAs for overall survival (OS) and recurrence free survival (RFS) prognosis of stage II patients. (A) Stage II OS Kaplan-Meier curve based on miR-320b - miR-326 - miR-331-3p - miR-339-5p - miR-484 (p-value = 0.032, Log rank test; HR= 5.23) (B) Stage II recurrence free survival RFS Kaplan-Meier curve based on miR-320a - miR-324-5p - miR-324-3p - miR-423-3p - miR-484 (p-value = 0.043, Log rank test). Time is represented in years. Higher (in red) and Lower (in blue) expression groups represent the group of patients with miRNA expression above and below miRNAs median expression, respectively. Censored data is represented by small plus signs in each group. The number of patients at risk for each group and per time point is shown in the table below each graph. HR, hazard ratio.

Figure 10

Panels of 5 miRNAs for overall survival (OS) and recurrence free survival (RFS) prognosis of stage III patients. (A) Stage III OS Kaplan-Meier curve based on miR-324-5p - miR-324-3p - miR-331-3p - miR-484 - miR-486-5p (p-value = 0.020, Log rank test; HR= 8.15). (B) Stage III recurrence free survival RFS Kaplan-Meier curve based on miR-299-5p - miR-324-5p - miR-324-3p - miR-331-3p - miR-484 (p = 0.030, Log rank test). Time is represented in years. Higher (in red) and Lower (in blue) expression groups represent the group of patients with miRNA expression above and below miRNAs median expression, respectively. Censored data is represented by small plus signs in each group. The number of patients at risk for each group and per time point is shown in the table below each graph. HR, hazard ratio.

Figure 11

Validation analysis for the 25 downregulated miRNAs in CRC. (A) The log₂(FC) values calculated for each dataset are reported with red scale boxes for upregulated miRNAs and blue scale boxes for the downregulated miRNAs. White boxes represent the inexistence of the miRNA on the dataset. (B) The miRNAs AUC values in each of the datasets GSE115513, GSE41655 and GSE71008 are reported as blue scale boxes. MiRNAs with AUC = 1 were considered perfect diagnostic biomarkers, 0.9 < AUC < 1 highly accurate, 0.7 < AUC ≤ 0.9 moderately accurate and 0.5 < AUC ≤ 0.7 less accurate [76]. (C) Stage III OS Kaplan-Meier curve based on miR-486-5p - miR-330-3p - miR-375 (p-value = 0.025, Log rank test; HR= 4.01). Time is represented in years. Higher (in red) and Lower (in blue) expression groups represent the group of patients with miRNA expression above and below miRNAs median expression, respectively. Censored data is represented by small plus signs in each group. The number of patients at risk for each group and per time point is shown in the table below each graph. HR, hazard ratio.