Mercurial > repos > greg > kaks_analysis_barplot

comparison kaks_analysis_barplot.R @ 3:9ec2b94ff094 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Wed, 08 Mar 2017 13:57:03 -0500
parents b4f599423810
children
comparison
equal deleted inserted replaced
2:b4f599423810 3:9ec2b94ff094
1 #!/usr/bin/env Rscript 1 #!/usr/bin/env Rscript
2 2
3 suppressPackageStartupMessages(library("optparse")) 3 suppressPackageStartupMessages(library("optparse"))
4 4
5 option_list <- list( 5 option_list <- list(
6 make_option(c("-c", "--components", action="store", dest="components", 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("-o", "--output", action="store", dest="output", default=NULL, help="Output 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 ) 9 )
9 10
10 parser <- OptionParser(usage="%prog [options] file", option_list=option_list) 11 parser <- OptionParser(usage="%prog [options] file", option_list=option_list)
11 args <- parse_args(parser, positional_arguments=TRUE) 12 args <- parse_args(parser, positional_arguments=TRUE)
12 opt <- args$options 13 opt <- args$options
13 14
14 15
15 get_num_components=function(components_dataset) 16 get_num_components = function(components_data)
16 { 17 {
17 # Read in the components data.
18 components_data <- read.delim(components_dataset, header=TRUE);
19 # Get the max of the number_comp column. 18 # Get the max of the number_comp column.
20 num_components <- max(components_data[3, ], na.rm=TRUE); 19 number_comp = components_data[, 3]
21 return = c(components_data, num_components) 20 num_components <- max(number_comp, na.rm=TRUE);
22 return 21 num_components
23 } 22 }
24 23
25 get_pi_mu_var = function(components_data, num_components) { 24 get_pi_mu_var = function(components_data, num_components)
25 {
26 # FixMe: enhance this to generically handle any integer value for num_components. 26 # FixMe: enhance this to generically handle any integer value for num_components.
27 if (num_components == 1) { 27 if (num_components == 1)
28 {
28 pi <- c(components_data[1, 9]); 29 pi <- c(components_data[1, 9]);
29 mu <- c(components_data[1, 7]); 30 mu <- c(components_data[1, 7]);
30 var <- c(components_data[1, 8]); 31 var <- c(components_data[1, 8]);
31 } else if (num_components == 2) { 32 }
33 else if (num_components == 2)
34 {
32 pi <- c(components_data[2, 9], components_data[3, 9]); 35 pi <- c(components_data[2, 9], components_data[3, 9]);
33 mu <- c(components_data[2, 7], components_data[3, 7]); 36 mu <- c(components_data[2, 7], components_data[3, 7]);
34 var <- c(components_data[2, 8], components_data[3, 8]); 37 var <- c(components_data[2, 8], components_data[3, 8]);
35 } else if (num_components == 3) { 38 }
39 else if (num_components == 3)
40 {
36 pi <- c(components_data[4, 9], components_data[5, 9], components_data[6, 9]); 41 pi <- c(components_data[4, 9], components_data[5, 9], components_data[6, 9]);
37 mu <- c(components_data[4, 7], components_data[5, 7], components_data[6, 7]); 42 mu <- c(components_data[4, 7], components_data[5, 7], components_data[6, 7]);
38 var <- c(components_data[4, 8], components_data[5, 8], components_data[6, 8]); 43 var <- c(components_data[4, 8], components_data[5, 8], components_data[6, 8]);
39 } else if (num_components == 4) { 44 }
45 else if (num_components == 4)
46 {
40 pi <- c(components_data[7, 9], components_data[8, 9], components_data[9, 9], components_data[10, 9]); 47 pi <- c(components_data[7, 9], components_data[8, 9], components_data[9, 9], components_data[10, 9]);
41 mu <- c(components_data[7, 7], components_data[8, 7], components_data[9, 7], components_data[10, 7]); 48 mu <- c(components_data[7, 7], components_data[8, 7], components_data[9, 7], components_data[10, 7]);
42 var <- c(components_data[7, 8], components_data[8, 8], components_data[9, 8], components_data[10, 8]); 49 var <- c(components_data[7, 8], components_data[8, 8], components_data[9, 8], components_data[10, 8]);
50 }
43 return = c(pi, mu, var) 51 return = c(pi, mu, var)
44 return 52 return
45 } 53 }
46 54
47 plot_ks<-function(ksfile, pi, mu, var) { 55 plot_ks<-function(kaks_input, output, pi, mu, var)
48 #change bin width 56 {
49 bin <- 0.05 * seq(0, 40); 57 # Start PDF device driver to save charts to output.
50 kaks <- read.table(file=ksfile, header=T); 58 pdf(file=output, bg="white")
51 kaks <- kaks[kaks$Ks<2,]; 59 # Change bin width
52 h.kst <- hist(kaks$Ks, breaks=bin, plot=F); 60 bin <- 0.05 * seq(0, 40);
53 nc <- h.kst$counts; 61 kaks <- read.table(file=kaks_input, header=T);
54 vx <- h.kst$mids; 62 kaks <- kaks[kaks$Ks<2,];
55 ntot <- sum(nc); 63 h.kst <- hist(kaks$Ks, breaks=bin, plot=F);
56 # Set margin for plot bottom, left top, right. 64 nc <- h.kst$counts;
57 par(mai=c(0.5, 0.5, 0, 0)); 65 vx <- h.kst$mids;
58 # Plot dimension in inches. 66 ntot <- sum(nc);
59 par(pin=c(2.5, 2.5)); 67 # Set margin for plot bottom, left top, right.
60 g <- calculate_fitted_density(pi, mu, var); 68 par(mai=c(0.5, 0.5, 0, 0));
61 h <- ntot * 2.5 / sum(g); 69 # Plot dimension in inches.
62 vx <- seq(1, 100) * 0.02; 70 par(pin=c(2.5, 2.5));
63 ymax <- max(nc) + 5; 71 g <- calculate_fitted_density(pi, mu, var);
64 barplot(nc, space=0.25, offset=0, width=0.04, xlim=c(0,2), ylim=c(0, ymax)); 72 h <- ntot * 2.5 / sum(g);
65 # Add x-axis. 73 vx <- seq(1, 100) * 0.02;
66 axis(1); 74 ymax <- max(nc) + 5;
67 color <- c('green', 'blue', 'black', 'red'); 75 barplot(nc, space=0.25, offset=0, width=0.04, xlim=c(0,2), ylim=c(0, ymax));
68 for (i in 1:length(mu)) { 76 # Add x-axis.
69 lines(vx, g[,i] * h, lwd=2, col=color[i]); 77 axis(1);
70 } 78 color <- c('green', 'blue', 'black', 'red');
79 for (i in 1:length(mu))
80 {
81 lines(vx, g[,i] * h, lwd=2, col=color[i]);
82 }
71 }; 83 };
72 84
73 calculate_fitted_density <- function(pi, mu, var) { 85 calculate_fitted_density <- function(pi, mu, var)
74 comp <- length(pi); 86 {
75 var <- var/mu^2; 87 comp <- length(pi);
76 mu <- log(mu);; 88 var <- var/mu^2;
77 #calculate lognormal density 89 mu <- log(mu);
78 vx <- seq(1, 100) * 0.02; 90 #calculate lognormal density
79 fx <- matrix(0, 100, comp); 91 vx <- seq(1, 100) * 0.02;
80 for (i in 1:100) { 92 fx <- matrix(0, 100, comp);
81 for (j in 1:comp) { 93 for (i in 1:100)
82 fx[i, j] <- pi[j] * dlnorm(vx[i], meanlog=mu[j], sdlog=(sqrt(var[j]))); 94 {
83 }; 95 for (j in 1:comp)
84 }; 96 {
85 fx; 97 fx[i, j] <- pi[j] * dlnorm(vx[i], meanlog=mu[j], sdlog=(sqrt(var[j])));
98 };
99 };
100 fx;
86 } 101 }
87 102
88 # Read in the components data and get the number of components. 103 # Read in the components data.
89 items <- get_num_components(opt$components) 104 components_data <- read.delim(opt$components_input, header=TRUE);
90 components_data <- items[1] 105 # Get the number of components.
91 num_components <- items[2] 106 num_components <- get_num_components(components_data)
92
93 # Set output file name.
94 if (is.null(opt$output)) {
95 # Name the output file based on the name of the
96 # input file, properly handling full path if passed.
97 input_filename = basename(opt$components)
98 items = strsplit(input_filename, ".")
99 output_filename <- paste(items[1], ".components.", num_components, ".pdf")
100 } else {
101 output_filename <- opt$output
102 }
103 107
104 # Set pi, mu, var. 108 # Set pi, mu, var.
105 items <- get_pi_mu_var(components_data, num_components) 109 items <- get_pi_mu_var(components_data, num_components);
106 pi <- items[1] 110 pi <- items[1];
107 mu <- items[2] 111 mu <- items[2];
108 var <- items[3] 112 var <- items[3];
109 113
110 # Plot the output. 114 # Plot the output.
111 plot_ks(output_filename, pi, mu, var) 115 plot_ks(opt$kaks_input, opt$output, pi, mu, var);