|
0
|
1 #!/usr/bin/env Rscript
|
|
|
2
|
|
|
3 suppressPackageStartupMessages(library("optparse"))
|
|
|
4
|
|
|
5 option_list <- list(
|
|
|
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"),
|
|
1
|
8 make_option(c("-n", "--number_comp"), action="store", dest="number_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"),
|
|
|
10 make_option(c("-s", "--specified_colors"), action="store", dest="specified_colors", default=NULL, help="List of component colors")
|
|
0
|
11 )
|
|
|
12
|
|
|
13 parser <- OptionParser(usage="%prog [options] file", option_list=option_list)
|
|
|
14 args <- parse_args(parser, positional_arguments=TRUE)
|
|
|
15 opt <- args$options
|
|
|
16
|
|
1
|
17 set_component_colors = function(colors, number_comp) {
|
|
|
18 # Handle colors for components.
|
|
|
19 if (is.null(colors)) {
|
|
|
20 # Randomly specify colors for components.
|
|
|
21 component_colors <- c("red", "yellow", "green", "black", "blue", "darkorange")
|
|
|
22 } else {
|
|
|
23 # Handle selected colors for components.
|
|
|
24 component_colors <- c()
|
|
|
25 colors <- as.character(colors)
|
|
|
26 items <- strsplit(colors, ",")
|
|
|
27 for (item in items) {
|
|
|
28 component_colors <- c(component_colors, item)
|
|
|
29 }
|
|
|
30 num_colors_specified <- length(component_colors)
|
|
|
31 if (num_colors_specified < number_comp) {
|
|
|
32 # The number of selected colors is less than the number of
|
|
|
33 # components, so we'll add random colors that were not
|
|
|
34 # selected to the set of component colors until we have a
|
|
|
35 # color for each component.
|
|
|
36 loop_count <- number_comp - num_colors_specified
|
|
|
37 for (i in 1:loop_count) {
|
|
|
38 if (!(is.element("red", component_colors))) {
|
|
|
39 component_colors <- c(component_colors, "red")
|
|
|
40 } else if (!(is.element("yellow", component_colors))) {
|
|
|
41 component_colors <- c(component_colors, "yellow")
|
|
|
42 } else if (!(is.element("green", component_colors))) {
|
|
|
43 component_colors <- c(component_colors, "green")
|
|
|
44 } else if (!(is.element("black", component_colors))) {
|
|
|
45 component_colors <- c(component_colors, "black")
|
|
|
46 } else if (!(is.element("blue", component_colors))) {
|
|
|
47 component_colors <- c(component_colors, "blue")
|
|
|
48 } else if (!(is.element("darkorange", component_colors))) {
|
|
|
49 component_colors <- c(component_colors, "darkorange")
|
|
|
50 }
|
|
|
51 }
|
|
|
52 }
|
|
|
53 }
|
|
|
54 return(component_colors)
|
|
0
|
55 }
|
|
|
56
|
|
1
|
57 get_pi_mu_var = function(components_data, number_comp) {
|
|
|
58 if (number_comp == 1) {
|
|
0
|
59 pi <- c(components_data[1, 9])
|
|
|
60 mu <- c(components_data[1, 7])
|
|
|
61 var <- c(components_data[1, 8])
|
|
1
|
62 } else if (number_comp == 2) {
|
|
0
|
63 pi <- c(components_data[2, 9], components_data[3, 9])
|
|
|
64 mu <- c(components_data[2, 7], components_data[3, 7])
|
|
|
65 var <- c(components_data[2, 8], components_data[3, 8])
|
|
1
|
66 } else if (number_comp == 3) {
|
|
0
|
67 pi <- c(components_data[4, 9], components_data[5, 9], components_data[6, 9])
|
|
|
68 mu <- c(components_data[4, 7], components_data[5, 7], components_data[6, 7])
|
|
|
69 var <- c(components_data[4, 8], components_data[5, 8], components_data[6, 8])
|
|
1
|
70 } else if (number_comp == 4) {
|
|
0
|
71 pi <- c(components_data[7, 9], components_data[8, 9], components_data[9, 9], components_data[10, 9])
|
|
|
72 mu <- c(components_data[7, 7], components_data[8, 7], components_data[9, 7], components_data[10, 7])
|
|
|
73 var <- c(components_data[7, 8], components_data[8, 8], components_data[9, 8], components_data[10, 8])
|
|
1
|
74 } else if (number_comp == 5) {
|
|
0
|
75 pi <- c(components_data[11, 9], components_data[12, 9], components_data[13, 9], components_data[14, 9], components_data[15, 9])
|
|
|
76 mu <- c(components_data[11, 7], components_data[12, 7], components_data[13, 7], components_data[14, 7], components_data[15, 7])
|
|
|
77 var <- c(components_data[11, 8], components_data[12, 8], components_data[13, 8], components_data[14, 8], components_data[15, 8])
|
|
1
|
78 } else if (number_comp == 6) {
|
|
0
|
79 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])
|
|
|
80 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])
|
|
|
81 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])
|
|
|
82 }
|
|
1
|
83 results <- c(pi, mu, var)
|
|
0
|
84 return(results)
|
|
|
85 }
|
|
|
86
|
|
1
|
87 plot_ks<-function(kaks_input, component_colors, output, pi, mu, var) {
|
|
0
|
88 # Start PDF device driver to save charts to output.
|
|
|
89 pdf(file=output, bg="white")
|
|
|
90 kaks <- read.table(file=kaks_input, header=T)
|
|
|
91 max_ks <- max(kaks$Ks, na.rm=TRUE)
|
|
|
92 # Change bin width
|
|
|
93 max_bin_range <- as.integer(max_ks / 0.05)
|
|
|
94 bin <- 0.05 * seq(0, (max_bin_range + 1 ))
|
|
|
95 kaks <- kaks[kaks$Ks<max_ks,]
|
|
|
96 h.kst <- hist(kaks$Ks, breaks=bin, plot=F)
|
|
|
97 nc <- h.kst$counts
|
|
|
98 vx <- h.kst$mids
|
|
|
99 ntot <- sum(nc)
|
|
|
100 # Set margin for plot bottom, left top, right.
|
|
|
101 par(mai=c(0.5, 0.5, 0, 0))
|
|
|
102 # Plot dimension in inches.
|
|
|
103 par(pin=c(3.0, 3.0))
|
|
|
104 g <- calculate_fitted_density(pi, mu, var, max_ks)
|
|
|
105 h <- ntot * 1.5 / sum(g)
|
|
|
106 vx <- seq(1, 100) * (max_ks / 100)
|
|
|
107 ymax <- max(nc)
|
|
|
108 barplot(nc, space=0.25, offset=0, width=0.04, xlim=c(0, max_ks), ylim=c(0, ymax), col="lightpink1", border="lightpink3")
|
|
|
109 # Add x-axis.
|
|
|
110 axis(1)
|
|
1
|
111 for (i in 1:length(mu)) {
|
|
|
112 lines(vx, g[,i] * h, lwd=2, col=component_colors[i])
|
|
0
|
113 }
|
|
|
114 }
|
|
|
115
|
|
1
|
116 calculate_fitted_density <- function(pi, mu, var, max_ks) {
|
|
0
|
117 comp <- length(pi)
|
|
|
118 var <- var/mu^2
|
|
|
119 mu <- log(mu)
|
|
|
120 # Calculate lognormal density.
|
|
|
121 vx <- seq(1, 100) * (max_ks / 100)
|
|
|
122 fx <- matrix(0, 100, comp)
|
|
1
|
123 for (i in 1:100) {
|
|
|
124 for (j in 1:comp) {
|
|
|
125 fx[i, j] <- pi[j] * dlnorm(vx[i], meanlog=mu[j], sdlog=(sqrt(var[j])))
|
|
|
126 if (is.nan(fx[i,j])) {
|
|
|
127 fx[i,j]<-0
|
|
|
128 }
|
|
0
|
129 }
|
|
|
130 }
|
|
|
131 return(fx)
|
|
|
132 }
|
|
|
133
|
|
|
134 # Read in the components data.
|
|
|
135 components_data <- read.delim(opt$components_input, header=TRUE)
|
|
1
|
136 number_comp <- as.integer(opt$number_comp)
|
|
|
137 if (number_comp == 0) {
|
|
|
138 # Default to 1 component to allow functional testing.
|
|
|
139 number_comp <- 1
|
|
|
140 }
|
|
|
141
|
|
|
142 # Set component colors.
|
|
|
143 component_colors <- set_component_colors(opt$specified_colors, number_comp)
|
|
0
|
144
|
|
|
145 # Set pi, mu, var.
|
|
1
|
146 items <- get_pi_mu_var(components_data, number_comp)
|
|
|
147 if (number_comp == 1) {
|
|
|
148 pi <- items[1]
|
|
|
149 mu <- items[2]
|
|
|
150 var <- items[3]
|
|
|
151 } else if (number_comp == 2) {
|
|
|
152 pi <- items[1:2]
|
|
|
153 mu <- items[3:4]
|
|
|
154 var <- items[5:6]
|
|
|
155 } else if (number_comp == 3) {
|
|
|
156 pi <- items[1:3]
|
|
|
157 mu <- items[4:6]
|
|
|
158 var <- items[7:9]
|
|
|
159 } else if (number_comp == 4) {
|
|
|
160 pi <- items[1:4]
|
|
|
161 mu <- items[5:8]
|
|
|
162 var <- items[9:12]
|
|
|
163 } else if (number_comp == 5) {
|
|
|
164 pi <- items[1:5]
|
|
|
165 mu <- items[6:10]
|
|
|
166 var <- items[11:15]
|
|
|
167 } else if (number_comp == 6) {
|
|
|
168 pi <- items[1:6]
|
|
|
169 mu <- items[7:12]
|
|
|
170 var <- items[13:18]
|
|
0
|
171 }
|
|
|
172
|
|
|
173 # Plot the output.
|
|
1
|
174 plot_ks(opt$kaks_input, component_colors, opt$output, pi, mu, var)
|
