Cell-type-resolved genetic variation shapes inflammatory bowel disease risk

Patient recruitment and single-cell sequencing

This study was approved by the National Health Service (NHS) Research Ethics Committee (Cambridge South, REC ID 17/EE/0338). Written informed consent was given by all participants.

296 individuals without IBD and 125 individuals with CD were recruited at Addenbrooke’s Hospital, Cambridge, UK. Patients with CD were required to have a history of ileal disease involvement. During routine endoscopy, gastrointestinal biopsies were collected from the rectum and terminal ileum, and blood samples were drawn after the procedure. Terminal ileal biopsies were obtained from 243 individuals without IBD and 119 individuals with CD. Ileal biopsies from patients with CD included both inflamed and uninflamed tissue. Rectal biopsies were collected from 275 individuals without IBD. Where possible, rectal and ileal samples were collected from the same healthy individuals, and blood and ileal samples were collected from the same patients. For gastrointestinal biopsies, we have developed a gentle single-cell dissociation protocol that minimizes transcriptional changes and cell stress or death during cell isolation. Early samples were processed using version 1 of the protocol, whereas most of the samples were processed using version 2^90,91. Whole blood was collected for scRNA-seq from 95 individuals with CD. Peripheral blood mononuclear blood cells (PBMCs) were isolated from the whole blood using the EasySep Direct Human PBMC Isolation Kit (StemCell Technologies).

Because IBD tissue displays increased cellular heterogeneity due to immune infiltration, we targeted 6,000 and 3,000 cells for CD and healthy biopsies, respectively. This is to ensure adequate representation of diverse cell populations and to maintain statistical power for expression quantification and eQTL mapping across annotations. scRNA-seq was performed using 3′ 10X Genomics Chromium kits (v3.0 and v3.1) according to the manufacturer’s instructions. Libraries were sequenced to a target depth of 50,000 reads per cell. CellRanger (v7.2.0) was used to align reads to the GRCh38 human genome reference with Ensembl v93 transcript definitions (GRCh38-3.0.0 reference file provided by 10X Genomics), and to generate cell-by-gene count matrices. The greater number of cells targeted for patients with CD is confirmed in the terminal ileum samples from the ‘atlasing’ dataset (see ‘Single-cell RNA sequencing quality control and cell clustering’) (Wilcoxon P value < 2.2 × 10⁻¹⁶), with no discernable influence on sequencing depth, as indicated by median number of counts per sample (P = 0.27) (Supplementary Fig. 6).

Single-cell RNA-sequencing quality control and cell clustering

Cell count matrices were initially subset to retain droplets likely to contain cells and corrected for ambient RNA using CellBender v2.1⁹², as previously described⁷⁶. Scrublet v0.2.1⁹³ was then used to identify and remove potential multiplets. All subsequent quality control, normalization, and clustering steps were performed with Scanpy v1.9.6⁹⁴, including normalizing to library sizes of 10,000 total counts using the function sc.pp.normalise_per_cell and log transformation using sc.pp.log1p to produce the ln[cp10k + 1] count matrix.

A lineage-sensitive, relative quality control strategy was applied to account for technical and biological variation across tissues and lineages (Supplementary Fig. 1). First, the input dataset (4,099,804 cells) was aggregated and subjected to an initial round of quality control. Cells were removed if they had fewer than 250 genes expressed, fewer than 500 total counts, or more than 20% or 50% of reads originating from the mitochondrial genome, for blood-derived and non-blood-derived cells, respectively. This left a total of 2,725,568 cells. Genes were then removed if they were not expressed in at least five cells. A total of 5,000 highly variable genes were then classified within each sample and merged using sc.pp.highly_variable_genes and flavour = seurat. Ribosomal, mitochondrial and immunoglobulin genes were then excluded from the highly variable set.

These data were then subject to an initial round of integration and clustering. Data were integrated using the scVI v0.16.3⁹⁵ function scvi.model.SCVI on raw count data of highly variable genes. This was trained on 2 layers with 30 latent variables and gene likelihood modelled as a negative binomial distribution. A total of 50 nearest neighbours was calculated per cell using the function sc.pp.neighbors to generate an adjacency matrix, which was used for calculation of both an initial uniform manifold approximation and projection (UMAP) embedding using function sc.tl.umap with both a minimum distance and spread of 0.5, and Leiden clusters at a resolution of 0.03 using sc.tl.leiden. Initial clusters were then grouped into cell lineages, based on the expression of canonical marker genes (Supplementary Fig. 7), including: B cells, epithelial cells, mast cells, mesenchymal cells, myeloid cells, platelets and T/ILC cells. Markers for these annotations are described in ‘Annotation of cell clusters’.

Within each cell lineage (B cells, epithelial cells, mesenchymal cells, myeloid cells (not including platelets or mast cells) and T/ILC cells), we then further subset cells based on relative thresholds for quality control parameters. Cells with number of unique genes expressed, total counts and percentage of genes originating from the mitochondrial genome that lay outside the median ± 2.5 median absolute deviations were removed. For the epithelial lineage, the median deviation was far greater than that for non-epithelial, and so a window of 2 median absolute deviations was utilized. Additionally, epithelial cells from samples with lower overall depth (indicated by median number of genes expressed per cell <2,000) were removed from the epithelial population only, and any infiltration of blood-derived cells were also removed. Highly variable genes were then re-selected from the full initial set within each lineage independently, and re-integrated using scVI as before, computing 20 nearest neighbours for non-epithelial cells and 50 nearest neighbours for epithelial cells.

Cell types were then identified by clustering within each lineage independently. As described previously⁷⁶, the resolution at which final clusters were derived was determined in a data-driven manner. Cells from each lineage were clustered at resolutions of 0.1 to 1.5 (interval of 0.1) using the sc.tl.leiden function. To determine the reproducibility of the clusters detected at each resolution, two-thirds of the data were then used to train a single layer dense neural network model using keras v2.4.3, and the per-cluster matthews correlation coefficient (MCC) of the fit was assessed in the remaining one-third of data. After clusters that were either comprising >40% cells from a single sample or showed outlying quality control parameters were removed, we selected the resolution at which all remaining clusters had a MCC > 0.75. Additionally, reproducible (MCC > 0.75) clusters demarcated by the expression of genes indicative of poor quality cells, such as MALAT1⁹⁶, and cross-lineage marker genes were also removed. The resolutions utilized for the final cell clusters were: epithelial cells, 0.9; mesenchymal cells, 1.0; myeloid cells, 0.7; B cells, 0.8; and T/ILC cells, 1.1.

Finally, the additionally filtered and clustered per-lineage clusters were then recombined with the mast and platelet cells. This population, known as the ‘atlasing’ dataset, comprised 1,852,681 cells and was re-integrated as described above. With scope to auto-annotate cells with reliable transcriptional profiles that may have been removed by the strict, within-lineage quality control, a CellTypist v1.6.2²⁴ model was trained on the atlasing dataset. This was assessed by quantifying the proportion of cells that were correctly assigned to the same annotations in the initial, manual annotation. Overall, this was achieved at a rate of 97% (median 98%), with only 4 clusters achieving an accuracy <90% and only one <80% (tissue plasma B CD38^hi IgG IgA, 87.8%; tissue plasma B CD38^low, 79%; blood atypical B cells, 88.7%; tissue memory B HIVEP3⁺, 80%). This model was then used to auto-annotate the full set of cells that passed the first round of quality control, and filtered for those with confidence <0.5, resulting in a final dataset that after intersection with genotyped individuals, results in 2,196,874 cells.

Annotation of cell clusters

As described above, cells from all tissues were initially clustered into cell lineages after the first round of quality control and integration. These were grouped and annotated on the basis of canonical marker genes (Supplementary Fig. 8); Epithelial cells were defined by high expression of EPCAM, CDH1 and KRT19; mesenchymal by expression of COL1A1, COL1A2, COL6A2 and VWF; B cells by expression of CD79A, MS4A1 and CD79B; plasma B cells by expression of MZB1 and JCHAIN; T and T/ILC cells by expression of CD3D, CD3E, CD3G, CCR7, IL7R, TRAC and NCAM1; myeloid cells by expression of ITGAM, CD14, CSF1R and TYROBP; mast cells by expression of TPSAB1, TPSB2 and CPA3; and platelets by expression of RUNX1, GATA2, HDC and SLC24A3 (but low expressors of mast cell genes). For within-lineage clustering, B cells and plasma B cells were combined into a single lineage, with mast cells and platelets excluded.

Following clustering within lineage, cell types were first annotated by major population using a series of canonical markers, detailed in subsections below. In the case of T/ILC and mesenchymal populations, this matched their clustering lineage, however; epithelial cells were further subdivided into enterocyte, colonocyte, stem cells and secretory cells; the B cell lineage was re-divided into B cells and plasma B cells; and the myeloid cell lineage were combined with mast and platelet cells. To interpret individual cell clusters, marker genes were calculated within major populations or across targeted subsets of clusters using scanpy’s sc.tl.rank_genes_groups function (method = ‘wilcoxon’). The significant results (FDR < 0.05), along with a discriminative subset that were filtered using scanpy’s sc.tl.filter_rank_genes_groups (min_in_group_fraction=0.5, max_out_group_fraction=0.3, min_fold_change=1.5), were used to annotate and label clusters using expert knowledge. Enrichment of specific populations in a single tissue, or combination of tissues, was determined by the overall proportion of cells from each site (Extended Data Fig. 2).

Epithelial cell lineages

Genes used to demarcate epithelial lineages are presented in Supplementary Fig. 8a. The stem cell major population was defined by expression of LGR5, OLFM4, ASCL2, SMOC2 and RGMB. One of these populations was specific to the rectum, while two were specific to the terminal ileum, and thus are described as ‘rectal’ and ‘intestinal’, respectively. Within the intestinal subsets, one showed greater expression of MKI67, and the other had greater expression of OLFM4 and LGR5, which were thus used to discriminate these populations.

Absorptive epithelial cells were also largely terminal ileum- or rectum-specific (Extended Data Fig. 2). The majority show clear separation in dimensionality reduced space (Fig. 1b), and were grouped into the enterocyte (high expression of RBP2, ANPEP and FABP2) and colonocyte (high expression of CA2, SLC26A2 and FABP1) categories, respectively. The one exception was a population composed of cells from both predominantly diseased terminal ileum and rectal origin, with high expression of OLFM4, thus referred to as ‘epithelial progenitor OLFM4^hi’, grouped into the colonocyte major population due to their embedding in dimensionality reduced space (Fig. 1b). Both enterocyte and colonocyte populations were characterized by the relative expression of marker genes, with discriminative ones being added to their label, including OLFM4 and MKI67 (indicative of proliferation), BEST4, CEACAM7, KRT20, IFI27, ALDOB and RBFOX1.

Secretory cells were broadly more shared across the terminal ileum and rectum, but with some site enriched clusters (Extended Data Fig. 2); goblet cells were identified by expression of MUC2, CLCA1 and FCGBP; tuft cells by expression of POU2F3, SH2D6, LRMP and TRPM5; Paneth cells by expression of DEFA5 and DEFA6; enteroendocrine cells by expression of NTS, PYY, GCG, ISL1; and enterochromaffin cells by expression of TPS1 and CES1.

Mesenchymal cell lineages

Markers used to demarcate mesenchymal cell subsets are displayed in Supplementary Fig. 8b. Fibroblast clusters were identified by their expression of ADAMDEC1, PDGFRA, BMP4, LUM, COL1A2, COL3A1 and COL1A1, subsets of which were then delineated by their expression of PI16, C3, CPM, NPY, GRIN2A and ADAMDEC1 (lamina propria fibroblast), and tissue specificity. Glial cells were identified by expression of NRXN1; smooth muscle cells by expression of ACTA2; pericytes by expression of PDGFRB, RGS5 and NDUFA4L2; endothelial cells by expression of PECAM1, VWF, PLVAP, ACKR1 and CD36. Endothelial cells were further divided into inferred residency by expression of high TLL1, ACKR1 and PECAM1 (venous), modest TLL1 and PECAM1 (lymphatic) and high CD36 greatest PECAM1 and PLVAP (capillary). Myofibroblasts were identified by expression of ACTA2 and MYH11, which were then delineated by relative expression of GREM2 and ADGRL3.

B cell lineages

Markers calculated within plasma B cells were not very discriminatory, however these clusters show a gradient of CD38 expression, which is therefore outlined in the label (Supplementary Fig. 8c). Each of these include a majority of IgA-expressing cells, with one cluster containing a minority of high IgG-expressing cells. For the non-plasma B cells; atypical B cells were identified by their expression of FCRL2, ITGAX and TBX21; naive B cells by expression of IGHD and FCER2; germinal centre/plasmablast B cells by expression of CD19 and TCL1A (for which we also identified a MKI67-expressing, proliferative subset); memory B cells by expression of BANK1 and FCRL2. The several clusters of memory B cells were also distinguished by the relative expression of marker genes; CD27, HIVEP3, PTPRJ, PDE4D and TEX9.

Myeloid cell lineages

Across the myeloid cell lineage, we identified several distinct clusters of monocytes, macrophage and dendritic cells (Supplementary Fig. 8d). Monocytes were defined by high expression of S100A8 and S100A9 and were broadly divisible by anatomical location. Blood monocytes were distinguished by expression of CD14 and FCGR3A (encoding CD16) and tissue monocytes by expression of SOD2, CXCL9 and CXCL10.

Macrophages were predominantly identified in the terminal ileum and rectum and were identified by high expression of ITGA4, CSF1R and MAF. Macrophage subpopulations could thereafter be distinguished by relative expression of C1QA/B/C, CD163, MRC1, SELENOP and HLA-DRA.

Three separate dendritic cell subsets were identified: cDC1s were identifiable from expression of CLEC9A and XCR1, cDC2s were identified by expression of CLEC10A, and pDCs identified by IL3RA.

T cell and ILC lineages

Within the T cell and ILC lineage, T cells were defined by expression of CD3D, CD3E, CD3G, TRAC, TRBC1 and CD96, refined by CD4, CD8A, CD8B, CD3 and GZMK expression (Supplementary Fig. 8e). Naive subsets were identified by high expression of SELL; γδ (GD) by expression of TRGC2 and TRDC (for which we also defined gut- and blood-enriched subsets referred to as ‘GD T KLRG1’ and ‘Tissue resident GD T’; see Extended Data Fig. 2); T_reg cells by expression of FOXP3and CTLA4; Effector memory subsets divided by expression of TIGIT and THEMIS; tissue resident memory with high expression of TRGC2. Two T cell clusters that did not express CD8 or CD4 were also identified, one of which strongly expressed TRDC, and the other that strongly expressed GZMK. NK1 cells were identified by expression of CX3CR1 and FCGR3A, and NK2 cells by XCL1, XCL2 and CD44. ILCs were identified by expression of IL1R1, ALDOC, LST1 and RORC.

Pseudobulking

Pseudobulk expression data were obtained by calculating the mean ln[cp10k + 1] expression across cells of a given annotation, as per previous sc-eQTL mapping benchmarks⁹⁷. This was performed at ‘all cells’, major population and cell-type resolutions, both within and across anatomical sites (Fig. 1a). Expression data were limited to samples with 5 or more cells in the given annotation and only taken forward for further analysis if there were 30 or more individuals contributing to the annotation. In addition, the genes present in each annotation matrix were limited to those with a total count greater than 0 in at least 20% of pseudobulked samples, or 40% of pseudobulked samples for ieQTL mapping. Pseudobulked expression matrices were then inverse normal transformed at the gene level prior to eQTL mapping.

DNA collection and quality control

On the day of the endoscopy, blood (or gastrointestinal tissue, where blood was unavailable) was collected from individuals for DNA array genotyping. DNA was extracted from blood or tissue samples and genotyped with the Applied Biosystems UK Biobank Axiom Array in five batches of samples by YourGene Health. Genotypes were called jointly across all five batches of samples to mitigate potential batch effects. Genotype data were filtered to remove low missingness, duplicate samples, variants not in Hardy–Weinberg equilibrium, individuals with non-European ancestry (as predicted by KING⁹⁸) and first-degree relatives. Variants were lifted over from human genome build 37 to human genome build 38 with LiftOver prior to imputation being performed using the Michigan Imputation Server with the TOPMed reference panel^99,100. Imputed variants with imputation quality lower than 0.3 were removed.

cis-eQTL mapping

eQTL mapping was carried out using TensorQTL¹⁰¹ permutation and nominal modes with default parameters. cis-eQTLs were mapped for variants with minor allele frequency > 0.05 that were within a 2 Mb window (1 Mb in either direction) of the TSS of a given gene. The following model was used for mapping cis-eQTLs:

where Ex is gene expression, G denotes a gene and C is an annotation. GT denotes the genotype alternative allele dosage of the individual (typed or imputed) and PC refers to principal components of genotype and annotation expression. To detect conditionally independent signals, a forwards regression framework was applied, as described previously¹³.

For interaction cis-eQTL mapping, a stricter minor allele frequency filter of 0.1 was applied to each variable level and an expanded model was used:

where I is the interaction variable and GT:I indicates the coefficient for the interaction between genotype and the interaction variable, reported in results. We tested for interactions with sex, age (in bins of 5 years), disease status (CD or healthy), smoking status, and inflammation status (TI-SES-CD score, as defined previously⁷⁶). The reported sex of each sample was verified by Y chromosome carriage and pseudobulk expression of XIST versus all genes on the Y chromosome.

Five genetic principal components were included to control for population structure between individuals. The number of expression principal components to include was determined by repeating the analysis several times, increasing numbers of principal components from 2 to 100 in intervals of 2. The optimal model was selected based on the number of expression principal components that returned the greatest number of eGenes (genome-wide FDR < 0.05), with ties broken by the model with fewest covariates (see ‘Multiple hypothesis testing correction’). For the interaction analysis, any expression principal components which were correlated with the interaction variable were removed from the model (Pearson’s R > 0.25).

Multiple hypothesis testing correction

Correction of P values to account for multiple hypothesis testing in eQTL mapping was done within each cellular annotation independently. This varied slightly for primary, conditional and interaction effects and is detailed below.

Multiple testing correction of primary signals

1. 1)

   Within each gene, multiple variants are being tested for their association with expression. To account for this, expression is permuted 1,000 times, and associations with each variant are recalculated. This procedure allows estimation of a null distribution of associations, and calculation of beta-approximated empirical P values (beta-P). The variant with the smallest beta-P value is then selected for each gene.

2. 2)

   Multiple genes are being tested for eQTL associations within each annotation. A q-value correction¹⁰² is then applied to the variants with the smallest beta-P value per gene, allowing estimation of a genome-wide FDR. We then apply a threshold of FDR < 0.05 to account for both the number of variants being tested per gene and the number of genes, while accommodating for the varied underlying null distributions.

Multiple testing correction of conditionally independent effects

For genes with a significant eQTL association in a cellular annotation (genome-wide FDR < 0.05), subsequent signals are then tested in a forward regression framework by iteratively incorporating variants into the tested model. Firstly, the minimum beta-P value with genome-wide FDR < 0.05 is defined across all genes (min-beta-P). During each round of conditional analysis, associations between variants and gene expression level are re-tested and permuted, and the variant with the smallest beta-P value is selected for each gene. If this variant has a beta-P < min-beta-P from the primary analysis, it is declared a conditionally independent signal and included into the model being tested for that gene in the next round. Because the beta-P must always be less than min-beta-P, conditional signals are iteratively mapped per gene until they would fail to pass the genome-wide FDR < 0.05 from the primary analysis.

Multiple testing correction of interaction effects

Following ieQTL testing, the number of effectively independent variants being tested is estimated for each gene using eigenMT¹⁰³. Nominal P values for ieQTL effects are then adjusted using the Bonferroni method to account for the number of effectively independent tests being performed (emt-p). The variant with the minimum emt-p value is then selected for each gene, and a Benjamini-Hochberg adjustment is then applied to account for the number of genes being tested (adjusted P). Reported ieQTL effects all have adjusted P < 0.05.

Aggregation of QTL effects

Lead variants for the same genetic association may differ across conditions or QTL types (primary, conditional or interaction) due to minor changes in significance. To group lead variants across all tests and define comparable QTL signals, linkage disequilibrium clumping was performed using the ‘–clump’ flag in plink v1.90b6.27¹⁰⁴ and a threshold of r² > 0.5 for variants within 1 Mb with no additional significance filter. This was performed using lead variants from conditionally independent eQTLs (genome-wide FDR < 0.05), ieQTLs (adjusted P < 0.05) and lead variants after intersection with GWAS summary statistics to help map colocalization events to eQTLs.

Inferring capture of our eQTLs in other studies

Harmonized nominal summary statistics for eQTLs from every GTEx v8 bulk RNA sequenced tissue(GTEx Consortium 2020) and sc-eQTLs for peripheral blood mononuclear cells of Yazar et al.¹⁸ were downloaded from eQTL catalogue (https://www.ebi.ac.uk/eqtl/, study IDs QTS000015 and QTS000038, respectively). For eQTLs detected in each of our annotations, we then calculated the enrichment of these associations among low P values, π1 (ref. ¹⁰²), in each external dataset using qvalue v2.34.0¹⁰⁵. This was also performed between eQTLs detected at each resolution in our dataset with eQTLs from GTEx tissues. For comparison with Yazar et al.¹⁸, cell types were grouped into major populations based on nomenclature (B – intermediate, memory, naive; Myeloid – CD14_Mono, CD16_Mono, Platelet, cDC2, pDC; T – CD4_CTL, CD4_Naive, CD4_TCM, CD4_TEM, CD8_Naive, CD8_TCM, CD8_TEM, MAIT, NK, NK_CD56bright, NK_Proliferating, Treg, dnT, gdT; unknown – HSPC, Haematopoietic stem; Plasma – Plasmablast). A detailed summary of these comparisons is shown in the Supplementary Note.

Colocalization

Colocalization was performed in a pairwise fashion between sc-eQTL summary statistics and GWAS summary statistics using coloc v5.2.2 and the coloc.abf function³². GWAS summary statistics from de Lange et al.³ were downloaded from GWAS Catalog (https://www.ebi.ac.uk/gwas/home) (see ‘Data availability’ and ‘Code availability’). Analyses were conducted in 2 Mb cis windows centred around the lead single-nucleotide polymorphism (SNP) of each eGene with a permuted FDR-corrected P value < 0.05 and using GWAS loci which were significant at the genome-wide suggestive threshold (P value < 5 × 10⁻⁵). Analysis was conducted using default prior probabilities of p1 = p2 = 10⁻⁴ and p12 = 10⁻⁵ where p1 and p2 are the prior probabilities of true genetic associations in each input dataset and p12 is the probability of true associations in both datasets. For colocalization between GWAS and ieQTLs, approximate Bayes factors were computed from the effect sizes and standard errors of the GT:I term, and the prior standard deviation on the true effect size was estimated as 0.3 based on the distribution of GT:I effect sizes across all nominal results on chromosome 21. A posterior probability of colocalization (PP_h4) greater than 0.75 was considered to be sufficient evidence for colocalization of genetic signals. Colocalizations which fell within the HLA on chromosome 6 (build 38, positions 28510120–33480577) were discarded as likely false positives due to the complicated patterns of linkage disequilibrium in this region. Colocalization events were assigned to the closest existing IBD locus, if the lead colocalized GWAS variant was within 500 kb of a known IBD locus as defined by the interval comprising ±500 kb containing: (1) the most likely causal variant from published fine-mapping analysis^3,106; or alternatively, if the locus has now been previously fine mapped, (2) from the GWAS lead variant as defined in refs. ^3,34.

We used the R package otargen (v1.1.5) to download all colocalizations between genes with molecular QTLs and IBD, CD and UC GWAS loci from the Open Targets genetics portal (accessed 22 January 2025). We retained any colocalizations with PP_h4 > 0.75 and assigned to IBD loci as previously described for our colocalizations. For comparison of disease effector genes nominated in this work with that of Cuomo et al.³⁵, as their colocalization is limited to IBD and with a threshold of PP_h4 > 0.8, we applied the same filtration to our own associations.

To check if disease effector genes identified at high resolutions were more specifically expressed, genes were grouped into the lowest resolution at which a colocalization was detected on their expression. The maximum specificity of gene expression at the cell-type level was then calculated using CELLEX¹⁰⁷, and compared across resolutions using a Wilcoxon test.

Detection of colocalization events between two phenotypes is highly dependent on the significance of those effects. This is therefore underpinned by the statistical power in those tests, driven by the total number of observations. As this varies substantially for sc-eQTL tests (Fig. 1b and Supplementary Fig. 3), inherent power to determine eQTL and colocalization events also varies. To circumvent the influence of varied statistical power in the nomination of cell types within which colocalization events are active, for each colocalization, we quantified the proportion of variance in disease effector gene expression (r²) that was explained by the lead eQTL variant after intersection with the GWAS variants across each condition. This was calculated by linear regression of disease effector gene expression in the inverse normal transformed pseudobulk data on genotype.

To assess the evidence for enrichment of Notch and Wnt pathway regulators in our disease effector genes, we performed a gene ontology over-representation test of ‘Notch signalling pathway’ (GO:0007219) and ‘Wnt signalling pathway’ (GO:0016055) in our disease effector genes against a background of all eGenes using clusterProfiler v4.10.1¹⁰⁸. This was performed using default parameters except minGSSize=3, pvalueCutoff = 1, qvalueCutoff = 1. These gene sets were also used to quantify the number of disease effector genes previously and presently nominated in these pathways.

Enrichment of eQTLs in genomic annotations

To calculate the enrichment of lead variants from eQTLs in ENCODE and FANTOM genomic annotations, we followed previous approaches¹², available on GitHub https://github.com/hakha-most/gwas_eqtl. In brief, eQTLs were first grouped by the minimum resolution at which they were detected, and their presence in functional elements was calculated. To calculate bootstrapped confidence intervals, we then sampled with replacement from independent linkage disequilibrium blocks (calculated previously¹⁰⁹) that contained eQTLs detected at each resolution one thousand times. Enrichment values in different annotations were then first normalized to that of random variants, before values in enhancers being normalized to those in promoters, in order to identify relative enrichment in these annotations specifically. Distal and proximal enhancer-like annotations from ENCODE were grouped together. P values were calculated by two-tailed t-tests comparing the bootstrapped distributions, with the 25th, 50th and 75th percentiles displayed in the enrichment plot.

Therapeutic evidence analysis

We downloaded data on existing therapeutics and their clinical approval stage from ChEMBL and Open Targets^68,110. Therapeutics were grouped into broad disease terms using the Ontology Lookup Service¹¹¹ and filtered for those which directly targeted genes which we had colocalized with UC, CD or IBD (PP_h4 > 0.75, within 0.5 Mb of existing IBD locus) but had not yet been trialed as IBD treatments according to ChEMBL.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.