annotate MatrixEQTL/demo/sample.all.r @ 4:cf4e9238644c draft default tip

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