3
+ − 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 }
