|
0
|
1 \name{modelLINEAR_CROSS}
|
|
|
2 \alias{modelLINEAR_CROSS}
|
|
|
3 \docType{data}
|
|
|
4 \title{
|
|
|
5 Constant for \code{\link{Matrix_eQTL_engine}}.
|
|
|
6 }
|
|
|
7 \description{
|
|
|
8 Set parameter \code{useModel = modelLINEAR_CROSS} in the call of \code{\link{Matrix_eQTL_main}} to indicate that Matrix eQTL should include the interaction of SNP and last covariate in the model and test for its significance.
|
|
|
9 }
|
|
|
10 \examples{
|
|
|
11 library('MatrixEQTL')
|
|
|
12
|
|
|
13 # Number of columns (samples)
|
|
|
14 n = 25;
|
|
|
15
|
|
|
16 # Number of covariates
|
|
|
17 nc = 10;
|
|
|
18
|
|
|
19 # Generate the standard deviation of the noise
|
|
|
20 noise.std = 0.1 + rnorm(n)^2;
|
|
|
21
|
|
|
22 # Generate the covariates
|
|
|
23 cvrt.mat = 2 + matrix(rnorm(n*nc), ncol = nc);
|
|
|
24
|
|
|
25 # Generate the vectors with single genotype and expression variables
|
|
|
26 snps.mat = cvrt.mat \%*\% rnorm(nc) + rnorm(n);
|
|
|
27 gene.mat = cvrt.mat \%*\% rnorm(nc) + rnorm(n) * noise.std +
|
|
|
28 1 + 0.5 * snps.mat + snps.mat * cvrt.mat[,nc];
|
|
|
29
|
|
|
30 # Create 3 SlicedData objects for the analysis
|
|
|
31 snps1 = SlicedData$new( matrix( snps.mat, nrow = 1 ) );
|
|
|
32 gene1 = SlicedData$new( matrix( gene.mat, nrow = 1 ) );
|
|
|
33 cvrt1 = SlicedData$new( t(cvrt.mat) );
|
|
|
34
|
|
|
35 # name of temporary output file
|
|
|
36 filename = tempfile();
|
|
|
37
|
|
|
38 # Call the main analysis function
|
|
|
39 me = Matrix_eQTL_main(
|
|
|
40 snps = snps1,
|
|
|
41 gene = gene1,
|
|
|
42 cvrt = cvrt1,
|
|
|
43 output_file_name = filename,
|
|
|
44 pvOutputThreshold = 1,
|
|
|
45 useModel = modelLINEAR_CROSS,
|
|
|
46 errorCovariance = diag(noise.std^2),
|
|
|
47 verbose = TRUE,
|
|
|
48 pvalue.hist = FALSE );
|
|
|
49 # remove the output file
|
|
|
50 unlink( filename );
|
|
|
51
|
|
|
52 # Pull Matrix eQTL results - t-statistic and p-value
|
|
|
53 beta = me$all$eqtls$beta;
|
|
|
54 tstat = me$all$eqtls$statistic;
|
|
|
55 pvalue = me$all$eqtls$pvalue;
|
|
|
56 rez = c(beta = beta, tstat = tstat, pvalue = pvalue)
|
|
|
57 # And compare to those from the linear regression in R
|
|
|
58 {
|
|
|
59 cat('\n\n Matrix eQTL: \n');
|
|
|
60 print(rez);
|
|
|
61 cat('\n R summary(lm()) output: \n')
|
|
|
62 lmodel = lm( gene.mat ~ snps.mat + cvrt.mat + snps.mat*cvrt.mat[,nc],
|
|
|
63 weights = 1/noise.std^2 );
|
|
|
64 lmout = tail(summary( lmodel )$coefficients,1)[,c(1,3,4)];
|
|
|
65 print( tail(lmout) );
|
|
|
66 }
|
|
|
67
|
|
|
68 # Results from Matrix eQTL and 'lm' must agree
|
|
|
69 stopifnot(all.equal(lmout, rez, check.attributes=FALSE))
|
|
|
70 }
|
|
|
71 \references{
|
|
|
72 The package website: \url{http://www.bios.unc.edu/research/genomic_software/Matrix_eQTL/}
|
|
|
73 }
|
|
|
74 \seealso{
|
|
|
75 See \code{\link{Matrix_eQTL_engine}} for reference and sample code.
|
|
|
76 }
|
