coher: coheritability estimation

Description Usage Arguments Details Value Author(s) References Examples

View source: R/coher.R

Description

coheritability estimation

Usage

1
coher(Y, X, a.quantile = 0.25, alpha.enet = 1, para = TRUE)

Arguments

Y

a named list of phenotypes

X

genotypes matrix

a.quantile

a threshold to remove week correlation covariates, only use with large data set.

alpha.enet

alpha value in elastic net

para

parallel options: TRUE or FALSE

Details

na

Value

coher.matrix

matrix of coheritabilities

coefficients.B

list of estimated coefficients from elastic-net for each phenotype

Author(s)

The Tien Mai

References

to be filled

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
##---- Should be DIRECTLY executable !! ----
##-- ==>  Define data, use random,
##--	or do  help(data=index)  for the standard data sets.

## The function is currently defined as
function (Y, X, a.quantile = 0.25, alpha.enet = 1, para = TRUE)
{
    library(parallel)
    require(doMC)
    registerDoMC(cores = detectCores() - 1)
    library(glmnet)
    q = length(Y)
    p = ncol(X)
    bhat = list()
    nameSNP = colnames(X)
    nameSAMPLe = rownames(X)
    selected.snps = c()
    for (j in 1:q) {
        matSNP = X[names(Y[[j]]), ]
        if (a.quantile > 0) {
            sam.cor <- apply(matSNP, 2, function(snp) cor.scale(snp,
                Y[[j]]))
            sam.cor[is.na(sam.cor)] = 0
            xsele = matSNP[, sam.cor > quantile(sam.cor, a.quantile)]
            sk = Matrix(xsele, sparse = T)
        }
        else if (a.quantile == 0) {
            sk = Matrix(X, sparse = T)
        }
        if (checkBinaryTrait(Y[[j]]) == "binomial") {
            enet <- cv.glmnet(x = sk, y = Y[[j]], parallel = para,
                alpha = alpha.enet, family = "binomial")
        }
        else {
            enet <- cv.glmnet(x = sk, y = Y[[j]], parallel = para,
                alpha = alpha.enet)
        }
        tam = coef(enet, s = "lambda.min")[-1]
        names(tam) = colnames(xsele)
        bhat[[j]] = rep(0, p)
        names(bhat[[j]]) = nameSNP
        bhat[[j]][names(tam)] <- tam
        selected.snps = c(selected.snps, names(tam[tam != 0]))
    }
    selected.snps = unique(selected.snps)
    covariance.selesnps = cov2(X[, selected.snps])
    mat.coher = matrix(0, nr = q, nc = q, dimnames = list(names(Y),
        names(Y)))
    for (i in 1:(q - 1)) {
        for (j in (i + 1):q) {
            mat.coher[i, j] = bhat[[i]][selected.snps] %*% covariance.selesnps %*%
                bhat[[j]][selected.snps]
            mat.coher[j, i] = bhat[[i]][selected.snps] %*% covariance.selesnps %*%
                bhat[[j]][selected.snps]/sqrt(bhat[[i]][selected.snps] %*%
                covariance.selesnps %*% bhat[[i]][selected.snps] *
                bhat[[j]][selected.snps] %*% covariance.selesnps %*%
                  bhat[[j]][selected.snps])
        }
    }
    return(coher.matrix = mat.coher, coefficients.B = bhat)
  }

tienmt/coher documentation built on Dec. 15, 2021, 6:59 a.m.