NeuroBooster Array: A Genome-Wide Genotyping Platform to Study Neurological Disorders Across Diverse Populations

NeuroBooster array design

The NBA contains a backbone of 1,914,934 genetic markers from the Infinium Global Diversity Array-8 v1.0 (GDA) complemented with custom content of 95,273 disease-associated variants involved in a wide range of neurological conditions. The custom content on the NBA was purposely designed to complement the base content of the GDA series of arrays from Illumina (Table 1).

The GDA is built on a high-density SNP global backbone optimized for cross-population imputation coverage of the genome [https://www.illumina.com/products]. The GDA contains a robust genome-wide scaffold designed to tag both common and low frequency variants (minor allele frequency > 1%) in global populations. This scaffold was designed through collaborations with the Consortium on Asthma among African-ancestry Populations in the Americas (CAAPA) [https://www.caapa-project.org/] and the Population Architecture using Genomics and Epidemiology (PAGE) study [https://www.pagestudy.org/]. For GDA design, more than 400,000 markers of exome content were gathered from 36,000 individuals of diverse ethnic groups, including African Americans, Hispanics, Pacific Islanders, East Asians, and individuals of mixed ancestry. The GDA includes content not found in the 1K Genomes Project, with more than 1000 whole-genome sequences of African ancestry and populations throughout the Americas, including the US, Caribbean, and Latin and South America. The GDA is designed to contain diverse exonic content from, Genome Aggregation Database version 2.1 (https://gnomad.broadinstitute.org/news/2018-10-gnomad-v2-1/, RRID:SCR_014964) including both cross-population and population-specific markers with either functionality or strong evidence for association with certain populations (Illumina Inc).

We designed the custom content for the NBA with four goals in mind. These goals included: 1) identify rare coding variants associated with disease of interest to researchers in the neurology field; 2) improve imputation of known GWAS loci for neurodegenerative diseases across populations of diverse ancestries; 3) generally improve imputation quality across populations of diverse ancestries to facilitate risk locus discovery; and 4) facilitate fine-mapping of risk loci.

To accomplish these aims, we aggregated custom content from a variety of neurodegenerative disease-relevant sources. In brief, building the array was carried out in five phases (described further below), consisting of sequential design components. The first design component was a systematic review as well as consultation with key opinion leaders in the field of neurology to identify putative rare variants associated with disease. The second design component consisted of annotating known GWAS risk loci and identifying multi-population tag SNPs to support improved imputation and locus tagging across ancestry groups in follow-up studies. Based on this information, the third component consisted of saturation of genotyping of risk loci from large GWAS regions. Next, we carried out probe design due diligence in collaboration with Illumina Inc. Finally, we used diverse samples and rare-variant-enriched cohorts to build a custom cluster file to improve genotype quality. The culmination of these phases of work allowed us to densely genotype and impute higher variant coverage into regions of interest related to neurological diseases.

Systematic review of publicly annotated databases and key opinion leaders

Whole-exome and whole-genome sequencing variant prioritization analyses were conducted in families and cohort studies. This includes some from diverse ancestries and population isolates, provided by International Parkinson’s Disease Genomics Consortium (IPDGC, https://pdgenetics.org) members and collaborators and used as the first layer of derived content nominated by key opinion leaders. We followed this up with systematic reviews of publicly annotated databases, including the Human Gene Mutation Database (HGMD) version 2019 (http://www.hgmd.org/, RRID:SCR_001621) ¹¹ and the Genomics England panel (accessed in May 2019) (https://panelapp.genomicsengland.co.uk/).

From HGMD, we extracted any nucleotide substitutions (missense/nonsense, regulatory, splicing) and also variants labeled as nucleotide deletions, insertions and indels for any disease or phenotype combinations suggested by key opinion leaders as potentially being implicated in the etiology of neurological diseases or conditions presenting with neurological complications. Our exhaustive list of search terms are highlighted in Supplementary Table 1. The database cross-referencing resulted in nearly 90,000 candidate variants initially. Due to array design and space concerns on the product, the variant list from this first portion of the design was reduced to 54,907 variants prior to Illumina due diligence (described later, Supplementary Table 5).

We also included an allocation of 1,809 polygenic risk score (PRS)-related variants for Parkinson’s disease risk nominated by PRSice-2 Version 2.2.13 ¹² (https://choishingwan.github.io/PRSice/, DOI: 10.5281/zenodo.3703335, RRID:SCR_017057). This includes variants that did not reach genome-wide significance but positively impacted PRS modeling efforts in previous reports. For more details on PRS construction, please refer to Nalls et al. 2019.

Identification of multi-ancestry tag SNPs from GWAS loci

Loci identified as genome-wide significant in large-scale GWAS meta-analyses of neurodegenerative diseases ^{13, 14, 15, 4, 16, 17, 18, 19} were the next target of our design process. Genome-wide significant loci from eight neurodegenerative disease-related conditions, including Alzheimer’s disease, Parkinson’s disease, amyotrophic lateral sclerosis, multiple sclerosis, progressive supranuclear palsy, corticobasal degeneration, multiple system atrophy, frontotemporal dementia and dementia with Lewy bodies were included with either multiple proxies for the top SNP at every locus, or technical replicates when proxies were not available. GWAS tagging SNPs per locus were defined as any SNP reaching a genome-wide significant p-value and correlated at r² < 0.3 with any other significant SNPs within 250 kilobases. This process was done for each disease of interest based on the combined European reference series from 1000 Genomes phase 1v5 data (https://www.internationalgenome.org/data-portal/data-collection/phase-1,RRID:SCR_0068 28) after filtering all variants for “PASS” annotation and minor allele count of at least three. This led to the successful inclusion of ∼180 GWAS-related variants (Supplementary Table 2).

To build multi-population references, the 1000 Genomes sequence data was again subsetted into ancestry reference “sub-populations”. These included conglomerate datasets from population labels as follows: MXL (Mexican Ancestry), CLM (Colombian), PEL (Peruvian), and PUR (Puerto Rican) were combined to form the AMR (American admixed) ancestry group; JPT (Japanese), CDX (Chinese Dai in Xishuangbanna), CHB (Han Chinese in Beijing), CHS (Han Chinese South), KHV (Kinh in Ho Chi Minh City, Vietnam) and CHD (Han Chinese in Beijing) were combined to form the EAS (East Asian) ancestry group; TSI (Toscani in Italy), IBS (Iberian Populations in Spain), GBR (British From England and Scotland), and CEU (Central European) were combined to form the EUR (European) ancestry group; PJL (Punjabi in Lahore, Pakistan), ITU (Indian Telugu in the U.K), STU (Sri Lankan Tamil in the UK), GIH (Gujarati Indians in Houston, Texas, USA), BEB (Bengali in Bangladesh) were combined to form the SAS (South Asians) ancestry group; GWD (Gambian in Western Division – Mandinka), MSL (Mende in Sierra Leone), ESN (Esan in Nigeria), YRI (Yoruba in Ibadan, Nigeria), LWK (Luhya in Webuye, Kenya), GWJ (Gambian Jola), GWF (Gambian Fula), and GWW (Gambian Wolof) were combined to form the AFR (African) ancestry group; ASW (African Ancestry in SW USA) and ACB (African Caribbean in Barbados) were combined to form the AAC (African admixed) ancestry group. Due to sample size issues for tagging analysis, a total of 115 samples from Finland were excluded from the references leaving us with 3,395 reference samples across all groups. Based on these reference samples and the 180 GWAS loci discussed earlier, the software TagIt version 1.0.8 https://bioinformaticshome.com/tools/imputation/descriptions/TagIt.html²⁰ was run using default settings per locus stratified by ancestry group to identify diverse tag SNP series (r2 > 0.3).

Saturation genotyping of known risk loci5

For each of the 180 GWAS variants of interest tagging unique loci, we identified upstream and downstream tag SNPs across all ancestry groups. For each region, all tag SNPs from diverse populations were included. To fill out content after the Illumina design triage described below in the due diligence content section, a surplus of EUR ancestry tag SNPs per locus were added. This phase of design included over 10,000 SNPs tagging diverse ancestral linkage structure at all loci. Regions tagged by the multi-population variants of interest covered over 38MB of the genome.

NeuroBooster array due diligence content

Illumina design scoring was undertaken at Illumina Inc to screen potential array content. This included a scoring of variants for quality and their ability to be included on arrays. Variants whose probes were problematic and for which low scoring could not reliably be built were excluded from the beadpool. Additionally, after clustering 2,793 diverse samples from the Global Parkinson’s Genetics Program ²¹ [https://gp2.org/training/] across multiple populations (Supplementary Table 3), we reviewed variants that were not accurately genotyped, with high rates of missingness across batches resulting in the exclusion of 49,590 variants (Supplementary Table 4). After iterative design scoring at Illumina Inc, custom content for 95,273 variants (and the additional 1,914,934 standard content variants) were annotated using the Variant Effect Predictor Release 100 and made available as part of the standard series of GP2 releases [https://useast.ensembl.org/info/docs/tools/vep/index.html].

NeuroBooster array data processing and custom clustering

Automated genotype data processing was conducted on GenoTools, a Python pipeline built for quality control and ancestry estimation of data. Additional details can be found at (https://github.com/GP2code/GenoTools). Quality control (QC) was performed according to standard protocols. Samples with a call rate below 95%, sex mismatches, or high heterozygosity (estimated by an |F| statistics of > 0.25) were excluded from analyses.

Further QC measures included the removal of SNPs with missingness above 5%, variants with significant deviations from Hardy-Weinberg Equilibrium (HWE P value < 1E-4), variants with non-random missingness and variants with missing data patterns (haplotype at P ≤ 1E-4 per ancestry). A total of 2,793 diverse case samples across six ancestry groups were clustered using the Illumina’s GenomeStudio Software v2.0.5 package (https://support.illumina.com/downloads/genomestudio-2-0.html, RRID:SCR_010973) (Figure 1, Supplementary Table 3).

• Download figure
• Open in new tab

Figure 1. Ancestry prediction clustering for samples genotyped on the NeuroBooster array

3 dimensional principal components analysis to group individuals based on their genetic makeup. A total of 2,793 samples from the Global Parkinson’s Genetics Program were included, including 2,373 Parkinson’s disease cases and 420 Gaucher disease cases. Each point represents a sample and the colors depict the ancestral background as shown in the color legend: EUR = orange, EAS = lemon green, AMR = yellow, AJ = lapis blue, AFR = teal blue, AAC = purple, SAS = magenta, CAS = dark pink, MDE = cerulean blue. EUR (Europe), EAS(East Asian), AMR (Latino/Admixed American), AJ (Ashkenazi Jewish), AAC (African-Admixed), AFR (African), SAS (South-Asian), CAS (Central-Asian), MDE (Middle East).

• Download figure
• Open in new tab

Figure 2. Brisbane plot showing the genomic density of SNPs on GP2 raw genotyped data generated on the NeuroBooster array

The clustering file is available for download via the GP2 Github repository (https://github.com/GP2code; DOI 10.5281/zenodo.7904832, release 5; https://gp2.org/). Clusters can be viewed using the GP2 data browser. Among the 2,793 samples, a total of 2,373 Parkinson’s disease patients were included in addition to 420 Gaucher disease cases to capture variants of interest in the GBA1 gene.

Clustering was done following Illumina Inc guidelines. Specifically, the robustness and reliability of genotype clustering was ensured by implementing a meticulous quality control protocol on the genotyping data. The protocol started with an evaluation of BeadChip performance in GenomeStudio using the Controls Dashboard. Any Beadchips with profiles suggestive of technical issues were promptly excluded from the analysis to negate the influence of technical inaccuracies on the results. Next, an assessment of sample call rates was conducted utilizing the Samples Table. To uphold the standard of genotyping data, samples presenting call rates lower than 0.98 were excluded from the subsequent analyses.

Following the assessment of call rates, an in-depth review of the Copy Number (CN) Metrics report was carried out to detect samples with outlier values for ‘LogRDev’ or ‘BAlleleDev’. Such samples, suggestive of potential genotyping anomalies, were excluded from further investigation to ensure a dataset free from major genotyping errors and anomalies. By employing this stringent sample inclusion/exclusion protocol, we assured that the subsequent genotype clustering was based on high-quality samples free from technical errors. This approach ensures that the genotyping data used for clustering is accurate, reliable, and most importantly, trustworthy.

Passing samples were then reclustered using the default Illumina Inc cluster algorithm built into GenomeStudio, with autosomes clustered together and sex chromosomes done separately.

NeuroBooster array genotyping protocol

An overview of the protocol can be visualized in Supplementary Figure 1. Briefly, for each sample a total of 200 ng of high-quality genomic DNA is amplified and enzymatically fragmented. The resulting fragments are alcohol-precipitated and resuspended in a buffer. Next, the fragmented DNA solution is hybridized to the NBA using a Tecan Freedom EVO liquid-handling robot (Tecan, Research Triangle Park, NC). After hybridization, automated allele-specific, enzymatic base extension and fluorophore staining are performed. The stained genotyping arrays are washed, sealed and vacuum-dried prior to scanning them on the Illumina iScan system. Raw data files are imported into GenomeStudio (version 2.0, Genotyping Module, Illumina Inc) using a custom-generated sample sheet, and genotypes are called using a GenCall threshold of 0.15.

NeuroBooster content overview

In total, the NBA contains 1,873,290 autosomal variants, 79,994 sex chromosomal variants (X and Y chromosomes), and 1,509 mitochondrial variants (Figure 2). The overlap between the NBA and NeuroChip is n= 126,220 variants. The NBA includes over 10,000 multi-ancestry GWAS locus tagging variants to facilitate imputation and analyses of these neurodegenerative disease-related GWAS loci across diverse populations. Detailed pathogenic inferences and annotated NBA functionally variant content for the array in total is provided through GP2 ²¹. Figure 3, Figure 4, and Supplementary Table 5 display the abundance of disease-associated variants per gene covered across the most prevalent neurodegenerative diseases based on initial systematic review and cross-referencing of the HGMD and Genomics England database content.

• Download figure
• Open in new tab

Figure 3. Overview of the number of Human Gene Mutation Database disease associated variants that are present on the NeuroBooster array for the most prevalent neurodegenerative diseases

AD = Alzheimer’s disease, ALS = amyotrophic lateral sclerosis, FTD = frontotemporal dementia, and PD = Parkinson’s disease

• Download figure
• Open in new tab

Figure 4A.

Overview of NeuroBooster array content by variant category

• Download figure
• Open in new tab

Figure 4B.

Overview of NeuroBooster array content by pathogenicity predictors

Genotyping accuracy

GenTrain scores were calculated for all NeuroBooster variants using GenomeStudio version 2.0. The GenTrain score is a statistical score based on the shapes of the different allelic clusters and their relative distance to each other (Illumina). Typically, GenTrain scores above 0.7 are considered high quality genotypes. The mean (and SD) for the array is 0.83 (0.091) (Supplementary Figure 2). Due to improved technology from Illumina Inc compared to NeuroChip or NeuroX arrays over the past decade as well as the iterative design process the array underwent in collaboration with Illumina Inc, we are able to rescue many variants not previously able to be genotyped (see Supplementary Figure 3 for example comparison).

Genetic variants in APOE, APP, GBA1, LRRK2, PARK2, DJ-1, PINK1, PRKN, PSEN1, SNCA, TREM2 and VPS35 contribute to AD, LBD and PD etiology among other neurodegenerative diseases. Genotyping of some of these genes remains complicated due to high GC content (APOE) or the presence of a pseudogene (GBA1). We comprehensively assessed variant content for tagging and imputed variants across multiple ancestries (additional details can be found in Supplementary Table 6, Supplementary Table 8).

Quality control metrics for imputation accuracy and allele frequencies intervals across diverse populations

Variants with a minor allele frequency (MAF) < 0.005 and Hardy-Weinberg equilibrium (HWE) p-value < 1E-5 were excluded before submission to the TOPMed Imputation server. The utilized TOPMed reference panel version, known as r2, encompasses genetic information from 97,256 reference samples and over 300 million genetic variants across the 22 autosomes and the X chromosome. As of October 2022, the TOPMed panel includes approximately 180,000 participants, with 29% of African, 19% of Latino/Admixed American ancestry, 8% of Asian ancestry, and 40% of European ancestry (https://topmed.nhlbi.nih.gov/). Further details about the TOPMed Study, Imputation Server, and Minimac Imputation can be accessed at https://topmed.nhlbi.nih.gov/, https://imputation.biodatacatalyst.nhlbi.nih.gov/#!(RRID:SCR_009292), https://genome.sph.umich.edu/w/index.php?title=Minimac:_Tutorial&oldid=14562, (RRID:SCR_009292). Following imputation, the resulting files underwent pruning based on a minor allele count (MAC) threshold of 10 and an imputation Rsq value of 0.3.

Imputation of autosomal variants was performed across 2,793 individuals across ten different ancestries (Supplementary Table 7). In its application to recent GP2 dataset releases [https://gp2.org/], imputation yields were solid. After filtering for imputation quality (r2 > 0.3 and minor allele count > 10) across ancestry groups, minimum variant yields using TOPMed imputation included 7.4 M to 19.6 M variants in total across ancestry groups. Lower frequency variants at MAF < 5% were also well imputed in this new data resource, with yields after filtering ranging from 0.8 M to 7.1 M across populations.

In addition to data imputation using TOPMed, we next assessed imputation accuracy performance by extracting the variant content available on the NBA from whole genome sequencing data from six ancestral populations available on 1000 Genomes (AAC, AFR, AMR, EAS, EUR, SAS). Optimal array design depends not only on tag SNP selection but also on empirical evaluation of imputation performance ²⁰. We compared imputation performance for neurodegenerative diseases related GWAS loci with and without the custom NBA content. For each of the population-specific GWAS scaffolds, imputation accuracy was assessed by MAF bins versus various imputation reference panels, including Haplotype Reference Consortium (HRC), the Consortium on Asthma among African-Ancestry Populations in the Americas (CAAPA) the Genome Asia Pilot (GAsP) (Supplementary Figure 4A-H).