|
0
|
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 }
|
