estimGenRel: Genomic relatedness

In timflutre/rutilstimflutre: Timothee Flutre's personal R code

Genomic relatedness

Description

Estimate genetic relationships between genotypes from their SNP genotypes. Note that this function estimates "relationships" and not "coancestries". See Toro et al (2011): "for diploid individuals, twice the coancestry coefficient is the additive relationship coefficient, which describes the ratio between the genetic covariance between individuals and the genetic variance of the base population".

Usage

estimGenRel(
  X,
  afs = NULL,
  thresh = NULL,
  relationships = "additive",
  method = "vanraden1",
  theta = 0.5,
  alpha = NULL,
  verbose = 1
)

Arguments

X	matrix of bi-allelic SNP genotypes encoded in allele doses in 0,1,2, with genotypes in rows and SNPs in columns; missing values should be encoded as NA
afs	vector of allele frequencies, corresponding to the alleles whose copies are counted in X (if NULL, will be calculated with estimSnpAf)
thresh	threshold on minor allele frequencies below which SNPs are ignored (e.g. 0.01; NULL to skip this step)
relationships	relationship to estimate (additive/dominant/gaussian) where "gaussian" corresponds to the Gaussian kernel from Endelman (2011)
method	if additive relationships, can be "noia" (see Vitezica et al, 2017), "vanraden1" (first method in VanRaden, 2008), "toro2011_eq10" (equation 10 using molecular covariance from Toro et al, 2011), "habier" (similar to "vanraden1" but without centering; from Habier et al, 2007), "astle-balding" (two times equation 2.2 in Astle & Balding, 2009), "yang" (similar to 'astle-balding' but without ignoring sampling error per SNP; from Yang et al, 2010), "zhou" (centering the genotypes with scale and not assuming that rare variants have larger effects; from Zhou et al, 2013) or "center-std"; if dominant relationships, can be "noia" (see Vitezica et al, 2017), "vitezica" (classical/statistical parametrization from Vitezica et al, 2013) or "su" (from Su et al, 2012)
theta	smoothing parameter for "gauss"
alpha	if not NULL, weight to make the matrix invertible as it often is singular because two genotypes are clones or because the centered-coding of SNP genotypes makes the last row predictable from the other ones (see Stranden and Christensen, 2011); otherwise use nearPD
verbose	verbosity level (0/1)

Value

symmetric matrix with the number of SNPs used as an attribute

Author(s)

Timothee Flutre

Examples

## Not run: set.seed(1859)
nb.genos <- 200
Ne <- 10^4
chrom.len <- 10^5
mu <- 10^(-8)
c <- 10^(-8)
genomes <- simulCoalescent(nb.inds=nb.genos,
                           pop.mut.rate=4 * Ne * mu * chrom.len,
                           pop.recomb.rate=4 * Ne * c * chrom.len,
                           chrom.len=chrom.len)
X <- genomes$genos

A.vr <- estimGenRel(X=X, relationships="additive", method="vanraden1")
mean(diag(A.vr)) # should be 1 in a base population at HWE
mean(A.vr[upper.tri(A.vr)]) # should be 0 in a base population at HWE

D.v <- estimGenRel(X=X, relationships="dominant", method="vitezica")
mean(diag(D.v)) # should be 1 in a base population at HWE
mean(D.v[upper.tri(D.v)]) # should be 0 in a base population at HWE

## End(Not run)

timflutre/rutilstimflutre documentation built on Feb. 12, 2025, 11:35 p.m.