Mercurial > repos > greg > ks_distribution

comparison ks_distribution.R @ 23:492ef68d73f7 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Mon, 26 Jun 2017 09:56:24 -0400
parents e6bbc7426cd1
children da8615a486fe
comparison
equal deleted inserted replaced
22:f3b7455d6c04 23:492ef68d73f7
5 option_list <- list( 5 option_list <- list(
6 make_option(c("-c", "--components_input"), action="store", dest="components_input", help="Ks significant components input dataset"), 6 make_option(c("-c", "--components_input"), action="store", dest="components_input", help="Ks significant components input dataset"),
7 make_option(c("-k", "--kaks_input"), action="store", dest="kaks_input", help="KaKs analysis input dataset"), 7 make_option(c("-k", "--kaks_input"), action="store", dest="kaks_input", help="KaKs analysis input dataset"),
8 make_option(c("-n", "--num_comp"), action="store", dest="num_comp", type="integer", help="Number of significant components in the Ks distribution"), 8 make_option(c("-n", "--num_comp"), action="store", dest="num_comp", type="integer", help="Number of significant components in the Ks distribution"),
9 make_option(c("-o", "--output"), action="store", dest="output", help="Output dataset"), 9 make_option(c("-o", "--output"), action="store", dest="output", help="Output dataset"),
10 make_option(c("-r", "--colors"), action="store", default=[], help="List of component colors") 10 make_option(c("-r", "--colors"), action="store", default=NA, help="List of component colors")
11 ) 11 )
12 12
13 parser <- OptionParser(usage="%prog [options] file", option_list=option_list) 13 parser <- OptionParser(usage="%prog [options] file", option_list=option_list)
14 args <- parse_args(parser, positional_arguments=TRUE) 14 args <- parse_args(parser, positional_arguments=TRUE)
15 opt <- args$options 15 opt <- args$options
16 16
17 get_pi_mu_var = function(components_data, num_components) 17 get_pi_mu_var = function(components_data, number_comp)
18 { 18 {
19 # FixMe: enhance this to generically handle any integer value for num_components. 19 # FixMe: enhance this to generically handle any integer value for number_comp.
20 if (num_components == 1) 20 if (number_comp == 1)
21 { 21 {
22 pi <- c(components_data[1, 9]) 22 pi <- c(components_data[1, 9])
23 mu <- c(components_data[1, 7]) 23 mu <- c(components_data[1, 7])
24 var <- c(components_data[1, 8]) 24 var <- c(components_data[1, 8])
25 } 25 }
26 else if (num_components == 2) 26 else if (number_comp == 2)
27 { 27 {
28 pi <- c(components_data[2, 9], components_data[3, 9]) 28 pi <- c(components_data[2, 9], components_data[3, 9])
29 mu <- c(components_data[2, 7], components_data[3, 7]) 29 mu <- c(components_data[2, 7], components_data[3, 7])
30 var <- c(components_data[2, 8], components_data[3, 8]) 30 var <- c(components_data[2, 8], components_data[3, 8])
31 } 31 }
32 else if (num_components == 3) 32 else if (number_comp == 3)
33 { 33 {
34 pi <- c(components_data[4, 9], components_data[5, 9], components_data[6, 9]) 34 pi <- c(components_data[4, 9], components_data[5, 9], components_data[6, 9])
35 mu <- c(components_data[4, 7], components_data[5, 7], components_data[6, 7]) 35 mu <- c(components_data[4, 7], components_data[5, 7], components_data[6, 7])
36 var <- c(components_data[4, 8], components_data[5, 8], components_data[6, 8]) 36 var <- c(components_data[4, 8], components_data[5, 8], components_data[6, 8])
37 } 37 }
38 else if (num_components == 4) 38 else if (number_comp == 4)
39 { 39 {
40 pi <- c(components_data[7, 9], components_data[8, 9], components_data[9, 9], components_data[10, 9]) 40 pi <- c(components_data[7, 9], components_data[8, 9], components_data[9, 9], components_data[10, 9])
41 mu <- c(components_data[7, 7], components_data[8, 7], components_data[9, 7], components_data[10, 7]) 41 mu <- c(components_data[7, 7], components_data[8, 7], components_data[9, 7], components_data[10, 7])
42 var <- c(components_data[7, 8], components_data[8, 8], components_data[9, 8], components_data[10, 8]) 42 var <- c(components_data[7, 8], components_data[8, 8], components_data[9, 8], components_data[10, 8])
43 } 43 }
44 else if (num_components == 5) 44 else if (number_comp == 5)
45 { 45 {
46 pi <- c(components_data[11, 9], components_data[12, 9], components_data[13, 9], components_data[14, 9], components_data[15, 9]) 46 pi <- c(components_data[11, 9], components_data[12, 9], components_data[13, 9], components_data[14, 9], components_data[15, 9])
47 mu <- c(components_data[11, 7], components_data[12, 7], components_data[13, 7], components_data[14, 7], components_data[15, 7]) 47 mu <- c(components_data[11, 7], components_data[12, 7], components_data[13, 7], components_data[14, 7], components_data[15, 7])
48 var <- c(components_data[11, 8], components_data[12, 8], components_data[13, 8], components_data[14, 8], components_data[15, 8]) 48 var <- c(components_data[11, 8], components_data[12, 8], components_data[13, 8], components_data[14, 8], components_data[15, 8])
49 } 49 }
50 else if (num_components == 6) 50 else if (number_comp == 6)
51 { 51 {
52 pi <- c(components_data[16, 9], components_data[17, 9], components_data[18, 9], components_data[19, 9], components_data[20, 9], components_data[21, 9]) 52 pi <- c(components_data[16, 9], components_data[17, 9], components_data[18, 9], components_data[19, 9], components_data[20, 9], components_data[21, 9])
53 mu <- c(components_data[16, 7], components_data[17, 7], components_data[18, 7], components_data[19, 7], components_data[20, 7], components_data[21, 7]) 53 mu <- c(components_data[16, 7], components_data[17, 7], components_data[18, 7], components_data[19, 7], components_data[20, 7], components_data[21, 7])
54 var <- c(components_data[16, 8], components_data[17, 8], components_data[18, 8], components_data[19, 8], components_data[20, 8], components_data[21, 8]) 54 var <- c(components_data[16, 8], components_data[17, 8], components_data[18, 8], components_data[19, 8], components_data[20, 8], components_data[21, 8])
55 } 55 }
56 results = c(pi, mu, var) 56 results = c(pi, mu, var)
57 return(results) 57 return(results)
58 } 58 }
59 59
60 plot_ks<-function(kaks_input, num_components, colors, output, pi, mu, var) 60 plot_ks<-function(kaks_input, number_comp, colors, output, pi, mu, var)
61 { 61 {
62 # Start PDF device driver to save charts to output. 62 # Start PDF device driver to save charts to output.
63 pdf(file=output, bg="white") 63 pdf(file=output, bg="white")
64 kaks <- read.table(file=kaks_input, header=T) 64 kaks <- read.table(file=kaks_input, header=T)
65 max_ks <- max(kaks$Ks, na.rm=TRUE) 65 max_ks <- max(kaks$Ks, na.rm=TRUE)
94 items <- strsplit(cStr, ",") 94 items <- strsplit(cStr, ",")
95 for (item in items) { 95 for (item in items) {
96 color <- c(color, item) 96 color <- c(color, item)
97 } 97 }
98 num_colors_specified = length(color) 98 num_colors_specified = length(color)
99 if (num_colors_specified < num_components): 99 if (num_colors_specified < number_comp):
100 { 100 {
101 for (i in num_colors_specified:num_components) 101 for (i in num_colors_specified:number_comp)
102 { 102 {
103 if !(any(color=='red')) 103 if !(any(color=='red'))
104 { 104 {
105 color <- c(color, 'red') 105 color <- c(color, 'red')
106 } 106 }
150 } 150 }
151 } 151 }
152 return(fx) 152 return(fx)
153 } 153 }
154 154
155 # Handle colors for components.
156 if (is.na(opt$colors))
157 {
158 # Randomly specify colors for components.
159 specified_colors <- c('red', 'yellow', 'green', 'black', 'blue', 'darkorange')
160 }
161 else
162 {
163 # Handle selected colors for components.
164 parser <- newJSONParser()
165 parser$addData(opt$colors)
166 raw_colors <- parser$getObject()
167 specified_colors <- c()
168 for (raw_color in raw_colors) {
169 specified_colors <- c(specified_colors, raw_color)
170 }
171
155 # Read in the components data. 172 # Read in the components data.
156 components_data <- read.delim(opt$components_input, header=TRUE) 173 components_data <- read.delim(opt$components_input, header=TRUE)
157 num_components <- opt$num_comp 174 number_comp <- opt$number_comp
158 175
159 # Set pi, mu, var. 176 # Set pi, mu, var.
160 items <- get_pi_mu_var(components_data, num_components) 177 items <- get_pi_mu_var(components_data, number_comp)
161 if (num_components == 1) 178 if (number_comp == 1)
162 { 179 {
163 pi <- items[1] 180 pi <- items[1]
164 mu <- items[2] 181 mu <- items[2]
165 var <- items[3] 182 var <- items[3]
166 } 183 }
167 if (num_components == 2) 184 if (number_comp == 2)
168 { 185 {
169 pi <- items[1:2] 186 pi <- items[1:2]
170 mu <- items[3:4] 187 mu <- items[3:4]
171 var <- items[5:6] 188 var <- items[5:6]
172 } 189 }
173 if (num_components == 3) 190 if (number_comp == 3)
174 { 191 {
175 pi <- items[1:3] 192 pi <- items[1:3]
176 mu <- items[4:6] 193 mu <- items[4:6]
177 var <- items[7:9] 194 var <- items[7:9]
178 } 195 }
179 if (num_components == 4) 196 if (number_comp == 4)
180 { 197 {
181 pi <- items[1:4] 198 pi <- items[1:4]
182 mu <- items[5:8] 199 mu <- items[5:8]
183 var <- items[9:12] 200 var <- items[9:12]
184 } 201 }
185 if (num_components == 5) 202 if (number_comp == 5)
186 { 203 {
187 pi <- items[1:5] 204 pi <- items[1:5]
188 mu <- items[6:10] 205 mu <- items[6:10]
189 var <- items[11:15] 206 var <- items[11:15]
190 } 207 }
191 if (num_components == 6) 208 if (number_comp == 6)
192 { 209 {
193 pi <- items[1:6] 210 pi <- items[1:6]
194 mu <- items[7:12] 211 mu <- items[7:12]
195 var <- items[13:18] 212 var <- items[13:18]
196 } 213 }
197 214
198 # Plot the output. 215 # Plot the output.
199 plot_ks(opt$kaks_input, num_components, opt$colors, opt$output, pi, mu, var) 216 plot_ks(opt$kaks_input, number_comp, specified_colors, opt$output, pi, mu, var)