|
0
|
1 #!/usr/bin/env Rscript
|
|
|
2
|
|
|
3 suppressPackageStartupMessages(library("optparse"))
|
|
26
|
4 suppressPackageStartupMessages(library("rjson"))
|
|
0
|
5
|
|
|
6 option_list <- list(
|
|
|
7 make_option(c("-c", "--components_input"), action="store", dest="components_input", help="Ks significant components input dataset"),
|
|
|
8 make_option(c("-k", "--kaks_input"), action="store", dest="kaks_input", help="KaKs analysis input dataset"),
|
|
9
|
9 make_option(c("-n", "--num_comp"), action="store", dest="num_comp", type="integer", help="Number of significant components in the Ks distribution"),
|
|
|
10 make_option(c("-o", "--output"), action="store", dest="output", help="Output dataset"),
|
|
23
|
11 make_option(c("-r", "--colors"), action="store", default=NA, help="List of component colors")
|
|
0
|
12 )
|
|
|
13
|
|
|
14 parser <- OptionParser(usage="%prog [options] file", option_list=option_list)
|
|
|
15 args <- parse_args(parser, positional_arguments=TRUE)
|
|
|
16 opt <- args$options
|
|
|
17
|
|
23
|
18 get_pi_mu_var = function(components_data, number_comp)
|
|
0
|
19 {
|
|
23
|
20 # FixMe: enhance this to generically handle any integer value for number_comp.
|
|
|
21 if (number_comp == 1)
|
|
0
|
22 {
|
|
4
|
23 pi <- c(components_data[1, 9])
|
|
|
24 mu <- c(components_data[1, 7])
|
|
|
25 var <- c(components_data[1, 8])
|
|
0
|
26 }
|
|
23
|
27 else if (number_comp == 2)
|
|
0
|
28 {
|
|
4
|
29 pi <- c(components_data[2, 9], components_data[3, 9])
|
|
|
30 mu <- c(components_data[2, 7], components_data[3, 7])
|
|
|
31 var <- c(components_data[2, 8], components_data[3, 8])
|
|
0
|
32 }
|
|
23
|
33 else if (number_comp == 3)
|
|
0
|
34 {
|
|
4
|
35 pi <- c(components_data[4, 9], components_data[5, 9], components_data[6, 9])
|
|
|
36 mu <- c(components_data[4, 7], components_data[5, 7], components_data[6, 7])
|
|
|
37 var <- c(components_data[4, 8], components_data[5, 8], components_data[6, 8])
|
|
0
|
38 }
|
|
23
|
39 else if (number_comp == 4)
|
|
0
|
40 {
|
|
4
|
41 pi <- c(components_data[7, 9], components_data[8, 9], components_data[9, 9], components_data[10, 9])
|
|
|
42 mu <- c(components_data[7, 7], components_data[8, 7], components_data[9, 7], components_data[10, 7])
|
|
|
43 var <- c(components_data[7, 8], components_data[8, 8], components_data[9, 8], components_data[10, 8])
|
|
0
|
44 }
|
|
23
|
45 else if (number_comp == 5)
|
|
4
|
46 {
|
|
|
47 pi <- c(components_data[11, 9], components_data[12, 9], components_data[13, 9], components_data[14, 9], components_data[15, 9])
|
|
|
48 mu <- c(components_data[11, 7], components_data[12, 7], components_data[13, 7], components_data[14, 7], components_data[15, 7])
|
|
|
49 var <- c(components_data[11, 8], components_data[12, 8], components_data[13, 8], components_data[14, 8], components_data[15, 8])
|
|
|
50 }
|
|
23
|
51 else if (number_comp == 6)
|
|
4
|
52 {
|
|
|
53 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])
|
|
|
54 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])
|
|
|
55 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])
|
|
|
56 }
|
|
|
57 results = c(pi, mu, var)
|
|
|
58 return(results)
|
|
0
|
59 }
|
|
|
60
|
|
23
|
61 plot_ks<-function(kaks_input, number_comp, colors, output, pi, mu, var)
|
|
0
|
62 {
|
|
|
63 # Start PDF device driver to save charts to output.
|
|
|
64 pdf(file=output, bg="white")
|
|
4
|
65 kaks <- read.table(file=kaks_input, header=T)
|
|
|
66 max_ks <- max(kaks$Ks, na.rm=TRUE)
|
|
0
|
67 # Change bin width
|
|
4
|
68 max_bin_range <- as.integer(max_ks / 0.05)
|
|
7
|
69 bin <- 0.05 * seq(0, (max_bin_range + 1 ))
|
|
|
70 kaks <- kaks[kaks$Ks<max_ks,]
|
|
4
|
71 h.kst <- hist(kaks$Ks, breaks=bin, plot=F)
|
|
|
72 nc <- h.kst$counts
|
|
|
73 vx <- h.kst$mids
|
|
|
74 ntot <- sum(nc)
|
|
0
|
75 # Set margin for plot bottom, left top, right.
|
|
4
|
76 par(mai=c(0.5, 0.5, 0, 0))
|
|
0
|
77 # Plot dimension in inches.
|
|
7
|
78 par(pin=c(3.0, 3.0))
|
|
|
79 g <- calculate_fitted_density(pi, mu, var, max_ks)
|
|
|
80 h <- ntot * 1.5 / sum(g)
|
|
|
81 vx <- seq(1, 100) * (max_ks / 100)
|
|
|
82 ymax <- max(nc)
|
|
4
|
83 barplot(nc, space=0.25, offset=0, width=0.04, xlim=c(0, max_ks), ylim=c(0, ymax), col="lightpink1", border="lightpink3")
|
|
0
|
84 # Add x-axis.
|
|
4
|
85 axis(1)
|
|
31
|
86 if (length(colors) == 0))
|
|
9
|
87 {
|
|
|
88 color <- c('red', 'yellow', 'green', 'black', 'blue', 'darkorange')
|
|
|
89 }
|
|
|
90 else
|
|
|
91 {
|
|
|
92 # Handle specified colors for components.
|
|
|
93 cStr <- unlist(colors)
|
|
|
94 color <- c()
|
|
|
95 items <- strsplit(cStr, ",")
|
|
|
96 for (item in items) {
|
|
|
97 color <- c(color, item)
|
|
|
98 }
|
|
|
99 num_colors_specified = length(color)
|
|
27
|
100 if (num_colors_specified < number_comp)
|
|
9
|
101 {
|
|
23
|
102 for (i in num_colors_specified:number_comp)
|
|
9
|
103 {
|
|
27
|
104 if (!any(color=='red'))
|
|
9
|
105 {
|
|
|
106 color <- c(color, 'red')
|
|
|
107 }
|
|
27
|
108 else if (!any(color=='yellow'))
|
|
9
|
109 {
|
|
|
110 color <- c(color, 'yellow')
|
|
|
111 }
|
|
27
|
112 else if (!any(color=='green'))
|
|
9
|
113 {
|
|
|
114 color <- c(color, 'green')
|
|
|
115 }
|
|
27
|
116 else if (!any(color=='black'))
|
|
9
|
117 {
|
|
|
118 color <- c(color, 'black')
|
|
|
119 }
|
|
27
|
120 else if (!any(color=='blue'))
|
|
9
|
121 {
|
|
|
122 color <- c(color, 'blue')
|
|
|
123 }
|
|
|
124 else
|
|
|
125 {
|
|
|
126 color <- c(color, 'darkorange')
|
|
|
127 }
|
|
|
128 }
|
|
|
129 }
|
|
|
130 }
|
|
0
|
131 for (i in 1:length(mu))
|
|
|
132 {
|
|
4
|
133 lines(vx, g[,i] * h, lwd=2, col=color[i])
|
|
0
|
134 }
|
|
4
|
135 }
|
|
0
|
136
|
|
7
|
137 calculate_fitted_density <- function(pi, mu, var, max_ks)
|
|
0
|
138 {
|
|
4
|
139 comp <- length(pi)
|
|
|
140 var <- var/mu^2
|
|
|
141 mu <- log(mu)
|
|
|
142 # Calculate lognormal density.
|
|
7
|
143 vx <- seq(1, 100) * (max_ks / 100)
|
|
4
|
144 fx <- matrix(0, 100, comp)
|
|
0
|
145 for (i in 1:100)
|
|
|
146 {
|
|
|
147 for (j in 1:comp)
|
|
|
148 {
|
|
4
|
149 fx[i, j] <- pi[j] * dlnorm(vx[i], meanlog=mu[j], sdlog=(sqrt(var[j])))
|
|
|
150 if (is.nan(fx[i,j])) fx[i,j]<-0
|
|
|
151 }
|
|
|
152 }
|
|
|
153 return(fx)
|
|
0
|
154 }
|
|
|
155
|
|
23
|
156 # Handle colors for components.
|
|
|
157 if (is.na(opt$colors))
|
|
28
|
158 {
|
|
23
|
159 # Randomly specify colors for components.
|
|
|
160 specified_colors <- c('red', 'yellow', 'green', 'black', 'blue', 'darkorange')
|
|
|
161 }
|
|
|
162 else
|
|
|
163 {
|
|
|
164 # Handle selected colors for components.
|
|
|
165 parser <- newJSONParser()
|
|
|
166 parser$addData(opt$colors)
|
|
|
167 raw_colors <- parser$getObject()
|
|
|
168 specified_colors <- c()
|
|
28
|
169 for (raw_color in raw_colors)
|
|
|
170 {
|
|
23
|
171 specified_colors <- c(specified_colors, raw_color)
|
|
28
|
172 }
|
|
23
|
173 }
|
|
|
174
|
|
0
|
175 # Read in the components data.
|
|
4
|
176 components_data <- read.delim(opt$components_input, header=TRUE)
|
|
23
|
177 number_comp <- opt$number_comp
|
|
0
|
178
|
|
|
179 # Set pi, mu, var.
|
|
23
|
180 items <- get_pi_mu_var(components_data, number_comp)
|
|
|
181 if (number_comp == 1)
|
|
7
|
182 {
|
|
|
183 pi <- items[1]
|
|
|
184 mu <- items[2]
|
|
|
185 var <- items[3]
|
|
|
186 }
|
|
30
|
187 else if (number_comp == 2)
|
|
7
|
188 {
|
|
|
189 pi <- items[1:2]
|
|
|
190 mu <- items[3:4]
|
|
|
191 var <- items[5:6]
|
|
|
192 }
|
|
30
|
193 else if (number_comp == 3)
|
|
7
|
194 {
|
|
|
195 pi <- items[1:3]
|
|
|
196 mu <- items[4:6]
|
|
|
197 var <- items[7:9]
|
|
|
198 }
|
|
30
|
199 else if (number_comp == 4)
|
|
7
|
200 {
|
|
|
201 pi <- items[1:4]
|
|
|
202 mu <- items[5:8]
|
|
|
203 var <- items[9:12]
|
|
|
204 }
|
|
30
|
205 else if (number_comp == 5)
|
|
7
|
206 {
|
|
|
207 pi <- items[1:5]
|
|
|
208 mu <- items[6:10]
|
|
|
209 var <- items[11:15]
|
|
|
210 }
|
|
30
|
211 else if (number_comp == 6)
|
|
7
|
212 {
|
|
|
213 pi <- items[1:6]
|
|
|
214 mu <- items[7:12]
|
|
|
215 var <- items[13:18]
|
|
|
216 }
|
|
0
|
217
|
|
|
218 # Plot the output.
|
|
23
|
219 plot_ks(opt$kaks_input, number_comp, specified_colors, opt$output, pi, mu, var)
|
