Mercurial > repos > greg > plant_tribes_ks_distribution

comparison ks_distribution.R @ 0:c5846258c458 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Thu, 08 Jun 2017 12:55:49 -0400
parents
children 56f42cc1dd58
comparison
equal deleted inserted replaced
-1:000000000000 0:c5846258c458
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 return(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 else if (num_components == 5)
52 {
53 pi <- c(components_data[11, 9], components_data[12, 9], components_data[13, 9], components_data[14, 9], components_data[15, 9])
54 mu <- c(components_data[11, 7], components_data[12, 7], components_data[13, 7], components_data[14, 7], components_data[15, 7])
55 var <- c(components_data[11, 8], components_data[12, 8], components_data[13, 8], components_data[14, 8], components_data[15, 8])
56 }
57 else if (num_components == 6)
58 {
59 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])
60 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])
61 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])
62 }
63 results = c(pi, mu, var)
64 return(results)
65 }
66
67 plot_ks<-function(kaks_input, output, pi, mu, var)
68 {
69 # Start PDF device driver to save charts to output.
70 pdf(file=output, bg="white")
71 kaks <- read.table(file=kaks_input, header=T)
72 max_ks <- max(kaks$Ks, na.rm=TRUE)
73 # Change bin width
74 max_bin_range <- as.integer(max_ks / 0.05)
75 bin <- 0.05 * seq(0, (max_bin_range + 1 ))
76 kaks <- kaks[kaks$Ks<max_ks,]
77 h.kst <- hist(kaks$Ks, breaks=bin, plot=F)
78 nc <- h.kst$counts
79 vx <- h.kst$mids
80 ntot <- sum(nc)
81 # Set margin for plot bottom, left top, right.
82 par(mai=c(0.5, 0.5, 0, 0))
83 # Plot dimension in inches.
84 par(pin=c(3.0, 3.0))
85 g <- calculate_fitted_density(pi, mu, var, max_ks)
86 h <- ntot * 1.5 / sum(g)
87 vx <- seq(1, 100) * (max_ks / 100)
88 ymax <- max(nc)
89 barplot(nc, space=0.25, offset=0, width=0.04, xlim=c(0, max_ks), ylim=c(0, ymax), col="lightpink1", border="lightpink3")
90 # Add x-axis.
91 axis(1)
92 color <- c('red', 'yellow','green','black','blue', 'darkorange' )
93 for (i in 1:length(mu))
94 {
95 lines(vx, g[,i] * h, lwd=2, col=color[i])
96 }
97 }
98
99 calculate_fitted_density <- function(pi, mu, var, max_ks)
100 {
101 comp <- length(pi)
102 var <- var/mu^2
103 mu <- log(mu)
104 # Calculate lognormal density.
105 vx <- seq(1, 100) * (max_ks / 100)
106 fx <- matrix(0, 100, comp)
107 for (i in 1:100)
108 {
109 for (j in 1:comp)
110 {
111 fx[i, j] <- pi[j] * dlnorm(vx[i], meanlog=mu[j], sdlog=(sqrt(var[j])))
112 if (is.nan(fx[i,j])) fx[i,j]<-0
113 }
114 }
115 return(fx)
116 }
117
118 # Read in the components data.
119 components_data <- read.delim(opt$components_input, header=TRUE)
120 # Get the number of components.
121 num_components <- get_num_components(components_data)
122
123 # Set pi, mu, var.
124 items <- get_pi_mu_var(components_data, num_components)
125 if (num_components == 1)
126 {
127 pi <- items[1]
128 mu <- items[2]
129 var <- items[3]
130 }
131 if (num_components == 2)
132 {
133 pi <- items[1:2]
134 mu <- items[3:4]
135 var <- items[5:6]
136 }
137 if (num_components == 3)
138 {
139 pi <- items[1:3]
140 mu <- items[4:6]
141 var <- items[7:9]
142 }
143 if (num_components == 4)
144 {
145 pi <- items[1:4]
146 mu <- items[5:8]
147 var <- items[9:12]
148 }
149 if (num_components == 5)
150 {
151 pi <- items[1:5]
152 mu <- items[6:10]
153 var <- items[11:15]
154 }
155 if (num_components == 6)
156 {
157 pi <- items[1:6]
158 mu <- items[7:12]
159 var <- items[13:18]
160 }
161
162 # Plot the output.
163 plot_ks(opt$kaks_input, opt$output, pi, mu, var)