Haplotype blocks of target loci

Groups a set of target loci into haplotype blocks by one of seven methods and reports the block assignment. Two calling conventions are supported.

Usage

hbr_haplotype_blocks(
  plink_prefix = NULL,
  units = NULL,
  cm_snp = NULL,
  chr_snp = NULL,
  method = c("recomb", "ld", "either", "fourgamete", "gabriel", "recomb+gabriel",
    "recomb+fourgamete"),
  flank_cM = 1,
  max_snps = 400L,
  max_span_cM = 5,
  min_gamete = 0.01,
  maf_min = 0.05,
  rmax = 0.01,
  r2 = 0.5,
  ld_max_cM = 3,
  r2_min = 0.25,
  gabriel_grid = 21L,
  carry = NULL,
  genetic_map = NULL,
  cM_per_Mb = NULL,
  verbose = FALSE
)

Arguments

plink_prefix: PLINK prefix, or an `hbr_run` object (in which case the run-level API is used). When the first positional argument is an `hbr_run`, this argument receives it.
units: Low-level API: data.frame with `chr` and `cm` per target locus, and `snp_col` when using the SNP-level methods. Alternatively, an `hbr_run` object (in which case the run-level API is used).
cm_snp, chr_snp: Low-level API: genetic positions and chromosomes of every SNP in the PLINK file.
method: One of the seven methods above.
flank_cM, max_snps, min_gamete: SNP-level methods only: window around the genes, SNPs kept per cluster, and gamete frequency that counts as recombination.
max_span_cM: Largest distance between the first and last locus of a block.
maf_min: Minimum minor-allele frequency.
rmax: Maximum Haldane recombination fraction between neighbours for `"recomb"`.
r2: Minimum carrier-indicator r-squared for `"ld"`.
ld_max_cM: Largest distance over which LD may join neighbours.
r2_min, gabriel_grid, carry: Reserved for future versions; kept for backward compatibility with the low-level API.
genetic_map, cM_per_Mb: Genetic positions for the run-level API: a data.frame with `chr`, `bp`, `cm` (for example [hbr_rice_genetic_map()]), or a constant rate. A `.bim` with an informative cM column is used when neither is supplied.
verbose: Print progress.

Value

For the low-level API, an integer vector of block IDs per row of `units`. For the run-level API, an object of class `hbr_blocks` with elements `blocks`, `units`, `lines`, `full`, `snp_clusters`, `settings` and `map_source`.

Details

**Low-level.** Supply `plink_prefix`, `units`, `cm_snp` and `chr_snp` directly. `units` must carry `chr`, `cm` and (for the SNP-level methods) `snp_col`; `cm_snp` and `chr_snp` give the genetic and physical position of every SNP in the PLINK file. The function returns an integer vector of block IDs per row of `units`. This is the legacy API and is fully backward compatible.

**Run-level.** Pass an `hbr_run` (from [hbr_read_results()]) as the first positional argument or as `plink_prefix`. The function reads the superior haplotypes, rebuilds the codes, replays the discovery gate to pick the best haplotype per gene, places each at its genetic map position, and returns a classed `hbr_blocks` object with the block table, the per-line carrier matrix, and (for SNP-level methods) the SNP windows used.

The seven block methods are:

`"recomb"`: Join neighbours whose Haldane recombination fraction is at most `rmax`. Describes the recombination expected in a new cross.
`"ld"`: Join neighbours whose carrier indicators have r-squared at least `r2` within `ld_max_cM`. Describes LD in the panel.
`"either"`: Union of `"recomb"` and `"ld"`.
`"fourgamete"`: Four-gamete rule on the SNPs around the genes (Wang et al. 2002; Haploview default).
`"gabriel"`: Gabriel-type D-prime confidence intervals (Gabriel et al. 2002).
`"recomb+gabriel"`: Union of `"recomb"` and `"gabriel"`, then capped at `max_span_cM`.
`"recomb+fourgamete"`: Union of `"recomb"` and `"fourgamete"`, then capped.

Blocks never span chromosomes. `"recomb"`, `"ld"` and `"either"` operate on the loci directly; the SNP-level methods read a window of `flank_cM` around each cluster of target genes.

Examples

# Legacy low-level API — no genotypes needed for the map methods
units <- data.frame(chr = c("1", "1", "1", "2", "2"),
                    cm  = c(10, 10.5, 50, 5, 5.2))
hbr_haplotype_blocks(method = "recomb", units = units)
#> [1] 1 1 2 3 3

# Run-level API on the bundled toy data
toy <- system.file("extdata", "hgsi_toy", package = "HaploBreedR")
if (nzchar(toy) && dir.exists(toy) &&
    requireNamespace("BEDMatrix", quietly = TRUE)) {
  run <- hbr_read_results(file.path(toy, "trait_YLD"))
  blk <- hbr_haplotype_blocks(
    run, plink_prefix = file.path(toy, "geno_toy"),
    method = "recomb")
  blk$blocks
}
#>   using bundled rice genetic map (Harushima 1998 / RGP); pass `genetic_map = "none"` to disable
#>   block chr start_cM end_cM n_target_haplotypes n_genes  genes traits
#> 1 HB001   2      0.0    0.0                   1       1 GENE04    YLD
#> 2 HB002   3      1.1    1.1                   1       1 GENE07    YLD
#>   panel_freq
#> 1     0.1667
#> 2     0.3333

Usage

Arguments

Value

Details

See also

Examples