Mercurial > repos > jasonxu > matrixeqtl

comparison MatrixEQTL/demo/sample.all.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 # Matrix eQTL by Andrey A. Shabalin
2 # http://www.bios.unc.edu/research/genomic_software/Matrix_eQTL/
3 #
4 # Be sure to use an up to date version of R and Matrix eQTL.
5
6 # source("Matrix_eQTL_R/Matrix_eQTL_engine.r");
7 library(MatrixEQTL)
8
9 ## Location of the package with the data files.
10 base.dir = find.package('MatrixEQTL');
11 # base.dir = '.';
12
13 ## Settings
14
15 # Linear model to use, modelANOVA, modelLINEAR, or modelLINEAR_CROSS
16 useModel = modelLINEAR; # modelANOVA, modelLINEAR, or modelLINEAR_CROSS
17
18 # Genotype file name
19 SNP_file_name = paste(base.dir, "/data/SNP.txt", sep="");
20
21 # Gene expression file name
22 expression_file_name = paste(base.dir, "/data/GE.txt", sep="");
23
24 # Covariates file name
25 # Set to character() for no covariates
26 covariates_file_name = paste(base.dir, "/data/Covariates.txt", sep="");
27
28 # Output file name
29 output_file_name = tempfile();
30
31 # Only associations significant at this level will be saved
32 pvOutputThreshold = 1e-2;
33
34 # Error covariance matrix
35 # Set to numeric() for identity.
36 errorCovariance = numeric();
37 # errorCovariance = read.table("Sample_Data/errorCovariance.txt");
38
39
40 ## Load genotype data
41
42 snps = SlicedData$new();
43 snps$fileDelimiter = "\t"; # the TAB character
44 snps$fileOmitCharacters = "NA"; # denote missing values;
45 snps$fileSkipRows = 1; # one row of column labels
46 snps$fileSkipColumns = 1; # one column of row labels
47 snps$fileSliceSize = 2000; # read file in slices of 2,000 rows
48 snps$LoadFile(SNP_file_name);
49
50 ## Load gene expression data
51
52 gene = SlicedData$new();
53 gene$fileDelimiter = "\t"; # the TAB character
54 gene$fileOmitCharacters = "NA"; # denote missing values;
55 gene$fileSkipRows = 1; # one row of column labels
56 gene$fileSkipColumns = 1; # one column of row labels
57 gene$fileSliceSize = 2000; # read file in slices of 2,000 rows
58 gene$LoadFile(expression_file_name);
59
60 ## Load covariates
61
62 cvrt = SlicedData$new();
63 cvrt$fileDelimiter = "\t"; # the TAB character
64 cvrt$fileOmitCharacters = "NA"; # denote missing values;
65 cvrt$fileSkipRows = 1; # one row of column labels
66 cvrt$fileSkipColumns = 1; # one column of row labels
67 if(length(covariates_file_name)>0) {
68 cvrt$LoadFile(covariates_file_name);
69 }
70
71 ## Run the analysis
72
73 me = Matrix_eQTL_engine(
74 snps = snps,
75 gene = gene,
76 cvrt = cvrt,
77 output_file_name = output_file_name,
78 pvOutputThreshold = pvOutputThreshold,
79 useModel = useModel,
80 errorCovariance = errorCovariance,
81 verbose = TRUE,
82 pvalue.hist = TRUE,
83 min.pv.by.genesnp = FALSE,
84 noFDRsaveMemory = FALSE);
85
86 unlink(output_file_name);
87
88 ## Results:
89
90 cat('Analysis done in: ', me$time.in.sec, ' seconds', '\n');
91 cat('Detected eQTLs:', '\n');
92 show(me$all$eqtls)
93
94 ## Plot the histogram of all p-values
95
96 plot(me)