Mercurial > repos > dereeper > blup_calculator
comparison blupcal_wrapper.R @ 0:45d215f2be74 draft
Uploaded
author | dereeper |
---|---|
date | Sat, 29 Dec 2018 18:44:05 -0500 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:45d215f2be74 |
---|---|
1 #!/usr/bin/Rscript | |
2 | |
3 options(show.error.messages = F, | |
4 error = function() { | |
5 cat(geterrmessage(), file = stderr()) | |
6 q("no", 1, F) } ) | |
7 loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8") | |
8 library("getopt") | |
9 | |
10 options(stringAsfactors = FALSE, useFancyQuotes = FALSE) | |
11 options(warn = -1) | |
12 args <- commandArgs(trailingOnly = TRUE) | |
13 option_specification <- matrix(c( | |
14 "tool_directory", "a", 1, "character" | |
15 ,"tabular_file", "b", 1, "character" | |
16 ,"replication_vector_column_index", "c", 1, "integer" | |
17 ,"genotype_vector_column_index", "d", 1, "integer" | |
18 ,"y_vector_column_index", "e", 1, "integer" | |
19 ,"design", "f", 1, "character" | |
20 ,"block_vector_column_index", "g", 1, "integer" | |
21 ,"summarize_by", "h", 1, "character" | |
22 ,"summarize_by_vector_column_index", "i", 1, "integer" | |
23 ,"group_variable_1", "j", 1, "character" | |
24 ,"group_variable_1_vector_column_index", "k", 1, "integer" | |
25 ,"group_variable_2", "l", 1, "character" | |
26 ,"group_variable_2_vector_column_index", "m", 1, "integer" | |
27 ,"output_file_path", "n", 1, "character" | |
28 ,"png_plots_file_path", "o", 1, "character" | |
29 ,"png_histograms_file_path", "p", 1, "character" | |
30 ,"pdf_plots_file_path", "q", 1, "character" | |
31 ,"pdf_histograms_file_path", "r", 1, "character" | |
32 ), byrow = TRUE, ncol = 4) | |
33 options <- getopt(option_specification) | |
34 | |
35 cat("\nTool Directory: ", options$tool_directory) | |
36 cat("\nTabular File: ", options$tabular_file) | |
37 cat("\nReplication Vector Column Index: ", options$replication_vector_column_index) | |
38 cat("\nGenotype Vector Column Index: ", options$genotype_vector_column_index) | |
39 cat("\nY Vector Column Index: ", options$y_vector_column_index) | |
40 cat("\nDesign: ", options$design) | |
41 cat("\nBlock Vector Column Index: ", options$block_vector_column_index) | |
42 cat("\nSummarize By: ", options$summarize_by) | |
43 cat("\nSummarize By Vector Column Index: ", options$summarize_by_vector_column_index) | |
44 cat("\nGroup Variable #1: ", options$group_variable_1) | |
45 cat("\nGroup Variable #1 Vector Column Index: ", options$group_variable_1_vector_column_index) | |
46 cat("\nGroup Variable #2: ", options$group_variable_2) | |
47 cat("\nGroup Variable #2 Vector Column Index: ", options$group_variable_2_vector_column_index) | |
48 cat("\nOutput file path: ", options$output_file_path) | |
49 cat("\nPNG plots file path: ", options$png_plots_file_path) | |
50 cat("\nPNG histograms file path: ", options$png_histograms_file_path) | |
51 cat("\nPDF plots file path: ", options$pdf_plots_file_path) | |
52 cat("\nPDF histograms file path: ", options$pdf_histograms_file_path) | |
53 cat("\n\n") | |
54 | |
55 tabular_data <- read.table(file = options$tabular_file, | |
56 header = TRUE, | |
57 sep = "\t", | |
58 stringsAsFactors = FALSE, | |
59 strip.white = TRUE, | |
60 na.strings = ".") | |
61 #tabular_data | |
62 | |
63 column_names <- colnames(tabular_data) | |
64 cat("Column names:\n") | |
65 column_names | |
66 cat("\n\n") | |
67 | |
68 # Replication | |
69 replication_vector_header <- column_names[options$replication_vector_column_index] | |
70 cat("\nreplication vector header: ", replication_vector_header) | |
71 | |
72 # Genotype | |
73 genotype_vector_header <- column_names[options$genotype_vector_column_index] | |
74 cat("\ngenotype vector header: ", genotype_vector_header) | |
75 | |
76 # Y | |
77 y_vector_header <- column_names[options$y_vector_column_index] | |
78 cat("\ny vector header: ", y_vector_header) | |
79 | |
80 # Design | |
81 design <- options$design | |
82 cat("\ndesign: ", design) | |
83 | |
84 # Block | |
85 if (design == "rcbd") { | |
86 block_vector_header <- NULL | |
87 } else { | |
88 block_vector_header <- column_names[options$block_vector_column_index] | |
89 } | |
90 cat("\nblock vector header: ", block_vector_header) | |
91 | |
92 # Summarize By | |
93 if (options$summarize_by == "false") { | |
94 summarize_by_vector_header <- NULL | |
95 } else { | |
96 summarize_by_vector_header <- column_names[options$summarize_by_vector_column_index] | |
97 } | |
98 cat("\nsummarize by vector header: ", summarize_by_vector_header) | |
99 | |
100 # Group Variable #1 | |
101 if (options$group_variable_1 == "false") { | |
102 group_variable_1_vector_header <- NULL | |
103 } else { | |
104 group_variable_1_vector_header <- column_names[options$group_variable_1_vector_column_index] | |
105 } | |
106 cat("\ngroup variable 1 vector header: ", group_variable_1_vector_header) | |
107 | |
108 # Group Variable #2 | |
109 if (options$group_variable_2 == "false") { | |
110 group_variable_2_vector_header <- NULL | |
111 } else { | |
112 group_variable_2_vector_header <- column_names[options$group_variable_2_vector_column_index] | |
113 } | |
114 cat("\ngroup variable 2 vector header: ", group_variable_2_vector_header) | |
115 | |
116 cat("\n\n") | |
117 | |
118 source(paste(options$tool_directory, "/blupcal.R", sep = "")) | |
119 | |
120 fit <- blupcal(data = tabular_data, | |
121 Replication = replication_vector_header, | |
122 Genotype = genotype_vector_header, | |
123 y = y_vector_header, | |
124 design = design, | |
125 Block = block_vector_header, | |
126 summarizeby = summarize_by_vector_header, | |
127 groupvar1 = group_variable_1_vector_header, | |
128 groupvar2 = group_variable_2_vector_header) | |
129 | |
130 cat("\n\n") | |
131 | |
132 output_file_path <- options$output_file_path | |
133 | |
134 cat("\nfit: ") | |
135 cat("\n") | |
136 fit | |
137 cat("\n") | |
138 | |
139 #fit$gid <- as.numeric(levels(fit$gid))[fit$gid] | |
140 if (options$summarize_by == "false") { | |
141 fit_view <- fit | |
142 } else { | |
143 fit_view <- fit[c(1, 7, 2, 3, 4, 5, 6)] | |
144 } | |
145 | |
146 formatted_table <- lapply(fit_view, function(.col) { | |
147 if (is.numeric(.col)) { | |
148 return(as.numeric(sprintf("%.3f", .col))) | |
149 } else { | |
150 return(.col) | |
151 } | |
152 }) | |
153 | |
154 write.table(formatted_table, | |
155 file = output_file_path, | |
156 sep = "\t", | |
157 row.names = FALSE, | |
158 quote = FALSE) | |
159 | |
160 png_plots_file_path <-options$png_plots_file_path | |
161 if (!(is.null(png_plots_file_path))) { | |
162 png(png_plots_file_path) | |
163 plot(fit) | |
164 dev.off() | |
165 } | |
166 | |
167 png_histograms_file_path <-options$png_histograms_file_path | |
168 if (!(is.null(png_histograms_file_path))) { | |
169 png(png_histograms_file_path) | |
170 plot2.blupcal(fit) | |
171 dev.off() | |
172 } | |
173 | |
174 pdf_plots_file_path <-options$pdf_plots_file_path | |
175 if (!(is.null(pdf_plots_file_path))) { | |
176 pdf(pdf_plots_file_path) | |
177 plot(fit) | |
178 dev.off() | |
179 } | |
180 | |
181 pdf_histograms_file_path <-options$pdf_histograms_file_path | |
182 if (!(is.null(pdf_histograms_file_path))) { | |
183 pdf(pdf_histograms_file_path) | |
184 plot2.blupcal(fit) | |
185 dev.off() | |
186 } |