Mercurial > repos > jasonxu > matrixeqtl
diff MatrixEQTL/demo/a.nocvrt.r @ 3:ae74f8fb3aef draft
Uploaded
author | jasonxu |
---|---|
date | Fri, 12 Mar 2021 08:20:57 +0000 |
parents | cd4c8e4a4b5b |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/MatrixEQTL/demo/a.nocvrt.r Fri Mar 12 08:20:57 2021 +0000 @@ -0,0 +1,56 @@ +library("MatrixEQTL"); + +# Number of columns (samples) +n = 100; + + + + + + + + + + +# Generate the vectors with genotype and expression variables +snps.mat = rnorm(n); +gene.mat = rnorm(n) + 0.5 * snps.mat; + +# Create 3 SlicedData objects for the analysis +snps1 = SlicedData$new( matrix( snps.mat, nrow = 1 ) ); +gene1 = SlicedData$new( matrix( gene.mat, nrow = 1 ) ); +cvrt1 = SlicedData$new(); + +# Produce no output files +filename = NULL; # tempfile() + +# Call the main analysis function +me = Matrix_eQTL_main( + snps = snps1, + gene = gene1, + cvrt = cvrt1, + output_file_name = filename, + pvOutputThreshold = 1, + useModel = modelLINEAR, + errorCovariance = numeric(), + verbose = TRUE, + pvalue.hist = FALSE ); + +# Pull Matrix eQTL results - t-statistic and p-value +beta = me$all$eqtls$beta; +tstat = me$all$eqtls$statistic; +pvalue = me$all$eqtls$pvalue; +rez = c(beta = beta, tstat = tstat, pvalue = pvalue); +# And compare to those from the linear regression in R +{ + cat("\n\n Matrix eQTL: \n"); + print(rez); + cat("\n R summary(lm()) output: \n"); + lmdl = lm( gene.mat ~ snps.mat ); + + lmout = summary(lmdl)$coefficients[2,c("Estimate","t value","Pr(>|t|)")]; + print( lmout ); +} + +# Results from Matrix eQTL and "lm" must agree +stopifnot(all.equal(lmout, rez, check.attributes=FALSE));