Linear modelling at the probe-level:
The number of significant genes present in each stage-control pair from the Tier-I contrasts is shown in Figure 1(i). Using the top 100 DM genes of the linear model (given in Supplementary Information S3), we found a clear separation between controls and stage samples (Figure 1(ii)). The top genes in each stage (by adjusted p-value of contrast with control) are shown in Table 2, with |lfc M-value| and inferred regulation status. Figure 2 shows boxplots of stagewise methylation levels for two representative genes: (1) TMEM179, mutations in which could cause MSI [40]; and (2) MEOX2 whose promoter methylation status is a known CRC marker [41].The top four genes of each stage were used to construct a stagewise methylation heatmap(Figure 3). The stagewise methylation patterns of the top five linear model genes are also shown, in Figure 4. It is notable that a naturally occuring read-through fusion protein GPR75-ASB3 is the top linear model gene with significant differential expression in all stages relative to the control. GPR75-ASB3 is positively differentially expressed in the lung as well as different keratinocyte cell types, and evidence is emerging of its role in other cancers [42]. In this light, GPR75-ASB3 could play a significant role in colorectal cancers which are of epithelial origin. The top 100 significant stage-specific genes, listed in S3, were used in the consensus analysis.
Table 2. Top ten genes of the linear model at the probe level. The log fold-change of M-value of the probe in each stage relative to the controls, followed by p-value adjusted for the false discovery rate, and the methylation status of the gene in the cancer stages with respect to the control
ID
|
Stage I
lfc (β1)
|
Stage II lfc (β2)
|
Stage III lfc (β3)
|
Stage IV lfc (β4)
|
Adj. p-val
|
Methylation status
|
GPR75-ASB3
|
2.28
|
2.19
|
2.16
|
2.32
|
1E-82
|
Hyper
|
TM4SF19
|
-3.63
|
-3.58
|
-3.72
|
-3.71
|
1E-82
|
Hypo
|
CNRIP1
|
2.74
|
2.60
|
2.68
|
2.97
|
1E-78
|
Hyper
|
PDE4A
|
1.68
|
1.58
|
1.60
|
1.71
|
1E-71
|
Hyper
|
KRTAP11-1
|
-2.36
|
-2.30
|
-2.37
|
-2.40
|
1E-70
|
Hypo
|
ADHFE1
|
3.15
|
2.97
|
3.00
|
3.43
|
1E-69
|
Hyper
|
FAM123A
|
3.56
|
3.18
|
3.43
|
3.90
|
1E-69
|
Hyper
|
KHDRBS2
|
2.30
|
2.16
|
2.10
|
2.34
|
1E-68
|
Hyper
|
AJAP1
|
2.52
|
2.44
|
2.46
|
2.64
|
1E-68
|
Hyper
|
NALCN
|
2.96
|
2.80
|
2.94
|
3.25
|
1E-68
|
Hyper
|
Linear modelling at the gene-level (averep):
The genes with more than one probe were averaged to a single methylation value, which was then further analyzed. The number of genes present in each stage-control pair from the Tier-I contrasts is shown in Figure 5(i). Using the top 100 genes of the linear model (given in Supplementary Information S4), we found a clear separation between controls and stage samples (Figure 5(ii)). The top genes in each stage (by adjusted p-value of contrast with control) are shown in Table 3, with |lfc M-value| and inferred regulation status. Figure 6 shows the boxplots of stagewise methylation levels for two representative genes, NALCN and GLRX. Mutations in NALCN have been reported in sporadic CRC [43]; here NALCN is seen to be significantly hypermethylated, indicating the same outcome (loss of function) could be effected in multiple ways. GLRX is a target of the activating transcription factor MEOX2 [44]. The top four genes of each stage were used to construct a stagewise methylation heatmap (Figure 7). The stagewise methylation patterns of the top five linear model genes are also shown, in Figure 8. It is observed that LY6H showed both hypermethylation and hypomethylation when compared to the controls, indicating the role of experimentation necessary to clarify its role in colorectal cancer progression. The top significant 100 genes of each stage, listed in S4, were used for the consensus analysis.
Table 3. Top ten genes of the linear model at the gene level, using average values of methylation. The log fold-change of M-value of the gene in each stage (relative to the control) is given, followed by p-value adjusted for the false discovery rate and the methylation status of the gene in the cancer stages with respect to the control. A consistent methylation pattern is observed for all the top genes.
ID
|
Stage I
lfc (β1)
|
Stage II lfc (β2)
|
Stage III lfc (β3)
|
Stage IV lfc (β4)
|
Adj. p-val
|
Methylation status
|
TM4SF19
|
-3.63
|
-3.57
|
-3.72
|
-3.71
|
1E-82
|
Hypo
|
GPR75-ASB3
|
2.28
|
2.19
|
2.15
|
2.32
|
1E-82
|
Hyper
|
CNRIP1
|
2.74
|
2.60
|
2.67
|
2.97
|
1E-77
|
Hyper
|
KRTAP11-1
|
-2.36
|
-2.30
|
-2.38
|
-2.40
|
1E-70
|
Hypo
|
ADHFE1
|
3.15
|
2.96
|
2.99
|
3.43
|
1E-69
|
Hyper
|
FAM123A
|
3.56
|
3.18
|
3.42
|
3.89
|
1E-68
|
Hyper
|
AJAP1
|
2.53
|
2.44
|
2.46
|
2.64
|
1E-67
|
Hyper
|
NALCN
|
2.96
|
2.79
|
2.95
|
3.25
|
1E-65
|
Hyper
|
IRF4
|
1.99
|
1.83
|
1.89
|
2.13
|
1E-65
|
Hyper
|
PRKAR1B
|
3.38
|
3.13
|
3.24
|
3.50
|
1E-65
|
Hyper
|
Association with phenotype
The ANOVA from CpGassoc yielded p-values and log fold-changes, which were used to identify significant genes for each stage using the criteria given in Methods (Figure 9). The top 100 genes of each stage from this analysis (given in Supplementary Information S5) were used for the consensus investigation.
DMP analysis with ChAMP
The summary features of the β matrix dataset were evaluated using ChAMP (Figure 10). The DMPs were identified using CHAMP analysis from the β matrix. All the inter-stage contrasts yielded null results (i.e, no significant genes), except for stageII – stageIV contrast. Due to this, the top 100 DMPs from the stage vs control contrasts were used for the consensus analysis directly. Contrasts that showed significant DMPs were subjected to a further DMR analysis, to enable identification of DM genes. The stage-salient DMR regions (genes) determined are provided in Supplementary Information S6, and summarized in Table 4. The stage-II vs stage-IV DMR contrast yielded three genes, namely PLAG1, SOCS2, and NNAT. It is observed that these genes might be critical players in the transition to malignancy. Interestingly, some genes were differentially methylated in all the stagewise contrasts with the control; such genes are differentially methylated agnostic of stage, and could serve as valuable drug targets for CRC therapy. The top such genes included EYA4, WT1, DCC, RP11, GATA4, MSX1, DLX5, BNC1, WT1-AS, and ZIM2. A total of 31 such genes were identified and tabulated in Supplementary Information S7. The DMPs and DMRs from the analysis were subjected to GSEA and these results could also be found in Supplementary Information S6. Figure 11 shows representative DMP and DMR plots using MethylMix.
Table 4. Contrast-wise counts of DM probes and DM regions.
Contrast
|
DMPs
|
DMRs
|
Control and Stage 1
|
11045
|
34
|
Control and Stage 2
|
11254
|
35
|
Control and Stage 3
|
11254
|
36
|
Control and Stage 4
|
11108
|
34
|
Stage 2 and Stage 4
|
404
|
3
|
Methylation and Gene Expression Correlation analysis
Mixture models of genes, indicative of the number of methylation states, were constructed using MethylMix, and the top three genes from an overall cancer vs control comparison are shown in Figure 12. The estimated correlation between the methylation levels and actual gene expression for the same genes is depicted in Figure 13. Genes were differentially methylated and designated as ‘driver’ genes if the p-value of the contrast being studied was significant. The calculated differential methylation (DM) values from stage vs control contrasts ranged from -0.7 to +0.8, and genes were classified as hyper- or hypo-methylated based on the DM value. There were 209, 441, 275, and 134 driver genes in each of the contrasts with the controls (stage-I, stage-II, stage-III and stage-IV, respectively). All between-stages contrasts yielded null DM genes. The results from this analysis, including driver genes for all the contrasts, are provided in Supplementary Information S8. Top 100 genes from each comparison were taken forward for the consensus analysis. Certain genes emerged common to all the four comparisons, indicating stage-agnostic differential methylation events. The top such genes included CCDC88B, C1orf59, CHFR, ZP2, HOXA9, ELF5, FAM50B, MUC17, TBX20, and VSIG2. Stage-agnostic genes hold promise as therapeutic targets for the treatment of colorectal cancer; the complete list of 56 stage-agnostic genes arising out of the MethylMix analysis is provided in Supplementary File S9.
BioMethyl analysis
The significant stage-specific DEGs identified by this BioMethyl are shown in Figure 14. Top 100 genes of each stage from this analysis were taken for consensus analysis. The stage-specific genes from this analysis are presented in the Supplementary Information S10.
Stage-salient consensus biomarkers
The top 100 significantly differentially-expressed genes of each stage from all the methods discussed above (collated in Supplementary Information S11) were used for the consensus determination. The consensus analysis yielded seven stage-salient DMGs: one stage-I gene (FBN1), one stage-II gene (FOXG1), one stage-III gene (HCN1) and four stage-IV genes (NELL1, ZNF135, FAM123A, LAMA1). Each of these stage-salient genes presented an |lfc M-value| > 0.4 with respect to the other stages, validating their salience. Figures 15,16 represent boxplots of the consensus biomarkers, and Table 5 presents a summary of the consensus analysis. Gene ontology (GO) analysis [45] of the consensus biomarkers yielded processes related to structural integrity of cell division processes, immunity dysfunction, and cell migration (Table 6). Detailed GO results are presented in the Supplementary Information S12.
Table 5. Stage-salient biomarkers. The results of the consensus analysis and univariate survival analysis are summarized. All the biomarkers showed hypermethylation, indicating that an inferred common down-regulation of gene expression. Univariate survival analysis does not show significance for NELL1 by itself.
HGNC ID
|
Gene Name
|
Methods in agreement
|
Salience
|
Meth. status
|
P-values from modelling
|
P-values from univariate survival analysis
|
M value
|
Averep
|
Cox analysis
|
KaplanMeier
|
3603
|
FBN1
|
Averep, CHAMP
|
I
|
Hyper
|
0.310
|
0.040
|
0.036
|
0.025
|
3811
|
FOXG1
|
Averep, Mvalue, CHAMP, Methylmix
|
II
|
Hyper
|
1E-16
|
0.003
|
0.019
|
0.037
|
4845
|
HCN1
|
Averep, Mvalue, CHAMP
|
III
|
Hyper
|
1E-17
|
0.022
|
0.031
|
0.059
|
7756
|
NELL1
|
Mvalue, CHAMP, Methylmix
|
IV
|
Hyper
|
1E-68
|
0.061
|
0.283
|
0.27
|
12919
|
ZNF135
|
Mvalue, CHAMP, Methylmix
|
IV
|
Hyper
|
1E-76
|
0.062
|
0.096
|
0.084
|
26360
|
FAM123A
|
Mvalue, CHAMP, Methylmix
|
IV
|
Hyper
|
1E-115
|
0.097
|
0.30
|
0.28
|
6481
|
LAMA1
|
Mvalue, CHAMP, Methylmix
|
IV
|
Hyper
|
1E-86
|
0.297
|
0.052
|
0.051
|
Table 6. GO analysis of stage-salient genes in the order of decreasing significance (i.e, increasing p – value). Ontology: Cellular Compartment (CC), Molecular Function (MF), Biological Process (BP).
GO ID
|
Term
|
Ontology
|
p-value
|
GO:1990047
|
spindle matrix
|
CC
|
0.0001
|
GO:0030109
|
HLA-B specific inhibitory MHC class I receptor activity
|
MF
|
0.0003
|
GO:0032396
|
inhibitory MHC class I receptor activity
|
MF
|
0.0006
|
GO:0042609
|
CD4 receptor binding
|
MF
|
0.0012
|
GO:0032393
|
MHC class I receptor activity
|
MF
|
0.0013
|
GO:0050930
|
induction of positive chemotaxis
|
BP
|
0.0016
|
GO:0050927
|
positive regulation of positive chemotaxis
|
BP
|
0.0033
|
GO:0050926
|
regulation of positive chemotaxis
|
BP
|
0.0034
|
GO:0008608
|
attachment of spindle microtubules to kinetochore
|
BP
|
0.0043
|
GO:0007094
|
mitotic spindle assembly checkpoint
|
BP
|
0.0044
|
Survival analysis:
We constructed independent prognostic models of the stage-salient genes and identified the prognostically significant biomarkers as FBN1, FOXG1, HCN1, and LAMA1. The corresponding univariate Kaplan-Meier plots are shown in Fig. 17. Rational combinations of stage-salient genes, termed ColoRectal cancer Signatures (CRS), were modelled using multivariate Kaplan-Meier regression, to yield a risk score. Risk scores were then used to estimate survival-effect significance, as described in Methods. The results of this exercise are summarised in Table 7. We found that CRS12 signature (consisting of FBN1 and FOXG1) yielded significant risk scores in the multivariate Kaplan-Meier analysis, and both CRS12 and CRS34 (consisting of HCN1, NELL1, ZNF135, FAM123A, LAMA1) were significant in estimating overall survival (prognosis p-value ≤ 0.02) (Figure 18). Supplementary Information S13 provides survival plots of all possible signatures. At the end of our analysis pipeline, CRS12 passed all the filters and emerged as a significant early-stage panel for CRC prognosis.
Table 7. Summary of selected multivariate prognostic models. Weight denotes the coefficient in the multivariate model. The ultimate significant signature is highlighted.
Signature
|
Stages
|
Biomarker
|
Weight
|
P-value
|
Multivariate model
|
Prognosis
|
CRS12
|
I+II
|
FBN1
|
-0.62
|
0.015
|
0.005
|
FOXG1
|
-1.05
|
CRS34
|
III+IV
|
NELL1
|
0.1
|
0.172
|
0.02
|
ZNF135
|
-0.21
|
FAM123A
|
-0.23
|
LAMA1
|
-0.39
|
HCN1
|
-1.1
|
CRS234
|
II+III+IV
|
FOXG1
|
-0.99
|
0.087
|
0.032
|
HCN1
|
-1.07
|
NELL1
|
-0.1
|
ZNF135
|
-0.22
|
FAM123A
|
-0.37
|
LAMA1
|
-0.27
|