Mercurial > repos > greg > kaks_analysis_barplot

comparison kaks_distribution.R @ 5:8d18cb8396a7 draft default tip

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