Mercurial > repos > greg > ks_distribution
changeset 44:69ac2edffdde draft
Uploaded
author | greg |
---|---|
date | Mon, 26 Jun 2017 13:58:37 -0400 |
parents | 5074cb26597a |
children | f8b303b9f455 |
files | ks_distribution.R |
diffstat | 1 files changed, 28 insertions(+), 71 deletions(-) [+] |
line wrap: on
line diff
--- a/ks_distribution.R Mon Jun 26 12:13:52 2017 -0400 +++ b/ks_distribution.R Mon Jun 26 13:58:37 2017 -0400 @@ -7,7 +7,7 @@ make_option(c("-k", "--kaks_input"), action="store", dest="kaks_input", help="KaKs analysis input dataset"), make_option(c("-n", "--number_comp"), action="store", dest="number_comp", type="integer", help="Number of significant components in the Ks distribution"), make_option(c("-o", "--output"), action="store", dest="output", help="Output dataset"), - make_option(c("-r", "--colors"), action="store", default="", help="List of component colors") + make_option(c("-r", "--colors"), action="store", default=NA, help="List of component colors") ) parser <- OptionParser(usage="%prog [options] file", option_list=option_list) @@ -17,47 +17,30 @@ set_component_colors = function(colors, number_comp) { # Handle colors for components. - if (nchar(colors) == 0) - { + if (is.na(colors) == 0) { # Randomly specify colors for components. specified_colors <- c('red', 'yellow', 'green', 'black', 'blue', 'darkorange') - } - else - { + } else { # Handle selected colors for components. specified_colors <- c() items <- strsplit(colors, ",") - for (item in items) - { + for (item in items) { specified_colors <- c(specified_colors, item) } num_colors_specified <- length(specified_colors) - if (num_colors_specified < number_comp) - { - for (i in num_colors_specified:number_comp) - { - if (!any(specified_colors == 'red')) - { + if (num_colors_specified < number_comp) { + for (i in num_colors_specified:number_comp) { + if (!any(specified_colors == 'red')) { specified_colors <- c(specified_colors, 'red') - } - else if (!any(specified_colors == 'yellow')) - { + } else if (!any(specified_colors == 'yellow')) { specified_colors <- c(specified_colors, 'yellow') - } - else if (!any(specified_colors == 'green')) - { + } else if (!any(specified_colors == 'green')) { specified_colors <- c(specified_colors, 'green') - } - else if (!any(specified_colors == 'black')) - { + } else if (!any(specified_colors == 'black')) { specified_colors <- c(specified_colors, 'black') - } - else if (!any(specified_colors == 'blue')) - { + } else if (!any(specified_colors == 'blue')) { specified_colors <- c(specified_colors, 'blue') - } - else - { + } else { specified_colors <- c(specified_colors, 'darkorange') } } @@ -69,38 +52,27 @@ get_pi_mu_var = function(components_data, number_comp) { # FixMe: enhance this to generically handle any integer value for number_comp. - if (number_comp == 1) - { + if (number_comp == 1) { pi <- c(components_data[1, 9]) mu <- c(components_data[1, 7]) var <- c(components_data[1, 8]) - } - else if (number_comp == 2) - { + } else if (number_comp == 2) { pi <- c(components_data[2, 9], components_data[3, 9]) mu <- c(components_data[2, 7], components_data[3, 7]) var <- c(components_data[2, 8], components_data[3, 8]) - } - else if (number_comp == 3) - { + } else if (number_comp == 3) { pi <- c(components_data[4, 9], components_data[5, 9], components_data[6, 9]) mu <- c(components_data[4, 7], components_data[5, 7], components_data[6, 7]) var <- c(components_data[4, 8], components_data[5, 8], components_data[6, 8]) - } - else if (number_comp == 4) - { + } else if (number_comp == 4) { pi <- c(components_data[7, 9], components_data[8, 9], components_data[9, 9], components_data[10, 9]) mu <- c(components_data[7, 7], components_data[8, 7], components_data[9, 7], components_data[10, 7]) var <- c(components_data[7, 8], components_data[8, 8], components_data[9, 8], components_data[10, 8]) - } - else if (number_comp == 5) - { + } else if (number_comp == 5) { pi <- c(components_data[11, 9], components_data[12, 9], components_data[13, 9], components_data[14, 9], components_data[15, 9]) mu <- c(components_data[11, 7], components_data[12, 7], components_data[13, 7], components_data[14, 7], components_data[15, 7]) var <- c(components_data[11, 8], components_data[12, 8], components_data[13, 8], components_data[14, 8], components_data[15, 8]) - } - else if (number_comp == 6) - { + } else if (number_comp == 6) { 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]) 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]) 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]) @@ -134,8 +106,7 @@ barplot(nc, space=0.25, offset=0, width=0.04, xlim=c(0, max_ks), ylim=c(0, ymax), col="lightpink1", border="lightpink3") # Add x-axis. axis(1) - for (i in 1:length(mu)) - { + for (i in 1:length(mu)) { lines(vx, g[,i] * h, lwd=2, col=specified_colors[i]) } } @@ -148,13 +119,10 @@ # Calculate lognormal density. vx <- seq(1, 100) * (max_ks / 100) fx <- matrix(0, 100, comp) - for (i in 1:100) - { - for (j in 1:comp) - { + for (i in 1:100) { + for (j in 1:comp) { fx[i, j] <- pi[j] * dlnorm(vx[i], meanlog=mu[j], sdlog=(sqrt(var[j]))) - if (is.nan(fx[i,j])) - { + if (is.nan(fx[i,j])) { fx[i,j]<-0 } } @@ -171,38 +139,27 @@ # Set pi, mu, var. items <- get_pi_mu_var(components_data, number_comp) -if (number_comp == 1) -{ +if (number_comp == 1) { pi <- items[1] mu <- items[2] var <- items[3] -} -else if (number_comp == 2) -{ +} else if (number_comp == 2) { pi <- items[1:2] mu <- items[3:4] var <- items[5:6] -} -else if (number_comp == 3) -{ +} else if (number_comp == 3) { pi <- items[1:3] mu <- items[4:6] var <- items[7:9] -} -else if (number_comp == 4) -{ +} else if (number_comp == 4) { pi <- items[1:4] mu <- items[5:8] var <- items[9:12] -} -else if (number_comp == 5) -{ +} else if (number_comp == 5) { pi <- items[1:5] mu <- items[6:10] var <- items[11:15] -} -else if (number_comp == 6) -{ +} else if (number_comp == 6) { pi <- items[1:6] mu <- items[7:12] var <- items[13:18]