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]