CoMM

CoMM_S2 to dissecting genetic contributions to complex traits by leveraging regulatory information using GWAS summary data and eQTL individual level data

CoMM_S2(x, y, w, hatmu, hats, R, opts = NULL, px = 1)

Arguments

x	centered genotype (cis-SNPs) matrix for eQTL.
y	gene expression vector.
w	covariates file for eQTL data.
hatmu	coefficients from summary GWAS statistics
hats	standard deviance from summary GWAS statistics
R	correlation matrix for LD information
opts	is a list consists of "max_iter", "dispF", "display_gap", "epsStopLogLik", "fix_alphag". max_iter is an positive integer which denotes the maximal iterations, default value 1e5. dispF is a logic value either 1 or 0, controlling whether or not showing interations information, default value 1. display_gap is a positive integer which denotes the gap showing interations information, default value 10. epsStopLogLik is a positive value which control interation digits, default value 1e-5. constraintalpha is a logic value either 1 or 0, controlling whether alpha is constrained, default value 0. Noted that for separate analsis, fix_alphag does't exist in opts.
px	indicator for using PX-VBEM (default is 1). px = 1 for using PX-VBEM and VBEM otherwise.

Value

List of model parameters

Details

CoMM_S2 fits the CoMM_S2 model using GWAS summary data and eQTL individual level data.

Author

Jin Liu, jin.liu@duke-nus.edu.sg

Examples

set.seed(100)
L = 1;
M = 100;
rho =0.5;
n1 = 400;
n2 = 5000;
lam = 0.8;

X <- matrix(rnorm((n1+n2)*M),nrow=n1+n2,ncol=M);
beta_prop = 0.2;
b = matrix(0, M, 1); m = M * beta_prop;
b[sample(M,m)] = rnorm(m); h2y = 0.05;
y0 <- X%*%b + 6;
y  <- y0 + (as.vector(var(y0)*(1-h2y)/h2y))^0.5*rnorm(n1+n2);

h2 = 0.001;
y1 <- y[1:n1]
X1 <- X[1:n1,]
y2 <- y0[(n1+1):(n1+n2)]
X2 <- X[(n1+1):(n1+n2),]
alpha0 <- 3
alpha <- 0.3
sz2 <- var(y2*alpha) * ((1-h2)/h2)
z <- alpha0 + y2*alpha + rnorm(n2,0,sqrt(sz2))
y = y1;

mean.x1 = apply(X1,2,mean);
x1m = sweep(X1,2,mean.x1);
x1p = x1m;

mean.x2 = apply(X2,2,mean);
x2m = sweep(X2,2,mean.x2);
x2p = x2m;

w1 = matrix(rep(1,n1),ncol=1);


hatmu = matrix(0, M, 1)
hats = matrix(0, M, 1)

for (m in 1:M){
  fm = lm(z~1+x2p[,m]);
  hatmu[m] = summary(fm)$coefficients[2,1]
  hats[m] = summary(fm)$coefficients[2,2];
}

x3p = x1p;              
sumx3p = apply(x3p*x3p, 2, sum)
R = matrix(0, M, M);
for (i1 in 1:M){
  for (j1 in 1:M){
    R[i1,j1] = t(x3p[,i1])%*%x3p[,j1]/sqrt(sumx3p[i1]*sumx3p[j1])                        
  }
}
R = R*lam + (1 - lam)*diag(M)           

opts = list(max_iter = 10000, dispF = 1, display_gap = 10, epsStopLogLik = 1e-5, fix_alphag = 0);
px = 1

fm = CoMM_S2(x1p, y, w1, hatmu, hats, R, opts, px);
#> ***Iteration*******Fnew********Fold**********Diff***
#>    1.0000e+01  -1.4437e+03  -1.4437e+03   1.8819e-02
#> ***Iteration*******Fnew********Fold**********Diff***
#>    2.0000e+01  -1.4436e+03  -1.4437e+03   7.9145e-04
#> ***Iteration*******Fnew********Fold**********Diff***
#>    3.0000e+01  -1.4436e+03  -1.4436e+03   5.1346e-05