annotate MatrixEQTL/demo/d.ANOVA5.r @ 0:cd4c8e4a4b5b draft

Uploaded
author jasonxu
date Fri, 12 Mar 2021 08:12:46 +0000
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
1 library("MatrixEQTL");
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
2
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
3 anova.groups = 5;
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
4 options(MatrixEQTL.ANOVA.categories = anova.groups);
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
5
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
6 # Number of columns (samples)
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
7 n = 100;
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
8 # Number of covariates
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
9 nc = 10;
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
10 # Generate the standard deviation of the noise
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
11 noise.std = 0.1 + rnorm(n)^2;
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
12 # Generate the covariates
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
13 cvrt.mat = 2 + matrix(rnorm(n*nc), ncol = nc);
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
14
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
15 # Generate the vectors with single genotype and expression variables
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
16 snps.mat = floor(runif(n, min = 0, max = anova.groups));
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
17 gene.mat = 1 + (snps.mat==1) + cvrt.mat %*% rnorm(nc) + rnorm(n) * noise.std;
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
18
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
19 # Create 3 SlicedData objects for the analysis
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
20 snps1 = SlicedData$new( matrix( snps.mat, nrow = 1 ) );
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
21 gene1 = SlicedData$new( matrix( gene.mat, nrow = 1 ) );
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
22 cvrt1 = SlicedData$new( t(cvrt.mat) );
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
23
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
24 # Produce no output files
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
25 filename = NULL; # tempfile()
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
26
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
27 # Call the main analysis function
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
28 me = Matrix_eQTL_main(
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
29 snps = snps1,
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
30 gene = gene1,
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
31 cvrt = cvrt1,
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
32 output_file_name = filename,
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
33 pvOutputThreshold = 1,
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
34 useModel = modelANOVA,
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
35 errorCovariance = diag(noise.std^2),
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
36 verbose = TRUE,
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
37 pvalue.hist = FALSE );
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
38
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
39 # Pull Matrix eQTL results - t-statistic and p-value
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
40
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
41 fstat = me$all$eqtls$statistic;
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
42 pvalue = me$all$eqtls$pvalue;
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
43 rez = c( Fstat = fstat, pvalue = pvalue);
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
44 # And compare to those from ANOVA in R
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
45 {
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
46 cat("\n\n Matrix eQTL: \n");
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
47 print(rez);
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
48 cat("\n R anova(lm()) output: \n");
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
49 lmdl = lm( gene.mat ~ cvrt.mat + factor(snps.mat),
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
50 weights = 1/noise.std^2 );
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
51 lmout = anova(lmdl)[2, c("F value","Pr(>F)")];
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
52 print( lmout );
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
53 }
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
54
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
55 # Results from Matrix eQTL and "lm" must agree
cd4c8e4a4b5b Uploaded
jasonxu
parents:
diff changeset
56 stopifnot(all.equal(lmout, rez, check.attributes=FALSE));