--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/functions-all-clayton-12-13.r	Thu Jan 17 15:59:06 2013 -0500
@@ -0,0 +1,3099 @@
+# revised on 2-20-10
+# - fix error in pass.structure: reverse rank.combined, so that big sig.value
+#  are ranked with small numbers (1, 2, ...)
+# - fix error on get.ez.tt.all: get ez.cutoff from sorted e.z
+
+#
+# modified EM procedure to compute empirical CDF more precisely - 09/2009
+
+
+
+# this file contains the functions for  
+# 1. computing the correspondence profile (upper rank intersection and derivatives)
+# 2. inference of copula mixture model
+#
+# It also has functions for
+# 1. reading peak caller results
+# 2. processing and matching called peaks
+# 3. plotting results
+
+
+################ read peak caller results
+
+# process narrow peak format
+# some peak callers may not report q-values, p-values or fold of enrichment
+# need further process before comparison
+#
+# stop.exclusive: Is the basepair of peak.list$stop exclusive? In narrowpeak and broadpeak format they are exclusive.
+# 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 
+# adjacent peaks, which creates trouble for finding correct intersect  
+process.narrowpeak <- function(narrow.file, chr.size, half.width=NULL, summit="offset", stop.exclusive=T, broadpeak=F){
+
+
+  aa <- read.table(narrow.file)
+
+  if(broadpeak){
+    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)
+  }else{
+    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)
+  }
+    
+  if(summit=="summit"){
+    bb.ori$summit <- bb.ori$summit-bb.ori$start # change summit to offset to avoid error when concatenating chromosomes
+  }
+ 
+  bb <- concatenate.chr(bb.ori, chr.size)
+
+  #bb <- bb.ori
+
+  # remove the peaks that has the same start and stop value
+  bb <- bb[bb$start != bb$stop,]
+
+  if(stop.exclusive==T){
+    bb$stop <- bb$stop-1
+  }
+
+  if(!is.null(half.width)){
+    bb$start.ori <- bb$start    
+    bb$stop.ori <- bb$stop 
+
+    # if peak is narrower than the specified window, stay with its width
+    # otherwise chop wider peaks to specified width
+    width <- bb$stop-bb$start +1
+    is.wider <- width > 2*half.width
+
+    if(summit=="offset" | summit=="summit"){ # if summit is offset from start
+      bb$start[is.wider] <- bb$start.ori[is.wider] + bb$summit[is.wider]-half.width
+      bb$stop[is.wider] <- bb$start.ori[is.wider] + bb$summit[is.wider]+half.width
+    } else { 
+      if(summit=="unknown"){
+        bb$start[is.wider] <- bb$start.ori[is.wider]+round(width[is.wider]/2) - half.width
+        bb$stop[is.wider] <- bb$start.ori[is.wider]+round(width[is.wider]/2) + half.width
+      }
+    }
+  }
+
+  bb <- clean.data(bb)
+  invisible(list(data.ori=bb.ori, data.cleaned=bb))
+}
+
+# clean data 
+# and concatenate chromosomes if needed
+clean.data <- function(adata){
+
+  # remove the peaks that has the same start and stop value
+  adata <- adata[adata$start != adata$stop,]
+
+  # if some stops and starts are the same, need fix them
+  stop.in.start <- is.element(adata$stop, adata$start)
+  n.fix <- sum(stop.in.start)
+  if(n.fix >0){
+    print(paste("Fix", n.fix, "stops\n"))
+    adata$stop[stop.in.start] <- adata$stop[stop.in.start]-1 
+  }  
+ 
+  return(adata) 
+}
+
+# concatenate peaks
+# peaks: the dataframe to have all the peaks
+# chr.file: the file to keep the length of each chromosome 
+# chr files should come from the species that the data is from
+#concatenate.chr <- function(peaks, chr.size){
+
+ # chr.size <- read.table(chr.file)
+#  chr.o <- order(chr.size[,1])
+#  chr.size <- chr.size[chr.o,]
+#
+#  chr.shift <- cumsum(c(0, chr.size[-nrow(chr.size),2]))
+#  chr.size.cum <- data.frame(chr=chr.size[,1], shift=chr.shift)  
+#
+#  for(i in 1:nrow(chr.size)){
+#    is.in <- as.character(peaks$chr) == as.character(chr.size.cum$chr[i])
+#    if(sum(is.in)>0){
+#      peaks[is.in,]$start <- peaks[is.in,]$start + chr.size.cum$shift[i]
+#      peaks[is.in,]$stop <- peaks[is.in,]$stop + chr.size.cum$shift[i]
+#    }
+#  }
+#
+#  invisible(peaks)
+#}
+
+
+
+
+# concatenate peaks
+# peaks: the dataframe to have all the peaks
+# chr.file: the file to keep the length of each chromosome 
+# chr files should come from the species that the data is from
+concatenate.chr <- function(peaks, chr.size){
+
+ # chr.size <- read.table(chr.file)
+  chr.o <- order(chr.size[,1])
+  chr.size <- chr.size[chr.o,]
+
+  chr.shift <- cumsum(c(0, chr.size[-nrow(chr.size),2]))
+  chr.size.cum <- data.frame(chr=chr.size[,1], shift=chr.shift)  
+
+  peaks$start.ori <- peaks$start
+  peaks$stop.ori <- peaks$stop
+  
+  for(i in 1:nrow(chr.size)){
+    is.in <- as.character(peaks$chr) == as.character(chr.size.cum$chr[i])
+    if(sum(is.in)>0){
+      peaks[is.in,]$start <- peaks[is.in,]$start + chr.size.cum$shift[i]
+      peaks[is.in,]$stop <- peaks[is.in,]$stop + chr.size.cum$shift[i]
+    }
+  }
+
+  invisible(peaks)
+}
+
+
+deconcatenate.chr <- function(peaks, chr.size){
+
+  chr.o <- order(chr.size[,1])
+  chr.size <- chr.size[chr.o,]
+
+  chr.shift <- cumsum(c(0, chr.size[-nrow(chr.size),2]))
+  chr.size.cum <- data.frame(chr=chr.size[,1], shift=chr.shift)  
+
+  peaks$chr <- rep(NA, nrow(peaks))
+  
+  for(i in 1:(nrow(chr.size.cum)-1)){
+    is.in <- peaks$start > chr.size.cum[i,2] & peaks$start <= chr.size.cum[i+1, 2]
+    if(sum(is.in)>0){
+      peaks[is.in,]$start <- peaks[is.in,]$start - chr.size.cum[i,2]
+      peaks[is.in,]$stop <- peaks[is.in,]$stop - chr.size.cum[i,2]+1    
+      peaks[is.in,]$chr <- chr.size[i,1]
+    }
+  }
+
+  if(i == nrow(chr.size.cum)){
+    is.in <- peaks$start > chr.size.cum[i, 2]
+    if(sum(is.in)>0){
+      peaks[is.in,]$start <- peaks[is.in,]$start - chr.size.cum[i,2]
+      peaks[is.in,]$stop <- peaks[is.in,]$stop - chr.size.cum[i,2]+1    
+      peaks[is.in,]$chr <- chr.size[i,1]
+    }
+  }
+  
+  invisible(peaks)
+}
+
+################ preprocessing peak calling output
+
+
+# 
+# read two calling results and sort by peak starting locations, 
+# then find overlap between peaks
+# INPUT:
+#   rep1: the 1st replicate
+#   rep2: the 2nd replicate
+# OUTPUT:
+#   id1, id2: the labels for the identified peaks on the replicates
+find.overlap <- function(rep1, rep2){
+
+  o1 <- order(rep1$start)
+  rep1 <- rep1[o1,]
+    
+  o2 <- order(rep2$start)
+  rep2 <- rep2[o2,]
+
+  n1 <- length(o1)
+  n2 <- length(o2)
+  
+  # assign common ID to peaks
+  id1 <- rep(0, n1) # ID assigned on rep1
+  id2 <- rep(0, n2) # ID assigned on rep2
+  id <- 1 # keep track common id's
+  
+  # check if two replicates overlap with each other
+  i <- 1
+  j <- 1
+
+  while(i <= n1|| j <= n2){
+
+    # && (id1[n1] ==0 || id2[n2] ==0)
+    
+    # if one list runs out
+    if(i > n1 && j < n2){
+      
+      j <- j+1
+      id2[j] <- id
+      id <- id +1
+      next
+    } else{
+      if(j > n2 && i < n1){
+        i <- i+1        
+        id1[i] <- id
+        id <- id +1
+        next
+      } else {
+        if(i >= n1 && j >=n2)
+          break
+      }
+    }
+
+    # if not overlap
+
+    if(!(rep1$start[i] <= rep2$stop[j] && rep2$start[j] <= rep1$stop[i])){
+
+      # at the start of loop, when both are not assigned an ID
+      # the one locates in front is assigned first
+      if(id1[i] ==0 && id2[j]==0){
+        if(rep1$stop[i] < rep2$stop[j]){
+          id1[i] <- id
+        } else {
+          id2[j] <- id
+        }
+      } else { # in the middle of the loop, when one is already assigned
+      # The one that has not assigned gets assigned
+      #  if(id1[i] ==0){ # id1[i] is not assigned
+      #    id1[i] <- id
+      #  } else { # id2[i] is not assigned
+      #    id2[j] <- id 
+      #  }
+
+        # order the id according to location
+        if(rep1$stop[i] <= rep2$stop[j]){
+          id1[i] <- max(id2[j], id1[i])
+          id2[j] <- id  
+        } else {
+          if(rep1$stop[i] > rep2$stop[j]){
+            id2[j] <- max(id1[i], id2[j])
+            id1[i] <- id
+          }
+        }
+          
+      }
+      
+      id <- id +1
+      
+    } else { # if overlap
+    
+      if(id1[i] == 0 && id2[j] == 0){ # not assign label yet
+        id1[i] <- id 
+        id2[j] <- id
+        id <- id +1
+      } else { # one peak is already assigned label, the other is 0
+        
+        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
+        id2[j] <- id1[i] # syncronize the labels        
+      }
+      
+    }
+    
+    if(rep1$stop[i] < rep2$stop[j]){
+      i <- i+1
+    } else {
+      j <- j+1
+    } 
+    
+  }
+
+  invisible(list(id1=id1, id2=id2))
+  
+}
+
+# Impute the missing significant value for the peaks called only on one replicate.
+# value 
+# INPUT:
+#   rep1, rep2: the two peak calling output 
+#   id1, id2: the IDs assigned by function find.overlap, vectors
+#        If id1[i]==id2[j], peak i on rep1 overlaps with peak j on rep2
+#   p.value.impute: the significant value to impute for the missing peaks 
+# OUTPUT:   
+#   rep1, rep2: peaks ordered by the start locations with imputed peaks
+#   id1, id2: the IDs with imputed peaks
+fill.missing.peaks <- function(rep1, rep2, id1, id2, p.value.impute){
+
+#   rep1 <- data.frame(chr=rep1$chr, start=rep1$start, stop=rep1$stop, sig.value=rep1$sig.value)
+#   rep2 <- data.frame(chr=rep2$chr, start=rep2$start, stop=rep2$stop, sig.value=rep2$sig.value)   
+   
+   o1 <- order(rep1$start)
+   rep1 <- rep1[o1,]
+    
+   o2 <- order(rep2$start)
+   rep2 <- rep2[o2,]  
+     
+   entry.in1.not2 <- !is.element(id1, id2)
+   entry.in2.not1 <- !is.element(id2, id1)
+
+   if(sum(entry.in1.not2) > 0){
+   
+     temp1 <- rep1[entry.in1.not2, ]
+
+     # impute sig.value
+     temp1$sig.value <- p.value.impute
+     temp1$signal.value <- p.value.impute
+     temp1$p.value <- p.value.impute
+     temp1$q.value <- p.value.impute
+     
+     rep2.filled <- rbind(rep2, temp1)
+     id2.filled <- c(id2, id1[entry.in1.not2])
+   } else {
+     id2.filled <- id2
+     rep2.filled <- rep2
+   }
+
+   if(sum(entry.in2.not1) > 0){
+
+     temp2 <- rep2[entry.in2.not1, ]
+
+     # fill in p.values to 1
+     temp2$sig.value <- p.value.impute
+     temp2$signal.value <- p.value.impute
+     temp2$p.value <- p.value.impute
+     temp2$q.value <- p.value.impute
+   
+
+     # append to the end
+     rep1.filled <- rbind(rep1, temp2)
+
+     id1.filled <- c(id1, id2[entry.in2.not1])
+   } else {
+     id1.filled <- id1
+     rep1.filled <- rep1
+   }
+
+   # sort rep1 and rep2 by the same id
+   o1 <- order(id1.filled)
+   rep1.ordered <- rep1.filled[o1, ]
+
+   o2 <- order(id2.filled)
+   rep2.ordered <- rep2.filled[o2, ]   
+   
+   invisible(list(rep1=rep1.ordered, rep2=rep2.ordered,
+                  id1=id1.filled[o1], id2=id2.filled[o2]))
+ }
+
+# Merge peaks with same ID on the same replicates 
+# (They are generated if two peaks on rep1 map to the same peak on rep2)
+# need peak.list have 3 columns: start, stop and sig.value 
+merge.peaks <- function(peak.list, id){
+
+  i <- 1
+  j <- 1
+  dup.index <- c()
+  sig.value <- c()
+  start.new <- c()
+  stop.new <- c()
+  id.new <- c()
+
+  # original data
+  chr <- c()
+  start.ori <- c()
+  stop.ori <- c()
+  
+  signal.value <- c()
+  p.value <- c()
+  q.value <- c()
+
+  while(i < length(id)){
+    
+    if(id[i] == id[i+1]){
+      dup.index <- c(dup.index, i, i+1) # push on dup.index
+    } else {
+      if(length(dup.index)>0){ # pop from dup.index
+        sig.value[j] <- mean(peak.list$sig.value[unique(dup.index)]) # mean of -log(pvalue)
+        start.new[j] <- peak.list$start[min(dup.index)]
+        stop.new[j] <- peak.list$stop[max(dup.index)]
+        id.new[j] <- id[max(dup.index)]
+        
+        signal.value[j] <- mean(peak.list$signal.value[unique(dup.index)]) # mean of -log(pvalue)
+        p.value[j] <- mean(peak.list$p.value[unique(dup.index)]) # mean of -log(pvalue)
+        q.value[j] <- mean(peak.list$q.value[unique(dup.index)]) # mean of -log(pvalue)
+
+        chr[j] <- as.character(peak.list$chr[min(dup.index)])
+        start.ori[j] <- peak.list$start.ori[min(dup.index)]
+        stop.ori[j] <- peak.list$stop.ori[max(dup.index)]
+        
+        dup.index <- c()
+      } else { # nothing to pop
+        sig.value[j] <- peak.list$sig.value[i]
+        start.new[j] <- peak.list$start[i]
+        stop.new[j] <- peak.list$stop[i]
+        id.new[j] <- id[i]
+
+        signal.value[j] <- peak.list$signal.value[i] 
+        p.value[j] <- peak.list$p.value[i] 
+        q.value[j] <- peak.list$q.value[i] 
+
+        chr[j] <- as.character(peak.list$chr[i])
+        start.ori[j] <- peak.list$start.ori[i]
+        stop.ori[j] <- peak.list$stop.ori[i]
+        
+      }
+      j <- j+1
+    }
+    i <- i+1
+  }
+
+  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)
+  invisible(data.new)
+}
+
+
+
+
+
+# a wrap function to fill in missing peaks, merge peaks and impute significant values
+# out1 and out2 are two peak calling outputs
+pair.peaks <- function(out1, out2, p.value.impute=0){
+
+  aa <- find.overlap(out1, out2)
+  bb <- fill.missing.peaks(out1, out2, aa$id1, aa$id2, p.value.impute=0)
+
+  cc1 <- merge.peaks(bb$rep1, bb$id1)
+  cc2 <- merge.peaks(bb$rep2, bb$id2)
+
+  invisible(list(merge1=cc1, merge2=cc2))
+}
+
+
+
+# overlap.ratio is a parameter to define the percentage of overlap
+# if overlap.ratio =0, 1 basepair overlap is counted as overlap
+# if overlap.ratio between 0 and 1, it is the minimum proportion of
+# overlap required to be called as a match
+# it is computed as the overlap part/min(peak1.length, peak2.length)
+pair.peaks.filter <- function(out1, out2, p.value.impute=0, overlap.ratio=0){
+
+  aa <- find.overlap(out1, out2)
+  bb <- fill.missing.peaks(out1, out2, aa$id1, aa$id2, p.value.impute=0)
+
+  cc1 <- merge.peaks(bb$rep1, bb$id1)
+  cc2 <- merge.peaks(bb$rep2, bb$id2)
+
+  frag12 <- cbind(cc1$start, cc1$stop, cc2$start, cc2$stop)
+  
+  frag.ratio <- apply(frag12, 1, overlap.middle)
+
+  frag.ratio[cc1$sig.value==p.value.impute | cc2$sig.value==p.value.impute] <- 0
+
+  cc1$frag.ratio <- frag.ratio
+  cc2$frag.ratio <- frag.ratio
+
+  merge1 <- cc1[cc1$frag.ratio >= overlap.ratio,]
+  merge2 <- cc2[cc2$frag.ratio >= overlap.ratio,]
+  
+  invisible(list(merge1=merge1, merge2=merge2))
+}
+
+# x[1], x[2] are the start and end of the first fragment
+# and x[3] and x[4] are the start and end of the 2nd fragment 
+# If there are two fragments, we can find the overlap by ordering the
+# start and stop of all the ends and find the difference between the middle two
+overlap.middle  <- function(x){
+
+  x.o <- x[order(x)]
+  f1 <- x[2]-x[1]
+  f2 <- x[4]-x[3]
+  
+  f.overlap <- abs(x.o[3]-x.o[2])
+  f.overlap.ratio <- f.overlap/min(f1, f2)
+
+  return(f.overlap.ratio)
+}
+
+
+
+#######
+####### compute correspondence profile
+#######
+
+# compute upper rank intersection for one t
+# tv: the upper percentile
+# x is sorted by the order of paired variable
+comp.uri <- function(tv, x){
+  n <- length(x)
+  qt <- quantile(x, prob=1-tv[1]) # tv[1] is t
+#  sum(x[1:ceiling(n*tv[2])] >= qt)/n/tv[2]- tv[1]*tv[2] #tv[2] is v
+  sum(x[1:ceiling(n*tv[2])] >= qt)/n
+
+}
+
+# compute the correspondence profile
+# tt, vv: vector between (0, 1) for percentages
+get.uri.2d <- function(x1, x2, tt, vv, spline.df=NULL){
+
+  o <- order(x1, x2, decreasing=T)
+  
+  # sort x2 by the order of x1
+  x2.ordered <- x2[o]
+  
+  tv <- cbind(tt, vv)
+  ntotal <- length(x1) # number of peaks    
+
+  uri <- apply(tv, 1, comp.uri, x=x2.ordered)
+
+  # compute the derivative of URI vs t using small bins
+  uri.binned <- uri[seq(1, length(uri), by=4)]
+  tt.binned <- tt[seq(1, length(uri), by=4)]
+  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)])
+
+  # smooth uri using spline
+  # first find where the jump is and don't fit the jump
+  # this is the index on the left
+  # jump.left.old  <- which.max(uri[-1]-uri[-length(uri)])
+  short.list.length <- min(sum(x1>0)/length(x1), sum(x2>0)/length(x2))
+
+  if(short.list.length < max(tt)){
+    jump.left <- which(tt>short.list.length)[1]-1
+  } else {
+    jump.left <- which.max(tt)
+  }
+
+#  reversed.index <- seq(length(tt), 1, by=-1)
+#  nequal <- sum(uri[reversed.index]== tt[reversed.index])
+#  temp  <- which(uri[reversed.index]== tt[reversed.index])[nequal]
+#  jump.left <- length(tt)-temp
+ 
+  if(jump.left < 6){
+   jump.left <- length(tt)
+  }
+    
+ 
+  if(is.null(spline.df))
+    uri.spl <- smooth.spline(tt[1:jump.left], uri[1:jump.left], df=6.4)
+  else{
+    uri.spl <- smooth.spline(tt[1:jump.left], uri[1:jump.left], df=spline.df)
+  }
+  # predict the first derivative
+  uri.der <- predict(uri.spl, tt[1:jump.left], deriv=1)
+
+  invisible(list(tv=tv, uri=uri, 
+                 uri.slope=uri.slope, t.binned=tt.binned[2:length(uri.binned)], 
+                 uri.spl=uri.spl, uri.der=uri.der, jump.left=jump.left,
+                 ntotal=ntotal))
+ }
+
+
+# change the scale of uri from based on t (percentage) to n (number of peaks or basepairs)
+# this is for plotting multiple pairwise URI's on the same plot 
+scale.t2n <- function(uri){
+
+  ntotal <- uri$ntotal
+  tv <- uri$tv*uri$ntotal
+  uri.uri <- uri$uri*uri$ntotal
+  jump.left <- uri$jump.left
+  uri.spl <- uri$uri.spl
+  uri.spl$x <- uri$uri.spl$x*uri$ntotal 
+  uri.spl$y <- uri$uri.spl$y*uri$ntotal
+
+  t.binned <- uri$t.binned*uri$ntotal
+  uri.slope <- uri$uri.slope
+  uri.der <- uri$uri.der
+  uri.der$x <- uri$uri.der$x*uri$ntotal
+  uri.der$y <- uri$uri.der$y
+
+  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)
+  return(uri.n)
+} 
+
+
+
+
+# a wrapper for running URI for peaks from peak calling results
+# both data1 and data2 are calling results in narrowpeak format
+compute.pair.uri <- function(data.1, data.2, sig.value1="signal.value", sig.value2="signal.value", spline.df=NULL, overlap.ratio=0){
+
+  tt <- seq(0.01, 1, by=0.01)
+  vv <- tt
+
+  if(sig.value1=="signal.value"){
+    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)
+  } else {
+    if(sig.value1=="p.value"){ 
+      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)
+    } else {
+      if(sig.value1=="q.value"){
+        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)
+      }
+    }
+  }
+
+  if(sig.value2=="signal.value"){
+    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)
+  } else {
+    if(sig.value2=="p.value"){ 
+      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)
+    } else {
+      if(sig.value2=="q.value"){
+        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)
+      }
+    }
+  }
+
+  ### by peaks
+  # data12.enrich <- pair.peaks(data.1.enrich, data.2.enrich)
+  data12.enrich <- pair.peaks.filter(data.1.enrich, data.2.enrich, p.value.impute=0, overlap.ratio)
+  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)
+  uri.n <- scale.t2n(uri)
+
+  return(list(uri=uri, uri.n=uri.n, data12.enrich=data12.enrich, sig.value1=sig.value1, sig.value2=sig.value2))
+
+
+}
+
+
+
+# compute uri for matched sample
+get.uri.matched <- function(data12, df=10){
+
+  tt <- seq(0.01, 1, by=0.01)
+  vv <- tt
+  uri <- get.uri.2d(data12$sample1$sig.value, data12$sample2$sig.value, tt, vv, spline.df=df)
+
+  # change scale from t to n
+  uri.n <- scale.t2n(uri)
+
+  return(list(uri=uri, uri.n=uri.n))
+  
+}
+
+# map.uv is a pair of significant values corresponding to specified consistency FDR
+# assuming values in map.uv and qvalue are linearly related
+# data.set is the original data set
+# sig.value is the name of the significant value in map.uv, say enrichment
+# nominal.value is the one we want to map to, say q-value
+# 
+map.sig.value <- function(data.set, map.uv, nominal.value){
+
+  index.nominal <- which(names(data.set$merge1)==nominal.value)
+  nentry <- nrow(map.uv)  
+  map.nominal <- rbind(map.uv[, c("sig.value1", "sig.value2")])
+
+  for(i in 1:nentry){
+
+    map.nominal[i, "sig.value1"] <- data.set$merge1[unique(which.min(abs(data.set$merge1$sig.value-map.uv[i, "sig.value1"]))), index.nominal]
+    map.nominal[i, "sig.value2"] <- data.set$merge2[unique(which.min(abs(data.set$merge2$sig.value-map.uv[i, "sig.value2"]))), index.nominal]
+  }
+
+  invisible(map.nominal)
+}
+
+
+############### plot correspondence profile
+
+# plot multiple comparison wrt one template
+# uri.list contains the total number of peaks
+# plot.missing=F: not plot the missing points on the right 
+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){
+
+  if(is.null(col.txt))
+    col.txt <- c("black", "red", "purple", "green", "blue", "cyan", "magenta", "orange", "grey")
+
+  n <- length(uri.n.list)  
+
+  ntotal <- c()
+  for(i in 1:n)
+    ntotal[i] <- uri.n.list[[i]]$ntotal
+
+  jump.left <- c()
+  jump.left.der <- c()
+  ncommon <- c()
+  for(i in 1:n){
+#    jump.left[i]  <- which.max(uri.n.list[[i]]$uri[-1]-uri.n.list[[i]]$uri[-length(uri.n.list[[i]]$uri)])
+#    if(jump.left[i] < 6)
+#      jump.left[i] <- length(uri.n.list[[i]]$uri)
+
+##  reversed.index <- seq(length(uri.n.list[[i]]$tv[,1]), 1, by=-1)
+##  nequal <- sum(uri.n.list[[i]]$uri[reversed.index]== uri.n.list[[i]]$tv[reversed.index,1])
+##  temp  <- which(uri.n.list[[i]]$uri[reversed.index]== uri.n.list[[i]]$tv[reversed.index,1])[nequal]
+##  jump.left[i] <- length(uri.n.list[[i]]$tv[,1])-temp
+##print(uri.n.list[[i]]$uri)
+##print(uri.n.list[[i]]$tv[,1])
+##   jump.left[i] <- uri.n.list[[i]]$jump.left
+
+#    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)])
+
+    jump.left[i] <- uri.n.list[[i]]$jump.left
+    jump.left.der[i] <- jump.left[i]
+    ncommon[i] <- uri.n.list[[i]]$tv[jump.left[i],1]
+  }
+
+
+  if(plot.missing){
+    max.peak <- max(ntotal)
+  } else {
+    max.peak <- max(ncommon)*1.05
+  }
+
+  if(!is.null(file.name)){
+    postscript(paste(plot.dir, "uri.", file.name, sep=""))
+    par(mfrow=c(1,1), mar=c(5,5,4,2))
+  }
+
+  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)
+
+  for(i in 1:n){
+
+    if(plot.missing){ 
+      points(uri.n.list[[i]]$tv[,1], uri.n.list[[i]]$uri, col=col.txt[i+col.start], cex=0.5 )
+    } else {
+      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)
+    }
+    lines(uri.n.list[[i]]$uri.spl, col=col.txt[i+col.start], lwd=4)
+  }
+  abline(coef=c(0,1), lty=3)
+  legend(0, max.peak, legend=legend.txt, col=col.txt[(col.start+1):length(col.txt)], lty=1, lwd=3, cex=2)
+  if(!is.null(title))
+    title(title.txt)
+
+  if(!is.null(file.name)){
+    dev.off()
+  }
+
+  if(!is.null(file.name)){
+    postscript(paste(plot.dir, "duri.", file.name, sep=""))
+    par(mfrow=c(1,1), mar=c(5,5,4,2))
+  } 
+  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)
+
+  for(i in 1:n){
+#    if(plot.missing){ 
+#      points(uri.n.list[[i]]$t.binned, uri.n.list[[i]]$uri.slope, col=col.txt[i+col.start], cex=0.5)
+#    } else {
+#      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)
+#    }
+    lines(uri.n.list[[i]]$uri.der, col=col.txt[i+col.start], lwd=4)
+  }
+  abline(h=1, lty=3)
+  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)
+
+  if(!is.null(title))
+    title(title.txt)
+
+  if(!is.null(file.name)){
+    dev.off()
+  }
+  
+}
+
+
+
+#######################
+####################### copula fitting for matched peaks
+#######################
+
+# estimation from mixed copula model 
+
+# 4-5-09
+# A nonparametric estimation of mixed copula model
+
+
+# updated
+
+# c1, c2, f1, f2, g1, g2 are vectors
+# c1*f1*g1 and c2*f2*g2 are copula densities for the two components
+# xd1 and yd1 are the values of marginals for the first component
+# xd2 and yd2 are the values of marginals for the 2nd component
+#
+# ez is the prob for being in the consistent group
+get.ez <- function(p, c1, c2, xd1, yd1, xd2, yd2){
+
+  return(p*c1*xd1*yd1/(p*c1*xd1*yd1 + (1-p)*c2*xd2*yd2))
+}
+
+# checked
+
+# this is C_12 not the copula density function c=C_12 * f1* f2
+# since nonparametric estimation is used here for f1 and f2, which
+# are constant throughout the iterations, we don't need them for optimization
+# 
+# bivariate gaussian copula function
+# t and s are vectors of same length, both are percentiles 
+# return a vector
+gaussian.cop.den <- function(t, s, rho){
+
+  A <- qnorm(t)^2 + qnorm(s)^2
+  B <- qnorm(t)*qnorm(s)
+
+  loglik <-  -log(1-rho^2)/2 - rho/(2*(1-rho^2))*(rho*A-2*B)
+
+  return(exp(loglik))
+}
+
+clayton.cop.den <- function(t, s, rho){
+
+  if(rho > 0)
+    return(exp(log(rho+1)-(rho+1)*(log(t)+log(s))-(2+1/rho)*log(t^(-rho) + s^(-rho)-1)))
+
+  if(rho==0)
+    return(1)
+
+  if(rho<0)
+    stop("Incorrect Clayton copula coefficient")
+  
+}
+
+
+# checked
+# estimate rho from Gaussian copula
+mle.gaussian.copula <- function(t, s, e.z){
+
+  # reparameterize to bound from rho=+-1
+  l.c <- function(rho, t, s, e.z){
+#    cat("rho=", rho, "\n")
+    sum(e.z*log(gaussian.cop.den(t, s, rho)))}
+
+  rho.max <- optimize(f=l.c, c(-0.998, 0.998), maximum=T, tol=0.00001, t=t, s=s, e.z=e.z)
+
+#print(rho.max$m)
+
+#cat("cor=", cor(qnorm(t)*e.z, qnorm(s)*e.z), "\t", "rho.max=", rho.max$m, "\n")
+#  return(sign(rho.max$m)/(1+rho.max$m))
+  return(rho.max$m)
+}
+
+
+# estimate mle from Clayton copula, 
+mle.clayton.copula <- function(t, s, e.z){
+
+  l.c <- function(rho, t, s, e.z){
+    lc <- sum(e.z*log(clayton.cop.den(t, s, rho)))
+#    cat("rho=", rho, "\t", "l.c=", lc, "\n")
+    return(lc)
+  }
+
+  rho.max <- optimize(f=l.c, c(0.1, 20), maximum=T, tol=0.00001, t=t, s=s, e.z=e.z)
+
+  return(rho.max$m)
+}
+
+
+
+# updated
+# mixture likelihood of two gaussian copula
+# nonparametric and ranked transformed
+loglik.2gaussian.copula <- function(x, y, p, rho1, rho2, x.mar, y.mar){
+ 
+  px.1 <- get.pdf.cdf(x, x.mar$f1)
+  px.2 <- get.pdf.cdf(x, x.mar$f2)
+  py.1 <- get.pdf.cdf(y, y.mar$f1)
+  py.2 <- get.pdf.cdf(y, y.mar$f2)
+
+  c1 <- gaussian.cop.den(px.1$cdf, py.1$cdf, rho1)
+  c2 <- gaussian.cop.den(px.2$cdf, py.2$cdf, rho2)
+
+  sum(log(p*c1*px.1$pdf*py.1$pdf + (1-p)*c2*px.2$pdf*py.2$pdf))
+}
+
+loglik.2copula <- function(x, y, p, rho1, rho2, x.mar, y.mar, copula.txt){
+
+  px.1 <- pdf.cdf$px.1
+  px.2 <- pdf.cdf$px.2
+  py.1 <- pdf.cdf$py.1
+  py.2 <- pdf.cdf$py.2
+
+  if(copula.txt=="gaussian"){
+    c1 <- gaussian.cop.den(px.1$cdf, py.1$cdf, rho1)
+    c2 <- gaussian.cop.den(px.2$cdf, py.2$cdf, rho2)
+  } else {
+    if(copula.txt=="clayton"){
+      c1 <- clayton.cop.den(px.1$cdf, py.1$cdf, rho1)
+      c2 <- clayton.cop.den(px.2$cdf, py.2$cdf, rho2)
+    }
+  }  
+  sum(log(p*c1*px.1$pdf*py.1$pdf + (1-p)*c2*px.2$pdf*py.2$pdf))
+}
+
+
+# estimate the marginals of each component using histogram estimator in EM
+# return the density, breaks, and cdf of the histogram estimator 
+est.mar.hist <- function(x, e.z, breaks){
+
+  binwidth <- c()
+  nbin <- length(breaks)-1
+  nx <- length(x) 
+
+  # the histogram
+  x1.pdf <- c()
+  x2.pdf <- c()
+  x1.cdf <- c()
+  x2.cdf <- c()
+
+  # the pdf for each point
+  x1.pdf.value <- rep(NA, nx)
+  x2.pdf.value <- rep(NA, nx)
+
+  x1.cdf.value <- rep(NA, nx)
+  x2.cdf.value <- rep(NA, nx) 
+
+  for(i in 1:nbin){
+
+    binwidth[i] <- breaks[i+1] - breaks[i]
+    if(i < nbin)
+      in.bin <- x>= breaks[i] & x < breaks[i+1]
+    else    # last bin
+      in.bin <- x>= breaks[i] & x <=breaks[i+1]
+
+    # each bin add one observation to avoid empty bins
+    # multiple (nx+nbin)/(nx+nbin+1) to avoid blowup when looking up for
+    # quantiles 
+    x1.pdf[i] <- (sum(e.z[in.bin])+1)/(sum(e.z)+nbin)/binwidth[i]*(nx+nbin)/(nx+nbin+1)        
+    x2.pdf[i] <- (sum(1-e.z[in.bin])+1)/(sum(1-e.z)+nbin)/binwidth[i]*(nx+nbin)/(nx+nbin+1) 
+
+
+#    x1.pdf[i] <- sum(e.z[in.bin])/sum(e.z)/binwidth[i]*nx/(nx+1)        
+#    x2.pdf[i] <- sum(1-e.z[in.bin])/sum(1-e.z)/binwidth[i]*nx/(nx+1) 
+    
+# treat each bin as a value for a discrete variable    
+#    x1.cdf[i] <- sum(x1.pdf[1:i]*binwidth[1:i])
+#    x2.cdf[i] <- sum(x2.pdf[1:i]*binwidth[1:i])
+
+
+    # cumulative density before reaching i
+    if(i>1){
+      x1.cdf[i] <- sum(x1.pdf[1:(i-1)]*binwidth[1:(i-1)])
+      x2.cdf[i] <- sum(x2.pdf[1:(i-1)]*binwidth[1:(i-1)])    
+    } else{
+      x1.cdf[i] <- 0
+      x2.cdf[i] <- 0
+    }
+
+    # make a vector of nx to store the values of pdf and cdf for each x
+    # this will speed up the computation dramatically
+    x1.pdf.value[in.bin] <- x1.pdf[i]
+    x2.pdf.value[in.bin] <- x2.pdf[i]
+
+    x1.cdf.value[in.bin] <- x1.cdf[i] + x1.pdf[i]*(x[in.bin]-breaks[i])
+    x2.cdf.value[in.bin] <- x2.cdf[i] + x2.pdf[i]*(x[in.bin]-breaks[i])      
+  }
+
+#  x1.cdf <- cumsum(x1.pdf*binwidth)
+#  x2.cdf <- cumsum(x2.pdf*binwidth)
+
+  f1 <-list(breaks=breaks, density=x1.pdf, cdf=x1.cdf)
+  f2 <-list(breaks=breaks, density=x2.pdf, cdf=x2.cdf)
+
+  f1.value <- list(pdf=x1.pdf.value, cdf=x1.cdf.value)
+  f2.value <- list(pdf=x2.pdf.value, cdf=x2.cdf.value)
+
+  return(list(f1=f1, f2=f2, f1.value=f1.value, f2.value=f2.value))
+}
+
+# estimate the marginal cdf from rank
+est.cdf.rank <- function(x, conf.z){
+
+  # add 1 to prevent blow up
+  x1.cdf <- rank(x[conf.z==1])/(length(x[conf.z==1])+1)
+
+  x2.cdf <- rank(x[conf.z==0])/(length(x[conf.z==0])+1)
+
+  return(list(cdf1=x1.cdf, cdf2=x2.cdf))
+}
+
+# df is a density function with fields: density, cdf and breaks, x is a scalar
+get.pdf <- function(x, df){
+
+  if(x < df$breaks[1])
+    cat("x is out of the range of df\n")
+
+  index <- which(df$breaks >= x)[1]
+
+  if(index==1)
+    index <- index +1
+  return(df$density[index-1])  
+}
+
+# get cdf from histgram estimator for a single value
+get.cdf <- function(x, df){
+
+  index <- which(df$breaks >= x)[1]
+  if(index==1)
+    index <- index +1
+  return(df$cdf[index-1])   
+}
+
+# df is a density function with fields: density, cdf and breaks
+get.pdf.cdf <- function(x.vec, df){
+
+  x.pdf <- sapply(x.vec, get.pdf, df=df)
+  x.cdf <- sapply(x.vec, get.cdf, df=df) 
+  return(list(cdf=x.cdf, pdf=x.pdf))
+}
+
+# E-step
+# x and y are the original observations or ranks
+# rho1 and rho2 are the parameters of each copula
+# f1, f2, g1, g2 are functions, each is a histogram 
+e.step.2gaussian <- function(x, y, p, rho1, rho2, x.mar, y.mar){
+
+  # get pdf and cdf of each component from functions in the corresponding component 
+  px.1 <- get.pdf.cdf(x, x.mar$f1)
+  px.2 <- get.pdf.cdf(x, x.mar$f2)
+  py.1 <- get.pdf.cdf(y, y.mar$f1)
+  py.2 <- get.pdf.cdf(y, y.mar$f2)
+
+  c1 <- gaussian.cop.den(px.1$cdf, py.1$cdf, rho1)
+  c2 <- gaussian.cop.den(px.2$cdf, py.2$cdf, rho2)
+  
+  return(get.ez(p, c1, c2, px.1$pdf, py.1$pdf, px.2$pdf, py.2$pdf))
+}
+
+# E-step
+# rho1 and rho2 are the parameters of each copula 
+e.step.2copula <- function(x, y, p, rho1, rho2, x.mar, y.mar, copula.txt){
+
+  # get pdf and cdf of each component from functions in the corresponding component 
+  px.1 <- get.pdf.cdf(x, x.mar$f1)
+  px.2 <- get.pdf.cdf(x, x.mar$f2)
+  py.1 <- get.pdf.cdf(y, y.mar$f1)
+  py.2 <- get.pdf.cdf(y, y.mar$f2)
+
+  if(copula.txt=="gaussian"){
+    c1 <- gaussian.cop.den(px.1$cdf, py.1$cdf, rho1)
+    c2 <- gaussian.cop.den(px.2$cdf, py.2$cdf, rho2)
+  } else {
+    if(copula.txt=="clayton"){
+      c1 <- clayton.cop.den(px.1$cdf, py.1$cdf, rho1)
+      c2 <- clayton.cop.den(px.2$cdf, py.2$cdf, rho2)
+    } 
+  }
+  return(get.ez(p, c1, c2, px.1$pdf, py.1$pdf, px.2$pdf, py.2$pdf))
+}
+
+
+
+
+# M-step
+m.step.2gaussian <- function(x, y, e.z, breaks){
+
+  # compute f1, f2, g1 and g2
+  x.mar <- est.mar.hist(x, e.z, breaks)
+  y.mar <- est.mar.hist(y, e.z, breaks)  
+
+  px.1 <- get.pdf.cdf(x, x.mar$f1)
+  px.2 <- get.pdf.cdf(x, x.mar$f2)
+  py.1 <- get.pdf.cdf(y, y.mar$f1)
+  py.2 <- get.pdf.cdf(y, y.mar$f2)
+
+  rho1 <- mle.gaussian.copula(px.1$cdf, py.1$cdf, e.z)  
+  rho2 <- mle.gaussian.copula(px.2$cdf, py.2$cdf, 1-e.z) 
+
+  p <- sum(e.z)/length(e.z) 
+
+  return(list(p=p, rho1=rho1, rho2=rho2, x.mar=x.mar, y.mar=y.mar))
+}
+
+m.step.2copula <- function(x, y, e.z, breaks, copula.txt){
+
+  # compute f1, f2, g1 and g2
+  x.mar <- est.mar.hist(x, e.z, breaks)
+  y.mar <- est.mar.hist(y, e.z, breaks)  
+
+  px.1 <- get.pdf.cdf(x, x.mar$f1)
+  px.2 <- get.pdf.cdf(x, x.mar$f2)
+  py.1 <- get.pdf.cdf(y, y.mar$f1)
+  py.2 <- get.pdf.cdf(y, y.mar$f2)
+
+  if(copula.txt=="gaussian"){
+    rho1 <- mle.gaussian.copula(px.1$cdf, py.1$cdf, e.z)  
+    rho2 <- mle.gaussian.copula(px.2$cdf, py.2$cdf, 1-e.z) 
+  } else {
+    if(copula.txt=="clayton"){
+      rho1 <- mle.clayton.copula(px.1$cdf, py.1$cdf, e.z)  
+      rho2 <- mle.clayton.copula(px.2$cdf, py.2$cdf, 1-e.z)      
+    }
+  }
+  
+  p <- sum(e.z)/length(e.z) 
+
+  return(list(p=p, rho1=rho1, rho2=rho2, x.mar=x.mar, y.mar=y.mar))
+}
+
+
+
+# E-step: pass values
+# x and y are the original observations or ranks
+# rho1 and rho2 are the parameters of each copula
+# f1, f2, g1, g2 are functions, each is a histogram 
+e.step.2gaussian.value <- function(x, y, p, rho1, rho2, pdf.cdf){
+
+  c1 <- gaussian.cop.den(pdf.cdf$px.1$cdf, pdf.cdf$py.1$cdf, rho1)
+  c2 <- gaussian.cop.den(pdf.cdf$px.2$cdf, pdf.cdf$py.2$cdf, rho2)
+  
+  e.z <- get.ez(p, c1, c2, pdf.cdf$px.1$pdf, pdf.cdf$py.1$pdf, 
+               pdf.cdf$px.2$pdf, pdf.cdf$py.2$pdf)
+  return(e.z)
+}
+
+
+e.step.2copula.value <- function(x, y, p, rho1, rho2, pdf.cdf, copula.txt){
+
+  if(copula.txt =="gaussian"){
+    c1 <- gaussian.cop.den(pdf.cdf$px.1$cdf, pdf.cdf$py.1$cdf, rho1)
+    c2 <- gaussian.cop.den(pdf.cdf$px.2$cdf, pdf.cdf$py.2$cdf, rho2)
+  } else {
+    if(copula.txt =="clayton"){
+      c1 <- clayton.cop.den(pdf.cdf$px.1$cdf, pdf.cdf$py.1$cdf, rho1)
+      c2 <- clayton.cop.den(pdf.cdf$px.2$cdf, pdf.cdf$py.2$cdf, rho2)      
+    }
+  }
+  
+  e.z <- get.ez(p, c1, c2, pdf.cdf$px.1$pdf, pdf.cdf$py.1$pdf, 
+               pdf.cdf$px.2$pdf, pdf.cdf$py.2$pdf)
+  return(e.z)
+}
+
+
+# M-step: pass values
+m.step.2gaussian.value <- function(x, y, e.z, breaks, fix.rho2){
+
+  # compute f1, f2, g1 and g2
+  x.mar <- est.mar.hist(x, e.z, breaks)
+  y.mar <- est.mar.hist(y, e.z, breaks)  
+
+#  px.1 <- get.pdf.cdf(x, x.mar$f1)
+#  px.2 <- get.pdf.cdf(x, x.mar$f2)
+#  py.1 <- get.pdf.cdf(y, y.mar$f1)
+#  py.2 <- get.pdf.cdf(y, y.mar$f2)
+
+  px.1 <- x.mar$f1.value
+  px.2 <- x.mar$f2.value
+  py.1 <- y.mar$f1.value
+  py.2 <- y.mar$f2.value
+
+  rho1 <- mle.gaussian.copula(px.1$cdf, py.1$cdf, e.z)  
+
+  if(!fix.rho2)
+    rho2 <- mle.gaussian.copula(px.2$cdf, py.2$cdf, 1-e.z) 
+  else
+    rho2 <- 0
+
+  p <- sum(e.z)/length(e.z) 
+
+  pdf.cdf <- list(px.1=px.1, px.2=px.2, py.1=py.1, py.2=py.2)
+
+  return(list(p=p, rho1=rho1, rho2=rho2, x.mar=x.mar, y.mar=y.mar,
+              pdf.cdf=pdf.cdf))
+}
+
+m.step.2gaussian.value2 <- function(x, y, e.z, breaks, fix.rho2, x.mar, y.mar){
+
+  # compute f1, f2, g1 and g2
+#  x.mar <- est.mar.hist(x, e.z, breaks)
+#  y.mar <- est.mar.hist(y, e.z, breaks)  
+
+#  px.1 <- get.pdf.cdf(x, x.mar$f1)
+#  px.2 <- get.pdf.cdf(x, x.mar$f2)
+#  py.1 <- get.pdf.cdf(y, y.mar$f1)
+#  py.2 <- get.pdf.cdf(y, y.mar$f2)
+
+  px.1 <- x.mar$f1.value
+  px.2 <- x.mar$f2.value
+  py.1 <- y.mar$f1.value
+  py.2 <- y.mar$f2.value
+
+  rho1 <- mle.gaussian.copula(px.1$cdf, py.1$cdf, e.z)  
+
+  if(!fix.rho2)
+    rho2 <- mle.gaussian.copula(px.2$cdf, py.2$cdf, 1-e.z) 
+  else
+    rho2 <- 0
+
+  p <- sum(e.z)/length(e.z) 
+
+  pdf.cdf <- list(px.1=px.1, px.2=px.2, py.1=py.1, py.2=py.2)
+
+  return(list(p=p, rho1=rho1, rho2=rho2, x.mar=x.mar, y.mar=y.mar,
+              pdf.cdf=pdf.cdf))
+}
+
+
+
+m.step.2copula.value <- function(x, y, e.z, breaks, fix.rho2, copula.txt){
+
+  # compute f1, f2, g1 and g2
+  x.mar <- est.mar.hist(x, e.z, breaks)
+  y.mar <- est.mar.hist(y, e.z, breaks)  
+
+#  px.1 <- get.pdf.cdf(x, x.mar$f1)
+#  px.2 <- get.pdf.cdf(x, x.mar$f2)
+#  py.1 <- get.pdf.cdf(y, y.mar$f1)
+#  py.2 <- get.pdf.cdf(y, y.mar$f2)
+
+  px.1 <- x.mar$f1.value
+  px.2 <- x.mar$f2.value
+  py.1 <- y.mar$f1.value
+  py.2 <- y.mar$f2.value
+
+  if(copula.txt=="gaussian"){
+    rho1 <- mle.gaussian.copula(px.1$cdf, py.1$cdf, e.z)  
+    
+    if(!fix.rho2)
+      rho2 <- mle.gaussian.copula(px.2$cdf, py.2$cdf, 1-e.z) 
+    else
+      rho2 <- 0
+  } else {
+
+    if(copula.txt=="clayton"){
+      rho1 <- mle.clayton.copula(px.1$cdf, py.1$cdf, e.z)  
+    
+      if(!fix.rho2)
+        rho2 <- mle.clayton.copula(px.2$cdf, py.2$cdf, 1-e.z) 
+      else
+        rho2 <- 0
+    }    
+  }
+    
+  p <- sum(e.z)/length(e.z) 
+
+  pdf.cdf <- list(px.1=px.1, px.2=px.2, py.1=py.1, py.2=py.2)
+
+  return(list(p=p, rho1=rho1, rho2=rho2, x.mar=x.mar, y.mar=y.mar,
+              pdf.cdf=pdf.cdf))
+}
+
+
+
+
+# updated
+# mixture likelihood of two gaussian copula
+# nonparametric and ranked transformed
+loglik.2gaussian.copula.value <- function(x, y, p, rho1, rho2, pdf.cdf){
+
+  px.1 <- pdf.cdf$px.1
+  px.2 <- pdf.cdf$px.2
+  py.1 <- pdf.cdf$py.1
+  py.2 <- pdf.cdf$py.2
+
+  c1 <- gaussian.cop.den(px.1$cdf, py.1$cdf, rho1)
+  c2 <- gaussian.cop.den(px.2$cdf, py.2$cdf, rho2)
+
+  sum(log(p*c1*px.1$pdf*py.1$pdf + (1-p)*c2*px.2$pdf*py.2$pdf))
+}
+
+
+
+# updated
+# mixture likelihood of two gaussian copula
+# nonparametric and ranked transformed
+loglik.2copula.value <- function(x, y, p, rho1, rho2, pdf.cdf, copula.txt){
+
+  px.1 <- pdf.cdf$px.1
+  px.2 <- pdf.cdf$px.2
+  py.1 <- pdf.cdf$py.1
+  py.2 <- pdf.cdf$py.2
+
+  if(copula.txt=="gaussian"){
+    c1 <- gaussian.cop.den(px.1$cdf, py.1$cdf, rho1)
+    c2 <- gaussian.cop.den(px.2$cdf, py.2$cdf, rho2)
+  } else {
+    if(copula.txt=="clayton"){
+      c1 <- clayton.cop.den(px.1$cdf, py.1$cdf, rho1)
+      c2 <- clayton.cop.den(px.2$cdf, py.2$cdf, rho2)
+    }
+  }
+
+  sum(log(p*c1*px.1$pdf*py.1$pdf + (1-p)*c2*px.2$pdf*py.2$pdf))
+}
+
+
+
+# EM for 2 Gaussian, speed up computation, unfinished
+
+em.2gaussian.quick <- function(x, y, p0, rho1.0, rho2.0, eps, fix.p=F, stoc=T, fix.rho2=T){
+
+  x <- rank(x, tie="random")
+  y <- rank(y, tie="random")
+
+#  x <- rank(x, tie="average")
+#  y <- rank(y, tie="average")
+
+  # nbin=20
+  xy.min <- min(x, y)
+  xy.max <- max(x, y)
+  binwidth <- (xy.max-xy.min)/50
+  breaks <- seq(xy.min-binwidth/100, xy.max+binwidth/100, by=(xy.max-xy.min+binwidth/50)/50)
+#  breaks <- seq(xy.min, xy.max, by=binwidth)
+  
+
+  # initiate marginals 
+  # initialization: first p0 data has 
+#  e.z <- e.step.2gaussian(x, y, p0, rho1.0, rho2.0, x0.mar, y0.mar) # this starting point assumes two components are overlapped
+
+  e.z <- c(rep(0.9, round(length(x)*p0)), rep(0.1, length(x)-round(length(x)*p0)))
+
+  if(!stoc)
+    para <- m.step.2gaussian.value(x, y, e.z, breaks, fix.rho2)
+  else 
+    para <- m.step.2gaussian.stoc.value(x, y, e.z, breaks, fix.rho2)
+
+
+  if(fix.p){
+    p <- p0
+  } else {
+    p <- para$p  
+  }
+
+  if(fix.rho2){
+    rho2 <- rho2.0
+  } else {
+    rho2 <- para$rho2
+  }
+
+#  rho1 <- 0.8
+  rho1 <- para$rho1
+
+  l0 <- loglik.2gaussian.copula.value(x, y, p, rho1, rho2, para$pdf.cdf)
+
+  loglik.trace <- c()
+  loglik.trace[1] <- l0
+#  loglik.trace[2] <- l1
+  to.run <- T
+
+  i <- 2
+
+  # this two lines to remove
+#  x.mar <- est.mar.hist(x, e.z, breaks)
+#  y.mar <- est.mar.hist(y, e.z, breaks)  
+  
+  while(to.run){
+
+    e.z <- e.step.2gaussian.value(x, y, p, rho1, rho2, para$pdf.cdf) 
+    if(!stoc)
+      para <- m.step.2gaussian.value(x, y, e.z, breaks, fix.rho2)
+    else
+      para <- m.step.2gaussian.stoc.value(x, y, e.z, breaks, fix.rho2)
+
+    # fix x.mar and y.mar : to remove
+#    if(!stoc)
+#      para <- m.step.2gaussian.value2(x, y, e.z, breaks, fix.rho2, x.mar, y.mar)
+#    else
+#      para <- m.step.2gaussian.stoc.value(x, y, e.z, breaks, fix.rho2)
+
+    
+    if(fix.p){
+      p <- p0
+    } else {
+      p <- para$p  
+    }
+
+    if(fix.rho2){
+      rho2 <- rho2.0
+    } else {
+      rho2 <- para$rho2
+    }
+
+#    rho1 <- 0.8
+    rho1 <- para$rho1
+
+  #  l0 <- l1
+    l1 <- loglik.2gaussian.copula.value(x, y, p, rho1, rho2, para$pdf.cdf)
+    loglik.trace[i] <- l1
+
+#cat("l1=", l1, "\n") 
+
+    # Aitken acceleration criterion
+    if(i > 2){
+      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])) 
+      to.run <- abs(l.inf - loglik.trace[i]) > eps 
+#cat("para=", "p=", para$p, " rho1=", rho1, " rho2=", rho2, "\n")
+#cat("l.inf=", l.inf, "\n")
+#cat(l.inf-loglik.trace[i], "\n")     
+    }
+
+    i <- i+1
+  }
+
+  bic <- -2*l1 + (2*(length(breaks)-1+1)+1-fix.p-fix.rho2)*log(length(x)) # parameters
+  return(list(para=list(p=para$p, rho1=rho1, rho2=rho2), 
+              loglik=l1, bic=bic, e.z=e.z, conf.z = para$conf.z, 
+              loglik.trace=loglik.trace, x.mar=para$x.mar, y.mar=para$y.mar,
+              breaks=breaks))
+}
+
+
+
+em.2copula.quick <- function(x, y, p0, rho1.0, rho2.0, eps, fix.p=F, stoc=T, fix.rho2=T, copula.txt, nbin=50){
+
+  x <- rank(x, tie="random")
+  y <- rank(y, tie="random")
+
+#  x <- rank(x, tie="first")
+#  y <- rank(y, tie="first")
+
+  # nbin=50
+  xy.min <- min(x, y)
+  xy.max <- max(x, y)
+  binwidth <- (xy.max-xy.min)/50
+  breaks <- seq(xy.min-binwidth/100, xy.max+binwidth/100, by=(xy.max-xy.min+binwidth/50)/nbin)  
+#  breaks <- seq(xy.min, xy.max, by=binwidth)
+  
+  # initiate marginals 
+  # initialization: first p0 data has 
+#  e.z <- e.step.2gaussian(x, y, p0, rho1.0, rho2.0, x0.mar, y0.mar) # this starting point assumes two components are overlapped
+
+  e.z <- c(rep(0.9, round(length(x)*p0)), rep(0.1, length(x)-round(length(x)*p0)))
+
+
+  if(!stoc)
+    para <- m.step.2copula.value(x, y, e.z, breaks, fix.rho2, copula.txt)
+  else 
+    para <- m.step.2copula.stoc.value(x, y, e.z, breaks, fix.rho2, copula.txt)
+
+  if(fix.p){
+    p <- p0
+  } else {
+    p <- para$p  
+  }
+
+  if(fix.rho2){
+    rho2 <- rho2.0
+  } else {
+    rho2 <- para$rho2
+  }
+
+  l0 <- loglik.2copula.value(x, y, p, para$rho1, rho2, para$pdf.cdf, copula.txt)
+
+  loglik.trace <- c()
+  loglik.trace[1] <- l0
+#  loglik.trace[2] <- l1
+  to.run <- T
+
+  i <- 2
+
+  while(to.run){
+
+    e.z <- e.step.2copula.value(x, y, p, para$rho1, rho2, para$pdf.cdf, copula.txt) 
+    if(!stoc)
+      para <- m.step.2copula.value(x, y, e.z, breaks, fix.rho2, copula.txt)
+    else
+      para <- m.step.2copula.stoc.value(x, y, e.z, breaks, fix.rho2, copula.txt)
+
+    if(fix.p){
+      p <- p0
+    } else {
+      p <- para$p  
+    }
+
+    if(fix.rho2){
+      rho2 <- rho2.0
+    } else {
+      rho2 <- para$rho2
+    }
+
+
+  #  l0 <- l1
+    l1 <- loglik.2copula.value(x, y, p, para$rho1, rho2, para$pdf.cdf, copula.txt)
+    loglik.trace[i] <- l1
+
+cat("l1=", l1, "\n") 
+
+    # Aitken acceleration criterion
+    if(i > 2){
+      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])) 
+      to.run <- abs(l.inf - loglik.trace[i]) > eps 
+cat("para=", "p=", para$p, " rho1=", para$rho1, " rho2=", rho2, "\n")
+#cat("l.inf=", l.inf, "\n")
+#cat(l.inf-loglik.trace[i], "\n")     
+    }
+
+    i <- i+1
+  }
+
+  bic <- -2*l1 + (2*(length(breaks)-1+1)+1-fix.p-fix.rho2)*log(length(x)) # parameters
+  return(list(para=list(p=para$p, rho1=para$rho1, rho2=rho2), 
+              loglik=l1, bic=bic, e.z=e.z, conf.z = para$conf.z, 
+              loglik.trace=loglik.trace, x.mar=para$x.mar, y.mar=para$y.mar,
+              breaks=breaks))
+}
+
+
+#######################
+####################### fit EM procedure for the matched peaks
+#######################
+
+# remove the unmatched ones
+#rm.unmatch <- function(sample1, sample2, p.value.impute=0){
+#
+#  sample1.prune <- sample1[sample1$sig.value > p.value.impute & sample2$sig.value > p.value.impute,]
+#  sample2.prune <- sample2[sample1$sig.value > p.value.impute & sample2$sig.value > p.value.impute,]
+# 
+#  invisible(list(sample1=sample1.prune$sig.value, sample2=sample2.prune$sig.value))
+#}
+
+
+# fit 2-component model
+#fit.em <- function(sample12, fix.rho2=T){
+#
+#  prune.sample <- rm.unmatch(sample12$merge1, sample12$merge2)
+#
+#  em.fit <- em.2gaussian.quick(-prune.sample$sample1, -prune.sample$sample2,
+# p0=0.5, rho1.0=0.7, rho2.0=0, eps=0.01, fix.p=F, stoc=F, fix.rho2)
+#
+#  invisible(list(em.fit=em.fit, data.pruned=prune.sample))
+#}
+
+
+rm.unmatch <- function(sample1, sample2, p.value.impute=0){
+
+  sample1.prune <- sample1[sample1$sig.value > p.value.impute & sample2$sig.value > p.value.impute,]
+  sample2.prune <- sample2[sample1$sig.value > p.value.impute & sample2$sig.value > p.value.impute,]
+ 
+  invisible(list(sample1=sample1.prune, sample2=sample2.prune))
+}
+
+
+# fit 2-component model
+fit.em <- function(sample12, fix.rho2=T){
+
+  prune.sample <- rm.unmatch(sample12$merge1, sample12$merge2)
+
+  em.fit <- em.2gaussian.quick(-prune.sample$sample1$sig.value, -prune.sample$sample2$sig.value,
+ p0=0.5, rho1.0=0.7, rho2.0=0, eps=0.01, fix.p=F, stoc=F, fix.rho2)
+
+  invisible(list(em.fit=em.fit, data.pruned=prune.sample))
+}
+
+
+
+fit.2copula.em <- function(sample12, fix.rho2=T, copula.txt){
+
+  prune.sample <- rm.unmatch(sample12$merge1, sample12$merge2)
+
+#  o <- order(prune.sample$sample1)
+#  n <- length(prune.sample$sample1)
+    
+#  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))))
+
+#  temp <- init.dist(f0, f1)
+  para <- list()
+  para$rho <- 0.6
+  para$p <- 0.3
+  para$mu <- 2.5
+  para$sigma <- 1
+##  para$mu <- -temp$mu
+##  para$sigma <- temp$sigma
+#cat("mu=", para$mu, "sigma=", para$sigma, "\n")
+  
+#  em.fit <- em.transform.1loop(-prune.sample$sample1, -prune.sample$sample2,
+  cat("EM is running")
+  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)
+
+  invisible(list(em.fit=em.fit, data.pruned=prune.sample))
+}
+
+
+
+
+# fit 1-component model
+fit.1.component <- function(data.pruned, breaks){
+
+#  gaussian.1 <- fit.gaussian.1(-data.pruned$sample1$sig.value, -data.pruned$sample2$sig.value, breaks)
+#  clayton.1 <- fit.clayton.1(-data.pruned$sample1$sig.value, -data.pruned$sample2$sig.value, breaks)
+
+  gaussian.1 <- fit.gaussian.1(-data.pruned$sample1, -data.pruned$sample2, breaks)
+  clayton.1 <- fit.clayton.1(-data.pruned$sample1, -data.pruned$sample2, breaks)
+
+  return(list(gaussian.1=gaussian.1, clayton.1=clayton.1))
+}
+
+
+
+#################
+# Fit a single component  
+#################
+
+# a single gaussian copula
+# if breaks=NULL, use empirical pdf, otherwise use histogram estimate
+fit.gaussian.1 <- function(x, y, breaks=NULL){
+
+  # rank transformed and compute the empirical cdf
+  t <- emp.mar.cdf.rank(x)
+  s <- emp.mar.cdf.rank(y)
+
+  mle.rho <- mle.gaussian.copula(t, s, rep(1, length(t)))
+
+  c1 <- gaussian.cop.den(t, s, mle.rho)
+cat("c1", sum(log(c1)), "\n")
+
+  if(is.null(breaks)){
+    f1 <- emp.mar.pdf.rank(t)
+    f2 <- emp.mar.pdf.rank(s)
+  } else {
+    x.mar <- est.mar.hist(rank(x), rep(1, length(x)), breaks)
+    y.mar <- est.mar.hist(rank(y), rep(1, length(y)), breaks)
+
+    f1 <- x.mar$f1.value$pdf  # only one component
+    f2 <- y.mar$f1.value$pdf
+  }
+
+
+cat("f1", sum(log(f1)), "\n")
+cat("f2", sum(log(f2)), "\n")
+
+  loglik <- sum(log(c1)+log(f1)+log(f2))
+
+  bic <- -2*loglik + log(length(t))*(1+length(breaks)-1)
+
+  return(list(rho=mle.rho, loglik=loglik, bic=bic))
+}
+
+
+# a single Clayton copula
+fit.clayton.1 <- function(x, y, breaks=NULL){
+
+  # rank transformed and compute the empirical cdf
+  t <- emp.mar.cdf.rank(x)
+  s <- emp.mar.cdf.rank(y)
+
+  mle.rho <- mle.clayton.copula(t, s, rep(1, length(t)))
+
+  c1 <- clayton.cop.den(t, s, mle.rho)
+
+  if(is.null(breaks)){
+    f1 <- emp.mar.pdf.rank(t)
+    f2 <- emp.mar.pdf.rank(s)
+  } else {
+    x.mar <- est.mar.hist(rank(x), rep(1, length(x)), breaks)
+    y.mar <- est.mar.hist(rank(y), rep(1, length(y)), breaks)
+
+    f1 <- x.mar$f1.value$pdf  # only one component
+    f2 <- y.mar$f1.value$pdf
+  }
+
+  loglik <- sum(log(c1)+log(f1)+log(f2))
+
+  bic <- -2*loglik + log(length(t))*(1+length(breaks)-1)
+
+  return(list(rho=mle.rho, tau=rho/(rho+2), loglik=loglik, bic=bic)) 
+}
+
+## obsolete function (01-06-2010)
+## compute the average posterior probability to belong to the random component
+## for peaks selected at different cutoffs 
+comp.uri.ez <- function(tt, u, v, e.z){
+
+   u.t <- quantile(u, prob=(1-tt))
+   v.t <- quantile(v, prob=(1-tt))
+
+ #  ez <- mean(e.z[u >= u.t & v >=u.t]) Is this wrong?
+   ez <- mean(e.z[u >= u.t & v >=v.t])
+
+   return(ez)
+}
+
+## obsolete function (01-06-2010)
+# compute the largest posterior error probability corresponding to
+# the square centered at the origin and spanned top tt% on both coordinates
+# so the consistent low rank ones are excluded
+# boundary.txt: either "max" or "min", if it is error prob, use "max"
+comp.ez.cutoff <- function(tt, u, v, e.z, boundary.txt){
+
+   u.t <- quantile(u, prob=(1-tt))
+   v.t <- quantile(v, prob=(1-tt))
+
+   if(boundary.txt == "max"){
+ #    ez.bound <- max(e.z[u >= u.t & v >=u.t])
+     ez.bound <- max(e.z[u >= u.t & v >=v.t])
+   } else {
+ #    ez.bound <- min(e.z[u >= u.t & v >=u.t])
+     ez.bound <- min(e.z[u >= u.t & v >=v.t])     
+   }
+
+   return(ez.bound)
+
+}
+
+# obsolete functions: 01-06-2010
+# compute the error rate
+# u.t and v.t are the quantiles
+# this one is used for the plots generated initially in the brief writeup  
+# and it was used for processing merged data in July before the IDR definition
+# is formalized
+# It does not implement the current definition of IDR
+get.ez.tt.old  <- function(em.fit, reverse=T, fdr.level=c(0.01, 0.05, 0.1)){
+
+  u <- em.fit$data.pruned$sample1
+  v <- em.fit$data.pruned$sample2
+
+  tt <- seq(0.01, 0.99, by=0.01)
+  if(reverse){ 
+    e.z <-  1-em.fit$em.fit$e.z # this is the error prob
+    uri.ez <- sapply(tt, comp.uri.ez, u=u, v=v, e.z=e.z)
+    ez.bound <- sapply(tt, comp.ez.cutoff, u=u, v=v, e.z=e.z, boundary.txt="max") 
+  } else {
+    e.z <-  em.fit$em.fit$e.z
+    uri.ez <- sapply(tt, comp.uri.ez, u=u, v=v, e.z=e.z)
+    ez.bound <- sapply(tt, comp.ez.cutoff, u=u, v=v, e.z=e.z, boundary.txt="min") 
+  }
+
+  u.t <- quantile(u, prob=(1-tt))
+  v.t <- quantile(v, prob=(1-tt))  
+
+  # find the levels on the two replicates
+  sig.value1 <- c()
+  sig.value2 <- c()
+  error.prob.cutoff <- c()
+  n.selected.match <- c()
+
+  for(i in 1:length(fdr.level)){
+
+    # find which uri.ez is closet to fdr.level
+    index <- which.min(abs(uri.ez - fdr.level[i]))
+    sig.value1[i] <- u.t[index]
+    sig.value2[i] <- v.t[index]
+    error.prob.cutoff[i] <- ez.bound[index]  
+    if(reverse){
+      n.selected.match[i] <- sum(e.z<=ez.bound[index])    
+    } else {
+      n.selected.match[i] <- sum(e.z>=ez.bound[index])    
+    }
+  }   
+
+  # output the cutoff of posterior probability, signal values on two replicates
+  map.uv <- cbind(error.prob.cutoff, sig.value1, sig.value2, n.selected.match)
+
+  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))
+}
+
+# created: 01-06-2010
+# Output IDR at various number of selected peaks
+# Find cutoff (idr cutoff, sig.value cutoff on each replicate) for specified IDR level
+# IDR definition is similar to FDR
+get.ez.tt <- function(em.fit, idr.level=c(0.01, 0.05, 0.1)){
+
+#  u <- em.fit$data.pruned$sample1$sig.value
+#  v <- em.fit$data.pruned$sample2$sig.value
+  u <- em.fit$data.pruned$sample1
+  v <- em.fit$data.pruned$sample2
+  
+  e.z <-  1-em.fit$em.fit$e.z # this is the error prob
+  
+  o <- order(e.z)
+  e.z.ordered <- e.z[o]
+  n.select <- c(1:length(e.z))
+  IDR <- cumsum(e.z.ordered)/n.select
+
+  u.o <- u[o]
+  v.o <- v[o]
+
+  n.level <- length(idr.level)
+#  sig.value1 <- rep(NA, n.level)
+#  sig.value2 <- rep(NA, n.level)
+  ez.cutoff <- rep(NA, n.level)
+  n.selected <- rep(NA, n.level)
+  
+  for(i in 1:length(idr.level)){
+
+    # find which uri.ez is closet to fdr.level
+    index <- which.min(abs(IDR - idr.level[i]))
+#    sig.value1[i] <- min(u.o[1:index])
+#    sig.value2[i] <- min(v.o[1:index])
+    ez.cutoff[i] <- e.z[index]      
+    n.selected[i] <- sum(e.z<=ez.cutoff[i])    
+  }   
+
+  # output the cutoff of posterior probability, number of selected overlapped peaks 
+#  map.uv <- cbind(ez.cutoff, sig.value1, sig.value2, n.selected)
+
+  map.uv <- cbind(ez.cutoff, n.selected)
+
+  return(list(n=n.select, IDR=IDR, idr.level=idr.level, map.uv=map.uv))
+}   
+  
+#  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))  
+  
+
+
+
+
+### compute the mean of the marginals
+get.mar.mean <- function(em.out){
+
+  x.f1 <- em.out$x.mar$f1
+  x.f2 <- em.out$x.mar$f2
+
+  y.f1 <- em.out$y.mar$f1
+  y.f2 <- em.out$y.mar$f2
+
+  x.stat1 <- get.hist.mean(x.f1)
+  x.stat2 <- get.hist.mean(x.f2)
+  y.stat1 <- get.hist.mean(y.f1)
+  y.stat2 <- get.hist.mean(y.f2)
+
+  return(list(x.mean1=x.stat1$mean, x.mean2=x.stat2$mean, 
+              y.mean1=y.stat1$mean, y.mean2=y.stat2$mean,
+              x.sd1=x.stat1$sd, x.sd2=x.stat2$sd, 
+              y.sd1=y.stat1$sd, y.sd2=y.stat2$sd
+              ))
+
+}
+
+
+# compute the mean of marginals
+get.hist.mean  <- function(x.f){
+
+  nbreaks <- length(x.f$breaks)
+  x.bin <- x.f$breaks[-1]-x.f$breaks[-nbreaks]
+
+  x.mid <- (x.f$breaks[-nbreaks]+x.f$breaks[-1])/2
+  x.mean <- sum(x.mid*x.f$density*x.bin)
+  x.sd <- sqrt(sum(x.mid*x.mid*x.f$density*x.bin)-x.mean^2)
+  
+  return(list(mean=x.mean, sd=x.sd))
+}
+
+get.hist.var <- function(x.f){
+
+  nbreaks <- length(x.f$breaks)
+  x.bin <- x.f$breaks[-1]-x.f$breaks[-nbreaks]
+
+  x.mid <- (x.f$breaks[-nbreaks]+x.f$breaks[-1])/2
+  x.mean <- sum(x.mid*x.f$density*x.bin)
+
+  return(mean=x.mean)  
+}
+
+# obsolete function (01-06-2010)
+# plot 
+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){
+
+  if(is.null(col.txt))
+    col.txt <- c("black", "red", "purple", "green", "blue", "cyan", "magenta", "orange", "grey")
+
+  x <- c()
+  y <- c()
+
+  for(i in 1:length(ez.list)){
+    x <- c(x, ez.list[[i]]$n)
+      
+    y <- c(y, ez.list[[i]]$uri.ez)
+  }
+
+  if(is.null(y.lim))
+    y.lim <- c(0, max(y))
+
+  if(!is.null(file.name)){
+    postscript(paste(plot.dir, "ez.", file.name, sep=""))
+    par(mfrow=c(1,1), mar=c(5,5,4,2))
+  }
+
+  plot(x, y, ylim=y.lim, type="n", xlab=xlab.txt, ylab=ylab.txt, lwd=5, cex=5, cex.axis=2, cex.lab=2)
+
+  for(i in 1:length(ez.list)){
+    lines(ez.list[[i]]$n, ez.list[[i]]$uri.ez, col=col.txt[i], cex=2, lwd=5)    
+  }
+
+#   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)
+#   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)
+  
+
+  legend(0, y.lim[2], legend=legend.txt, col=col.txt[1:length(col.txt)], lty=1, lwd=5, cex=2)
+
+  if(!is.null(title))
+    title(title.txt)
+
+  if(!is.null(file.name)){
+    dev.off()
+  }
+  
+}
+
+
+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){
+
+  if(is.null(col.txt))
+    col.txt <- c("black", "red", "purple", "green", "blue", "cyan", "magenta", "orange", "grey")
+
+  n.entry <- length(ez.list)
+  x <- rep(NA, n.entry)
+  y.max <- rep(NA, n.entry)
+
+  for(i in 1:n.entry){
+    x[i] <- max(ez.list[[i]]$n)
+      
+    y.max[i] <- max(ez.list[[i]]$IDR)
+  
+  }
+
+  if(is.null(y.lim))
+    y.lim <- c(0, max(y.max))
+
+  if(!is.null(file.name)){
+    postscript(paste(plot.dir, "ez.", file.name, sep=""))
+    par(mfrow=c(1,1), mar=c(5,5,4,2))
+  }
+
+
+  
+  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)
+
+  q <- seq(0.01, 0.99, by=0.01)
+  
+  for(i in 1:length(ez.list)){
+
+    n.plot <- round(quantile(ez.list[[i]]$n, prob=q))
+    IDR.plot <- ez.list[[i]]$IDR[n.plot]
+    lines(n.plot, IDR.plot, col=col.txt[i], cex=2, lwd=5)    
+  }
+
+
+  legend(0, y.lim[2], legend=legend.txt, col=col.txt[1:length(col.txt)], lty=1, lwd=5, cex=2)
+
+  if(!is.null(title))
+    title(title.txt)
+
+  if(!is.null(file.name)){
+    dev.off()
+  }
+  
+}
+
+
+
+#############################################################################
+#############################################################################
+# statistics about peaks selected on the individual replicates
+#
+# idr.level: the consistency cutoff, say 0.05
+# uri.output: a list of uri.output from consistency analysis generated by batch-consistency-analysis.r
+# ez.list : a list of IDRs computed from get.ez.tt using the same idr.level
+#
+##################
+
+
+# obsolete?
+# compute the error rate
+# u.t and v.t are the quantiles
+# 
+# map back to all peaks and report the number of peaks selected
+get.ez.tt.all.old  <- function(em.fit, all.data1, all.data2, idr.level){
+
+  u <- em.fit$data.pruned$sample1
+  v <- em.fit$data.pruned$sample2
+
+  tt <- seq(0.01, 0.99, by=0.01)
+#  if(reverse){ 
+    e.z <-  1-em.fit$em.fit$e.z # this is the error prob
+    uri.ez <- sapply(tt, comp.uri.ez, u=u, v=v, e.z=e.z)
+    ez.bound <- sapply(tt, comp.ez.cutoff, u=u, v=v, e.z=e.z, boundary.txt="max") 
+#  } else {
+#    e.z <-  em.fit$em.fit$e.z
+#    uri.ez <- sapply(tt, comp.uri.ez, u=u, v=v, e.z=e.z)
+#    ez.bound <- sapply(tt, comp.ez.cutoff, u=u, v=v, e.z=e.z, boundary.txt="min") 
+#  }
+
+  u.t <- quantile(u, prob=(1-tt))
+  v.t <- quantile(v, prob=(1-tt))  
+
+  # find the levels on the two replicates
+  sig.value1 <- c()
+  sig.value2 <- c()
+  error.prob.cutoff <- c()
+  n.selected.match <- c()
+  npeak.rep1 <- c()
+  npeak.rep2 <- c()
+
+  for(i in 1:length(idr.level)){
+
+    # find which uri.ez is closet to idr.level
+    index <- which.min(abs(uri.ez - as.numeric(idr.level[i])))
+
+    sig.value1[i] <- u.t[index]
+    sig.value2[i] <- v.t[index]
+    error.prob.cutoff[i] <- ez.bound[index]  
+    n.selected.match[i] <- sum(u>= u.t[index] & v>=v.t[index])
+
+    npeak.rep1[i] <- sum(all.data1["sig.value"] >= sig.value1[i])
+    npeak.rep2[i] <- sum(all.data2["sig.value"] >= sig.value2[i])    
+  }   
+
+
+  # output the cutoff of posterior probability, signal values on two replicates
+  map.uv <- cbind(error.prob.cutoff, sig.value1, sig.value2, n.selected.match, npeak.rep1, npeak.rep2)
+
+  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))
+}
+
+
+get.ez.tt.all <- function(em.fit, all.data1, all.data2, idr.level=c(0.01, 0.05, 0.1)){
+
+  u <- em.fit$data.pruned$sample1$sig.value
+  v <- em.fit$data.pruned$sample2$sig.value
+#  u <- em.fit$data.pruned$sample1
+#  v <- em.fit$data.pruned$sample2
+  
+  e.z <-  1-em.fit$em.fit$e.z # this is the error prob
+  
+  o <- order(e.z)
+  e.z.ordered <- e.z[o]
+  n.select <- c(1:length(e.z))
+  IDR <- cumsum(e.z.ordered)/n.select
+
+  u.o <- u[o]
+  v.o <- v[o]
+
+  n.level <- length(idr.level)
+#  sig.value1 <- rep(NA, n.level)
+#  sig.value2 <- rep(NA, n.level)
+  ez.cutoff <- rep(NA, n.level)
+  n.selected <- rep(NA, n.level)
+  npeak.rep1 <- rep(NA, n.level)
+  npeak.rep2 <- rep(NA, n.level)
+  
+  for(i in 1:length(idr.level)){
+
+    # find which uri.ez is closet to fdr.level
+    index <- which.min(abs(IDR - idr.level[i]))
+#    sig.value1[i] <- min(u.o[1:index])
+#    sig.value2[i] <- min(v.o[1:index])
+    ez.cutoff[i] <- e.z.ordered[index]      # fixed on 02/20/10
+    n.selected[i] <- sum(e.z<=ez.cutoff[i])
+#    npeak.rep1[i] <- sum(all.data1["sig.value"] >= sig.value1[i])
+#    npeak.rep2[i] <- sum(all.data2["sig.value"] >= sig.value2[i])     
+  }   
+
+  # output the cutoff of posterior probability, number of selected overlapped peaks 
+  map.uv <- cbind(ez.cutoff, n.selected)
+
+  return(list(n=n.select, IDR=IDR, idr.level=idr.level, map.uv=map.uv))
+}   
+  
+#  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))  
+  
+
+
+
+
+
+####### the following is for determining thresholds for merged dataset
+
+############# select peaks above a given threshold
+#
+# pass.threshold: a simple method, passing the threshold on the threshold on the individual replicate to the pooled sample 
+#                 
+# 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
+#          the 2nd column is the cutoff of ez, the rest of columns are consistency analysis for other replicates 
+# sig.value.name: the name of the sig.value column
+# combined: combined dataset
+# nrep: number of pairs of comparisons
+#
+# Procedure:
+# 1. Find the significant threshold corresponding to the idr cutoff on the matched peaks. 
+# 2. Each time we will get two or more (if >2 replicates) cutoffs and will report the most stringent and the least stringent
+#    cutoff and the number of peaks selected at those two cutoffs
+#############
+
+pass.threshold <- function(sig.map.list, sig.value.name, combined, idr.level, nrep, chr.size){
+
+  sig.map <- c()
+
+  # choose idr.level
+  idr.index <- which(rbind(sig.map.list[[1]])[,1] == idr.level)
+  if(length(i) ==0){
+    print("no level matches specified idr.level")
+    return(-1)
+  }
+
+  for(i in 1:length(sig.map.list))
+    sig.map <- c(sig.map, rbind(sig.map.list[[i]])[idr.index, c("sig.value1", "sig.value2")])
+  
+  
+  npeak.tight <- c()  
+  npeak.loose <- c()
+
+
+  max.sig <- max(sig.map)
+  min.sig <- min(sig.map)
+  selected.sig.tight <- combined[combined[,sig.value.name]>=max.sig, ]
+  selected.sig.loose <- combined[combined[,sig.value.name]>=min.sig, ]
+
+  selected.sig.tight <- deconcatenate.chr(selected.sig.tight, chr.size)[,c("chr", "start", "stop", "signal.value", "p.value", "q.value")]
+  selected.sig.loose <- deconcatenate.chr(selected.sig.loose, chr.size)[,c("chr", "start", "stop", "signal.value", "p.value", "q.value")]
+  
+  npeak.tight <- nrow(selected.sig.tight)
+  npeak.loose <- nrow(selected.sig.loose)
+  
+  
+  npeak.stat <- list(idr.level=idr.level, max.sig=max.sig, min.sig=min.sig, npeak.tight=npeak.tight, npeak.loose=npeak.loose)
+
+  invisible(list(npeak.stat=npeak.stat, combined.selected.tight=selected.sig.tight, combined.selected.loose=selected.sig.loose))  
+}
+
+#################
+# pass the regions selected from consistency analysis to combined data
+# Threshold is determined on the replicates, the regions above the threshold are selected
+# then peaks on the combined data are selected from the selected regions
+#
+# To avoid being too stringent, regions satisfying the following conditions are selected
+# 1. regions above the significant threshold determined by consistency analysis on either replicate
+# 2. regions that have consistent low peaks, i.e. posterior prob > threshold but not passing the significant threshold
+#
+# This method doesn't make a difference when using different thresholds
+#################
+
+pass.region <- function(sig.map.list, uri.output, ez.list, em.output, combined, idr.level, sig.value.impute=0, chr.size){
+ 
+  combined <- combined[, c("start", "stop", "sig.value", "signal.value", "p.value", "q.value")]
+  npair <- length(uri.output) # number of pairs of consistency analysis
+  combined.region <- c()
+
+  # choose idr.level
+  idr.index <- which(rbind(sig.map.list[[1]])[,1] == idr.level)
+  if(length(idr.index) ==0){
+    print("no level matches specified idr.level")
+    return(-1)
+  }
+
+    for(j in 1:npair){
+      # select peaks from individual replicates using individual cutoff
+      above.1 <- uri.output[[j]]$data12.enrich$merge1["sig.value"] >= ez.list[[j]]$map.uv[idr.index,"sig.value1"]
+      above.2 <- uri.output[[j]]$data12.enrich$merge1["sig.value"] >= ez.list[[j]]$map.uv[idr.index,"sig.value2"]
+      selected.sig.rep1 <- uri.output[[j]]$data12.enrich$merge1[above.1, c("start", "stop", "sig.value", "signal.value", "p.value", "q.value")]
+      selected.sig.rep2 <- uri.output[[j]]$data12.enrich$merge2[above.2, c("start", "stop", "sig.value", "signal.value", "p.value", "q.value")] 
+      
+      # find the peaks that are overlapped with reliable peaks in the individual replicates
+      overlap.1 <- pair.peaks(selected.sig.rep1, combined)$merge2
+      overlap.2 <- pair.peaks(selected.sig.rep2, combined)$merge2
+
+      # choose the ones with significant value > 0, which are the overlapped ones
+
+      combined.in1 <- overlap.1[overlap.1$sig.value > sig.value.impute, c("start", "stop", "sig.value", "signal.value", "p.value", "q.value")]
+      combined.in2 <- overlap.2[overlap.2$sig.value > sig.value.impute, c("start", "stop", "sig.value", "signal.value", "p.value", "q.value")]
+
+      ## consistent low significant ones
+      ## first find consistenct ones, ie. high posterior prob
+      # is.consistent <- ez.list[[j]]$e.z < ez.list[[j]]$ez.cutoff 
+
+      # data.matched <- keep.match(uri.output[[j]]$data12.enrich$merge1[!above.1, ], uri.output[[j]]$data12.enrich$merge2[!above.2, ], sig.value.impute=0)
+      # data.matched$sample1 <- data.matched$sample1[, c("start", "stop", "sig.value", "signal.value", "p.value", "q.value")]
+      # data.matched$sample2 <- data.matched$sample2[, c("start", "stop", "sig.value", "signal.value", "p.value", "q.value")]
+
+      # consistent.in1 <- data.matched$sample1[is.consistent, ]
+      # consistent.in2 <- data.matched$sample2[is.consistent, ]
+
+      # overlap.consistent.1 <- pair.peaks(consistent.in1, combined)$merge2
+      # overlap.consistent.2 <- pair.peaks(consistent.in2, combined)$merge2
+
+      ## choose the ones with significant value > 0, which are the overlapped ones
+
+      # combined.consistent.in1 <- overlap.consistent.1[overlap.consistent.1$sig.value > sig.value.impute, ]
+      # combined.consistent.in2 <- overlap.consistent.2[overlap.consistent.2$sig.value > sig.value.impute, ]
+
+      # combined.region <- rbind(combined.region, combined.in1, combined.in2, combined.consistent.in1, combined.consistent.in2)
+
+       combined.region <- rbind(combined.region, combined.in1, combined.in2)
+
+      is.repeated <- duplicated(combined.region$start)
+      combined.region <- combined.region[!is.repeated, c("start", "stop", "sig.value", "signal.value", "p.value", "q.value")]
+      
+    }
+    npeak <- nrow(combined.region)
+  
+   sig.combined <- c(min(combined.region[,"sig.value"], na.rm=T), max(combined.region[,"sig.value"], na.rm=T))
+
+  # idr.combined <- c(min(combined.region[,"q.value"], na.rm=T), max(combined.region[,"q.value"], na.rm=T))
+
+   npeak.stat <- list(idr.level=idr.level, npeak=npeak)
+
+   combined.region <- deconcatenate.chr(combined.region, chr.size)[,c("chr", "start", "stop", "signal.value", "p.value", "q.value")]
+
+  invisible(list(npeak.stat=npeak.stat, combined.selected=combined.region, sig.combined=sig.combined))
+}
+
+################
+# pass structure: this method does another round of inference on the combined data
+#
+# To make the mixture structure comparable on the replicates and the combined data, the 2nd inference is done on the peaks
+# at the reliable regions on the combined data, using rank transformed significant values. The mixture structure is estimated using my consistency analysis, which 
+# estimates marginal distributions of ranks using nonparametric ways. Then the significant values are found out.
+# There are several advantages to do it this way:  
+# 1. The premise of passing structure is that the means and variance (i.e. distribution) of two replicates should be the same
+#    The significant values on the two replicates clearly have different distributions. The structure estimated from consistency
+#    analysis will generate similar rank distribution on two replicates by its setup (i.e. same number of peaks are paired up).  
+# 2. Because pooled sample is a black box, the structure is more likely to be followed in the matched regions than other locations,
+#    after all, we don't know what other things are. If even the structure doesn't hold on the matched regions, 
+#    which is possible, let alone the other regions. Focusing on the reliable regions helps to get rid of those unknown noises.
+#  
+# 
+# modified on 2-20-10: reverse rank.combined, make big sig.value with small
+# ranks, to be consistent with f1 and f2
+################ 
+
+pass.structure <- function(uri.output, em.output, combined, idr.level, sig.value.impute, chr.size, overlap.ratio=0){
+
+  columns.keep <- c("sig.value", "start", "stop", "signal.value", "p.value", "q.value", "chr", "start.ori", "stop.ori")
+  combined <- combined[, columns.keep]
+  combined.selected.all <- c()
+
+  for(j in 1:npair){
+
+    sample1 <- uri.output[[j]]$data12.enrich$merge1[, columns.keep]
+    sample2 <- uri.output[[j]]$data12.enrich$merge2[, columns.keep]
+        
+    # find peaks on the matched region on the combined one
+    data.matched <- keep.match(sample1, sample2, sig.value.impute=sig.value.impute)
+
+    data.matched$sample1 <- data.matched$sample1[, columns.keep]
+    data.matched$sample2 <- data.matched$sample2[, columns.keep]
+
+    overlap.1 <- pair.peaks.filter(data.matched$sample1, combined, p.value.impute=sig.value.impute, overlap.ratio)$merge2
+    overlap.2 <- pair.peaks.filter(data.matched$sample2, combined, p.value.impute=sig.value.impute, overlap.ratio)$merge2
+
+    # choose the ones with significant value > sig.value.impute, which are the overlapped ones
+
+    combined.in1 <- overlap.1[overlap.1$sig.value > sig.value.impute, ]
+    combined.in2 <- overlap.2[overlap.2$sig.value > sig.value.impute, ]
+
+    combined.region <- rbind(combined.in1, combined.in2)
+  
+    is.repeated <- duplicated(combined.region$start)
+    combined.region <- combined.region[!is.repeated,]
+
+    # now rank the peaks in matched region
+    rank.combined <- rank(-combined.region$sig.value)
+    
+    # now transform the parameters estimated into the new scale
+    npeaks.overlap <- nrow(combined.region)
+    npeaks.consistent <- nrow(cbind(em.output[[j]]$data.pruned$sample1))
+
+    
+    # the breaks are the same for x and y
+    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)
+    # the first break boundary goes up when changing scale, need set it back to be a bit smaller than 1
+    f1$breaks[1] <- min(f1$breaks[1], 0.95)
+    
+    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)
+    # the first break boundary goes up when changing scale, need set it back to be a bit smaller than 1
+    f2$breaks[1] <- min(f2$breaks[1], 0.95)
+    
+    p <- em.output[[j]]$em.fit$para$p
+ 
+    # find the posterior probability
+    errorprob.combined <- get.comp2.prob(rank.combined, p, f1, f2)
+
+    # compute the FDR and find cutoff of posterior prob and the sig value
+    o <- order(errorprob.combined)
+    idr <- cumsum(errorprob.combined[o])/c(1:length(o))
+    idr.index <- which(idr > idr.level)[1]
+    errorprob.cutoff <- errorprob.combined[o][idr.index]
+
+    # find the minimum significant measure among selected peaks
+    sig.value <- min(combined.region$sig.value[o][1:idr.index])
+  #  sig.value <- quantile(combined.region$sig.value[o][1:idr.index], prob=0.05)
+#sig.value <- quantile(combined.region$sig.value[errorprob.combined<=errorprob.cutoff], prob=0.05)
+    
+    # apply the significant value on the whole pooled list
+    combined.selected <- combined[combined$sig.value >= sig.value,]
+
+    combined.selected.all <- rbind(combined.selected.all, combined.selected)
+  }
+
+  is.repeated <- duplicated(combined.selected.all$start)
+  combined.selected.all <- combined.selected.all[!is.repeated,]
+
+  npeak <- nrow(combined.selected.all)  
+  
+  npeak.stat <- list(idr.level=idr.level, npeak=npeak)
+
+  sig.combined <- c(min(combined.selected.all[,"sig.value"], na.rm=T), max(combined.selected.all[,"sig.value"], na.rm=T))
+
+ #  idr.combined <- c(min(combined.selected.all[,"q.value"], na.rm=T), max(combined.selected.all[,"q.value"], na.rm=T))
+ # combined.selected.all <- deconcatenate.chr(combined.selected.all, chr.size)[,c("chr", "start", "stop", "signal.value", "p.value", "q.value")]
+
+  combined.selected.all <- combined.selected.all[,  c("chr", "start.ori", "stop.ori", "signal.value", "p.value", "q.value")]
+  colnames(combined.selected.all) <- c("chr", "start", "stop", "signal.value", "p.value", "q.value")
+  
+  invisible(list(npeak.stat=npeak.stat, combined.selected=combined.selected.all, sig.combined=sig.combined))
+}
+
+
+
+# get the posterior probability of the 2nd component
+get.comp2.prob <- function(x, p, f1, f2){
+
+  # get pdf and cdf of each component from functions in the corresponding component 
+  px.1 <- sapply(x, get.pdf, df=f1)
+  px.2 <- sapply(x, get.pdf, df=f2)
+
+  comp2prob <- 1 - p*px.1/(p*px.1+(1-p)*px.2)
+  
+  return(comp2prob)
+}
+
+keep.match <- function(sample1, sample2, sig.value.impute=0){
+
+  sample1.prune <- sample1[sample1$sig.value > sig.value.impute & sample2$sig.value > sig.value.impute,]
+  sample2.prune <- sample2[sample1$sig.value > sig.value.impute & sample2$sig.value > sig.value.impute,]
+ 
+  invisible(list(sample1=sample1.prune, sample2=sample2.prune))
+}
+
+
+##############################################
+#
+# The following is for simulation
+#
+##############################################
+
+
+# simulate gaussian copula
+# u is the uniform random variable and rho is correlation coefficient 
+simu.gaussian.copula <- function(u, rho){
+
+  n <- length(u)
+
+  # simulate y given x=qnorm(u)
+  y <- qnorm(u)*rho + rnorm(n)*sqrt(1-rho^2)
+
+  v <- pnorm(y)
+
+  invisible(v)
+}
+
+## simulate Clayton copula from its generating function
+## Genest and MacKay (1986)
+
+phi.ori <- function(t, s){
+
+ (t^(-s) -1)/s
+}
+
+
+phi.inv <- function(y, s){
+
+  exp(-log(s*y+1)/s)
+}
+
+phi.der <- function(t, s){
+
+  -t^(-s-1)
+}
+
+phi.der.inv <- function(y, s){
+
+  exp(log(-y)/(-s-1))
+}
+
+get.w <- function(u, t, s){
+
+  phi.der.inv(phi.der(u, s)/t, s)
+} 
+
+get.v <- function(w, u, s){
+
+  phi.inv(phi.ori(w, s) - phi.ori(u, s), s) 
+}
+
+# u is a uniform random variable, s is the association parameter
+simu.clayton.copula <- function(u, s){
+
+  t <- runif(length(u))
+
+  if(s>0){
+    w <- get.w(u, t, s)
+    v <- get.v(w, u, s)
+    return(v)    
+  }
+
+  if(s==0){
+    return(t)
+  }
+
+  if(s <0){
+    print("Invalid association parameters for clayton copula")
+  }
+  
+}
+
+
+
+###### 09-09-09
+
+# simulate a two-component copula mixture:
+# - marginal distributions for the two variables in each component are both 
+#   normal and with the same parameters 
+# p is the mixing proportion of component 1
+# n is the total sample size
+simu.copula.2mix <- function(s1, s2, p, n, mu1, mu2, sd1, sd2, copula.txt){
+
+  n1 <- round(n*p)
+  n2 <- n-n1
+
+  u1 <- runif(n1)
+  
+  if(copula.txt =="clayton")
+    v1 <- simu.clayton.copula(u1, s1)
+  else{
+    if(copula.txt =="gaussian")
+      v1 <- simu.gaussian.copula(u1, s1)
+  }
+
+  u2 <- runif(n2)
+
+  if(copula.txt =="clayton")
+    v2 <- simu.clayton.copula(u2, s2)
+  else{
+    if(copula.txt =="gaussian")
+      v2 <- simu.gaussian.copula(u2, s2)
+  }
+
+  # generate test statistics
+  sample1.1 <- qnorm(u1, mu1, sd1)
+  sample1.2 <- qnorm(v1, mu1, sd1)
+
+  sample2.1 <- qnorm(u2, mu2, sd2)
+  sample2.2 <- qnorm(v2, mu2, sd2)
+
+  return(list(u=c(u1, u2), v=c(v1, v2), 
+              u.inv=c(sample1.1, sample2.1), v.inv=c(sample1.2, sample2.2),
+              label=c(rep(1, n1), rep(2, n2))))
+}
+
+# using inverse of the cdf to generate original observations 
+
+simu.copula.2mix.inv <- function(s1, s2, p, n, cdf1.x, cdf1.y, cdf2.x, cdf2.y, copula.txt){
+
+  n1 <- round(n*p)
+  n2 <- n-n1
+
+  u1 <- runif(n1)
+  
+  if(copula.txt =="clayton")
+    v1 <- simu.clayton.copula(u1, s1)
+  else{
+    if(copula.txt =="gaussian")
+      v1 <- simu.gaussian.copula(u1, s1)
+  }
+
+  u2 <- runif(n2)
+
+  if(copula.txt =="clayton")
+    v2 <- simu.clayton.copula(u2, s2)
+  else{
+    if(copula.txt =="gaussian")
+      v2 <- simu.gaussian.copula(u2, s2)
+  }
+
+  # generate test statistics
+#  sample1.1 <- qnorm(u1, mu1, sd1)
+#  sample1.2 <- qnorm(v1, mu1, sd1)
+
+#  sample2.1 <- qnorm(u2, mu2, sd2)
+#  sample2.2 <- qnorm(v2, mu2, sd2)
+  
+  sample1.x <- inv.cdf.vec(u1, cdf1.x)
+  sample1.y <- inv.cdf.vec(v1, cdf1.y)
+
+  sample2.x <- inv.cdf.vec(u2, cdf2.x)
+  sample2.y <- inv.cdf.vec(v2, cdf2.y)
+  
+  
+  return(list(u=c(u1, u2), v=c(v1, v2), 
+              u.inv=c(sample1.x, sample2.x), v.inv=c(sample1.y, sample2.y),
+              label=c(rep(1, n1), rep(2, n2))))
+}
+
+# obtain original observation by converting cdf into quantiles
+# u is one cdf
+# u.cdf is a cdf (assuming it is a histogram) and has the break points (cdf$cdf and cdf$breaks)
+# the smallest value of cdf=0 and the largest =1 
+inv.cdf <- function(u, u.cdf){
+
+  # which bin it falls into
+  i <- which(u.cdf$cdf> u)[1]
+  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]
+
+  return(q.u)
+}
+
+inv.cdf.vec <- function(u, u.cdf){
+
+  # check if cdf has the right range (0, 1)  
+  ncdf <- length(u.cdf$cdf)
+  nbreaks <- length(u.cdf$breaks)
+  
+  if(ncdf == nbreaks-1 & u.cdf$cdf[ncdf]< 1)
+    u.cdf[ncdf] <- 1
+    
+  q.u <- sapply(u, inv.cdf, u.cdf)
+
+  return(q.u) 
+}
+
+# here we simulate a likely real situation
+# the test statistics from two normal distributions
+# according to their labels, then convert them into p-values w.r.t H0 using
+# one-sided test.
+# The test statistics are correlated for the signal component and independent
+# for the noise component
+# For the signal component, Y = X + eps, where eps ~ N(0, sigma^2)
+simu.test.stat <- function(p, n, mu1, sd1, mu0, sd0, sd.e){
+
+  # first component - signal
+  n.signal <- round(n*p)
+  n.noise <- n - n.signal
+
+  # labels
+  labels <- c(rep(1, n.signal), rep(0, n.noise))
+  
+  # test statistics for signal and noise
+  mu.signal <- rnorm(n.signal, mu1, sd1)
+  x.signal <- mu.signal + rnorm(n.signal, 0, sd.e)
+  x.noise <- rnorm(n.noise, mu0, sd0) + rnorm(n.noise, 0, sd.e)
+
+  y.signal <- mu.signal + rnorm(n.signal, 0, sd.e)
+                                        # sd.e can be dependent on signal
+  y.noise <- rnorm(n.noise, mu0, sd0) + rnorm(n.noise, 0, sd.e)
+
+  # concatenate
+  x <- c(x.signal, x.noise)
+  y <- c(y.signal, y.noise)
+  
+  # convert to p-values based on H0
+  p.x <- 1-pnorm(x, mu0, sqrt(sd0^2+sd.e^2))
+  p.y <- 1-pnorm(y, mu0, sqrt(sd0^2+sd.e^2))
+
+  return(list(p.x=p.x, p.y=p.y, x=x, y=y, labels=labels))
+  
+}
+
+# compute the tradeoff and calibration
+forward.decoy.tradeoff.ndecoy <- function(xx, labels, ndecoy){
+
+  xx <- round(xx, 5) 
+  o <- order(xx, decreasing=T)
+
+  rand <- 1-labels # if rand==0, consistent
+  # order the random indicator in the same order
+  rand.o <- rand[o]
+
+  if(sum(rand.o) > ndecoy){
+    index.decoy <- which(cumsum(rand.o)==ndecoy)
+  } else {
+    index.decoy <- which(cumsum(rand.o)==sum(rand.o))
+  }
+  
+  cutoff.decoy <- xx[o][index.decoy]
+    
+  # only consider the unique ones
+  cutoff.unique <- unique(xx[o])
+
+  cutoff <- cutoff.unique[cutoff.unique >= cutoff.decoy[length(cutoff.decoy)]]
+
+   get.decoy.count <- function(cut.off){
+     above <- rep(0, length(xx))
+     above[xx >= cut.off] <- 1
+     decoy.count <- sum(above==1 & rand==1)     
+     return(decoy.count)
+   }
+
+   get.forward.count <- function(cut.off){
+     above <- rep(0, length(xx))
+     above[xx >= cut.off] <- 1
+     forward.count <- sum(above==1 & rand==0)
+     return(forward.count)
+   }
+
+   get.est.fdr <- function(cut.off){
+     above <- rep(0, length(xx))
+     above[xx >= cut.off] <- 1
+     est.fdr <- 1-mean(xx[above==1])
+     return(est.fdr)
+   }
+
+  # assuming rand=0 is right
+   get.false.neg.count <- function(cut.off){
+     below <- rep(0, length(xx))
+     below[xx < cut.off] <- 1
+     false.neg.count <- sum(below==1 & rand==0)
+     return(false.neg.count)
+   }
+
+  get.false.pos.count <- function(cut.off){
+     above <- rep(0, length(xx))
+     above[xx >= cut.off] <- 1
+     false.pos.count <- sum(above==1 & rand==1)
+     return(false.pos.count)
+   } 
+
+   decoy <- sapply(cutoff, get.decoy.count)
+   forward <- sapply(cutoff, get.forward.count)
+
+   est.fdr <- sapply(cutoff, get.est.fdr)
+   emp.fdr <- decoy/(decoy+forward)
+  
+   # compute specificity and sensitivity
+   # assuming rand=1 is wrong and rand=0 is right
+   false.neg <- sapply(cutoff, get.false.neg.count)
+   false.pos <- sapply(cutoff, get.false.pos.count)
+  
+   true.pos <- sum(rand==0)-false.neg
+   true.neg <- sum(rand==1)-false.pos
+  
+   sensitivity <- true.pos/(true.pos+false.neg)
+   specificity <- true.neg/(true.neg+false.pos)
+  
+   return(list(decoy=decoy, forward=forward, cutoff=cutoff, est.fdr=est.fdr, emp.fdr=emp.fdr, sensitivity=sensitivity, specificity=specificity))  
+}
+
+
+# compute the em for jackknife and all data, and find FDR
+get.emp.jack <- function(a, p0){
+
+  nobs <- length(a$labels)
+  est <- list()
+  est.all <- list()
+
+  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)
+#  temp.all <- em.2copula.quick(a$p.x, a$p.y, p0=p0, rho1.0=0.7,
+#      rho2.0=0, eps=0.01, fix.p=T, stoc=F, fix.rho2=T, "gaussian")
+
+  est.all$p <- temp.all$para$p
+  est.all$rho1 <- temp.all$para$rho1
+  est.all$FDR <- get.FDR(temp.all$e.z)
+
+  FDR <- list()
+  p <- c()
+  rho1 <- c()
+
+
+  for(i in 1:nobs){
+
+    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)    
+#    temp <- em.2copula.quick(a$p.x[-i], a$p.y[-i], p0=p0, rho1.0=0.7,
+#      rho2.0=0, eps=0.01, fix.p=T, stoc=F, fix.rho2=T, "gaussian")
+
+    est[[i]] <- list(p=temp$para$p, rho1=temp$para$rho1, FDR=get.FDR(temp$e.z))
+
+    FDR[[i]] <- est[[i]]$FDR # this is the FDR for top n peaks
+    p[i] <- est[[i]]$p
+    rho1[i] <- est[[i]]$rho1 
+  }
+
+  est.jack <- list(FDR=FDR, p=p, rho1=rho1) 
+  return(list(est.jack=est.jack, est.all=est.all))
+}
+
+
+# get the npeaks corresponding to the nominal FDR estimated from the sample
+# and find the corresponding FDR from the entire data
+get.FDR.jack <- function(est, FDR.nominal){
+  
+  nobs <- length(est$est.jack$FDR)
+  FDR.all <- c()  
+  top.n <- c()
+  
+  for(i in 1:nobs){
+    top.n[i] <- max(which(est$est.jack$FDR[[i]] <= FDR.nominal))
+    FDR.all[i] <- est$est.all$FDR[top.n[i]]
+  }
+  
+  invisible(list(FDR.all=FDR.all, top.n=top.n))
+}
+
+# compute Jackknife peudonumber
+# a is the dataset
+get.emp.IF <- function(a, p0){
+
+  nobs <- length(a$labels)
+  est <- list()
+  est.all <- list()
+
+  temp.all <- em.2copula.quick(a$p.x, a$p.y, p0=p0, rho1.0=0.7,
+      rho2.0=0, eps=0.01, fix.p=T, stoc=F, fix.rho2=T, "gaussian")
+
+  est.all$p <- temp.all$para$p
+  est.all$rho1 <- temp.all$para$rho1
+  est.all$FDR <- get.FDR(temp.all$e.z)
+
+  IF.FDR <- list()
+  IF.p <- c()
+  IF.rho1 <- c()
+  
+  for(i in 1:nobs){
+    
+    temp <- em.2copula.quick(a$p.x[-i], a$p.y[-i], p0=p0, rho1.0=0.7,
+      rho2.0=0, eps=0.01, fix.p=T, stoc=F, fix.rho2=T, "gaussian")
+
+    est[[i]] <- list(p=temp$para$p, rho1=temp$para$rho1, FDR=get.FDR(temp$e.z))
+
+    IF.FDR[[i]] <- (nobs-1)*(est.all$FDR[-nobs] - est[[i]]$FDR) # this is the FDR for top n peaks
+    IF.p[i] <- (nobs-1)*(est.all$p - est[[i]]$p)
+    IF.rho1[i] <- (nobs-1)*(est.all$rho1 - est[[i]]$rho1) 
+  }
+
+  emp.IF <- list(FDR=IF.FDR, p=IF.p, rho1=IF.rho1) 
+
+  invisible(list(emp.IF=emp.IF, est.all=est.all, est=est))
+}
+
+# e.z is the posterior probability of being in signal component
+get.FDR <- function(e.z){
+
+  e.z.o <- order(1-e.z)
+  FDR <- cumsum(1-e.z[e.z.o])/c(1:length(e.z.o))
+
+  invisible(FDR)
+}
+
+# get the FDR of selecting the top n peaks
+# IF.est is the sample influence function
+# top.n
+get.IF.FDR <- function(IF.est, top.n){
+
+  nobs <- length(IF.est$emp.IF$FDR)
+  FDR <- c()  
+  
+  # influence function of p
+  for(i in 1:nobs)
+    FDR[i] <- IF.est$emp.IF$FDR[[i]][top.n]
+  
+  invisible(FDR)
+}
+
+# get the sample influence function for FDR at a given FDR size
+# 1. find the number of peaks selected at a given FDR computed from all obs
+# 2. use the number to find the sample influence function for FDR
+# IF.est$est.all is the FDR with all peaks
+get.IF.FDR.all <- function(IF.est, FDR.size){
+
+  top.n <- which.min(abs(IF.est$est.all$FDR -FDR.size))
+  nobs <- length(IF.est$est.all$FDR)
+  FDR <- c()  
+  
+  # influence function of p
+  for(i in 1:nobs)
+    FDR[i] <- IF.est$emp.IF$FDR[[i]][top.n]
+  
+  invisible(list(FDR=FDR, top.n=top.n))
+}
+
+plot.simu.uri <- function(x, y){
+
+  tt <- seq(0.01, 0.99, by=0.01)
+  uri <- sapply(tt, comp.uri.prob, u=x, v=y)
+  uri.thin <- uri[seq(1, length(tt), by=3)]
+  tt.thin <- tt[seq(1, length(tt), by=3)]
+  duri <- (uri.thin[-1]-uri.thin[-length(uri.thin)])/(tt.thin[-1]-tt.thin[-length(tt.thin)])
+  uri.spl <- smooth.spline(tt, uri, df=6.4)
+  uri.der <- predict(uri.spl, tt, deriv=1)
+
+  par(mfrow=c(2,2))
+  plot(x[1:n0], y[1:n0])
+  points(x[(n0+1):n], y[(n0+1):n], col=2)
+  plot(rank(-x)[1:n0], rank(-y)[1:n0])
+  points(rank(-x)[(1+n0):n], rank(-y)[(1+n0):n])
+  plot(tt, uri)
+  lines(c(0,1), c(0,1), lty=2)
+  title(paste("rho1=", rho1, " rho2=", rho2, "p=", p, sep=""))
+  plot(tt.thin[-1], duri)
+  lines(uri.der)
+  abline(h=1)
+  invisible(list(x=x, y=y, uri=uri, tt=tt, duri=duri, tt.thin=tt.thin, uri.der=uri.der))
+
+}
+
+
+###### new fitting procedure
+
+
+
+
+# 1. rank pairs
+
+# 2. initialization
+# 3. convert to pseudo-number
+
+# 4. EM
+
+# need plugin and test
+# find the middle point between the bins
+get.pseudo.mix <- function(x, mu, sigma, rho, p){
+
+  
+  # first compute cdf for points on the grid
+  # generate 200 points between [-3, mu+3*sigma]
+  nw <- 1000
+  w <- seq(min(-3, mu-3*sigma), max(mu+3*sigma, 3), length=nw) 
+  w.cdf <- p*pnorm(w, mean=mu, sd=sigma) + (1-p)*pnorm(w, mean=0, sd=1)
+
+  i <- 1
+
+  quan.x <- rep(NA, length(x))
+
+  for(i in c(1:nw)){
+    index <- which(x >= w.cdf[i] & x < w.cdf[i+1])
+    quan.x[index] <- (x[index]-w.cdf[i])*(w[i+1]-w[i])/(w.cdf[i+1]-w.cdf[i]) +w[i]
+  }
+
+  index <- which(x < w.cdf[1])
+  if(length(index)>0)
+    quan.x[index] <- w[1]
+
+  index <- which(x > w.cdf[nw])
+  if(length(index)>0)
+    quan.x[index] <- w[nw]  
+  
+#  linear.ext <- function(x, w, w.cdf){
+  # linear interpolation
+#    index.up <- which(w.cdf>= x)[1]
+#    left.index <- which(w.cdf <=x)
+#    index.down <- left.index[length(left.index)]
+#    quan.x <- (w[index.up] + w[index.down])/2  
+#  }
+  
+#  x.pseudo <- sapply(x, linear.ext, w=w, w.cdf=w.cdf)
+
+#  invisible(x.pseudo)
+  invisible(quan.x)
+}
+
+
+# EM to compute the latent structure
+# steps:
+# 1. raw values are first transformed into pseudovalues
+# 2. EM is used to compute the underlining structure, which is a mixture
+#    of two normals
+em.transform <- function(x, y, mu, sigma, rho, p, eps){
+  
+  x.cdf.func <- ecdf(x)
+  y.cdf.func <- ecdf(y)
+  afactor <- length(x)/(length(x)+1)
+  x.cdf <- x.cdf.func(x)*afactor
+  y.cdf <- y.cdf.func(y)*afactor
+  
+  # initialization
+  para <- list()
+  para$mu <- mu
+  para$sigma <- sigma
+  para$rho <- rho
+  para$p <- p  
+
+  j <- 1
+  to.run <- T
+  loglik.trace <- c()
+  loglik.inner.trace <- c()
+  
+  #to.run.inner <- T
+  z.1 <- get.pseudo.mix(x.cdf, para$mu, para$sigma, para$rho, para$p)
+  z.2 <- get.pseudo.mix(y.cdf, para$mu, para$sigma, para$rho, para$p)
+
+#  cat("length(z1)", length(z.1), "\n")
+  while(to.run){
+    
+    # get pseudo value in each cycle
+#    z.1 <- get.pseudo.mix(x.cdf, para$mu, para$sigma, para$rho, para$p)
+#    z.2 <- get.pseudo.mix(y.cdf, para$mu, para$sigma, para$rho, para$p)
+
+    i <- 1
+    while(to.run){
+      
+      # EM for latent structure
+      e.z <- e.step.2normal(z.1, z.2, para$mu, para$sigma, para$rho, para$p)
+      para <- m.step.2normal(z.1, z.2, e.z)
+#para$rho <- rho
+#para$p <- p    
+#para$mu <- mu
+#para$sigma <- sigma    
+      if(i > 1)
+        l.old <- l.new
+    
+      # this is just the mixture likelihood of two-component Gaussian
+      l.new <- loglik.2binormal(z.1, z.2, para$mu, para$sigma, para$rho, para$p)
+
+      loglik.inner.trace[i] <- l.new 
+
+      if(i > 1){
+        to.run <- loglik.inner.trace[i]-loglik.inner.trace[i-1]>eps         
+      }
+        
+    
+#      if(i > 2){
+#        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]))
+
+#        if(loglik.inner.trace[i-1]!=loglik.inner.trace[i-2])
+#          to.run <- abs(l.inf - loglik.inner.trace[i]) > eps
+#        else
+#          to.run <- F
+          
+#      }
+
+      cat("loglik.inner.trace[", i, "]=", loglik.inner.trace[i], "\n")
+    cat("mu=", para$mu, "sigma=", para$sigma, "p=", para$p, "rho=", para$rho, "\n\n")
+      
+      i <- i+1
+    }
+    
+
+    # get pseudo value in each cycle
+    z.1 <- get.pseudo.mix(x.cdf, para$mu, para$sigma, para$rho, para$p)
+    z.2 <- get.pseudo.mix(y.cdf, para$mu, para$sigma, para$rho, para$p)
+
+    if(j > 1)
+      l.old.outer <- l.new.outer
+
+    l.new.outer <- loglik.2binormal(z.1, z.2, para$mu, para$sigma, para$rho, para$p)
+
+    loglik.trace[j] <- l.new.outer
+    
+    if(j == 1)
+      to.run <- T
+    else{ # stop when iteration>100
+      if(j > 100)
+        to.run <- F
+      else
+        to.run <- l.new.outer - l.old.outer > eps
+    }
+
+#    if(j %% 10==0)
+      cat("loglik.trace[", j, "]=", loglik.trace[j], "\n")
+    cat("mu=", para$mu, "sigma=", para$sigma, "p=", para$p, "rho=", para$rho, "\n")
+    
+    j <- j+1
+  }
+
+  bic <- -2*l.new + 4*log(length(z.1))
+  
+  return(list(para=list(p=para$p, rho=para$rho, mu=para$mu, sigma=para$sigma),
+              loglik=l.new, bic=bic, e.z=e.z, loglik.trace=loglik.trace))
+}  
+
+
+
+
+# compute log-likelihood for mixture of two bivariate normals
+loglik.2binormal <- function(z.1, z.2, mu, sigma, rho, p){
+
+  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)))
+  
+#  l.m <- sum((p*d.binormal(z.1, z.2, mu, sigma, rho) + (1-p)*d.binormal(z.1, z.2, 0, 1, 0)))
+  return(l.m) 
+}
+
+# check this when rho=1
+
+# density of binomial distribution with equal mean and sigma on both dimensions
+d.binormal <- function(z.1, z.2, mu, sigma, rho){
+
+  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))
+
+  return(loglik)
+}
+
+# E-step for computing the latent strucutre
+# e.z is the prob to be in the consistent group
+# e.step for estimating posterior prob
+# z.1 and z.2 can be vectors or scalars
+e.step.2normal <- function(z.1, z.2, mu, sigma, rho, p){
+
+  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)
+  
+  invisible(e.z)
+}
+
+# M-step for computing the latent structure
+# m.step for estimating proportion, mean, sd and correlation coefficient
+m.step.2normal <- function(z.1, z.2, e.z){
+
+  p <- mean(e.z)
+  mu <- sum((z.1+z.2)*e.z)/2/sum(e.z) 
+  sigma <- sqrt(sum(e.z*((z.1-mu)^2+(z.2-mu)^2))/2/sum(e.z))
+  rho <- 2*sum(e.z*(z.1-mu)*(z.2-mu))/(sum(e.z*((z.1-mu)^2+(z.2-mu)^2)))
+
+  return(list(p=p, mu=mu, sigma=sigma, rho=rho))
+}
+
+
+# assume top p percent of observations are true
+# x and y are ranks, estimate
+init <- function(x, y, x.label){
+  
+  x.o <- order(x)
+
+  x.ordered <- x[x.o]
+  y.ordered <- y[x.o]
+  x.label.ordered <- x.label[x.o]
+  
+  n <- length(x)
+  p <- sum(x.label)/n
+  
+  rho <- cor(x.ordered[1:ceiling(p*n)], y.ordered[1:ceiling(p*n)])
+
+  temp <- find.mu.sigma(x.ordered, x.label.ordered)
+  mu <- temp$mu
+  sigma <- temp$sigma
+  
+  invisible(list(mu=mu, sigma=sigma, rho=rho, p=p))
+
+}
+
+# find mu and sigma if the distributions of marginal ranks are known
+# take the medians of the two dist and map back to the original
+init.dist <- function(f0, f1){
+
+  # take the median in f0
+  index.median.0 <- which(f0$cdf>0.5)[1]
+  q.0.small <- f0$cdf[index.median.0] # because f0 and f1 have the same bins
+  q.1.small <- f1$cdf[index.median.0]
+
+  # take the median in f1
+  index.median.1 <- which(f1$cdf>0.5)[1]
+  q.0.big <- f0$cdf[index.median.1] # because f0 and f1 have the same bins
+  q.1.big <- f1$cdf[index.median.1]
+
+  # find pseudo value for x.middle[1] on normal(0,1) 
+  pseudo.small.0 <- qnorm(q.0.small, mean=0, sd=1)
+  pseudo.small.1 <- qnorm(q.1.small, mean=0, sd=1)
+
+  # find pseudo value for x.middle[2] on normal(0,1) 
+  pseudo.big.0 <- qnorm(q.0.big, mean=0, sd=1)
+  pseudo.big.1 <- qnorm(q.1.big, mean=0, sd=1)
+
+  mu <- (pseudo.small.0*pseudo.big.1 - pseudo.small.1*pseudo.big.0)/(pseudo.big.1-pseudo.small.1) 
+
+  sigma <- (pseudo.small.0-mu)/pseudo.small.1
+
+  return(list(mu=mu, sigma=sigma))  
+}
+
+# generate labels
+
+# find the part of data with overlap
+
+# find the percentile on noise and signal
+
+# Suppose there are signal and noise components, with mean=0 and sd=1 for noise
+# x and x.label are the rank of the observations and their labels,
+# find the mean and sd of the other component
+# x.label takes values of 0 and 1
+find.mu.sigma <- function(x, x.label){
+
+  x.0 <- x[x.label==0]
+  x.1 <- x[x.label==1]
+
+  n.x0 <- length(x.0)
+  n.x1 <- length(x.1)
+
+  x.end <- c(min(x.0), min(x.1), max(x.0), max(x.1))
+  o <- order(x.end)
+  x.middle <- x.end[o][c(2,3)]
+
+  # the smaller end of the overlap
+  q.1.small <- mean(x.1 <= x.middle[1])*n.x1/(n.x1+1)
+  q.0.small <- mean(x.0 <= x.middle[1])*n.x0/(n.x0+1)
+
+  # the bigger end of the overlap
+  q.1.big <- mean(x.1 <= x.middle[2])*n.x1/(n.x1+1)
+  q.0.big <- mean(x.0 <= x.middle[2])*n.x0/(n.x0+1)
+
+  # find pseudo value for x.middle[1] on normal(0,1) 
+  pseudo.small.0 <- qnorm(q.0.small, mean=0, sd=1)
+  pseudo.small.1 <- qnorm(q.1.small, mean=0, sd=1)
+
+  # find pseudo value for x.middle[2] on normal(0,1) 
+  pseudo.big.0 <- qnorm(q.0.big, mean=0, sd=1)
+  pseudo.big.1 <- qnorm(q.1.big, mean=0, sd=1)
+
+  mu <- (pseudo.small.0*pseudo.big.1 - pseudo.small.1*pseudo.big.0)/(pseudo.big.1-pseudo.small.1) 
+
+  sigma <- (pseudo.small.0-mu)/pseudo.small.1
+
+  return(list(mu=mu, sigma=sigma))
+}