# blup: Convert breeding values into BLUP solutions of HapAllele... In GHap: Genome-Wide Haplotyping

### 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 <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 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 19, 2017, 8:43 p.m.

Search within the GHap package
Search all R packages, documentation and source code