# HG changeset patch
# User greg
# Date 1489689699 14400
# Node ID 8d18cb8396a7b076eeb887d850270b178ff7394d
# Parent a419970e9c195713174e1aee365dfe969cd3d5da
Uploaded
diff -r a419970e9c19 -r 8d18cb8396a7 kaks_analysis_barplot.R
--- a/kaks_analysis_barplot.R Wed Mar 08 14:07:30 2017 -0500
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,115 +0,0 @@
-#!/usr/bin/env Rscript
-
-suppressPackageStartupMessages(library("optparse"))
-
-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")
-)
-
-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);
- num_components
-}
-
-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)
- {
- pi <- c(components_data[1, 9]);
- mu <- c(components_data[1, 7]);
- var <- c(components_data[1, 8]);
- }
- else if (num_components == 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)
- {
- 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)
- {
- 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]);
- }
- return = c(pi, mu, var)
- return
-}
-
-plot_ks<-function(kaks_input, output, pi, mu, var)
-{
- # Start PDF device driver to save charts to output.
- pdf(file=output, bg="white")
- # Change bin width
- bin <- 0.05 * seq(0, 40);
- kaks <- read.table(file=kaks_input, header=T);
- kaks <- kaks[kaks$Ks<2,];
- h.kst <- hist(kaks$Ks, breaks=bin, plot=F);
- nc <- h.kst$counts;
- vx <- h.kst$mids;
- ntot <- sum(nc);
- # Set margin for plot bottom, left top, right.
- par(mai=c(0.5, 0.5, 0, 0));
- # Plot dimension in inches.
- par(pin=c(2.5, 2.5));
- g <- calculate_fitted_density(pi, mu, var);
- h <- ntot * 2.5 / sum(g);
- vx <- seq(1, 100) * 0.02;
- ymax <- max(nc) + 5;
- barplot(nc, space=0.25, offset=0, width=0.04, xlim=c(0,2), ylim=c(0, ymax));
- # Add x-axis.
- axis(1);
- color <- c('green', 'blue', 'black', 'red');
- for (i in 1:length(mu))
- {
- lines(vx, g[,i] * h, lwd=2, col=color[i]);
- }
-};
-
-calculate_fitted_density <- function(pi, mu, var)
-{
- comp <- length(pi);
- var <- var/mu^2;
- mu <- log(mu);
- #calculate lognormal density
- vx <- seq(1, 100) * 0.02;
- 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])));
- };
- };
- fx;
-}
-
-# 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)
-
-# Set pi, mu, var.
-items <- get_pi_mu_var(components_data, num_components);
-pi <- items[1];
-mu <- items[2];
-var <- items[3];
-
-# Plot the output.
-plot_ks(opt$kaks_input, opt$output, pi, mu, var);
diff -r a419970e9c19 -r 8d18cb8396a7 kaks_analysis_barplot.xml
--- a/kaks_analysis_barplot.xml Wed Mar 08 14:07:30 2017 -0500
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,43 +0,0 @@
-
-
-
- r-optparse
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-**What it does**
-
-Draws a barplot of the significant components in the ks distribution.
-
-
-
- @unpublished{None,
- author = {Eric Wafula},
- title = {None},
- year = {None},
- url = {https://github.com/dePamphilis/PlantTribes}
- }
-
-
-
diff -r a419970e9c19 -r 8d18cb8396a7 kaks_distribution.R
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/kaks_distribution.R Thu Mar 16 14:41:39 2017 -0400
@@ -0,0 +1,115 @@
+#!/usr/bin/env Rscript
+
+suppressPackageStartupMessages(library("optparse"))
+
+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")
+)
+
+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);
+ num_components
+}
+
+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)
+ {
+ pi <- c(components_data[1, 9]);
+ mu <- c(components_data[1, 7]);
+ var <- c(components_data[1, 8]);
+ }
+ else if (num_components == 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)
+ {
+ 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)
+ {
+ 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]);
+ }
+ return = c(pi, mu, var)
+ return
+}
+
+plot_ks<-function(kaks_input, output, pi, mu, var)
+{
+ # Start PDF device driver to save charts to output.
+ pdf(file=output, bg="white")
+ # Change bin width
+ bin <- 0.05 * seq(0, 40);
+ kaks <- read.table(file=kaks_input, header=T);
+ kaks <- kaks[kaks$Ks<2,];
+ h.kst <- hist(kaks$Ks, breaks=bin, plot=F);
+ nc <- h.kst$counts;
+ vx <- h.kst$mids;
+ ntot <- sum(nc);
+ # Set margin for plot bottom, left top, right.
+ par(mai=c(0.5, 0.5, 0, 0));
+ # Plot dimension in inches.
+ par(pin=c(2.5, 2.5));
+ g <- calculate_fitted_density(pi, mu, var);
+ h <- ntot * 2.5 / sum(g);
+ vx <- seq(1, 100) * 0.02;
+ ymax <- max(nc) + 5;
+ barplot(nc, space=0.25, offset=0, width=0.04, xlim=c(0,2), ylim=c(0, ymax));
+ # Add x-axis.
+ axis(1);
+ color <- c('green', 'blue', 'black', 'red');
+ for (i in 1:length(mu))
+ {
+ lines(vx, g[,i] * h, lwd=2, col=color[i]);
+ }
+};
+
+calculate_fitted_density <- function(pi, mu, var)
+{
+ comp <- length(pi);
+ var <- var/mu^2;
+ mu <- log(mu);
+ #calculate lognormal density
+ vx <- seq(1, 100) * 0.02;
+ 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])));
+ };
+ };
+ fx;
+}
+
+# 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)
+
+# Set pi, mu, var.
+items <- get_pi_mu_var(components_data, num_components);
+pi <- items[1];
+mu <- items[2];
+var <- items[3];
+
+# Plot the output.
+plot_ks(opt$kaks_input, opt$output, pi, mu, var);
diff -r a419970e9c19 -r 8d18cb8396a7 kaks_distribution.xml
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/kaks_distribution.xml Thu Mar 16 14:41:39 2017 -0400
@@ -0,0 +1,43 @@
+
+ plots synonymous (ks) distribution and fit significant component
+
+ r-optparse
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+**What it does**
+
+Draws a barplot of the significant components in the ks distribution.
+
+
+
+ @article{None,
+ journal = {None},
+ author = {1. Wafula EK},
+ title = {Manuscript in preparation},
+ year = {None},
+ url = {https://github.com/dePamphilis/PlantTribes},}
+
+
+