annotate Functions/Kernel_function_form.R @ 1:cb8bac9d0d37 draft default tip

planemo upload commit 4ec21642211c1fb7427e4c98fdf0f4b9a3f8a185-dirty
author cristian
date Thu, 07 Sep 2017 10:21:45 -0400
parents 1535ffddeff4
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
1 ##############################################
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
2 # KDE function for single observation sequence #
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
3 ##############################################
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
4
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
5 # load packages
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
6 require(moments, quietly = TRUE)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
7 require(doParallel, quietly = TRUE)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
8
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
9
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
10 # get names for the columns of the table characterizing the peaks
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
11 col.names.peaks <- function()
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
12 {
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
13 v <- c("Number of modes", "Number of modes (5% excluded)",
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
14 "Number of modes (10% excluded)", "Skewness",
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
15 "Mode skewness", "Nonparametric skewness", "Q50 skewness",
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
16 "Absolute Q50 mode skewness", "Absolute Q80 mode skewness",
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
17 "Peak 1", "Probability Mass 1",
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
18 "Peak 2", "Probability Mass 2",
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
19 "Peak 3", "Probability Mass 3",
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
20 "Peak 4", "Probability Mass 4",
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
21 "Peak 5", "Probability Mass 5",
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
22 "Peak 6", "Probability Mass 6",
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
23 "Peak 7", "Probability Mass 7",
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
24 "Peak 8", "Probability Mass 8",
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
25 "Peak 9", "Probability Mass 9",
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
26 "Peak 10", "Probability Mass 10",
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
27 "Warning close modes",
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
28 "Number close modes",
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
29 "Modes (close modes excluded)",
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
30 "SD", "IQR 80", "IQR 90",
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
31 "Total number of sequences")
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
32 return(v)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
33 }
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
34
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
35
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
36 # get names for the columns of the table containing the boostrap results
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
37 col.names.bs <- function()
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
38 {
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
39 v <- c("Number of modes (NM)",
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
40 "% of samples with same NM",
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
41 "% of samples with more NM",
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
42 "% of samples with less NM",
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
43 "no. of samples with same NM",
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
44 "% BS samples excluded by prob. mass crit.",
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
45 "Warning CI")
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
46 return(v)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
47 }
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
48
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
49
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
50 # plot KDE for one set up CpGo/e ratios
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
51 # obs: CpGo/e ratios
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
52 # bs.cis: Is bootstrap done?
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
53 # t.name: Title of the plot
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
54 # t.sub: Text that is added below the title
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
55 # t.legend: Is a legend printed?
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
56 plot.KDE <- function(obs, t.name, bs.cis = FALSE, bstrp.reps = 1500, conf.lev = 0.95, t.sub = NULL, t.legend = TRUE, min.dist = 0.2, mode.mass = 0.01, band.width = 1.06)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
57 {
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
58 # determine directory where functions are located
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
59 cmdArgs <- commandArgs(trailingOnly = FALSE)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
60 str <- "--file="
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
61 match <- grep(str, cmdArgs)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
62 if (length(match) == 0) {
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
63 FCTN.DIR <- "../Galaxy/Functions"
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
64 } else {
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
65 path <- normalizePath( sub(str, "", cmdArgs[match]) )
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
66 FCTN.DIR <- file.path(dirname(path), "Functions")
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
67 }
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
68
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
69 # part 1: initialize parameters etc
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
70 # ---------------------------------
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
71
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
72 # table with names and number of peaks
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
73 v <- col.names.peaks()
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
74 tab1.m <- data.frame(matrix(NA, nrow = 1, ncol = length(v)))
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
75 names(tab1.m) <- col.names.peaks()
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
76 # parameters to set in any case
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
77 num.points <- 10 ^ 4 # number of points where to estimate the density
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
78 p.bw <- "nrd" # algorithm for the bandwidth selection, "nrd" for Scott's bandwith
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
79 use.seed <- TRUE
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
80 threshold.modes <- mode.mass
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
81 threshold.bs.ci <- 0.2 # only changes of +- threshold.bs.ci% in prob. mass is allowed for entering the CI calculation
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
82 # bootstrap with optional parameters
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
83 if (bs.cis) {
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
84 ncpus = max(1, detectCores()) # number of processsors
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
85 cl <- makeCluster(ncpus)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
86 registerDoParallel(cl)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
87 # table for the bootstrap
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
88 tab2.m <- data.frame(matrix(NA, nrow = 1, ncol = 7))
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
89 names(tab2.m) <- col.names.bs()
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
90 }
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
91
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
92
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
93 # part 2: estimate the KD and calculate probability masses
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
94 # --------------------------------------------------------
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
95
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
96 source( file.path(FCTN.DIR, "density_pm.R") )
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
97 estimate <- density_pm(obs, num.points, p.bw = band.width * sd(obs) / length(obs)^0.2, threshold.modes = threshold.modes)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
98 ker <- estimate$ker # kernel density
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
99 p <- estimate$peaks
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
100 v <- estimate$valleys
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
101 pm <- estimate$pm
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
102
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
103 # select for each peak a line type
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
104 num.pm <- length(p) # number of peaks
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
105 p5 <- estimate$p5
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
106 p10 <- estimate$p10
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
107 lty = rep(1, num.pm) #pm > 0.10
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
108 if (length(p10) > 0) { #pm < 0.10
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
109 lty[p10] <- 2
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
110 }
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
111 if (length(p5) > 0) { #pm < 0.05
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
112 lty[p5] <- 4
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
113 }
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
114 p.lty <- lty
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
115
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
116 # part 3: bootstrap
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
117 # -----------------
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
118
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
119 if (bs.cis == TRUE) {
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
120 # do bootstrapping
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
121 estimateB = foreach(j = 1 : bstrp.reps ) %dopar% {
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
122 if (use.seed == TRUE) {set.seed(j)}
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
123 source( file.path(FCTN.DIR, "density_pm.R") ) # As doParallel is used, the source code has to be included here
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
124 obs.boot = obs[sample(seq(obs), replace = TRUE)]
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
125 density_pm_boot(obs.boot, num.points = num.points, p.bw = band.width * sd(obs) / length(obs)^0.2, threshold.modes = threshold.modes)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
126 }
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
127 stopCluster(cl)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
128
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
129 # calculate CIs based on samples where the number of peaks of the KDE coincide with the original data
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
130 # and the probability mass of the peaks of the sample don not divert to strongly from the original data
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
131 # ... extract modes and pm for samples
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
132 num.pmB <- sapply(lapply(estimateB, "[[", "peaks"), length) # no of peaks before cleaning
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
133 num.peaks.ok <- num.pmB == num.pm # samples with same number of peaks before cleaning
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
134 pB <- t(sapply( estimateB[num.peaks.ok], "[[", "peaks") )
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
135 pmB <- t(sapply( estimateB[num.peaks.ok], "[[", "pm") )
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
136 if (num.pm == 1) { # put into right matrix form
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
137 pB <- t(pB)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
138 pmB <- t(pmB)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
139 }
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
140
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
141 # ... check if prob. mass changes too strongly for any mode
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
142 t.pmB <- t(pmB)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
143 keep.ub <- !apply(pm * (1 + threshold.bs.ci) < t.pmB, 2, any)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
144 keep.lb <- !apply(pm * (1 - threshold.bs.ci) > t.pmB, 2, any)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
145 mass.ok <- keep.ub & keep.lb
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
146 pB.cl <- as.matrix( pB[mass.ok, ] )
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
147
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
148 # ... determine CIs
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
149 q <- (1 - conf.lev) / 2
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
150 p.CI <- apply(pB.cl, 2, quantile, probs = c(q, 1 - q))
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
151 }
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
152
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
153
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
154 # part 4: plots
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
155 # -------------
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
156
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
157 #t.breaks <- seq(0, max(obs)*1.05, by = 0.03)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
158 t.breaks = 50
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
159 hist_data <- hist(obs, breaks = t.breaks, plot = FALSE)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
160 hist(obs, breaks = t.breaks, prob = TRUE, main = t.name,
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
161 # sub = paste("Gaussian kernel with band width", band.width),
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
162 xlab = "CpG o/e ratio",
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
163 col = grey(0.9), border = grey(0.6)) #, xlim = c(-0.05, max(obs)*1.1))
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
164 if (!is.null(t.sub)) {
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
165 mtext(t.sub)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
166 }
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
167 if (bs.cis == TRUE) { # CI
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
168 j <- 1 : ncol(p.CI)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
169 rect(ker$x[p.CI[1, j]], 0, ker$x[p.CI[2, j]],
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
170 15, density = 20, angle = 45 + (j - 1) * 90, col = "blue")
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
171 }
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
172 lines(ker, col = "red", lwd = 2) # density
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
173
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
174 # vertical lines at peaks
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
175 x.pos = ker$x[p]
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
176 ok <- c()
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
177 close <- c()
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
178 for (i in 1:length(x.pos)) {
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
179 b <- TRUE
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
180 if (i > 1) {
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
181 if (x.pos[i] - x.pos[i - 1] < min.dist) {
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
182 b <- FALSE
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
183 }
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
184 }
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
185 if (i < length(x.pos)) {
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
186 if (x.pos[i + 1] - x.pos[i] < min.dist) {
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
187 b <- FALSE
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
188 }
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
189 }
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
190
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
191 if (b) {
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
192 ok <- c(ok, i)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
193 } else {
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
194 close <- c(close, i)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
195 }
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
196 }
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
197 # peaks
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
198 if (length(ok) > 0) {
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
199 abline(v = ker$x[p][ok], col = "blue", lwd = 3, lty = p.lty)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
200 }
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
201 if (length(close) > 0) {
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
202 abline(v = ker$x[p][close], col = "orange", lwd = 3, lty = p.lty)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
203 }
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
204 # valleys
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
205 vals = ""
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
206 if (length(x.pos)>1) {
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
207 # ker$x[v][c(-1, -length(ker$x[v]))]
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
208 vals = ker$x[v][c(-1, -length(ker$x[v]))]
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
209 abline(v = vals, col = "black", lwd = 1)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
210 }
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
211
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
212 # legend
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
213 if(t.legend == TRUE) {
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
214 t.sym <- expression(""<="", ""<"", "">="")
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
215 thr = threshold.modes
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
216 if (thr >= 0.1) {
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
217 md.labs <- substitute(paste("Mode with PM ", sym3, " ", thr, sep = ""), list(thr = thr, sym3 = t.sym[[3]]))
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
218 } else {
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
219 md.labs <- substitute(paste("Mode with PM ", sym3, " 0.1", sep = ""), list(sym3 = t.sym[[3]]))
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
220 if (thr >= 0.05) {
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
221 md.labs <- c( md.labs, substitute(paste("Mode with ", thr, " ", sym1, " PM ", sym2, " 0.1", sep = ""), list(thr = thr, sym1 = t.sym[[1]], sym2 = t.sym[[2]])) )
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
222 } else {
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
223 md.labs <- c( md.labs, substitute(paste("Mode with 0.05 ", sym1, " PM ", sym2, " 0.1", sep = ""), list(sym1 = t.sym[[1]], sym2 = t.sym[[2]])),
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
224 substitute(paste("Mode with ", thr, sym1, " PM ", sym2, " 0.05"), list(thr = thr, sym1 = t.sym[[1]], sym2 = t.sym[[2]])))
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
225 }
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
226 }
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
227 legend("topright",
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
228 c(expression("Estimated density"), md.labs),
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
229 lty = c(1, 1, 2, 4), lwd = c(2, 3, 3, 3),
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
230 col = c("red", "blue", "blue", "blue"), bg = "white")
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
231 }
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
232
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
233
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
234 # part 5: return results
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
235 # ----------------------
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
236
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
237 # part 5 a): results table 1
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
238 # tbd (maybe): introduce maximum number of modes (10 right now)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
239 tab1.m[1, "Number of modes"] <- num.pm
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
240 tab1.m[1, 2] = num.pm - length(estimate$p5)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
241 tab1.m[1, 3] = num.pm - length(estimate$p10)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
242 for(j in 1 : 10){
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
243 if(num.pm < j){
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
244 tab1.m[1, j * 2 + 8] = c(" ")
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
245 tab1.m[1, j * 2 + 9] = c(" ")
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
246 } else{
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
247 tab1.m[1, j * 2 + 8] = ker$x[p][j]
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
248 tab1.m[1, j * 2 + 9] = pm[j]
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
249 }
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
250 }
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
251
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
252 # fill table 1 with descriptives
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
253 tab1.m[1, "Skewness"] <- skewness(obs)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
254 mode <- ker$x[ which.max(ker$y) ]
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
255 tab1.m[1, "Mode skewness"] <- (mean(obs) - mode) / sd(obs)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
256 tab1.m[1, "Nonparametric skewness"] <- (mean(obs) - median(obs)) / sd(obs)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
257 q <- quantile(obs, c(0.25, 0.5, 0.75))
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
258 tab1.m[1, "Q50 skewness"] <- (q[3] + q[1] - 2 * q[2]) / (q[3] - q[1])
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
259 tab1.m[1, "Absolute Q50 mode skewness"] <- (q[3] + q[1]) / 2 - mode
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
260 q <- quantile(obs, c(0.1, 0.5, 0.9))
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
261 tab1.m[1, "Absolute Q80 mode skewness"] <- (q[3] + q[1]) / 2 - mode
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
262 tab1.m[1, "SD"] <- sd(obs)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
263 tab1.m[1, "IQR 80"] <- diff(quantile(obs, c(0.1, 0.9)))
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
264 tab1.m[1, "IQR 90"] <- diff(quantile(obs, c(0.05, 0.95)))
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
265 tab1.m[1, "Total number of sequences"] = length(obs)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
266
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
267 # check if any peak is closer than a given threshold to any other
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
268 num.close.modes <- sum(diff(ker$x[p]) < min.dist)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
269 if ( any(diff(ker$x[p]) < min.dist) && (num.pm > 1) ) {
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
270 tab1.m[1, "Warning close modes"] <- "Modes too close"
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
271 tab1.m[1, "Number close modes"] <- num.close.modes
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
272 tab1.m[1, "Modes (close modes excluded)"] <- num.pm - num.close.modes
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
273 } else {
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
274 tab1.m[1, "Modes (close modes excluded)"] <- num.pm
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
275 }
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
276
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
277 # part 5 b): results table 2
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
278 if (bs.cis == TRUE) {
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
279 ker <- lapply( estimateB, "[[", "ker")
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
280 peaks <- lapply( estimateB, "[[", "peaks")
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
281 num.peaks <- c()
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
282 for (i in 1:length(peaks)) {
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
283 curr.peaks <- ker[[i]]$x[ peaks[[i]] ]
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
284 num.cl <- sum(diff(curr.peaks) < min.dist)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
285 num.peaks <- c(num.peaks, length(curr.peaks) - num.cl)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
286 }
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
287
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
288 # fill table 2 with stats on number of modes in bs samples
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
289 num <- num.pm - num.close.modes
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
290 tab2.m[1, "Number of modes (NM)"] <- num
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
291 tab2.m[1, "% of samples with same NM"] <- 100 * sum(num.peaks == num) / bstrp.reps # equal
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
292 tab2.m[1, "% of samples with more NM"] <- 100 * sum(num.peaks > num) / bstrp.reps # more
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
293 tab2.m[1, "% of samples with less NM"] <- 100 * sum(num.peaks < num) / bstrp.reps # less
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
294 if (num.pm > 1) {
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
295 tab2.m[1, "Warning CI"] <- "CI's may be unreliable"
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
296 }
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
297
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
298 tab2.m[1, "no. of samples with same NM"] <- sum(num.peaks == num)
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
299 tab2.m[1, "% BS samples excluded by prob. mass crit."] <- (1 - sum(mass.ok) / sum(num.peaks.ok)) * 100
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
300 }
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
301
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
302 # return the results
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
303 if (bs.cis){
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
304 return(list(tab.des = tab1.m, tab.bs = tab2.m, valleys = vals))
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
305 } else {
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
306 return(list(tab.des = tab1.m, valleys = vals))
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
307 }
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
308 }
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
309
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
310
1535ffddeff4 planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
cristian
parents:
diff changeset
311