Bioconductor-style infrastructure for breeding and quantitative genomics: genotypes, phenotypes, pedigree and marker coordinates in one object, with the relationship matrices that quantitative genetics runs on.
Bioconductor has excellent containers for expression, variants and single-cell data. It has nothing for the shape of data a breeding programme actually produces: a marker matrix, a phenotype table, a pedigree that reaches back further than the genotyping did, and a marker map. This package provides that container, and the matrices computed from it.
Keeping genotypes, phenotypes and pedigree in separate objects means every
filter or reordering is a chance to desynchronise them. BreedingExperiment
extends RangedSummarizedExperiment, so:
GRanges, so a data set can be cut to a region,
overlapped with annotation, or handled with any Bioconductor tool that speaks
genomic ranges;From R-universe:
install.packages("BreedingExperiment",
repos = c("https://mqfarooqi1.r-universe.dev",
"https://bioconductor.org/packages/release/bioc",
"https://cloud.r-project.org"))
Or from GitHub:
# install.packages("BiocManager")
BiocManager::install("SummarizedExperiment")
remotes::install_github("mqfarooqi1/BreedingExperiment")
The package is under review for Bioconductor;
once accepted it will install with BiocManager::install("BreedingExperiment").
library(BreedingExperiment)
data(demoBreeding)
demoBreeding
#> class: BreedingExperiment
#> markers: 900 individuals genotyped: 120
#> assays(1): genotype
#> phenotypes(6): generation sex yield stature trueBV_yield trueBV_stature
#> sequences(10): chr1 chr2 chr3 chr4 chr5 chr6
#> pedigree: 180 individuals (120 genotyped, 60 not; 30 founders)
demoBreeding is a simulated five-generation population shaped like a real
programme: the two earliest generations were born before genotyping began, so
60 of the 180 individuals have no markers; one trait is recorded on females
only; a few sires dominate; and about one per cent of calls are missing. Use
simulateBreeding() to generate your own.
Your own data goes in the same way:
be <- BreedingExperiment(
genotypes = dosage_matrix, # markers x individuals, coded 0/1/2
phenotypes = phenotype_table, # one row per individual
pedigree = pedigree_table, # id, sire, dam
rowRanges = marker_granges) # marker positions
| Function | Matrix | Source |
|---|---|---|
Amatrix() |
Pedigree numerator relationships | Henderson (1976) |
Gmatrix() |
Genomic relationships from markers | VanRaden (2008) |
Dmatrix() |
Dominance relationships | Vitezica et al. (2013) |
Hmatrix() |
Single step: pedigree and markers combined | Legarra et al. (2009); Christensen & Lund (2010) |
A <- Amatrix(be) # every individual in the pedigree
G <- Gmatrix(be) # the genotyped ones
H <- Hmatrix(be) # all of them, using both sources of information
Real programmes genotype only part of the population. Using only G throws
away every ungenotyped relative; using only A throws away the markers.
Hmatrix() keeps both, so genotyped and ungenotyped individuals can be
analysed in one model — the relationship structure behind single-step genomic
prediction (ssGBLUP).
G is first tuned to the scale of A22 and then blended with it, which also
guarantees the result is invertible when markers are fewer than individuals.
be <- filterMarkers(be, min_maf = 0.05, min_call_rate = 0.9)
be <- imputeMarkers(be)
alleleFrequency(be) # freq, maf, call rate, heterozygosity per marker
inbreeding(be) # from the diagonal of A
The relationship matrices are checked against values that follow from theory,
not against a previous run: parent–offspring 0.5, grandparent 0.25, full sibs
0.5, the offspring of a full-sib mating inbred at F = 0.25, H positive
definite, and Hmatrix(blend = 1) reproducing A22 exactly.
O(n²) in memory, which is
fine for tens of thousands of individuals but not for a national evaluation
with millions. A sparse inverse is the obvious next step.MIT licensed. Muhammad Farooqi.