Mercurial > repos > greg > ks_distribution

comparison ks_distribution.R @ 9:214e2710c51e draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Fri, 23 Jun 2017 14:12:59 -0400
parents 22cae2172406
children 1f6943662833
comparison
equal deleted inserted replaced
8:1650842a90ba 9:214e2710c51e
3 suppressPackageStartupMessages(library("optparse")) 3 suppressPackageStartupMessages(library("optparse"))
4 4
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("-o", "--output"), action="store", dest="output", help="Output dataset") 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"),
10 make_option(c("-r", "--colors"), action="store", default=NA, help="List of component colors"),
9 ) 11 )
10 12
11 parser <- OptionParser(usage="%prog [options] file", option_list=option_list) 13 parser <- OptionParser(usage="%prog [options] file", option_list=option_list)
12 args <- parse_args(parser, positional_arguments=TRUE) 14 args <- parse_args(parser, positional_arguments=TRUE)
13 opt <- args$options 15 opt <- args$options
14
15
16 get_num_components = function(components_data)
17 {
18 # Get the max of the number_comp column.
19 number_comp = components_data[, 3]
20 num_components <- max(number_comp, na.rm=TRUE)
21 return(num_components)
22 }
23 16
24 get_pi_mu_var = function(components_data, num_components) 17 get_pi_mu_var = function(components_data, num_components)
25 { 18 {
26 # FixMe: enhance this to generically handle any integer value for num_components. 19 # FixMe: enhance this to generically handle any integer value for num_components.
27 if (num_components == 1) 20 if (num_components == 1)
62 } 55 }
63 results = c(pi, mu, var) 56 results = c(pi, mu, var)
64 return(results) 57 return(results)
65 } 58 }
66 59
67 plot_ks<-function(kaks_input, output, pi, mu, var) 60 plot_ks<-function(kaks_input, num_components, colors, output, pi, mu, var)
68 { 61 {
69 # Start PDF device driver to save charts to output. 62 # Start PDF device driver to save charts to output.
70 pdf(file=output, bg="white") 63 pdf(file=output, bg="white")
71 kaks <- read.table(file=kaks_input, header=T) 64 kaks <- read.table(file=kaks_input, header=T)
72 max_ks <- max(kaks$Ks, na.rm=TRUE) 65 max_ks <- max(kaks$Ks, na.rm=TRUE)
87 vx <- seq(1, 100) * (max_ks / 100) 80 vx <- seq(1, 100) * (max_ks / 100)
88 ymax <- max(nc) 81 ymax <- max(nc)
89 barplot(nc, space=0.25, offset=0, width=0.04, xlim=c(0, max_ks), ylim=c(0, ymax), col="lightpink1", border="lightpink3") 82 barplot(nc, space=0.25, offset=0, width=0.04, xlim=c(0, max_ks), ylim=c(0, ymax), col="lightpink1", border="lightpink3")
90 # Add x-axis. 83 # Add x-axis.
91 axis(1) 84 axis(1)
92 color <- c('red', 'yellow','green','black','blue', 'darkorange' ) 85 if(is.na(colors)
86 {
87 color <- c('red', 'yellow', 'green', 'black', 'blue', 'darkorange')
88 }
89 else
90 {
91 # Handle specified colors for components.
92 cStr <- unlist(colors)
93 color <- c()
94 items <- strsplit(cStr, ",")
95 for (item in items) {
96 color <- c(color, item)
97 }
98 num_colors_specified = length(color)
99 if num_colors_specified < num_components:
100 {
101 for (i in num_colors_specified:num_components)
102 {
103 if !(any(color=='red')
104 {
105 color <- c(color, 'red')
106 }
107 else if !(any(color=='yellow')
108 {
109 color <- c(color, 'yellow')
110 }
111 else if !(any(color=='green')
112 {
113 color <- c(color, 'green')
114 }
115 else if !(any(color=='black')
116 {
117 color <- c(color, 'black')
118 }
119 else if !(any(color=='blue')
120 {
121 color <- c(color, 'blue')
122 }
123 else
124 {
125 color <- c(color, 'darkorange')
126 }
127 }
128 }
129 }
93 for (i in 1:length(mu)) 130 for (i in 1:length(mu))
94 { 131 {
95 lines(vx, g[,i] * h, lwd=2, col=color[i]) 132 lines(vx, g[,i] * h, lwd=2, col=color[i])
96 } 133 }
97 } 134 }
115 return(fx) 152 return(fx)
116 } 153 }
117 154
118 # Read in the components data. 155 # Read in the components data.
119 components_data <- read.delim(opt$components_input, header=TRUE) 156 components_data <- read.delim(opt$components_input, header=TRUE)
120 # Get the number of components. 157 num_components <- opt$num_comp
121 num_components <- get_num_components(components_data)
122 158
123 # Set pi, mu, var. 159 # Set pi, mu, var.
124 items <- get_pi_mu_var(components_data, num_components) 160 items <- get_pi_mu_var(components_data, num_components)
125 if (num_components == 1) 161 if (num_components == 1)
126 { 162 {
158 mu <- items[7:12] 194 mu <- items[7:12]
159 var <- items[13:18] 195 var <- items[13:18]
160 } 196 }
161 197
162 # Plot the output. 198 # Plot the output.
163 plot_ks(opt$kaks_input, opt$output, pi, mu, var) 199 plot_ks(opt$kaks_input, num_components, opt$colors, opt$output, pi, mu, var)