Mercurial > repos > greg > plant_tribes_ks_distribution

changeset 1:56f42cc1dd58 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Wed, 28 Jun 2017 11:20:17 -0400
parents c5846258c458
children 2c8e564adede
files components.tabular kaks_input1.tabular ks_distribution.R ks_distribution.xml macros.xml output.pdf test-data/components.ptkscmp test-data/kaks_input1.tabular test-data/output.pdf
diffstat 9 files changed, 186 insertions(+), 144 deletions(-) [+]
line wrap: on
line diff
--- a/components.tabular	Thu Jun 08 12:55:49 2017 -0400
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,2 +0,0 @@
-species	n	number_comp	lnL	AIC	BIC	mean	variance	porportion
-species1	3	1	-3.4750	6.95	6.95	3.1183	5.7732	1.00
--- a/kaks_input1.tabular	Thu Jun 08 12:55:49 2017 -0400
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,4 +0,0 @@
-SEQ1	SEQ2	Ka	Ks	Ka\Ks
-contig_241;	contig_241	1.5312	7.1619	0.2138
-contig_300	contig_300;	0.8653	3.7872	0.2285
-contig_586	contig_586;	1.7791	1.1181	1.5912
--- a/ks_distribution.R	Thu Jun 08 12:55:49 2017 -0400
+++ b/ks_distribution.R	Wed Jun 28 11:20:17 2017 -0400
@@ -5,67 +5,86 @@
 option_list <- list(
     make_option(c("-c", "--components_input"), action="store", dest="components_input", help="Ks significant components input dataset"),
     make_option(c("-k", "--kaks_input"), action="store", dest="kaks_input", help="KaKs analysis input dataset"),
-    make_option(c("-o", "--output"), action="store", dest="output", help="Output 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("-s", "--specified_colors"), action="store", dest="specified_colors", default=NULL, help="List of component colors")
 )
 
 parser <- OptionParser(usage="%prog [options] file", option_list=option_list)
 args <- parse_args(parser, positional_arguments=TRUE)
 opt <- args$options
 
-
-get_num_components = function(components_data)
-{
-    # Get the max of the number_comp column.
-    number_comp = components_data[, 3]
-    num_components <- max(number_comp, na.rm=TRUE)
-    return(num_components)
+set_component_colors = function(colors, number_comp) {
+    # Handle colors for components.
+    if (is.null(colors)) {
+        # Randomly specify colors for components.
+        component_colors <- c("red", "yellow", "green", "black", "blue", "darkorange")
+    } else {
+        # Handle selected colors for components.
+        component_colors <- c()
+        colors <- as.character(colors) 
+        items <- strsplit(colors, ",") 
+        for (item in items) { 
+            component_colors <- c(component_colors, item)
+        }
+        num_colors_specified <- length(component_colors)
+        if (num_colors_specified < number_comp) {
+            # The number of selected colors is less than the number of
+            # components, so we'll add random colors that were not
+            # selected to the set of component colors until we have a
+            # color for each component.
+            loop_count <- number_comp - num_colors_specified
+            for (i in 1:loop_count) {
+                if (!(is.element("red", component_colors))) {
+                    component_colors <- c(component_colors, "red")
+                } else if (!(is.element("yellow", component_colors))) {
+                    component_colors <- c(component_colors, "yellow")
+                } else if (!(is.element("green", component_colors))) {
+                    component_colors <- c(component_colors, "green")
+                } else if (!(is.element("black", component_colors))) {
+                    component_colors <- c(component_colors, "black")
+                } else if (!(is.element("blue", component_colors))) {
+                    component_colors <- c(component_colors, "blue")
+                } else if (!(is.element("darkorange", component_colors))) {
+                    component_colors <- c(component_colors, "darkorange")
+                }
+            }
+        }
+    }
+    return(component_colors)
 }
 
-get_pi_mu_var = function(components_data, num_components)
-{
-    # FixMe: enhance this to generically handle any integer value for num_components.
-    if (num_components == 1)
-    {
+get_pi_mu_var = function(components_data, number_comp) {
+    if (number_comp == 1) {
         pi <- c(components_data[1, 9])
         mu <- c(components_data[1, 7])
         var <- c(components_data[1, 8])
-    }
-    else if (num_components == 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 (num_components == 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 (num_components == 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 (num_components == 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 (num_components == 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])
     }
-    results = c(pi, mu, var)
+    results <- c(pi, mu, var)
     return(results)
 }
 
-plot_ks<-function(kaks_input, output, pi, mu, var)
-{
+plot_ks<-function(kaks_input, component_colors, output, pi, mu, var) {
     # Start PDF device driver to save charts to output.
     pdf(file=output, bg="white")
     kaks <- read.table(file=kaks_input, header=T)
@@ -89,27 +108,24 @@
     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)
-    color <- c('red', 'yellow','green','black','blue', 'darkorange' )
-    for (i in 1:length(mu))
-    {
-       lines(vx, g[,i] * h, lwd=2, col=color[i])
+    for (i in 1:length(mu)) {
+       lines(vx, g[,i] * h, lwd=2, col=component_colors[i])
     }
 }
 
-calculate_fitted_density <- function(pi, mu, var, max_ks)
-{
+calculate_fitted_density <- function(pi, mu, var, max_ks) {
     comp <- length(pi)
     var <- var/mu^2
     mu <- log(mu)
     # 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)
-        {
-           fx[i, j] <- pi[j] * dlnorm(vx[i], meanlog=mu[j], sdlog=(sqrt(var[j])))
-           if (is.nan(fx[i,j])) fx[i,j]<-0
+    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])) {
+                fx[i,j]<-0
+            }
         }
      }
     return(fx)
@@ -117,47 +133,42 @@
 
 # Read in the components data.
 components_data <- read.delim(opt$components_input, header=TRUE)
-# Get the number of components.
-num_components <- get_num_components(components_data)
+number_comp <- as.integer(opt$number_comp)
+if (number_comp == 0) {
+    # Default to 1 component to allow functional testing.
+    number_comp <- 1
+}
+
+# Set component colors.
+component_colors <- set_component_colors(opt$specified_colors, number_comp)
 
 # Set pi, mu, var.
-items <- get_pi_mu_var(components_data, num_components)
-if (num_components == 1)
-{
-	pi <- items[1]
-	mu <- items[2]
-	var <- items[3]
-}
-if (num_components == 2)
-{
-	pi <- items[1:2]
-	mu <- items[3:4]
-	var <- items[5:6]
-}
-if (num_components == 3)
-{
-	pi <- items[1:3]
-	mu <- items[4:6]
-	var <- items[7:9]
-}
-if (num_components == 4)
-{
-	pi <- items[1:4]
-	mu <- items[5:8]
-	var <- items[9:12]
-}
-if (num_components == 5)
-{
-	pi <- items[1:5]
-	mu <- items[6:10]
-	var <- items[11:15]
-}
-if (num_components == 6)
-{
-	pi <- items[1:6]
-	mu <- items[7:12]
-	var <- items[13:18]
+items <- get_pi_mu_var(components_data, number_comp)
+if (number_comp == 1) {
+    pi <- items[1]
+    mu <- items[2]
+    var <- items[3]
+} else if (number_comp == 2) {
+    pi <- items[1:2]
+    mu <- items[3:4]
+    var <- items[5:6]
+} else if (number_comp == 3) {
+    pi <- items[1:3]
+    mu <- items[4:6]
+    var <- items[7:9]
+} else if (number_comp == 4) {
+    pi <- items[1:4]
+    mu <- items[5:8]
+    var <- items[9:12]
+} else if (number_comp == 5) {
+    pi <- items[1:5]
+    mu <- items[6:10]
+    var <- items[11:15]
+} else if (number_comp == 6) {
+    pi <- items[1:6]
+    mu <- items[7:12]
+    var <- items[13:18]
 }
 
 # Plot the output.
-plot_ks(opt$kaks_input, opt$output, pi, mu, var)
+plot_ks(opt$kaks_input, component_colors, opt$output, pi, mu, var)
--- a/ks_distribution.xml	Thu Jun 08 12:55:49 2017 -0400
+++ b/ks_distribution.xml	Wed Jun 28 11:20:17 2017 -0400
@@ -1,4 +1,4 @@
-<tool id="ks_distribution" name="KsDistribution" version="1.0.0">
+<tool id="ks_distribution" name="KsDistribution" version="1.0.1">
     <description>plots the distribution of synonymous substitution (Ks) rates and fits significant component(s)</description>
     <macros>
         <import>macros.xml</import>
@@ -8,11 +8,39 @@
 Rscript $__tool_directory__/ks_distribution.R
 -k '$input'
 -c '$components'
+-n $components.metadata.number_comp
 -o '$output'
+#if str($choose_colors_cond.choose_colors) == 'yes':
+    #set specified_colors = list()
+    #for $selection in $choose_colors_cond.colors:
+        $specified_colors.append(str($selection.color))
+    #end for
+    #set specified_colors = ",".join($specified_colors)
+    -s '$specified_colors'
+#end if
     ]]></command>
     <inputs>
         <param name="input" format="tabular" type="data" label="KaKsAnalysis tabular file" />
-        <param name="components" format="tabular" type="data" label="Significant components" />
+        <param name="components" format="ptkscmp" type="data" label="Significant components" />
+        <conditional name="choose_colors_cond">
+            <param name="choose_colors" type="select" label="Choose colors for significant components">
+                <option value="no" selected="true">No</option>
+                <option value="yes">Yes</option>
+            </param>
+            <when value="no" />
+            <when value="yes">
+                <repeat name="colors" title="Component colors" min="1">
+                    <param name="color" type="select" label="Color">
+                        <option value="red" selected="True">red</option>
+                        <option value="green">green</option>
+                        <option value="blue">blue</option>
+                        <option value="yellow">yellow</option>
+                        <option value="darkorange">dark orange</option>
+                        <option value="black">black</option>
+                    </param>
+                </repeat>
+            </when>
+        </conditional>
     </inputs>
     <outputs>
         <data name="output" format="pdf"/>
@@ -20,7 +48,7 @@
     <tests>
         <test>
             <param name="input" value="kaks_input1.tabular" ftype="tabular" />
-            <param name="components" value="components.tabular" ftype="tabular" />
+            <param name="components" value="components.ptkscmp" ftype="ptkscmp" />
             <output name="output" file="output.pdf" ftype="pdf" compare="contains" />
         </test>
     </tests>
@@ -37,6 +65,9 @@
 
  * **Synonymous substitution rates** - estimated synonymous substitution (Ks) rates output file produced by the KaKsAnalysis tool selected from your history.
  * **Synonymous components** - estimated significant component(s) output file produced by the KaKsAnalysis tool selected from your history.
+ * **Choose colors for significant components** - select 'Yes' to specify component colors or 'No' for colors chosen randlomly.
+
+   * **Component colors** - select a color from the palette for each component (colors will be chosen randomly for unspecified components).
     </help>
     <citations>
         <expand macro="citation1" />
--- a/macros.xml	Thu Jun 08 12:55:49 2017 -0400
+++ b/macros.xml	Wed Jun 28 11:20:17 2017 -0400
@@ -3,7 +3,7 @@
     <token name="@WRAPPER_VERSION@">1.0</token>
     <xml name="requirements_assembly_post_processor">
         <requirements>
-            <requirement type="package" version="1.0.0">plant_tribes_assembly_post_processor</requirement>
+            <requirement type="package" version="1.0.1">plant_tribes_assembly_post_processor</requirement>
         </requirements>
     </xml>
     <xml name="requirements_gene_family_aligner">
--- a/output.pdf	Thu Jun 08 12:55:49 2017 -0400
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,52 +0,0 @@
-1 0 obj
-<<
-/Title (R Graphics Output)
-/Creator (R)
->>
-endobj
-2 0 obj
-<< /Type /Catalog /Pages 3 0 R >>
-endobj
-7 0 obj
-<< /Type /Page /Parent 3 0 R /Contents 8 0 R /Resources 4 0 R >>
-endobj
-8 0 obj
-<<
->>
-stream
-endobj
-3 0 obj
-<< /Type /Pages /Kids [ 7 0 R ] /Count 1 /MediaBox [0 0 504 504] >>
-endobj
-4 0 obj
-<<
-/ProcSet [/PDF /Text]
-/Font <</F2 10 0 R >>
-/ExtGState << >>
-/ColorSpace << /sRGB 5 0 R >>
->>
-endobj
-5 0 obj
-[/ICCBased 6 0 R]
-endobj
-6 0 obj
-<< /Alternate /DeviceRGB /N 3 /Length 2596 /Filter /FlateDecode >>
-stream
-9 0 obj
-<<
-/Type /Encoding /BaseEncoding /WinAnsiEncoding
-/Differences [ 45/minus 96/quoteleft
-144/dotlessi /grave /acute /circumflex /tilde /macron /breve /dotaccent
-/dieresis /.notdef /ring /cedilla /.notdef /hungarumlaut /ogonek /caron /space]
->>
-endobj
-10 0 obj
-<< /Type /Font /Subtype /Type1 /Name /F2 /BaseFont /Helvetica
-/Encoding 9 0 R >>
-endobj
-xref
-0 11
-trailer
-<< /Size 11 /Info 1 0 R /Root 2 0 R >>
-startxref
-%%EOF
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/components.ptkscmp	Wed Jun 28 11:20:17 2017 -0400
@@ -0,0 +1,2 @@
+species	n	number_comp	lnL	AIC	BIC	mean	variance	porportion
+species1	3	1	-3.4750	6.95	6.95	3.1183	5.7732	1.00
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/kaks_input1.tabular	Wed Jun 28 11:20:17 2017 -0400
@@ -0,0 +1,4 @@
+SEQ1	SEQ2	Ka	Ks	Ka\Ks
+contig_241;	contig_241	1.5312	7.1619	0.2138
+contig_300	contig_300;	0.8653	3.7872	0.2285
+contig_586	contig_586;	1.7791	1.1181	1.5912
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/output.pdf	Wed Jun 28 11:20:17 2017 -0400
@@ -0,0 +1,52 @@
+1 0 obj
+<<
+/Title (R Graphics Output)
+/Creator (R)
+>>
+endobj
+2 0 obj
+<< /Type /Catalog /Pages 3 0 R >>
+endobj
+7 0 obj
+<< /Type /Page /Parent 3 0 R /Contents 8 0 R /Resources 4 0 R >>
+endobj
+8 0 obj
+<<
+>>
+stream
+endobj
+3 0 obj
+<< /Type /Pages /Kids [ 7 0 R ] /Count 1 /MediaBox [0 0 504 504] >>
+endobj
+4 0 obj
+<<
+/ProcSet [/PDF /Text]
+/Font <</F2 10 0 R >>
+/ExtGState << >>
+/ColorSpace << /sRGB 5 0 R >>
+>>
+endobj
+5 0 obj
+[/ICCBased 6 0 R]
+endobj
+6 0 obj
+<< /Alternate /DeviceRGB /N 3 /Length 2596 /Filter /FlateDecode >>
+stream
+9 0 obj
+<<
+/Type /Encoding /BaseEncoding /WinAnsiEncoding
+/Differences [ 45/minus 96/quoteleft
+144/dotlessi /grave /acute /circumflex /tilde /macron /breve /dotaccent
+/dieresis /.notdef /ring /cedilla /.notdef /hungarumlaut /ogonek /caron /space]
+>>
+endobj
+10 0 obj
+<< /Type /Font /Subtype /Type1 /Name /F2 /BaseFont /Helvetica
+/Encoding 9 0 R >>
+endobj
+xref
+0 11
+trailer
+<< /Size 11 /Info 1 0 R /Root 2 0 R >>
+startxref
+%%EOF