Linkage disequilibrium analysis

Introduction

library(GeneTrackR)

GeneTrackR separates linkage disequilibrium analysis into two stages. compute_ld_block() calculates pairwise LD from VCF genotypes and stores the complete result in an LDTrack; plot_ld_block() then adds a triangular heatmap and optionally a gene/variant region track to the same object. Keeping calculation and plotting separate makes the numerical LD table, matrix, region definition, sample settings, and final figure available for downstream analysis.

The bundled demo now contains three deterministic LD scenarios:

varLD01-varLD12: twelve SNPs forming the primary LD-gradient example. varLD01-varLD06 are a perfect high-LD core, while varLD07-varLD12 progressively decorrelate from that core so the heatmap contains high, intermediate, and low r2 values;
varLow01-varLow04: four variants designed to have little linkage, with most pairwise r2 = 0 and a maximum of approximately 0.111;
varPair01-varPair02: exactly two variants, retained to test the single-diamond LD plot.

The added gradient variants are outside the GeneA gene body, so the canonical GeneA haplotype example still contains 11 variants and four balanced haplotypes.

This module follows nine steps:

prepare genotype-rich VCF data and gene annotation;
compute the primary twelve-variant LD-gradient example;
inspect the LDTrack object and pairwise statistics;
inspect the LD matrix and quantify high/intermediate/low LD pairs;
draw the triangular LD heatmap using the default Reds palette;
add gene structure, natural-variant markers, and connector lines;
compare the LD-gradient, low-LD, and two-variant regions;
control samples, variant types, LD methods, and retained genotype matrices;
interpret LD results and pass a selected region to haplotype analysis.

Package documentation: https://renscq.github.io/GeneTrackR/.

Installation

After Bioconductor acceptance, install GeneTrackR with BiocManager:

if (!requireNamespace("BiocManager", quietly = TRUE)) {
  install.packages("BiocManager")
}
BiocManager::install("GeneTrackR")
library(GeneTrackR)

The development version remains available from the GeneTrackR GitHub repository.

Step 1. Prepare VCF genotypes and annotation

Locate and read the same demo VCF and GenePred annotation used in the variant and haplotype workflows:

vcf_file <- system.file(
  "extdata",
  "gtr_demo_variants.vcf",
  package = "GeneTrackR"
)
gp_file <- system.file(
  "extdata",
  "gtr_demo.genePredExt",
  package = "GeneTrackR"
)
vcf <- read_vcf(
  vcf_file,
  mode = "memory",
  keep_genotype = TRUE,
  verbose = FALSE,
  progress = FALSE
)
gp <- read_genepred(
  gp_file,
  format = "genePredExt",
  verbose = FALSE,
  progress = FALSE
)

LD calculation requires VCF genotype columns. If the VCF was previously read with keep_genotype = FALSE, read it again with keep_genotype = TRUE before calling compute_ld_block().

For large bgzip-compressed and indexed VCF files, a lazy VariantTrack can be used, but lazy LD queries require a complete chrom, start, and end interval so only the requested region is retrieved.

Step 2. Compute the primary twelve-variant LD-gradient example

The main demo interval starts inside GeneA and extends into the short downstream interval before GeneB. It contains 12 designed LD SNPs plus two GeneA indels. Use variant_type = "snp" so the LD example contains exactly varLD01-varLD12:

ld_gradient <- compute_ld_block(
  vcf,
  chrom = "chr1",
  start = 12342620,
  end = 12355500,
  variant_type = "snp",
  method = "r2",
  min_pair_samples = 3,
  verbose = FALSE
)
ld_gradient
#> <LDTrack>
#>   region    : chr1:12342620-12355500
#>   method    : r2
#>   variants  : 12
#>   pairs     : 66
#>   samples   : 36
#>   figure    : no

The expected result is:

12 retained SNPs;
66 unique variant pairs (12 * 11 / 2);
36 available samples per pair;
a perfect six-variant core (varLD01-varLD06, r2 = 1 within the core);
a progressive decline in LD toward varLD12, including weak pairs with r2 < 0.2 and pairs reaching r2 = 0.

Check the retained variants before interpreting the heatmap:

ld_gradient$variants
#>     variant_index  chrom      pos variant_id    ref    alt variant_type
#>             <int> <char>    <int>     <char> <char> <char>       <char>
#>  1:             1   chr1 12342620    varLD01      G      A          SNP
#>  2:             2   chr1 12342710    varLD02      A      T          SNP
#>  3:             3   chr1 12342805    varLD03      C      T          SNP
#>  4:             4   chr1 12342920    varLD04      G      C          SNP
#>  5:             5   chr1 12343040    varLD05      T      G          SNP
#>  6:             6   chr1 12343180    varLD06      A      C          SNP
#>  7:             7   chr1 12352500    varLD07      G      T          SNP
#>  8:             8   chr1 12353100    varLD08      A      C          SNP
#>  9:             9   chr1 12353700    varLD09      C      A          SNP
#> 10:            10   chr1 12354300    varLD10      T      G          SNP
#> 11:            11   chr1 12354900    varLD11      G      C          SNP
#> 12:            12   chr1 12355500    varLD12      A      G          SNP

The first six variants remain inside GeneA and use the same genotype pattern required by the GeneA haplotype demo. The six added variants lie downstream of the GeneA gene body, so they increase LD-plot density without changing the 11-variant GeneA haplotype definition.

Step 3. Inspect the LDTrack object and pairwise statistics

LDTrack keeps the complete state of an LD calculation:

names(ld_gradient)
#> [1] "data"     "matrix"   "variants" "region"   "genotype" "figure"   "meta"

The main components are:

Component | Content | Typical use |

— | — | — |

data | one row per unique variant pair | inspect r2, Dprime, distance, sample number, and allele frequencies |

matrix | symmetric matrix of the selected LD metric | matrix-based inspection or export |

variants | ordered metadata for retained variants | confirm positions and IDs |

region | chromosome and queried boundaries | record the analyzed interval |

genotype | optional variant-by-sample dosage matrix | detailed QC when explicitly retained |

figure | figure created by plot_ld_block() | plotting/export |

meta | method, sample count, ploidy, and calculation notes | reproduce the analysis settings |

Inspect the pairwise table directly:

ld_gradient$data[, .(
  variant_i,
  variant_j,
  distance_bp,
  n_samples,
  r,
  r2,
  Dprime,
  ld
)]
#>     variant_i variant_j distance_bp n_samples         r         r2    Dprime
#>        <char>    <char>       <int>     <int>     <num>      <num>     <num>
#>  1:   varLD01   varLD02          90        36 1.0000000 1.00000000 1.0000000
#>  2:   varLD01   varLD03         185        36 1.0000000 1.00000000 1.0000000
#>  3:   varLD01   varLD04         300        36 1.0000000 1.00000000 1.0000000
#>  4:   varLD01   varLD05         420        36 1.0000000 1.00000000 1.0000000
#>  5:   varLD01   varLD06         560        36 1.0000000 1.00000000 1.0000000
#>  6:   varLD01   varLD07        9880        36 0.8888889 0.79012346 1.0000000
#>  7:   varLD01   varLD08       10480        36 0.7777778 0.60493827 1.0000000
#>  8:   varLD01   varLD09       11080        36 0.6666667 0.44444444 1.0000000
#>  9:   varLD01   varLD10       11680        36 0.5007734 0.25077399 1.0000000
#> 10:   varLD01   varLD11       12280        36 0.3333333 0.11111111 0.6666667
#> 11:   varLD01   varLD12       12880        36 0.0000000 0.00000000 0.0000000
#> 12:   varLD02   varLD03          95        36 1.0000000 1.00000000 1.0000000
#> 13:   varLD02   varLD04         210        36 1.0000000 1.00000000 1.0000000
#> 14:   varLD02   varLD05         330        36 1.0000000 1.00000000 1.0000000
#> 15:   varLD02   varLD06         470        36 1.0000000 1.00000000 1.0000000
#> 16:   varLD02   varLD07        9790        36 0.8888889 0.79012346 1.0000000
#> 17:   varLD02   varLD08       10390        36 0.7777778 0.60493827 1.0000000
#> 18:   varLD02   varLD09       10990        36 0.6666667 0.44444444 1.0000000
#> 19:   varLD02   varLD10       11590        36 0.5007734 0.25077399 1.0000000
#> 20:   varLD02   varLD11       12190        36 0.3333333 0.11111111 0.6666667
#> 21:   varLD02   varLD12       12790        36 0.0000000 0.00000000 0.0000000
#> 22:   varLD03   varLD04         115        36 1.0000000 1.00000000 1.0000000
#> 23:   varLD03   varLD05         235        36 1.0000000 1.00000000 1.0000000
#> 24:   varLD03   varLD06         375        36 1.0000000 1.00000000 1.0000000
#> 25:   varLD03   varLD07        9695        36 0.8888889 0.79012346 1.0000000
#> 26:   varLD03   varLD08       10295        36 0.7777778 0.60493827 1.0000000
#> 27:   varLD03   varLD09       10895        36 0.6666667 0.44444444 1.0000000
#> 28:   varLD03   varLD10       11495        36 0.5007734 0.25077399 1.0000000
#> 29:   varLD03   varLD11       12095        36 0.3333333 0.11111111 0.6666667
#> 30:   varLD03   varLD12       12695        36 0.0000000 0.00000000 0.0000000
#> 31:   varLD04   varLD05         120        36 1.0000000 1.00000000 1.0000000
#> 32:   varLD04   varLD06         260        36 1.0000000 1.00000000 1.0000000
#> 33:   varLD04   varLD07        9580        36 0.8888889 0.79012346 1.0000000
#> 34:   varLD04   varLD08       10180        36 0.7777778 0.60493827 1.0000000
#> 35:   varLD04   varLD09       10780        36 0.6666667 0.44444444 1.0000000
#> 36:   varLD04   varLD10       11380        36 0.5007734 0.25077399 1.0000000
#> 37:   varLD04   varLD11       11980        36 0.3333333 0.11111111 0.6666667
#> 38:   varLD04   varLD12       12580        36 0.0000000 0.00000000 0.0000000
#> 39:   varLD05   varLD06         140        36 1.0000000 1.00000000 1.0000000
#> 40:   varLD05   varLD07        9460        36 0.8888889 0.79012346 1.0000000
#> 41:   varLD05   varLD08       10060        36 0.7777778 0.60493827 1.0000000
#> 42:   varLD05   varLD09       10660        36 0.6666667 0.44444444 1.0000000
#> 43:   varLD05   varLD10       11260        36 0.5007734 0.25077399 1.0000000
#> 44:   varLD05   varLD11       11860        36 0.3333333 0.11111111 0.6666667
#> 45:   varLD05   varLD12       12460        36 0.0000000 0.00000000 0.0000000
#> 46:   varLD06   varLD07        9320        36 0.8888889 0.79012346 1.0000000
#> 47:   varLD06   varLD08        9920        36 0.7777778 0.60493827 1.0000000
#> 48:   varLD06   varLD09       10520        36 0.6666667 0.44444444 1.0000000
#> 49:   varLD06   varLD10       11120        36 0.5007734 0.25077399 1.0000000
#> 50:   varLD06   varLD11       11720        36 0.3333333 0.11111111 0.6666667
#> 51:   varLD06   varLD12       12320        36 0.0000000 0.00000000 0.0000000
#> 52:   varLD07   varLD08         600        36 0.8888889 0.79012346 1.0000000
#> 53:   varLD07   varLD09        1200        36 0.7777778 0.60493827 1.0000000
#> 54:   varLD07   varLD10        1800        36 0.6120564 0.37461300 1.0000000
#> 55:   varLD07   varLD11        2400        36 0.4444444 0.19753086 0.8888889
#> 56:   varLD07   varLD12        3000        36 0.1111111 0.01234568 0.2222222
#> 57:   varLD08   varLD09         600        36 0.8888889 0.79012346 1.0000000
#> 58:   varLD08   varLD10        1200        36 0.7233393 0.52321981 1.0000000
#> 59:   varLD08   varLD11        1800        36 0.5555556 0.30864198 1.0000000
#> 60:   varLD08   varLD12        2400        36 0.2222222 0.04938272 0.4444444
#> 61:   varLD09   varLD10         600        36 0.8346223 0.69659443 1.0000000
#> 62:   varLD09   varLD11        1200        36 0.6666667 0.44444444 1.0000000
#> 63:   varLD09   varLD12        1800        36 0.3333333 0.11111111 0.6666667
#> 64:   varLD10   varLD11         600        36 0.8346223 0.69659443 1.0000000
#> 65:   varLD10   varLD12        1200        36 0.5007734 0.25077399 1.0000000
#> 66:   varLD11   varLD12         600        36 0.6666667 0.44444444 1.0000000
#>     variant_i variant_j distance_bp n_samples         r         r2    Dprime
#>        <char>    <char>       <int>     <int>     <num>      <num>     <num>
#>             ld
#>          <num>
#>  1: 1.00000000
#>  2: 1.00000000
#>  3: 1.00000000
#>  4: 1.00000000
#>  5: 1.00000000
#>  6: 0.79012346
#>  7: 0.60493827
#>  8: 0.44444444
#>  9: 0.25077399
#> 10: 0.11111111
#> 11: 0.00000000
#> 12: 1.00000000
#> 13: 1.00000000
#> 14: 1.00000000
#> 15: 1.00000000
#> 16: 0.79012346
#> 17: 0.60493827
#> 18: 0.44444444
#> 19: 0.25077399
#> 20: 0.11111111
#> 21: 0.00000000
#> 22: 1.00000000
#> 23: 1.00000000
#> 24: 1.00000000
#> 25: 0.79012346
#> 26: 0.60493827
#> 27: 0.44444444
#> 28: 0.25077399
#> 29: 0.11111111
#> 30: 0.00000000
#> 31: 1.00000000
#> 32: 1.00000000
#> 33: 0.79012346
#> 34: 0.60493827
#> 35: 0.44444444
#> 36: 0.25077399
#> 37: 0.11111111
#> 38: 0.00000000
#> 39: 1.00000000
#> 40: 0.79012346
#> 41: 0.60493827
#> 42: 0.44444444
#> 43: 0.25077399
#> 44: 0.11111111
#> 45: 0.00000000
#> 46: 0.79012346
#> 47: 0.60493827
#> 48: 0.44444444
#> 49: 0.25077399
#> 50: 0.11111111
#> 51: 0.00000000
#> 52: 0.79012346
#> 53: 0.60493827
#> 54: 0.37461300
#> 55: 0.19753086
#> 56: 0.01234568
#> 57: 0.79012346
#> 58: 0.52321981
#> 59: 0.30864198
#> 60: 0.04938272
#> 61: 0.69659443
#> 62: 0.44444444
#> 63: 0.11111111
#> 64: 0.69659443
#> 65: 0.25077399
#> 66: 0.44444444
#>             ld
#>          <num>

Important columns are:

distance_bp: physical distance between the two variants;
n_samples: paired non-missing sample number used for that comparison;
r: dosage correlation;
r2: squared dosage correlation;
Dprime: dosage-based D-prime approximation;
ld: the metric selected by method; for method = "r2", ld is identical to r2.

The region and settings used for calculation are preserved in the object:

ld_gradient$region
#> $chrom
#> [1] "chr1"
#> 
#> $start
#> [1] 12342620
#> 
#> $end
#> [1] 12355500
ld_gradient$meta
#> $method
#> [1] "r2"
#> 
#> $variant_type
#> [1] "snp"
#> 
#> $sample_n
#> [1] 36
#> 
#> $min_pair_samples
#> [1] 3
#> 
#> $ploidy
#> [1] 2
#> 
#> $dosage_definition
#> [1] "number of non-reference allele copies parsed from GT"
#> 
#> $dprime_note
#> [1] "Dprime is approximated from genotype dosage for unphased VCFs. Use phased haplotypes for exact haplotype-based Dprime."

Step 4. Inspect the LD matrix and quantify the LD gradient

The symmetric LD matrix uses variant IDs as row and column names and contains 1 on the diagonal:

round(ld_gradient$matrix, 3)
#>         varLD01 varLD02 varLD03 varLD04 varLD05 varLD06 varLD07 varLD08 varLD09
#> varLD01   1.000   1.000   1.000   1.000   1.000   1.000   0.790   0.605   0.444
#> varLD02   1.000   1.000   1.000   1.000   1.000   1.000   0.790   0.605   0.444
#> varLD03   1.000   1.000   1.000   1.000   1.000   1.000   0.790   0.605   0.444
#> varLD04   1.000   1.000   1.000   1.000   1.000   1.000   0.790   0.605   0.444
#> varLD05   1.000   1.000   1.000   1.000   1.000   1.000   0.790   0.605   0.444
#> varLD06   1.000   1.000   1.000   1.000   1.000   1.000   0.790   0.605   0.444
#> varLD07   0.790   0.790   0.790   0.790   0.790   0.790   1.000   0.790   0.605
#> varLD08   0.605   0.605   0.605   0.605   0.605   0.605   0.790   1.000   0.790
#> varLD09   0.444   0.444   0.444   0.444   0.444   0.444   0.605   0.790   1.000
#> varLD10   0.251   0.251   0.251   0.251   0.251   0.251   0.375   0.523   0.697
#> varLD11   0.111   0.111   0.111   0.111   0.111   0.111   0.198   0.309   0.444
#> varLD12   0.000   0.000   0.000   0.000   0.000   0.000   0.012   0.049   0.111
#>         varLD10 varLD11 varLD12
#> varLD01   0.251   0.111   0.000
#> varLD02   0.251   0.111   0.000
#> varLD03   0.251   0.111   0.000
#> varLD04   0.251   0.111   0.000
#> varLD05   0.251   0.111   0.000
#> varLD06   0.251   0.111   0.000
#> varLD07   0.375   0.198   0.012
#> varLD08   0.523   0.309   0.049
#> varLD09   0.697   0.444   0.111
#> varLD10   1.000   0.697   0.251
#> varLD11   0.697   1.000   0.444
#> varLD12   0.251   0.444   1.000

Unlike the previous all-identical six-variant example, the new matrix is intentionally heterogeneous. The upper-left core is strongly linked, whereas LD decreases as progressively decorrelated downstream variants are compared with the core.

Summarize the distribution directly:

range(ld_gradient$data$r2, na.rm = TRUE)
#> [1] 0 1
summary(ld_gradient$data$r2)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>  0.0000  0.2508  0.4838  0.5228  0.7901  1.0000

Define descriptive groups explicitly when useful:

ld_pair_summary <- ld_gradient$data[, .(
  high_ld = sum(r2 >= 0.8, na.rm = TRUE),
  intermediate_ld = sum(r2 >= 0.2 & r2 < 0.8, na.rm = TRUE),
  low_ld = sum(r2 < 0.2, na.rm = TRUE)
)]
ld_pair_summary
#>    high_ld intermediate_ld low_ld
#>      <int>           <int>  <int>
#> 1:      15              35     16

For the deterministic demo, 15 pairs have r2 >= 0.8 and 16 pairs have r2 < 0.2. The remaining pairs form the intermediate gradient.

GeneTrackR calculates pairwise LD but does not automatically call biological LD-block boundaries from a fixed threshold. A rule such as r2 >= 0.8 is therefore a user-defined interpretation threshold, not an implicit block-calling rule inside compute_ld_block().

Step 5. Draw the triangular LD heatmap

plot_ld_block() now defaults to the sequential Reds palette, which is well suited to an r2 scale from 0 to 1:

ld_gradient <- plot_ld_block(
  ld_gradient,
  label_by = "variant_id",
  show_variant_labels = TRUE,
  show_region = FALSE,
  return_object = TRUE,
  font = 10
)
ld_gradient$figure

The default color interpretation is straightforward: pale cells indicate low LD and progressively darker red cells indicate stronger LD.

The palette can still be changed explicitly when needed. Any sequential or diverging RColorBrewer palette accepted by GeneTrackR can be supplied; for example:

ld_gradient_blue <- plot_ld_block(
  ld_gradient,
  color_palette = "Blues",
  show_variant_labels = FALSE,
  return_object = TRUE,
  font = 10
)
ld_gradient_blue$figure

plot_ld_block() returns the updated LDTrack by default. The calculation results remain unchanged and the generated ggplot/patchwork figure is stored in $figure.

Step 6. Add gene structure and natural-variant markers

For locus interpretation, add the compact genomic region track above the triangular heatmap:

ld_gradient_region <- plot_ld_block(
  ld_gradient,
  show_region = TRUE,
  annotation = gp,
  connect_region = TRUE,
  show_variant_marker = TRUE,
  variant_marker_size = 3,
  label_by = "variant_id",
  show_variant_labels = TRUE,
  region_height = 1.25,
  connector_height = 0.35,
  return_object = TRUE,
  font = 10
)
ld_gradient_region$figure

The combined plot contains three aligned components when connect_region = TRUE:

compact gene/variant region track;
connector lines from genomic variant positions to LD columns;
triangular LD heatmap.

With 12 variants, the region panel is now dense enough to show the relationship between physical position and the LD gradient. Variant triangle markers can be hidden with show_variant_marker = FALSE, or resized with variant_marker_size.

Step 7. Compare the LD-gradient, low-LD, and two-variant regions

The dedicated low-LD interval remains useful as an extreme comparison:

ld_low <- compute_ld_block(
  vcf,
  chrom = "chr2",
  start = 1999000,
  end = 2008500,
  method = "r2",
  verbose = FALSE
)
ld_low$data[, .(
  variant_i,
  variant_j,
  distance_bp,
  r2
)]
#>    variant_i variant_j distance_bp        r2
#>       <char>    <char>       <int>     <num>
#> 1:  varLow01  varLow02        1300 0.0000000
#> 2:  varLow01  varLow03        3000 0.0000000
#> 3:  varLow01  varLow04        5300 0.1111111
#> 4:  varLow02  varLow03        1700 0.0000000
#> 5:  varLow02  varLow04        4000 0.0000000
#> 6:  varLow03  varLow04        2300 0.0000000

The six pairwise values are designed so that most are zero and the maximum is approximately 1 / 9:

range(ld_low$data$r2, na.rm = TRUE)
#> [1] 0.0000000 0.1111111
max(ld_low$data$r2, na.rm = TRUE)
#> [1] 0.1111111

Plot it using the same default Reds scale:

ld_low <- plot_ld_block(
  ld_low,
  label_by = "variant_id",
  return_object = TRUE,
  font = 10
)
ld_low$figure

The two-variant demo region contains exactly one pairwise comparison:

ld_pair <- compute_ld_block(
  vcf,
  chrom = "chr2",
  start = 16995001,
  end = 17006000,
  method = "r2",
  verbose = FALSE
)
ld_pair$data
#>    variant_i variant_j index_i index_j    pos_i    pos_j distance_bp n_samples
#>       <char>    <char>   <int>   <int>    <int>    <int>       <int>     <int>
#> 1: varPair01 varPair02       1       2 16996000 17004500        8500        36
#>        r    r2     D Dprime_signed Dprime   p_i   p_j    ld method
#>    <num> <num> <num>         <num>  <num> <num> <num> <num> <char>
#> 1:     1     1   0.5             1      1   0.5   0.5     1     r2

plot_ld_block() draws this case as one centered diamond:

ld_pair <- plot_ld_block(
  ld_pair,
  label_by = "variant_id",
  return_object = TRUE,
  font = 10
)
ld_pair$figure

When only the figure is needed for a temporary plotting call, compatibility mode can return it directly:

pair_figure <- plot_ld_block(
  ld_pair,
  return_object = FALSE,
  font = 10
)
pair_figure

Step 8. Control samples, variant types, LD methods, and genotype retention

Select samples

LD is population dependent. Use samples when LD should be estimated in a defined population subset:

subset_samples <- sprintf("S%02d", 1:24)
ld_subset <- compute_ld_block(
  vcf,
  chrom = "chr1",
  start = 12342620,
  end = 12355500,
  variant_type = "snp",
  samples = subset_samples,
  min_pair_samples = 10,
  method = "r2",
  verbose = FALSE
)
ld_subset$meta$sample_n
#> [1] 24

Select variant types

variant_type controls the records retained before pairwise calculation. It is important in the primary demo because the genomic interval also contains GeneA indels that are not part of the 12-SNP LD truth set:

ld_snp <- compute_ld_block(
  vcf,
  chrom = "chr1",
  start = 12342620,
  end = 12355500,
  variant_type = "snp",
  method = "r2",
  verbose = FALSE
)

Available choices are both, snp, and ind. At least two retained variants are required.

Retain the dosage matrix when needed

The genotype dosage matrix is omitted by default to reduce object size. Retain it only when detailed genotype QC or downstream matrix analysis is needed:

ld_with_gt <- compute_ld_block(
  vcf,
  chrom = "chr1",
  start = 12342620,
  end = 12355500,
  variant_type = "snp",
  method = "r2",
  keep_genotype_matrix = TRUE,
  verbose = FALSE
)
dim(ld_with_gt$genotype)
#> [1] 12 36
ld_with_gt$genotype[, 1:6, drop = FALSE]
#>         S01 S02 S03 S04 S05 S06
#> varLD01   0   0   0   0   0   0
#> varLD02   0   0   0   0   0   0
#> varLD03   0   0   0   0   0   0
#> varLD04   0   0   0   0   0   0
#> varLD05   0   0   0   0   0   0
#> varLD06   0   0   0   0   0   0
#> varLD07   2   0   0   0   0   0
#> varLD08   2   2   0   0   0   0
#> varLD09   2   2   2   0   0   0
#> varLD10   2   2   2   2   2   0
#> varLD11   2   2   2   2   2   2
#> varLD12   2   2   2   2   2   2

For the demo, the retained genotype matrix is 12 x 36.

Use D-prime cautiously

method = "Dprime" is also available:

ld_dprime <- compute_ld_block(
  vcf,
  chrom = "chr2",
  start = 1999000,
  end = 2008500,
  method = "Dprime",
  verbose = FALSE
)
ld_dprime$data[, .(
  variant_i,
  variant_j,
  r2,
  Dprime,
  ld
)]
#>    variant_i variant_j        r2    Dprime        ld
#>       <char>    <char>     <num>     <num>     <num>
#> 1:  varLow01  varLow02 0.0000000 0.0000000 0.0000000
#> 2:  varLow01  varLow03 0.0000000 0.0000000 0.0000000
#> 3:  varLow01  varLow04 0.1111111 0.6666667 0.6666667
#> 4:  varLow02  varLow03 0.0000000 0.0000000 0.0000000
#> 5:  varLow02  varLow04 0.0000000 0.0000000 0.0000000
#> 6:  varLow03  varLow04 0.0000000 0.0000000 0.0000000

The current Dprime implementation is a dosage-based approximation for unphased VCF genotypes. It should not be described as an exact haplotype-based D-prime estimate. For ordinary unphased genotype VCF workflows, r2 is the recommended primary GeneTrackR LD metric.

Step 9. Interpret LD and pass the selected region downstream

LD measures statistical association between variants; it does not by itself identify a causal variant or prove that all variants in a visually strong region should be treated as one biological unit. Interpret the pattern together with the focal variant, physical distance, gene structure, population, variant quality, and phenotype/association evidence.

Once an interval has been selected biologically, the same region can be passed to the haplotype workflow. To reproduce the 12-SNP LD definition, retain SNPs only:

hap_ld_region <- hap_region_variant(
  vcf,
  chrom = "chr1",
  start = 12342620,
  end = 12355500,
  variant_type = "SNP",
  genotype_mode = "string"
)
#> [GeneTrackR] Retrieved variants: 12.
hap_ld_region
#> <HapVariant>
#>   region    : chr1:12342620-12355500
#>   variants  : 12
#>   samples   : 36
#>   haplotypes: 14

This keeps the analysis chain explicit:


VCF genotypes

    -> pairwise LD calculation

    -> inspect the LD gradient and regional gene context

    -> define the biological interval of interest

    -> construct regional haplotypes

    -> phenotype/refinement analysis if appropriate

For real GWAS loci, it is usually preferable to define the LD query around a lead SNP or candidate interval first, inspect the pairwise pattern, and only then decide which linked variants or genes should move into the haplotype analysis.

Session information

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] GeneTrackR_0.99.2 rmarkdown_2.32   
#> 
#> loaded via a namespace (and not attached):
#>  [1] gtable_0.3.6         jsonlite_2.0.0       crayon_1.5.3        
#>  [4] compiler_4.6.1       zip_3.0.2            Rcpp_1.1.2          
#>  [7] Biostrings_2.81.9    bitops_1.1-0         Rsamtools_2.29.0    
#> [10] GenomicRanges_1.65.4 parallel_4.6.1       jquerylib_0.1.4     
#> [13] IRanges_2.47.5       Seqinfo_1.3.2        scales_1.4.0        
#> [16] BiocParallel_1.47.0  yaml_2.3.12          fastmap_1.2.0       
#> [19] XVector_0.53.0       ggplot2_4.0.3        R6_2.6.1            
#> [22] patchwork_1.3.2      labeling_0.4.3       generics_0.1.4      
#> [25] openxlsx_4.2.9       knitr_1.52           BiocGenerics_0.59.12
#> [28] maketools_1.3.2      bslib_0.12.0         RColorBrewer_1.1-3  
#> [31] rlang_1.3.0          stringi_1.8.9        cachem_1.1.0        
#> [34] xfun_0.61            sass_0.4.10          sys_3.4.3           
#> [37] S7_0.2.2             otel_0.2.0           cli_3.6.6           
#> [40] withr_3.0.3          digest_0.6.39        grid_4.6.1          
#> [43] lifecycle_1.0.5      S4Vectors_0.51.10    vctrs_0.7.3         
#> [46] evaluate_1.0.5       glue_1.8.1           data.table_1.18.6.1 
#> [49] farver_2.1.2         codetools_0.2-20     buildtools_1.0.0    
#> [52] stats4_4.6.1         tools_4.6.1          htmltools_0.5.9