Patients with previously diagnosed chronic granulomatous disease (CGD) (n = 45) and first-degree relatives (n = 3) followed at, or referred to, the National Institutes of Health were included in the study after being consented to IRB-approved protocols NCT00001355, NCT00404560, NCT00001467, or NCT00128973. Routine functional analysis included DHR assay for ROS production, immunoblot and/or flow cytometry to determine the presence of specific NADPH oxidase component proteins, and ddPCR for quantification of ΔGT as previously published (2). Illumina short-read whole exome (n = 21) or whole genome (n = 21) sequencing was performed at commercial laboratories (Johns Hopkins Genomics, Baltimore, MD; Baylor College of Medicine Human Genome Sequencing Center, Houston, TX). Oxford Nanopore long-read whole genome sequencing (n = 6) was performed at Johns Hopkins Genomics, Genetics Resource Core Facility, Baltimore, MD.

Non NCF1-CGD patients included 48 patients with X-linked CYBB mutations (including 4 skewed female patients), 4 bi-allelic CYBA, and 1 bi-allelic NCF2, recruited under the same protocols.

NCBI reference sequences for NCF1 (NM_00236.7; NP_000256.4), NCF1B (NR_003186), and NCF1C (NR_003187) were aligned using Sequencher (GeneCodes, Ann Arbor, MI).

To effectively map sequencing reads to this region we created an hg38 reference fasta file by masking the pseudogene sequence (chr7:73220639–73235945 and chr7:75156639-75172044) with “N”s using “bedtools maskfasta” (10, 11). Reads were mapped to our masked-reference using bwa-mem (https://bio-bwa.sourceforge.net/bwa.shtml). Variant calling was performed using GATK Best Practices (12), altering the ploidy to 6 to adjust for the tripling of reads in our region of interest. A custom script (https://github.com/niaid/NCF1_variant_calling), adapted from Almeida de Jesus, et al (13), was used to emit putative variant sites when 2 or more reads with alternate alleles were present, regardless of read balance. These sites were then fed to GATK for genotyping. The result of these methods is a list of predicted variants, most of which we were unable to discover using standard methods. We also calculate the alternate allele balance (AltAB) defined as the frequency of variant allele reads versus total read depth from NCF1 plus NCF1B/NCF1C at any given nucleotide. While AltAB is equivalent to variant allele frequency (VAF) used for high throughput sequencing data, given the 6 alleles present in the NCF1/NCF1B/NCF1C locus, we use AltAB to distinguish from normal, heterozygous loci.

Variants identified using publicly accessible 1000G short read sequence data (https://ftp.sra.ebi.ac.uk/vol1/run/) were parsed by ancestry to provide a reference dataset of masked NCF1 variant frequency in a healthy control population ( Supplementary Table 1 , 1000_Genomes_NCF1_variants.xlsx).

Primers were selected from outside the deleted region and used to amplify DNA from P5, P6, and the father of P6 (P6-F). Primers used were NCF1del 169F 5’- AAGATAAACCCAAACTAAGGGACATTCTACAAGG– 3’ and NCF1del 5960R 5’ – ATTTTATTTTGAGATGGAGTTTTGTCCTTGTTGC – 3’. Amplification was performed in 15 μl reactions using Platinum Taq HiFi (ThermoFisher) and cycling conditions 95°C 3 min, (95°C 20 sec, 62.3°C 10 sec, 68°C 30 sec) x35, 68°C 3 min; product was visualized on 1.2% agarose gel.

Patients were referred to the NIH for clinical evaluation or functional testing for suspected chronic granulomatous disease. All patients (n=48; 45 NCF1-CGD plus 3 NCF1-carriers) ( Table 1 ) were characterized by functional assays including flow cytometric DHR, immunoblots and/or flow cytometry for gp91^phox (CYBB), p22^phox (CYBA), p47^phox (NCF1), p67^phox (NCF2), and p40^phox (NCF4). All NCF1-CGD patients had diminished PMA-induced neutrophil superoxide production and DHR ( Supplementary Table 2 ) and undetectable p47^phox by immunoblot or flow cytometry. Frequency of ΔGT for each individual was performed using ddPCR ( Table 1 ). Genetic sequencing included short-read whole exome (WES) (n=21) or whole genome (WGS) (n=21) or long-read WGS (n=3 NCF1-CGD patients and 3 NCF1 carriers).

Laboratory testing among the NCF1-CGD cohort resembled previously reported patients (2, 14) with all patients displaying residual ROS after PMA stimulation and absent p47^phox protein (2). Among NCF1-CGD patients, ddPCR results demonstrated 37/45 (82.2%) had no intact GTGT alleles, 6/45 (13.3%) had 1 GTGT allele, while 2/45 (4.4%) had 2 intact GTGT alleles.

The presence of NCF1B/NCF1C (ΨNCF1) prevents unique alignment of sequencing reads. Therefore, we developed a bioinformatic pipeline to mask ΨNCF1, thereby aligning all sequence reads to NCF1 ( Figure 2A ). Utilizing this pipeline we then performed variant calling using standard methods. Adjusting the ploidy parameter from 2 to 6 allowed identification of variants from NCF1 reference sequence. For each variant, the ratio of variant calls to read depth was used to establish the alternate allele balance (AltAB) from NCF1B and NCF1C. Given that 3 unique autosomal genomic regions were included, AltAB increments approximate 1/6, with variants occurring in both ΨNCF1 loci having AltAB ~0.66 and those occurring in only NCF1B or NCF1C having AltAB ~0.33. We then used AltAB of ΔGT to compare 1000G, NCF1, and non-NCF1 CGD patients. Both 1000G and non-NCF1 patients had median AltAB of 0.66 indicating 4 copies of ΔGT from ΨNCF1 and two intact NCF1 GTGT alleles; by contrast, NCF1 patients had a median AltAB of 1.0 (indicating no remaining GTGT alleles) with only 11 individuals deviating from that ( Figure 2B ). Since we had previously performed ddPCR to determine ΔGT copies in these patients, we compared AltAB to ddPCR results ( Figure 2C ). Each of the 11 samples with ΔGT AltAB<1 (P1 through P11) was from patients with 1 (n=6) or 2 (n=2) copies of GTGT; also included were 3 NCF1 mutation carriers, with 2 (n=2) or 3 (n=1) copies of GTGT. Plotting ddPCR-determined GTGT copy number versus AltAB demonstrated concordance between the two methods (R² = 0.980), validating the use of AltAB to detect the most common NCF1 mutation.

Using NCBI reference sequences for the 3 loci we established a variant table, allowing assignment of variants to NCF1B, NCF1C or both (ΨNCF1) ( Supplementary Table 3 ). We then examined the AltAB distribution of these variants among 1000G. ΨNCF1 variants ( Supplementary Figure 1 ) mostly displayed a median AltAB (mAltAB) of 0.66, consistent with 4/6 copies of ΔGT, although the distribution was broad with 3/6 or 5/6 copies not uncommon and several other variants displaying 2/6 copies. For 9 loci, mAltAB was 1.0 indicating incorporation of the variant in NCF1 or incorrect NCF1 reference sequence. It is noteworthy that, despite a normal distribution and mAltAB ~0.66, 12 individuals have AltAB=1.0 for c.269G>A encoding p.R90H ( Supplementary Figure 1 , blue).

Examining NCF1C variants, in which mAltAB should be ~0.33, 3/14 variants have mAltAB>0.33 with multimodal distribution; an additional 4 variants have mAltAB<0.33, 2 of which have bimodal distributions of 1/6 or 2/6 while the remaining 2 variants have mAltAB=0.26 (Supplementary Figure 2 ). A similar NCF1B variant analysis reveals 2/14 variants with mAltAB~0.66, 4/14 variants with multimodal distributions and mAltAB>0.33 ( Supplementary Figure 3 ). Collectively, this analysis identifies a subset of variants having frequencies consistent with their presence in one or both pseudogenes (noted in Supplementary Table 3 ) available for further analyses of the locus.

Using these data, we next examined distribution of AltAB across NCF1 among NCF1 patients, 1000G, and non-NCF1 CGD patients. It was previously suggested that inclusion of ΔGT arises from meiotic crossover between NCF1 and one of the pseudogenes (7). Supporting this, NCF1-CGD patients have higher AltAB values throughout the NCF1 locus compared to 1000G controls and non-NCF1 patients ( Figure 3A ) indicating a higher proportion of ΨNCF1 variants among patients. NCF1 patients displayed a clear skewing of AltAB with mAltAB=0.66 while both 1000G and non-NCF1 patients had mAltAB of 0.52 and 0.49 respectively. This suggests larger incorporation of ΨNCF1 than solely the ΔGT-containing exon 2 ( Figure 2B ). To explore this, we plotted AltAB across 4 variants specific to NCF1C or NCF1B. NCF1 patients display increased AltAB across NCF1C variants compared to 1000G or non-NCF1 patients ( Figure 3B ); by contrast, there is no significant increase in AltAB across NCF1B SNPs ( Figure 3C ), although some individual patients have increased NCF1B SNP AltAB. Together, these data suggest NCF1 replacement with NCF1C occurs more frequently than with NCF1B.

These data suggested full or partial replacement of NCF1 with NCF1C in the majority of ΔGT NCF1 alleles, which we sought to confirm. Using variants demonstrated to have normal distribution among 1000G, we normalized AltAB at each site to 1000G mAltAB and plotted the normalized value across the locus. For patients with no intact GTGT, there would be 6 copies of ΔGT rather than 4, giving a normalized AltAB of 1.5. In the majority of patients (27/37 homozygous ΔGT, 73.0%), we observe normalized AltAB of 1.5 at ΨNCF1 SNPs across the locus indicating full replacement of NCF1 with a pseudogene ( Figures 4A, B ); additionally, 2/4 heterozygous ΔGT patients replaced one copy of NCF1 with a pseudogene resulting in normalized AltAB of ~1.25. Using the NCF1B or NCF1C specific variants, we observed cases of full NCF1 replacement by NCF1C as demonstrated by a normalized NCF1C AltAB of 2 across the locus indicating 4 copies of NCF1C and 2 copies of NCF1B (P43, Figure 4A ). Alternative arrangements are present as well, including replacement of NCF1 in the setting of 5 copies of NCF1C and only 1 NCF1B (P46, Figure 4B ).

Not all patients replaced the full gene, there was a frequent crossover point in intron 5 ( Figures 4C, D red arrows) corresponding to consecutive AluJr and AluSx1 repeats present in all three loci. Among NCF1-CGD patients with short-read sequencing (n=42), we could identify recurrent breakpoints between exons 2-4 (5/42, 11.9%) ( Figure 4D , blue arrow) and within intron 5 (10/42, 23.8%); 4 patients within these two groups exhibited both breakpoints ( Figure 4D ), one patient each had breaks between introns 4–5 and within intron 8; one patient was indeterminate due to lack of informative SNPs ( Table 1 ). In P21, NCF1 was replaced by NCF1C on one allele and only the 5’ portion of NCF1B on the other ( Figure 4C ). Lastly, P2 carried a missense mutation, (c.839T>C, p.L280P) on one allele, and ΔGT on the other. In this patient the mAltAB was ~1.25 with increased NCF1C in the 5’ region, and normalization of AltAB to 1 after intron 5 ( Figure 4D ). We validated this approach in two NCF1-CGD siblings, in whom the pattern of normalized AltAB for NCF1B, NCF1C, and ΨNCF1 was similar, as would be expected given the same parental alleles ( Supplementary Figure 4 ).

Having confirmed NCF1 replacement by NCF1C or NCF1B in ΔGT patients we looked for the inverse allele in a control population – that is, replacement of NCF1C or NCF1B by NCF1. Examination of ΔGT frequency in 1000G revealed 24 individuals (24/2504; 1%) with ΔGT AltAB<0.4, indicating fewer than 4 copies of ΔGT. In contrast to the rearrangements seen in CGD patients, these individuals had had NCF1 replacement of a pseudogene. One individual had AltAB = 0.38 for ΔGT, indicating only 2 remaining copies of a pseudogene. SNP analysis across NCF1B and NCF1C revealed a total loss of NCF1B but maintenance of NCF1C SNPs, suggesting replacement of NCF1B by NCF1 ( Figure 4E ). Another individual displayed apparent loss of at least one copy of NCF1C ( Figure 4F ), although exact determination was difficult due to the lack of NCF1C specific variants with normal distribution in the 5’ region.

Having demonstrated AltAB as a valid metric for ΔGT, we used AltAB to screen for non-ΔGT mutations among patients with available short-read sequencing data. Similar to examination of other genes, we first established a variation reference. Analyzing 1000G data with our pipeline identified all variants with an AltAB frequency > 0.08. We determined the number of individuals with each variant as well as minimum, maximum, mean, median and standard deviation of AltAB for each variant; those data were further split by ancestry ( Supplementary Table 1 ). Using the variants identified in 1000G, we examined NCF1-CGD patients with 1 or 2 intact GTGT alleles indicating non-ΔGT variations. We considered variants rare or unique among 1000G and the human variation database, Genome Aggregation Database (gnomAD) (15), and predicted deleterious by the bioinformatic algorithm, Combined Annotation Dependent Depletion (16), using a threshold of CADD>20, identifying 6 variants ( Figure 5 , upper; Table 1 ). Each patient with 1 intact GTGT copy had a single mutation, while P1, with 2 GTGT alleles, carried 2 separate mutations; each variant had a frequency approximating 1/6, suggesting that patients were heterozygous for the variants. P1 had two variants, c.72+3G>A at the start of intron 1 and c.125G>A (p.R42Q) within exon 2, both previously reported (14). P2 also carried a previously reported mutation, c.574G>A (p.G192S), which occurs at the last base of exon 6 leading to impaired splicing of the exon (14). P3, P4, and P5 carried novel variants: c.892_905+11del spanning the exon 9/intron 9 boundary, c.839T>C (p.L280P), and c.107C>T (p.S36L), respectively. Three of the variants, c.107C>T, c.125G>A, and c.574G>A, are present in gnomAD at low frequency (~1/100,000) however only c.107C>T occurs in 1000G (n=2). While it is not always possible to phase variants with short-read sequencing, c.107C>T is within 35 bases of ΔGT; using the Integrative Genomics Viewer (17) (IGV), a platform allowing visualization of high throughput sequencing alignments, we determined this variant was allelic with ΔGT ( Supplementary Figure 5 ), indicating the variant was not causative of CGD.

Since we were unable to identify a second mutation in the P5/P6 family, we performed Oxford Nanopore long-read sequencing. With an average read length of 25kb, we confidently spanned the full NCF1 locus and were able to phase variants across the region. In P6, the niece of P5, we identified a 5394 bp deletion fully encompassing exon 1 and part of intron 1 (chr7:74,770,624-74,776,017), ( Supplementary Figure 6 ). Examination of sequencing reads in IGV revealed the large deletion was not allelic with ΔGT ( Supplementary Figure 6 ). Using primers spanning the deleted region, the presence of the 5394 bp deletion was confirmed in P5 and the father of P6 (P6-F, an obligate carrier) ( Supplementary Figure 7A ). Additionally, the c.107C>T variant in P5, allelic with ΔGT in P5, was not present in P6 suggesting a commonly inherited large deletion within this family and different ΔGT alleles present in P5 and P6. ( Supplementary Figure 7B ). It is notable that the deceased brother of P5 (uncle of P6) was diagnosed with CGD after Mycobacterium fortuitum infection at age 27 (18).

Two additional NCF1-CGD patients were diagnosed using long read sequencing. P7 carried one variant occurring at the last base of exon 6, (c.574G>A; p.G192S/splice) which was non-allelic with ΔGT. P8, with two copies of intact GTGT, had three damaging variants identified, c.72+3G>T at the start of intron 1, c.500A>C encoding p.Y167S, and c.579G>A encoding p.W193*. Examination of reads in IGV revealed that c.500A>C was allelic with the c.579G>A mutation but not the splice mutation ( Supplementary Figure 8 ). Lastly, three presumed carriers of non-ΔGT mutations were sequenced. In two cases, a single, known pathogenic mutation was identified which was allelic with intact GTGT; P9, mother of P1, carried one copy of c.125G>A (p.R42Q), while P10 carried c.574G>A (p.G192S/splice). No mutation was identified in P11, the father of an NCF1-CGD patient. It is possible the patient carries a de novo mutation on the paternal allele or P11 is a germline mosaic.