blup: Convert breeding values into BLUP solutions of HapAllele...

Description Usage Arguments Details Value Author(s) References Examples

Description

Given genomic estimated breeding values (GEBVs), compute Best Linear Unbiased Predictor (BLUP) solutions for HapAllele effects.

Usage

1
2
ghap.blup(gebvs,haplo,invcov, gebvsweights = NULL, haploweights = NULL, nperm = 1,
          only.active.alleles = TRUE, ncores = 1)

Arguments

gebvs

A vector of GEBVs. The vector must be named and all names must be present in the GHap.haplo object.

haplo

A GHap.haplo object.

invcov

The inverse covariance (i.e., genomic kinship) matrix for GEBVs.

gebvsweights

A numeric vector providing individual-specific weights.

haploweights

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{DM}'\mathbf{K}^{-1}\mathbf{\hat{u}}

where \mathbf{M} is the N x H centered matrix of HapGenotypes observed for N individuals and H 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{M}, \mathbf{K} is the haplotype-based kinship matrix and \hat{u} is the vector of GEBVs. 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 <[email protected]>

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
65
66
67
68
69
70
71
72
73
74
# #### 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 - markers with maf > 0.05
# maf <- ghap.maf(phase, ncores = 2)
# markers <- phase$marker[maf > 0.05]
# phase <- ghap.subsetphase(phase, unique(phase$id), markers)
# 
# # Generate blocks of 5 markers sliding 5 markers at a time
# blocks.mkr <- ghap.blockgen(phase, windowsize = 5, slide = 5, unit = "marker")
#
# # Generate matrix of haplotype genotypes
# ghap.haplotyping(phase, blocks.mkr, batchsize = 100, ncores = 2, outfile = "human")
#
# # Load haplotype genotypes
# haplo <- ghap.loadhaplo("human.hapsamples", "human.hapalleles", "human.hapgenotypes")
#
#
# ### RUN ###
# # Subset common haplotypes in Europeans
# EUR.ids <- haplo$id[haplo$pop %in% c("TSI","CEU")]
# haplo <- ghap.subsethaplo(haplo,EUR.ids,rep(TRUE,times=haplo$nalleles))
# hapstats <- ghap.hapstats(haplo, ncores = 2)
# common <- hapstats$TYPE %in% c("REGULAR","MAJOR") &
#  hapstats$FREQ > 0.05 &
#  hapstats$FREQ < 0.95
# haplo <- ghap.subsethaplo(haplo,EUR.ids,common)
#
# #Compute relationship matrix
# K <- ghap.kinship(haplo, batchsize = 100)
# 
# # Quantitative trait with 50% heritability
# # Unbalanced repeated measurements (0 to 30)
# # Two major haplotypes accounting for 50% of the genetic variance
# myseed <- 123456789
# set.seed(myseed)
# major <- sample(which(haplo$allele.in == TRUE),size = 2)
# g2 <- runif(n = 2, min = 0, max = 1)
# g2 <- (g2/sum(g2))*0.5
# sim <- ghap.simpheno(haplo, kinship = K, h2 = 0.5, g2 = g2, nrep = 30,
#                      balanced = FALSE, major = major, seed = myseed)
# 
# #Fit model using REML
# model <- ghap.lmm(fixed = phenotype ~ 1, random = ~ individual,
#                   covmat = list(individual = K), data = sim$data)
# 
#
# ### RUN ###
# 
# #BLUP GWAS
# gebvs <- model$random$individual
# gebvsw <- table(sim$data$individual)
# gebvsw <- gebvsw + mean(gebvsw)
# gebvsw <- gebvsw[names(gebvs)]
# Kinv <- ghap.kinv(K)
# gwas.blup <- ghap.blup(gebvs = gebvs, haplo = haplo, gebvsweights = gebvsw,
#                        ncores = 4, invcov = Kinv)
# plot(gwas.blup$BP1/1e+6,gwas.blup$pVAR*100,pch=20,
#      xlab="Position (in Mb)",ylab="Variance explained (%)")
# abline(v=haplo$bp1[major]/1e+6)
# #BLUP with one update
# w <- gwas.blup$VAR*nrow(gwas.blup)
# K2 <- ghap.kinship(haplo=haplo,weights = w)
# Kinv2 <- ghap.kinv(K2)
# gwas.blup2 <- ghap.blup(gebvs = gebvs, haplo = haplo, invcov = Kinv2, ncores = 2,
#                         gebvsweights = gebvsw, haploweights = w)
# plot(gwas.blup2$BP1/1e+6,gwas.blup2$pVAR*100,pch=20,
#      xlab="Position (in Mb)",ylab="Variance explained (%)")
# abline(v=haplo$bp1[major]/1e+6)

GHap documentation built on May 29, 2017, 9:56 p.m.