Mercurial > repos > jasonxu > matrixeqtl

comparison MatrixEQTL/demo/sample.cis.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 snps_location_file_name = paste(base.dir, "/data/snpsloc.txt", sep="");
21
22 # Gene expression file name
23 expression_file_name = paste(base.dir, "/data/GE.txt", sep="");
24 gene_location_file_name = paste(base.dir, "/data/geneloc.txt", sep="");
25
26 # Covariates file name
27 # Set to character() for no covariates
28 covariates_file_name = paste(base.dir, "/data/Covariates.txt", sep="");
29
30 # Output file name
31 output_file_name_cis = tempfile();
32 output_file_name_tra = tempfile();
33
34 # Only associations significant at this level will be saved
35 pvOutputThreshold_cis = 2e-2;
36 pvOutputThreshold_tra = 1e-2;
37
38 # Error covariance matrix
39 # Set to numeric() for identity.
40 errorCovariance = numeric();
41 # errorCovariance = read.table("Sample_Data/errorCovariance.txt");
42
43 # Distance for local gene-SNP pairs
44 cisDist = 1e6;
45
46 ## Load genotype data
47
48 snps = SlicedData$new();
49 snps$fileDelimiter = "\t"; # the TAB character
50 snps$fileOmitCharacters = "NA"; # denote missing values;
51 snps$fileSkipRows = 1; # one row of column labels
52 snps$fileSkipColumns = 1; # one column of row labels
53 snps$fileSliceSize = 2000; # read file in slices of 2,000 rows
54 snps$LoadFile(SNP_file_name);
55
56 ## Load gene expression data
57
58 gene = SlicedData$new();
59 gene$fileDelimiter = "\t"; # the TAB character
60 gene$fileOmitCharacters = "NA"; # denote missing values;
61 gene$fileSkipRows = 1; # one row of column labels
62 gene$fileSkipColumns = 1; # one column of row labels
63 gene$fileSliceSize = 2000; # read file in slices of 2,000 rows
64 gene$LoadFile(expression_file_name);
65
66 ## Load covariates
67
68 cvrt = SlicedData$new();
69 cvrt$fileDelimiter = "\t"; # the TAB character
70 cvrt$fileOmitCharacters = "NA"; # denote missing values;
71 cvrt$fileSkipRows = 1; # one row of column labels
72 cvrt$fileSkipColumns = 1; # one column of row labels
73 if(length(covariates_file_name)>0) {
74 cvrt$LoadFile(covariates_file_name);
75 }
76
77 ## Run the analysis
78 snpspos = read.table(snps_location_file_name, header = TRUE, stringsAsFactors = FALSE);
79 genepos = read.table(gene_location_file_name, header = TRUE, stringsAsFactors = FALSE);
80
81 me = Matrix_eQTL_main(
82 snps = snps,
83 gene = gene,
84 cvrt = cvrt,
85 output_file_name = output_file_name_tra,
86 pvOutputThreshold = pvOutputThreshold_tra,
87 useModel = useModel,
88 errorCovariance = errorCovariance,
89 verbose = TRUE,
90 output_file_name.cis = output_file_name_cis,
91 pvOutputThreshold.cis = pvOutputThreshold_cis,
92 snpspos = snpspos,
93 genepos = genepos,
94 cisDist = cisDist,
95 pvalue.hist = TRUE,
96 min.pv.by.genesnp = TRUE,
97 noFDRsaveMemory = FALSE);
98
99 unlink(output_file_name_tra);
100 unlink(output_file_name_cis);
101
102 ## Results:
103
104 cat('Analysis done in: ', me$time.in.sec, ' seconds', '\n');
105 cat('Detected local eQTLs:', '\n');
106 show(me$cis$eqtls)
107 cat('Detected distant eQTLs:', '\n');
108 show(me$trans$eqtls)
109
110 ## Plot the histogram of local and distant p-values
111
112 plot(me)