blup: Convert breeding values into BLUP solutions of HapAllele effects

Description

Given solutions for a mixed linear model, compute Best Linear Unbiased Predictor (BLUP) solutions for HapAllele effects.

Usage

1
ghap.blup(blmm,haplo,weights=NULL,nperm=1,only.active.alleles = TRUE,ncores = 1)

Arguments

blmm

A GHap.blmm object, such as supplied by the ghap.blmm function.

haplo

A GHap.haplo object.

weights

A numeric vector providing HapAllele-specific weights.

nperm

Number of permutations to be performed for significance assessment (default = 1).

only.active.alleles

A logical value specifying whether only active haplotype alleles should be included in the calculations (default = TRUE).

ncores

A numeric value specifying the number of cores to be used in parallel computing (default = 1).

Details

The function uses the equation:

\mathbf{\hat{a}} = q\mathbf{DH}'\mathbf{K}^{-1}\mathbf{\hat{u}}

where \mathbf{H} is the n x m centered matrix of HapGenotypes observed for n individuals and m HapAlleles, \mathbf{D} = diag(d_i), d_i is the weight of HapAllele i (default d_i = 1), q is the inverse weighted sum of variances in the columns of \mathbf{H}, \mathbf{K} is the haplotype-based kinship matrix and \hat{u} is the vector of estimated breeding values. The permutation procedure consists in randomizing the vector \hat{\mathbf{u}} and computing the null statistic max(\mathbf{\hat{a}}). The permutation p-value is computed as the number of times the HapAllele effect was smaller than the null statistic.

Value

The function returns a dataframe with columns:

BLOCK

Block alias.

CHR

Chromosome name.

BP1

Block start position.

BP2

Block end position.

ALLELE

Haplotype allele identity.

SCORE

BLUP for the random effect of the haplotype allele.

FREQ

Frequency of the haplotype allele.

VAR

Variance in allele-specific breeding values.

pVAR

Proportion of variance explained by the haplotype allele.

CENTER

Average genotype (meaningful only for predictions with ghap.profile).

SCALE

A constant set to 1 (meaningful only for predictions with ghap.profile).

P

P-value for the permutation test. This column is suppressed if nperm < 1.

Author(s)

Yuri Tani Utsunomiya <ytutsunomiya@gmail.com>

References

I. Stranden and D.J. Garrick. Technical note: derivation of equivalent computing algorithms for genomic predictions and reliabilities of animal merit. J Dairy Sci. 2009. 92:2971-2975.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
# #### DO NOT RUN IF NOT NECESSARY ###
# 
# # Copy the example data in the current working directory
# ghap.makefile()
# 
# # Load data
# phase <- ghap.loadphase("human.samples", "human.markers", "human.phase")
# 
# # Subset data - randomly select 3000 markers with maf > 0.02
# maf <- ghap.maf(phase, ncores = 2)
# set.seed(1988)
# markers <- sample(phase$marker[maf > 0.02], 3000, replace = FALSE)
# phase <- ghap.subsetphase(phase, unique(phase$id), markers)
# rm(maf,markers)
# 
# # Generate block coordinates based on windows of 10 markers, sliding 5 marker at a time
# blocks <- ghap.blockgen(phase, 10, 5, "marker")
# 
# # Generate matrix of haplotype genotypes
# ghap.haplotyping(phase, blocks, batchsize = 100, ncores = 2, freq = 0.05, outfile = "example")
# 
# # Load haplotype genotypes
# haplo <- ghap.loadhaplo("example.hapsamples", "example.hapalleles", "example.hapgenotypes")
# 
# # Compute kinship matrix
# K <- ghap.kinship(haplo, batchsize = 100)
# 
# # Quantitative trait with 50% heritability
# # One major haplotype accounting for 30% of the genetic variance
# sim <- ghap.simpheno(haplo = haplo, K = K, h2 = 0.5, g2 = 0.3, major = 1000,seed=1988)
# 
# # Continuous model
# model <- ghap.blmm(fixed = phenotype ~ 1, random = "individual", data = sim$data, K = K)
# 
# ### RUN ###
# 
# 
# # BLUP solutions for HapAllele effects
# blup <- ghap.blup(blmm = model, haplo = haplo, ncores=2)
# blup$POS <- (blup$BP1+blup$BP2)/2e+6
# plot(blup$POS,100*blup$pVAR, xlab="Chromosome 2 position (in Mb)",
#      ylab="Haplotype variance explained (%)", pch=20, col="#471FAA99")
# abline(v=(haplo$bp1[1000]+haplo$bp2[1000])/2e+6,lty=3)
# 
# # Permutation test
# blup <- ghap.blup(blmm = model, haplo = haplo, ncores=2, nperm = 1000)
# blup$POS <- (blup$BP1+blup$BP2)/2e+6
# plot(blup$POS,-log10(blup$P), xlab="Chromosome 2 position (in Mb)",
#      ylab=expression(-log[10](p)), pch=20, col="#471FAA99")
# abline(v=(haplo$bp1[1000]+haplo$bp2[1000])/2e+6,lty=3)
# abline(h=-log10(0.05))
# 
# # Re-weighted solutions
# w <- nrow(blup)*blup$pVAR
# K <- ghap.kinship(haplo, weights = w, batchsize = 100)
# model2 <- ghap.mme(fixed = phenotype ~ 1, random = "individual", data = sim$data, K = K, 
#                    varcomp=c(model$vare,model$varu))
# blup <- ghap.blup(blmm = model2, haplo = haplo, weights = w, ncores=2)
# blup$POS <- (blup$BP1+blup$BP2)/2e+6
# plot(blup$POS,100*blup$pVAR, xlab="Chromosome 2 position (in Mb)",
#      ylab="Haplotype variance explained (%)", pch=20, col="#471FAA99")
# abline(v=(haplo$bp1[1000]+haplo$bp2[1000])/2e+6,lty=3)
# sim$major.effect
# blup$SCORE[1000]

Questions? Problems? Suggestions? or email at ian@mutexlabs.com.

All documentation is copyright its authors; we didn't write any of that.