Systematic identification of disease-causing promoter and untranslated region variants in 8,040 undiagnosed individuals with rare disease

Identifying known disease genes using PanelApp

PanelApp gene panels were obtained from panelapp.genomicsengland.co.uk using lynx v2.8.9(16), on 12/09/2022. These were filtered to include only the 6,504 genes where the strength of association for one or more gene panel was “green”, corresponding to those with a confident link to the phenotype.

We further filtered to only include genes known to cause a disorder with a dominant mode of inheritance (MOI), inclusive of any genes associated with both dominant and recessive phenotypes. Finally, we selected only genes with transcripts in the MANE v1.0 dataset (17). In total we included 1,536 genes/1,567 transcripts.

Annotating non-coding regions of interest

Transcripts were defined using MANE v1.0, inclusive of 19,062 MANE ‘select’ and 58 ‘plus clinical’ transcripts(17). UTR exon and intron coordinates were taken directly from the MANE .gff file.

Proximal promoter regions were defined using candidate cis regulatory elements (cCREs) obtained from ENCODE(18). Accurate promoter definition is hampered by their tissue specificity. In tissues where a promoter is inactive, it is often marked by a minimal nucleosome free region, but this region may be expanded when the promoter is active. To account for this, as well as promoters that are not annotated at all in the ENCODE dataset, we calculated the average size of all ‘promoter-like’ cCREs that overlap with TSS of MANE transcripts. We calculated the 25th and 75th percentiles of the distribution of distances these cCREs extend up-(25%=181bp; 75%=266bp) and down-stream (25%=67bp; 75%=139bp) of the TSS (Supplementary Figure 1). The 25th percentiles (-181bp to +67bp from TSS) were used to define a ‘minimal’ promoter region.

For transcripts with a cCRE that overlaps the TSS:

• ● If the cCRE extends ≥ 181bp upstream and ≥ 67bp downstream of the TSS (i.e. at least the minimal 25th percentile definition) the exact region defined by the cCRE is used (Supplementary Figure 1d; n=7,368).

• ● If the cCRE falls short in either (or both) direction(s), it is extended to reach the 25th percentile distance in that/those direction(s) (Supplementary Figure 1e; extended upstream n=2,953; extended downstream n=2,918, extended in both directions n=464).

For transcripts with no TSS overlapping cCRE (n=5,417), the 75th percentiles are used to define a promoter region that stretches 266bp up- and 139bp down-stream of the TSS.

To ensure identified variants don’t have a protein-coding impact(19) we used bedtools(19) to exclude any positions that intersect with a CDS position in any MANE transcript. In total, we defined 20,417,669 near-coding bases across the 1,567 green PanelApp genes, for an average of 13,030 bases per transcript (min=264, max=791,387), and between 17 and 18,786 per region (Supplementary Table 1). The final set of near-coding regions defined across all green PanelApp genes is in Supplementary Table 2.

Identifying and filtering de novo variants

We used a dataset of previously identified and filtered de novo variants (DNVs) within GEL(20), accessed using the RLabKey API(21). We filtered individuals to remove any with subsequently withdrawn consent, and to only include those with a ‘participant type’ of ‘proband’, where neither parent was classified as ‘affected’ or had any associated Human Phenotype Ontology (HPO) terms, and for whom variant calls were on the GRCh38 reference genome. This resulted in an initial set of 9,665 trio probands.

Variants were filtered to only those that passed the most stringent set of GEL filters(22). We removed variants with allele frequency (AF) ≥ 0.00005 and allele count (AC) >=5 in the GEL defined set of 55,603 unrelated individuals or with AF ≥ 0.0005 in any of the major population groups in gnomAD v3.1.1.(23). We restricted our analyses to DNVs within our defined near-coding regions of PanelApp genes with high confidence phenotypic associations (flagged as ‘green’ genes) for the individual’s phenotype. Finally, we excluded participants with an identified coding diagnosis (see below). This resulted in a set of 1,278 DNVs, in 1,094 probands.

Identifying individuals with existing diagnoses

A list of all participants for whom a confirmed diagnosis was recorded was obtained from the Genomics England ‘exit questionnaire’ table, identified with those for whom the family case was flagged as “solved”. The associated variants were cross referenced with MANE v1.0 coding regions and our promoter, UTR, and UTR intronic regions, with variants mapped onto GRCh38 by the ‘LiftOver’(24) tool, where required.

Region level variant annotations

Variants for both the de novo and burden testing arms of our analysis underwent initial annotation using Ensembl’s variant effect predictor (VEP) v99.1(25) with UTRannotator, SpliceAI v1.3, and CADD v.1.6 plugins(26), as well as custom annotations for PhyloP 100 way vertebrate conservation scores(27), and ClinVar (28) clinical significance annotations (accessed 2022/08/12).

5’UTR variants annotated by UTRannotator(29) were filtered to identify those with the highest likelihood of disrupting translation. To this end, we extracted all variants annotated as:

• ● uAUG gain resulting in creation of an overlapping open reading frame (oORF) with a strong or moderate Kozak consensus sequence;

• ● uSTOP loss, with no alternate stop prior to the CDS start (i.e. also resulting in an oORF) and a strong or moderate Kozak consensus sequence;

• ● uAUG loss, with strong Kozak consensus sequence;

• ● uFrameshift resulting in an oORF with a strong or moderate Kozak consensus sequence.

For SpliceAI scores, we took the highest delta value across all predictions, with a cutoff of 0.2 for the de novo analysis and 0.5 for the burden testing.

Variants across all regions with a ClinVar annotation indicating a benign or protective effect (benign, likely benign, benign/likely benign, protective) were excluded.

Cut-offs for CADD PHRED(26) (25.3) and PhyloP(27) (7.367) scores were taken as the supporting evidence thresholds from Pejaver et al(30). For variants with multiple recorded scores, the maximum was taken. We note that these scores were calibrated for missense variants, but no alternative exists for non-coding region variants due to the paucity of variants available to benchmark against. Due to this, CADD PHRED and PhyloP scores were not used to annotate variants in deep intronic regions (>20bp from the end of the exon), to reduce noise.

Internal ribosome entry site (IRES) data were obtained from IRESbase(31) on 23/08/2022, microRNA (miRNA) binding sites were obtained from the literature(32–35), and downstream open reading frame (dORF) coordinates were obtained from Chothani et al(36). For each, locations were cross referenced with our variant positions, and any intersecting variant flagged.

Given the large proportion of variants that fall into miRNA and IRES sites, we excluded any variants that also had a CADD ≤ 22.7 or PhyloP ≤ 1.879. These scores were suggested by Pejaver et al as supporting evidence for a benign classification(30).

Kozak consensus sequence variants in the -3 position were identified with reference to MANE v1.0 CDS start positions (i.e. the R of the gccRccAUGG motif). Any variant that changed a reference A or G to a C or T was annotated as potentially Kozak disrupting(37).

RNA binding protein predictions were generated using the methods detailed in Findlay et al(38) for all possible variants within motifs that are proximal to ENCODE eCLIP sites, that are also high affinity sites as predicted by RBPamp(39). These were intersected with our variants and filtered to retain only those with a reference affinity of ≥ 0.1 and with an impact of ‘loss of binding’ predicted by the RBP binding affinity model (defined as alternative allele affinity / reference allele affinity < ⅓).

Using MANE v1.0 mRNA sequences, we identified the locations of all 3’UTR AAUAAA and AUUAAA polyA signal motifs. We filtered intersecting variants to those that did not result in the creation of an alternative known motif (AAUAAA, AUUAAA, AGUAAA, UAUAAA, CAUAAA, GAUAAA, AAUAUA, AAUACA, AAUAGA, AAAAAG, ACUAAA, or AAAAAA)(40).

Transcription factor binding site (TFBS) locations were obtained from ENCODE(41) and converted using bigBedToBed(24) on the command line, resulting in 4,465,728 TFBS footprints. Any variant not within a footprint identified by ENCODE as falling within the ‘core’ region of a DNase I hypersensitive (DHS) peak was excluded. Remaining variants were annotated using FABIAN(42), limiting to only transcription factor flexible models as these have been shown to outperform positional weight matrices(43). The resultant data was transformed to produce one score per transcription factor, per variant: Where A is each model’s predicted change in binding affinity and ‘N’ is the total number of these predictions provided for that transcription factor. Scores ≥ 0.04 were recorded as predicted gain and those ≤ -0.04 as predicted loss. For each variant we then calculated the mean gain/loss/total scores and retained any variant with a loss score ≤ -0.4.

Clinical Review of candidate variants

For each participant carrying a candidate diagnostic de novo variant, we compared the similarity between the HPO terms assigned at recruitment with the phenotype expected for a heterozygous loss of function variant in the gene. Variants were interpreted under the assumption that they caused loss-of-function (LoF) and were of high penetrance. Expected phenotypes for each gene were sought from OMIM and the published literature. Where we identified a plausible phenotype match, we raised a clinical collaboration request with Genomics England to confirm or refute our findings by collaboration with the recruiting clinician.

Defining case and control sets for burden testing

From GEL version 15, we selected participants meeting all of the following criteria:

1. Variants called on genome build ‘GRCh38 and with delivery version ‘V4’

2. Consent not subsequently withdrawn

3. Karyotype one of ‘XX’, ‘XY’, ‘NA’, ‘Other’ and karyotypic and phenotypic sex not in conflict

Cases were defined as:

1. Individuals with a participant type of ‘proband’

2. With at least one ‘green’ PanelApp gene in a virtual gene panel assigned to them

3. Without an existing coding diagnosis (see above).

Controls were taken as the unaffected parents of participants with rare disease. Defined as:

1. Participant type is ‘Mother’ or ‘Father’

2. Affected status is ‘Unaffected’

3. No recorded HPO terms

The genetically-inferred ancestry of each participant, as calculated by GEL, was obtained from LabKey. Participants with a single origin ancestry match of 99% or greater were retained and defined as that ancestry(44). Through this approach, we defined a total of 19,220 cases and 20,683 controls.

Filtering aggregated variants

Variants within MANEv1.0 transcripts, for all potential case and control participants that passed all internal QC filters were extracted from the aggregated variant VCF files in GEL(45).

In line with recommendations from Pedersen et al(46) we filtered variants to those with genotype quality (GQ) ≥ 20, read depth (DP) ≥ 10, missingness ≤ 5% heterozygous allele balance (AB) 0.2 ≤ AB ≤ 0.8, and homozygous AB ≤ 0.02. If a variant call failed one or more of these filters in 25% or more cases that call was excluded. We further filtered to only those with GEL internal and gnomAD (v3.1.1; in any population) AF ≤ 0.0001. We retained a total of 18,498,584 variants, a mean variant count per individual of 463.59 (461.74 in cases and 465.74 in controls).

Participant Matching

To exclude any individuals with very high numbers of called variants (suggestive of systematic error), we calculated a population-specific Z-score per participant as follows: Where ‘x’ is the variant count in that participant, across all MANE transcripts, from the start of the promoter to the transcript end, ‘μ ‘is the population mean, and ‘σ’ is the population standard deviation, where the population is all individuals defined of the same genetic ancestry (see above). Participants with a Z-score of ± 2 were dropped (N=1,560, 728 probands, 832 controls) resulting in a set of 18,492 probands and 19,851 control participants.

Within each cohort, we removed individuals with KING score(47) ≥ 0.0442 within the Genomics England relatedness data(44), indicative of being a 3rd degree relative, by randomly selecting one participant for removal in an iterative process until no further relatedness in individuals was detected.

We then matched each proband 1:1 with a single control participant by sex and genetically-inferred ancestry, ensuring that the matched proband and control did not share a family ID. The resultant matched cohort consisted of 18,304 probands, paired with 11,641 unique controls. To avoid potential biases when matching participants caused by low population numbers, we limited to genetically-inferred ancestries where the number of both case and control participants was > 200. This resulted in a cohort of 17,641 case probands, and 11,227 control participants with either European or South Asian genetically-inferred ancestry (Supplementary Table 3). To limit bias due to the presence of large gene panels, this cohort was then limited to probands who had 100 or fewer green dominant genes assigned (Supplementary Figure 2), resulting in a final cohort of 7,862 probands and 6,371 matched controls.

Burden testing

Aggregated variants filtered as above were further restricted to those with AF ≤ 0.00005 for both internal and gnomAD major population frequencies and to exclude any with an allele count (AC) across the entirety of AggV2 of >= 5. These 1,079,616 variants were annotated and filtered with reference to the annotations described above, with the addition of a more stringent SpliceAI threshold of 0.5 (Supplementary Figure 3).

A simple burden test was performed across all defined near-coding region and variant annotations comparing individuals that had one or more annotated variants meeting our criteria in any near coding region to those that did not, using a Fisher’s exact test performed in R (Supplementary Table 3). The test was repeated for each region annotation separately, using Bonferroni correction for multiple testing.

To estimate the number of participants required to see a significant enrichment across all region and variant annotations, we iteratively increased the number of case and control participants by 1, while maintaining the proportion of observed cases and controls with candidate variants. Fisher’s tests were performed for each iteration, until the resulting P-value was ≤0.0031, a Bonferroni adjusted threshold accounting for 16 tests.

RNA sequencing

Blood was collected from a subset of 100,000 Genomes Project probands in PaxGene tubes to preserve RNA at the time of recruitment. RNA was extracted, depleted of globin and ribosomal RNAs, and subjected to sequencing by Illumina using 100bp paired-end reads, with a mean of 102M mapped reads per individual. Alignment was performed using Illumina’s DRAGEN pipeline. IGV(48) was used to inspect sequencing reads and generate Sashimi plots to show splicing junctions supported by 5 or more reads in areas of interest. FRASER2(49) and OUTRIDER(50) were used to detect abnormal splicing events and expression differences with 499 samples used as controls.

DNA methylation analyses

DNA methylation array testing was performed on a diagnostic basis as described previously(51,52).

Availability of data and materials

The datasets supporting the conclusions of this article are available in the near-coding annotation github repository (https://github.com/Computational-Rare-Disease-Genomics-WHG/Near_coding_annotation).

Unless otherwise stated, all analyses were performed using R Statistical software version 4.0.2(53), with the packages ‘dplyr’(54), ‘tidyr’(55), ‘stringr’(56), ‘Rlabkey’(21), ’UpSetR’(57), and ‘ggplot2’(58).

Strict region-specific filtering prioritises likely deleterious de novo promoter and UTR variants

We identified 767,063 rare (AF ≤ 0.005%) high-confidence DNVs in 10,665 trio probands in GEL (71.9 per proband), 685,438 of these variants were in participants whose parents were classed as unaffected and had no HPO terms (9,665 probands). Of the remaining probands, 8,040 did not have an existing confirmed diagnosis attributed to a coding variant (576,030 variants, 71.6 per proband). We filtered to include only DNVs in UTR exons and introns (both defined using MANEv1.0 transcripts), and promoter regions (defined by ENCODE candidate cis regulatory elements; see Methods). Henceforth, we collectively refer to these regions as ‘near-coding’. We limited our analysis to variants which fell in or near known monogenic disorder genes (3,316 variants) and filtered these to genes which could be plausibly associated with the participant’s phenotype. Accordingly, we filtered for DNVs of genes flagged as ‘green’ in one or more PanelApp(59) gene panel(s) assigned to the individual, and which were associated with disorders with a dominant mode of inheritance. Of note, 309 probands did not have any assigned green dominant genes. In total, we proceeded with 1,311 candidate DNVs in 1,118 probands.

To identify likely disease-causing DNVs we used a region-specific annotation and filtering approach. We prioritised 5’UTR variants that create uAUGs or disrupt uORFs using UTRannotator(29), that overlap IRES defined by IRESbase(31), or that lead to disruption of the Kozak consensus sequence(37). 3’UTR variants were prioritised if they disrupt a polyadenylation site or signal sequence, disrupt a miRNA binding site, disrupt an RBP motif, or if they disrupt the start/stop of a dORF with evidence of translation from ribosome profiling (from Chothani et al(36)). Given the large numbers of variants annotated as within IRES or miRNA binding sites, these variants were further filtered to remove any with CADD (<22.7) or PhyloP (<1.879) scores in support of being benign(30). Across all UTR exons, and in 5’ and 3’ UTR introns, variants with SpliceAI masked delta scores ≥ 0.2 were prioritised. Promoter variants were prioritised if they are predicted to disrupt a transcription factor binding site using FABIAN(42). Finally, across all region annotations, variants with a CADD score ≥ 25.3 and/or a PhyloP score ≥ 7.367 were flagged (thresholds taken from Pejaver et al(30)). After filtering to only include variants with one or more of these variant annotations, we retained eleven candidate DNVs (0.8% of the initial 1,311 DNVs) each found in a different individual (Figure 1; Table 1).

• Download figure
• Open in new tab

Figure 1: Prioritised de novo variants split by region and variant annotations.

DNVs were identified from the Genomics England de novo dataset in the following regions: Promoter (mustard), UTR exons (raspberry), UTR/Promoter overlapping region (mustard and raspberry stripes), and UTR introns (teal). The gene names corresponding to identified DNVs are written above the corresponding bar. Those in black represent likely diagnoses (nine probands), with those in grey not being a good phenotypic match (two probands). Novel potential diagnoses are marked by an asterisk. Vertical bars in the top panel denote the number of variants identified with specific region and variant annotations that are represented by the bar colour (region annotations), and in the upset plot below (variant annotations). The total number of DNVs with each variant annotation is shown by the horizontal bars to the left of the upset.

Promoter and UTR DNVs provide a diagnosis for undiagnosed individuals with rare disease

Of the eleven remaining candidate variants, nine (82%) were assessed to be a good match for the individual’s phenotype after detailed clinical review (see methods). Three of these had been flagged as diagnostic variants in GEL (in the ‘exit questionnaire’ table) prior to starting this work: two 5’UTR splicing variants in PAX6 in two individuals with aniridia (OMIM:617141) and one 5’UTR splicing variant in RPL26 in an individual with a previously undiagnosed monogenic disorder. A further variant, a 5’UTR variant that creates an upstream start codon in MEF2C, we previously identified as occurring de novo in three unrelated individuals with severe developmental disorders(2). Our approach successfully prioritised all rare DNVs within our candidate regions that had previously been identified as likely diagnostic in GEL. Together, these data demonstrate that our pipeline effectively identifies known diagnostic variants.

Four of the remaining five variants represent likely new diagnoses: (1) a 5’UTR uAUG-creating variant in SLC2A1 in a patient with GLUT1 deficiency syndrome (OMIM:606777) that was not flagged by GEL as diagnostic, but that has been published previously(3) (Figure 2A). This uAUG is created into a strong start codon context and functional studies support its translation(3). Translation from this uAUG will prevent translation of the downstream CDS, leading to loss-of-function (Figure 2A). After returning this diagnosis to the recruiting clinical team it was classified as Likely Pathogenic and the individual is now on treatment; (2) A NIPBL splice disrupting (SpliceAI=0.24) variant 17 bp upstream of the final 5’UTR acceptor site in a participant with a phenotype closely related to Cornelia de Lange syndrome (OMIM:122470, Figure 2B). This variant introduces an AG dinucleotide which is predicted to result in a premature acceptor, however, the positioning of this within the ‘AG exclusion zone’ may also cause skipping of the exon containing the CDS start codon or other splice defects(60) (Figure 2B). The exact impact of this variant will need to be confirmed through RNA studies, but RNA was not available for the patient; (3) A promoter variant that is located in a highly evolutionarily conserved position (PhyloP=7.426) 13 bp upstream of the TSS of ZBTB18 in a participant with Intellectual disability; (4) A 5’UTR splice-site variant in SETD5 in an individual originally suspected to have Silver Russell Syndrome (OMIM:180860). This variant is predicted to result in loss of the splice donor (SpliceAI=0.97) of the first 5’UTR exon at the canonical +1 position. DNA methylation signature analysis in this patient revealed an episignature consistent with SETD5-related neurodevelopmental disorder (Figure 2C) and no other candidate variants were identified after screening the protein-coding regions of SETD5.

• Download figure
• Open in new tab

Figure 2: Candidate diagnostic de novo variants.

A. Gene diagram showing the creation of an out of frame overlapping ORF (oORF; in red) in the SLC2A1 gene in the proband. B. Illustration of the AG exclusion zone in the NIPBL gene. The T>A variant at the -17 position is marked in red, the most strongly predicted branch point (Branchpointer(62) 0.48), directly upstream of the AG exclusion zone is shown in blue. C. Multidimensional scaling plot showing differential methylation in SETD5. The position of both variants found in this gene are shown as red dotted lines. D. Sashimi plot showing aberrant splicing in the MANE Plus clinical transcript ENST00000371085. The proband shows some retention of the intron containing the variant (which is marked by a red dotted line) and increased skipping of the following exon compared to the controls (6.06X% vs 0.65X% and 1X%).

We also prioritised a cryptic splice variant in GNAS (SpliceAI=0.67) in a participant with hypothyroidism. Whilst we originally annotated this variant as within a 3’UTR intron for the MANE Plus Clinical transcript, the intron is between two CDS exons of the MANE Select transcript. Blood RNA-sequencing from the patient showed evidence of abnormal splicing of the MANE Select transcript, including intron retention (FRASER2 adjusted P=1.67×10^-23), and significantly reduced expression (OUTRIDER adjusted P=0.0019; fold-change=0.66; Figure 2D), however the exact mechanism through which this variant could lead to disease is unclear.

For all candidate variants, we checked whether they were found in any other individuals across the full GEL cohort (i.e. not limiting to full trios or DNVs). Whilst we did not observe recurrence of any of the exact variants identified, we did identify a second participant with a different SETD5 variant at the same genomic position (chr3:9397974 CAAGGT>C, hg38). On closer investigation, this variant is consistent with a germline de novo, but it fell just below the required coverage in one parent so it was excluded from the conservative high confidence de novo callset. DNA methylation signature analysis also confirmed SETD5 as the diagnosis in this individual (Figure 2C). In total, we identified a likely disease-causing ‘near-coding’ DNV in ten of 8,040 individuals (0.0012%; nine initially prioritised variants and one additional SETD5 variant) who did not previously have a coding diagnosis. We classified all six newly identified variants as Likely Pathogenic following the ACMG/AMP guidelines (Supplementary Table 4)(8,61).

Burden testing does not detect a significant enrichment of variants with any collective region or variant annotation

Given we were able to identify disease-causing near-coding variants using our region-based filtering pipeline, we sought to further use this approach to quantify the enrichment of potentially damaging promoter and UTR variants. However, given the small number of trios within GEL with an unaffected child, and the fact that mutational models to directly assess enrichment of de novo variants (by comparing observed to expected numbers) have not been well calibrated for non-coding regions, specifically struggling with the 5’ end of genes(63), we instead used the full aggregated set of inherited and de novo variants for our analysis. We matched each of 19,220 rare disease probands without a recorded protein-coding diagnosis, with replacement, to an unrelated unaffected individual from within the rare disease arm (unaffected parents) of GEL as a control on sex and genetically-inferred ancestry (see methods). The control individual was assigned the same dominant, green panelApp genes as had been assigned to the rare disease proband by GEL, allowing us to control for gene and region level differences in mutability (Figure 3). Given the disparity in panel size, with many probands having over 500 assigned green dominant genes (Supplementary Figure 2), to reduce noise we filtered probands to include only those to whom 100 or fewer green dominant PanelApp genes had been assigned (28% of all probands).

• Download figure
• Open in new tab

Figure 3: Burden testing results.

Counts of variants and odd ratios (log10) testing for an enrichment of variants in cases compared to matched control participants (Fisher’s test), collectively by A. region annotation, and B. variant annotation. Annotation groups with fewer than 10 participants are omitted. Error bars represent 95% confidence intervals. Variants in 5’UTRs (P=0.032) and variants with SpliceAI ≥0.5 (P=0.008) are enriched in cases over matched controls, but neither is significant after correcting for multiple testing (Bonferroni threshold adjusting for 16 tests =0.0031). Full results are in Supplementary Table 5.

After participant matching, we analysed a final set of 7,862 probands and 6,371 matched controls (1,295 matched controls were partnered with more than one proband). For all individuals, we extracted variants from GEL’s aggregated variant dataset (AggV2) and filtered these using the same region-specific criteria applied to the de novo variants. Given that we used a high sensitivity SpliceAI threshold to prioritise DNVs with a high prior of pathogenicity, this was raised to a stricter cutoff of 0.5 for this analysis (Supplementary Figure 3). As we are not powered to analyse individual genes or gene-regions, we performed burden testing collectively across all prioritised variants with the same regional (e.g., 5’UTR) or variant-level (e.g. SpliceAI) annotations, across all participants and their assigned green genes. Whilst we observed a greater number of probands with prioritised variants compared with matched controls for the majority of regional and variant-level annotations we identified a greater number of probands with prioritised variants compared with matched controls, no specific annotation was significantly enriched for variants in cases after correcting for multiple testing (Figure 3; Supplementary Table 5). We also did not observe a significant enrichment when combining across all regions and variant annotations (Fisher’s P=0.109, OR=1.09, 95% CI=[0.981,1.210]).

Assuming a constant ratio between case and control variants, and hence ORs, we would need an estimated 11,610 cases and controls for significance at P<0.05 in the combined test across all region and variant annotations, and 26,066 for significance at a study-wide P-value threshold of <0.0031, correcting for 16 tests (Supplementary Figure 4).