|
0
|
1 #!/usr/bin/env Rscript
|
|
|
2
|
|
|
3 suppressPackageStartupMessages(library("optparse"))
|
|
|
4
|
|
|
5 option_list <- list(
|
|
|
6 make_option(c("-c", "--components"), action="store", dest="components", help="Ks significant components input dataset"),
|
|
|
7 make_option(c("-o", "--output"), action="store", dest="output", default=NULL, help="Output dataset"),
|
|
|
8 )
|
|
|
9
|
|
|
10 parser <- OptionParser(usage="%prog [options] file", option_list=option_list)
|
|
|
11 args <- parse_args(parser, positional_arguments=TRUE)
|
|
|
12 opt <- args$options
|
|
|
13
|
|
|
14
|
|
|
15 get_num_components=function(components_dataset)
|
|
|
16 {
|
|
|
17 # Read in the components data.
|
|
|
18 components_data <- read.delim(components_dataset, header=TRUE);
|
|
|
19 # Get the max of the number_comp column.
|
|
|
20 num_components <- max(components_data[3, ], na.rm=TRUE);
|
|
|
21 return = c(components_data, num_components)
|
|
|
22 return
|
|
|
23 }
|
|
|
24
|
|
|
25 get_pi_mu_var = function(components_data, num_components) {
|
|
|
26 # FixMe: enhance this to generically handle any integer value for num_components.
|
|
|
27 if (num_components == 1) {
|
|
|
28 pi <- c(components_data[1, 9]);
|
|
|
29 mu <- c(components_data[1, 7]);
|
|
|
30 var <- c(components_data[1, 8]);
|
|
|
31 } else if (num_components == 2) {
|
|
|
32 pi <- c(components_data[2, 9], components_data[3, 9]);
|
|
|
33 mu <- c(components_data[2, 7], components_data[3, 7]);
|
|
|
34 var <- c(components_data[2, 8], components_data[3, 8]);
|
|
|
35 } else if (num_components == 3) {
|
|
|
36 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]);
|
|
|
38 var <- c(components_data[4, 8], components_data[5, 8], components_data[6, 8]);
|
|
|
39 } else if (num_components == 4) {
|
|
|
40 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]);
|
|
|
42 var <- c(components_data[7, 8], components_data[8, 8], components_data[9, 8], components_data[10, 8]);
|
|
|
43 return = c(pi, mu, var)
|
|
|
44 return
|
|
|
45 }
|
|
|
46
|
|
|
47 plot_ks<-function(ksfile, pi, mu, var) {
|
|
|
48 #change bin width
|
|
|
49 bin <- 0.05 * seq(0, 40);
|
|
|
50 kaks <- read.table(file=ksfile, header=T);
|
|
|
51 kaks <- kaks[kaks$Ks<2,];
|
|
|
52 h.kst <- hist(kaks$Ks, breaks=bin, plot=F);
|
|
|
53 nc <- h.kst$counts;
|
|
|
54 vx <- h.kst$mids;
|
|
|
55 ntot <- sum(nc);
|
|
|
56 # Set margin for plot bottom, left top, right.
|
|
|
57 par(mai=c(0.5, 0.5, 0, 0));
|
|
|
58 # Plot dimension in inches.
|
|
|
59 par(pin=c(2.5, 2.5));
|
|
|
60 g <- calculate_fitted_density(pi, mu, var);
|
|
|
61 h <- ntot * 2.5 / sum(g);
|
|
|
62 vx <- seq(1, 100) * 0.02;
|
|
|
63 ymax <- max(nc) + 5;
|
|
|
64 barplot(nc, space=0.25, offset=0, width=0.04, xlim=c(0,2), ylim=c(0, ymax));
|
|
|
65 # Add x-axis.
|
|
|
66 axis(1);
|
|
|
67 color <- c('green', 'blue', 'black', 'red');
|
|
|
68 for (i in 1:length(mu)) {
|
|
|
69 lines(vx, g[,i] * h, lwd=2, col=color[i]);
|
|
|
70 }
|
|
|
71 };
|
|
|
72
|
|
|
73 calculate_fitted_density <- function(pi, mu, var) {
|
|
|
74 comp <- length(pi);
|
|
|
75 var <- var/mu^2;
|
|
|
76 mu <- log(mu);;
|
|
|
77 #calculate lognormal density
|
|
|
78 vx <- seq(1, 100) * 0.02;
|
|
|
79 fx <- matrix(0, 100, comp);
|
|
|
80 for (i in 1:100) {
|
|
|
81 for (j in 1:comp) {
|
|
|
82 fx[i, j] <- pi[j] * dlnorm(vx[i], meanlog=mu[j], sdlog=(sqrt(var[j])));
|
|
|
83 };
|
|
|
84 };
|
|
|
85 fx;
|
|
|
86 }
|
|
|
87
|
|
|
88 # Read in the components data and get the number of components.
|
|
|
89 items <- get_num_components(opt$components)
|
|
|
90 components_data <- items[1]
|
|
|
91 num_components <- items[2]
|
|
|
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
|
|
|
104 # Set pi, mu, var.
|
|
|
105 items <- get_pi_mu_var(components_data, num_components)
|
|
|
106 pi <- items[1]
|
|
|
107 mu <- items[2]
|
|
|
108 var <- items[3]
|
|
|
109
|
|
|
110 # Plot the output.
|
|
|
111 plot_ks(output_filename, pi, mu, var)
|
