Four haplotype terms, linkage structure and allele frequencies

Published

November 11, 2026

Code 
source("_setup.R")
spec <- read.csv("data/final_panel_spec.csv")
hap  <- read.csv("data/haplotype_frequencies.csv")
af   <- read.csv("data/panel_allele_frequencies.csv")

Why haplotype terms

The unfinished predecessor analysis simulated three PTX3 SNPs as independent multiplicative risks. They are in near-complete linkage disequilibrium and tag one h1/h2 haplotype, so independent simulation triple-counts a single signal. The same error would recur within CLEC7A and CD209 if their two SNPs each entered as additive terms. The panel therefore enters as four haplotype terms, with IL10 and IFNG as optional single SNPs. TLR9 rs352140 is dropped (one usable study of three).

Code 
spec |>
  mutate(variants = gsub(";", ", ", variants), study_ids = gsub(";", ", ", study_ids)) |>
  transmute(Term = term, Variants = variants, `Studies (union)` = n_studies_union,
            `Independent cohorts` = n_independent_cohorts, `Study IDs` = study_ids) |>
  gt()
Table 1: Final panel specification.
Term Variants Studies (union) Independent cohorts Study IDs
PTX3 haplotype rs2305619, rs3816527, rs1840680 7 7 S026, S029, S035, S045, S056, S077, S080
TLR4 haplotype rs4986790, rs4986791 5 5 S046, S056, S061, S066, S077
CLEC7A haplotype rs16910526, rs7309123 11 9 S023, S035, S045, S046, S056, S070, S071, S073, S074, S077, S080
CD209 haplotype rs4804800, rs7248637 5 4 S021, S045, S071, S077, S080
IL10 rs1800896 (optional) rs1800896 3 3 S003, S061, S077
IFNG rs2069705 (optional) rs2069705 3 3 S046, S072, S080

Herrero-Sanchez 2018 is excluded from the primary analysis because its endpoint is Candida or Aspergillus detection or colonisation and not proven or probable invasive mould disease. The cohort counts above include it, and excluding it takes PTX3 from 7 to 6 and CLEC7A from 9 to 8.

Linkage disequilibrium

Pairwise values in 1000 Genomes European samples.

Code 
hap |>
  distinct(gene, locus_1, locus_2, r2, d_prime) |>
  transmute(Gene = gene, `Locus 1` = locus_1, `Locus 2` = locus_2, `r²` = r2, `D′` = d_prime) |>
  gt() |>
  fmt_number(c(`r²`, `D′`), decimals = 3)
Table 2: Pairwise linkage between the two loci of each block.
Gene Locus 1 Locus 2 r² D′
PTX3 rs2305619 rs3816527 0.787 1.000
PTX3 rs2305619 rs1840680 0.996 1.000
PTX3 rs3816527 rs1840680 0.790 1.000
TLR4 rs4986790 rs4986791 0.945 0.981
CLEC7A rs16910526 rs7309123 0.104 1.000
CD209 rs4804800 rs7248637 0.413 0.671

For CLEC7A, D′ = 1 means no recombination between the two sites. The rare Y238X allele (rs16910526) therefore sits entirely on one rs7309123 background, and the low r² reflects the frequency difference and not independence. For CD209, D′ = 0.67 and all four haplotypes are present.

Haplotype phase check

Haplotype frequencies were derived from allele frequencies and D, assuming the minor alleles are in coupling (D positive). D′ near 1 fixes the size of D but not its sign. The reported r² and the allele frequencies fix |D| independently of that assumption, so the sign can be checked arithmetically: if even the largest possible negative D is smaller than the |D| that r² requires, negative phase is impossible.

Code 
phase <- hap |>
  group_by(gene, locus_1, locus_2) |>
  summarise(
    r2 = first(r2),
    p1 = sum(haplotype_freq_1000G_EUR[!allele_1 %in% c("major", "other")]),
    p2 = sum(haplotype_freq_1000G_EUR[!allele_2 %in% c("major", "other")]),
    .groups = "drop"
  ) |>
  mutate(
    required_D = sqrt(r2 * p1 * (1 - p1) * p2 * (1 - p2)),
    ceiling_if_negative = pmin(p1 * p2, (1 - p1) * (1 - p2)),
    negative_possible = required_D <= ceiling_if_negative + 1e-3
  )

phase |>
  transmute(Gene = gene, `Locus 1` = locus_1, `Locus 2` = locus_2,
            `Required |D|` = required_D, `Ceiling if D < 0` = ceiling_if_negative,
            `Negative phase possible` = ifelse(negative_possible, "yes", "no")) |>
  gt() |>
  fmt_number(c(`Required |D|`, `Ceiling if D < 0`), decimals = 4)
Table 3: Phase check. A negative D is possible only if the required |D| is no larger than the negative-phase ceiling.
Gene Locus 1 Locus 2 Required |D| Ceiling if D < 0 Negative phase possible
CD209 rs4804800 rs7248637 0.0678 0.0143 no
CLEC7A rs16910526 rs7309123 0.0408 0.0298 no
PTX3 rs2305619 rs1840680 0.2475 0.2078 no
PTX3 rs2305619 rs3816527 0.2162 0.1814 no
PTX3 rs3816527 rs1840680 0.2166 0.1810 no
TLR4 rs4986790 rs4986791 0.0524 0.0033 no

All six pairs are consistent only with coupling. The check assumes r² and D′ were taken from the 1000 Genomes source and not back-calculated from the assumed D. The phased chromosomes below confirm it directly.

Haplotype frequencies

Code 
hap |>
  mutate(pair = paste0(gene, ": ", locus_1, " x ", locus_2)) |>
  transmute(pair, `Allele 1` = allele_1, `Allele 2` = allele_2, Frequency = haplotype_freq_1000G_EUR) |>
  gt(groupname_col = "pair") |>
  fmt_number(Frequency, decimals = 4)
Table 4: Two-locus haplotype frequencies, 1000 Genomes European. ‘major’ and ‘other’ denote the non-minor allele.
Allele 1 Allele 2 Frequency
PTX3: rs2305619 x rs3816527
A C 0.3976
A other 0.0587
other C 0.0000
other other 0.5437
PTX3: rs2305619 x rs1840680
A A 0.4553
A other 0.0010
other A 0.0000
other other 0.5437
PTX3: rs3816527 x rs1840680
C A 0.3976
C other 0.0000
other A 0.0577
other other 0.5447
TLR4: rs4986790 x rs4986791
G T 0.0557
G other 0.0010
other T 0.0020
other other 0.9413
CLEC7A: rs16910526 x rs7309123
C G 0.0706
C major 0.0000
major G 0.3519
major major 0.5775
CD209: rs4804800 x rs7248637
G A 0.0821
G major 0.0422
major A 0.0332
major major 0.8425

That file holds pairwise frequencies derived from allele frequencies and D. The phased frequencies in the next section supersede it.

Phased haplotype frequencies

These come directly from phased 1000 Genomes phase 3 European genotypes (503 samples, 1,006 chromosomes, via the Ensembl REST API), then are rescaled to the gnomAD v4 non-Finnish-European allele frequencies that the project uses as primary. Rescaling is by iterative proportional fitting, which keeps the odds-ratio structure of the observed table, moves only the allele frequencies, and works for three loci. A haplotype that is absent in 1000G stays absent.

Code 
hp <- read.csv("data/haplotype_frequencies_phased.csv")
ldp <- read.csv("data/haplotype_ld_phased.csv")

PTX3, three loci

Alleles are in the order rs2305619, rs3816527, rs1840680.

Code 
hp |>
  filter(block == "PTX3_trio", n_chromosomes > 0) |>
  arrange(desc(freq_1000G_EUR)) |>
  transmute(Haplotype = haplotype, Chromosomes = n_chromosomes,
            `1000G EUR` = freq_1000G_EUR, SE = se_1000G_EUR, `gnomAD NFE (rescaled)` = freq_gnomAD_NFE_rescaled) |>
  gt() |>
  fmt_number(c(`1000G EUR`, SE, `gnomAD NFE (rescaled)`), decimals = 4)
Table 5: PTX3 three-locus haplotypes. SE is the binomial standard error on 1,006 chromosomes.
Haplotype Chromosomes 1000G EUR SE gnomAD NFE (rescaled)
G-A-G 547 0.5437 0.0157 0.5180
A-C-A 400 0.3976 0.0154 0.4242
A-A-A 58 0.0577 0.0073 0.0564
A-A-G 1 0.0010 0.0010 0.0014

There are three common haplotypes, not two. rs2305619 A and rs1840680 A mark the same class of chromosome (r² = 0.996), which is A-C-A plus A-A-A. rs3816527 C splits that class, and A-C-A is about 87% of it. The single A-A-G chromosome is too rare to estimate; its rescaled frequency is fixed by the margins alone.

CLEC7A, TLR4 and CD209

Code 
hp |>
  filter(block %in% c("CLEC7A", "TLR4", "CD209")) |>
  mutate(Block = paste0(block, ": ", gsub(";", " x ", loci))) |>
  transmute(Block, Haplotype = haplotype, Chromosomes = n_chromosomes,
            `1000G EUR` = freq_1000G_EUR, `gnomAD NFE (rescaled)` = freq_gnomAD_NFE_rescaled,
            `gnomAD NFE (D′ kept)` = freq_gnomAD_NFE_dprime) |>
  gt(groupname_col = "Block") |>
  fmt_number(c(`1000G EUR`, `gnomAD NFE (rescaled)`, `gnomAD NFE (D′ kept)`), decimals = 4)
Table 6: Two-locus haplotypes. ‘D′ kept’ rescales by preserving D′; the fitted values use iterative proportional fitting.
Haplotype Chromosomes 1000G EUR gnomAD NFE (rescaled) gnomAD NFE (D′ kept)
CLEC7A: rs16910526 x rs7309123
C-G 71 0.0706 0.0771 0.0771
A-G 354 0.3519 0.3492 0.3492
C-C 0 0.0000 0.0000 0.0000
A-C 581 0.5775 0.5737 0.5737
TLR4: rs4986790 x rs4986791
G-T 56 0.0557 0.0608 0.0603
A-T 2 0.0020 0.0034 0.0039
G-C 1 0.0010 0.0006 0.0011
A-C 947 0.9414 0.9352 0.9347
CD209: rs4804800 x rs7248637
G-A 84 0.0835 0.0778 0.0750
A-A 32 0.0318 0.0257 0.0285
G-G 41 0.0408 0.0473 0.0501
A-G 849 0.8439 0.8492 0.8464

The two rescaling methods agree closely. The largest difference is about 0.003, in CD209, where D′ is 0.69. For CLEC7A they are identical because D′ = 1.

Code 
ldp |>
  filter(!is.na(r2_1000G)) |>
  transmute(Block = block, Loci = gsub(";", " x ", loci),
            `r² (1000G)` = r2_1000G, `D′ (1000G)` = dprime_1000G,
            `r² (rescaled)` = r2_rescaled, `D′ (rescaled)` = dprime_rescaled) |>
  gt() |>
  fmt_number(where(is.numeric), decimals = 3)
Table 7: Linkage in the phased data and after rescaling to gnomAD NFE margins.
Block Loci r² (1000G) D′ (1000G) r² (rescaled) D′ (rescaled)
PTX3_a rs2305619 x rs3816527 0.787 1.000 0.792 1.000
PTX3_b rs2305619 x rs1840680 0.996 1.000 0.994 1.000
PTX3_c rs3816527 x rs1840680 0.790 1.000 0.796 1.000
CLEC7A rs16910526 x rs7309123 0.104 1.000 0.112 1.000
TLR4 rs4986790 x rs4986791 0.946 0.981 0.933 0.989
CD209 rs4804800 x rs7248637 0.431 0.685 0.414 0.716

The phased r² reproduces the earlier table for PTX3, TLR4 and CLEC7A. For CD209 it is 0.431 against 0.414 before. In every pair the minor alleles are in coupling.

NoteAssumptions in the rescaling

1000 Genomes European samples include a Finnish group (FIN), while gnomAD NFE excludes Finnish individuals. Rescaling assumes the linkage structure carries over and moves only the allele frequencies. Haplotypes seen on a handful of chromosomes are poorly estimated, and the table gives each standard error.

Allele frequencies

Code 
clean <- function(x) gsub(",", ", ", gsub("[{}\"]", "", x))
af |>
  transmute(rsID = rsid, Gene = gene, Position = chr_pos,
            `gnomAD NFE` = clean(gnomad_nfe_freqs), `1000G EUR` = clean(tg_eur_freqs)) |>
  gt()
Table 8: Allele frequencies for the 13 variants with three or more meta-analysable studies.
rsID Gene Position gnomAD NFE 1000G EUR
rs16910526 CLEC7A 12:10118488 A: 0.9229, C: 0.0771 C: 0.0706, A: 0.9294
rs7309123 CLEC7A 12:10119994 C: 0.5737, G: 0.4263 G: 0.4225, C: 0.5775
rs2305619 PTX3 3:157437072 A: 0.482, G: 0.518 A: 0.4563, G: 0.5437
rs4986790 TLR4 9:117713024 G: 0.0614, A: 0.9386 A: 0.9433, G: 0.0567
rs4986791 TLR4 9:117713324 C: 0.9358, T: 0.0642 T: 0.0577, C: 0.9423
rs3816527 PTX3 3:157437525 A: 0.5758, C: 0.4242 A: 0.6024, C: 0.3976
rs4804800 CD209 19:7740242 A: 0.8749, G: 0.1251 G: 0.1243, A: 0.8757
rs7248637 CD209 19:7742141 G: 0.8965, A: 0.1035 G: 0.8847, A: 0.1153
rs1840680 PTX3 3:157438240 A: 0.4806, G: 0.5194 A: 0.4553, G: 0.5447
rs5743708 TLR2 4:153705165 A: 0.0294, G: 0.9706 A: 0.0239, G: 0.9761
rs1800896 IL10 1:206773552 T: 0.5134, C: 0.4866 T: 0.5467, C: 0.4533
rs2069705 IFNG 12:68161231 A: 0.6624, G: 0.3379, C: 0.0 A: 0.673, G: 0.327
rs5743810 TLR6 4:38828729 A: 0.4217, G: 0.5783 A: 0.4085, G: 0.5915
Noters1800896 is reported on the opposite strand

gnomAD gives rs1800896 as T/C where the literature writes A/G. It is the same variant on the other strand, and it is the most-cited IL10 promoter SNP in this field, so a natural place for an orientation error.