Harnessing the 100,000 Genomes Project whole genome sequencing data - an unbiased systematic tool to filter by biologically validated regions of functionality

Defining test and validation datasets through the 100,000 Genomes Project

The 100,000 Genomes Project was set up by the UK Department of Health and Social Security in 2013, to sequence whole genomes from National Health Service (NHS) patients recruited through one of 13 Genomic Medicine Centres (GMCs). For successfully nominated Rare Diseases and Cancers, patient/family DNA sampling was accompanied by phenotypic data entries and integration with NHS health records. Anonymised raw sequencing data were made available in the Genomics England Research Environment.[12] Separately, Genomics England performed data alignments, and variant classifications for the clinical diagnostic pipeline: For known disease-causing genes, loss of function variants were defined as Tier 1, other variants as Tier 2 variants. Variants in other protein coding genes not in the virtual gene panel for the disease or any listed patient phenotype were classified as Tier 3 variants.

In order to provide test and validation sets, we focused on WGS data from patients with one particular rare disease, hereditary hemorrhagic telangiectasia (HHT).[13] Due to specific circumstances in NHS commissioning, a group of HHT patients who had not previously undergone genetic testing had been recruited, and a subset of these provided a validation dataset, as the Tiering process identified variants in known disease-causing genes.

Generating GROFFFY

Genomic coordinates of regions included in the GROFFFY were generated from publicly available databases using the Imperial College High Performing Computing service. As our goal was to miss no functional information, wherever possible, experimentally-derived biological data were used in preference to computational predicted files with potential for false negatives. Genomic coordinates were extracted from data aligned to GRCh38 genome build[14], and excluded data originating in cancer cells. Full details are provide in the online Data Supplement (Supplementary Methods and Tables 1,2). Briefly: Genomic coordinates for transcribed loci including protein-coding genes, pseudogenes and long non-coding RNAs (lncRNA) were extracted from GENCODE Version 31[10]. LncRNA coordinates were also downloaded from lncipedia version 5.2,[15] and micro RNAs (miRNA) coordinates downloaded from miRbase release 22.1.[16]Genomic coordinates for candidate regulatory element (cRE) regions were generated by merging 18,828 bed files identifying regions where DNA had interaction with different elements, such as histone modifications, protein binding, and nuclease accessibility. Thus, DNA binding regions from 2,823 experiments performed in 9 different labs, 4,409 bed files identifying representative DNase hypersensitivity sites from 631 experiments, contained within the ENCODE Encyclopedia registry [7] were downloaded, together with coordinates for CpG islands in the UCSC database.[9]Genomic coordinates of regions excluded were based on the ENCODE blacklist (accession ENCSR636HFF) [5] downloaded from ENCODE,[6,7] and RepeatMasker downloaded from the UCSC database [8,9].

Analysis of Whole Genome Sequencing Data

The final genomics coordinates for GROFFFY were transferred into the Genomics England Research Environment. The full pipeline employed for the current manuscript is illustrated as Supplementary Figure 1. At all stages, the number of target set DNA variants after exclusion was recorded and outputted to log files. Stepwise, filter by value functions were written to apply allele frequency filters (excluding variants where general allele frequency exceeded 0.0002[3] in either the 1000 Genome Project[17] or gnomAD database[18]); to exclude synonymous variants (note that synonymous variants within splice regions were retained as they were annotated as “splice region” variants); and to remove additional variants in annotated WGS files for the 34 “Validation Set” HHT DNAs where a causative variant had already been identified in known HHT genes through Tiering. Combined Annotation-Dependent depletion score annotation, using version 1.5 of the CADD database,[19] was quite slow, and the process was put towards the end of the analysis pipeline, so that there were fewer variants that needed to be annotated.

• Download figure
• Open in new tab

Figure 1: GROFFFY categorization of the human genome:

A) Heatmap displaying 61,799 data points at 500kb resolution mapped against GROFFFY. Colors represent the percentage of each 500kb region included in GROFFFY with red highest (> 90%), orange ∼70%, yellow ∼40%, and blue <30%. All chromosomes are left-aligned, and as all are shorter than chromosome 1, the right-aligned blue regions indicate the absence of chromosomal material for the indicated chromosome. Other broad blue regions highlight regions particularly devoid of genes and regulatory elements, for example the pericentromeric long arms of chromosomes 1 and 9, and the short arms of chromosomes 13, 14, 15. B) The human karyotype for comparison, plotted in R using sourced ideograms with centromeres highlighted in red. (Note consistency with sizes in A). C) Randomly chosen region of chromosome 3 (chr3:25,597,986-2-25,783,443). The top 4 tracks illustrate GENCODE gene annotations, CpG islands, long non-coding RNAs, and miRNAs. The final filter is shown as the lowest track. Note this contains both intra and intergenic regions for the region, and that the raw data were not subjected to any processed annotation tracks.

GROFFFY defines biologically validated regions of functionality in the human genome

GROFFFY essentially excludes biologically less important regions of the genome by only including regions where biological experiments have generated evidence in favour of functional roles. Due to many overlapping regions, merged bed files were approximately 10-fold the size of the human genome (e.g. for DNA-binding regions, 49,941,695 Kb) before sorting. Following sorting, the GROFFFY filter region based on positive selection of transcribed loci and candidate regulatory element (cREs), and masking of repetitive regions, included 44.4% of the human genome. Genome-wide illustration of GROFFFY by heat mapping datapoints at 500kb resolution highlighted genomic regions that are enriched or devoid of genes and regulatory elements (Figure 1A, 1B). To illustrate GROFFFY at higher resolution, we randomly selected a 185kb region containing two protein-coding genes. Figure 1C illustrates the raw experimental data of the GROFFFY region as the lowest track.

GROFFFY prioritisation of variants within WGS data

GROFFFY was applied to WGS-derived DNA variants in individuals with HHT [13] recruited to the 100,000 Genomes Project. The experimental target dataset comprised 98 whole genomes. The validation dataset included 34 whole genomes where a disease-causing variant had already been identified through clinical diagnostic pipelines.

We focused first on the experimental target dataset of 98 whole genomes. The scale of the bioinformatics challenge was emphasized by a total of 26,271,838 unique changes compared to the human genome reference sequence (GRCh38 genome build [14]). Pre filtration, the number of DNA variants per individual target dataset DNA ranged from 4,786,039 to 5,070,340 (mean 4,867,167, SD 48,868), and approximately one third (32.28%), lay within intergenic areas.

Applying GROFFFY as a first filter reduced the mean number of variants by 2,812,015 per DNA, while retaining variants within intergenic areas (Supplementary Table 3A, Supplementary Figure 2A). Post filtration, intergenic variants represented 33.02% of the total variant number.

• Download figure
• Open in new tab

Figure 2: GROFFFY validation by retention of pathogenic variants:

Comparison of variants in the 34 validation dataset whole genomes by Combined Annotation-Dependent depletion (CADD) scores indicative of likely pathogenicity pre, and post filtration. A): Total number of variants with indicated CADD scores (note logarithmic scale pre filtration versus linear scale post filtration. B): Non-pathogenic variants by CADD score categories, and molecular subtype as noted in key. C): Pathogenic variant scores by CADD score categories, and molecular subtype as noted in key. Note plots in C are identical, as all variants were still present post filtration.

We next tested as if seeking a single rare, disease-causing, pathogenic variant in each target dataset DNA. In keeping with recommendations for rare disease pathogenicity assignments,[3] 2×10⁻⁴ was set as the population allele frequency threshold to be included: Referencing allele frequencies in the 1000 Genome Project[17] and gnomAD[18] removed a mean of 2,476,589 and 2,483,377 variants per DNA respectively (Supplementary Figure 2A). Further rare-variant specific filters were then employed to exclude (i) the non-pathogenic variants present in the validation dataset DNAs (“white listed variants”), (ii) synonymous variants that were not in splice or regulatory regions, and (iii) variants where deleteriousness was not predicted based on low Combined Annotation-Dependent depletion (CADD) scores.[19] These filters removed a mean of 307,152 further variants. Overall, following application of all filters, a mean of 21,486 (0.44%) of variants were left per DNA, totaling 873,554 different variants across all DNAs. Comparing the pre and post filter variant numbers generated a Wilcoxon rank-sum p-value of <2.2e^-16. Irrespective of other filters applied, GROFFFY, and its individual components, significantly reduced the number of variants compared to the other tested filter sets (Supplementary Table 3B, Supplementary Figure 2B).

Validating GROFFFY

Finally, we demonstrated that applying GROFFFY did not delete key variants by examining the validation dataset. These 34 DNAs each harbored a single variant in one of the known genes for HHT, a vascular dyplasia caused by a single null allele most commonly in ENG, ACVRL1 or SMAD4.[13] With the exception of a mosaic case identified only by the clinical diagnostic pipelines (and hence not present in the unfiltered dataset[20]), all pathogenic variants were retained post filtration (Figure 2). identicality of plots as all variants were still present post filtration.