comparison functions-all-clayton-12-13.r @ 9:c8f47481b704 draft

Uploaded
author modencode-dcc
date Thu, 17 Jan 2013 15:59:06 -0500
parents 5e6efd5f3567
children
comparison
equal deleted inserted replaced
8:6a4297fa2df1 9:c8f47481b704
1 # revised on 2-20-10
2 # - fix error in pass.structure: reverse rank.combined, so that big sig.value
3 # are ranked with small numbers (1, 2, ...)
4 # - fix error on get.ez.tt.all: get ez.cutoff from sorted e.z
5
6 #
7 # modified EM procedure to compute empirical CDF more precisely - 09/2009
8
9
10
11 # this file contains the functions for
12 # 1. computing the correspondence profile (upper rank intersection and derivatives)
13 # 2. inference of copula mixture model
14 #
15 # It also has functions for
16 # 1. reading peak caller results
17 # 2. processing and matching called peaks
18 # 3. plotting results
19
20
21 ################ read peak caller results
22
23 # process narrow peak format
24 # some peak callers may not report q-values, p-values or fold of enrichment
25 # need further process before comparison
26 #
27 # stop.exclusive: Is the basepair of peak.list$stop exclusive? In narrowpeak and broadpeak format they are exclusive.
28 # If it is exclusive, we need subtract peak.list$stop by 1 to avoid the same basepair being both a start and a stop of two
29 # adjacent peaks, which creates trouble for finding correct intersect
30 process.narrowpeak <- function(narrow.file, chr.size, half.width=NULL, summit="offset", stop.exclusive=T, broadpeak=F){
31
32
33 aa <- read.table(narrow.file)
34
35 if(broadpeak){
36 bb.ori <- data.frame(chr=aa$V1, start=aa$V2, stop=aa$V3, signal.value=aa$V7, p.value=aa$V8, q.value=aa$V9)
37 }else{
38 bb.ori <- data.frame(chr=aa$V1, start=aa$V2, stop=aa$V3, signal.value=aa$V7, p.value=aa$V8, q.value=aa$V9, summit=aa$V10)
39 }
40
41 if(summit=="summit"){
42 bb.ori$summit <- bb.ori$summit-bb.ori$start # change summit to offset to avoid error when concatenating chromosomes
43 }
44
45 bb <- concatenate.chr(bb.ori, chr.size)
46
47 #bb <- bb.ori
48
49 # remove the peaks that has the same start and stop value
50 bb <- bb[bb$start != bb$stop,]
51
52 if(stop.exclusive==T){
53 bb$stop <- bb$stop-1
54 }
55
56 if(!is.null(half.width)){
57 bb$start.ori <- bb$start
58 bb$stop.ori <- bb$stop
59
60 # if peak is narrower than the specified window, stay with its width
61 # otherwise chop wider peaks to specified width
62 width <- bb$stop-bb$start +1
63 is.wider <- width > 2*half.width
64
65 if(summit=="offset" | summit=="summit"){ # if summit is offset from start
66 bb$start[is.wider] <- bb$start.ori[is.wider] + bb$summit[is.wider]-half.width
67 bb$stop[is.wider] <- bb$start.ori[is.wider] + bb$summit[is.wider]+half.width
68 } else {
69 if(summit=="unknown"){
70 bb$start[is.wider] <- bb$start.ori[is.wider]+round(width[is.wider]/2) - half.width
71 bb$stop[is.wider] <- bb$start.ori[is.wider]+round(width[is.wider]/2) + half.width
72 }
73 }
74 }
75
76 bb <- clean.data(bb)
77 invisible(list(data.ori=bb.ori, data.cleaned=bb))
78 }
79
80 # clean data
81 # and concatenate chromosomes if needed
82 clean.data <- function(adata){
83
84 # remove the peaks that has the same start and stop value
85 adata <- adata[adata$start != adata$stop,]
86
87 # if some stops and starts are the same, need fix them
88 stop.in.start <- is.element(adata$stop, adata$start)
89 n.fix <- sum(stop.in.start)
90 if(n.fix >0){
91 print(paste("Fix", n.fix, "stops\n"))
92 adata$stop[stop.in.start] <- adata$stop[stop.in.start]-1
93 }
94
95 return(adata)
96 }
97
98 # concatenate peaks
99 # peaks: the dataframe to have all the peaks
100 # chr.file: the file to keep the length of each chromosome
101 # chr files should come from the species that the data is from
102 #concatenate.chr <- function(peaks, chr.size){
103
104 # chr.size <- read.table(chr.file)
105 # chr.o <- order(chr.size[,1])
106 # chr.size <- chr.size[chr.o,]
107 #
108 # chr.shift <- cumsum(c(0, chr.size[-nrow(chr.size),2]))
109 # chr.size.cum <- data.frame(chr=chr.size[,1], shift=chr.shift)
110 #
111 # for(i in 1:nrow(chr.size)){
112 # is.in <- as.character(peaks$chr) == as.character(chr.size.cum$chr[i])
113 # if(sum(is.in)>0){
114 # peaks[is.in,]$start <- peaks[is.in,]$start + chr.size.cum$shift[i]
115 # peaks[is.in,]$stop <- peaks[is.in,]$stop + chr.size.cum$shift[i]
116 # }
117 # }
118 #
119 # invisible(peaks)
120 #}
121
122
123
124
125 # concatenate peaks
126 # peaks: the dataframe to have all the peaks
127 # chr.file: the file to keep the length of each chromosome
128 # chr files should come from the species that the data is from
129 concatenate.chr <- function(peaks, chr.size){
130
131 # chr.size <- read.table(chr.file)
132 chr.o <- order(chr.size[,1])
133 chr.size <- chr.size[chr.o,]
134
135 chr.shift <- cumsum(c(0, chr.size[-nrow(chr.size),2]))
136 chr.size.cum <- data.frame(chr=chr.size[,1], shift=chr.shift)
137
138 peaks$start.ori <- peaks$start
139 peaks$stop.ori <- peaks$stop
140
141 for(i in 1:nrow(chr.size)){
142 is.in <- as.character(peaks$chr) == as.character(chr.size.cum$chr[i])
143 if(sum(is.in)>0){
144 peaks[is.in,]$start <- peaks[is.in,]$start + chr.size.cum$shift[i]
145 peaks[is.in,]$stop <- peaks[is.in,]$stop + chr.size.cum$shift[i]
146 }
147 }
148
149 invisible(peaks)
150 }
151
152
153 deconcatenate.chr <- function(peaks, chr.size){
154
155 chr.o <- order(chr.size[,1])
156 chr.size <- chr.size[chr.o,]
157
158 chr.shift <- cumsum(c(0, chr.size[-nrow(chr.size),2]))
159 chr.size.cum <- data.frame(chr=chr.size[,1], shift=chr.shift)
160
161 peaks$chr <- rep(NA, nrow(peaks))
162
163 for(i in 1:(nrow(chr.size.cum)-1)){
164 is.in <- peaks$start > chr.size.cum[i,2] & peaks$start <= chr.size.cum[i+1, 2]
165 if(sum(is.in)>0){
166 peaks[is.in,]$start <- peaks[is.in,]$start - chr.size.cum[i,2]
167 peaks[is.in,]$stop <- peaks[is.in,]$stop - chr.size.cum[i,2]+1
168 peaks[is.in,]$chr <- chr.size[i,1]
169 }
170 }
171
172 if(i == nrow(chr.size.cum)){
173 is.in <- peaks$start > chr.size.cum[i, 2]
174 if(sum(is.in)>0){
175 peaks[is.in,]$start <- peaks[is.in,]$start - chr.size.cum[i,2]
176 peaks[is.in,]$stop <- peaks[is.in,]$stop - chr.size.cum[i,2]+1
177 peaks[is.in,]$chr <- chr.size[i,1]
178 }
179 }
180
181 invisible(peaks)
182 }
183
184 ################ preprocessing peak calling output
185
186
187 #
188 # read two calling results and sort by peak starting locations,
189 # then find overlap between peaks
190 # INPUT:
191 # rep1: the 1st replicate
192 # rep2: the 2nd replicate
193 # OUTPUT:
194 # id1, id2: the labels for the identified peaks on the replicates
195 find.overlap <- function(rep1, rep2){
196
197 o1 <- order(rep1$start)
198 rep1 <- rep1[o1,]
199
200 o2 <- order(rep2$start)
201 rep2 <- rep2[o2,]
202
203 n1 <- length(o1)
204 n2 <- length(o2)
205
206 # assign common ID to peaks
207 id1 <- rep(0, n1) # ID assigned on rep1
208 id2 <- rep(0, n2) # ID assigned on rep2
209 id <- 1 # keep track common id's
210
211 # check if two replicates overlap with each other
212 i <- 1
213 j <- 1
214
215 while(i <= n1|| j <= n2){
216
217 # && (id1[n1] ==0 || id2[n2] ==0)
218
219 # if one list runs out
220 if(i > n1 && j < n2){
221
222 j <- j+1
223 id2[j] <- id
224 id <- id +1
225 next
226 } else{
227 if(j > n2 && i < n1){
228 i <- i+1
229 id1[i] <- id
230 id <- id +1
231 next
232 } else {
233 if(i >= n1 && j >=n2)
234 break
235 }
236 }
237
238 # if not overlap
239
240 if(!(rep1$start[i] <= rep2$stop[j] && rep2$start[j] <= rep1$stop[i])){
241
242 # at the start of loop, when both are not assigned an ID
243 # the one locates in front is assigned first
244 if(id1[i] ==0 && id2[j]==0){
245 if(rep1$stop[i] < rep2$stop[j]){
246 id1[i] <- id
247 } else {
248 id2[j] <- id
249 }
250 } else { # in the middle of the loop, when one is already assigned
251 # The one that has not assigned gets assigned
252 # if(id1[i] ==0){ # id1[i] is not assigned
253 # id1[i] <- id
254 # } else { # id2[i] is not assigned
255 # id2[j] <- id
256 # }
257
258 # order the id according to location
259 if(rep1$stop[i] <= rep2$stop[j]){
260 id1[i] <- max(id2[j], id1[i])
261 id2[j] <- id
262 } else {
263 if(rep1$stop[i] > rep2$stop[j]){
264 id2[j] <- max(id1[i], id2[j])
265 id1[i] <- id
266 }
267 }
268
269 }
270
271 id <- id +1
272
273 } else { # if overlap
274
275 if(id1[i] == 0 && id2[j] == 0){ # not assign label yet
276 id1[i] <- id
277 id2[j] <- id
278 id <- id +1
279 } else { # one peak is already assigned label, the other is 0
280
281 id1[i] <- max(id1[i], id2[j]) # this is a way to copy the label of the assigned peak without knowing which one is already assigned
282 id2[j] <- id1[i] # syncronize the labels
283 }
284
285 }
286
287 if(rep1$stop[i] < rep2$stop[j]){
288 i <- i+1
289 } else {
290 j <- j+1
291 }
292
293 }
294
295 invisible(list(id1=id1, id2=id2))
296
297 }
298
299 # Impute the missing significant value for the peaks called only on one replicate.
300 # value
301 # INPUT:
302 # rep1, rep2: the two peak calling output
303 # id1, id2: the IDs assigned by function find.overlap, vectors
304 # If id1[i]==id2[j], peak i on rep1 overlaps with peak j on rep2
305 # p.value.impute: the significant value to impute for the missing peaks
306 # OUTPUT:
307 # rep1, rep2: peaks ordered by the start locations with imputed peaks
308 # id1, id2: the IDs with imputed peaks
309 fill.missing.peaks <- function(rep1, rep2, id1, id2, p.value.impute){
310
311 # rep1 <- data.frame(chr=rep1$chr, start=rep1$start, stop=rep1$stop, sig.value=rep1$sig.value)
312 # rep2 <- data.frame(chr=rep2$chr, start=rep2$start, stop=rep2$stop, sig.value=rep2$sig.value)
313
314 o1 <- order(rep1$start)
315 rep1 <- rep1[o1,]
316
317 o2 <- order(rep2$start)
318 rep2 <- rep2[o2,]
319
320 entry.in1.not2 <- !is.element(id1, id2)
321 entry.in2.not1 <- !is.element(id2, id1)
322
323 if(sum(entry.in1.not2) > 0){
324
325 temp1 <- rep1[entry.in1.not2, ]
326
327 # impute sig.value
328 temp1$sig.value <- p.value.impute
329 temp1$signal.value <- p.value.impute
330 temp1$p.value <- p.value.impute
331 temp1$q.value <- p.value.impute
332
333 rep2.filled <- rbind(rep2, temp1)
334 id2.filled <- c(id2, id1[entry.in1.not2])
335 } else {
336 id2.filled <- id2
337 rep2.filled <- rep2
338 }
339
340 if(sum(entry.in2.not1) > 0){
341
342 temp2 <- rep2[entry.in2.not1, ]
343
344 # fill in p.values to 1
345 temp2$sig.value <- p.value.impute
346 temp2$signal.value <- p.value.impute
347 temp2$p.value <- p.value.impute
348 temp2$q.value <- p.value.impute
349
350
351 # append to the end
352 rep1.filled <- rbind(rep1, temp2)
353
354 id1.filled <- c(id1, id2[entry.in2.not1])
355 } else {
356 id1.filled <- id1
357 rep1.filled <- rep1
358 }
359
360 # sort rep1 and rep2 by the same id
361 o1 <- order(id1.filled)
362 rep1.ordered <- rep1.filled[o1, ]
363
364 o2 <- order(id2.filled)
365 rep2.ordered <- rep2.filled[o2, ]
366
367 invisible(list(rep1=rep1.ordered, rep2=rep2.ordered,
368 id1=id1.filled[o1], id2=id2.filled[o2]))
369 }
370
371 # Merge peaks with same ID on the same replicates
372 # (They are generated if two peaks on rep1 map to the same peak on rep2)
373 # need peak.list have 3 columns: start, stop and sig.value
374 merge.peaks <- function(peak.list, id){
375
376 i <- 1
377 j <- 1
378 dup.index <- c()
379 sig.value <- c()
380 start.new <- c()
381 stop.new <- c()
382 id.new <- c()
383
384 # original data
385 chr <- c()
386 start.ori <- c()
387 stop.ori <- c()
388
389 signal.value <- c()
390 p.value <- c()
391 q.value <- c()
392
393 while(i < length(id)){
394
395 if(id[i] == id[i+1]){
396 dup.index <- c(dup.index, i, i+1) # push on dup.index
397 } else {
398 if(length(dup.index)>0){ # pop from dup.index
399 sig.value[j] <- mean(peak.list$sig.value[unique(dup.index)]) # mean of -log(pvalue)
400 start.new[j] <- peak.list$start[min(dup.index)]
401 stop.new[j] <- peak.list$stop[max(dup.index)]
402 id.new[j] <- id[max(dup.index)]
403
404 signal.value[j] <- mean(peak.list$signal.value[unique(dup.index)]) # mean of -log(pvalue)
405 p.value[j] <- mean(peak.list$p.value[unique(dup.index)]) # mean of -log(pvalue)
406 q.value[j] <- mean(peak.list$q.value[unique(dup.index)]) # mean of -log(pvalue)
407
408 chr[j] <- as.character(peak.list$chr[min(dup.index)])
409 start.ori[j] <- peak.list$start.ori[min(dup.index)]
410 stop.ori[j] <- peak.list$stop.ori[max(dup.index)]
411
412 dup.index <- c()
413 } else { # nothing to pop
414 sig.value[j] <- peak.list$sig.value[i]
415 start.new[j] <- peak.list$start[i]
416 stop.new[j] <- peak.list$stop[i]
417 id.new[j] <- id[i]
418
419 signal.value[j] <- peak.list$signal.value[i]
420 p.value[j] <- peak.list$p.value[i]
421 q.value[j] <- peak.list$q.value[i]
422
423 chr[j] <- as.character(peak.list$chr[i])
424 start.ori[j] <- peak.list$start.ori[i]
425 stop.ori[j] <- peak.list$stop.ori[i]
426
427 }
428 j <- j+1
429 }
430 i <- i+1
431 }
432
433 data.new <- data.frame(id=id.new, sig.value=sig.value, start=start.new, stop=stop.new, signal.value=signal.value, p.value=p.value, q.value=q.value, chr=chr, start.ori=start.ori, stop.ori=stop.ori)
434 invisible(data.new)
435 }
436
437
438
439
440
441 # a wrap function to fill in missing peaks, merge peaks and impute significant values
442 # out1 and out2 are two peak calling outputs
443 pair.peaks <- function(out1, out2, p.value.impute=0){
444
445 aa <- find.overlap(out1, out2)
446 bb <- fill.missing.peaks(out1, out2, aa$id1, aa$id2, p.value.impute=0)
447
448 cc1 <- merge.peaks(bb$rep1, bb$id1)
449 cc2 <- merge.peaks(bb$rep2, bb$id2)
450
451 invisible(list(merge1=cc1, merge2=cc2))
452 }
453
454
455
456 # overlap.ratio is a parameter to define the percentage of overlap
457 # if overlap.ratio =0, 1 basepair overlap is counted as overlap
458 # if overlap.ratio between 0 and 1, it is the minimum proportion of
459 # overlap required to be called as a match
460 # it is computed as the overlap part/min(peak1.length, peak2.length)
461 pair.peaks.filter <- function(out1, out2, p.value.impute=0, overlap.ratio=0){
462
463 aa <- find.overlap(out1, out2)
464 bb <- fill.missing.peaks(out1, out2, aa$id1, aa$id2, p.value.impute=0)
465
466 cc1 <- merge.peaks(bb$rep1, bb$id1)
467 cc2 <- merge.peaks(bb$rep2, bb$id2)
468
469 frag12 <- cbind(cc1$start, cc1$stop, cc2$start, cc2$stop)
470
471 frag.ratio <- apply(frag12, 1, overlap.middle)
472
473 frag.ratio[cc1$sig.value==p.value.impute | cc2$sig.value==p.value.impute] <- 0
474
475 cc1$frag.ratio <- frag.ratio
476 cc2$frag.ratio <- frag.ratio
477
478 merge1 <- cc1[cc1$frag.ratio >= overlap.ratio,]
479 merge2 <- cc2[cc2$frag.ratio >= overlap.ratio,]
480
481 invisible(list(merge1=merge1, merge2=merge2))
482 }
483
484 # x[1], x[2] are the start and end of the first fragment
485 # and x[3] and x[4] are the start and end of the 2nd fragment
486 # If there are two fragments, we can find the overlap by ordering the
487 # start and stop of all the ends and find the difference between the middle two
488 overlap.middle <- function(x){
489
490 x.o <- x[order(x)]
491 f1 <- x[2]-x[1]
492 f2 <- x[4]-x[3]
493
494 f.overlap <- abs(x.o[3]-x.o[2])
495 f.overlap.ratio <- f.overlap/min(f1, f2)
496
497 return(f.overlap.ratio)
498 }
499
500
501
502 #######
503 ####### compute correspondence profile
504 #######
505
506 # compute upper rank intersection for one t
507 # tv: the upper percentile
508 # x is sorted by the order of paired variable
509 comp.uri <- function(tv, x){
510 n <- length(x)
511 qt <- quantile(x, prob=1-tv[1]) # tv[1] is t
512 # sum(x[1:ceiling(n*tv[2])] >= qt)/n/tv[2]- tv[1]*tv[2] #tv[2] is v
513 sum(x[1:ceiling(n*tv[2])] >= qt)/n
514
515 }
516
517 # compute the correspondence profile
518 # tt, vv: vector between (0, 1) for percentages
519 get.uri.2d <- function(x1, x2, tt, vv, spline.df=NULL){
520
521 o <- order(x1, x2, decreasing=T)
522
523 # sort x2 by the order of x1
524 x2.ordered <- x2[o]
525
526 tv <- cbind(tt, vv)
527 ntotal <- length(x1) # number of peaks
528
529 uri <- apply(tv, 1, comp.uri, x=x2.ordered)
530
531 # compute the derivative of URI vs t using small bins
532 uri.binned <- uri[seq(1, length(uri), by=4)]
533 tt.binned <- tt[seq(1, length(uri), by=4)]
534 uri.slope <- (uri.binned[2:(length(uri.binned))] - uri.binned[1:(length(uri.binned)-1)])/(tt.binned[2:(length(uri.binned))] - tt.binned[1:(length(tt.binned)-1)])
535
536 # smooth uri using spline
537 # first find where the jump is and don't fit the jump
538 # this is the index on the left
539 # jump.left.old <- which.max(uri[-1]-uri[-length(uri)])
540 short.list.length <- min(sum(x1>0)/length(x1), sum(x2>0)/length(x2))
541
542 if(short.list.length < max(tt)){
543 jump.left <- which(tt>short.list.length)[1]-1
544 } else {
545 jump.left <- which.max(tt)
546 }
547
548 # reversed.index <- seq(length(tt), 1, by=-1)
549 # nequal <- sum(uri[reversed.index]== tt[reversed.index])
550 # temp <- which(uri[reversed.index]== tt[reversed.index])[nequal]
551 # jump.left <- length(tt)-temp
552
553 if(jump.left < 6){
554 jump.left <- length(tt)
555 }
556
557
558 if(is.null(spline.df))
559 uri.spl <- smooth.spline(tt[1:jump.left], uri[1:jump.left], df=6.4)
560 else{
561 uri.spl <- smooth.spline(tt[1:jump.left], uri[1:jump.left], df=spline.df)
562 }
563 # predict the first derivative
564 uri.der <- predict(uri.spl, tt[1:jump.left], deriv=1)
565
566 invisible(list(tv=tv, uri=uri,
567 uri.slope=uri.slope, t.binned=tt.binned[2:length(uri.binned)],
568 uri.spl=uri.spl, uri.der=uri.der, jump.left=jump.left,
569 ntotal=ntotal))
570 }
571
572
573 # change the scale of uri from based on t (percentage) to n (number of peaks or basepairs)
574 # this is for plotting multiple pairwise URI's on the same plot
575 scale.t2n <- function(uri){
576
577 ntotal <- uri$ntotal
578 tv <- uri$tv*uri$ntotal
579 uri.uri <- uri$uri*uri$ntotal
580 jump.left <- uri$jump.left
581 uri.spl <- uri$uri.spl
582 uri.spl$x <- uri$uri.spl$x*uri$ntotal
583 uri.spl$y <- uri$uri.spl$y*uri$ntotal
584
585 t.binned <- uri$t.binned*uri$ntotal
586 uri.slope <- uri$uri.slope
587 uri.der <- uri$uri.der
588 uri.der$x <- uri$uri.der$x*uri$ntotal
589 uri.der$y <- uri$uri.der$y
590
591 uri.n <- list(tv=tv, uri=uri.uri, t.binned=t.binned, uri.slope=uri.slope, uri.spl=uri.spl, uri.der=uri.der, ntotal=ntotal, jump.left=jump.left)
592 return(uri.n)
593 }
594
595
596
597
598 # a wrapper for running URI for peaks from peak calling results
599 # both data1 and data2 are calling results in narrowpeak format
600 compute.pair.uri <- function(data.1, data.2, sig.value1="signal.value", sig.value2="signal.value", spline.df=NULL, overlap.ratio=0){
601
602 tt <- seq(0.01, 1, by=0.01)
603 vv <- tt
604
605 if(sig.value1=="signal.value"){
606 data.1.enrich <- data.frame(chr=data.1$chr, start.ori=data.1$start.ori, stop.ori=data.1$stop.ori, start=data.1$start, stop=data.1$stop, sig.value=data.1$signal.value, signal.value=data.1$signal.value, p.value=data.1$p.value, q.value=data.1$q.value)
607 } else {
608 if(sig.value1=="p.value"){
609 data.1.enrich <- data.frame(chr=data.1$chr, start.ori=data.1$start.ori, stop.ori=data.1$stop.ori, start=data.1$start, stop=data.1$stop, sig.value=data.1$p.value, signal.value=data.1$signal.value, p.value=data.1$p.value, q.value=data.1$q.value)
610 } else {
611 if(sig.value1=="q.value"){
612 data.1.enrich <- data.frame(chr=data.1$chr, start.ori=data.1$start.ori, stop.ori=data.1$stop.ori, start=data.1$start, stop=data.1$stop, sig.value=data.1$q.value, signal.value=data.1$signal.value, p.value=data.1$p.value, q.value=data.1$q.value)
613 }
614 }
615 }
616
617 if(sig.value2=="signal.value"){
618 data.2.enrich <- data.frame(chr=data.2$chr, start.ori=data.2$start.ori, stop.ori=data.2$stop.ori, start=data.2$start, stop=data.2$stop, sig.value=data.2$signal.value, signal.value=data.2$signal.value, p.value=data.2$p.value, q.value=data.2$q.value)
619 } else {
620 if(sig.value2=="p.value"){
621 data.2.enrich <- data.frame(chr=data.2$chr, start.ori=data.2$start.ori, stop.ori=data.2$stop.ori, start=data.2$start, stop=data.2$stop, sig.value=data.2$p.value, signal.value=data.2$signal.value, p.value=data.2$p.value, q.value=data.2$q.value)
622 } else {
623 if(sig.value2=="q.value"){
624 data.2.enrich <- data.frame(chr=data.2$chr, start.ori=data.2$start.ori, stop.ori=data.2$stop.ori, start=data.2$start, stop=data.2$stop, sig.value=data.2$q.value, signal.value=data.2$signal.value, p.value=data.2$p.value, q.value=data.2$q.value)
625 }
626 }
627 }
628
629 ### by peaks
630 # data12.enrich <- pair.peaks(data.1.enrich, data.2.enrich)
631 data12.enrich <- pair.peaks.filter(data.1.enrich, data.2.enrich, p.value.impute=0, overlap.ratio)
632 uri <- get.uri.2d(as.numeric(as.character(data12.enrich$merge1$sig.value)), as.numeric(as.character(data12.enrich$merge2$sig.value)), tt, vv, spline.df=spline.df)
633 uri.n <- scale.t2n(uri)
634
635 return(list(uri=uri, uri.n=uri.n, data12.enrich=data12.enrich, sig.value1=sig.value1, sig.value2=sig.value2))
636
637
638 }
639
640
641
642 # compute uri for matched sample
643 get.uri.matched <- function(data12, df=10){
644
645 tt <- seq(0.01, 1, by=0.01)
646 vv <- tt
647 uri <- get.uri.2d(data12$sample1$sig.value, data12$sample2$sig.value, tt, vv, spline.df=df)
648
649 # change scale from t to n
650 uri.n <- scale.t2n(uri)
651
652 return(list(uri=uri, uri.n=uri.n))
653
654 }
655
656 # map.uv is a pair of significant values corresponding to specified consistency FDR
657 # assuming values in map.uv and qvalue are linearly related
658 # data.set is the original data set
659 # sig.value is the name of the significant value in map.uv, say enrichment
660 # nominal.value is the one we want to map to, say q-value
661 #
662 map.sig.value <- function(data.set, map.uv, nominal.value){
663
664 index.nominal <- which(names(data.set$merge1)==nominal.value)
665 nentry <- nrow(map.uv)
666 map.nominal <- rbind(map.uv[, c("sig.value1", "sig.value2")])
667
668 for(i in 1:nentry){
669
670 map.nominal[i, "sig.value1"] <- data.set$merge1[unique(which.min(abs(data.set$merge1$sig.value-map.uv[i, "sig.value1"]))), index.nominal]
671 map.nominal[i, "sig.value2"] <- data.set$merge2[unique(which.min(abs(data.set$merge2$sig.value-map.uv[i, "sig.value2"]))), index.nominal]
672 }
673
674 invisible(map.nominal)
675 }
676
677
678 ############### plot correspondence profile
679
680 # plot multiple comparison wrt one template
681 # uri.list contains the total number of peaks
682 # plot.missing=F: not plot the missing points on the right
683 plot.uri.group <- function(uri.n.list, plot.dir, file.name=NULL, legend.txt, xlab.txt="num of significant peaks", ylab.txt="num of peaks in common", col.start=0, col.txt=NULL, plot.missing=F, title.txt=NULL){
684
685 if(is.null(col.txt))
686 col.txt <- c("black", "red", "purple", "green", "blue", "cyan", "magenta", "orange", "grey")
687
688 n <- length(uri.n.list)
689
690 ntotal <- c()
691 for(i in 1:n)
692 ntotal[i] <- uri.n.list[[i]]$ntotal
693
694 jump.left <- c()
695 jump.left.der <- c()
696 ncommon <- c()
697 for(i in 1:n){
698 # jump.left[i] <- which.max(uri.n.list[[i]]$uri[-1]-uri.n.list[[i]]$uri[-length(uri.n.list[[i]]$uri)])
699 # if(jump.left[i] < 6)
700 # jump.left[i] <- length(uri.n.list[[i]]$uri)
701
702 ## reversed.index <- seq(length(uri.n.list[[i]]$tv[,1]), 1, by=-1)
703 ## nequal <- sum(uri.n.list[[i]]$uri[reversed.index]== uri.n.list[[i]]$tv[reversed.index,1])
704 ## temp <- which(uri.n.list[[i]]$uri[reversed.index]== uri.n.list[[i]]$tv[reversed.index,1])[nequal]
705 ## jump.left[i] <- length(uri.n.list[[i]]$tv[,1])-temp
706 ##print(uri.n.list[[i]]$uri)
707 ##print(uri.n.list[[i]]$tv[,1])
708 ## jump.left[i] <- uri.n.list[[i]]$jump.left
709
710 # jump.left.der[i] <- sum(uri.n.list[[i]]$t.binned < uri.n.list[[i]]$uri.der$x[length(uri.n.list[[i]]$uri.der$x)])
711
712 jump.left[i] <- uri.n.list[[i]]$jump.left
713 jump.left.der[i] <- jump.left[i]
714 ncommon[i] <- uri.n.list[[i]]$tv[jump.left[i],1]
715 }
716
717
718 if(plot.missing){
719 max.peak <- max(ntotal)
720 } else {
721 max.peak <- max(ncommon)*1.05
722 }
723
724 if(!is.null(file.name)){
725 postscript(paste(plot.dir, "uri.", file.name, sep=""))
726 par(mfrow=c(1,1), mar=c(5,5,4,2))
727 }
728
729 plot(uri.n.list[[1]]$tv[,1], uri.n.list[[1]]$uri, type="n", xlab=xlab.txt, ylab=ylab.txt, xlim=c(0, max.peak), ylim=c(0, max.peak), cex.lab=2)
730
731 for(i in 1:n){
732
733 if(plot.missing){
734 points(uri.n.list[[i]]$tv[,1], uri.n.list[[i]]$uri, col=col.txt[i+col.start], cex=0.5 )
735 } else {
736 points(uri.n.list[[i]]$tv[1:jump.left[i],1], uri.n.list[[i]]$uri[1:jump.left[i]], col=col.txt[i+col.start], cex=0.5)
737 }
738 lines(uri.n.list[[i]]$uri.spl, col=col.txt[i+col.start], lwd=4)
739 }
740 abline(coef=c(0,1), lty=3)
741 legend(0, max.peak, legend=legend.txt, col=col.txt[(col.start+1):length(col.txt)], lty=1, lwd=3, cex=2)
742 if(!is.null(title))
743 title(title.txt)
744
745 if(!is.null(file.name)){
746 dev.off()
747 }
748
749 if(!is.null(file.name)){
750 postscript(paste(plot.dir, "duri.", file.name, sep=""))
751 par(mfrow=c(1,1), mar=c(5,5,4,2))
752 }
753 plot(uri.n.list[[1]]$t.binned, uri.n.list[[1]]$uri.slope, type="n", xlab=xlab.txt, ylab="slope", xlim=c(0, max.peak), ylim=c(0, 1.5), cex.lab=2)
754
755 for(i in 1:n){
756 # if(plot.missing){
757 # points(uri.n.list[[i]]$t.binned, uri.n.list[[i]]$uri.slope, col=col.txt[i+col.start], cex=0.5)
758 # } else {
759 # points(uri.n.list[[i]]$t.binned[1:jump.left.der[i]], uri.n.list[[i]]$uri.slope[1:jump.left.der[i]], col=col.txt[i+col.start], cex=0.5)
760 # }
761 lines(uri.n.list[[i]]$uri.der, col=col.txt[i+col.start], lwd=4)
762 }
763 abline(h=1, lty=3)
764 legend(0.5*max.peak, 1.5, legend=legend.txt, col=col.txt[(col.start+1):length(col.txt)], lty=1, lwd=3, cex=2)
765
766 if(!is.null(title))
767 title(title.txt)
768
769 if(!is.null(file.name)){
770 dev.off()
771 }
772
773 }
774
775
776
777 #######################
778 ####################### copula fitting for matched peaks
779 #######################
780
781 # estimation from mixed copula model
782
783 # 4-5-09
784 # A nonparametric estimation of mixed copula model
785
786
787 # updated
788
789 # c1, c2, f1, f2, g1, g2 are vectors
790 # c1*f1*g1 and c2*f2*g2 are copula densities for the two components
791 # xd1 and yd1 are the values of marginals for the first component
792 # xd2 and yd2 are the values of marginals for the 2nd component
793 #
794 # ez is the prob for being in the consistent group
795 get.ez <- function(p, c1, c2, xd1, yd1, xd2, yd2){
796
797 return(p*c1*xd1*yd1/(p*c1*xd1*yd1 + (1-p)*c2*xd2*yd2))
798 }
799
800 # checked
801
802 # this is C_12 not the copula density function c=C_12 * f1* f2
803 # since nonparametric estimation is used here for f1 and f2, which
804 # are constant throughout the iterations, we don't need them for optimization
805 #
806 # bivariate gaussian copula function
807 # t and s are vectors of same length, both are percentiles
808 # return a vector
809 gaussian.cop.den <- function(t, s, rho){
810
811 A <- qnorm(t)^2 + qnorm(s)^2
812 B <- qnorm(t)*qnorm(s)
813
814 loglik <- -log(1-rho^2)/2 - rho/(2*(1-rho^2))*(rho*A-2*B)
815
816 return(exp(loglik))
817 }
818
819 clayton.cop.den <- function(t, s, rho){
820
821 if(rho > 0)
822 return(exp(log(rho+1)-(rho+1)*(log(t)+log(s))-(2+1/rho)*log(t^(-rho) + s^(-rho)-1)))
823
824 if(rho==0)
825 return(1)
826
827 if(rho<0)
828 stop("Incorrect Clayton copula coefficient")
829
830 }
831
832
833 # checked
834 # estimate rho from Gaussian copula
835 mle.gaussian.copula <- function(t, s, e.z){
836
837 # reparameterize to bound from rho=+-1
838 l.c <- function(rho, t, s, e.z){
839 # cat("rho=", rho, "\n")
840 sum(e.z*log(gaussian.cop.den(t, s, rho)))}
841
842 rho.max <- optimize(f=l.c, c(-0.998, 0.998), maximum=T, tol=0.00001, t=t, s=s, e.z=e.z)
843
844 #print(rho.max$m)
845
846 #cat("cor=", cor(qnorm(t)*e.z, qnorm(s)*e.z), "\t", "rho.max=", rho.max$m, "\n")
847 # return(sign(rho.max$m)/(1+rho.max$m))
848 return(rho.max$m)
849 }
850
851
852 # estimate mle from Clayton copula,
853 mle.clayton.copula <- function(t, s, e.z){
854
855 l.c <- function(rho, t, s, e.z){
856 lc <- sum(e.z*log(clayton.cop.den(t, s, rho)))
857 # cat("rho=", rho, "\t", "l.c=", lc, "\n")
858 return(lc)
859 }
860
861 rho.max <- optimize(f=l.c, c(0.1, 20), maximum=T, tol=0.00001, t=t, s=s, e.z=e.z)
862
863 return(rho.max$m)
864 }
865
866
867
868 # updated
869 # mixture likelihood of two gaussian copula
870 # nonparametric and ranked transformed
871 loglik.2gaussian.copula <- function(x, y, p, rho1, rho2, x.mar, y.mar){
872
873 px.1 <- get.pdf.cdf(x, x.mar$f1)
874 px.2 <- get.pdf.cdf(x, x.mar$f2)
875 py.1 <- get.pdf.cdf(y, y.mar$f1)
876 py.2 <- get.pdf.cdf(y, y.mar$f2)
877
878 c1 <- gaussian.cop.den(px.1$cdf, py.1$cdf, rho1)
879 c2 <- gaussian.cop.den(px.2$cdf, py.2$cdf, rho2)
880
881 sum(log(p*c1*px.1$pdf*py.1$pdf + (1-p)*c2*px.2$pdf*py.2$pdf))
882 }
883
884 loglik.2copula <- function(x, y, p, rho1, rho2, x.mar, y.mar, copula.txt){
885
886 px.1 <- pdf.cdf$px.1
887 px.2 <- pdf.cdf$px.2
888 py.1 <- pdf.cdf$py.1
889 py.2 <- pdf.cdf$py.2
890
891 if(copula.txt=="gaussian"){
892 c1 <- gaussian.cop.den(px.1$cdf, py.1$cdf, rho1)
893 c2 <- gaussian.cop.den(px.2$cdf, py.2$cdf, rho2)
894 } else {
895 if(copula.txt=="clayton"){
896 c1 <- clayton.cop.den(px.1$cdf, py.1$cdf, rho1)
897 c2 <- clayton.cop.den(px.2$cdf, py.2$cdf, rho2)
898 }
899 }
900 sum(log(p*c1*px.1$pdf*py.1$pdf + (1-p)*c2*px.2$pdf*py.2$pdf))
901 }
902
903
904 # estimate the marginals of each component using histogram estimator in EM
905 # return the density, breaks, and cdf of the histogram estimator
906 est.mar.hist <- function(x, e.z, breaks){
907
908 binwidth <- c()
909 nbin <- length(breaks)-1
910 nx <- length(x)
911
912 # the histogram
913 x1.pdf <- c()
914 x2.pdf <- c()
915 x1.cdf <- c()
916 x2.cdf <- c()
917
918 # the pdf for each point
919 x1.pdf.value <- rep(NA, nx)
920 x2.pdf.value <- rep(NA, nx)
921
922 x1.cdf.value <- rep(NA, nx)
923 x2.cdf.value <- rep(NA, nx)
924
925 for(i in 1:nbin){
926
927 binwidth[i] <- breaks[i+1] - breaks[i]
928 if(i < nbin)
929 in.bin <- x>= breaks[i] & x < breaks[i+1]
930 else # last bin
931 in.bin <- x>= breaks[i] & x <=breaks[i+1]
932
933 # each bin add one observation to avoid empty bins
934 # multiple (nx+nbin)/(nx+nbin+1) to avoid blowup when looking up for
935 # quantiles
936 x1.pdf[i] <- (sum(e.z[in.bin])+1)/(sum(e.z)+nbin)/binwidth[i]*(nx+nbin)/(nx+nbin+1)
937 x2.pdf[i] <- (sum(1-e.z[in.bin])+1)/(sum(1-e.z)+nbin)/binwidth[i]*(nx+nbin)/(nx+nbin+1)
938
939
940 # x1.pdf[i] <- sum(e.z[in.bin])/sum(e.z)/binwidth[i]*nx/(nx+1)
941 # x2.pdf[i] <- sum(1-e.z[in.bin])/sum(1-e.z)/binwidth[i]*nx/(nx+1)
942
943 # treat each bin as a value for a discrete variable
944 # x1.cdf[i] <- sum(x1.pdf[1:i]*binwidth[1:i])
945 # x2.cdf[i] <- sum(x2.pdf[1:i]*binwidth[1:i])
946
947
948 # cumulative density before reaching i
949 if(i>1){
950 x1.cdf[i] <- sum(x1.pdf[1:(i-1)]*binwidth[1:(i-1)])
951 x2.cdf[i] <- sum(x2.pdf[1:(i-1)]*binwidth[1:(i-1)])
952 } else{
953 x1.cdf[i] <- 0
954 x2.cdf[i] <- 0
955 }
956
957 # make a vector of nx to store the values of pdf and cdf for each x
958 # this will speed up the computation dramatically
959 x1.pdf.value[in.bin] <- x1.pdf[i]
960 x2.pdf.value[in.bin] <- x2.pdf[i]
961
962 x1.cdf.value[in.bin] <- x1.cdf[i] + x1.pdf[i]*(x[in.bin]-breaks[i])
963 x2.cdf.value[in.bin] <- x2.cdf[i] + x2.pdf[i]*(x[in.bin]-breaks[i])
964 }
965
966 # x1.cdf <- cumsum(x1.pdf*binwidth)
967 # x2.cdf <- cumsum(x2.pdf*binwidth)
968
969 f1 <-list(breaks=breaks, density=x1.pdf, cdf=x1.cdf)
970 f2 <-list(breaks=breaks, density=x2.pdf, cdf=x2.cdf)
971
972 f1.value <- list(pdf=x1.pdf.value, cdf=x1.cdf.value)
973 f2.value <- list(pdf=x2.pdf.value, cdf=x2.cdf.value)
974
975 return(list(f1=f1, f2=f2, f1.value=f1.value, f2.value=f2.value))
976 }
977
978 # estimate the marginal cdf from rank
979 est.cdf.rank <- function(x, conf.z){
980
981 # add 1 to prevent blow up
982 x1.cdf <- rank(x[conf.z==1])/(length(x[conf.z==1])+1)
983
984 x2.cdf <- rank(x[conf.z==0])/(length(x[conf.z==0])+1)
985
986 return(list(cdf1=x1.cdf, cdf2=x2.cdf))
987 }
988
989 # df is a density function with fields: density, cdf and breaks, x is a scalar
990 get.pdf <- function(x, df){
991
992 if(x < df$breaks[1])
993 cat("x is out of the range of df\n")
994
995 index <- which(df$breaks >= x)[1]
996
997 if(index==1)
998 index <- index +1
999 return(df$density[index-1])
1000 }
1001
1002 # get cdf from histgram estimator for a single value
1003 get.cdf <- function(x, df){
1004
1005 index <- which(df$breaks >= x)[1]
1006 if(index==1)
1007 index <- index +1
1008 return(df$cdf[index-1])
1009 }
1010
1011 # df is a density function with fields: density, cdf and breaks
1012 get.pdf.cdf <- function(x.vec, df){
1013
1014 x.pdf <- sapply(x.vec, get.pdf, df=df)
1015 x.cdf <- sapply(x.vec, get.cdf, df=df)
1016 return(list(cdf=x.cdf, pdf=x.pdf))
1017 }
1018
1019 # E-step
1020 # x and y are the original observations or ranks
1021 # rho1 and rho2 are the parameters of each copula
1022 # f1, f2, g1, g2 are functions, each is a histogram
1023 e.step.2gaussian <- function(x, y, p, rho1, rho2, x.mar, y.mar){
1024
1025 # get pdf and cdf of each component from functions in the corresponding component
1026 px.1 <- get.pdf.cdf(x, x.mar$f1)
1027 px.2 <- get.pdf.cdf(x, x.mar$f2)
1028 py.1 <- get.pdf.cdf(y, y.mar$f1)
1029 py.2 <- get.pdf.cdf(y, y.mar$f2)
1030
1031 c1 <- gaussian.cop.den(px.1$cdf, py.1$cdf, rho1)
1032 c2 <- gaussian.cop.den(px.2$cdf, py.2$cdf, rho2)
1033
1034 return(get.ez(p, c1, c2, px.1$pdf, py.1$pdf, px.2$pdf, py.2$pdf))
1035 }
1036
1037 # E-step
1038 # rho1 and rho2 are the parameters of each copula
1039 e.step.2copula <- function(x, y, p, rho1, rho2, x.mar, y.mar, copula.txt){
1040
1041 # get pdf and cdf of each component from functions in the corresponding component
1042 px.1 <- get.pdf.cdf(x, x.mar$f1)
1043 px.2 <- get.pdf.cdf(x, x.mar$f2)
1044 py.1 <- get.pdf.cdf(y, y.mar$f1)
1045 py.2 <- get.pdf.cdf(y, y.mar$f2)
1046
1047 if(copula.txt=="gaussian"){
1048 c1 <- gaussian.cop.den(px.1$cdf, py.1$cdf, rho1)
1049 c2 <- gaussian.cop.den(px.2$cdf, py.2$cdf, rho2)
1050 } else {
1051 if(copula.txt=="clayton"){
1052 c1 <- clayton.cop.den(px.1$cdf, py.1$cdf, rho1)
1053 c2 <- clayton.cop.den(px.2$cdf, py.2$cdf, rho2)
1054 }
1055 }
1056 return(get.ez(p, c1, c2, px.1$pdf, py.1$pdf, px.2$pdf, py.2$pdf))
1057 }
1058
1059
1060
1061
1062 # M-step
1063 m.step.2gaussian <- function(x, y, e.z, breaks){
1064
1065 # compute f1, f2, g1 and g2
1066 x.mar <- est.mar.hist(x, e.z, breaks)
1067 y.mar <- est.mar.hist(y, e.z, breaks)
1068
1069 px.1 <- get.pdf.cdf(x, x.mar$f1)
1070 px.2 <- get.pdf.cdf(x, x.mar$f2)
1071 py.1 <- get.pdf.cdf(y, y.mar$f1)
1072 py.2 <- get.pdf.cdf(y, y.mar$f2)
1073
1074 rho1 <- mle.gaussian.copula(px.1$cdf, py.1$cdf, e.z)
1075 rho2 <- mle.gaussian.copula(px.2$cdf, py.2$cdf, 1-e.z)
1076
1077 p <- sum(e.z)/length(e.z)
1078
1079 return(list(p=p, rho1=rho1, rho2=rho2, x.mar=x.mar, y.mar=y.mar))
1080 }
1081
1082 m.step.2copula <- function(x, y, e.z, breaks, copula.txt){
1083
1084 # compute f1, f2, g1 and g2
1085 x.mar <- est.mar.hist(x, e.z, breaks)
1086 y.mar <- est.mar.hist(y, e.z, breaks)
1087
1088 px.1 <- get.pdf.cdf(x, x.mar$f1)
1089 px.2 <- get.pdf.cdf(x, x.mar$f2)
1090 py.1 <- get.pdf.cdf(y, y.mar$f1)
1091 py.2 <- get.pdf.cdf(y, y.mar$f2)
1092
1093 if(copula.txt=="gaussian"){
1094 rho1 <- mle.gaussian.copula(px.1$cdf, py.1$cdf, e.z)
1095 rho2 <- mle.gaussian.copula(px.2$cdf, py.2$cdf, 1-e.z)
1096 } else {
1097 if(copula.txt=="clayton"){
1098 rho1 <- mle.clayton.copula(px.1$cdf, py.1$cdf, e.z)
1099 rho2 <- mle.clayton.copula(px.2$cdf, py.2$cdf, 1-e.z)
1100 }
1101 }
1102
1103 p <- sum(e.z)/length(e.z)
1104
1105 return(list(p=p, rho1=rho1, rho2=rho2, x.mar=x.mar, y.mar=y.mar))
1106 }
1107
1108
1109
1110 # E-step: pass values
1111 # x and y are the original observations or ranks
1112 # rho1 and rho2 are the parameters of each copula
1113 # f1, f2, g1, g2 are functions, each is a histogram
1114 e.step.2gaussian.value <- function(x, y, p, rho1, rho2, pdf.cdf){
1115
1116 c1 <- gaussian.cop.den(pdf.cdf$px.1$cdf, pdf.cdf$py.1$cdf, rho1)
1117 c2 <- gaussian.cop.den(pdf.cdf$px.2$cdf, pdf.cdf$py.2$cdf, rho2)
1118
1119 e.z <- get.ez(p, c1, c2, pdf.cdf$px.1$pdf, pdf.cdf$py.1$pdf,
1120 pdf.cdf$px.2$pdf, pdf.cdf$py.2$pdf)
1121 return(e.z)
1122 }
1123
1124
1125 e.step.2copula.value <- function(x, y, p, rho1, rho2, pdf.cdf, copula.txt){
1126
1127 if(copula.txt =="gaussian"){
1128 c1 <- gaussian.cop.den(pdf.cdf$px.1$cdf, pdf.cdf$py.1$cdf, rho1)
1129 c2 <- gaussian.cop.den(pdf.cdf$px.2$cdf, pdf.cdf$py.2$cdf, rho2)
1130 } else {
1131 if(copula.txt =="clayton"){
1132 c1 <- clayton.cop.den(pdf.cdf$px.1$cdf, pdf.cdf$py.1$cdf, rho1)
1133 c2 <- clayton.cop.den(pdf.cdf$px.2$cdf, pdf.cdf$py.2$cdf, rho2)
1134 }
1135 }
1136
1137 e.z <- get.ez(p, c1, c2, pdf.cdf$px.1$pdf, pdf.cdf$py.1$pdf,
1138 pdf.cdf$px.2$pdf, pdf.cdf$py.2$pdf)
1139 return(e.z)
1140 }
1141
1142
1143 # M-step: pass values
1144 m.step.2gaussian.value <- function(x, y, e.z, breaks, fix.rho2){
1145
1146 # compute f1, f2, g1 and g2
1147 x.mar <- est.mar.hist(x, e.z, breaks)
1148 y.mar <- est.mar.hist(y, e.z, breaks)
1149
1150 # px.1 <- get.pdf.cdf(x, x.mar$f1)
1151 # px.2 <- get.pdf.cdf(x, x.mar$f2)
1152 # py.1 <- get.pdf.cdf(y, y.mar$f1)
1153 # py.2 <- get.pdf.cdf(y, y.mar$f2)
1154
1155 px.1 <- x.mar$f1.value
1156 px.2 <- x.mar$f2.value
1157 py.1 <- y.mar$f1.value
1158 py.2 <- y.mar$f2.value
1159
1160 rho1 <- mle.gaussian.copula(px.1$cdf, py.1$cdf, e.z)
1161
1162 if(!fix.rho2)
1163 rho2 <- mle.gaussian.copula(px.2$cdf, py.2$cdf, 1-e.z)
1164 else
1165 rho2 <- 0
1166
1167 p <- sum(e.z)/length(e.z)
1168
1169 pdf.cdf <- list(px.1=px.1, px.2=px.2, py.1=py.1, py.2=py.2)
1170
1171 return(list(p=p, rho1=rho1, rho2=rho2, x.mar=x.mar, y.mar=y.mar,
1172 pdf.cdf=pdf.cdf))
1173 }
1174
1175 m.step.2gaussian.value2 <- function(x, y, e.z, breaks, fix.rho2, x.mar, y.mar){
1176
1177 # compute f1, f2, g1 and g2
1178 # x.mar <- est.mar.hist(x, e.z, breaks)
1179 # y.mar <- est.mar.hist(y, e.z, breaks)
1180
1181 # px.1 <- get.pdf.cdf(x, x.mar$f1)
1182 # px.2 <- get.pdf.cdf(x, x.mar$f2)
1183 # py.1 <- get.pdf.cdf(y, y.mar$f1)
1184 # py.2 <- get.pdf.cdf(y, y.mar$f2)
1185
1186 px.1 <- x.mar$f1.value
1187 px.2 <- x.mar$f2.value
1188 py.1 <- y.mar$f1.value
1189 py.2 <- y.mar$f2.value
1190
1191 rho1 <- mle.gaussian.copula(px.1$cdf, py.1$cdf, e.z)
1192
1193 if(!fix.rho2)
1194 rho2 <- mle.gaussian.copula(px.2$cdf, py.2$cdf, 1-e.z)
1195 else
1196 rho2 <- 0
1197
1198 p <- sum(e.z)/length(e.z)
1199
1200 pdf.cdf <- list(px.1=px.1, px.2=px.2, py.1=py.1, py.2=py.2)
1201
1202 return(list(p=p, rho1=rho1, rho2=rho2, x.mar=x.mar, y.mar=y.mar,
1203 pdf.cdf=pdf.cdf))
1204 }
1205
1206
1207
1208 m.step.2copula.value <- function(x, y, e.z, breaks, fix.rho2, copula.txt){
1209
1210 # compute f1, f2, g1 and g2
1211 x.mar <- est.mar.hist(x, e.z, breaks)
1212 y.mar <- est.mar.hist(y, e.z, breaks)
1213
1214 # px.1 <- get.pdf.cdf(x, x.mar$f1)
1215 # px.2 <- get.pdf.cdf(x, x.mar$f2)
1216 # py.1 <- get.pdf.cdf(y, y.mar$f1)
1217 # py.2 <- get.pdf.cdf(y, y.mar$f2)
1218
1219 px.1 <- x.mar$f1.value
1220 px.2 <- x.mar$f2.value
1221 py.1 <- y.mar$f1.value
1222 py.2 <- y.mar$f2.value
1223
1224 if(copula.txt=="gaussian"){
1225 rho1 <- mle.gaussian.copula(px.1$cdf, py.1$cdf, e.z)
1226
1227 if(!fix.rho2)
1228 rho2 <- mle.gaussian.copula(px.2$cdf, py.2$cdf, 1-e.z)
1229 else
1230 rho2 <- 0
1231 } else {
1232
1233 if(copula.txt=="clayton"){
1234 rho1 <- mle.clayton.copula(px.1$cdf, py.1$cdf, e.z)
1235
1236 if(!fix.rho2)
1237 rho2 <- mle.clayton.copula(px.2$cdf, py.2$cdf, 1-e.z)
1238 else
1239 rho2 <- 0
1240 }
1241 }
1242
1243 p <- sum(e.z)/length(e.z)
1244
1245 pdf.cdf <- list(px.1=px.1, px.2=px.2, py.1=py.1, py.2=py.2)
1246
1247 return(list(p=p, rho1=rho1, rho2=rho2, x.mar=x.mar, y.mar=y.mar,
1248 pdf.cdf=pdf.cdf))
1249 }
1250
1251
1252
1253
1254 # updated
1255 # mixture likelihood of two gaussian copula
1256 # nonparametric and ranked transformed
1257 loglik.2gaussian.copula.value <- function(x, y, p, rho1, rho2, pdf.cdf){
1258
1259 px.1 <- pdf.cdf$px.1
1260 px.2 <- pdf.cdf$px.2
1261 py.1 <- pdf.cdf$py.1
1262 py.2 <- pdf.cdf$py.2
1263
1264 c1 <- gaussian.cop.den(px.1$cdf, py.1$cdf, rho1)
1265 c2 <- gaussian.cop.den(px.2$cdf, py.2$cdf, rho2)
1266
1267 sum(log(p*c1*px.1$pdf*py.1$pdf + (1-p)*c2*px.2$pdf*py.2$pdf))
1268 }
1269
1270
1271
1272 # updated
1273 # mixture likelihood of two gaussian copula
1274 # nonparametric and ranked transformed
1275 loglik.2copula.value <- function(x, y, p, rho1, rho2, pdf.cdf, copula.txt){
1276
1277 px.1 <- pdf.cdf$px.1
1278 px.2 <- pdf.cdf$px.2
1279 py.1 <- pdf.cdf$py.1
1280 py.2 <- pdf.cdf$py.2
1281
1282 if(copula.txt=="gaussian"){
1283 c1 <- gaussian.cop.den(px.1$cdf, py.1$cdf, rho1)
1284 c2 <- gaussian.cop.den(px.2$cdf, py.2$cdf, rho2)
1285 } else {
1286 if(copula.txt=="clayton"){
1287 c1 <- clayton.cop.den(px.1$cdf, py.1$cdf, rho1)
1288 c2 <- clayton.cop.den(px.2$cdf, py.2$cdf, rho2)
1289 }
1290 }
1291
1292 sum(log(p*c1*px.1$pdf*py.1$pdf + (1-p)*c2*px.2$pdf*py.2$pdf))
1293 }
1294
1295
1296
1297 # EM for 2 Gaussian, speed up computation, unfinished
1298
1299 em.2gaussian.quick <- function(x, y, p0, rho1.0, rho2.0, eps, fix.p=F, stoc=T, fix.rho2=T){
1300
1301 x <- rank(x, tie="random")
1302 y <- rank(y, tie="random")
1303
1304 # x <- rank(x, tie="average")
1305 # y <- rank(y, tie="average")
1306
1307 # nbin=20
1308 xy.min <- min(x, y)
1309 xy.max <- max(x, y)
1310 binwidth <- (xy.max-xy.min)/50
1311 breaks <- seq(xy.min-binwidth/100, xy.max+binwidth/100, by=(xy.max-xy.min+binwidth/50)/50)
1312 # breaks <- seq(xy.min, xy.max, by=binwidth)
1313
1314
1315 # initiate marginals
1316 # initialization: first p0 data has
1317 # e.z <- e.step.2gaussian(x, y, p0, rho1.0, rho2.0, x0.mar, y0.mar) # this starting point assumes two components are overlapped
1318
1319 e.z <- c(rep(0.9, round(length(x)*p0)), rep(0.1, length(x)-round(length(x)*p0)))
1320
1321 if(!stoc)
1322 para <- m.step.2gaussian.value(x, y, e.z, breaks, fix.rho2)
1323 else
1324 para <- m.step.2gaussian.stoc.value(x, y, e.z, breaks, fix.rho2)
1325
1326
1327 if(fix.p){
1328 p <- p0
1329 } else {
1330 p <- para$p
1331 }
1332
1333 if(fix.rho2){
1334 rho2 <- rho2.0
1335 } else {
1336 rho2 <- para$rho2
1337 }
1338
1339 # rho1 <- 0.8
1340 rho1 <- para$rho1
1341
1342 l0 <- loglik.2gaussian.copula.value(x, y, p, rho1, rho2, para$pdf.cdf)
1343
1344 loglik.trace <- c()
1345 loglik.trace[1] <- l0
1346 # loglik.trace[2] <- l1
1347 to.run <- T
1348
1349 i <- 2
1350
1351 # this two lines to remove
1352 # x.mar <- est.mar.hist(x, e.z, breaks)
1353 # y.mar <- est.mar.hist(y, e.z, breaks)
1354
1355 while(to.run){
1356
1357 e.z <- e.step.2gaussian.value(x, y, p, rho1, rho2, para$pdf.cdf)
1358 if(!stoc)
1359 para <- m.step.2gaussian.value(x, y, e.z, breaks, fix.rho2)
1360 else
1361 para <- m.step.2gaussian.stoc.value(x, y, e.z, breaks, fix.rho2)
1362
1363 # fix x.mar and y.mar : to remove
1364 # if(!stoc)
1365 # para <- m.step.2gaussian.value2(x, y, e.z, breaks, fix.rho2, x.mar, y.mar)
1366 # else
1367 # para <- m.step.2gaussian.stoc.value(x, y, e.z, breaks, fix.rho2)
1368
1369
1370 if(fix.p){
1371 p <- p0
1372 } else {
1373 p <- para$p
1374 }
1375
1376 if(fix.rho2){
1377 rho2 <- rho2.0
1378 } else {
1379 rho2 <- para$rho2
1380 }
1381
1382 # rho1 <- 0.8
1383 rho1 <- para$rho1
1384
1385 # l0 <- l1
1386 l1 <- loglik.2gaussian.copula.value(x, y, p, rho1, rho2, para$pdf.cdf)
1387 loglik.trace[i] <- l1
1388
1389 #cat("l1=", l1, "\n")
1390
1391 # Aitken acceleration criterion
1392 if(i > 2){
1393 l.inf <- loglik.trace[i-2] + (loglik.trace[i-1] - loglik.trace[i-2])/(1-(loglik.trace[i]-loglik.trace[i-1])/(loglik.trace[i-1]-loglik.trace[i-2]))
1394 to.run <- abs(l.inf - loglik.trace[i]) > eps
1395 #cat("para=", "p=", para$p, " rho1=", rho1, " rho2=", rho2, "\n")
1396 #cat("l.inf=", l.inf, "\n")
1397 #cat(l.inf-loglik.trace[i], "\n")
1398 }
1399
1400 i <- i+1
1401 }
1402
1403 bic <- -2*l1 + (2*(length(breaks)-1+1)+1-fix.p-fix.rho2)*log(length(x)) # parameters
1404 return(list(para=list(p=para$p, rho1=rho1, rho2=rho2),
1405 loglik=l1, bic=bic, e.z=e.z, conf.z = para$conf.z,
1406 loglik.trace=loglik.trace, x.mar=para$x.mar, y.mar=para$y.mar,
1407 breaks=breaks))
1408 }
1409
1410
1411
1412 em.2copula.quick <- function(x, y, p0, rho1.0, rho2.0, eps, fix.p=F, stoc=T, fix.rho2=T, copula.txt, nbin=50){
1413
1414 x <- rank(x, tie="random")
1415 y <- rank(y, tie="random")
1416
1417 # x <- rank(x, tie="first")
1418 # y <- rank(y, tie="first")
1419
1420 # nbin=50
1421 xy.min <- min(x, y)
1422 xy.max <- max(x, y)
1423 binwidth <- (xy.max-xy.min)/50
1424 breaks <- seq(xy.min-binwidth/100, xy.max+binwidth/100, by=(xy.max-xy.min+binwidth/50)/nbin)
1425 # breaks <- seq(xy.min, xy.max, by=binwidth)
1426
1427 # initiate marginals
1428 # initialization: first p0 data has
1429 # e.z <- e.step.2gaussian(x, y, p0, rho1.0, rho2.0, x0.mar, y0.mar) # this starting point assumes two components are overlapped
1430
1431 e.z <- c(rep(0.9, round(length(x)*p0)), rep(0.1, length(x)-round(length(x)*p0)))
1432
1433
1434 if(!stoc)
1435 para <- m.step.2copula.value(x, y, e.z, breaks, fix.rho2, copula.txt)
1436 else
1437 para <- m.step.2copula.stoc.value(x, y, e.z, breaks, fix.rho2, copula.txt)
1438
1439 if(fix.p){
1440 p <- p0
1441 } else {
1442 p <- para$p
1443 }
1444
1445 if(fix.rho2){
1446 rho2 <- rho2.0
1447 } else {
1448 rho2 <- para$rho2
1449 }
1450
1451 l0 <- loglik.2copula.value(x, y, p, para$rho1, rho2, para$pdf.cdf, copula.txt)
1452
1453 loglik.trace <- c()
1454 loglik.trace[1] <- l0
1455 # loglik.trace[2] <- l1
1456 to.run <- T
1457
1458 i <- 2
1459
1460 while(to.run){
1461
1462 e.z <- e.step.2copula.value(x, y, p, para$rho1, rho2, para$pdf.cdf, copula.txt)
1463 if(!stoc)
1464 para <- m.step.2copula.value(x, y, e.z, breaks, fix.rho2, copula.txt)
1465 else
1466 para <- m.step.2copula.stoc.value(x, y, e.z, breaks, fix.rho2, copula.txt)
1467
1468 if(fix.p){
1469 p <- p0
1470 } else {
1471 p <- para$p
1472 }
1473
1474 if(fix.rho2){
1475 rho2 <- rho2.0
1476 } else {
1477 rho2 <- para$rho2
1478 }
1479
1480
1481 # l0 <- l1
1482 l1 <- loglik.2copula.value(x, y, p, para$rho1, rho2, para$pdf.cdf, copula.txt)
1483 loglik.trace[i] <- l1
1484
1485 cat("l1=", l1, "\n")
1486
1487 # Aitken acceleration criterion
1488 if(i > 2){
1489 l.inf <- loglik.trace[i-2] + (loglik.trace[i-1] - loglik.trace[i-2])/(1-(loglik.trace[i]-loglik.trace[i-1])/(loglik.trace[i-1]-loglik.trace[i-2]))
1490 to.run <- abs(l.inf - loglik.trace[i]) > eps
1491 cat("para=", "p=", para$p, " rho1=", para$rho1, " rho2=", rho2, "\n")
1492 #cat("l.inf=", l.inf, "\n")
1493 #cat(l.inf-loglik.trace[i], "\n")
1494 }
1495
1496 i <- i+1
1497 }
1498
1499 bic <- -2*l1 + (2*(length(breaks)-1+1)+1-fix.p-fix.rho2)*log(length(x)) # parameters
1500 return(list(para=list(p=para$p, rho1=para$rho1, rho2=rho2),
1501 loglik=l1, bic=bic, e.z=e.z, conf.z = para$conf.z,
1502 loglik.trace=loglik.trace, x.mar=para$x.mar, y.mar=para$y.mar,
1503 breaks=breaks))
1504 }
1505
1506
1507 #######################
1508 ####################### fit EM procedure for the matched peaks
1509 #######################
1510
1511 # remove the unmatched ones
1512 #rm.unmatch <- function(sample1, sample2, p.value.impute=0){
1513 #
1514 # sample1.prune <- sample1[sample1$sig.value > p.value.impute & sample2$sig.value > p.value.impute,]
1515 # sample2.prune <- sample2[sample1$sig.value > p.value.impute & sample2$sig.value > p.value.impute,]
1516 #
1517 # invisible(list(sample1=sample1.prune$sig.value, sample2=sample2.prune$sig.value))
1518 #}
1519
1520
1521 # fit 2-component model
1522 #fit.em <- function(sample12, fix.rho2=T){
1523 #
1524 # prune.sample <- rm.unmatch(sample12$merge1, sample12$merge2)
1525 #
1526 # em.fit <- em.2gaussian.quick(-prune.sample$sample1, -prune.sample$sample2,
1527 # p0=0.5, rho1.0=0.7, rho2.0=0, eps=0.01, fix.p=F, stoc=F, fix.rho2)
1528 #
1529 # invisible(list(em.fit=em.fit, data.pruned=prune.sample))
1530 #}
1531
1532
1533 rm.unmatch <- function(sample1, sample2, p.value.impute=0){
1534
1535 sample1.prune <- sample1[sample1$sig.value > p.value.impute & sample2$sig.value > p.value.impute,]
1536 sample2.prune <- sample2[sample1$sig.value > p.value.impute & sample2$sig.value > p.value.impute,]
1537
1538 invisible(list(sample1=sample1.prune, sample2=sample2.prune))
1539 }
1540
1541
1542 # fit 2-component model
1543 fit.em <- function(sample12, fix.rho2=T){
1544
1545 prune.sample <- rm.unmatch(sample12$merge1, sample12$merge2)
1546
1547 em.fit <- em.2gaussian.quick(-prune.sample$sample1$sig.value, -prune.sample$sample2$sig.value,
1548 p0=0.5, rho1.0=0.7, rho2.0=0, eps=0.01, fix.p=F, stoc=F, fix.rho2)
1549
1550 invisible(list(em.fit=em.fit, data.pruned=prune.sample))
1551 }
1552
1553
1554
1555 fit.2copula.em <- function(sample12, fix.rho2=T, copula.txt){
1556
1557 prune.sample <- rm.unmatch(sample12$merge1, sample12$merge2)
1558
1559 # o <- order(prune.sample$sample1)
1560 # n <- length(prune.sample$sample1)
1561
1562 # para <- init(prune.sample$sample1$sig.value, prune.sample$sample2$sig.value, c(rep(0, round(n/3)), rep(c(0,1), round(n/6)), rep(1, n-round(n/3)-round(n/6))))
1563
1564 # temp <- init.dist(f0, f1)
1565 para <- list()
1566 para$rho <- 0.6
1567 para$p <- 0.3
1568 para$mu <- 2.5
1569 para$sigma <- 1
1570 ## para$mu <- -temp$mu
1571 ## para$sigma <- temp$sigma
1572 #cat("mu=", para$mu, "sigma=", para$sigma, "\n")
1573
1574 # em.fit <- em.transform.1loop(-prune.sample$sample1, -prune.sample$sample2,
1575 cat("EM is running")
1576 em.fit <- em.transform(prune.sample$sample1$sig.value, prune.sample$sample2$sig.value, para$mu, para$sigma, para$rho, para$p, eps=0.01)
1577
1578 invisible(list(em.fit=em.fit, data.pruned=prune.sample))
1579 }
1580
1581
1582
1583
1584 # fit 1-component model
1585 fit.1.component <- function(data.pruned, breaks){
1586
1587 # gaussian.1 <- fit.gaussian.1(-data.pruned$sample1$sig.value, -data.pruned$sample2$sig.value, breaks)
1588 # clayton.1 <- fit.clayton.1(-data.pruned$sample1$sig.value, -data.pruned$sample2$sig.value, breaks)
1589
1590 gaussian.1 <- fit.gaussian.1(-data.pruned$sample1, -data.pruned$sample2, breaks)
1591 clayton.1 <- fit.clayton.1(-data.pruned$sample1, -data.pruned$sample2, breaks)
1592
1593 return(list(gaussian.1=gaussian.1, clayton.1=clayton.1))
1594 }
1595
1596
1597
1598 #################
1599 # Fit a single component
1600 #################
1601
1602 # a single gaussian copula
1603 # if breaks=NULL, use empirical pdf, otherwise use histogram estimate
1604 fit.gaussian.1 <- function(x, y, breaks=NULL){
1605
1606 # rank transformed and compute the empirical cdf
1607 t <- emp.mar.cdf.rank(x)
1608 s <- emp.mar.cdf.rank(y)
1609
1610 mle.rho <- mle.gaussian.copula(t, s, rep(1, length(t)))
1611
1612 c1 <- gaussian.cop.den(t, s, mle.rho)
1613 cat("c1", sum(log(c1)), "\n")
1614
1615 if(is.null(breaks)){
1616 f1 <- emp.mar.pdf.rank(t)
1617 f2 <- emp.mar.pdf.rank(s)
1618 } else {
1619 x.mar <- est.mar.hist(rank(x), rep(1, length(x)), breaks)
1620 y.mar <- est.mar.hist(rank(y), rep(1, length(y)), breaks)
1621
1622 f1 <- x.mar$f1.value$pdf # only one component
1623 f2 <- y.mar$f1.value$pdf
1624 }
1625
1626
1627 cat("f1", sum(log(f1)), "\n")
1628 cat("f2", sum(log(f2)), "\n")
1629
1630 loglik <- sum(log(c1)+log(f1)+log(f2))
1631
1632 bic <- -2*loglik + log(length(t))*(1+length(breaks)-1)
1633
1634 return(list(rho=mle.rho, loglik=loglik, bic=bic))
1635 }
1636
1637
1638 # a single Clayton copula
1639 fit.clayton.1 <- function(x, y, breaks=NULL){
1640
1641 # rank transformed and compute the empirical cdf
1642 t <- emp.mar.cdf.rank(x)
1643 s <- emp.mar.cdf.rank(y)
1644
1645 mle.rho <- mle.clayton.copula(t, s, rep(1, length(t)))
1646
1647 c1 <- clayton.cop.den(t, s, mle.rho)
1648
1649 if(is.null(breaks)){
1650 f1 <- emp.mar.pdf.rank(t)
1651 f2 <- emp.mar.pdf.rank(s)
1652 } else {
1653 x.mar <- est.mar.hist(rank(x), rep(1, length(x)), breaks)
1654 y.mar <- est.mar.hist(rank(y), rep(1, length(y)), breaks)
1655
1656 f1 <- x.mar$f1.value$pdf # only one component
1657 f2 <- y.mar$f1.value$pdf
1658 }
1659
1660 loglik <- sum(log(c1)+log(f1)+log(f2))
1661
1662 bic <- -2*loglik + log(length(t))*(1+length(breaks)-1)
1663
1664 return(list(rho=mle.rho, tau=rho/(rho+2), loglik=loglik, bic=bic))
1665 }
1666
1667 ## obsolete function (01-06-2010)
1668 ## compute the average posterior probability to belong to the random component
1669 ## for peaks selected at different cutoffs
1670 comp.uri.ez <- function(tt, u, v, e.z){
1671
1672 u.t <- quantile(u, prob=(1-tt))
1673 v.t <- quantile(v, prob=(1-tt))
1674
1675 # ez <- mean(e.z[u >= u.t & v >=u.t]) Is this wrong?
1676 ez <- mean(e.z[u >= u.t & v >=v.t])
1677
1678 return(ez)
1679 }
1680
1681 ## obsolete function (01-06-2010)
1682 # compute the largest posterior error probability corresponding to
1683 # the square centered at the origin and spanned top tt% on both coordinates
1684 # so the consistent low rank ones are excluded
1685 # boundary.txt: either "max" or "min", if it is error prob, use "max"
1686 comp.ez.cutoff <- function(tt, u, v, e.z, boundary.txt){
1687
1688 u.t <- quantile(u, prob=(1-tt))
1689 v.t <- quantile(v, prob=(1-tt))
1690
1691 if(boundary.txt == "max"){
1692 # ez.bound <- max(e.z[u >= u.t & v >=u.t])
1693 ez.bound <- max(e.z[u >= u.t & v >=v.t])
1694 } else {
1695 # ez.bound <- min(e.z[u >= u.t & v >=u.t])
1696 ez.bound <- min(e.z[u >= u.t & v >=v.t])
1697 }
1698
1699 return(ez.bound)
1700
1701 }
1702
1703 # obsolete functions: 01-06-2010
1704 # compute the error rate
1705 # u.t and v.t are the quantiles
1706 # this one is used for the plots generated initially in the brief writeup
1707 # and it was used for processing merged data in July before the IDR definition
1708 # is formalized
1709 # It does not implement the current definition of IDR
1710 get.ez.tt.old <- function(em.fit, reverse=T, fdr.level=c(0.01, 0.05, 0.1)){
1711
1712 u <- em.fit$data.pruned$sample1
1713 v <- em.fit$data.pruned$sample2
1714
1715 tt <- seq(0.01, 0.99, by=0.01)
1716 if(reverse){
1717 e.z <- 1-em.fit$em.fit$e.z # this is the error prob
1718 uri.ez <- sapply(tt, comp.uri.ez, u=u, v=v, e.z=e.z)
1719 ez.bound <- sapply(tt, comp.ez.cutoff, u=u, v=v, e.z=e.z, boundary.txt="max")
1720 } else {
1721 e.z <- em.fit$em.fit$e.z
1722 uri.ez <- sapply(tt, comp.uri.ez, u=u, v=v, e.z=e.z)
1723 ez.bound <- sapply(tt, comp.ez.cutoff, u=u, v=v, e.z=e.z, boundary.txt="min")
1724 }
1725
1726 u.t <- quantile(u, prob=(1-tt))
1727 v.t <- quantile(v, prob=(1-tt))
1728
1729 # find the levels on the two replicates
1730 sig.value1 <- c()
1731 sig.value2 <- c()
1732 error.prob.cutoff <- c()
1733 n.selected.match <- c()
1734
1735 for(i in 1:length(fdr.level)){
1736
1737 # find which uri.ez is closet to fdr.level
1738 index <- which.min(abs(uri.ez - fdr.level[i]))
1739 sig.value1[i] <- u.t[index]
1740 sig.value2[i] <- v.t[index]
1741 error.prob.cutoff[i] <- ez.bound[index]
1742 if(reverse){
1743 n.selected.match[i] <- sum(e.z<=ez.bound[index])
1744 } else {
1745 n.selected.match[i] <- sum(e.z>=ez.bound[index])
1746 }
1747 }
1748
1749 # output the cutoff of posterior probability, signal values on two replicates
1750 map.uv <- cbind(error.prob.cutoff, sig.value1, sig.value2, n.selected.match)
1751
1752 return(list(n=tt*length(u), uri.ez=uri.ez, u.t=u.t, v.t=v.t, tt=tt, fdr.level=fdr.level, map.uv=map.uv, e.z=e.z, error.prob.cutoff=error.prob.cutoff))
1753 }
1754
1755 # created: 01-06-2010
1756 # Output IDR at various number of selected peaks
1757 # Find cutoff (idr cutoff, sig.value cutoff on each replicate) for specified IDR level
1758 # IDR definition is similar to FDR
1759 get.ez.tt <- function(em.fit, idr.level=c(0.01, 0.05, 0.1)){
1760
1761 # u <- em.fit$data.pruned$sample1$sig.value
1762 # v <- em.fit$data.pruned$sample2$sig.value
1763 u <- em.fit$data.pruned$sample1
1764 v <- em.fit$data.pruned$sample2
1765
1766 e.z <- 1-em.fit$em.fit$e.z # this is the error prob
1767
1768 o <- order(e.z)
1769 e.z.ordered <- e.z[o]
1770 n.select <- c(1:length(e.z))
1771 IDR <- cumsum(e.z.ordered)/n.select
1772
1773 u.o <- u[o]
1774 v.o <- v[o]
1775
1776 n.level <- length(idr.level)
1777 # sig.value1 <- rep(NA, n.level)
1778 # sig.value2 <- rep(NA, n.level)
1779 ez.cutoff <- rep(NA, n.level)
1780 n.selected <- rep(NA, n.level)
1781
1782 for(i in 1:length(idr.level)){
1783
1784 # find which uri.ez is closet to fdr.level
1785 index <- which.min(abs(IDR - idr.level[i]))
1786 # sig.value1[i] <- min(u.o[1:index])
1787 # sig.value2[i] <- min(v.o[1:index])
1788 ez.cutoff[i] <- e.z[index]
1789 n.selected[i] <- sum(e.z<=ez.cutoff[i])
1790 }
1791
1792 # output the cutoff of posterior probability, number of selected overlapped peaks
1793 # map.uv <- cbind(ez.cutoff, sig.value1, sig.value2, n.selected)
1794
1795 map.uv <- cbind(ez.cutoff, n.selected)
1796
1797 return(list(n=n.select, IDR=IDR, idr.level=idr.level, map.uv=map.uv))
1798 }
1799
1800 # return(list(n=tt*length(u), uri.ez=uri.ez, fdr.level=fdr.level, map.uv=map.uv, e.z=e.z, error.prob.cutoff=error.prob.cutoff))
1801
1802
1803
1804
1805
1806 ### compute the mean of the marginals
1807 get.mar.mean <- function(em.out){
1808
1809 x.f1 <- em.out$x.mar$f1
1810 x.f2 <- em.out$x.mar$f2
1811
1812 y.f1 <- em.out$y.mar$f1
1813 y.f2 <- em.out$y.mar$f2
1814
1815 x.stat1 <- get.hist.mean(x.f1)
1816 x.stat2 <- get.hist.mean(x.f2)
1817 y.stat1 <- get.hist.mean(y.f1)
1818 y.stat2 <- get.hist.mean(y.f2)
1819
1820 return(list(x.mean1=x.stat1$mean, x.mean2=x.stat2$mean,
1821 y.mean1=y.stat1$mean, y.mean2=y.stat2$mean,
1822 x.sd1=x.stat1$sd, x.sd2=x.stat2$sd,
1823 y.sd1=y.stat1$sd, y.sd2=y.stat2$sd
1824 ))
1825
1826 }
1827
1828
1829 # compute the mean of marginals
1830 get.hist.mean <- function(x.f){
1831
1832 nbreaks <- length(x.f$breaks)
1833 x.bin <- x.f$breaks[-1]-x.f$breaks[-nbreaks]
1834
1835 x.mid <- (x.f$breaks[-nbreaks]+x.f$breaks[-1])/2
1836 x.mean <- sum(x.mid*x.f$density*x.bin)
1837 x.sd <- sqrt(sum(x.mid*x.mid*x.f$density*x.bin)-x.mean^2)
1838
1839 return(list(mean=x.mean, sd=x.sd))
1840 }
1841
1842 get.hist.var <- function(x.f){
1843
1844 nbreaks <- length(x.f$breaks)
1845 x.bin <- x.f$breaks[-1]-x.f$breaks[-nbreaks]
1846
1847 x.mid <- (x.f$breaks[-nbreaks]+x.f$breaks[-1])/2
1848 x.mean <- sum(x.mid*x.f$density*x.bin)
1849
1850 return(mean=x.mean)
1851 }
1852
1853 # obsolete function (01-06-2010)
1854 # plot
1855 plot.ez.group.old <- function(ez.list, plot.dir, file.name=NULL, legend.txt, y.lim=NULL, xlab.txt="num of significant peaks", ylab.txt="avg posterior prob of being random", col.txt=NULL, title.txt=NULL){
1856
1857 if(is.null(col.txt))
1858 col.txt <- c("black", "red", "purple", "green", "blue", "cyan", "magenta", "orange", "grey")
1859
1860 x <- c()
1861 y <- c()
1862
1863 for(i in 1:length(ez.list)){
1864 x <- c(x, ez.list[[i]]$n)
1865
1866 y <- c(y, ez.list[[i]]$uri.ez)
1867 }
1868
1869 if(is.null(y.lim))
1870 y.lim <- c(0, max(y))
1871
1872 if(!is.null(file.name)){
1873 postscript(paste(plot.dir, "ez.", file.name, sep=""))
1874 par(mfrow=c(1,1), mar=c(5,5,4,2))
1875 }
1876
1877 plot(x, y, ylim=y.lim, type="n", xlab=xlab.txt, ylab=ylab.txt, lwd=5, cex=5, cex.axis=2, cex.lab=2)
1878
1879 for(i in 1:length(ez.list)){
1880 lines(ez.list[[i]]$n, ez.list[[i]]$uri.ez, col=col.txt[i], cex=2, lwd=5)
1881 }
1882
1883 # plot(ez.list[[1]]$u.t, y, ylim=y.lim, type="l", xlab="rep-sig", ylab=ylab.txt, lwd=5, cex=5, cex.axis=2, cex.lab=2)
1884 # plot(ez.list[[1]]$v.t, y, ylim=y.lim, type="l", xlab="rep-sig", ylab=ylab.txt, lwd=5, cex=5, cex.axis=2, cex.lab=2)
1885
1886
1887 legend(0, y.lim[2], legend=legend.txt, col=col.txt[1:length(col.txt)], lty=1, lwd=5, cex=2)
1888
1889 if(!is.null(title))
1890 title(title.txt)
1891
1892 if(!is.null(file.name)){
1893 dev.off()
1894 }
1895
1896 }
1897
1898
1899 plot.ez.group <- function(ez.list, plot.dir, file.name=NULL, legend.txt, y.lim=NULL, xlab.txt="num of significant peaks", ylab.txt="IDR", col.txt=NULL, title.txt=NULL){
1900
1901 if(is.null(col.txt))
1902 col.txt <- c("black", "red", "purple", "green", "blue", "cyan", "magenta", "orange", "grey")
1903
1904 n.entry <- length(ez.list)
1905 x <- rep(NA, n.entry)
1906 y.max <- rep(NA, n.entry)
1907
1908 for(i in 1:n.entry){
1909 x[i] <- max(ez.list[[i]]$n)
1910
1911 y.max[i] <- max(ez.list[[i]]$IDR)
1912
1913 }
1914
1915 if(is.null(y.lim))
1916 y.lim <- c(0, max(y.max))
1917
1918 if(!is.null(file.name)){
1919 postscript(paste(plot.dir, "ez.", file.name, sep=""))
1920 par(mfrow=c(1,1), mar=c(5,5,4,2))
1921 }
1922
1923
1924
1925 plot(c(0, max(x)), y.lim, ylim=y.lim, type="n", xlab=xlab.txt, ylab=ylab.txt, lwd=5, cex=5, cex.axis=2, cex.lab=2)
1926
1927 q <- seq(0.01, 0.99, by=0.01)
1928
1929 for(i in 1:length(ez.list)){
1930
1931 n.plot <- round(quantile(ez.list[[i]]$n, prob=q))
1932 IDR.plot <- ez.list[[i]]$IDR[n.plot]
1933 lines(n.plot, IDR.plot, col=col.txt[i], cex=2, lwd=5)
1934 }
1935
1936
1937 legend(0, y.lim[2], legend=legend.txt, col=col.txt[1:length(col.txt)], lty=1, lwd=5, cex=2)
1938
1939 if(!is.null(title))
1940 title(title.txt)
1941
1942 if(!is.null(file.name)){
1943 dev.off()
1944 }
1945
1946 }
1947
1948
1949
1950 #############################################################################
1951 #############################################################################
1952 # statistics about peaks selected on the individual replicates
1953 #
1954 # idr.level: the consistency cutoff, say 0.05
1955 # uri.output: a list of uri.output from consistency analysis generated by batch-consistency-analysis.r
1956 # ez.list : a list of IDRs computed from get.ez.tt using the same idr.level
1957 #
1958 ##################
1959
1960
1961 # obsolete?
1962 # compute the error rate
1963 # u.t and v.t are the quantiles
1964 #
1965 # map back to all peaks and report the number of peaks selected
1966 get.ez.tt.all.old <- function(em.fit, all.data1, all.data2, idr.level){
1967
1968 u <- em.fit$data.pruned$sample1
1969 v <- em.fit$data.pruned$sample2
1970
1971 tt <- seq(0.01, 0.99, by=0.01)
1972 # if(reverse){
1973 e.z <- 1-em.fit$em.fit$e.z # this is the error prob
1974 uri.ez <- sapply(tt, comp.uri.ez, u=u, v=v, e.z=e.z)
1975 ez.bound <- sapply(tt, comp.ez.cutoff, u=u, v=v, e.z=e.z, boundary.txt="max")
1976 # } else {
1977 # e.z <- em.fit$em.fit$e.z
1978 # uri.ez <- sapply(tt, comp.uri.ez, u=u, v=v, e.z=e.z)
1979 # ez.bound <- sapply(tt, comp.ez.cutoff, u=u, v=v, e.z=e.z, boundary.txt="min")
1980 # }
1981
1982 u.t <- quantile(u, prob=(1-tt))
1983 v.t <- quantile(v, prob=(1-tt))
1984
1985 # find the levels on the two replicates
1986 sig.value1 <- c()
1987 sig.value2 <- c()
1988 error.prob.cutoff <- c()
1989 n.selected.match <- c()
1990 npeak.rep1 <- c()
1991 npeak.rep2 <- c()
1992
1993 for(i in 1:length(idr.level)){
1994
1995 # find which uri.ez is closet to idr.level
1996 index <- which.min(abs(uri.ez - as.numeric(idr.level[i])))
1997
1998 sig.value1[i] <- u.t[index]
1999 sig.value2[i] <- v.t[index]
2000 error.prob.cutoff[i] <- ez.bound[index]
2001 n.selected.match[i] <- sum(u>= u.t[index] & v>=v.t[index])
2002
2003 npeak.rep1[i] <- sum(all.data1["sig.value"] >= sig.value1[i])
2004 npeak.rep2[i] <- sum(all.data2["sig.value"] >= sig.value2[i])
2005 }
2006
2007
2008 # output the cutoff of posterior probability, signal values on two replicates
2009 map.uv <- cbind(error.prob.cutoff, sig.value1, sig.value2, n.selected.match, npeak.rep1, npeak.rep2)
2010
2011 return(list(n=tt*length(u), uri.ez=uri.ez, u.t=u.t, v.t=v.t, tt=tt, idr.level=idr.level, map.uv=map.uv, e.z=e.z, error.prob.cutoff=error.prob.cutoff))
2012 }
2013
2014
2015 get.ez.tt.all <- function(em.fit, all.data1, all.data2, idr.level=c(0.01, 0.05, 0.1)){
2016
2017 u <- em.fit$data.pruned$sample1$sig.value
2018 v <- em.fit$data.pruned$sample2$sig.value
2019 # u <- em.fit$data.pruned$sample1
2020 # v <- em.fit$data.pruned$sample2
2021
2022 e.z <- 1-em.fit$em.fit$e.z # this is the error prob
2023
2024 o <- order(e.z)
2025 e.z.ordered <- e.z[o]
2026 n.select <- c(1:length(e.z))
2027 IDR <- cumsum(e.z.ordered)/n.select
2028
2029 u.o <- u[o]
2030 v.o <- v[o]
2031
2032 n.level <- length(idr.level)
2033 # sig.value1 <- rep(NA, n.level)
2034 # sig.value2 <- rep(NA, n.level)
2035 ez.cutoff <- rep(NA, n.level)
2036 n.selected <- rep(NA, n.level)
2037 npeak.rep1 <- rep(NA, n.level)
2038 npeak.rep2 <- rep(NA, n.level)
2039
2040 for(i in 1:length(idr.level)){
2041
2042 # find which uri.ez is closet to fdr.level
2043 index <- which.min(abs(IDR - idr.level[i]))
2044 # sig.value1[i] <- min(u.o[1:index])
2045 # sig.value2[i] <- min(v.o[1:index])
2046 ez.cutoff[i] <- e.z.ordered[index] # fixed on 02/20/10
2047 n.selected[i] <- sum(e.z<=ez.cutoff[i])
2048 # npeak.rep1[i] <- sum(all.data1["sig.value"] >= sig.value1[i])
2049 # npeak.rep2[i] <- sum(all.data2["sig.value"] >= sig.value2[i])
2050 }
2051
2052 # output the cutoff of posterior probability, number of selected overlapped peaks
2053 map.uv <- cbind(ez.cutoff, n.selected)
2054
2055 return(list(n=n.select, IDR=IDR, idr.level=idr.level, map.uv=map.uv))
2056 }
2057
2058 # return(list(n=tt*length(u), uri.ez=uri.ez, fdr.level=fdr.level, map.uv=map.uv, e.z=e.z, error.prob.cutoff=error.prob.cutoff))
2059
2060
2061
2062
2063
2064
2065 ####### the following is for determining thresholds for merged dataset
2066
2067 ############# select peaks above a given threshold
2068 #
2069 # pass.threshold: a simple method, passing the threshold on the threshold on the individual replicate to the pooled sample
2070 #
2071 # sig.map.list: a list of matrix to include all the cutoff values, each row corresponds to a cutoff. The first column is idr.level
2072 # the 2nd column is the cutoff of ez, the rest of columns are consistency analysis for other replicates
2073 # sig.value.name: the name of the sig.value column
2074 # combined: combined dataset
2075 # nrep: number of pairs of comparisons
2076 #
2077 # Procedure:
2078 # 1. Find the significant threshold corresponding to the idr cutoff on the matched peaks.
2079 # 2. Each time we will get two or more (if >2 replicates) cutoffs and will report the most stringent and the least stringent
2080 # cutoff and the number of peaks selected at those two cutoffs
2081 #############
2082
2083 pass.threshold <- function(sig.map.list, sig.value.name, combined, idr.level, nrep, chr.size){
2084
2085 sig.map <- c()
2086
2087 # choose idr.level
2088 idr.index <- which(rbind(sig.map.list[[1]])[,1] == idr.level)
2089 if(length(i) ==0){
2090 print("no level matches specified idr.level")
2091 return(-1)
2092 }
2093
2094 for(i in 1:length(sig.map.list))
2095 sig.map <- c(sig.map, rbind(sig.map.list[[i]])[idr.index, c("sig.value1", "sig.value2")])
2096
2097
2098 npeak.tight <- c()
2099 npeak.loose <- c()
2100
2101
2102 max.sig <- max(sig.map)
2103 min.sig <- min(sig.map)
2104 selected.sig.tight <- combined[combined[,sig.value.name]>=max.sig, ]
2105 selected.sig.loose <- combined[combined[,sig.value.name]>=min.sig, ]
2106
2107 selected.sig.tight <- deconcatenate.chr(selected.sig.tight, chr.size)[,c("chr", "start", "stop", "signal.value", "p.value", "q.value")]
2108 selected.sig.loose <- deconcatenate.chr(selected.sig.loose, chr.size)[,c("chr", "start", "stop", "signal.value", "p.value", "q.value")]
2109
2110 npeak.tight <- nrow(selected.sig.tight)
2111 npeak.loose <- nrow(selected.sig.loose)
2112
2113
2114 npeak.stat <- list(idr.level=idr.level, max.sig=max.sig, min.sig=min.sig, npeak.tight=npeak.tight, npeak.loose=npeak.loose)
2115
2116 invisible(list(npeak.stat=npeak.stat, combined.selected.tight=selected.sig.tight, combined.selected.loose=selected.sig.loose))
2117 }
2118
2119 #################
2120 # pass the regions selected from consistency analysis to combined data
2121 # Threshold is determined on the replicates, the regions above the threshold are selected
2122 # then peaks on the combined data are selected from the selected regions
2123 #
2124 # To avoid being too stringent, regions satisfying the following conditions are selected
2125 # 1. regions above the significant threshold determined by consistency analysis on either replicate
2126 # 2. regions that have consistent low peaks, i.e. posterior prob > threshold but not passing the significant threshold
2127 #
2128 # This method doesn't make a difference when using different thresholds
2129 #################
2130
2131 pass.region <- function(sig.map.list, uri.output, ez.list, em.output, combined, idr.level, sig.value.impute=0, chr.size){
2132
2133 combined <- combined[, c("start", "stop", "sig.value", "signal.value", "p.value", "q.value")]
2134 npair <- length(uri.output) # number of pairs of consistency analysis
2135 combined.region <- c()
2136
2137 # choose idr.level
2138 idr.index <- which(rbind(sig.map.list[[1]])[,1] == idr.level)
2139 if(length(idr.index) ==0){
2140 print("no level matches specified idr.level")
2141 return(-1)
2142 }
2143
2144 for(j in 1:npair){
2145 # select peaks from individual replicates using individual cutoff
2146 above.1 <- uri.output[[j]]$data12.enrich$merge1["sig.value"] >= ez.list[[j]]$map.uv[idr.index,"sig.value1"]
2147 above.2 <- uri.output[[j]]$data12.enrich$merge1["sig.value"] >= ez.list[[j]]$map.uv[idr.index,"sig.value2"]
2148 selected.sig.rep1 <- uri.output[[j]]$data12.enrich$merge1[above.1, c("start", "stop", "sig.value", "signal.value", "p.value", "q.value")]
2149 selected.sig.rep2 <- uri.output[[j]]$data12.enrich$merge2[above.2, c("start", "stop", "sig.value", "signal.value", "p.value", "q.value")]
2150
2151 # find the peaks that are overlapped with reliable peaks in the individual replicates
2152 overlap.1 <- pair.peaks(selected.sig.rep1, combined)$merge2
2153 overlap.2 <- pair.peaks(selected.sig.rep2, combined)$merge2
2154
2155 # choose the ones with significant value > 0, which are the overlapped ones
2156
2157 combined.in1 <- overlap.1[overlap.1$sig.value > sig.value.impute, c("start", "stop", "sig.value", "signal.value", "p.value", "q.value")]
2158 combined.in2 <- overlap.2[overlap.2$sig.value > sig.value.impute, c("start", "stop", "sig.value", "signal.value", "p.value", "q.value")]
2159
2160 ## consistent low significant ones
2161 ## first find consistenct ones, ie. high posterior prob
2162 # is.consistent <- ez.list[[j]]$e.z < ez.list[[j]]$ez.cutoff
2163
2164 # data.matched <- keep.match(uri.output[[j]]$data12.enrich$merge1[!above.1, ], uri.output[[j]]$data12.enrich$merge2[!above.2, ], sig.value.impute=0)
2165 # data.matched$sample1 <- data.matched$sample1[, c("start", "stop", "sig.value", "signal.value", "p.value", "q.value")]
2166 # data.matched$sample2 <- data.matched$sample2[, c("start", "stop", "sig.value", "signal.value", "p.value", "q.value")]
2167
2168 # consistent.in1 <- data.matched$sample1[is.consistent, ]
2169 # consistent.in2 <- data.matched$sample2[is.consistent, ]
2170
2171 # overlap.consistent.1 <- pair.peaks(consistent.in1, combined)$merge2
2172 # overlap.consistent.2 <- pair.peaks(consistent.in2, combined)$merge2
2173
2174 ## choose the ones with significant value > 0, which are the overlapped ones
2175
2176 # combined.consistent.in1 <- overlap.consistent.1[overlap.consistent.1$sig.value > sig.value.impute, ]
2177 # combined.consistent.in2 <- overlap.consistent.2[overlap.consistent.2$sig.value > sig.value.impute, ]
2178
2179 # combined.region <- rbind(combined.region, combined.in1, combined.in2, combined.consistent.in1, combined.consistent.in2)
2180
2181 combined.region <- rbind(combined.region, combined.in1, combined.in2)
2182
2183 is.repeated <- duplicated(combined.region$start)
2184 combined.region <- combined.region[!is.repeated, c("start", "stop", "sig.value", "signal.value", "p.value", "q.value")]
2185
2186 }
2187 npeak <- nrow(combined.region)
2188
2189 sig.combined <- c(min(combined.region[,"sig.value"], na.rm=T), max(combined.region[,"sig.value"], na.rm=T))
2190
2191 # idr.combined <- c(min(combined.region[,"q.value"], na.rm=T), max(combined.region[,"q.value"], na.rm=T))
2192
2193 npeak.stat <- list(idr.level=idr.level, npeak=npeak)
2194
2195 combined.region <- deconcatenate.chr(combined.region, chr.size)[,c("chr", "start", "stop", "signal.value", "p.value", "q.value")]
2196
2197 invisible(list(npeak.stat=npeak.stat, combined.selected=combined.region, sig.combined=sig.combined))
2198 }
2199
2200 ################
2201 # pass structure: this method does another round of inference on the combined data
2202 #
2203 # To make the mixture structure comparable on the replicates and the combined data, the 2nd inference is done on the peaks
2204 # at the reliable regions on the combined data, using rank transformed significant values. The mixture structure is estimated using my consistency analysis, which
2205 # estimates marginal distributions of ranks using nonparametric ways. Then the significant values are found out.
2206 # There are several advantages to do it this way:
2207 # 1. The premise of passing structure is that the means and variance (i.e. distribution) of two replicates should be the same
2208 # The significant values on the two replicates clearly have different distributions. The structure estimated from consistency
2209 # analysis will generate similar rank distribution on two replicates by its setup (i.e. same number of peaks are paired up).
2210 # 2. Because pooled sample is a black box, the structure is more likely to be followed in the matched regions than other locations,
2211 # after all, we don't know what other things are. If even the structure doesn't hold on the matched regions,
2212 # which is possible, let alone the other regions. Focusing on the reliable regions helps to get rid of those unknown noises.
2213 #
2214 #
2215 # modified on 2-20-10: reverse rank.combined, make big sig.value with small
2216 # ranks, to be consistent with f1 and f2
2217 ################
2218
2219 pass.structure <- function(uri.output, em.output, combined, idr.level, sig.value.impute, chr.size, overlap.ratio=0){
2220
2221 columns.keep <- c("sig.value", "start", "stop", "signal.value", "p.value", "q.value", "chr", "start.ori", "stop.ori")
2222 combined <- combined[, columns.keep]
2223 combined.selected.all <- c()
2224
2225 for(j in 1:npair){
2226
2227 sample1 <- uri.output[[j]]$data12.enrich$merge1[, columns.keep]
2228 sample2 <- uri.output[[j]]$data12.enrich$merge2[, columns.keep]
2229
2230 # find peaks on the matched region on the combined one
2231 data.matched <- keep.match(sample1, sample2, sig.value.impute=sig.value.impute)
2232
2233 data.matched$sample1 <- data.matched$sample1[, columns.keep]
2234 data.matched$sample2 <- data.matched$sample2[, columns.keep]
2235
2236 overlap.1 <- pair.peaks.filter(data.matched$sample1, combined, p.value.impute=sig.value.impute, overlap.ratio)$merge2
2237 overlap.2 <- pair.peaks.filter(data.matched$sample2, combined, p.value.impute=sig.value.impute, overlap.ratio)$merge2
2238
2239 # choose the ones with significant value > sig.value.impute, which are the overlapped ones
2240
2241 combined.in1 <- overlap.1[overlap.1$sig.value > sig.value.impute, ]
2242 combined.in2 <- overlap.2[overlap.2$sig.value > sig.value.impute, ]
2243
2244 combined.region <- rbind(combined.in1, combined.in2)
2245
2246 is.repeated <- duplicated(combined.region$start)
2247 combined.region <- combined.region[!is.repeated,]
2248
2249 # now rank the peaks in matched region
2250 rank.combined <- rank(-combined.region$sig.value)
2251
2252 # now transform the parameters estimated into the new scale
2253 npeaks.overlap <- nrow(combined.region)
2254 npeaks.consistent <- nrow(cbind(em.output[[j]]$data.pruned$sample1))
2255
2256
2257 # the breaks are the same for x and y
2258 f1 <- list(breaks=em.output[[j]]$em.fit$x.mar$f1$breaks*npeaks.overlap/npeaks.consistent, density=(em.output[[j]]$em.fit$x.mar$f1$density+em.output[[j]]$em.fit$y.mar$f1$density)/2)
2259 # the first break boundary goes up when changing scale, need set it back to be a bit smaller than 1
2260 f1$breaks[1] <- min(f1$breaks[1], 0.95)
2261
2262 f2 <- list(breaks=em.output[[j]]$em.fit$x.mar$f2$breaks*npeaks.overlap/npeaks.consistent, density=(em.output[[j]]$em.fit$x.mar$f2$density+em.output[[j]]$em.fit$y.mar$f2$density)/2)
2263 # the first break boundary goes up when changing scale, need set it back to be a bit smaller than 1
2264 f2$breaks[1] <- min(f2$breaks[1], 0.95)
2265
2266 p <- em.output[[j]]$em.fit$para$p
2267
2268 # find the posterior probability
2269 errorprob.combined <- get.comp2.prob(rank.combined, p, f1, f2)
2270
2271 # compute the FDR and find cutoff of posterior prob and the sig value
2272 o <- order(errorprob.combined)
2273 idr <- cumsum(errorprob.combined[o])/c(1:length(o))
2274 idr.index <- which(idr > idr.level)[1]
2275 errorprob.cutoff <- errorprob.combined[o][idr.index]
2276
2277 # find the minimum significant measure among selected peaks
2278 sig.value <- min(combined.region$sig.value[o][1:idr.index])
2279 # sig.value <- quantile(combined.region$sig.value[o][1:idr.index], prob=0.05)
2280 #sig.value <- quantile(combined.region$sig.value[errorprob.combined<=errorprob.cutoff], prob=0.05)
2281
2282 # apply the significant value on the whole pooled list
2283 combined.selected <- combined[combined$sig.value >= sig.value,]
2284
2285 combined.selected.all <- rbind(combined.selected.all, combined.selected)
2286 }
2287
2288 is.repeated <- duplicated(combined.selected.all$start)
2289 combined.selected.all <- combined.selected.all[!is.repeated,]
2290
2291 npeak <- nrow(combined.selected.all)
2292
2293 npeak.stat <- list(idr.level=idr.level, npeak=npeak)
2294
2295 sig.combined <- c(min(combined.selected.all[,"sig.value"], na.rm=T), max(combined.selected.all[,"sig.value"], na.rm=T))
2296
2297 # idr.combined <- c(min(combined.selected.all[,"q.value"], na.rm=T), max(combined.selected.all[,"q.value"], na.rm=T))
2298 # combined.selected.all <- deconcatenate.chr(combined.selected.all, chr.size)[,c("chr", "start", "stop", "signal.value", "p.value", "q.value")]
2299
2300 combined.selected.all <- combined.selected.all[, c("chr", "start.ori", "stop.ori", "signal.value", "p.value", "q.value")]
2301 colnames(combined.selected.all) <- c("chr", "start", "stop", "signal.value", "p.value", "q.value")
2302
2303 invisible(list(npeak.stat=npeak.stat, combined.selected=combined.selected.all, sig.combined=sig.combined))
2304 }
2305
2306
2307
2308 # get the posterior probability of the 2nd component
2309 get.comp2.prob <- function(x, p, f1, f2){
2310
2311 # get pdf and cdf of each component from functions in the corresponding component
2312 px.1 <- sapply(x, get.pdf, df=f1)
2313 px.2 <- sapply(x, get.pdf, df=f2)
2314
2315 comp2prob <- 1 - p*px.1/(p*px.1+(1-p)*px.2)
2316
2317 return(comp2prob)
2318 }
2319
2320 keep.match <- function(sample1, sample2, sig.value.impute=0){
2321
2322 sample1.prune <- sample1[sample1$sig.value > sig.value.impute & sample2$sig.value > sig.value.impute,]
2323 sample2.prune <- sample2[sample1$sig.value > sig.value.impute & sample2$sig.value > sig.value.impute,]
2324
2325 invisible(list(sample1=sample1.prune, sample2=sample2.prune))
2326 }
2327
2328
2329 ##############################################
2330 #
2331 # The following is for simulation
2332 #
2333 ##############################################
2334
2335
2336 # simulate gaussian copula
2337 # u is the uniform random variable and rho is correlation coefficient
2338 simu.gaussian.copula <- function(u, rho){
2339
2340 n <- length(u)
2341
2342 # simulate y given x=qnorm(u)
2343 y <- qnorm(u)*rho + rnorm(n)*sqrt(1-rho^2)
2344
2345 v <- pnorm(y)
2346
2347 invisible(v)
2348 }
2349
2350 ## simulate Clayton copula from its generating function
2351 ## Genest and MacKay (1986)
2352
2353 phi.ori <- function(t, s){
2354
2355 (t^(-s) -1)/s
2356 }
2357
2358
2359 phi.inv <- function(y, s){
2360
2361 exp(-log(s*y+1)/s)
2362 }
2363
2364 phi.der <- function(t, s){
2365
2366 -t^(-s-1)
2367 }
2368
2369 phi.der.inv <- function(y, s){
2370
2371 exp(log(-y)/(-s-1))
2372 }
2373
2374 get.w <- function(u, t, s){
2375
2376 phi.der.inv(phi.der(u, s)/t, s)
2377 }
2378
2379 get.v <- function(w, u, s){
2380
2381 phi.inv(phi.ori(w, s) - phi.ori(u, s), s)
2382 }
2383
2384 # u is a uniform random variable, s is the association parameter
2385 simu.clayton.copula <- function(u, s){
2386
2387 t <- runif(length(u))
2388
2389 if(s>0){
2390 w <- get.w(u, t, s)
2391 v <- get.v(w, u, s)
2392 return(v)
2393 }
2394
2395 if(s==0){
2396 return(t)
2397 }
2398
2399 if(s <0){
2400 print("Invalid association parameters for clayton copula")
2401 }
2402
2403 }
2404
2405
2406
2407 ###### 09-09-09
2408
2409 # simulate a two-component copula mixture:
2410 # - marginal distributions for the two variables in each component are both
2411 # normal and with the same parameters
2412 # p is the mixing proportion of component 1
2413 # n is the total sample size
2414 simu.copula.2mix <- function(s1, s2, p, n, mu1, mu2, sd1, sd2, copula.txt){
2415
2416 n1 <- round(n*p)
2417 n2 <- n-n1
2418
2419 u1 <- runif(n1)
2420
2421 if(copula.txt =="clayton")
2422 v1 <- simu.clayton.copula(u1, s1)
2423 else{
2424 if(copula.txt =="gaussian")
2425 v1 <- simu.gaussian.copula(u1, s1)
2426 }
2427
2428 u2 <- runif(n2)
2429
2430 if(copula.txt =="clayton")
2431 v2 <- simu.clayton.copula(u2, s2)
2432 else{
2433 if(copula.txt =="gaussian")
2434 v2 <- simu.gaussian.copula(u2, s2)
2435 }
2436
2437 # generate test statistics
2438 sample1.1 <- qnorm(u1, mu1, sd1)
2439 sample1.2 <- qnorm(v1, mu1, sd1)
2440
2441 sample2.1 <- qnorm(u2, mu2, sd2)
2442 sample2.2 <- qnorm(v2, mu2, sd2)
2443
2444 return(list(u=c(u1, u2), v=c(v1, v2),
2445 u.inv=c(sample1.1, sample2.1), v.inv=c(sample1.2, sample2.2),
2446 label=c(rep(1, n1), rep(2, n2))))
2447 }
2448
2449 # using inverse of the cdf to generate original observations
2450
2451 simu.copula.2mix.inv <- function(s1, s2, p, n, cdf1.x, cdf1.y, cdf2.x, cdf2.y, copula.txt){
2452
2453 n1 <- round(n*p)
2454 n2 <- n-n1
2455
2456 u1 <- runif(n1)
2457
2458 if(copula.txt =="clayton")
2459 v1 <- simu.clayton.copula(u1, s1)
2460 else{
2461 if(copula.txt =="gaussian")
2462 v1 <- simu.gaussian.copula(u1, s1)
2463 }
2464
2465 u2 <- runif(n2)
2466
2467 if(copula.txt =="clayton")
2468 v2 <- simu.clayton.copula(u2, s2)
2469 else{
2470 if(copula.txt =="gaussian")
2471 v2 <- simu.gaussian.copula(u2, s2)
2472 }
2473
2474 # generate test statistics
2475 # sample1.1 <- qnorm(u1, mu1, sd1)
2476 # sample1.2 <- qnorm(v1, mu1, sd1)
2477
2478 # sample2.1 <- qnorm(u2, mu2, sd2)
2479 # sample2.2 <- qnorm(v2, mu2, sd2)
2480
2481 sample1.x <- inv.cdf.vec(u1, cdf1.x)
2482 sample1.y <- inv.cdf.vec(v1, cdf1.y)
2483
2484 sample2.x <- inv.cdf.vec(u2, cdf2.x)
2485 sample2.y <- inv.cdf.vec(v2, cdf2.y)
2486
2487
2488 return(list(u=c(u1, u2), v=c(v1, v2),
2489 u.inv=c(sample1.x, sample2.x), v.inv=c(sample1.y, sample2.y),
2490 label=c(rep(1, n1), rep(2, n2))))
2491 }
2492
2493 # obtain original observation by converting cdf into quantiles
2494 # u is one cdf
2495 # u.cdf is a cdf (assuming it is a histogram) and has the break points (cdf$cdf and cdf$breaks)
2496 # the smallest value of cdf=0 and the largest =1
2497 inv.cdf <- function(u, u.cdf){
2498
2499 # which bin it falls into
2500 i <- which(u.cdf$cdf> u)[1]
2501 q.u <- (u - u.cdf$cdf[i-1])/(u.cdf$cdf[i] - u.cdf$cdf[i-1])* (u.cdf$breaks[i]-u.cdf$breaks[i-1]) + u.cdf$breaks[i-1]
2502
2503 return(q.u)
2504 }
2505
2506 inv.cdf.vec <- function(u, u.cdf){
2507
2508 # check if cdf has the right range (0, 1)
2509 ncdf <- length(u.cdf$cdf)
2510 nbreaks <- length(u.cdf$breaks)
2511
2512 if(ncdf == nbreaks-1 & u.cdf$cdf[ncdf]< 1)
2513 u.cdf[ncdf] <- 1
2514
2515 q.u <- sapply(u, inv.cdf, u.cdf)
2516
2517 return(q.u)
2518 }
2519
2520 # here we simulate a likely real situation
2521 # the test statistics from two normal distributions
2522 # according to their labels, then convert them into p-values w.r.t H0 using
2523 # one-sided test.
2524 # The test statistics are correlated for the signal component and independent
2525 # for the noise component
2526 # For the signal component, Y = X + eps, where eps ~ N(0, sigma^2)
2527 simu.test.stat <- function(p, n, mu1, sd1, mu0, sd0, sd.e){
2528
2529 # first component - signal
2530 n.signal <- round(n*p)
2531 n.noise <- n - n.signal
2532
2533 # labels
2534 labels <- c(rep(1, n.signal), rep(0, n.noise))
2535
2536 # test statistics for signal and noise
2537 mu.signal <- rnorm(n.signal, mu1, sd1)
2538 x.signal <- mu.signal + rnorm(n.signal, 0, sd.e)
2539 x.noise <- rnorm(n.noise, mu0, sd0) + rnorm(n.noise, 0, sd.e)
2540
2541 y.signal <- mu.signal + rnorm(n.signal, 0, sd.e)
2542 # sd.e can be dependent on signal
2543 y.noise <- rnorm(n.noise, mu0, sd0) + rnorm(n.noise, 0, sd.e)
2544
2545 # concatenate
2546 x <- c(x.signal, x.noise)
2547 y <- c(y.signal, y.noise)
2548
2549 # convert to p-values based on H0
2550 p.x <- 1-pnorm(x, mu0, sqrt(sd0^2+sd.e^2))
2551 p.y <- 1-pnorm(y, mu0, sqrt(sd0^2+sd.e^2))
2552
2553 return(list(p.x=p.x, p.y=p.y, x=x, y=y, labels=labels))
2554
2555 }
2556
2557 # compute the tradeoff and calibration
2558 forward.decoy.tradeoff.ndecoy <- function(xx, labels, ndecoy){
2559
2560 xx <- round(xx, 5)
2561 o <- order(xx, decreasing=T)
2562
2563 rand <- 1-labels # if rand==0, consistent
2564 # order the random indicator in the same order
2565 rand.o <- rand[o]
2566
2567 if(sum(rand.o) > ndecoy){
2568 index.decoy <- which(cumsum(rand.o)==ndecoy)
2569 } else {
2570 index.decoy <- which(cumsum(rand.o)==sum(rand.o))
2571 }
2572
2573 cutoff.decoy <- xx[o][index.decoy]
2574
2575 # only consider the unique ones
2576 cutoff.unique <- unique(xx[o])
2577
2578 cutoff <- cutoff.unique[cutoff.unique >= cutoff.decoy[length(cutoff.decoy)]]
2579
2580 get.decoy.count <- function(cut.off){
2581 above <- rep(0, length(xx))
2582 above[xx >= cut.off] <- 1
2583 decoy.count <- sum(above==1 & rand==1)
2584 return(decoy.count)
2585 }
2586
2587 get.forward.count <- function(cut.off){
2588 above <- rep(0, length(xx))
2589 above[xx >= cut.off] <- 1
2590 forward.count <- sum(above==1 & rand==0)
2591 return(forward.count)
2592 }
2593
2594 get.est.fdr <- function(cut.off){
2595 above <- rep(0, length(xx))
2596 above[xx >= cut.off] <- 1
2597 est.fdr <- 1-mean(xx[above==1])
2598 return(est.fdr)
2599 }
2600
2601 # assuming rand=0 is right
2602 get.false.neg.count <- function(cut.off){
2603 below <- rep(0, length(xx))
2604 below[xx < cut.off] <- 1
2605 false.neg.count <- sum(below==1 & rand==0)
2606 return(false.neg.count)
2607 }
2608
2609 get.false.pos.count <- function(cut.off){
2610 above <- rep(0, length(xx))
2611 above[xx >= cut.off] <- 1
2612 false.pos.count <- sum(above==1 & rand==1)
2613 return(false.pos.count)
2614 }
2615
2616 decoy <- sapply(cutoff, get.decoy.count)
2617 forward <- sapply(cutoff, get.forward.count)
2618
2619 est.fdr <- sapply(cutoff, get.est.fdr)
2620 emp.fdr <- decoy/(decoy+forward)
2621
2622 # compute specificity and sensitivity
2623 # assuming rand=1 is wrong and rand=0 is right
2624 false.neg <- sapply(cutoff, get.false.neg.count)
2625 false.pos <- sapply(cutoff, get.false.pos.count)
2626
2627 true.pos <- sum(rand==0)-false.neg
2628 true.neg <- sum(rand==1)-false.pos
2629
2630 sensitivity <- true.pos/(true.pos+false.neg)
2631 specificity <- true.neg/(true.neg+false.pos)
2632
2633 return(list(decoy=decoy, forward=forward, cutoff=cutoff, est.fdr=est.fdr, emp.fdr=emp.fdr, sensitivity=sensitivity, specificity=specificity))
2634 }
2635
2636
2637 # compute the em for jackknife and all data, and find FDR
2638 get.emp.jack <- function(a, p0){
2639
2640 nobs <- length(a$labels)
2641 est <- list()
2642 est.all <- list()
2643
2644 temp.all <- em.transform(-a$p.x, -a$p.y, mu=1.5, sigma=1.4, rho=0.4, p=0.7, eps=0.01)
2645 # temp.all <- em.2copula.quick(a$p.x, a$p.y, p0=p0, rho1.0=0.7,
2646 # rho2.0=0, eps=0.01, fix.p=T, stoc=F, fix.rho2=T, "gaussian")
2647
2648 est.all$p <- temp.all$para$p
2649 est.all$rho1 <- temp.all$para$rho1
2650 est.all$FDR <- get.FDR(temp.all$e.z)
2651
2652 FDR <- list()
2653 p <- c()
2654 rho1 <- c()
2655
2656
2657 for(i in 1:nobs){
2658
2659 temp <- em.transform(-a$p.x[-i], -a$p.y[-i], mu=1.5, sigma=1.4, rho=0.4, p=0.7, eps=0.01)
2660 # temp <- em.2copula.quick(a$p.x[-i], a$p.y[-i], p0=p0, rho1.0=0.7,
2661 # rho2.0=0, eps=0.01, fix.p=T, stoc=F, fix.rho2=T, "gaussian")
2662
2663 est[[i]] <- list(p=temp$para$p, rho1=temp$para$rho1, FDR=get.FDR(temp$e.z))
2664
2665 FDR[[i]] <- est[[i]]$FDR # this is the FDR for top n peaks
2666 p[i] <- est[[i]]$p
2667 rho1[i] <- est[[i]]$rho1
2668 }
2669
2670 est.jack <- list(FDR=FDR, p=p, rho1=rho1)
2671 return(list(est.jack=est.jack, est.all=est.all))
2672 }
2673
2674
2675 # get the npeaks corresponding to the nominal FDR estimated from the sample
2676 # and find the corresponding FDR from the entire data
2677 get.FDR.jack <- function(est, FDR.nominal){
2678
2679 nobs <- length(est$est.jack$FDR)
2680 FDR.all <- c()
2681 top.n <- c()
2682
2683 for(i in 1:nobs){
2684 top.n[i] <- max(which(est$est.jack$FDR[[i]] <= FDR.nominal))
2685 FDR.all[i] <- est$est.all$FDR[top.n[i]]
2686 }
2687
2688 invisible(list(FDR.all=FDR.all, top.n=top.n))
2689 }
2690
2691 # compute Jackknife peudonumber
2692 # a is the dataset
2693 get.emp.IF <- function(a, p0){
2694
2695 nobs <- length(a$labels)
2696 est <- list()
2697 est.all <- list()
2698
2699 temp.all <- em.2copula.quick(a$p.x, a$p.y, p0=p0, rho1.0=0.7,
2700 rho2.0=0, eps=0.01, fix.p=T, stoc=F, fix.rho2=T, "gaussian")
2701
2702 est.all$p <- temp.all$para$p
2703 est.all$rho1 <- temp.all$para$rho1
2704 est.all$FDR <- get.FDR(temp.all$e.z)
2705
2706 IF.FDR <- list()
2707 IF.p <- c()
2708 IF.rho1 <- c()
2709
2710 for(i in 1:nobs){
2711
2712 temp <- em.2copula.quick(a$p.x[-i], a$p.y[-i], p0=p0, rho1.0=0.7,
2713 rho2.0=0, eps=0.01, fix.p=T, stoc=F, fix.rho2=T, "gaussian")
2714
2715 est[[i]] <- list(p=temp$para$p, rho1=temp$para$rho1, FDR=get.FDR(temp$e.z))
2716
2717 IF.FDR[[i]] <- (nobs-1)*(est.all$FDR[-nobs] - est[[i]]$FDR) # this is the FDR for top n peaks
2718 IF.p[i] <- (nobs-1)*(est.all$p - est[[i]]$p)
2719 IF.rho1[i] <- (nobs-1)*(est.all$rho1 - est[[i]]$rho1)
2720 }
2721
2722 emp.IF <- list(FDR=IF.FDR, p=IF.p, rho1=IF.rho1)
2723
2724 invisible(list(emp.IF=emp.IF, est.all=est.all, est=est))
2725 }
2726
2727 # e.z is the posterior probability of being in signal component
2728 get.FDR <- function(e.z){
2729
2730 e.z.o <- order(1-e.z)
2731 FDR <- cumsum(1-e.z[e.z.o])/c(1:length(e.z.o))
2732
2733 invisible(FDR)
2734 }
2735
2736 # get the FDR of selecting the top n peaks
2737 # IF.est is the sample influence function
2738 # top.n
2739 get.IF.FDR <- function(IF.est, top.n){
2740
2741 nobs <- length(IF.est$emp.IF$FDR)
2742 FDR <- c()
2743
2744 # influence function of p
2745 for(i in 1:nobs)
2746 FDR[i] <- IF.est$emp.IF$FDR[[i]][top.n]
2747
2748 invisible(FDR)
2749 }
2750
2751 # get the sample influence function for FDR at a given FDR size
2752 # 1. find the number of peaks selected at a given FDR computed from all obs
2753 # 2. use the number to find the sample influence function for FDR
2754 # IF.est$est.all is the FDR with all peaks
2755 get.IF.FDR.all <- function(IF.est, FDR.size){
2756
2757 top.n <- which.min(abs(IF.est$est.all$FDR -FDR.size))
2758 nobs <- length(IF.est$est.all$FDR)
2759 FDR <- c()
2760
2761 # influence function of p
2762 for(i in 1:nobs)
2763 FDR[i] <- IF.est$emp.IF$FDR[[i]][top.n]
2764
2765 invisible(list(FDR=FDR, top.n=top.n))
2766 }
2767
2768 plot.simu.uri <- function(x, y){
2769
2770 tt <- seq(0.01, 0.99, by=0.01)
2771 uri <- sapply(tt, comp.uri.prob, u=x, v=y)
2772 uri.thin <- uri[seq(1, length(tt), by=3)]
2773 tt.thin <- tt[seq(1, length(tt), by=3)]
2774 duri <- (uri.thin[-1]-uri.thin[-length(uri.thin)])/(tt.thin[-1]-tt.thin[-length(tt.thin)])
2775 uri.spl <- smooth.spline(tt, uri, df=6.4)
2776 uri.der <- predict(uri.spl, tt, deriv=1)
2777
2778 par(mfrow=c(2,2))
2779 plot(x[1:n0], y[1:n0])
2780 points(x[(n0+1):n], y[(n0+1):n], col=2)
2781 plot(rank(-x)[1:n0], rank(-y)[1:n0])
2782 points(rank(-x)[(1+n0):n], rank(-y)[(1+n0):n])
2783 plot(tt, uri)
2784 lines(c(0,1), c(0,1), lty=2)
2785 title(paste("rho1=", rho1, " rho2=", rho2, "p=", p, sep=""))
2786 plot(tt.thin[-1], duri)
2787 lines(uri.der)
2788 abline(h=1)
2789 invisible(list(x=x, y=y, uri=uri, tt=tt, duri=duri, tt.thin=tt.thin, uri.der=uri.der))
2790
2791 }
2792
2793
2794 ###### new fitting procedure
2795
2796
2797
2798
2799 # 1. rank pairs
2800
2801 # 2. initialization
2802 # 3. convert to pseudo-number
2803
2804 # 4. EM
2805
2806 # need plugin and test
2807 # find the middle point between the bins
2808 get.pseudo.mix <- function(x, mu, sigma, rho, p){
2809
2810
2811 # first compute cdf for points on the grid
2812 # generate 200 points between [-3, mu+3*sigma]
2813 nw <- 1000
2814 w <- seq(min(-3, mu-3*sigma), max(mu+3*sigma, 3), length=nw)
2815 w.cdf <- p*pnorm(w, mean=mu, sd=sigma) + (1-p)*pnorm(w, mean=0, sd=1)
2816
2817 i <- 1
2818
2819 quan.x <- rep(NA, length(x))
2820
2821 for(i in c(1:nw)){
2822 index <- which(x >= w.cdf[i] & x < w.cdf[i+1])
2823 quan.x[index] <- (x[index]-w.cdf[i])*(w[i+1]-w[i])/(w.cdf[i+1]-w.cdf[i]) +w[i]
2824 }
2825
2826 index <- which(x < w.cdf[1])
2827 if(length(index)>0)
2828 quan.x[index] <- w[1]
2829
2830 index <- which(x > w.cdf[nw])
2831 if(length(index)>0)
2832 quan.x[index] <- w[nw]
2833
2834 # linear.ext <- function(x, w, w.cdf){
2835 # linear interpolation
2836 # index.up <- which(w.cdf>= x)[1]
2837 # left.index <- which(w.cdf <=x)
2838 # index.down <- left.index[length(left.index)]
2839 # quan.x <- (w[index.up] + w[index.down])/2
2840 # }
2841
2842 # x.pseudo <- sapply(x, linear.ext, w=w, w.cdf=w.cdf)
2843
2844 # invisible(x.pseudo)
2845 invisible(quan.x)
2846 }
2847
2848
2849 # EM to compute the latent structure
2850 # steps:
2851 # 1. raw values are first transformed into pseudovalues
2852 # 2. EM is used to compute the underlining structure, which is a mixture
2853 # of two normals
2854 em.transform <- function(x, y, mu, sigma, rho, p, eps){
2855
2856 x.cdf.func <- ecdf(x)
2857 y.cdf.func <- ecdf(y)
2858 afactor <- length(x)/(length(x)+1)
2859 x.cdf <- x.cdf.func(x)*afactor
2860 y.cdf <- y.cdf.func(y)*afactor
2861
2862 # initialization
2863 para <- list()
2864 para$mu <- mu
2865 para$sigma <- sigma
2866 para$rho <- rho
2867 para$p <- p
2868
2869 j <- 1
2870 to.run <- T
2871 loglik.trace <- c()
2872 loglik.inner.trace <- c()
2873
2874 #to.run.inner <- T
2875 z.1 <- get.pseudo.mix(x.cdf, para$mu, para$sigma, para$rho, para$p)
2876 z.2 <- get.pseudo.mix(y.cdf, para$mu, para$sigma, para$rho, para$p)
2877
2878 # cat("length(z1)", length(z.1), "\n")
2879 while(to.run){
2880
2881 # get pseudo value in each cycle
2882 # z.1 <- get.pseudo.mix(x.cdf, para$mu, para$sigma, para$rho, para$p)
2883 # z.2 <- get.pseudo.mix(y.cdf, para$mu, para$sigma, para$rho, para$p)
2884
2885 i <- 1
2886 while(to.run){
2887
2888 # EM for latent structure
2889 e.z <- e.step.2normal(z.1, z.2, para$mu, para$sigma, para$rho, para$p)
2890 para <- m.step.2normal(z.1, z.2, e.z)
2891 #para$rho <- rho
2892 #para$p <- p
2893 #para$mu <- mu
2894 #para$sigma <- sigma
2895 if(i > 1)
2896 l.old <- l.new
2897
2898 # this is just the mixture likelihood of two-component Gaussian
2899 l.new <- loglik.2binormal(z.1, z.2, para$mu, para$sigma, para$rho, para$p)
2900
2901 loglik.inner.trace[i] <- l.new
2902
2903 if(i > 1){
2904 to.run <- loglik.inner.trace[i]-loglik.inner.trace[i-1]>eps
2905 }
2906
2907
2908 # if(i > 2){
2909 # l.inf <- loglik.inner.trace[i-2] + (loglik.inner.trace[i-1] - loglik.inner.trace[i-2])/(1-(loglik.inner.trace[i]-loglik.inner.trace[i-1])/(loglik.inner.trace[i-1]-loglik.inner.trace[i-2]))
2910
2911 # if(loglik.inner.trace[i-1]!=loglik.inner.trace[i-2])
2912 # to.run <- abs(l.inf - loglik.inner.trace[i]) > eps
2913 # else
2914 # to.run <- F
2915
2916 # }
2917
2918 cat("loglik.inner.trace[", i, "]=", loglik.inner.trace[i], "\n")
2919 cat("mu=", para$mu, "sigma=", para$sigma, "p=", para$p, "rho=", para$rho, "\n\n")
2920
2921 i <- i+1
2922 }
2923
2924
2925 # get pseudo value in each cycle
2926 z.1 <- get.pseudo.mix(x.cdf, para$mu, para$sigma, para$rho, para$p)
2927 z.2 <- get.pseudo.mix(y.cdf, para$mu, para$sigma, para$rho, para$p)
2928
2929 if(j > 1)
2930 l.old.outer <- l.new.outer
2931
2932 l.new.outer <- loglik.2binormal(z.1, z.2, para$mu, para$sigma, para$rho, para$p)
2933
2934 loglik.trace[j] <- l.new.outer
2935
2936 if(j == 1)
2937 to.run <- T
2938 else{ # stop when iteration>100
2939 if(j > 100)
2940 to.run <- F
2941 else
2942 to.run <- l.new.outer - l.old.outer > eps
2943 }
2944
2945 # if(j %% 10==0)
2946 cat("loglik.trace[", j, "]=", loglik.trace[j], "\n")
2947 cat("mu=", para$mu, "sigma=", para$sigma, "p=", para$p, "rho=", para$rho, "\n")
2948
2949 j <- j+1
2950 }
2951
2952 bic <- -2*l.new + 4*log(length(z.1))
2953
2954 return(list(para=list(p=para$p, rho=para$rho, mu=para$mu, sigma=para$sigma),
2955 loglik=l.new, bic=bic, e.z=e.z, loglik.trace=loglik.trace))
2956 }
2957
2958
2959
2960
2961 # compute log-likelihood for mixture of two bivariate normals
2962 loglik.2binormal <- function(z.1, z.2, mu, sigma, rho, p){
2963
2964 l.m <- sum(d.binormal(z.1, z.2, 0, 1, 0)+log(p*exp(d.binormal(z.1, z.2, mu, sigma, rho)-d.binormal(z.1, z.2, 0, 1, 0))+(1-p)))
2965
2966 # l.m <- sum((p*d.binormal(z.1, z.2, mu, sigma, rho) + (1-p)*d.binormal(z.1, z.2, 0, 1, 0)))
2967 return(l.m)
2968 }
2969
2970 # check this when rho=1
2971
2972 # density of binomial distribution with equal mean and sigma on both dimensions
2973 d.binormal <- function(z.1, z.2, mu, sigma, rho){
2974
2975 loglik <- (-log(2)-log(pi)-2*log(sigma) - log(1-rho^2)/2 - (0.5/(1-rho^2)/sigma^2)*((z.1-mu)^2 -2*rho*(z.1-mu)*(z.2-mu) + (z.2-mu)^2))
2976
2977 return(loglik)
2978 }
2979
2980 # E-step for computing the latent strucutre
2981 # e.z is the prob to be in the consistent group
2982 # e.step for estimating posterior prob
2983 # z.1 and z.2 can be vectors or scalars
2984 e.step.2normal <- function(z.1, z.2, mu, sigma, rho, p){
2985
2986 e.z <- p/((1-p)*exp(d.binormal(z.1, z.2, 0, 1, 0)-d.binormal(z.1, z.2, mu, sigma, rho))+ p)
2987
2988 invisible(e.z)
2989 }
2990
2991 # M-step for computing the latent structure
2992 # m.step for estimating proportion, mean, sd and correlation coefficient
2993 m.step.2normal <- function(z.1, z.2, e.z){
2994
2995 p <- mean(e.z)
2996 mu <- sum((z.1+z.2)*e.z)/2/sum(e.z)
2997 sigma <- sqrt(sum(e.z*((z.1-mu)^2+(z.2-mu)^2))/2/sum(e.z))
2998 rho <- 2*sum(e.z*(z.1-mu)*(z.2-mu))/(sum(e.z*((z.1-mu)^2+(z.2-mu)^2)))
2999
3000 return(list(p=p, mu=mu, sigma=sigma, rho=rho))
3001 }
3002
3003
3004 # assume top p percent of observations are true
3005 # x and y are ranks, estimate
3006 init <- function(x, y, x.label){
3007
3008 x.o <- order(x)
3009
3010 x.ordered <- x[x.o]
3011 y.ordered <- y[x.o]
3012 x.label.ordered <- x.label[x.o]
3013
3014 n <- length(x)
3015 p <- sum(x.label)/n
3016
3017 rho <- cor(x.ordered[1:ceiling(p*n)], y.ordered[1:ceiling(p*n)])
3018
3019 temp <- find.mu.sigma(x.ordered, x.label.ordered)
3020 mu <- temp$mu
3021 sigma <- temp$sigma
3022
3023 invisible(list(mu=mu, sigma=sigma, rho=rho, p=p))
3024
3025 }
3026
3027 # find mu and sigma if the distributions of marginal ranks are known
3028 # take the medians of the two dist and map back to the original
3029 init.dist <- function(f0, f1){
3030
3031 # take the median in f0
3032 index.median.0 <- which(f0$cdf>0.5)[1]
3033 q.0.small <- f0$cdf[index.median.0] # because f0 and f1 have the same bins
3034 q.1.small <- f1$cdf[index.median.0]
3035
3036 # take the median in f1
3037 index.median.1 <- which(f1$cdf>0.5)[1]
3038 q.0.big <- f0$cdf[index.median.1] # because f0 and f1 have the same bins
3039 q.1.big <- f1$cdf[index.median.1]
3040
3041 # find pseudo value for x.middle[1] on normal(0,1)
3042 pseudo.small.0 <- qnorm(q.0.small, mean=0, sd=1)
3043 pseudo.small.1 <- qnorm(q.1.small, mean=0, sd=1)
3044
3045 # find pseudo value for x.middle[2] on normal(0,1)
3046 pseudo.big.0 <- qnorm(q.0.big, mean=0, sd=1)
3047 pseudo.big.1 <- qnorm(q.1.big, mean=0, sd=1)
3048
3049 mu <- (pseudo.small.0*pseudo.big.1 - pseudo.small.1*pseudo.big.0)/(pseudo.big.1-pseudo.small.1)
3050
3051 sigma <- (pseudo.small.0-mu)/pseudo.small.1
3052
3053 return(list(mu=mu, sigma=sigma))
3054 }
3055
3056 # generate labels
3057
3058 # find the part of data with overlap
3059
3060 # find the percentile on noise and signal
3061
3062 # Suppose there are signal and noise components, with mean=0 and sd=1 for noise
3063 # x and x.label are the rank of the observations and their labels,
3064 # find the mean and sd of the other component
3065 # x.label takes values of 0 and 1
3066 find.mu.sigma <- function(x, x.label){
3067
3068 x.0 <- x[x.label==0]
3069 x.1 <- x[x.label==1]
3070
3071 n.x0 <- length(x.0)
3072 n.x1 <- length(x.1)
3073
3074 x.end <- c(min(x.0), min(x.1), max(x.0), max(x.1))
3075 o <- order(x.end)
3076 x.middle <- x.end[o][c(2,3)]
3077
3078 # the smaller end of the overlap
3079 q.1.small <- mean(x.1 <= x.middle[1])*n.x1/(n.x1+1)
3080 q.0.small <- mean(x.0 <= x.middle[1])*n.x0/(n.x0+1)
3081
3082 # the bigger end of the overlap
3083 q.1.big <- mean(x.1 <= x.middle[2])*n.x1/(n.x1+1)
3084 q.0.big <- mean(x.0 <= x.middle[2])*n.x0/(n.x0+1)
3085
3086 # find pseudo value for x.middle[1] on normal(0,1)
3087 pseudo.small.0 <- qnorm(q.0.small, mean=0, sd=1)
3088 pseudo.small.1 <- qnorm(q.1.small, mean=0, sd=1)
3089
3090 # find pseudo value for x.middle[2] on normal(0,1)
3091 pseudo.big.0 <- qnorm(q.0.big, mean=0, sd=1)
3092 pseudo.big.1 <- qnorm(q.1.big, mean=0, sd=1)
3093
3094 mu <- (pseudo.small.0*pseudo.big.1 - pseudo.small.1*pseudo.big.0)/(pseudo.big.1-pseudo.small.1)
3095
3096 sigma <- (pseudo.small.0-mu)/pseudo.small.1
3097
3098 return(list(mu=mu, sigma=sigma))
3099 }