kinship: Relationship matrix based on genomic data

In GHap: Genome-Wide Haplotyping

Relationship matrix based on genomic data

Description

This function computes marker-based and HapAllele-based relationship matrices.

Usage

ghap.kinship(object, weights = NULL,
             sparsity = NULL,
             type = 1, batchsize = NULL,
             only.active.samples = TRUE,
             only.active.variants = TRUE,
             ncores = 1, verbose = TRUE)

Arguments

object	A valid GHap object (phase, haplo or plink).
weights	A numeric vector providing variant-specific weights.
sparsity	A numeric value specifying a relationship cut-off (default = NULL). All relationships below the specified cut-off will be set to zero, inducing sparsity into the relationship matrix.
type	A numeric value indicating the type of relationship matrix (see details).
batchsize	A numeric value controlling the number of variants to be processed at a time (default = nalleles/10).
only.active.samples	A logical value specifying whether only active samples should be included in the output (default = TRUE).
only.active.variants	A logical value specifying whether only active variants should be included in the output (default = TRUE).
ncores	A numeric value specifying the number of cores to be used in parallel computations (default = 1).
verbose	A logical value specfying whether log messages should be printed (default = TRUE).

Details

Let \mathbf{M} be the n x m matrix of genotypes, where n is the number of individuals and m is the number of variants (i.e, markers or HapAlleles). Prior to computation, genotypes in matrix \mathbf{M} are transformed according to the desired relationship type. After that transformation, the relationship matrix is computed as:

\mathbf{K} = q^{-1}\mathbf{MDM}'

where \mathbf{D} = diag(d_i), d_i is the weight of variant i (default d_i = 1), and q is a scaling factor. The argument type controls the genotype transformation and the scaling factor, and includes the following option:
Type = 1 (General additive 1)
Genotype transformation: m - mean(m)
Scaling factor: tr(\mathbf{MDM}')^{-1}n

Type = 2 (General additive 2)
Genotype transformation: (m - mean(m))/sd(m)
Scaling factor: m

Type = 3 (VanRaden, 2008)
Genotype transformation: m - 2*p[j]
Scaling factor: 2*sum(p*(1-p))

Type = 4 (GCTA)
Genotype transformation: (m - 2*p[j])/sqrt(2*p[j]*(1-p[j]))
Scaling factor: m

Type = 5 (Dominance 1)
Genotype transformation:
0 = -2*p[j]^2
1 = 2*p[j]*(1-p[j])
2 = -2*(1-p[j])^2
Scaling factor: 4*sum(p^2*(1-p)^2).

Type = 6 (Dominance 2)
Genotype transformation:
0 = -2*p[j]^2
1 = 2*p[j]*(1-p[j])
2 = -2*(1-p[j])^2
Scaling factor:
tr(\mathbf{MDM}')^{-1}n.

Value

The function returns a n x n relationship matrix, where n is the number of individuals.

Author(s)

Yuri Tani Utsunomiya <[e-mail]>
Marco Milanesi <[e-mail]>

References

P. M. VanRaden. Efficient methods to compute genomic predictions. J. Dairy. Sci. 2008. 91:4414-4423.

Examples

# #### DO NOT RUN IF NOT NECESSARY ###
# 
# # Copy plink data in the current working directory
# exfiles <- ghap.makefile(dataset = "example",
#                          format = "plink",
#                          verbose = TRUE)
# file.copy(from = exfiles, to = "./")
# 
# # Copy metadata in the current working directory
# exfiles <- ghap.makefile(dataset = "example",
#                          format = "meta",
#                          verbose = TRUE)
# file.copy(from = exfiles, to = "./")
# 
# # Load plink data
# plink <- ghap.loadplink("example")
# 
# # Load phenotype data
# df <- read.table(file = "example.phenotypes", header=T)
# 
# ### RUN ###
# 
# # Subset pure1 population
# pure1 <- plink$id[which(plink$pop == "Pure1")]
# plink <- ghap.subset(object = plink, ids = pure1, variants = plink$marker)
# 
# # Compute different types of relationship matrices
# K1 <- ghap.kinship(plink, type = 1) # General additive 1
# K2 <- ghap.kinship(plink, type = 2) # General additive 2
# K3 <- ghap.kinship(plink, type = 3) # VanRaden 2008
# K4 <- ghap.kinship(plink, type = 4) # GCTA GRM
# K5 <- ghap.kinship(plink, type = 5) # Dominance 1
# K6 <- ghap.kinship(plink, type = 6) # Dominance 2

GHap documentation built on July 2, 2022, 1:07 a.m.