|
0
|
1 library("MatrixEQTL");
|
|
|
2
|
|
|
3 # Number of columns (samples)
|
|
|
4 n = 100;
|
|
|
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 genotype and expression variables
|
|
|
16 snps.mat = cvrt.mat %*% rnorm(nc) + rnorm(n);
|
|
|
17 gene.mat = cvrt.mat %*% rnorm(nc) + rnorm(n) * noise.std + 0.5 * snps.mat + 1;
|
|
|
18
|
|
|
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,
|
|
|
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,
|
|
|
50 weights = 1/noise.std^2 );
|
|
|
51 lmout = summary(lmdl)$coefficients[2,c("Estimate","t value","Pr(>|t|)")];
|
|
|
52 print( lmout );
|
|
|
53 }
|
|
|
54
|
|
|
55 # Results from Matrix eQTL and "lm" must agree
|
|
|
56 stopifnot(all.equal(lmout, rez, check.attributes=FALSE));
|
