Association analyses identify 31 new risk loci for colorectal cancer susceptibility.

Colorectal cancer (CRC) is a leading cause of cancer-related death worldwide, and has a strong heritable basis. We report a genome-wide association analysis of 34,627 CRC cases and 71,379 controls of European ancestry that identifies SNPs at 31 new CRC risk loci. We also identify eight independent risk SNPs at the new and previously reported European CRC loci, and a further nine CRC SNPs at loci previously only identified in Asian populations. We use in situ promoter capture Hi-C (CHi-C), gene expression, and in silico annotation methods to identify likely target genes of CRC SNPs. Whilst these new SNP associations implicate target genes that are enriched for known CRC pathways such as Wnt and BMP, they also highlight novel pathways with no prior links to colorectal tumourigenesis. These findings provide further insight into CRC susceptibility and enhance the prospects of applying genetic risk scores to personalised screening and prevention.

Introduction

Many colorectal cancers (CRC) develop in genetically susceptible individuals[1] and genome-wide association studies (GWAS) of CRC have thus far reported 43 SNPs mapping to 40 risk loci in European populations[2,3]. In Asians, 18 SNPs mapping to 16 risk loci have been identified[4,5], a number of which overlap with those reported in Europeans. Collectively across ethnicities GWAS has provided evidence for 53 unique CRC susceptibility loci. While much of the heritable risk of CRC remains unexplained, statistical modelling indicates that further common risk variants remain to be discovered[6].

To gain a more comprehensive insight into CRC aetiology, we conducted a GWAS meta-analysis that includes additional, unreported datasets. We examine the possible gene regulatory mechanisms underlying all GWAS risk loci by analysing in situ promoter Capture Hi-C (CHi-C) to characterise chromatin interactions between predisposition loci and target genes, examine gene expression data and integrate these data with chromatin immunoprecipitation-sequencing (ChIP-seq) data. Finally, we quantify the contribution of the loci identified in this study, together with previously identified loci to the heritable risk of CRC and estimate the sample sizes required to explain the remaining heritability.

Results

Association analysis

Thus far, studies have identified 61 SNPs that are associated with CRC risk in European and Asian populations (Supplementary Data 1). To identify additional CRC risk loci, we conducted five new CRC GWAS, followed by a meta-analysis with 10 published GWAS totalling 34,627 cases and 71,379 controls of European ancestry under the auspices of the COGENT (COlorectal cancer GENeTics) consortium[7] (Fig. 1, Supplementary Data 2). Following established quality control measures for each dataset[8] (Supplementary Data 3), the genotypes of over 10 million SNPs in each study were imputed, primarily using 1000 Genomes and UK10K data as reference (see Methods). After filtering out SNPs with a minor allele frequency <0.5% and imputation quality score <0.8, we assessed associations between CRC status and SNP genotype in each study using logistic regression. Risk estimates were combined through an inverse-variance weighted fixed-effects meta-analysis. We found little evidence of genomic inflation in any of the GWAS datasets (individual λGC values 1.01–1.11; meta-analysis λ1000 = 1.01, Supplementary Figure 1).

Fig. 1

Study design

Excluding flanking regions of 500 kb around each previously identified CRC risk SNP, we identified 623 SNPs associated with CRC at genome-wide significance (logistic regression, P < 5 × 10−8). After implementing a stepwise model selection, these SNPs were resolved into 31 novel risk loci, with 27 exhibiting Bayesian False Discovery Probabilities (BFDPs)[9] <0.1 (Table 1, Fig. 2, Supplementary Figure 2). The association at 20q13.13 (rs6066825) had only been previously identified as significant in a multi-ethnic study[10]. Two new associations (rs3131043 and rs9271770) were identified within the 6p21.33 major histocompatibility (MHC) region, with rs3131043 located 470 kb 5′ of HLA-C, and rs9271770 located between HLA-DRB1 and HLA-DQA1. Imputation of the MHC region using SNP2HLA[11] provided no evidence for additional MHC risk loci.

Table 1

Summary results for the new colorectal cancer risk loci in Europeans

SNP	Cytoband	Position (bp, GRCh37)	Risk/alt allele	RAF	OR	95% CI	P-value	BFDP	I2 (%)	P het	Average info score
rs61776719	1p34.3	38,461,319	C/A	0.45	1.07	(1.05; 1.10)	2.19 × 10−10	1.98 × 10−3	1	0.44	0.89
rs12143541	1p32.3	55,247,852	G/A	0.15	1.10	(1.06; 1.13)	9.44 × 10−10	7.44 × 10−3	16	0.28	0.95
rs11692435	2q11.2	98,275,354	G/A	0.90	1.12	(1.07; 1.16)	1.22 × 10−8	0.079	29	0.14	0.97
rs11893063	2q33.1	199,601,925	A/G	0.47	1.07	(1.04; 1.09)	9.34 × 10−9	0.069	43	0.04	0.96
rs7593422	2q33.1	200,131,695	T/A	0.55	1.07	(1.05; 1.10)	3.56 × 10−11	3.50 × 10−4	15	0.28	0.99
rs9831861	3p21.1	53,088,285	G/T	0.59	1.07	(1.05; 1.09)	4.17 × 10−10	3.72 × 10−3	0	0.87	0.99
rs12635946	3q13.2	112,916,918	C/T	0.62	1.08	(1.06; 1.10)	1.02 × 10−11	1.03 × 10−4	11	0.33	0.97
rs17035289	4q24	106,048,291	T/C	0.83	1.10	(1.07; 1.13)	2.73 × 10−10	2.30 × 10−3	0	0.95	1.00
rs75686861	4q31.21	145,621,328	A/G	0.10	1.12	(1.08; 1.16)	1.76 × 10−9	0.014	0	0.49	0.92
rs2070699	6p24.1	12,292,772	T/G	0.48	1.07	(1.04; 1.09)	3.88 × 10−9	0.031	29	0.14	0.95
rs3131043	6p21.33	30,758,466	G/A	0.43	1.07	(1.05; 1.1)	2.67 × 10−8	0.159	60	0.01	0.91
rs9271770	6p21.32	32,594,248	A/G	0.81	1.08	(1.05; 1.11)	3.60 × 10−8	0.192	0	0.91	0.93
rs6928864	6q21	105,966,894	C/A	0.91	1.13	(1.09; 1.19)	1.37 × 10−8	0.094	0	0.73	0.98
rs10951878	7p12.3	46,926,695	C/T	0.49	1.06	(1.04; 1.09)	1.10 × 10−8	0.080	0	0.65	0.99
rs3801081	7p12.3	47,511,161	G/A	0.68	1.08	(1.06; 1.11)	2.00 × 10−11	1.96 × 10−4	50	0.01	1.00
rs1412834	9p21.3	22,110,131	T/C	0.50	1.08	(1.06; 1.11)	4.13 × 10−14	5.05 × 10−7	14	0.30	1.00
rs4450168	11p15.4	10,286,755	C/A	0.17	1.10	(1.06; 1.13)	1.24 × 10−8	0.079	0	0.81	0.86
rs7398375	12q13.3	57,540,848	C/G	0.72	1.09	(1.06; 1.13)	3.91 × 10−10	3.23 × 10−3	0	0.93	0.83
rs12427600	13q13.3	37,460,648	C/T	0.24	1.09	(1.06; 1.11)	5.43 × 10−11	5.01 × 10−4	0	0.81	0.99
rs45597035	13q22.1	73,649,152	A/G	0.64	1.08	(1.05; 1.10)	2.16 × 10−10	1.94 × 10−3	0	0.53	0.96
rs1330889	13q22.3	78,609,615	C/T	0.87	1.11	(1.07; 1.14)	6.50 × 10−10	5.25 × 10−3	0	0.59	0.97
rs7993934	13q34	111,074,915	T/C	0.65	1.08	(1.05; 1.10)	3.03 × 10−11	2.94 × 10−4	0	0.55	0.99
rs4776316	15q22.31	67,007,813	A/G	0.73	1.08	(1.05; 1.10)	1.11 × 10−8	0.076	22	0.21	0.95
rs10152518	15q23	68,177,162	G/A	0.19	1.08	(1.05; 1.11)	3.24 × 10−8	0.180	0	0.84	0.97
rs7495132	15q26.1	91,172,901	T/C	0.12	1.11	(1.07; 1.14)	7.92 × 10−10	6.34 × 10−3	29	0.14	0.99
rs61336918	16q23.2	80,007,266	A/T	0.29	1.09	(1.06; 1.12)	2.04 × 10−12	2.14 × 10−5	0	0.90	0.99
rs1078643	17p12	10,707,241	A/G	0.77	1.09	(1.06; 1.12)	4.14 × 10−11	3.81 × 10−4	0	0.56	0.92
rs285245	19p13.11	16,420,817	T/C	0.11	1.11	(1.07; 1.15)	3.71 × 10−8	0.195	2	0.42	0.91
rs12979278	19q13.33	49,218,602	T/C	0.53	1.07	(1.05; 1.09)	6.11 × 10−10	5.35 × 10−3	15	0.28	0.96
rs6066825	20q13.13	47,340,117	A/G	0.65	1.10	(1.08; 1.13)	3.82 × 10−17	5.67 × 10−10	0	0.49	0.99
rs3787089	20q13.33	62,316,630	C/T	0.32	1.07	(1.05; 1.10)	5.80 × 10−9	0.043	0	0.80	0.96
Associations previously only identified in Asian populations
rs639933	5q31.1	134,467,751	C/A	0.38	1.07	(1.05; 1.10)	1.14 × 10−9	9.50 × 10−3	0	0.73	0.98
rs6933790	6p21.1	41,672,769	T/C	0.83	1.10	(1.07; 1.14)	3.65 × 10−10	3.03 × 10−3	21	0.23	0.91
rs704017	10q22.3	80,819,132	G/A	0.60	1.10	(1.08; 1.13)	2.96 × 10−16	4.15 × 10−9	23	0.21	0.95
rs12255141	10q25.2	114,294,892	G/A	0.10	1.11	(1.07; 1.15)	2.97 × 10−9	0.022	0	0.81	0.96
rs10849438	12p13.31	6,412,036	G/T	0.12	1.12	(1.08; 1.16)	1.04 × 10−10	9.49 × 10−4	21	0.23	0.95
rs73975588	17p13.3	816,741	A/C	0.87	1.10	(1.06; 1.13)	8.71 × 10−9	0.058	33	0.11	0.97
rs9797885	19q13.2	41,873,001	G/A	0.71	1.08	(1.05; 1.10)	2.77 × 10−10	2.43 × 10−3	0	0.70	0.99
rs6055286	20p12.3	7,718,045	A/G	0.15	1.11	(1.07; 1.14)	9.69 × 10−11	8.61 × 10−4	50	0.02	0.97
rs2179593	20q13.12	42,660,286	A/C	0.72	1.07	(1.05; 1.10)	4.62 × 10−9	0.035	0	0.67	0.97

BFDP calculated using prior = 10−5 and maximum relative risk = 1.2

RAF risk allele frequency in Europeans, OR odds ratio, CI confidence interval, BFDP Bayesian False Discovery Probability, I2 proportion of the total variation due to heterogeneity, PP-value for heterogeneity

Fig. 2

Manhattan plot showing all loci containing genetic risk variants independently associated with colorectal cancer risk at P < 5 × 10−8 in European populations. SNPs on the left of the plot are new SNPs identified in this study, and SNPs on the right were identified in previous studies and replicated at genome-wide significance in this study. The 79 risk SNPs consisted of 31 previously reported SNPs, 39 new risk SNPs, and nine SNPs previously identified in Asian but not in European populations (denoted in dark gold). Dotted lines indicate SNPs that were identified as independent through conditional analysis. Square brackets indicate the position of the sentinel SNP relative to nearest genes (“gene1-[]-gene2” for intergenic SNPs, “[gene]” for intragenic SNPs). The distance from the sentinel SNP to each gene is proportional to the number of dashes. The red line indicates the genome-wide significance threshold. The x-axis represents the −log10P-values of the SNPs, and the y-axis represents the chromosomal positions

We confirmed 28 of the 40 risk loci for CRC published as genome-wide significant in Europeans (i.e. P < 5 × 10−8) (Supplementary Data 1). For four previously reported risk loci[2,12-14], we observed associations that were just below genome-wide significance (3q26.2, rs10936599, P = 1.41 × 10−7; 12p13.32, rs3217810, P = 1.09 × 10−6; 16q22.1, rs9929218, P = 4.96 × 10−7; 16q24.1, rs2696839, 1.28 × 10−6). In contrast, there was limited support in our current study for eight of the associations previously reported by others[2,10,15-17] (2q32.3, rs11903757, P = 0.23; 3p14.1, rs812481, P = 0.44; 4q22.2, rs1370821, P = 3.41 × 10−5; 4q26, rs3987, P = 0.10; 4q32.2, rs35509282, P = 0.24; 10q11.23, rs10994860, P = 3.65 × 10−4; 12q24.22, rs73208120, P = 0.03; 20q11.22, rs2295444, P = 0.03), all having a BFDP >0.99 (Supplementary Data 1). Of the 16 reported Asian-specific loci[4,5], nine harboured genome-wide significant signals in the current study (all BFDP <0.06), albeit sometimes at SNPs with low r2 but high D′ with the original SNP in Europeans, consistent with differences in allele frequencies in the different populations (Supplementary Data 1). Conditioning on the reported Asian SNPs, five of the nine European risk SNPs were independent of the Asian SNP (Pconditional < 5 × 10−8, Supplementary Data 4). We found no evidence of association signals at the remaining previously reported Asian SNPs.

Next, we performed an analysis conditioned on the sentinel SNP (r2 < 0.1 and Pconditional < 5 × 10−8; Table 2) to search for further independent signals at these new and previously reported risk loci. We confirmed the presence of previously reported dual signals at 14q22.2, 15q13.3 and 20p12.3[18]. For the new risk loci, an additional independent signal was identified at 5p15.3. In addition, a further seven signals were found at five previously reported risk loci: 11q13.4, 12p13.32, 15q13.3, 16q24.1, 20q13.13. Two of these signals were at the 15q13.3 locus, of which one was 5′ of GREM1 and the other intronic to FMN1. A further two signals were proximal and distal of 20q13.13. At 12p13.32 and 16q24.1, genome-wide associations marked by rs12818766 and rs899244, respectively, were shown. These were independent of the previously reported associations[2,14] at rs3217810 and rs2696839 (pairwise r2 = 0.0).

Table 2

Colorectal cancer variants identified in analysis conditioning on the sentinel SNP at each risk locus

Conditional (Sentinel) SNPs	Cytoband (position (bp, GRCh37))	Risk/Alt Allele	RAF	OR (95% CI)	P-value	Conditional OR (95% CI)	Conditional P-value	BFDP	LD with sentinel SNP (r2;D’)	I2 (%)	P het	Average info score
rs77776598 (rs2735940)	5p15.33 (1,240,998)	C/T	0.06	1.14 (1.09;1.20)	7.90 × 10−9	1.16 (1.11;1.21)	2.84 × 10−10	0.003	0.00; 0.33	0	0.93	0.99
rs4944940 (rs3824999)	11q13.4 (74,415,252)	G/A	0.96	1.31 (1.24;1.39)	1.05 × 10−20	1.28 (1.21;1.35)	3.21 × 10−17	2.73 × 10−9	0.00; 0.19	6	0.38	0.95
rs12818766 (rs3217810)	12p13.32 (4,376,091)	A/G	0.18	1.10 (1.07;1.13)	2.15 × 10−9	1.10 (1.07;1.13)	5.29 × 10−9	0.037	0.00; 0.06	30	0.16	0.89
rs1570405a (rs4444235)	14q22.2 (54,554,234)	G/A	0.31	1.06 (1.03;1.08)	9.81 × 10−7	1.07 (1.04;1.09)	1.91 × 10−8	0.125	0.02; 0.19	0	0.46	1.00
rs16969681b (rs73376930)	15q13.3 (32,993,111)	T/C	0.09	1.22 (1.18;1.27)	2.97 × 10−27	1.21 (1.16;1.25)	2.85 × 10−24	1.33 × 10−16	0.01; 0.32	42	0.04	0.99
rs16959063 (rs73376930)	15q13.3 (33,105,730)	A/G	0.01	1.30 (1.18;1.42)	3.72 × 10−8	1.33 (1.21;1.45)	5.40 × 10−9	0.23	0.00; 0.40	30	0.13	0.96
rs17816465 (rs73376930)	15q13.3 (33,156,386)	A/G	0.20	1.11 (1.08;1.14)	1.12 × 10−14	1.12 (1.09;1.15)	8.36 × 10−15	1.07 × 10−7	0.00; 0.11	44	0.04	0.97
rs899244 (rs2696839)	16q24.1 (86,700,030)	T/C	0.21	1.09 (1.06;1.12)	1.11 × 10−10	1.09 (1.06;1.12)	1.13 × 10−10	4.06 × 10−3	0.00; 0.04	14	0.29	0.99
rs6085661c (rs961253)	20p12.3 (6,693,128)	T/C	0.39	1.09 (1.06;1.11)	1.63 × 10−14	1.09 (1.07;1.11)	2.95 × 10−15	3.88 × 10−8	0.00; 0.08	0	0.96	1.00
rs4811050 (rs1810502)	20q13.13 (48,980,670)	A/G	0.18	1.10 (1.07;1.13)	2.43 × 10−11	1.09 (1.06;1.12)	4.07 × 10−9	4.06 × 10−3	0.04; 0.45	20	0.23	0.99
rs6091213 (rs1810502)	20q13.13 (49,384,745)	C/T	0.26	1.08 (1.05;1.11)	4.35 × 10−10	1.08 (1.05;1.11)	5.68 × 10−10	4.76 × 10−3	0.00; 0.05	6	0.39	0.94

BFDP calculated using prior = 10−5 and maximum relative risk = 1.2. LD calculated based on European populations in the 1000 Genomes Project data. BFDP calculated using conditional analysis results, with prior = 10−5 and maximum relative risk = 1.2

RAF risk allele frequency, OR odds ratio, CI confidence interval, BFDP Bayesian False Discovery Probability, I2 proportion of the total variation due to heterogeneity, PP-value for heterogeneity

aTags to rs1957636 (r2 = 0.67, D′ = 1). Previously identified in Tomlinson IP, Nat Genet, 2008 (PMID:18372905)

bPreviously identified in Tomlinson IP, Nat Genet, 2008 (PMID:18372905)

cTags to rs4813802 (r2 = 0.75, D′ = 0.93). Previously identified in Tomlinson IP, Nat Genet, 2008 (PMID:18372905)

In total, we identified 39 new independent SNPs associated with CRC susceptibility at genome-wide significance in Europeans. Together with the nine associations previously identified in Asian populations, and the 31 previously identified SNPs that were confirmed here, this brought the number of identified CRC association signals in Europeans to 79. Several of these risk loci map to regions previously identified in other cancers. In particular, three regions harbour susceptibility loci for multiple cancers[19], specifically 5p15.33 (TERT-CLPTM1L), 9p21.3 (CDKN2A) and 20q13.33 (RTEL1) (Supplementary Data 5).

Functional annotation and biological inference of risk loci

To the extent that they have been deciphered, most GWAS risk loci map to non-coding regions of the genome influencing gene regulation[19]. Consistent with this, we found evidence that the CRC risk SNPs mapped to regions enriched for active enhancer marks (H3K4me1 and H3K27ac) in colonic crypts (permutation test, P = 0.034 and 0.033, respectively) and colorectal tumours (P = 4.2 × 10−3 and 4.0 × 10−5) (Supplementary Figure 3). To determine whether the CRC SNPs overlapped with active regulatory regions in a cell-type specific manner[20], we analysed the H3K4me3, H3K27ac, H3K4me1, H3K27me3, H3K9ac, H3K9me3 and H3K36me3 chromatin marks across multiple cell types from the NIH Roadmap Epigenomics project[21]. Colonic and rectal mucosa cells showed the strongest enrichment of risk SNPs at active enhancer and promoter regions (H3K4me3, H3K4me1 and H3K27ac marks, P < 5 × 10−4) (Supplementary Figure 3).

Given our observation that the risk loci map to putative regulatory regions, we examined both histone modifications and transcription factor (TF) binding sites in LoVo and HT29 CRC cells across the risk SNPs. Using variant set enrichment[22], we identified regions of strong LD (defined as r2 > 0.8 and D′ > 0.8) with each risk SNP and determined the overlap with ChIP-seq data from the Systems Biology of Colorectal cancer (SYSCOL) study and inhouse-generated histone data. We identified an over-representation of binding for MYC, ETS2, cohesin loading factor NIPBL and cohesin-related proteins RAD21, SMC1A and SMC3 (Supplementary Figure 4). About 87% (69/79) of the risk SNPs were predicted to disrupt binding motifs of specific TFs, notably CTCF, SOX and FOX, with 35% located within TF binding peaks from LoVo, HT29 or ENCODE ChIP-seq data (Supplementary Data 6).

The upstream mechanisms by which predisposition SNPs influence disease risk is often through effects on cis-regulatory transcriptional networks, specifically through chromatin-looping interactions that are fundamental for regulation of gene expression. Therefore, to link regulatory regions containing risk SNPs to promoters of candidate target genes, we applied in situ promoter capture Hi-C (CHi-C) data in LoVo and HT29 cells (Supplementary Data 9). About 38% of the risk SNPs mapped to regions that showed statistically significant chromatin-looping interactions with the promoters of respective target genes. Notably, as well as confirming the interaction between rs6983267 and MYC at 8q24.21 (Supplementary Figure 2), the looping interaction from an active enhancer region at 10q25.2 implicates TCF7L2 as the target gene of rs12255141 variation (Fig. 3). TCF7L2 (previously known as TCF4) is a key transcription factor in the Wnt pathway and plays an important role in the development and progression of CRC[23]. Intriguingly, TCF7L2 has been shown to bind to a MYC enhancer containing rs6983267[24] and to a GREM1 enhancer near rs16969681[25]. Based on ChromHMM, this region is annotated as a promoter in HCT116 cells, but not in normal colonic and rectal mucosa. Additionally this locus has been implicated in lung cancer[26] and low-grade glioma[27]. Similarly, the 9p21.3 chromatin interaction provides evidence to support CDKN2B as the target gene for rs1412834 variation, a region of somatic loss.

Fig. 3

Regional plots of exemplar colorectal cancer risk loci. In the main panel, −log10P-values (y-axis) of the SNPs are shown according to their chromosomal positions (x-axis). The colour intensity of each symbol reflects the extent of LD with the top SNP: white (r2 = 0) through to dark red (r2 = 1.0), with r2 estimated from EUR 1000 Genomes data. Genetic recombination rates (cM/Mb) are shown with a light blue line. Physical positions are based on GRCh37 of the human genome. Where available, the upper panel shows Hi-C contacts from HT29 or LoVo. The lower panel shows the chromatin state segmentation track from the Roadmap Epigenomics project (colonic mucosa, rectal mucosa, sigmoid colon), and HCT116. Also shown are the relative positions of genes and transcripts mapping to each region of association. a rs12255141 (10q25.2); b rs12979278 (19q13); c rs2735940 (5p15); d rs7398375 (12q13.3)

We sought to gain further insight into the target genes at each locus, and hence the biological mechanisms for the associations, by performing expression quantitative trait locus (eQTL) analysis in colorectal tissue. We analysed inhouse eQTL data generated from samples of normal colonic mucosa (INTERMPHEN study, n = 131 patients) and GTEx data from transverse colon (n = 246). For the previously identified risk loci, there were eQTLs for rs4546885 and LAMC1 (1q25.3), rs13020391 and lnc-RNA RP11–378A13.1 (2q35), and rs3087967 and COLCA1, COLCA2 and C11orf53 (11q23.1). Amongst the eQTL associations at the new risk loci, pre-eminent eQTLs were rs9831861 and SFMBT1 (3p21.1), rs12427600 and SMAD9 (13q13.3), and rs12979278 and FUT2 and MAMSTR (19q13.33) (Supplementary Data 7). However, while multiple nominally significant cis-eQTLs were present, nearly half of all loci had no evidence of cis-eQTLs in the sample sets used.

In addition to eQTL analysis, we performed Summary-data-based Mendelian Randomization (SMR) analysis[28] as a more stringent test for causal differences in gene transcription (Supplementary Data 8). There was support for the 11q23.1 locus SNP influencing CRC risk through differential expression of one or more of COLCA1, COLCA2 and C11orf53 transcripts (PSMR < 10−10). There was also evidence that the 3p21.1 and 19q13.33 SNPs acted through SFMBT1 and FUT2, respectively, (PSMR < 10−5), as well as the 6p21.31 SNP acted through class II HLA expression (PSMR < 5 × 10−4).

Based on genetic fine-mapping and functional annotation, our data indicated several candidate target genes with functions previously unconnected to colorectal tumourigenesis (Supplementary Data 9). The SFMBT1 protein (3p21.1) acts as a histone reader and a component of a transcriptional repressor complex[29]. TNS3 at 7p12.3 encodes the focal adhesion protein TENSIN3, to which the intestinal stem cell marker protein Musashi1 has been reported to bind. Tns3-null mice exhibit impaired intestinal epithelial development, probably because of defects in Rho GTPase signalling and cell adhesion[30]. LRP1 (12q13.3, LDL receptor-related protein 1) (Fig. 3) may be involved in Wnt-signalling[31], although its role in the intestines has not previously been conclusively demonstrated. FUT2 at 19q13.33 encodes fucosyltransferase II. Variation at this locus is associated with differential interactions with intestinal bacteria and viruses. Our data thus provide evidence for a role of the microbiome in CRC risk[32]. PTPN1 (20q13.13), also known as PTP1B, encodes a non-receptor tyrosine phosphatase involved in regulating JAK-signalling, IR, c-Src, CTNNB1, and EGFR.

We annotated all risk loci with five types of functional data: (i) presence of a CHi-C contact linking to a gene promoter, (ii) presence of an association from eQTL, (iii) presence of a regulatory state, (iv) evidence of TF binding, and (v) presence of a nonsynonymous coding change (Supplementary Data 9). Collectively this analysis suggested three primary candidate disease mechanisms across a number of risk loci: firstly, genes linked to BMP/TGF-β signalling (e.g. GREM1, BMP2, BMP4, SMAD7, SMAD9); secondly, genes with roles either directly or indirectly linked to MYC (e.g. MYC, TCF7L2); and thirdly genes with roles in maintenance of chromosome integrity (e.g. TERT, RTEL1) and DNA repair (e.g. POLD3) (Supplementary Figure 5).

Pathway gene set enrichment analyses[33] revealed several growth or development related pathways were enriched, notably TGF-β signalling and immune response pathways (Supplementary Figure 6, Supplementary Data 10). Other cancer-related themes included apoptosis and leukocyte differentiation pathways. We used Data-driven Expression-Prioritized Integration for Complex Traits (DEPICT)[34] to predict gene targets based on gene functions that are shared across genome-wide significant risk loci, as well as those associated at P < 10−5 as advocated to mitigate against type II error. Tissue-specificity with respect to colonic tissue was evident (permutation test, P < 5 × 10−3) and among the protein-coding genes predicted, there was enrichment for TGF-β and PI3K-signalling pathways, and abnormal intestinal crypt gene sub-networks (P < 10−5; Supplementary Data 11).

Contribution of risk SNPs to heritability

Using Linkage Disequilibrium Adjusted Kinships (LDAK)[35] in conjunction with the GWAS data generated on unselected CRC cases (i.e. COIN, CORSA, Croatia, DACHS, FIN, SCOT, Scotland1, SOCCS/LBC, SOCCS/GS, UKBB, VQ58 studies) we estimated that the heritability of CRC attributable to all common variation is 0.29 (95% confidence interval: 0.24–0.35). To estimate the sample size required to explain a greater proportion of the GWAS heritability, we implemented a likelihood-based approach using association statistics in combination with LD information to model the effect-size distribution[36], which was best represented by a three-component model (mixture of two normal distributions). Under this model, to identify SNPs explaining 80% of the GWAS heritability, it is likely to require effective sample sizes in excess of 300,000 if solely based on GWAS associations (Supplementary Figure 7).

After adjusting for winner’s curse[37], the 79 SNPs thus far shown to be associated with CRC susceptibility in Europeans explain 11% of the 2.2-fold familial relative risk (FRR)[38], whilst all common genetic variants identifiable through GWAS could explain 73% of the FRR. Thus, the identified susceptibility SNPs collectively account for approximately 15% of the FRR of CRC that can be explained by common genetic variation. We incorporated the newly identified SNPs into risk prediction models for CRC and derived a polygenic risk score (PRS) based on a total of 79 GWAS significant risk variants. Individuals in the top 1% have a 2.6-fold increased risk of CRC compared with the population average (Supplementary Figure 8). Risk re-classification using this PRS offers the prospect of optimising prevention programmes for CRC in the population, for example through targeting screening[6], and also preventative interventions. The identification of further risk loci through the analysis of even larger GWAS is likely to improve the performance of any PRS model.

Co-heritability with non-cancer traits

We implemented cross-trait LD score regression[39] to investigate co-heritability globally between CRC and 41 traits with publicly available GWAS summary statistics data. None of the genetic correlations remained significant after Bonferroni correction (two-sided Z-test, P-threshold: 0.05/41 = 1.2 × 10−3). However, nominally significant positive associations with CRC risk (Supplementary Data 12) included insulin resistance, comprising raised fasting insulin, glucose and HbA1c (positive), hyperlipidaemia, comprising raised total cholesterol and low-density lipoprotein cholesterol, and ulcerative colitis, all of which are traits or diseases previously reported in observational epidemiological studies to be associated with CRC risk[40,41].

Discussion

Here we report a comprehensive analysis that sheds new light on the molecular basis of genetic risk for a common cancer, and greatly increases the number of known CRC risk SNPs. To identify the most credible target genes at each site, we have performed detailed annotation using public databases, and have also acquired our own disease-specific data from ChIP-seq, promoter capture Hi-C and gene expression analyses.

Given that there remains significant missing common heritability for CRC, additional GWAS meta-analyses are likely to lead to discovery of more risk loci. Such an assertion is directly supported a contemporaneous study[42], which has reported the identification of 40 independent signals; 30 novel loci and 10 conditionally independent association signals at previously and newly identified CRC risk loci. Of these, 18 were replicated in our analysis, with an additional five exhibiting an independent signal present at the same locus (Supplementary Data 13).

Overall, our findings provide new insights into the biological basis of CRC, not only confirming the importance of established gene networks, but also providing evidence that point to a role for the gut microbiome in CRC causation, and identifying several functional mechanisms previously unsuspected of any involvement in colorectal tumourigenesis. Several of the gene pathways identified through GWAS may provide potential novel targets for chemoprevention and chemotherapeutic intervention.

Methods

Ethics

Collection of patient samples and associated clinico-pathological information was undertaken with written informed consent and relevant ethical review board approval at respective study centres in accordance with the tenets of the Declaration of Helsinki. Specifically: (i) UK National Cancer Research Network Multi-Research Ethics Committee (02/0/097 [NSCCG], 01/0/5) [SOCCS], 05/S1401/89 [GS:SFHS], LREC/1998/4/183 [LBC1921], 2003/2/29 [LBC1936], 17/SC/0079 [CORGI] and 07/S0703/136 [SCOT]); (ii) The research activities of UK Biobank were approved by the North West Multi-centre Research Ethics Committee (11/NW/0382) in relation to the process of participant invitation, assessment and follow-up procedures. Additionally, ethics approvals from the National Information Governance Board for Health & Social Care in England and Wales and approval from the Community Health Index Advisory Group in Scotland were also obtained to gain access to the information that would allow the invitation of participants. This study did not need to re-contact the participants, and no separate ethics approval was required according to the Ethics and Governance Framework (EGF) of UK Biobank; (iii) South East Ethics Committee MREC (03/1/014); (iv) Written informed consent was obtained from all participants of CORSA. The study was approved by the ethical review committee of the Medical University of Vienna (MUW, EK Nr. 703/2010) and the “Ethikkommission Burgenland” (KRAGES, 33/2010) and (v) Finnish National Supervisory Authority for Welfare and Health, National Institute for Health and Welfare (THL/151/5.05.00/2017), the Ethics Committee of the Hospital District of Helsinki and Uusimaa (HUS/408/13/03/03/09).

The diagnosis of colorectal cancer (ICD-9 153, 154; ICD-10 C18.9, C19, C20) was established in all cases in accordance with World Health Organization guidelines.

Primary GWAS

We analysed data from five primary GWAS (Supplementary Data 2 and Supplementary Data 3):

The NSCCG-OncoArray GWAS comprised 6240 cases ascertained through the National Study of Colorectal Cancer Genetics (NSCCG)[43] and 1041 cases collected through the CORGI consortium, genotyped using the Illumina OncoArray. Patients were selected for having a family history of CRC (at least one first-degree relative) or age of diagnosis below 58. Controls were also genotyped using the OncoArray and comprised (i) 3031 cancer-free men recruited by the PRACTICAL Consortium—the UK Genetic Prostate Cancer Study (UKGPCS) (age <65 years), a study conducted through the Royal Marsden NHS Foundation Trust and SEARCH (Study of Epidemiology & Risk Factors in Cancer), recruited via GP practices in East Anglia (2003–2009) and (ii) 4,488 cancer-free women across the UK, recruited via the Breast Cancer Association Consortium (BCAC).

The SCOT GWAS comprised 3076 cases from the Short Course Oncology Treatment (SCOT) trial—a study of adjuvant chemotherapy in colorectal cancer by the CACTUS and OCTO groups[44]. Controls comprised 4349 cancer-free individuals from The Heinz Nixdorf Recall study[45]. Both cases and controls were genotyped using the Illumina Global Screening Array.

SOCCS/Generation Scotland (SOCCS/GS) comprised 4772 cases from the Study of Colorectal Cancer in Scotland (SOCCS)[12,13] and 12,158 controls including 2221 population-based controls from SOCCS and additional 9937 population controls without prior history of colorectal cancer from Generation Scotland-Scottish Family Health Study (GS:SFHS)[46].

SOCCS/Lothian Birth Cohort (SOCCS/LBC) GWAS comprised 1037 cases from the Study of Colorectal Cancer in Scotland (SOCCS)[47] and 1522 population-based controls without prior history of malignant tumours from the Lothian Birth Cohorts (LBC) of 1921 and 1936[48].

UK Biobank (UKBB) GWAS comprised 6360 cases and 25,440 population-based control individuals. UK Biobank is a large cohort study with more than 500,000 individuals recruited. Biological samples of these participants were genotyped using the custom-designed Affymetrix UK BiLEVE Axiom array on an initial 50,000 participants and Affymetrix UK Biobank Axiom array on the remaining 450,000 participants. The two arrays had over 95% common content. Genotyping was done at the Affymetrix Research Services Laboratory in Santa Clara, California, USA. Details on genotyping and quality control were previously reported[49]. Self-reported cases of cancers of bowel, colon or rectum, if not confirmed by the ICD9 or ICD10 codes were excluded from the analysis. Healthy control individuals without history of cancer and/or colorectal adenoma were included in the analysis after matching one case to four controls by age, gender, date of blood draw, ethnicity and region of residence (two first letters of postal code).

Published GWAS

We made use of 10 previously published GWAS (Supplementary Data 2): (1) UK1 (CORGI study) comprised 940 cases with colorectal neoplasia and 965 controls[12]; (2) Scotland1 (COGS study) included 1012 CRC cases and 1012 controls[12]; (3) VQ58 comprised 1800 cases from the UK-based VICTOR and QUASAR2 adjuvant chemotherapy clinical trials and 2690 population control genotypes from the Wellcome Trust Case Control Consortium 2 (WTCCC2) 1958 birth cohort[50]; (4) CCFR1 comprised 1290 familial CRC cases and 1055 controls from the Colon Cancer Family Registry (CCFR)[15]; (5) CCFR2 included a further 796 cases from the CCFR and 2236 controls from the Cancer Genetic Markers of Susceptibility (CGEMS) studies of breast and prostate cancer[51,52]; (6) COIN was based on 2244 CRC cases ascertained through two independent Medical Research Council clinical trials of advanced/metastatic CRC (COIN and COIN-B)[53] and controls comprised 2162 individuals from the UK Blood Service Control Group genotyped as part of the WTCCC2; (7) Finnish GWAS (FIN)[3] was based on 1172 CRC cases and 8266 cancer-free controls ascertained through FINRISK, Health 2000, Finnish Twin Cohort and Helsinki Birth Cohort Studies; (8) CORSA (COloRectal cancer Study of Austria) a molecular epidemiological study of 978 cases and 855 colonoscopy-negative controls[54]; (9) DACHS (Darmkrebs: Chancen der Verhütung durch Screening)[55] based on 1105 cases and 700 controls and (10) Croatia consisted of 764 cases and 460 population-based controls[56].

The VQ58, UK1 and Scotland1 GWAS were genotyped using Illumina Hap300, Hap240S, Hap370, Hap550 or Omni2.5 M arrays. 1958BC genotyping was performed as part of the WTCCC2 study on Hap1.2M-Duo Custom arrays. The CCFR samples were genotyped using Illumina Hap1M, Hap1M-Duo or Omni-express arrays. CGEMS samples were genotyped using Illumina Hap300 and Hap240 or Hap550 arrays. The COIN cases were genotyped using Affymetrix Axiom Arrays and the Blood Service controls were genotyped using Affymetrix 6.0 arrays. FIN cases were genotyped using Illumina HumanOmni 2.5M8v1 and controls using Illumina HumanHap 670k and 610k arrays. DACHS study samples were genotyped using the Illumina OncoArray, CORSA study sampels were genotyped on the Affymetrix Axiom Genome-Wide CEU 1 Array, and Croatia study samples were genotyped on Illumina OmniExpressExome BeadChip 8v1.1 or 8v1.3.

Quality control

Standard quality control (QC) measures were applied to each GWAS[8]. Specifically, individuals with low SNP call rate (<95%) as well as individuals evaluated to be of non-European ancestry (using the HapMap version 2 CEU, JPT/CHB and YRI populations as a reference) were excluded (Supplementary Figure 9). For apparent first-degree relative pairs, we excluded the control from a case-control pair; otherwise, we excluded the individual with the lower call rate. SNPs with a call rate <95% were excluded as were those with a MAF <0.5% or displaying significant deviation from Hardy–Weinberg equilibrium (P < 10−5). QC details are provided in Supplementary Data 3. All genotype analyses were performed using PLINK v1.9[57].

Imputation and statistical analysis

Prediction of the untyped SNPs was carried out using SHAPEIT v2.837[58] and IMPUTE v2.3.2[59]. The CCFR1, CCFR2, COIN, CORSA, Croatia, NSCCG-OncoArray, SCOT, Scotland1, SOCCS/GS, SOCCS/LBC, UK1 and VQ58 samples used a merged reference panel using data from 1000 Genomes Project (phase 1, December 2013 release) and UK10K (April 2014 release). Imputation of UKBB was based on data from 1000 Genomes Project (phase 3), UK10K and Haplotype Reference Consortium. The FIN and DACHS GWAS were imputed using a reference panel comprised of 1000 Genomes Projects Project with an additional population matched reference panel: 3882 Sequencing Initiative Suomi (SISu) haplotypes for the FIN study, and 3000 sequenced CRC cases for the DACHS study. We imposed predefined thresholds for imputation quality to retain potential risk variants with MAF >0.5% for validation. Poorly imputed SNPs defined by an information measure <0.80 were excluded. Tests of association between imputed SNPs and CRC were performed under an additive genetic model in SNPTEST v2.5.2[60]. Principal components were added to adjust for population stratification where required (i.e. DACHS, FIN, NSCCG-OncoArray, SCOT and UKBB).

To determine whether specific coding variants within HLA genes contributed to the diverse association signals, we imputed the classical HLA alleles (A, B, C, DQA1, DQB1 and DRB1) and coding variants across the HLA region using SNP2HLA[11]. The imputation was based on a reference panel from the Type 1 Diabetes Genetics Consortium (T1DGC) consisting of genotype data from 5225 individuals of European descent with genotyping data of 8961 common SNPs and indel polymorphisms across the HLA region, and four digit genotyping data of the HLA class I and II molecules. For the X chromosome, genotypes were phased and imputed as for the autosomal chromosome, with the inclusion of the “chrX” flag. X chromosome association analysis was performed in SNPTEST using a maximum likelihood model, assuming complete inactivation of one allele in females and equal effect-size between males and females.

The adequacy of the case-control matching and possibility of differential genotyping of cases and controls was evaluated using a Q–Q plot of test statistics in individual studies (Supplementary Figure 1). Meta-analyses were performed using the fixed-effects inverse-variance method using META v1.7[61]. Cochran’s Q-statistic to test for heterogeneity and the I statistic to quantify the proportion of the total variation due to heterogeneity were calculated. A Q–Q plot of the meta-analysis test statistics was also performed (Supplementary Figure 1). None of the studies showed evidence of genomic inflation, where λGC values for the CCFR1, CCFR2, COIN, CORSA, Croatia, DACHS, FIN, NSCCG-OncoArray, SCOT, Scotland1, SOCCS/GS, SOCCS/LBC, UKBB, UK1 and VQ58 studies were 1.03, 1.08, 1.09, 1.11, 1.01, 1.01, 1.09, 1.10, 1.08, 1.02, 1.09, 1.04, 1.05, 1.02 and 1.06, respectively. Estimates were calculated using the regression method, as implemented in GenABEL.

Definition of known and new risk loci

We sought to identify all associations for CRC previously reported at a significance level P < 5 × 10−8 by referencing the NHGRI-EBI Catalog of published genome-wide association studies, and a literature search for the years 1998–2018 using PubMed (performed January 2018). Additional articles were ascertained through references cited in primary publications. Where multiple studies reported associations in the same region, we only considered the first reported genome-wide significant association. New loci were identified based on SNPs at P < 5 × 10−8 using the meta-analysis summary statistics, with LD correlations from a reference panel of the European 1000 Genomes Project samples combined with UK10K. We only included one SNP per 500 kb interval. To measure the probability of the hits being false positives, the Bayesian False-Discovery Probability (BFDP)[9] was calculated based on a plausible OR of 1.2 (based on the 95th percentile of the meta-analysis OR values) and a prior probability of association of 10−5. A conditional analysis was performed using Genome-wide Complex Trait Analysis (GCTA)[62], conditioning on the new and known SNPs, and SNPs with Pconditioned < 5 × 10−8 and r2 > 0.1 were clumped using PLINK. The NSCCG-Oncoarray data were used to provide the LD reference data.

Fidelity of imputation

The reliability of imputation of the novel risk SNPs identified (all with an IMPUTE2 r2 > 0.8) was assessed for 51 SNPs (comprising all new signals not directly genotyped) by examining the concordance between imputed and whole-genome sequenced genotypes in a subset of 201 samples from the CORGI and NSCCG studies. More than 98% concordance was found between the directly sequenced and imputed SNPs (Supplementary Data 14).

eQTL analysis

In the INTERMPHEN study, biopsies of normal colorectal mucosa (trios of rectum, proximal colon and distal colon) were obtained from 131 UK individuals with self-reported European ancestry without CRC. Genotyping was performed using the Illumina Infinium Human Core Exome array, with quality control and imputation as above. RNA-seq was performed and data analysed as per the GTEx Project pipeline v7 using the 1000 Genomes and UK10K data as reference. Gene-level expression quantification was based on the GENCODE 19 annotation, collapsed to a single transcript model for each gene using a custom isoform procedure. Gene-level quantification (read counts and TPM values) was performed with RNA-SeQC v1.1.8. Gene expression was normalised using the TMM algorithm, implemented in edgeR, with inverse normal transformation, based on gene expression thresholds of >0.1 Transcripts Per Million (TPM) in ≥20% of samples and ≥6 reads in ≥20% of samples. cis-eQTL mapping was performed separately for proximal colon, distal colon and rectum samples using FastQTL. Principal components for the SNP data and additional covariate factors were identified using Probabilistic Estimation of Expression Residuals (PEER). P-values were generated for each variant-gene pair testing alternative hypothesis that the slope of a linear regression model between genotype and expression deviates from 0. The mapping window was defined as 1 Mb either side of the transcription start site. Beta distribution-adjusted empirical P-values from FastQTL were used to calculate Q-values, and FDR threshold of ≤0.05 was applied to identify genes with a significant eQTL. The normalised effect size of the eQTLs was defined as the slope of the linear regression, and computed as the effect of the alternative allele relative to the reference allele in the human genome reference GRCh37/hg19). MetaTissue was used to generate a “pan-colonic” eQTL measure from the three individual RNA-seq datasets per patient.

To supplement this analysis, we performed SMR analysis[28] including all eQTLs with nominally significant associations (P < 0.05). We additionally examined for heterogeneity using the heterogeneity in dependent instruments (HEIDI) test, where PHEIDI < 0.05 were considered as reflective of heterogeneity and were excluded.

Promoter capture Hi-C

In situ promoter capture Hi-C (CHi-C) on LoVo and HT29 cell lines was performed as previously described[63]. Hi-C and CHi-C libraries were sequenced using HiSeq 2000 (Illumina). Reads were aligned to the GRCh37 build using bowtie2 v2.2.6 and identification of valid di-tags was performed using HiCUP v0.5.9. To declare significant contacts, HiCUP output was processed using CHiCAGO v1.1.8. For each cell line, data from three independent biological replicates were combined to obtain a definitive set of contacts. As advocated, interactions with a score ≥5.0 were considered to be statistically significant[64].

Chromatin state annotation

Colorectal cancer risk loci and SNPs in LD (r2 > 0.8) were annotated for putative functional effect based upon ChIP-seq H3K4me1 (C15410194), H3K9me3 (C15410193), H3K27me3 (C15410195) and H3K36me3 (C15410192) for LoVo, and H3K4me1 and H3K9me3 for HT29. ChIP libraries were sequenced using HiSeq 2000 (Illumina) with 100 bp single-ended reads. Generated raw reads were filtered for quality (Phred33 ≥ 30) and length (n ≥ 32), and adapter sequences were removed using Trimmomatic v0.22. Reads passing filters were then aligned to the human reference (hg19) using BWA v0.6.1. Peak calls are obtained using MACS2 v 2.0.10.07132012.

Histone mark and transcription factor enrichment analysis

ChIP-seq data from colon crypt and tumour samples was obtained for H3K27ac and H3K4me1[65]. Multiple samples of the same tissue type or tumour stage were merged together. Additional ChIP-seq data from the Roadmap Epigenomics project[21] was obtained for H3K4me3, H3K27ac, H3K4me1, H3K27me3, H3K9ac, H3K9me3 and H3K36me3 marks in up to 114 tissues. Overlap enrichment analysis of CRC risk SNPs with these peaks was performed using EPIGWAS, as described by Trynka et al.[20]. Briefly, we evaluated if CRC risk SNPs and SNPs in LD (r2 > 0.8) with the sentinel SNP, were enriched at ChIP-seq peaks in tissues by a permutation procedure with 105 iterations.

To examine enrichment in specific TF binding across risk loci, we adapted the variant set enrichment method of Cowper-Sal lari et al.[22]. Briefly, for each risk locus, a region of strong LD (defined as r2 > 0.8 and D′ > 0.8) was determined, and these SNPs were termed the associated variant set (AVS). ChIP-seq uniform peak data were obtained for LoVo and HT29 cell lines (198 and 29 experiments, respectively)[66] and the above described histone marks. For each of these marks, the overlap of the SNPs in the AVS and the binding sites was determined to produce a mapping tally. A null distribution was produced by randomly selecting SNPs with the same characteristics as the risk-associated SNPs, and the null mapping tally calculated. This process was repeated 105 times, and P-values calculated as the proportion of permutations where the null mapping tally was greater or equal to the AVS mapping tally. An enrichment score was calculated by normalising the tallies to the median of the null distribution. Thus, the enrichment score is the number of standard deviations of the AVS mapping tally from the median of the null distribution tallies.

Functional annotation

For the integrated functional annotation of risk loci, LD blocks were defined as all SNPs in r2 > 0.8 with the sentinel SNP. Risk loci were then annotated with five types of functional data: (i) presence of a CHi-C contact linking to a gene promoter, (ii) presence of an association from eQTL, (iii) presence of a regulatory state, (iv) evidence of TF binding, and (v) presence of a nonsynonymous coding change. Candidate causal genes were then assigned to CRC risk loci using the target genes implicated in annotation tracks (i), (ii), (iiii) and (iv). If the data supported multiple gene candidates, the gene with the highest number of individual functional data points was considered as the candidate. Where multiple genes had the same number of data points, all genes were listed. Direct nonsynonymous coding variants were allocated additional weighting. Competing mechanisms for the same gene (e.g. both coding and promoter variants) were allowed for. Finally, if no evidence was provided by these criteria, if the lead SNP was intronic we assigned candidacy on this basis, or if intergenic the nearest gene neighbour. Chromatin data were obtained from HaploReg v4 and regulatory regions from Ensembl.

Regional plots were created using visPIG[67], using the data described above. We used ChromHMM to integrate DNAse, H3K4me3, H3K4me1, H3K27ac, Pol2 and CTCF states from the CRC cell line HCT116 using a multivariate Hidden Markov Model[68]. Chromatin annotation tracks for colonic mucosa (E075), rectal mucosa (E101) and sigmoid colon (E106) were obtained from the Roadmap Epigenomics project[21], using the core 15-state model data based on H3K4me3, H3K4me1, H3K36me3, H3K27me3 and H3K9me3 marks.

Transcription factor binding disruption analysis

To determine if the risk variants or their proxies were disrupting motif binding sites, we used the motifbreakR package[69]. This tool predicts the effects of variants on TF binding motifs, using position probability matrices to determine the likelihood of observing a particular nucleotide at a specific position within a TF binding site. We tested the SNPs by estimating their effects on over 2,800 binding motifs as characterised by ENCODE, FactorBook, HOCOMOCO and HOMER. Scores were calculated using the relative entropy algorithm.

Heritability analysis

We used LDAK[35] to estimate the polygenic variance (i.e. heritability) ascribable to SNPs from summary statistic data for the GWAS datasets which were based on unselected cases (i.e. CORSA, COIN, Croatia, DACHS, FIN, SCOT, Scotland1, SOCCS/GS, SOCCS/LBC, UKBB and VQ58). SNP-specific expected heritability, adjusted for LD, MAF and genotype certainty, was calculated from the UK10K and 1000 Genomes data. Individuals were excluded if they were closely related, had divergent ancestry from CEU, or had a call rate <0.99. SNPs were excluded if they showed deviation from HWE with P < 1 × 10−5, genotype yield <95%, MAF <1%, SNP imputation score <0.99, and the absence of the SNP in the GWAS summary statistic data. This resulted in a total 6,024,731 SNPs used to estimate the heritability of CRC.

To estimate the sample size required to detect a given proportion of the GWAS heritability we implemented a likelihood-based approach to model the effect-size distribution[36], using association statistics from the meta-analysis, and LD information from individuals of European ancestry in the 1000 Genomes Project Phase 3. LD values were based on an r2 threshold of 0.1 and a window size of 1MB. The goodness of fit of the observed distribution of P-values against the expected from a two-component model (single normal distribution) and a three-component model (mixture of two normal distributions) were assessed, and a better fit was observed for the latter model. The percentage of GWAS heritability explained for a projected sample size was determined using this model, based on power calculations for the discovery of genome-wide significant SNPs. The genetic variance explained was calculated as the proportion of total GWAS heritability explained by SNPs reaching genome-wide significance at a given sample size. The 95% confidence intervals were determined using 105 simulations.

Cross-trait genetic correlation

LD score regression[39] was used to determine if any traits were correlated with CRC risk. GWAS summary data were obtained for allergy, asthma, coronary artery disease, fatty acids, lipids (total cholesterol, high density lipoprotein, low-density lipoprotein, triglycerides), auto-immune diseases (Crohn’s disease, rheumatoid arthritis, atopic dermatitis, celiac disease, multiple sclerosis, primary biliary cirrhosis, inflammatory bowel disease, ulcerative colitis, systemic lupus erythematosus), anthropometric measures (BMI, height, body fat), glucose sensitivity (fasting glucose, fasting insulin, HbA1c), childhood measures (birth weight, birth length, childhood obesity, childhood BMI), eGFR and type 2 diabetes. All data were obtained for European populations. Summary statistics were reformatted to be consistent, and constrained to HapMap3 SNPs as these have been found to generally impute well. LD Scores were determined using 1000 Genomes European data.

Familial risk explained by risk SNPs

Under a multiplicative model, the contribution of risk SNPs to the familial risk of CRC was calculated from , where λ0 is the familial risk to first-degree relatives of CRC cases, assumed to be 2.2[38], and λ is the familial relative risk associated with SNP k, calculated as , where pk is the risk allele frequency for SNP k, qk  1−pk, and rk is the estimated per-allele OR from the meta-analysis[70]. The OR estimates were adjusted for the winner’s curse using the FDR Inverse Quantile Transformation (FIQT) method[37]. We constructed a PRS including all 79 CRC risk SNPs discovered or validated by this GWAS in the risk-score modelling. The distribution of risk on an RR scale in the population is assumed to be log-normal with arbitrary population mean μ set to -σ2/2 and variance where β and p correspond to the log odds ratio and the risk allele frequency, respectively, for SNP k. The distribution of PRS among cases is right-shifted by σ2 so that the overall mean PRS is 1.0[71]. The risk distribution was also performed assuming all common variation, using , where λsib = 1.79, as determined using the heritability estimate from GCTA.

Pathway analysis

SNPs were assigned to genes as described in the functional annotation section. The genes that mapped to genome-wide significant CRC risk SNPs were analysed using InBio Map, a manually curated database of protein-protein interactions.

Gene set enrichment was calculated using GenGen. Enrichment scores were calculated using the meta-analysis results and were based on 103 permutations on the χ2 values between SNPs. Pathway definitions were obtained from the Bader Lab[33], University of Toronto, July 2018 release. This data contained pathway information from Gene Ontology (GO), Reactome, HumanCyc, MSigdb C2 (curated dataset), NCI Pathway, NetPath and PANTHER for a total of 7269 pathways. GO annotations that were inferred computationally were excluded. To avoid biasing the results, the meta-analysis SNPs were pruned to only those with an r2 < 0.1 and a distance greater than 500 kb. Pathways were visualised using Cytoscape v3.6.1, together with the EnrichmentMap v3.1.0 and AutoAnnotate v1.2 plugins. Only pathways with an FDR <0.05 and edges with a similarity coefficient (number of shared genes between pathways) >0.55 were displayed.

URLs

Bader Lab pathway data: http://download.baderlab.org/EM_Genesets/July_01_2018/Human/symbol/

FastQTL: https://github.com/francois-a/fastqtl

GTEx: https://www.gtexportal.org/home/

InBioMap: https://www.intomics.com/inbio/map/#home

LD scores: https://data.broadinstitute.org/alkesgroup/LDSCORE/

NHGRI-EBI GWAS Catalog: https://www.ebi.ac.uk/gwas/

PredictDB: http://predictdb.org/

Roadmap Epigenomics data: https://egg2.wustl.edu/roadmap/web_portal/chr_state_learning.html

SYSCOL: http://syscol-project.eu/

UK Biobank: http://www.ukbiobank.ac.uk/scientists-3/genetic-data/