Package 'simplePHENOTYPES'

Description

Subset the individuals of a Population

Usage

## S3 method for class 'Population'
x[i]

Arguments

Value

A Population with the selected individuals.

Examples

data("SNP55K_maize282_maf04")
pop <- as_population(SNP55K_maize282_maf04)
pop[1:2]

Description

The additive relationship matrix $A$ of a Population's individuals, built from its recorded pedigree (see parentage()) by the tabular method (Emik & Terrill 1949) on coefficients of relationship (Wright 1922): processing individuals parents-first, $A_{ij} = (A_{i,m(j)} + A_{i,f(j)}) / 2$ and $A_{jj} = 1 + F_j$ with $F_j = A_{m(j) f(j)} / 2$. A self has both parents equal, so $F = (1 + F_P) / 2$; a doubled haploid is fully inbred, $A_{jj} = 2$, and relates to others as its parent does. Ancestors not in pop still enter through the pedigree.

Usage

a_matrix(pop, ids = NULL, founder_f = 0)

Arguments

Details

Founders are the base of the pedigree: unrelated to each other, and by default non-inbred (founder_f = 0, $A_{jj} = 1$), which is right for outbred (random-mating) founders. Inbred founders – e.g. the bundled maize inbred lines – are not: their inbreeding coefficient enters the base as $A_{jj} = 1 + F_j$ and their relationships with descendants, so pass founder_f = 1 for fully inbred lines (a founder $P$ crossed to another inbred line then has $A_{P,F1} = 1$ and the two founders' $F1$ sibs $A = 1$, against $0.5$ with the default). Leaving the default on inbred founders halves those relationships and disagrees systematically with the genomic matrix of g_matrix(), whose founder diagonal is close to 2 for inbred lines. Downstream individuals follow from the recursion; selfing and doubled-haploid rules are not altered by founder_f. Note that a Population in which the same individual appears twice (e.g. c(pop, pop[1])) lists that individual under two ids with a relationship equal to its own diagonal; a_matrix() does not refuse it.

Value

A symmetric matrix with dimnames ids.

References

Wright S (1922) Coefficients of inbreeding and relationship. The American Naturalist 56:330–338. doi:10.1086/279872

Emik LO, Terrill CE (1949) Systematic procedures for calculating inbreeding coefficients. Journal of Heredity 40:51–55. doi:10.1093/oxfordjournals.jhered.a105986

See Also

parentage(), predict_ebv(), g_matrix()

Examples

data("SNP55K_maize282_maf04")
pop <- as_population(SNP55K_maize282_maf04, individuals = 1:2)
f1 <- cross(pop[1], pop[2], n = 2, seed = 1)
s1 <- selfcross(f1[1], n = 2, seed = 2)
a_matrix(c(f1, s1))
# the bundled lines are inbred: declare it (A between a line and its F1 is 1)
a_matrix(c(pop, f1), founder_f = 1)

Description

Add an additive variance component

Usage

additive(
  sim,
  prop = NULL,
  n_qtn = NULL,
  qtn = NULL,
  effect = NULL,
  phase = c("coupling", "repulsion"),
  dist = "geometric",
  orthogonal = FALSE,
  a = NULL,
  d = NULL
)

Arguments

Details

The additive value is the -1/0/1 dosage weighted by the QTN effects, centered, and scaled to prop of the phenotypic variance. This is a simulation convention, not Fisher's average-effect decomposition; for an additive-only model prop equals the narrow-sense h2 under Hardy-Weinberg.

Combining this layer with dominance() on the same loci (the default) makes the realized genetic variance differ from prop_A + prop_D by ⁠2Cov(c_A, c_D)⁠ (with one additive and one dominance layer; several layers of one type add their mutual covariance to Var(c_A) or Var(c_D)), a structural term that depends on the allele frequencies and on which allele is coded +1 (see simulate_phenotype()); the object reports the realized Var(A), Var(D) and 2Cov(A,D), and Var(c_A), Var(c_D) and 2Cov(c_A,c_D), in ⁠$ad_report⁠ and in print(). Use orthogonal = TRUE for a Fisher-orthogonal additive/dominance model.

Value

the updated phenotype_sim.

Orthogonal genotypic model (orthogonal = TRUE)

Instead of the dosage coding above, the layer builds each locus's genotypic value from an additive effect a and a dominance deviation d (value -a / +d / +a for gene content 0 / 1 / 2), and the whole genotypic value is scaled to prop of the phenotypic variance. Its additive and dominance variances then emerge from a, d and the allele frequencies rather than from separate layer proportions: the additive (breeding-value) component uses the average effect of substitution $\alpha_j = a_j + d_j(1 - 2 p_j)$, and the dominance deviation is the realized residual $D = g - A$. The additive and dominance components are orthogonal ($Cov(A, D) = 0$) in expectation under random mating, which puts each locus in Hardy-Weinberg proportions; $\alpha$ above is the one-generation transmitting-ability form. This does not require linkage equilibrium – between-locus LD is compatible with orthogonality in this no-epistasis model – but per-locus HWE alone is not sufficient under arbitrary nonrandom multilocus genotype association (e.g. selection or population structure), which can correlate one locus's gene content with another's heterozygosity. A finite or structured sample therefore carries a (usually small) $Cov(A, D)$ that breeders customarily ignore; this implementation does not, and the variance budget reports the realized shares $Var(A)/Var(g)$ and $Var(D)/Var(g)$ on separate additive and dominance rows, plus an add_dom_cov row for $2\,Cov(A, D)/Var(g)$, so the three sum to the layer prop (the add_dom_cov row is $\approx 0$ for a large random-mating or F2 sample and non-zero otherwise). The realized narrow-sense heritability follows from the degree of dominance d / abs(a) (magnitude 1 = complete dominance, > 1 = overdominance; abs because a may be negative – which the variance-partition grammar cannot express, because per-component scaling washes a global degree out), and the breeding value used by select_ind() / optimum_contribution() picks up the dominance-induced average effect automatically. Because the effects are fixed, this mode does not draw a separate dominance() layer, is incompatible with vary_qtn, and cannot be used with the correlation-controlling "pleiotropy" (multi-trait) or "ld" architectures (their effects come from the correlated draw; ⁠qtn =⁠ alone is accepted).

References

Fisher RA (1918) The correlation between relatives on the supposition of Mendelian inheritance. Trans. R. Soc. Edinb. 52:399-433; Falconer DS (1985) A note on Fisher's 'average effect' and 'average excess'. Genet. Res. 46(3):337-347 (the average effect of substitution $\alpha = a + d(q - p)$ and its dependence on random-mating / Hardy-Weinberg proportions); Falconer DS, Mackay TFC (1996) Introduction to Quantitative Genetics, 4th ed. Longman, Harlow; Lynch M, Walsh B (1998) Genetics and Analysis of Quantitative Traits. Sinauer, Sunderland, MA (additive / dominance decomposition).

See Also

dominance(), epistasis(), vqtl() for the other layers.

Examples

data("SNP55K_maize282_maf04")

# Additive trait: 3 QTNs explaining half the phenotypic variance,
# so h2 = 0.5 and the remaining 0.5 is residual.
ph <- simulate_phenotype(SNP55K_maize282_maf04, seed = 1) |>
  additive(prop = 0.5, n_qtn = 3)
ph

head(phenotypes_long(ph))

# A major QTN with a set share of the variance: stack two additive layers,
# one locus explaining 20% of V_P and 50 equal-effect loci another 30%
g <- SNP55K_maize282_maf04
major <- simulate_phenotype(g, seed = 1) |>
  additive(prop = 0.2, qtn = g$snp[100]) |>
  additive(prop = 0.3, n_qtn = 50, effect = rep(1, 50))
tb <- qtn_table(major)
head(tb[order(-tb$var_explained), c("snp", "var_explained")], 3)

# Effect sizes follow a geometric series by default. `effect` sets its base,
# so QTN effects here are 0.5, 0.25, 0.125, ...
simulate_phenotype(SNP55K_maize282_maf04, seed = 1) |>
  additive(prop = 0.5, n_qtn = 3, effect = 0.5)

# Or give the effects explicitly, one per QTN.
simulate_phenotype(SNP55K_maize282_maf04, seed = 1) |>
  additive(prop = 0.5, n_qtn = 3, effect = c(0.6, 0.3, 0.1))

# Fix the additive QTNs yourself by marker name.
simulate_phenotype(SNP55K_maize282_maf04, seed = 1) |>
  additive(prop = 0.5, qtn = c("ss196442916", "ss196439337", "ss196480535"))

# Orthogonal genotypic model: give per-locus a and d; the additive and
# dominance variances emerge (see the var_budget). Dominance needs
# heterozygotes, so simulate on a segregating F2.
pop <- as_population(SNP55K_maize282_maf04, individuals = 1:20)
f2  <- selfcross(cross(pop[1], pop[2], n = 1, seed = 1), n = 40, seed = 2)
og  <- simulate_phenotype(f2, h2 = 0.6, seed = 3) |>
  additive(orthogonal = TRUE, a = 0.5, d = 0.5, n_qtn = 8)
og$var_budget

Description

Scores each individual by $\sum_j \mathrm{dosage}_{ij}\,\mathrm{effect}_j$ over a given, frozen architecture (loci qtn and their effects), using the -1/0/1 dosages directly with no per-population centring or rescaling.

Usage

additive_value(x, qtn, effect)

Arguments

Details

This is the accessor to use when comparing populations across generations – e.g. to show a selection response as the mean additive value climbs. It differs from genetic_values(), which re-centres and re-scales every layer to its target prop on whatever population is being scored: that makes genetic_values() read roughly mean 0 / variance prop at every generation, so it cannot express a cross-generational trend. Here the loci and effects are fixed by the caller, so the scale is stable and the numbers are comparable from one generation to the next. Effects are on the same -1/0/1 dosage scale as an additive() layer's effect (and as qtn_table()'s effect column).

Value

a named numeric vector of additive values, one per individual.

See Also

genetic_values(), dosages(), qtn_table(), select_ind().

Examples

data("SNP55K_maize282_maf04")
pop <- as_population(SNP55K_maize282_maf04, individuals = 1:20)
# Score on a fixed 3-locus architecture; the scale does not depend on the set.
av  <- additive_value(pop, qtn = c(1, 5, 9), effect = c(0.5, -1, 2))
head(av)

Description

Converts character SNP genotypes to numeric coding, the format every simulation function in the package expects. HapMap, VCF, GDS and PLINK bed/ped input are recognized automatically, from a file path or from an object already in memory.

Usage

as_numeric(x, ...)

Arguments

Details

When to_file = TRUE and no file_name is given, the default output name is ⁠<label>_numeric.txt⁠, where the label is the name of the object passed as x (for a file input, the input file's name without its extension). A label that is not already a portable file name (other characters replaced, or longer than 100 characters) is sanitized and gets a short hash of the original label appended; this makes accidental clashes unlikely but does not make names unique. Different inputs can still map to the same default name (labels that differ only in case on a case-insensitive file system, files such as same.hmp.txt and same.txt in one directory, or two inline objects of the same class and size, whose generic label is built from class and dimensions). A default-named file that already exists is overwritten, with a warning naming the file; an explicit file_name is used as given, without a warning. Pass file_name to choose the output path.

Value

The genotypes in numeric format: five metadata columns (snp, allele, chr, pos, cm), plus the optional counted column when counted_column = TRUE, followed by one column per individual. Returned as a data frame when to_r = TRUE, otherwise written to file.

Which allele becomes 1

By default (method = "frequency") the most frequent allele is the reference: it is coded 1, the minor allele -1, and the heterozygote 0. The choice is made per marker from the observed genotype calls, so it does not depend on how the input file ordered its alleles – with one exception, the tie rule below.

Tie rule (minor-allele frequency exactly 0.5). When the two homozygote classes are equally frequent at a marker (for example 50 of 100 doubled-haploid lines carry each allele), no allele is more frequent, and the coding follows the allele order instead: the first-listed allele ("allele 1") is coded as the reference (1, or 2 under code_as = "012"), and the second-listed as the minor allele (this is the strict > comparison in the numericalization kernel: the alleles are swapped only when the second allele is strictly more frequent). "First-listed" means the first allele of the input's allele pair: REF for VCF, the first allele of the alleles column for HapMap, A1 (as SNPRelate lists it) for PLINK bed/ped, and the first allele of snp.allele in a GDS file. For inputs that declare no allele pair (nucleotide tables and Illumina FinalReport) the alphabetically first allele is allele 1. Every reader applies the same rule, so the same VCF gives identical dosage columns from a file path and from a data frame, tied markers included.

The alleles column of HapMap input is checked, not trusted. If the calls at a marker contain alleles that are not in its declared pair (a stale or opposite-strand alleles column), a warning names the affected rows and the marker is oriented from the calls; the returned allele column then shows the alleles actually observed.

The returned allele column lists the two alleles as "allele1/allele2" in the order above (for a marker whose second allele was never observed, only allele 1). The column itself records the raw orientation, not which allele was coded 1.

Which allele is counted. The allele coded 1 (or 2 under code_as = "012") at each marker is recorded in the "counted_allele" attribute of the returned data frame, a character vector with one allele per marker (NA where it is not known, and absent under model = "Dom", which counts no allele). No dosage value depends on it: it only lets cross(), c.Population() and the selection schemes see whether two separately converted panels count the same allele. Two panels that both carry it and count different alleles at a marker are refused when crossed or pooled (e.g. an all-AA panel and an all-GG panel would otherwise both encode +1 and cross to an F1 of +1 instead of the heterozygote 0). The attribute is a convenience of the R object: a numeric text file cannot carry it, and [ keeps it unchanged when it subsets or reorders the rows of the data frame (so it no longer lines up with them). The durable form is the optional counted column: with counted_column = TRUE the result, and the file written, gain a character column counted immediately after cm (before the individuals), holding the same allele per marker (NA where it is not known). It travels with its marker through files and row subsetting, as_population(), filter_geno() and as_numeric() itself accept it, and when present it is the record used, so two panels read back from files, or subsetted, are still checked against each other. The default (counted_column = FALSE) writes exactly the five metadata columns and changes no value. counted_column = TRUE is an error under model = "Dom" (no allele is counted) and for numeric input that carries neither the column nor the attribute; numeric input that already has the column keeps it. A numeric data frame without the column (or with a numeric sixth column, which is an individual) has no record and gets the older label-only check: convert panels that will be crossed jointly, or with the same ref_allele, to be safe.

Calls are matched case-insensitively. A marker with more than two alleles is set to missing with a warning that counts the affected markers; missing marker IDs in VCF / PLINK input are filled with chr:pos, and duplicated marker IDs or sample names are an error (downstream functions index by name).

VCF calls. Only complete, biallelic diploid genotype calls are used, from a file path and from an in-memory VCF table alike: haploid calls ("0", "1"), partially missing calls ("./1") and calls that name a second ALT allele ("0/2", "2/2") are set to missing, and one warning per kind counts them. The marker's orientation is then chosen from the remaining calls.

Set method = "reference" with ref_allele to fix the reference allele yourself instead — for example to keep the coding consistent with an external panel, where the major allele in your sample may not be the one you want coded 1.

code_as = "012" switches to counting the reference allele (2/1/0) rather than centering on the heterozygote.

The output, and the cm column

The result has five metadata columns — snp, allele, chr, pos, cm — followed by one column per individual (with counted_column = TRUE, a sixth metadata column counted follows cm; see "Which allele is counted"). The types are the same for every input format: snp, allele and chr are character, pos is integer (double only if a position is not a whole number below 2^31, and integer NA when the input carries no positions), cm is double, and the genotype columns are integer. Numeric-format input (a data frame or a file written by as_numeric()) is brought to this schema too (logical and whole-number double genotype columns become integer, values unchanged), so converting a written file again returns the same table.

cm is the genetic map in centiMorgans, and it is NA unless the input format carries one. HapMap and VCF have no genetic-distance field, so cm is NA for them; PLINK .bim / .map files do (column 3), and it is kept when populated. This matters because meiosis needs recombination distances: as_population() and the crossing functions refuse to run on an all-NA map. Fill one in from the physical positions when you have no real map:

num <- as_numeric("my_genotypes.hmp.txt", to_r = TRUE)
num$cm <- synthetic_map(num$chr, num$pos)

Note also that aggregate(cm ~ chr, ...) returns "no rows to aggregate" on an all-NA map, since the formula interface drops missing values.

See Also

synthetic_map() to supply a genetic map, as_population() and cross() for what that map is needed for.

Examples

# A small HapMap table: 11 metadata columns, then one column per individual.
hmp <- data.frame(
  `rs#`   = c("snp1", "snp2", "snp3"),
  alleles = c("C/T", "A/G", "G/T"),
  chrom   = c(1, 1, 2),
  pos     = c(100, 500, 200),
  strand  = "+", `assembly#` = NA, center = NA,
  protLSID = NA, assayLSID = NA, panelLSID = NA, QCcode = NA,
  L1 = c("C", "A", "G"),
  L2 = c("T", "G", "T"),
  L3 = c("C", "A", "T"),
  L4 = c("T", "G", "G"),
  check.names = FALSE, stringsAsFactors = FALSE
)

# to_r = TRUE returns the result instead of writing it to disk.
num <- as_numeric(hmp, to_r = TRUE, verbose = FALSE)
num

# 0 / 1 / 2 coding instead of -1 / 0 / 1:
as_numeric(hmp, to_r = TRUE, code_as = "012", verbose = FALSE)

# Orient the coding against a reference allele rather than by frequency:
as_numeric(hmp, to_r = TRUE, verbose = FALSE,
           method = "reference", ref_allele = c("C", "A", "G"))

# From a file, the format is detected from its contents:
## Not run: 
num <- as_numeric("my_genotypes.hmp.txt", to_r = TRUE)
num <- as_numeric("my_genotypes.vcf", to_r = TRUE)

## End(Not run)

Description

Meiosis needs to know which allele sits on which of an individual's two homologous chromosomes, information that a -1/0/1 dosage matrix does not carry. as_population() builds that phased representation so the result can be passed to cross(), selfcross() and double_haploid().

Usage

as_population(geno, individuals = NULL, pool = NA_character_)

Arguments

Details

A Population is also accepted anywhere simulate_phenotype() takes geno, so a simulated pedigree can be phenotyped directly.

Value

A Population. Its individuals are recorded as pedigree founders.

Phasing of heterozygotes

Homozygous genotypes (-1 and 1) determine both strands exactly. Heterozygotes do not, and are assigned arbitrarily as allele-1 on the first strand and allele-2 on the second. For inbred panels such as SNP55K_maize282_maf04 this is nearly lossless because heterozygotes are rare, but for an outbred sample the phase of each heterozygote is a guess, and linkage between heterozygous sites in the first generation will not be realistic: every heterozygous site of such a founder sits in perfect coupling on the first strand, so its first-generation gametes carry the maximum coupling-phase linkage disequilibrium. The current public API does not import external haplotype phase, so as_population() should not be used for multi-generation recombination studies of substantially heterozygous, unphased founders.

Genetic map, units and chromosome order

The cm column must be a genetic map in centiMorgans (it is divided by 100 to give Morgans for meiosis). A map whose largest position is at most 5 across 20 or more markers looks like Morgans (or a proportion) and draws a warning: crossovers would be 100 times too rare. cm may start anywhere on a chromosome: the expected number of crossovers on a chromosome (Poisson with interference = "poisson", gamma-spaced by default; see cross()) equals its last map position in Morgans (the isqg convention), which is not its span max(cm) - min(cm) when the first marker is not at 0. Chromosomes are processed, and their random draws are consumed, in a fixed canonical order that does not depend on the storage type of chr or on the locale: labels that are numbers first, in numeric order (⁠1, 2, 10⁠), then other labels by their non-numeric prefix in byte order (uppercase before lowercase) and then by their trailing number (⁠chr1, chr2, chr10, chrX⁠). So the same map with chr stored as integers or as text gives the same seeded progeny. (Text labels that used to sort as ⁠"1", "10", "2"⁠ are now ordered ⁠1, 2, 10⁠; a seeded run on such a map draws its chromosomes in the new order.)

Allele orientation (crossing populations built from separate files)

The -1/0/1 coding is relative: +1 is the allele that as_numeric() considered the reference (by default the most frequent allele of that data set). Two populations made from separately converted panels can therefore code opposite alleles as +1 at a marker, and crossing them (or pooling them with c.Population()) would then silently mix up the alleles. Convert the panels jointly, or give each the same reference alleles with as_numeric(method = "reference", ref_allele = ). as_numeric() records the allele it counted as +1 at every marker (the "counted_allele" attribute of its result, see as_numeric()), and as_population() keeps it; the optional counted column of as_numeric(counted_column = TRUE) (placed after cm) carries the same record through text files and row subsetting and, when present, is the record used. When both populations carry that record, cross() and c.Population() compare it per marker and stop if the two panels count different alleles as +1. Where a record is missing on either side (numeric data read back from a text file or subsetted without the counted column, other software), the allele label is compared instead: the check warns when the two panels list a marker's alleles in opposite order (a sign the orientation may differ) and stops when they share no allele. The check only sees what these records show: it cannot detect a difference that neither records.

See Also

cross(), selfcross(), double_haploid(), synthetic_map()

Examples

data("SNP55K_maize282_maf04")
pop <- as_population(SNP55K_maize282_maf04, individuals = c("33-16", "38-11"))
pop

Description

Each individual's expected fraction of genes from each founder pool (breed, heterotic group), traced through its recorded pedigree: a founder is 1 for its own pool (the pool given to as_population()), and a progeny is the mean of its parents (a self or doubled haploid equals its parent). These are expected fractions – the realized fraction of an individual's genome from each breed varies around them and is not tracked (a Population carries no per-locus founder origin).

Usage

breed_composition(pop)

Arguments

Value

A matrix, individuals by pools, rows summing to 1. Founders imported without a pool label count toward a column "<unassigned>".

See Also

crossbreed(), heterosis(), parentage()

Examples

g <- data.frame(snp = paste0("m", 1:10), allele = "A/G", chr = 1,
                pos = 1:10, cm = seq(0, 90, by = 10), P1 = 1L, P2 = -1L)
A <- as_population(g[, 1:6], pool = "A"); B <- as_population(g[, -6], pool = "B")
f1 <- cross(A, B, n = 2, seed = 1)
bc <- cross(f1[1], A, n = 2, seed = 2)
breed_composition(bc)                               # 3/4 A, 1/4 B

Description

Advances a population as an undivided bulk: each generation every plant is selfed and the progeny are pooled, then n seeds are drawn at random from the pooled progeny to form the next generation (Bernardo 2020). Every plant is assumed to set ample seed, so the draw is uniform over the current plants: the number of advanced seeds each plant contributes is multinomial(n, 1/N) (mean n / N, variance about n / N), not fixed – including when n is a multiple of the current size N. Line identity is not tracked.

Usage

bulk(x, generations = 5L, n = NULL, seed = NULL, interference = NULL)

Arguments

Value

a Population (the final bulk).

References

Bernardo R (2020) Breeding for Quantitative Traits in Plants, 3rd ed. Stemma Press, Woodbury, Minnesota.

See Also

single_seed_descent(), pedigree(), recurrent_selection().

Examples

data("SNP55K_maize282_maf04")
pop <- as_population(SNP55K_maize282_maf04, individuals = 1:10)
f1  <- cross(pop[1], pop[2], n = 30, seed = 1)
bk  <- bulk(f1, generations = 4, n = 30, seed = 2)
n_individuals(bk)

Description

Column-binds several Populations (e.g. the progeny of many crosses) into one. All inputs must share the same marker map: identical snp names and chromosomes, and pos and cm equal to a numerical tolerance (absolute 1e-8 for values below 1, relative 1e-8 above – so round-off noise is accepted, any real map difference is not). The map of the first input is kept. Individual ids are made unique across the pooled set.

Usage

## S3 method for class 'Population'
c(...)

Arguments

Value

a single Population.

Examples

data("SNP55K_maize282_maf04")
pop <- as_population(SNP55K_maize282_maf04, individuals = 1:6)
a <- selfcross(pop[1], n = 3, seed = 1)
b <- selfcross(pop[2], n = 3, seed = 2)
c(a, b)

Description

The realized co-heritability matrix of a multi-trait simulation: $\mathrm{Cov}(G_i, G_j) / \sqrt{V_{P_i} V_{P_j}} = r_{G_{ij}} h_i h_j$, the genetic covariance standardized by the phenotypic standard deviations. G are the total genetic values, so this is the broad-sense co-heritability and the diagonal the realized broad-sense heritability Var(G_i) / Var(P_i). Only for a purely additive model does it equal the narrow-sense co-heritability $r_{A_{ij}} h_i h_j$ that predicts the correlated response to phenotypic selection, $CR_j = i \, r_{A_{ij}} h_i h_j \sigma_{P_j}$ (Falconer & Mackay 1996); with dominance or epistasis the non-additive covariance is not transmitted (e.g. pure dominance at p = 0.5 gives a positive value here but no response). G are genetic_values(); for a liability_threshold() trait the liability is used, the scale on which its genetic model is defined.

Usage

coheritability(sim, rep = 1L)

Arguments

Value

a symmetric ⁠n_traits x n_traits⁠ matrix with Trait_k dimnames.

References

Falconer DS, Mackay TFC (1996) Introduction to Quantitative Genetics, 4th ed. Longman, Harlow (correlated response to selection).

See Also

genetic_values(), liability_threshold()

Examples

data("SNP55K_maize282_maf04")
sim <- simulate_phenotype(SNP55K_maize282_maf04, n_traits = 2,
                          architecture = "pleiotropy", cor = 0.6,
                          h2 = c(0.5, 0.3), n_qtn = 20, seed = 1)
coheritability(sim)

Description

Scores candidates by their merit in crosses – general combining ability (GCA, the average merit of a parent in its crosses) and specific combining ability (SCA, the deviation of a particular cross from what its parents' GCAs predict), following the two-way model of Sprague & Tatum (1942) and, for the diallel, Griffing (1956). This is the criterion of half-sib / reciprocal recurrent selection and hybrid development. It is not the per se value (genotypic_value()) and, in general, not the transmissible breeding value.

Usage

combining_ability(
  candidates,
  testers = NULL,
  qtn,
  a,
  d = 0,
  design = c("topcross", "factorial", "diallel"),
  method = c("expected", "simulated"),
  n_progeny = NULL,
  h2 = NULL,
  var_e = NULL,
  ref = NULL,
  seed = NULL,
  interference = NULL
)

Arguments

Details

The architecture is frozen, as for genotypic_value(): loci qtn with additive effects a and dominance deviations d on the realized scale, so a locus has genotypic value -a / d / +a at gene content 0 / 1 / 2 (use template_effects() to take them from a simulation). Epistasis is outside this model.

Value

A combining_ability object: a list with gca (named, over candidates), gca_testers (factorial only), sca (matrix: candidates x testers, or the symmetric candidate x candidate matrix for a diallel, NA on its diagonal), cross_means (the Y matrix), grand_mean (mu), design, method, and for "simulated" also n_progeny, var_e and progeny (the progeny Population, pedigree recorded).

Expected (no simulation)

method = "expected" gives the conditional expectation of each cross's progeny value given the parents' genotypes – the idealized truth a simulation knows. It uses no random numbers. At one locus a parent with gene content x transmits the counted allele with probability x / 2, so for parents with g_i = x_i / 2 and g_k = x_k / 2 the progeny is AA, Aa or aa with probabilities $g_i g_k$, $g_i (1 - g_k) + (1 - g_i) g_k$ and $(1 - g_i)(1 - g_k)$, and its expected value is (this package's derivation, summed over loci)

$E[G_{ik}] = g_i g_k a + [g_i (1 - g_k) + (1 - g_i) g_k] d - (1 - g_i)(1 - g_k) a.$

Averaged over a set of testers this is linear in the candidate's gene content, with slope half the tester-referenced average effect $\alpha_T = a + d (1 - 2 p_T)$ (p_T the testers' mean gamete frequency): Falconer & Mackay's (1996) average effect $\alpha = a + d (q - p)$ at the testers' allele frequencies. So a candidate's GCA depends on the tester, and with testers drawn from the candidates' own population (p_T = p) it equals half its breeding value; with d = 0 every expected SCA is exactly zero. Expectations depend only on the parents' gamete frequencies, so linkage does not enter them (it enters the variance of a realized family).

Simulated (an estimate)

method = "simulated" realizes n_progeny progeny of every cross with the crossing core (mate()), scores them with genotypic_value(), optionally adds an independent residual per progeny (exactly one of h2 / var_e, as in phenotype_value(), h2 being the heritability of the total genotypic value in ref, default the progeny), and averages per cross. It is a finite-sample estimate of the expected values, with Mendelian (and, with a residual, environmental) sampling error. Consequently a simulated SCA is not zero for d = 0 when the parents are heterozygous (segregation among the n_progeny progeny of a cross is a sampling deviation from the expectation, shrinking with n_progeny); it is zero for d = 0 only in expectation, and exactly in a simulation from fully inbred parents (identical gametes, so no Mendelian sampling) that has no residual: an independent residual per progeny is environmental sampling error, which inbred parents do not remove, so with h2 / var_e the simulated SCA is not zero even then (it shrinks with n_progeny). ref, like h2, only scales the residual and is an error without h2. With seed the caller's random-number stream is left as it was found.

Designs and centering

• "topcross": every candidate crossed to every tester; GCA over candidates.

• "factorial": the same crosses (North Carolina Design II), with GCA reported for candidates and for testers.

• "diallel": every unordered pair of distinct candidates (testers unused; Griffing's method 4 layout, at least three candidates). With Y the cross means and mu their mean, the factorial / topcross GCA is the row (column) mean minus mu and SCA = Y - mu - g_i - g_k. In the diallel each cross involves two candidates, so this package's least-squares solution is g_i = (m_i - mu)(p - 1)/(p - 2), with m_i the mean of candidate i's p - 1 crosses; again SCA = Y - mu - g_i - g_k. In every design GCAs sum to zero and each candidate's SCAs sum to zero. Only the diallel without reciprocals is modelled (Griffing's method 4, one cross per unordered pair): reciprocal or maternal effects are outside the model and the cross matrix is symmetric. Individuals are identified by pedigree key, which includes the founders' pool label: two populations built from the same genotypes under different as_population(pool = ) labels are different individuals (a cross between them is recorded as a cross, not a self, although its expected value is the same).

References

Falconer DS, Mackay TFC (1996) Introduction to Quantitative Genetics, 4th ed. Longman, Harlow – the genotypic values $-a, d, a$ and the average effect $\alpha = a + d(q - p)$.

Sprague GF, Tatum LA (1942) General vs. specific combining ability in single crosses of corn. Journal of the American Society of Agronomy (now Agronomy Journal) 34:923–932. doi:10.2134/agronj1942.00021962003400100008x

Griffing B (1956) Concept of general and specific combining ability in relation to diallel crossing systems. Australian Journal of Biological Sciences 9:463–493. doi:10.1071/BI9560463

See Also

template_effects(), genotypic_value(), mate(), select_ind()

Examples

data("SNP55K_maize282_maf04")
pop <- as_population(SNP55K_maize282_maf04, individuals = 1:8)
q <- c("ss196442916", "ss196439337", "ss196480535")
ca <- combining_ability(pop[1:5], pop[6:8], qtn = q, a = c(1, 0.5, 0.25),
                        d = c(0.5, 0.5, 0))
ca$gca
ca$sca

Description

complex_phenotypes() merges two or more realized phenotype_sim objects built on the same genotypes with the same number of traits. The inputs' genetic values are summed (so each contributes in proportion to its genetic variance), the sum is scaled to the requested genetic variance share h2, and a common residual is added; the inputs' individual residuals are discarded. Realized H2 is reported from the resulting sample and can differ slightly because the generated genetic and residual vectors need not be exactly orthogonal. This recreates partial pleiotropy by combining, for example, a "pleiotropy" model with an "independent" model.

Usage

complex_phenotypes(
  ...,
  h2,
  reps = 1,
  resid_cor = NULL,
  residual_mode = c("fixed", "random")
)

Arguments

Details

Every input must be a complete model, exactly as for phenotypes_long(): if an input set an h2, its layer proportions must fill it (otherwise complex_phenotypes() errors with the same "Incomplete h2 allocation" message). The result is terminal: it has no layers, so additive(), dominance(), epistasis(), vqtl() and transcriptome() refuse to extend it (build the full model first, then combine). It carries no per-input state: mediation_split() is NULL for a combined model (the environmental expression-mediated part of an input is discarded with the input's residual, so no split of the combined phenotype is defined), and breeding-value selection methods that need per-locus effects refuse it.

Value

a combined phenotype_sim (architecture "complex").

Examples

data("SNP55K_maize282_maf04")
pleio <- simulate_phenotype(SNP55K_maize282_maf04, architecture = "pleiotropy",
                            n_traits = 2, seed = 10, cor = 0.5)
pleio <- additive(pleio, prop = 0.4, n_qtn = 3)
indep <- simulate_phenotype(SNP55K_maize282_maf04, n_traits = 2, seed = 11)
indep <- additive(indep, prop = 0.3, n_qtn = 3)
both <- complex_phenotypes(pleio, indep, h2 = 0.5)

Description

The first-order autoregressive correlation matrix $R_{ij} = \rho^{|i - j|}$: adjacent traits correlate rho, and the correlation decays geometrically with the distance between them. Use it for repeated measurements of one trait over time points or environments (high-throughput phenotyping), as the genetic correlation target cor of architecture = "pleiotropy" and/or the residual correlation resid_cor of simulate_phenotype().

Usage

cor_ar1(n_traits, rho)

Arguments

Value

an ⁠n_traits x n_traits⁠ positive-definite correlation matrix.

See Also

simulate_phenotype()

Examples

cor_ar1(4, 0.8)
data("SNP55K_maize282_maf04")
# five time points: genetic and residual correlations decay with time lag
sim <- simulate_phenotype(SNP55K_maize282_maf04, n_traits = 5,
                          architecture = "pleiotropy", cor = cor_ar1(5, 0.9),
                          resid_cor = cor_ar1(5, 0.5), h2 = 0.4, n_qtn = 30,
                          seed = 1)

Description

Superseded. create_phenotypes() is the frozen original engine, retained unchanged for backward compatibility (full v1 signature, legacy seed arithmetic, and RNGversion("3.5.1")). It is bugfix-only and is not a wrapper over the v2 grammar (the two implementations coexist). For new work prefer the composable grammar: simulate_phenotype() with additive(), dominance(), epistasis(), vqtl() and complex_phenotypes().

Usage

create_phenotypes(
  geno_obj = NULL,
  geno_file = NULL,
  geno_path = NULL,
  QTN_list = list(add = list(NULL), dom = list(NULL), epi = list(NULL), var = list(NULL)),
  prefix = NULL,
  rep = NULL,
  ntraits = 1,
  h2 = NULL,
  mean = NULL,
  model = NULL,
  architecture = "pleiotropic",
  add_QTN_num = NULL,
  dom_QTN_num = NULL,
  epi_QTN_num = NULL,
  var_QTN_num = NULL,
  epi_type = NULL,
  epi_interaction = 2,
  pleio_a = NULL,
  pleio_d = NULL,
  pleio_e = NULL,
  trait_spec_a_QTN_num = NULL,
  trait_spec_d_QTN_num = NULL,
  trait_spec_e_QTN_num = NULL,
  add_effect = NULL,
  dom_effect = NULL,
  epi_effect = NULL,
  var_effect = NULL,
  remove_add_effect = FALSE,
  same_add_dom_QTN = FALSE,
  same_mv_QTN = FALSE,
  big_add_QTN_effect = NULL,
  degree_of_dom = 1,
  type_of_ld = "indirect",
  ld_min = 0.2,
  ld_max = 0.8,
  ld_method = "composite",
  sim_method = "geometric",
  vary_QTN = FALSE,
  cor = NULL,
  cor_res = NULL,
  QTN_variance = FALSE,
  seed = NULL,
  home_dir = NULL,
  output_dir = NULL,
  export_gt = FALSE,
  output_format = "long",
  to_r = FALSE,
  out_geno = NULL,
  chr_prefix = "chr",
  remove_QTN = FALSE,
  warning_file_saver = TRUE,
  constraints = list(maf_above = NULL, maf_below = NULL, hets = NULL),
  maf_cutoff = NULL,
  nrows = Inf,
  na_string = "NA",
  SNP_effect = "Add",
  SNP_impute = "Middle",
  quiet = FALSE,
  verbose = TRUE,
  RNGversion = "3.5.1"
)

Arguments

Value

Single or multi-trait phenotypes in one of many formats. Numericalized marker data set with or without the selected QTNs. Diagnostic files (log, QTN information, summary of LD between QTNs, proportion of phenotypic variance explained by each QTN).

Rejected inputs and unsupported combinations

The frozen engine validates its arguments up front and stops with a message that names the argument and the remedy, instead of failing cryptically or returning wrong numbers. Valid inputs produce exactly the numbers of earlier releases. Inner errors are re-signalled (earlier releases printed the message and returned NULL), and on error only the run folder created by the call is removed. The following are rejected:

• QTN_list with ntraits = 1 (a single trait cannot be simulated from a user-specified marker list; use random QTNs or simulate_phenotype() with ⁠qtn =⁠), with same_add_dom_QTN = TRUE, with architecture = "LD", or with a number of trait-specific marker vectors different from ntraits.

• Single-trait models combining dominance and epistasis ("DE", "ADE") and single-trait same_add_dom_QTN = TRUE.

• architecture = "LD" with ntraits > 2, with ld_max >= 1 (an absolute LD of 1 is met by a marker paired with itself) or ld_min > ld_max, and model = "D" with type_of_ld = "indirect". When the LD search cannot deliver distinct, linked marker pairs inside ⁠[ld_min, ld_max]⁠ for the chosen seed (indirect: a marker selected for both traits or a duplicated QTN; direct: a pair outside the window, on two chromosomes or paired with itself), the call stops with an "LD contract" error: change seed, the LD window or type_of_ld. For type_of_ld = "direct" the search is first run exactly as in earlier releases and, only when that attempt fails, repeated with derived seeds (section "Direct-LD search retries"), so the error is now rare.

• h2 outside ⁠[0, 1]⁠; h2 positive but not above 0.05 (that is, ⁠0 < h2 <= 0.05⁠, including exactly 0.05) with rep > 1 (see seed); to_r = TRUE with several rows of h2 unless a single trait is simulated with vary_QTN = FALSE.

• cor that is not a symmetric positive-definite ⁠ntraits x ntraits⁠ matrix (it is no longer repaired silently) and cor_res that is not a symmetric, positive semi-definite correlation matrix.

• Effect sizes given as a mixture of one value per QTN and a single geometric-series base, or with a length that is neither 1 nor the number of QTNs.

• Variance QTL ("V") with ntraits > 1, with h2 = 0, or with var_effect values that make the standard-deviation multiplier negative, or with a constant additive baseline (for example model = "AV" with add_effect = 0: the variance multiplier is undefined, so the call stops instead of returning NaN phenotypes).

• Unknown output_format, output_format = "wide" with ntraits > 1 and rep = 1, rep/model/seed that are missing or malformed, and more QTNs than markers.

Duplicated marker positions

Duplicated marker positions (chr_pos, the same chromosome and position on several rows) are accepted by the fully pleiotropic and the partially pleiotropic architectures: QTN columns are addressed by row position, never by the ⁠Chr_<chr>_<pos>⁠ name, so duplicated names cannot change which markers enter a genetic value (an epistatic interaction, for example, multiplies the intended columns). The architecture = "LD" searches are stricter: a duplicated QTN or a QTN pair sharing a chromosome position is reported as an "LD contract" error (see above).

Direct-LD search retries

With architecture = "LD" and type_of_ld = "direct" the marker search of every replicate first runs exactly as in earlier releases (same seeds, same random draws); whenever that attempt meets the LD contract (distinct markers on one chromosome with an absolute LD inside the inclusive window ⁠[ld_min, ld_max]⁠, and the reported LD equal to the recomputed one) its output is unchanged, bit for bit. Only when it does not (the contract check fails, the search stops with "None of the selected SNPs met the minimum/maximum LD threshold", or it runs past the first/last marker) the replicate is searched again, up to 50 attempts in total, from the seed seed - (a - 1) * 1000003 for attempt a (seed + (a - 1) * 1000003 when seed <= 0); the retry seeds are then combined with the replicate and QTN offsets exactly like the original one. Retry attempts also reset the neighbour pointers after a marker is re-drawn, which the earlier dominance walks (model containing "D") did not: they kept pairing the new marker with the neighbours of the old one, so most pairs failed the contract. The window itself is never relaxed: every accepted pair is verified against ⁠[ld_min, ld_max]⁠ and the chromosome constraint. After 50 failed attempts the "LD contract" error is raised as before; a retry attempt that uses up every candidate marker without finding a pair inside the window (a window too narrow for the data, which used to end in the cryptic "invalid first argument") stops the search at once with an informative error. The seed of the attempt that succeeded is the one written to the ⁠Seed_num_for_*⁠ file and, when verbose = TRUE, a message names the attempt. Indirect LD is not retried. A dominance-only model ("D") can still stop with the "All individuals are homozygote for the selected dominance QTNs" error, which is a separate check on the selected markers.

Author(s)

Samuel B Fernandes and Alexander E Lipka. Last update: Jan 19, 2021

References

Fernandes, S.B., and Lipka, A.E., 2020 simplePHENOTYPES: SIMulation of pleiotropic, linked and epistatic SIMulation of Pleiotropic, Linked and Epistatic PHENOTYPES. BMC Bioinformatics 21(1):491, doi:10.1186/s12859-020-03804-y

Examples

# Simulate 50 replications of a single phenotype.
data("SNP55K_maize282_maf04")
pheno <- 
  create_phenotypes(
    geno_obj = SNP55K_maize282_maf04,
    add_QTN_num = 3,
    add_effect = 0.2,
    big_add_QTN_effect = 0.9,
    rep = 10,
    h2 = 0.7,
    model = "A",
    to_r = TRUE,
    home_dir = tempdir(),
    quiet = T
    )
# For more examples, please run the following:
# vignette("simplePHENOTYPES")

Description

Produces n progeny from a biparental cross. Each progeny receives one recombinant gamete from each parent, so the two parents contribute one homologue apiece.

Usage

cross(mother, father, n = 1, seed = NULL, interference = NULL)

Arguments

Details

Each gamete carries on average one crossover per Morgan of a chromosome's length. By default their spacing follows the gamma model of crossover interference (nu = 2.6, p = 0, AlphaSimR's default; see the section "Crossover interference"). With interference = "poisson" recombination follows the count-location model of isqg: the number of crossovers on a chromosome is Poisson with mean equal to its length in Morgans, and their positions are uniform along it. That length is the chromosome's last map position (cm / 100), not its span max(cm) - min(cm): positions are used as given and are not rebased to a zero origin (the isqg convention). The two agree when a chromosome's first marker is at 0; when it is not, the extra crossovers fall upstream of the first marker and only swap the whole chromosome, so the recombination fraction between markers is still Haldane's, $(1 - e^{-2d})/2$. Chromosomes assort independently. The genetic map is taken from the cm column of the population's marker map, so it must not be missing - see synthetic_map() if you only have physical positions.

Both parents must carry the same marker map (marker names, chromosomes and positions; the allele column is not part of it). The dosages are relative to the allele coded +1, which as_numeric() chooses per data set; crossing populations built from separately converted panels can therefore mix up the alleles. Convert the panels together, or with as_numeric(method = "reference", ref_allele = ). Where both populations record the allele column, a marker whose alleles are listed in opposite order draws a warning and markers with no allele in common an error (see as_population()).

Chromosomes are processed, and their random draws consumed, in a canonical order that does not depend on the storage type of chr or on the locale (see as_population()).

Value

A Population of n progeny.

Crossover interference

By default (interference = NULL, with the option below unset) crossovers follow a two-pathway gamma model of interference with nu = 2.6 and p = 0: AlphaSimR's default (its SimParam fields v and p), which approximates Kosambi's map function. Give interference = list(nu = , p = ) for other values, or interference = "poisson" for Poisson crossovers (no interference), whose random draws are exactly isqg's. Before version 2.0.0.9003 Poisson was the default; seeded results of every function that runs meiosis changed with the new default (see NEWS.md). The gamma model:

• The four-strand bivalent carries chiasmata at an average of 2 per Morgan. A gamete takes part in each chiasma independently with probability 1/2 (no chromatid interference), so it carries one crossover per Morgan on average: the genetic map (cm) keeps its meaning for every nu and p, as under Haldane.

• A share p of the chiasmata (intensity $2p$ per Morgan) is non-interfering: a Poisson process.

• The rest (intensity $2(1 - p)$) is interfering: a stationary renewal process whose gaps are gamma distributed with shape nu and rate $2 \nu (1 - p)$ per Morgan (mean gap $1 / (2 (1 - p))$ Morgans). ⁠1 <= nu <= 1e6⁠ is the interference strength (larger values are rejected: the gamma gaps are then numerically a fixed spacing); nu = 1 is exponential gaps, i.e. no interference. The first chiasma is drawn from the stationary (equilibrium) distribution, so the process does not "start" at the chromosome end.

• The two pathways are independent; the gamete's crossovers are the union.

The recombination fraction between two loci $d$ Morgans apart is then $r(d) = [1 - P_0(d)] / 2$, with $P_0(d)$ the probability of no chiasma in the interval: $P_0(d) = e^{-2pd} \{1 - F_{\nu+1}(y) - (y / \nu) [1 - F_\nu(y)]\}$, $y = 2 \nu (1 - p) d$, $F_a$ the Gamma(a, 1) distribution function. With nu = 1, or p = 1, the model is Poisson in distribution (Haldane), although it then uses its own random stream, not isqg's. p is in ⁠[0, 1]⁠; nu must be in ⁠[1, 1e6]⁠ (negative interference is not modelled). The gamma model of interference is that of McPeek and Speed (1995) and the two-pathway extension that of Housworth and Stahl (2003) (see References; the formulas above were derived and are checked numerically in this package's tests, not copied from either paper or from AlphaSimR). The draws are made in R; the Rust core only applies the drawn crossovers and is unchanged.

To use one interference model for a whole session or scheme without passing ⁠interference =⁠ to every call, set the option simplePHENOTYPES.interference to "poisson" (the old default, e.g. to reproduce results of earlier versions) or to a list(nu = , p = ) (NULL, the default, is the gamma model above). An interference argument that is NULL is then replaced by the option, in cross(), selfcross(), double_haploid(), mate(), crossbreed() and in every function that forwards it (single_seed_descent(), bulk(), pedigree(), recurrent_selection(), cross_usefulness(), combining_ability(), progeny_test()), so the precedence is: an explicit non-NULL argument, then the option, then the default gamma model. The option is validated by the same rules as the argument, and an invalid value is an error that names the option. A call overrides the option with an explicit value (interference = "poisson" or a list). With the option set, the draws are exactly those made when the same value is passed to every call.

References

Toledo, F.H., Perez-Rodriguez, P., Crossa, J. and Burgueno, J. (2019). isqg: A Binary Framework for in Silico Quantitative Genetics. G3 Genes|Genomes|Genetics 9(8), 2425–2428. doi:10.1534/g3.119.400373

McPeek, M.S. and Speed, T.P. (1995). Modeling interference in genetic recombination. Genetics 139(2), 1031–1044. doi:10.1093/genetics/139.2.1031

Housworth, E.A. and Stahl, F.W. (2003). Crossover interference in humans. American Journal of Human Genetics 73(1), 188–197. doi:10.1086/376610

See Also

selfcross(), double_haploid(), as_population()

Examples

data("SNP55K_maize282_maf04")
pop <- as_population(SNP55K_maize282_maf04, individuals = c("33-16", "38-11"))
f1 <- cross(pop[1], pop[2], n = 5, seed = 1)
f1

Description

For each candidate biparental cross, simulates a progeny family and reports its additive genetic mean, genetic standard deviation, and usefulness $U = \mu + i\,\sigma$ (Zhong & Jannink 2007; Lehermeier et al. 2017). A cross with a lower mean but more genetic variance can outrank a safer one, which is the information parent-mean ranking misses.

Usage

cross_usefulness(
  sim,
  pairs = NULL,
  scheme = c("dh", "selfcross", "cross"),
  n_progeny = 100L,
  generations = 5L,
  select_top = 0.1,
  trait = 1L,
  direction = c("high", "low"),
  seed = NULL,
  interference = NULL
)

Arguments

Details

$U = \mu + i\,\sigma$ is a normal-theory criterion: it equals the expected genetic value of the best select_top fraction of progeny only when that family's genetic values are approximately normally distributed, with i the standardized selection intensity for the fraction. For a large polygenic family this holds well; for a small family or one driven by a few large-effect QTNs the distribution is discrete/skewed and U is an approximation that can exceed the realized top-fraction mean (and even the family maximum). Use it as a ranking criterion under that assumption; read the realized family from the returned mean/sd (or simulate the family directly with cross()/selfcross()) if you need the exact tail mean.

The intensity is the infinite-population normal value $i = \varphi(\Phi^{-1}(1-p))/p$ for select_top = p (1.755 at p = 0.1); the mean of the top p of a finite family is slightly smaller (about 1.73 for the top 10 of 100), and a family with n_progeny * select_top < 1 selects less than one progeny, so U is then only a normal-theory ranking score. Prefer n_progeny * select_top >= 1, in practice a few tens of progeny at least.

Heterozygous parents and "dh" / "selfcross". These schemes derive each family from a single F1 individual of the cross (one draw of the two parental gametes), then double or self it. For inbred (homozygous) parents that F1 is unique, so the family-generating distribution is that of the cross itself (no conditioning on a random F1). The returned mean and sd are nevertheless Monte Carlo estimates from n_progeny simulated progeny (meiosis and, for "selfcross", the selfing draws remain random), so they still vary with seed and converge to the cross's moments only as n_progeny grows. For a heterozygous parent the F1 is also one random draw, so the reported mean and sd are additionally conditional on which F1 was drawn: they change with seed for that reason as well, and the standard deviation sd (its square, the family variance, likewise) is on average lower than the F1-averaged (many-F1) value. The function warns once per call when a parent in pairs is heterozygous at a scoring locus under these schemes. Use inbred parents, or scheme = "cross" (which draws every progeny from an independent meiosis of the parents), when the parents are not inbred; the sampling scheme itself is unchanged.

Progeny are scored on the additive QTN effects already realized in sim (the template), applied on a common scale, so the means and standard deviations are directly comparable between crosses. Only the additive component enters the score (the breeding-value basis of the usefulness criterion); dominance and epistasis layers are ignored, and doubled-haploid / inbred families carry no dominance anyway.

Value

a data frame with parent1, parent2, mean, sd, usefulness, sorted best first, carrying attribute intensity (the i used).

References

Zhong S, Jannink J-L (2007) Using quantitative trait loci results to discriminate among crosses on the basis of their progeny mean and variance. Genetics 177(1):567–576. doi:10.1534/genetics.107.075358; Lehermeier C, Teyssedre S, Schon C-C (2017) Genetic gain increases by applying the usefulness criterion with improved variance prediction in selection of crosses. Genetics 207(4):1651–1661. doi:10.1534/genetics.117.300403

See Also

cross(), single_seed_descent(), select_ind().

Examples

data("SNP55K_maize282_maf04")
pop <- as_population(SNP55K_maize282_maf04, individuals = 1:6)
sim <- simulate_phenotype(pop, h2 = 0.5, seed = 1) |> additive(n_qtn = 30)
u <- cross_usefulness(sim, scheme = "dh", n_progeny = 30, seed = 2)
head(u)   # crosses ranked by expected best-progeny value

Description

Runs a classical crossbreeding system on pure-breed Populations and returns the final crossbred generation, its pedigree traced to the breeds (see breed_composition()). Each generation mates n_progeny random dam-sire pairs (one progeny each) with mate(); sex is not modelled, so "dam" and "sire" are the roles of the two pools in the plan, not sexes.

Usage

crossbreed(
  breeds,
  system = c("two_way", "backcross", "three_way", "terminal", "rotational"),
  n_progeny,
  generations = 1L,
  sire_breed = NULL,
  seed = NULL,
  interference = NULL
)

Arguments

Details

• "two_way": breed 1 dams x breed 2 sires (an F1).

• "backcross": F1 dams x breed-1 sires.

• "three_way": F1 (breed 1 x breed 2) dams x breed-3 sires.

• "terminal": F1 (breed 1 x breed 2) dams x sires of sire_breed (a breed not in the F1); the products are not bred further.

• "rotational": F1 dams, then each generation's crossbred dams x sires of the next breed in rotation (breed 1, 2, ... then back to 1), for generations generations after the F1. The expected composition converges (for two breeds, toward alternating 2/3 : 1/3).

Value

The final generation as a Population, with attribute history: a data frame with one row per generation (generation, sire_breed and the mean expected breed fraction per breed).

See Also

breed_composition(), heterosis(), mate()

Examples

g <- data.frame(snp = paste0("m", 1:10), allele = "A/G", chr = 1,
                pos = 1:10, cm = seq(0, 90, by = 10))
mk <- function(pool, v) as_population(cbind(g, matrix(v, 10, 4,
  dimnames = list(NULL, paste0(pool, 1:4)))), pool = pool)
br <- list(A = mk("A", 1L), B = mk("B", -1L))
rot <- crossbreed(br, "rotational", n_progeny = 20, generations = 6, seed = 1)
attr(rot, "history")

Description

Add a dominance variance component

Usage

dominance(
  sim,
  prop = NULL,
  same_as_add = TRUE,
  n_qtn = NULL,
  qtn = NULL,
  dist = "geometric",
  effect = NULL
)

Arguments

Details

Dominance is modelled as a deviation applied to heterozygotes (the het indicator), with its share of phenotypic variance set by prop. In this variance-partition grammar there is no separate "degree of dominance" argument (there is no ⁠degree =⁠): the ratio of the scaled components' variances is prop_dominance / prop_additive, set through the layer proportions. (A single degree-of-dominance scalar would be washed out by the per-component variance scaling and is therefore not offered; for a per-locus degree of dominance d / abs(a) use additive(orthogonal = TRUE, a =, d =).)

Additive + dominance on the same loci is not orthogonal. The dosage and the heterozygote indicator covary by $-(2p - 1) 2pq$ at each locus, so the realized genetic variance is ⁠prop_A + prop_D + 2Cov(c_A, c_D)⁠ for one additive and one dominance layer (several layers of a type add their mutual covariance to Var(c_A) or Var(c_D)) – a structural, allele-coding-dependent bias, not a finite-sample effect; its per-locus sign is that of ⁠1 - 2p⁠ for the counted-allele frequency p, so it is one-signed across loci only when those frequencies are on one side of 0.5 – and the classical Va:Vd of the block differs from prop_A:prop_D. The object reports the realized Var(A), Var(D) and 2Cov(A,D), and Var(c_A), Var(c_D) and 2Cov(c_A,c_D) (⁠$ad_report⁠, and the note printed by print()); use additive(orthogonal = TRUE, a =, d =) for a Fisher-orthogonal partition.

Dominance needs heterozygotes to be identifiable: it is identically zero at a locus with no heterozygous individuals. If the selected loci carry none (common on a near-inbred panel such as the bundled maize lines, ~0.4% heterozygous), dominance() errors rather than silently substituting loci; if only some of them carry none, it warns that the rest absorb prop. Pre-filter to heterozygous markers with filter_geno(hets = "include"), fix het-bearing loci with ⁠qtn =⁠, or simulate dominance on an outbred or F2-type population instead.

Under architecture = "pleiotropy" the dominance effects are drawn with the same covariance as the additive layer (DECISION-023): shared loci carry correlated effects and trait-specific loci independent ones, split by pi, so the dominance component targets cor – its realized correlation converges to cor as the loci and individuals grow, for loci in approximate linkage equilibrium; see simulate_phenotype() cor – rather than every trait getting one identical effect series. With same_as_add = TRUE the additive layer's shared and trait-specific loci are reused. Each locus's effect is scaled by the realized standard deviation of its heterozygote indicator. Fixing loci with ⁠qtn =⁠ works under every architecture: under "pleiotropy" every locus affects every trait (the same loci for each trait), with an explicit effect setting the effects when no correlation is controlled; under "ld" dominance reuses the additive layer's linked loci (same_as_add = TRUE) or takes qtn = list(trait1_loci, trait2_loci) of disjoint linked pairs (no marker may already be causal for the other trait in an earlier layer), since a fresh random draw could make a marker causal for both traits across layers.

Value

the updated phenotype_sim.

See Also

additive(), epistasis(), vqtl(), filter_geno().

Examples

data("SNP55K_maize282_maf04")

# Dominance reuses the additive QTNs by default, so the trait has
# h2 = 0.4 + 0.1 = 0.5 spread over one set of loci.
ph <- simulate_phenotype(SNP55K_maize282_maf04, seed = 1) |>
  additive(prop = 0.4, n_qtn = 5) |>
  dominance(prop = 0.1)
ph

# same_as_add = FALSE draws separate dominance-only loci instead.
simulate_phenotype(SNP55K_maize282_maf04, seed = 2) |>
  additive(prop = 0.4, n_qtn = 5) |>
  dominance(prop = 0.1, same_as_add = FALSE, n_qtn = 3)

# Give the dominance effects yourself: a geometric base, or one value per QTN.
simulate_phenotype(SNP55K_maize282_maf04, seed = 2) |>
  additive(prop = 0.4, n_qtn = 5) |>
  dominance(prop = 0.1, same_as_add = FALSE, n_qtn = 3,
            effect = c(0.6, 0.3, 0.1))

Description

Dosage matrix of a Population

Usage

dosages(x)

Arguments

Value

An integer matrix of markers by individuals, coded -1/0/1.

Examples

data("SNP55K_maize282_maf04")
pop <- as_population(SNP55K_maize282_maf04, individuals = 1:3)
dim(dosages(pop))

Description

Produces n doubled-haploid progeny: a single recombinant gamete is drawn and then doubled, so every individual is completely homozygous at every marker. This is the one-generation route to a fully inbred line, and the result contains no heterozygotes at all.

Usage

double_haploid(parent, n = 1, seed = NULL, interference = NULL)

Arguments

Value

A Population of n fully homozygous progeny.

Many doubled-haploid families in one call

double_haploid() takes one parent. To make doubled haploids from many parents (or many families from one) call mate() with a plan whose rows have design = "dh", e.g. mate(data.frame(mother = ids, father = ids, n = 100, design = "dh"), pop): all rows are executed in one pass through the Rust core, with the same random draws, in the same order, and the same progeny and pedigree as calling double_haploid() on each parent in turn under one seed, at a fraction of the cost per call (see mate()).

References

Toledo, F.H., Perez-Rodriguez, P., Crossa, J. and Burgueno, J. (2019). isqg: A Binary Framework for in Silico Quantitative Genetics. G3 Genes|Genomes|Genetics 9(8), 2425–2428. doi:10.1534/g3.119.400373

See Also

cross(), selfcross()

Examples

data("SNP55K_maize282_maf04")
pop <- as_population(SNP55K_maize282_maf04, individuals = c("33-16", "38-11"))
f1 <- cross(pop[1], pop[2], n = 1, seed = 1)
dh <- double_haploid(f1, n = 10, seed = 3)
# No heterozygotes by construction:
any(dosages(dh) == 0)

Description

Add an epistatic variance component

Usage

epistasis(
  sim,
  prop = NULL,
  n_pairs = NULL,
  interaction = 2,
  interaction_type = "a",
  qtn = NULL,
  effect = NULL,
  dist = "geometric"
)

Arguments

Details

The epistatic value is the product of the interacting markers' design terms times a per-set effect, then centered and scaled to prop. Each design term is centered per locus before multiplying, which removes each locus's mean so the product is not dominated by a constant offset. This does not orthogonalize the interaction against the additive/dominance main effects: the centered product can still be correlated with the constituent dosages – strongly so when the interacting loci are in LD or away from allele frequency 0.5 (with two perfectly linked loci it can be collinear with each additive term). It is therefore not a Fisher/NOIA-orthogonal variance component (that is a planned genotypic-model mode), so the "epistatic proportion" is a simulation quantity, not an orthogonal share of variance.

interaction_type values of "d" depend on heterozygotes and are near-degenerate on a fully inbred panel (see dominance()); use "a" (the default) for inbred data, or an outbred/F2-type population for "d".

Under architecture = "pleiotropy" the interacting sets are split by pi into shared sets (common to every trait, with correlated effects) and trait-specific sets (independent effects), drawn with the same covariance as the additive layer, so the epistatic component targets cor (its realized correlation converges as the sets and individuals grow, for loci in approximate linkage equilibrium; DECISION-023). Each set's effect is scaled by the realized standard deviation of its centered-product design column; ⁠effect =⁠ is rejected there while a correlation is controlled (without cor / pi it sets the effects, and fixed ⁠qtn =⁠ sets are accepted). epistasis() is not supported under architecture = "ld": that architecture makes each trait's causal loci distinct markers in LD within the r2 window, and epistatic sets have no such linked construction – they would be unlinked and could reuse a marker another layer made causal for the other trait.

Value

the updated phenotype_sim.

See Also

additive(), dominance(), vqtl().

Examples

data("SNP55K_maize282_maf04")

# Two pairs of interacting loci contributing 10% of the variance,
# on top of an additive background.
ph <- simulate_phenotype(SNP55K_maize282_maf04, seed = 1) |>
  additive(prop = 0.4, n_qtn = 5) |>
  epistasis(prop = 0.1, n_pairs = 2)
ph

# interaction = 3 makes each epistatic term a three-way interaction
# rather than the default pairwise one.
simulate_phenotype(SNP55K_maize282_maf04, seed = 3) |>
  additive(prop = 0.3, n_qtn = 4) |>
  epistasis(prop = 0.2, n_pairs = 2, interaction = 3)

# Additive-by-dominance and dominance-by-dominance interactions depend on
# heterozygotes, so run them on a segregating F2 (the maize282 panel is
# inbred and has none).
pop <- as_population(SNP55K_maize282_maf04, individuals = 1:20)
f2 <- selfcross(cross(pop[1], pop[2], n = 1, seed = 1), n = 40, seed = 2)
simulate_phenotype(f2, seed = 4) |>
  additive(prop = 0.3, n_qtn = 4) |>
  epistasis(prop = 0.1, n_pairs = 2, interaction_type = c("a", "d")) |>
  epistasis(prop = 0.1, n_pairs = 2, interaction_type = c("d", "d"))

Description

Groups a Population's individuals by shared parents, for select_ind()'s ⁠family =⁠ argument (method = "within_family", "among_family" or "combined") or for progeny testing.

Usage

families(
  x,
  by = c("full_sib", "maternal_half_sib", "paternal_half_sib", "selfed")
)

Arguments

Value

A factor with one element per individual of x, named by id; levels are named by the parent id(s). Founders, and for "selfed" any individual that is not a self or doubled haploid, are NA.

See Also

parentage(), select_ind()

Examples

data("SNP55K_maize282_maf04")
pop <- as_population(SNP55K_maize282_maf04, individuals = 1:3)
fs <- c(cross(pop[1], pop[2], n = 4, seed = 1),
        cross(pop[1], pop[3], n = 4, seed = 2))
families(fs, "full_sib")
families(fs, "maternal_half_sib")

Description

One place to apply every marker filter simplePHENOTYPES implements, so a genotype can be trimmed once and then reused across simulations. It takes the simplePHENOTYPES numeric format (a data frame whose first five columns are snp, allele, chr, pos, cm (and optionally the counted column of as_numeric(counted_column = TRUE)), followed by one column per individual, e.g. SNP55K_maize282_maf04), an individuals-by-markers numeric matrix coded -1/0/1, or a Population (founders or the progeny of cross(), selfcross() and double_haploid()), and returns the same kind of object with the failing markers dropped. Every genotype input that simulate_phenotype() accepts can therefore be filtered, so the same MAF, heterozygosity and LD constraints apply whatever the architecture or layer (additive, dominance, epistasis, vQTL) drawn from it afterwards.

Usage

filter_geno(
  geno,
  maf_above = NULL,
  maf_below = NULL,
  hets = c("any", "include", "remove"),
  remove_monomorphic = TRUE,
  indep_pairwise = NULL,
  indep_pairphase = NULL,
  indep = NULL,
  blocks = FALSE,
  block_max_kb = 500,
  window_unit = c("variants", "kb"),
  code_as = c("-101", "012"),
  verbose = TRUE
)

Arguments

Details

Filters are applied in order: monomorphic removal, minor-allele frequency, heterozygosity, then LD pruning. A marker must pass every requested filter to be kept. There is no call-rate filter (PLINK's --geno / --mind) and no Hardy-Weinberg test: a marker with no called genotype is kept unless remove_monomorphic = TRUE (its minor-allele frequency is then 0), so drop or impute such markers beforehand if they matter.

Value

the filtered genotype object, in the same format as geno (a Population stays a Population).

LD pruning and haplotype blocks

The three sliding-window pruners take the same three-number specification used by PLINK (c(window, step, threshold)); window_unit sets whether window counts markers (default) or kilobases ("kb", using the pos column). indep_pairwise, indep and indep_pairphase reproduce PLINK 1.9's --indep-pairwise, --indep and --indep-pairphase marker-for-marker for ordinary use: on complete-call autosomal genotypes, with step <= window, the kept marker set (and hence the removed set) matches PLINK 1.9. The committed evidence is a set of golden PLINK v1.9.0-b.8 results on the bundled SNP55K_maize282_maf04 panel (nine prune configurations genome-wide, plus ⁠--blocks no-pheno-req⁠ at 200 and 500 kb), compared in the same marker order by tests/testthat/test-filter-geno-plink-parity.R, together with hard-coded threshold-epsilon and monomorphic corner cases; there is no PLINK run on other panels or on randomly generated genotypes. This holds for variant windows (all three) and kb windows (indep_pairwise); blocks = TRUE likewise reproduces the block definitions of PLINK's ⁠--blocks no-pheno-req⁠ at a matching --blocks-max-kb (simplePHENOTYPES imposes no phenotype requirement, and its block_max_kb defaults to 500, versus PLINK's 200-kb default – so a bare PLINK --blocks will differ unless run with no-pheno-req and the same span). All loci are treated as autosomal diploid: simplePHENOTYPES has no sex chromosome, so PLINK's haploid chromosome-X weighting is neither modelled nor matched.

Beyond that envelope the match is best-effort, not guaranteed, in a few pathological corners that do not arise in normal genotype panels: non-overlapping windows (step > window) combined with markers pruned mid-window; indep (VIF) when markers outnumber individuals, so the marker correlation matrix is rank-deficient and pruning turns on PLINK's exact matrix-inversion tie handling; indep when a pair is near-perfectly negatively correlated in PLINK's minor-allele coding (the VIF collinearity scan compares ⁠|r|⁠, which matches every real panel tested but not that hand-built sign case); indep_pairphase at a threshold within floating-point tolerance of a realised haplotypic r^2 (that r^2 comes from a transcendental cubic solver – acos/cos/log – so, unlike the integer-exact composite r^2, it can differ from PLINK's in the last bit); and kb + indep. PLINK computes LD from founders only; filter_geno() does not apply that rule: it computes LD from every individual it is given, including the non-founders of a crossing Population (its pedigree is not consulted; the parity is verified on founder-only data frames) – pass a founder-only genotype set, e.g. the founders before crossing, to reproduce a PLINK run that had non-founders. Missing calls (NA) are handled per pair on the complete observations (a pairwise.complete.obs correlation for indep_pairwise/indep, complete pairs only for indep_pairphase) and are not claimed identical to PLINK's missing-data path; drop or impute missing genotypes first for PLINK parity. Variant windows (window_unit = "variants") tolerate missing chr/pos (for example genotypes converted from a nucleotide table): markers with a missing chr form one group, ordered as supplied. window_unit = "kb" and blocks = TRUE need a physical map, so they stop with an error when chr or pos is missing for any marker.

indep_pairwise = c(window, step, r2)

PLINK --indep-pairwise: drop one marker of every pair whose composite genotype r^2 (squared correlation of the -1/0/1 dosage) exceeds r2, removing the lower-MAF marker of the pair (ties keep the earlier one), e.g. c(50, 5, 0.2).

indep_pairphase = c(window, step, r2)

PLINK --indep-pairphase: like indep_pairwise, but r^2 is the haplotypic (phased) r^2 of Hill & Robertson (1968), with the double-heterozygote phase resolved per pair by PLINK's exact two-locus maximum-likelihood solver, rather than the composite genotype r^2, e.g. c(50, 5, 0.2).

indep = c(window, step, vif)

PLINK --indep: multi-marker pruning by variance inflation factor – within each window the largest-VIF marker (VIF = the diagonal of the inverse marker correlation matrix) is dropped while any VIF exceeds vif, e.g. c(50, 5, 2) (Purcell et al. 2007).

blocks = TRUE

PLINK --blocks: define haplotype blocks by the method of Gabriel et al. (2002) – markers are grouped where the D' confidence intervals show strong LD and few recombination-consistent pairs – then keep one tag marker (highest MAF) per block. The block definitions reproduce PLINK 1.9's ⁠--blocks no-pheno-req⁠ (Haploview algorithm) block-for-block on the MAF >= 0.05 markers Haploview considers, at a matching span; block_max_kb (default 500) is PLINK's --blocks-max-kb (default 200). The tag-per-block collapse is a simplePHENOTYPES convenience on top of PLINK's block definitions.

The VIF pruner emits a one-per-session citation for PLINK (Purcell et al. 2007), and the block method one for Gabriel et al. (2002); the two pairwise pruners need no separate citation.

References

Purcell, S. et al. (2007). PLINK: a tool set for whole-genome association and population-based linkage analyses. Am. J. Hum. Genet. 81, 559–575. doi:10.1086/519795 (the VIF pruner).
Gabriel, S.B. et al. (2002). The structure of haplotype blocks in the human genome. Science 296, 2225–2229. doi:10.1126/science.1069424 (the haplotype-block definition).
Hill, W.G. and Robertson, A. (1968). Linkage disequilibrium in finite populations. Theor. Appl. Genet. 38, 226–231. doi:10.1007/BF01245622 (the phased r^2 used by indep_pairphase).
Gaunt, T.R., Rodriguez, S. and Day, I.N.M. (2007). Cubic exact solutions for the estimation of pairwise haplotype frequencies. BMC Bioinformatics 8, 428. doi:10.1186/1471-2105-8-428 (the exact cubic solution for two-locus haplotype frequencies behind the phased r^2 of indep_pairphase).

See Also

as_numeric(), simulate_phenotype(), dominance().

Examples

data("SNP55K_maize282_maf04")

# Common-variant, heterozygote-bearing markers (e.g. before a dominance model)
g <- filter_geno(SNP55K_maize282_maf04, maf_above = 0.1, hets = "include")

# Pairwise r^2 pruning in 50-marker windows
p <- filter_geno(SNP55K_maize282_maf04, indep_pairwise = c(50, 5, 0.2))

Description

Simulates founders with historical linkage disequilibrium and mutation-drift allele frequencies, as AlphaSimR's runMacs() does with MaCS, natively and reproducibly. Each chromosome is walked with a sequential Markov coalescent (SMC', Marjoram & Wall 2006) – the model AlphaSimR runs, since it calls MaCS with a one-base history window – under a piecewise-constant population-size history, with infinite-sites mutation. The haplotypes become a Population ready for cross(), selfcross(), double_haploid() and the selection engine.

Usage

founders_coalescent(
  n_ind,
  n_chr = 1,
  seg_sites = NULL,
  inbred = FALSE,
  species = c("GENERIC", "MAIZE", "WHEAT", "CATTLE"),
  split = NULL,
  theta = NULL,
  rho = NULL,
  history = NULL,
  morgans = NULL,
  bp = NULL,
  pool = NA_character_,
  seed = NULL
)

Arguments

Value

a Population of n_ind founders. Allele 1 (the counted allele) is the derived allele. Its attribute "coalescent" records the parameters, the per-chromosome seeds and the number of segregating sites simulated before subsetting (n_total).

Model and units

As in ms and MaCS, time is in units of $4N_0$ generations, theta = $4N_0\mu$ and rho = $4N_0 r$ are per chromosome, and history gives the population size relative to $N_0$ from each time back (ms -eN). A pair of lineages coalesces at rate $2/\lambda(u)$; along the chromosome, mutations and recombination points arrive at rates $\theta L$ and $\rho L$ for a marginal tree of total length $L$. The genetic map is linear in physical position, morgans long, and starts at 0 cM at the first segregating site (as runMacs() rebases it).

Species presets

species reproduces the demographic commands runMacs() builds in AlphaSimR 2.1.0 (read from its source; parameter values chosen by AlphaSimR, not estimates made here):

GENERIC

$10^8$ bp, $\theta = 1000$, $\rho = 400$, 1 Morgan, $N_e = 100$, sizes 5 / 15 / 60 / 120 / 1000 from times 0.25 / 2.5 / 25 / 250 / 2500.

MAIZE

$2 \times 10^8$ bp, $\theta = 1000$, $\rho = 800$, 2 Morgans, $N_e = 100$, a 15-step history.

WHEAT

$8 \times 10^8$ bp, $\theta = 320$, $\rho = 288$, 1.43 Morgans, $N_e = 50$, a 26-step history.

CATTLE

$2.8 \times 10^9 / 30$ bp per chromosome, mutation rate $9.4 \times 10^{-9}$ and recombination rate $9.26 \times 10^{-9}$ per bp, $N_e = 90$, a 12-step history in generations.

theta, rho, history, morgans and bp override the preset (the preset's $N_e$ still converts split to coalescent time).

References

Marjoram P, Wall JD (2006) Fast "coalescent" simulation. BMC Genet 7:16 (PMID 16539698). Chen GK, Marjoram P, Wall JD (2009) Fast and flexible simulation of DNA sequence data. Genome Res 19:136-142 (PMID 19029539). Hudson RR (2002) Generating samples under a Wright-Fisher neutral model of genetic variation. Bioinformatics 18:337-338 (PMID 11847089). Gaynor RC, Gorjanc G, Hickey JM (2021) AlphaSimR: an R package for breeding program simulations. G3 11(2):jkaa017. doi:10.1093/g3journal/jkaa017 (the runMacs() presets).

See Also

population_from_haplotypes(), as_population(), cross().

Examples

# 50 founders, 2 chromosomes, 200 segregating sites each (GENERIC preset)
f <- founders_coalescent(50, n_chr = 2, seg_sites = 200, seed = 1)
f
# inbred lines, e.g. for doubled-haploid programs
lines <- founders_coalescent(20, seg_sites = 100, inbred = TRUE, seed = 2)
all(dosages(lines) != 0)

Description

Builds the additive genomic relationship matrix $G = ZZ'/(2\sum p_j(1-p_j))$ from marker dosages, where $Z = M - 2p$ centres the 0/1/2 genotype matrix by twice the allele frequency (VanRaden 2008, method 1).

Usage

g_matrix(x, ridge = 0, base_freq = NULL)

Arguments

Details

By default p is the allele frequency of the current individuals (base_freq = NULL), and markers monomorphic in this set are dropped. Then G is a moving-base matrix: relatedness and the diagonal are expressed relative to this sample, so $F_i = G_{ii} - 1$ is inbreeding relative to the current set, not to a fixed founder base – across cycles of selection the base shifts, and a marker fixed by selection (p -> 0 or 1) drops out rather than reading as fixation. For a fixed reference (e.g. recurrent selection or founder-relative inbreeding) pass base_freq: markers polymorphic in the base are kept even when currently monomorphic, and ⁠2p⁠ centring uses the base frequencies.

Value

an individuals x individuals matrix with id dimnames.

References

VanRaden PM (2008) Efficient methods to compute genomic predictions. Journal of Dairy Science 91(11):4414–4423. doi:10.3168/jds.2007-0980

See Also

optimum_contribution(), dosages().

Examples

data("SNP55K_maize282_maf04")
pop <- as_population(SNP55K_maize282_maf04, individuals = 1:20)
G <- g_matrix(pop)
dim(G)
# Genomic inbreeding of the first few lines:
round(diag(G)[1:3] - 1, 3)

Description

The genetic component of each individual's phenotype: the sum of the scaled additive, dominance and epistatic layers, before the residual is added. This is the quantity the v1 engine wrote to Genetic_values.txt.

Usage

genetic_values(sim, rep = 1L)

Arguments

Details

When a phenotype includes a genome-derived transcriptome() layer, the genetic-mediated part of that expression component (the share tracing to the genome through gene expression) is also included, so it counts toward heritability. A real expression source contributes nothing here: its genetic content is not asserted, so it is treated as an environmental predictor and excluded from the genetic value.

Variance QTL are deliberately excluded: a vqtl() layer modulates the spread of the residual rather than contributing a genetic value, so it has no place in this matrix. By default this returns replication 1. When the simulation used vary_qtn = TRUE, use rep to retrieve the architecture and genetic values belonging to another replication.

Value

An individuals-by-traits numeric matrix, with individual IDs as row names and trait names as column names.

See Also

qtn_table() for the loci and effects that produce these values, phenotypes_long() for the phenotypes themselves.

Examples

data("SNP55K_maize282_maf04")
ph <- simulate_phenotype(SNP55K_maize282_maf04, n_traits = 2, h2 = 0.5,
                         seed = 1) |>
  additive(n_qtn = 5)

g <- genetic_values(ph)
head(g)

# Realized genetic correlation between the two traits
round(cor(g[, 1], g[, 2]), 3)

# Realized heritability: genetic variance over phenotypic variance
pheno <- phenotypes_wide(ph)
round(var(g[, 1]) / var(pheno$Trait_1), 3)

Description

The additive-plus-dominance counterpart of additive_value(). Scores each individual by its total genotypic value

$G_i = \sum_j \left[ a_j\,\mathrm{dosage}_{ij} + d_j\,\mathbf{1}(\mathrm{dosage}_{ij} = 0) \right]$

over a given, frozen architecture: a per-locus additive effect a and dominance deviation d, on the -1/0/1 dosage scale (so the heterozygote is dosage 0), with no per-population centring or rescaling. The per-locus term is -a / +d / +a for gene content 0 / 1 / 2 – exactly the raw genotypic value an additive(orthogonal = TRUE, a =, d =) layer and a dominance() layer build – but held fixed so it is comparable across generations, like additive_value().

Usage

genotypic_value(x, qtn, a, d)

Arguments

Details

This scores each individual's own per se total genotypic value $G = A + D$. It is not a parental combining ability: a per se genotypic value is not a parent's testcross merit – at a pure-dominance locus (a = 0, d > 0) a heterozygous Aa parent scores above a homozygous AA parent, yet crossed to an aa tester it is the AA parent whose progeny mean is higher, so per se scores can reverse testcross ranking. A reciprocal-recurrent or hybrid program therefore uses this by scoring the realized testcross / hybrid progeny (score the progeny population and dominance drives its mean), not the parents per se.

It is also not the transmissible breeding value. For selection on breeding value use select_ind() with on = "bv", which scores the average-effect breeding value $A_i = \sum_j \alpha_j (x_{ij} - 2 p_j)$ with $\alpha_j = a_j + d_j(1 - 2 p_j)$. Passing those $\alpha$ to additive_value() yields only a ranking-equivalent score ($\alpha_j\,\mathrm{dosage}_{ij}$, which differs from $A_i$ by an additive constant), not the population-centred breeding value.

Value

a named numeric vector of total genotypic values, one per individual.

References

Falconer DS, Mackay TFC (1996) Introduction to Quantitative Genetics, 4th ed. Longman, Harlow – the genotypic-value parameters $a$ (departure of the homozygote from the midpoint) and $d$ (dominance deviation of the heterozygote), and the average effect $\alpha = a + d(q - p)$.

See Also

additive_value(), phenotype_value(), genetic_values(), select_ind().

Examples

data("SNP55K_maize282_maf04")
pop <- as_population(SNP55K_maize282_maf04, individuals = 1:20)
# a = additive effects, d = dominance deviations (complete dominance at locus 5).
gv <- genotypic_value(pop, qtn = c(1, 5, 9),
                      a = c(0.5, -1, 2), d = c(0, 1, 0))
head(gv)

Description

Scores each individual's G x E slope (its response to the environmental covariate) on a frozen architecture:

$s_i = \mathrm{intercept} + \sum_j \mathrm{dosage}_{ij}\,b_j,$

with the -1/0/1 dosages and no per-population centring or rescaling, like additive_value(). This is the additive-by-environment trait of AlphaSimR (SimParam$addTraitAG()): its per-locus gxeEff are the effect here and its gxeInt the intercept, and AlphaSimR's centred genotype x - 1 (x = 0/1/2 copies of the counted allele) is this package's dosage. phenotype_value() adds s_i * w to the phenotype, where w is the environmental covariate of the trial (see its gxe / env / var_env arguments).

Usage

gxe_value(x, qtn, effect, intercept = 0)

Arguments

Value

a named numeric vector of slopes, one per individual.

References

Gaynor RC, Gorjanc G, Hickey JM (2021) AlphaSimR: an R package for breeding program simulations. G3 11(2):jkaa017. doi:10.1093/g3journal/jkaa017 (the additive-by-environment trait additive-by-environment trait). The exact slope scaling and the phenotype formula reproduced here are those of the AlphaSimR 2.1.0 source (SimParam$addTraitAG(), calcPheno()), not equations printed in the paper.

See Also

phenotype_value(), additive_value().

Examples

data("SNP55K_maize282_maf04")
pop <- as_population(SNP55K_maize282_maf04, individuals = 1:20)
s <- gxe_value(pop, qtn = c(1, 5, 9), effect = c(0.2, -0.1, 0.3))
head(s)

Description

The inverse of population_from_haplotypes(): the two homologous strands of every individual, 1 marking the counted (+1) allele (so that dosages(pop) == cis + trans - 1).

Usage

haplotypes(x)

Arguments

Value

A list with the integer matrices cis and trans, markers by individuals (rows named by map$snp, columns by the individual ids).

See Also

population_from_haplotypes(), dosages().

Examples

data("SNP55K_maize282_maf04")
pop <- as_population(SNP55K_maize282_maf04, individuals = 1:3)
h <- haplotypes(pop)
dim(h$cis)
stopifnot(identical(h$cis + h$trans - 1L, dosages(pop)))

Description

Realized heterosis of a crossbred population: its mean genotypic value minus the mean of the pure breeds weighted by its expected breed composition (breed_composition()) – the mid-parent for an F1. Also returns, for every pair of breeds, the expected F1 heterosis: the expected mean of random matings between the two breed populations minus their mid-parent, computed from the breeds' own genotypes (this package's derivation; it uses only the breeds' gamete frequencies, so linkage does not enter). Per locus, with gamete frequencies $p_A$, $p_B$ of the +a allele, the F1's allele frequency is the breeds' average, so the additive part cancels against the mid-parent and the heterosis is $d [h_{AB} - (h_A + h_B) / 2]$, where $h_{AB} = p_A (1 - p_B) + (1 - p_A) p_B$ is the F1's heterozygote frequency and $h_A$, $h_B$ the breeds' own. When each breed is in Hardy-Weinberg proportions ($h = 2p(1 - p)$) this is $\sum_j d_j (p_{Aj} - p_{Bj})^2$ – Falconer & Mackay's (1996) $H_{F1} = \sum d y^2$, with $y$ the difference in allele frequency between the two populations – zero for a purely additive architecture; for inbred-line breeds ($h = 0$) it is $\sum_j d_j (p_{Aj} + p_{Bj} - 2 p_{Aj} p_{Bj})$. The classical fractions of the F1 heterosis retained in later crosses (one half in an F2 or a backcross, two thirds in a two-breed rotation at equilibrium) follow from the probability that an individual's two alleles come from different breeds, a dominance (breed-origin) model, and are exact only under conditions on the breeds' Hardy-Weinberg deviations $h - 2p(1 - p)$ at the locus (not, in general, on each breed separately). A single fixed inbred line has no deviation ($p = 0$ or $1$, so $h = 0 = 2p(1 - p)$), but a breed made of several fully inbred lines that differ at a locus does ($h = 0 < 2p(1 - p)$). Exactly, at one locus with non-zero F1 heterosis the retained fraction in a backcross to breed A is one half if and only if $h_A = 2 p_A (1 - p_A)$, i.e. only the recurrent breed A need be in Hardy-Weinberg proportions (whatever breed B is), and in the F2 if and only if $[h_A - 2 p_A (1 - p_A)] + [h_B - 2 p_B (1 - p_B)] = 0$, i.e. the two deviations sum to zero (for example $-0.12$ and $+0.12$, with neither breed in Hardy-Weinberg proportions); for breeds that violate this the F2 and backcross retain a different fraction. (The two-thirds rotation figure is the classical value for breeds in Hardy-Weinberg proportions; its condition for breeds with deviations is not derived here.) For example, an equal mixture of AA and aa lines crossed to an aa line (pure dominance) has F1 heterozygosity 1/2 and backcross heterozygosity 1/2 against baselines of 0, so the retention is 1, not 1/2. realized is computed from the actual genotypes, so it does not assume either fraction; epistasis is outside this function.

Usage

heterosis(pop, breeds, qtn, a, d = 0)