Mercurial > repos > dereeper > mlmm
annotate MLMM.pl @ 1:380b364980f9 draft default tip
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
author | dereeper |
---|---|
date | Mon, 16 Apr 2018 08:50:05 -0400 |
parents | |
children |
rev | line source |
---|---|
1
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
1 #!/usr/bin/perl |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
2 |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
3 use strict; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
4 use Getopt::Long; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
5 use Bio::SeqIO; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
6 |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
7 |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
8 my $usage = qq~Usage:$0 <args> [<opts>] |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
9 where <args> are: |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
10 -g, --geno <Genotype input> |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
11 -i, --info <SNP information. Genome position.> |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
12 -p, --pheno <Phenotype input> |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
13 -o, --out <output name> |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
14 -d, --directory <directory for MLMM R libraries> |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
15 -s, --step_number <number of steps. Maximum: 20. Default: 10> |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
16 -m, --method <Method: mbonf or extBIC. Default: mbonf> |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
17 ~; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
18 $usage .= "\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
19 |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
20 my ($geno,$map,$pheno,$out,$dir,$steps,$method); |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
21 |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
22 |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
23 GetOptions( |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
24 "geno=s" => \$geno, |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
25 "info=s" => \$map, |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
26 "pheno=s" => \$pheno, |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
27 "out=s" => \$out, |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
28 "dir=s" => \$dir, |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
29 "steps=s" => \$steps, |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
30 "method=s" => \$method |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
31 ); |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
32 |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
33 |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
34 die $usage |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
35 if ( !$geno || !$map || !$pheno || !$out || !$dir || !$steps || !$method); |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
36 |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
37 my $max_steps = 10; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
38 my $plot_opt = "mbonf"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
39 if ($method && $method ne 'mbonf' && $method ne 'extBIC') |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
40 { |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
41 print "Aborted: Method must be mbonf or extBIC.\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
42 exit; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
43 } |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
44 else |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
45 { |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
46 $plot_opt = $method; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
47 } |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
48 if ($steps && $steps !~/\d+/ && $steps > 20 && $steps < 2) |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
49 { |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
50 print "Aborted: Number of steps must be greater than 2 and lower than 20.\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
51 exit; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
52 } |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
53 else |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
54 { |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
55 $max_steps = $steps; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
56 } |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
57 |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
58 |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
59 my $chunk = 2; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
60 |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
61 my $RSCRIPT_EXE = "Rscript"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
62 my $R_DIR = $dir; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
63 |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
64 |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
65 my $head_trait = `head -1 $pheno`; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
66 my @headers_traits = split(/\t/,$head_trait); |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
67 my $trait_name = $headers_traits[1]; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
68 |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
69 |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
70 open( my $RCMD, ">rscript" ) or throw Error::Simple($!); |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
71 |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
72 |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
73 |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
74 |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
75 print $RCMD "Y_file <- \"" . $pheno . "\"\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
76 print $RCMD "X_file <- \"" . $geno . "\"\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
77 if($map) |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
78 { |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
79 print $RCMD "map_file <- \"$map\"\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
80 print $RCMD "map <- read.table(map_file, sep = \"\\t\", header = T)\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
81 } |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
82 print $RCMD "mlmm_data = list()\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
83 print $RCMD "mlmm_data\$chunk <- $chunk\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
84 print $RCMD "mlmm_data\$maxsteps <- $max_steps\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
85 print $RCMD "genot <- read.table(X_file, sep = \"\\t\", header = T)\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
86 print $RCMD "genot_mat <- as.matrix(genot[, 2:ncol(genot)])\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
87 print $RCMD "rownames(genot_mat) <- genot\$Ind_id\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
88 |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
89 print $RCMD "phenot <- read.table(Y_file, sep = \"\\t\", header = T)\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
90 |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
91 |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
92 |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
93 # missing data imputation |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
94 print $RCMD "genot_imp <- genot_mat\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
95 print $RCMD "average <- colMeans(genot_imp, na.rm = T)\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
96 print $RCMD "for (i in 1:ncol(genot_imp)){genot_imp[is.na(genot_imp[,i]), i] <- average[i]}\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
97 |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
98 # kinship matrix computation |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
99 print $RCMD "average <- colMeans(genot_imp, na.rm = T)\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
100 print $RCMD "stdev <- apply(genot_imp, 2, sd)\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
101 print $RCMD "genot_stand <- sweep(sweep(genot_imp, 2, average, \"-\"), 2, stdev, \"/\")\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
102 print $RCMD "K_mat <- (genot_stand %*% t(genot_stand)) / ncol(genot_stand)\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
103 print $RCMD "write.table(K_mat, '$out.kinship', sep='\\t', dec='.', quote=F, col.names=T, row.names=T)\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
104 |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
105 print $RCMD "source(\"" . $R_DIR. "/mlmm.r\")\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
106 print $RCMD "source(\"" . $R_DIR. "/emma.r\")\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
107 |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
108 # mlmm |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
109 print $RCMD "mygwas <- mlmm(Y = phenot\$$trait_name, X = genot_imp, K = K_mat, nbchunks=mlmm_data\$chunk, maxsteps=mlmm_data\$maxsteps)\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
110 |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
111 # plots |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
112 print $RCMD "pdf('$out.pdf')\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
113 print $RCMD "plot_step_table(mygwas, \"h2\")\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
114 print $RCMD "plot_step_table(mygwas, \"extBIC\")\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
115 print $RCMD "plot_step_table(mygwas, \"maxpval\")\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
116 print $RCMD "plot_step_RSS(mygwas)\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
117 # for (my $j = 1; $j <= ($max_steps - 1); $j++) |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
118 # { |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
119 # print $RCMD "plot_fwd_GWAS(mygwas, step = $j, snp_info = map, pval_filt = 0.1)\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
120 # } |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
121 print $RCMD "plot_opt_GWAS(mygwas, opt = \"extBIC\", snp_info = map, pval_filt = 0.1)\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
122 print $RCMD "plot_opt_GWAS(mygwas, opt = \"mbonf\", snp_info = map, pval_filt = 0.1)\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
123 #print $RCMD "qqplot_fwd_GWAS(mygwas, nsteps = mlmm_data\$maxsteps)\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
124 print $RCMD "qqplot_opt_GWAS(mygwas, opt = \"extBIC\")\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
125 print $RCMD "qqplot_opt_GWAS(mygwas, opt = \"mbonf\")\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
126 |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
127 # outputs |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
128 print $RCMD "write.table(mygwas\$RSSout, '$out.rss', sep='\\t', dec='.', quote=F, col.names=T, row.names=F)\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
129 print $RCMD "write.table(mygwas\$step_table, '$out.steptable', sep='\\t', dec='.', quote=F, col.names=T, row.names=F)\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
130 |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
131 $plot_opt = "\$opt_" . $plot_opt; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
132 print $RCMD "pval = mygwas" . $plot_opt . "\$out\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
133 print $RCMD "colnames(pval) = c(\"Marker_name\", \"Pvalue\")\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
134 print $RCMD "info_tmp = map\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
135 print $RCMD "colnames(info_tmp) = c(\"Marker_name\", \"Chr\", \"Pos\")\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
136 print $RCMD "res_asso = pval\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
137 print $RCMD qq~ |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
138 if(exists("info_tmp")){ |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
139 res_asso = merge(info_tmp, res_asso, by="Marker_name") |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
140 if( !is.element("Trait", colnames(info_tmp)) ){ |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
141 m = matrix(data="traitname", ncol=1, nrow=nrow(res_asso), dimnames=list(c(), c("Trait"))) |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
142 res_asso = cbind(m, res_asso) |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
143 } |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
144 } |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
145 ~; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
146 print $RCMD "res_asso = res_asso[order(res_asso[, \"Trait\"], res_asso[, \"Chr\"], res_asso[, \"Pos\"]), ]\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
147 print $RCMD "write.table(res_asso, '$out.res_asso', sep='\t', dec='.', quote=F, col.names=T, row.names=F)\n"; |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
148 close($RCMD); |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
149 |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
150 system("$RSCRIPT_EXE --vanilla rscript"); |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
151 |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
152 |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
153 |
380b364980f9
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
dereeper
parents:
diff
changeset
|
154 |