Mercurial > repos > jasonxu > matrixeqtl

comparison MatrixEQTL/demo/e.interaction.r @ 0:cd4c8e4a4b5b draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author jasonxu
date Fri, 12 Mar 2021 08:12:46 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:cd4c8e4a4b5b
1 library("MatrixEQTL");
2
3 # Number of columns (samples)
4 n = 25;
5
6 # Number of covariates
7 nc = 10;
8
9 # Generate the standard deviation of the noise
10 noise.std = 0.1 + rnorm(n)^2;
11
12 # Generate the covariates
13 cvrt.mat = 2 + matrix(rnorm(n*nc), ncol = nc);
14
15 # Generate the vectors with single genotype and expression variables
16 snps.mat = cvrt.mat %*% rnorm(nc) + rnorm(n);
17 gene.mat = cvrt.mat %*% rnorm(nc) + rnorm(n) * noise.std +
18 1 + 0.5 * snps.mat + snps.mat * cvrt.mat[,nc];
19 # Create 3 SlicedData objects for the analysis
20 snps1 = SlicedData$new( matrix( snps.mat, nrow = 1 ) );
21 gene1 = SlicedData$new( matrix( gene.mat, nrow = 1 ) );
22 cvrt1 = SlicedData$new( t(cvrt.mat) );
23
24 # Produce no output files
25 filename = NULL; # tempfile()
26
27 # Call the main analysis function
28 me = Matrix_eQTL_main(
29 snps = snps1,
30 gene = gene1,
31 cvrt = cvrt1,
32 output_file_name = filename,
33 pvOutputThreshold = 1,
34 useModel = modelLINEAR_CROSS,
35 errorCovariance = diag(noise.std^2),
36 verbose = TRUE,
37 pvalue.hist = FALSE );
38
39 # Pull Matrix eQTL results - t-statistic and p-value
40 beta = me$all$eqtls$beta;
41 tstat = me$all$eqtls$statistic;
42 pvalue = me$all$eqtls$pvalue;
43 rez = c(beta = beta, tstat = tstat, pvalue = pvalue);
44 # And compare to those from the linear regression in R
45 {
46 cat("\n\n Matrix eQTL: \n");
47 print(rez);
48 cat("\n R summary(lm()) output: \n");
49 lmdl = lm( gene.mat ~ snps.mat + cvrt.mat + snps.mat*cvrt.mat[,nc],
50 weights = 1/noise.std^2 );
51 lmout = tail(summary(lmdl)$coefficients,1)[,c(1,3,4)];
52 print( tail(lmout) );
53 }
54
55 # Results from Matrix eQTL and "lm" must agree
56 stopifnot(all.equal(lmout, rez, check.attributes=FALSE));