R/LSD_functions.R
af47efdf
 ._get_first_nonzero_bin <- function(Brick = NULL, chr = NULL, resolution = NA){
632a1a64
     RowSums <- Brick_get_matrix_mcols(Brick = Brick, chr1 = chr, chr2 = chr, 
974f4033
         resolution = resolution, what = "chr1_row_sums")
a47def9c
     return(min(which(RowSums > 0)))
86474139
 }
af47efdf
 ._get_sparsity_index <- function(Brick = NULL, chr = NULL, resolution = NA){
632a1a64
     Sparsity.index <- Brick_get_matrix_mcols(Brick = Brick, chr1 = chr, 
af47efdf
         chr2 = chr, resolution = resolution, what = "sparsity")
625325f1
     return(Sparsity.index)
86474139
 }
625325f1
 Backwards.Difference <- function(Vector=NULL,sparse=FALSE,sparsity.idx=NULL,
6a22d578
     sparsity_threshold=NULL){
86474139
         if(is.null(Vector)){
             stop("Vector variable cannot be empty.")
         }
         if(sparse){
             VectorDiff <- rep(NA,length(Vector))
625325f1
             temp.diff <- rev(diff(
6a22d578
                 rev(Vector[sparsity.idx > sparsity_threshold])))
86474139
             temp.diff[length(temp.diff)+1] <- 0
6a22d578
             VectorDiff[sparsity.idx > sparsity_threshold] <- temp.diff
86474139
         }else{
             VectorDiff<-rev(diff(rev(Vector)))
             VectorDiff[length(Vector)]=0
         }
         return(VectorDiff)
 }
625325f1
 Forwards.Difference <- function(Vector=NULL,sparse=FALSE,sparsity.idx=NULL,
6a22d578
     sparsity_threshold=NULL){
86474139
         if(is.null(Vector)){
             stop("Vector variable cannot be empty.")
         }
         if(sparse){
             VectorDiff <- rep(NA,length(Vector))
6a22d578
             temp.diff <- diff(Vector[sparsity.idx > sparsity_threshold])
86474139
             temp.diff <- c(0,temp.diff)
6a22d578
             VectorDiff[sparsity.idx > sparsity_threshold] <- temp.diff
86474139
         }else{
         VectorDiff<-diff(Vector)
         VectorDiff=c(0,VectorDiff)
         }
         return(VectorDiff)
 }
6a22d578
 ComputeOutlierOverIQRInsideWindow <- function(lookup_window=NULL,
625325f1
     diff.values=NULL,values=NULL, row.sums = NULL,
6a22d578
         min_sum = NULL, tukeys_constant=NULL,tail=NULL,
         sparse=FALSE,strict = FALSE,sparsity.idx=NULL,sparsity_threshold=NULL){
86474139
         seq.over.value <- seq_along(diff.values)
6a22d578
         Filter <- row.sums > min_sum
86474139
         if(sparse){
6a22d578
             Filter <- Filter & sparsity.idx > sparsity_threshold
86474139
         }
a47def9c
         seq.over.value <- seq.over.value[Filter]
86474139
         seq.over.seq <- seq_along(seq.over.value)
         outlier.list <- lapply(seq.over.seq,function(x.seq){
6a22d578
             lookup_window.range <- (
                 (x.seq - lookup_window) : (x.seq + lookup_window))
             lookup_window.range <- seq.over.value[lookup_window.range[
             lookup_window.range>0 & lookup_window.range<=max(seq.over.seq)]]
             offset <- (min(lookup_window.range)-1)
             diff.value.window <- diff.values[lookup_window.range]
             value.window <- values[lookup_window.range]
86474139
             value.quartile<-quantile(diff.value.window,na.rm=TRUE)
             InterQuartile <- value.quartile[4]-value.quartile[2]
             if(tail=="lower.tail"){
                 #Calculate Inner fences based on accepted formula
6a22d578
                 fences <- value.quartile[2] - (InterQuartile*tukeys_constant)
86474139
                 Outlier.Filter <- !is.na(value.window) &
                 diff.value.window <= fences &
                 diff.value.window < value.window
                 if(strict){
                     Outlier.Filter <- Outlier.Filter & (value.window < 0)
                 }
             }else if(tail=="upper.tail"){
6a22d578
                 fences <- value.quartile[4] + (InterQuartile*tukeys_constant)
86474139
                 Outlier.Filter <- !is.na(value.window) &
                 diff.value.window >= fences &
                 diff.value.window > value.window
                 if(strict){
                     Outlier.Filter <- Outlier.Filter & (value.window > 0)
                 }
             }
             outliers<-which(Outlier.Filter)
             if(length(outliers)>0){
6a22d578
                 outliers <- lookup_window.range[outliers]
86474139
             }
             outliers
625325f1
         })
86474139
         outlier.vector.dups <- do.call(c,outlier.list)
         outlier.vector.uniq.sorted <- sort(unique(outlier.vector.dups))
         return(outlier.vector.uniq.sorted)
 }
6a22d578
 CreateDomainlist <- function(start.vector=NULL,end.vector=NULL,fill_gaps=NULL){
a47def9c
     Domains.by.start.list <- lapply(start.vector,function(x){
         data.frame(startbin=x, endbin=min(end.vector[end.vector > x]))
     })
     Domains.by.start.df <- do.call(rbind,Domains.by.start.list)
     Domains.by.start.df$gap.fill <- as.numeric(FALSE)
     Domains.by.start.df$level <- 2
     Domains.by.end.df <- NULL
     Domains.by.assumption.df <- NULL
6a22d578
     if(fill_gaps){
625325f1
         uncurated.ends <- end.vector[
         !(end.vector %in% Domains.by.start.df[,"endbin"])]
a47def9c
         if(length(uncurated.ends) > 0){
             Domains.by.end.list <- lapply(uncurated.ends,function(x){
625325f1
                 data.frame(startbin=(
                     max(Domains.by.start.df$endbin[
                         Domains.by.start.df$endbin<x])+1),endbin=x)
a47def9c
             })
             Domains.by.end.df <- do.call(rbind,Domains.by.end.list)
             Domains.by.end.df$gap.fill <- as.numeric(TRUE)
             Domains.by.end.df$level <- 1
         }
86474139
     }
625325f1
     All.Domains <- rbind(Domains.by.start.df,
         Domains.by.end.df,
         Domains.by.assumption.df)
a47def9c
     All.Domains.sorted <- All.Domains[order(All.Domains$startbin),]
     return(All.Domains.sorted)
 }
625325f1
 ComputeDirectionalityIndex <- function(Matrix = NULL, Window.size=NULL, 
     filter = NULL, start = NULL, end = NULL){
     Sequence <- seq_len(nrow(Matrix))
01a05500
     Sequence <- Sequence[start:end]
a47def9c
     DI.Data <- rep(NA,length(Sequence))
01a05500
     Bins.to.process <- Sequence[filter[start:end]]
625325f1
     All.bins <- seq_len(nrow(Matrix))
01a05500
     All.bins <- All.bins[filter]
     DI.list <- vapply(Bins.to.process,function(i){
86474139
             Upstream<-0
             Downstream<-0
             My.DI.Data <- NA
01a05500
             Relative.mid <- which(All.bins == i)
625325f1
             Window.range <- c(
                 (Relative.mid - Window.size) : (Relative.mid + Window.size))
             Window.range <- All.bins[
             Window.range[Window.range >= 1 & Window.range <= length(All.bins)]]
86474139
             Upstream.range <- Window.range[Window.range < i]
             Downstream.range <- Window.range[Window.range > i]
01a05500
             Row.vector <- Matrix[i,]
             Row.vector[is.na(Row.vector) | is.infinite(Row.vector)] <- 0
             if(length(Upstream.range) > 0){
                 Upstream <- sum(Row.vector[Upstream.range])
             }
             if(length(Downstream.range) > 0){
                 Downstream <- sum(Row.vector[Downstream.range])
             }
             Expected <- (Upstream + Downstream)/2
             if( Expected == 0 | Upstream == Downstream ){
                 My.DI.Data <- 0
86474139
             }else{
625325f1
                 My.DI.Data <- (
                     (Downstream - Upstream)/abs(Downstream - Upstream)
                     ) * (
                     ((Upstream - Expected)^2)/Expected + (
                         (Downstream - Expected)^2)/Expected)
86474139
             }
01a05500
         },1)
     DI.Data[Sequence %in% Bins.to.process] <- DI.list
86474139
     return(DI.Data)
 }
af47efdf
 get_directionality_index_by_chunks <- function(Brick = NULL, chr = NULL,
     resolution = NA, di_window = NULL, distance = NULL, chunk_size = 500, 
     sparse = FALSE, sparsity_threshold = 0.8, min_sum = -1, force = FALSE){
     Ranges <- Brick_get_bintable(Brick = Brick, chr = chr, 
         resolution = resolution)
     First.non.zero.bin <- ._get_first_nonzero_bin(Brick = Brick, chr = chr,
         resolution = resolution)
a47def9c
     chr.length <- length(Ranges)
632a1a64
     RowSums <- Brick_get_matrix_mcols(Brick = Brick, chr1 = chr, chr2 = chr, 
974f4033
         resolution = resolution, what = "chr1_row_sums")
a47def9c
     if(sparse){
632a1a64
         SparsityIndex <- Brick_get_matrix_mcols(Brick = Brick, chr1 = chr, 
af47efdf
             chr2 = chr, resolution = resolution, what = "sparsity")
86474139
     }
6a22d578
     if((chunk_size - (di_window*2))/di_window < 10){
         stop("chunk_size is too small for this di_window\n")
86474139
     }
6a22d578
     if(any(di_window > distance)){
         stop("di_window cannot be larger than distance\n")
86474139
     }
c1a95915
     Span <- (chr.length - First.non.zero.bin)
6a22d578
     Iterations <- Span/chunk_size
     Starts <- seq(from = First.non.zero.bin, to = chr.length, by = chunk_size)
c1a95915
     Starts <- Starts[Starts != chr.length]
     Ends <- c(Starts[-1] -1, chr.length)
a47def9c
     DI.data.list <- lapply(seq_along(Starts), function(x){
c1a95915
         Start <- Starts[x]
         End <- Ends[x]
a47def9c
         Position.start <- Start
         Position.end <- End
6a22d578
         Start <- ifelse((Start - di_window) < First.non.zero.bin, 
             First.non.zero.bin, (Start - di_window))
         End <- ifelse((End + di_window) > chr.length, 
             chr.length, (End + di_window))
a47def9c
         RowSums.subset <- RowSums[Start:End]
6a22d578
         Filter <- RowSums.subset > min_sum
a47def9c
         if(sparse){
             Sparsity.index.subset <- SparsityIndex[Start:End]
6a22d578
             Filter <- Filter & (Sparsity.index.subset > sparsity_threshold)
a47def9c
         }
         Total.length <- length(Filter)
         Filter.extend.length <- length(which(!Filter))
         True.length <- length(which(Filter))
625325f1
         extend <- Filter.extend.length
a47def9c
         while(True.length < Total.length){
625325f1
             Start <- ifelse((Start - extend) < First.non.zero.bin, 
                 First.non.zero.bin, Start - extend)
             End <- ifelse((End + extend) > chr.length, 
                 chr.length, (End + extend))
a47def9c
             RowSums.subset <- RowSums[Start:End]
6a22d578
             Filter <- RowSums.subset > min_sum
a47def9c
             if(sparse){
                 Sparsity.index.subset <- SparsityIndex[Start:End]
6a22d578
                 Filter <- Filter & (Sparsity.index.subset > sparsity_threshold)
a47def9c
             }
             True.length <- length(which(Filter))
             extend <- extend + 1
         }
632a1a64
         Matrix <- Brick_get_vector_values(Brick = Brick, chr1 = chr, 
af47efdf
             resolution = resolution, chr2 = chr, xaxis=c(Start:End), 
             yaxis=c(Start:End), force = force)
c1a95915
         # cat((Start - 1),"\n")
60fdd192
         # message(Position.start," ",Position.end,"\n")
         # message(Position.start - (Start - 1),Position.end - (Start - 1),"\n")
625325f1
         DI.data <- ComputeDirectionalityIndex(Matrix = Matrix, 
6a22d578
             Window.size = di_window, filter = Filter, 
625325f1
             start = Position.start - (Start - 1), 
             end = Position.end - (Start - 1))
a47def9c
         return(DI.data)
     })
     DI.data <- do.call(c, DI.data.list)
     DI.data <- c(rep(NA,First.non.zero.bin - 1),DI.data)
     Ranges$DI.Data <- DI.data
     return(Ranges)
 }
01a05500
 MakeBoundaries <- function(chr = NULL, Ranges = NULL, Binsize = NULL){
     Ends <- end(Ranges) 
625325f1
     Ends <- Ends[seq_len(length(Ends)-1)]
01a05500
     Starts <- start(Ranges)
625325f1
     Starts <- Starts[seq_len(length(Starts))[-1]] - 1 
01a05500
     Domain.boundaries <- unique(Starts,Ends)
632a1a64
     Boundary.ranges <- Brick_make_ranges(
6a22d578
         chrom=rep(chr,length(Domain.boundaries)),
         start=(Domain.boundaries-(Binsize/2))+1,
         end=Domain.boundaries+(Binsize/2))
01a05500
 }
 
 #' Do TAD Calls with Local Score Differentiator on a Hi-C matrix
 #' 
 #' `Local_score_differentiator` calls topologically associated domains on Hi-C 
 #' matrices. Local score differentiator at the most fundamental level is a 
 #' change point detector, which detects change points in the directionality 
 #' index using various thresholds defined on a local directionality index 
 #' distributions.
 #' The directionality index (DI) is calculated as defined by Dixon et al., 2012 
 #' Nature. Next, the difference of DI is calculated between neighbouring bins to
 #' get the change in DI distribution in each bin. When a DI value goes from a
 #' highly negative value to a highly positive one as expected to occur at domain
 #' boundaries, the ensuing DI difference distribution becomes a very flat 
 #' distribution interjected by very large peaks signifying regions where such
 #' a change may take place. We use two difference vectors, one is the difference
 #' vector between a bin and its adjacent downstream bin and another is the 
 #' difference between a bin and its adjacent upstream bin. Using these vectors,
 #' and the original directionality index, we define domain borders as outliers.
 #' 
 #' To define an outlier, fences are first defined. The fences are defined using
6a22d578
 #' tukeys_constant x inter-quartile range of the directionality index. The upper
01a05500
 #' fence used for detecting domain starts is the 75th quartile + 
6a22d578
 #' (IQR x tukeys_constant), while the lower fence is the 
 #' 25th quartile - (IQR x tukeys_constant). For domain starts the DI difference
01a05500
 #' must be greater than or equal to the upper fence, it must be greater than the
 #' DI and the DI must be a finite real value. If strict is TRUE, DI will also
 #' be required to be greater than 0. Similarly, for domain ends the 
 #' DI difference must be lower than or equal to the lower fence, it must be 
 #' lower than the DI and the DI must be a finite real value. If strict is TRUE,
 #' DI will also be required to be lower than 0. 
 #' 
 #' After defining outliers, each domain start will be associated to its 
6a22d578
 #' nearest downstream domain end. If \emph{fill_gaps} is defined as TRUE and
01a05500
 #' there are domain ends which remain unassociated to a domain start, These 
 #' domain ends will be associated to the bin adjacent to their nearest upstream
 #' domain end. This associations will be marked by metadata columns, gap.fill= 1
 #' and level = 1.
 #' 
 #' This function provides the capability to call very accurante TAD definitions
 #' in a very fast way. 
 #' 
632a1a64
 #' @inheritParams Brick_get_chrominfo
af47efdf
 #' @inheritParams Brick_add_ranges
01a05500
 #' 
 #' @param chrs \strong{Optional}. Default NULL
 #' If present, only TAD calls for elements in \emph{chrs} will be done.
 #' 
6a22d578
 #' @param min_sum \strong{Optional}. Default -1
 #' Process bins in the matrix with row.sums greater than \emph{min_sum}.
01a05500
 #' 
6a22d578
 #' @param di_window \strong{Optional}. Default 200
 #' Use \emph{di_window} to define the directionality index.
01a05500
 #' 
6a22d578
 #' @param lookup_window \strong{Optional}. Default 200
 #' Use \emph{lookup_window} local window to call borders. At smaller 
 #' \emph{di_window} values we recommend setting this to 2*\emph{di_window}
01a05500
 #' 
6a22d578
 #' @param tukeys_constant \strong{Optional}. Default 1.5
 #' \emph{tukeys_constant}*IQR (inter-quartile range) defines the lower and upper
01a05500
 #' fence values.
 #' 
 #' @param strict \strong{Optional}. Default TRUE
 #' If TRUE, \emph{strict} creates an additional filter on the directionality 
 #' index requiring it to be either greater than or less than 0 on the right tail
 #' or left tail respectively.  
 #' 
6a22d578
 #' @param fill_gaps \strong{Optional}. Default TRUE
01a05500
 #' If TRUE, this will affect the TAD stiching process. All Border starts are 
 #' stiched to the next downstream border ends. Therefore, at times border ends 
 #' remain unassociated to a border start. These border ends are stiched to the 
6a22d578
 #' adjacent downstream bin from their upstream border end when \emph{fill_gaps} 
01a05500
 #' is true. 
 #' 
 #' TADs inferred in this way will be annotated with two metadata columns in the 
 #' GRanges object. \emph{gap.fill} will hold a value of 1 and \emph{level} will 
625325f1
 #' hold a value 1. TADs which were not filled in will hold a gap.fill value of
 #' 0 and a level value of 2.
01a05500
 #' 
6a22d578
 #' @param ignore_sparse \strong{Optional}. Default TRUE
01a05500
 #' If TRUE, a matrix which has been defined as sparse during the matrix loading
6a22d578
 #' process will be treated as a dense matrix. The \emph{sparsity_threshold} 
01a05500
 #' filter will not be applied. Please note, that if a matrix is defined as 
6a22d578
 #' sparse and fill_gaps is TRUE, fill_gaps will be turned off.
01a05500
 #' 
6a22d578
 #' @param sparsity_threshold \strong{Optional}. Default 0.8
01a05500
 #' Sparsity threshold relates to the sparsity index, which is computed as the 
 #' number of non-zero bins at a certain distance from the diagonal. If a matrix
6a22d578
 #' is sparse and ignore_sparse is FALSE, bins which have a sparsity index value
01a05500
 #' below this threshold will be discarded from DI computation.
 #' 
6a22d578
 #' @param remove_empty Not implemented.
01a05500
 #' After implementation, this will ensure that the presence of centromeric 
 #' regions is accounted for.
 #' 
6a22d578
 #' @param chunk_size \strong{Optional}. Default 500
01a05500
 #' The size of the matrix chunk to process. This value should be larger than 2x
6a22d578
 #' di_window.
01a05500
 #' 
6a22d578
 #' @param force_retrieve \strong{Optional}. Default TRUE
01a05500
 #' If TRUE, this will force the retrieval of a matrix chunk even when the 
632a1a64
 #' retrieval includes interaction points which were not loaded into a Brick 
01a05500
 #' store (larger chunks). Please note, that this does not mean that DI can be 
 #' computed at distances larger than max distance. Rather, this is meant to aid
 #' faster computation.
 #' 
 #' @return A ranges object containing domain definitions. The starts and ends
 #' of the ranges coincide with the starts and ends of their contained bins from 
 #' the bintable. 
 #' 
625325f1
 #' @examples
fb79dd23
 #' Bintable.path <- system.file(file.path("extdata", "Bintable_100kb.bins"), 
 #' package = "HiCBricks")
6a22d578
 #' 
fb79dd23
 #' out_dir <- file.path(tempdir(), "lsd_test")
 #' dir.create(out_dir)
 #' 
 #' My_BrickContainer <- Create_many_Bricks(BinTable = Bintable.path, 
 #'   bin_delim = " ", output_directory = out_dir, file_prefix = "Test",
e813fdca
 #'   experiment_name = "Vignette Test", resolution = 100000,
fb79dd23
 #'   remove_existing = TRUE)
 #' 
 #' Matrix_file <- system.file(file.path("extdata", 
 #' "Sexton2012_yaffetanay_CisTrans_100000_corrected_chr3R.txt.gz"), 
 #' package = "HiCBricks")
 #' 
 #' Brick_load_matrix(Brick = My_BrickContainer, chr1 = "chr3R", 
e813fdca
 #' chr2 = "chr3R", matrix_file = Matrix_file, delim = " ",
fb79dd23
 #' remove_prior = TRUE, resolution = 100000)
 #' 
 #' TAD_ranges <- Brick_local_score_differentiator(Brick = My_BrickContainer, 
 #' chrs = "chr3R", resolution = 100000, di_window = 10, lookup_window = 30, 
 #' strict = TRUE, fill_gaps = TRUE, chunk_size = 500)
2c40cea7
 Brick_local_score_differentiator <- function(Brick, chrs = NULL, 
af47efdf
     resolution = NA, all_resolutions = FALSE, min_sum = -1, di_window = 200L, 
     lookup_window = 200L, tukeys_constant=1.5, strict = TRUE, fill_gaps=TRUE, 
     ignore_sparse=TRUE, sparsity_threshold=0.8, remove_empty = NULL, 
     chunk_size = 500, force_retrieve = TRUE){
     BrickContainer_resolution_check(resolution, all_resolutions)
     ChromInfo <- Brick_get_chrominfo(Brick = Brick, 
         resolution = resolution)
a47def9c
     Chromosomes <- ChromInfo[,'chr']
01a05500
     if(!is.null(chrs)){
a47def9c
         Chromosomes <- ChromInfo[ChromInfo[,'chr'] %in% chrs,'chr']
86474139
     }
8ebacdac
     chr_done_filter <- vapply(Chromosomes, function(chr){
         Brick_matrix_isdone(Brick = Brick, chr1 = chr, chr2 = chr, 
             resolution = resolution)
aa29c477
     }, TRUE)
8ebacdac
     if(!all(chr_done_filter)){
         message("Skipping intra-chromosomal maps containing no data...")
         message(paste(Chromosomes[!chr_done_filter], collapse = ", "), 
             " will be skipped")
     }
a47def9c
     Chrom.domains.ranges.list <- lapply(Chromosomes, function(chr){
af47efdf
         Ranges <- Brick_get_bintable(Brick = Brick, chr = chr, 
             resolution = resolution)
         sparse <- Brick_matrix_issparse(Brick = Brick, chr1 = chr, chr2 = chr, 
             resolution = resolution)
632a1a64
         max.distance <- Brick_matrix_maxdist(Brick = Brick, chr1 = chr, 
af47efdf
             chr2 = chr, resolution = resolution)
6a22d578
         if(ignore_sparse){
a47def9c
             sparse=FALSE
         }
6a22d578
         if(sparse & fill_gaps){
             fill_gaps=FALSE
a47def9c
         }
60fdd192
         message("[1] Computing DI for ",chr,"\n")
af47efdf
         Ranges <- get_directionality_index_by_chunks(Brick = Brick, 
             chr = chr, 
             resolution = resolution,
6a22d578
             di_window = di_window, 
af47efdf
             distance = max.distance, 
             chunk_size = chunk_size, 
             sparse=sparse, 
6a22d578
             sparsity_threshold=sparsity_threshold,
             min_sum = min_sum, force = force_retrieve)
01a05500
 
632a1a64
         RowSums <- Brick_get_matrix_mcols(Brick = Brick, chr1 = chr, 
974f4033
             chr2 = chr, resolution = resolution, what = "chr1_row_sums")
a47def9c
         Ranges$row.sums <- RowSums
60fdd192
         message("[2] Computing DI Differences for ",chr,"\n")
a47def9c
         if(sparse){
af47efdf
             SparsityIndex <- Brick_get_matrix_mcols(Brick = Brick, 
                 chr1 = chr, 
                 chr2 = chr, 
                 resolution = resolution, 
                 what = "sparsity")
625325f1
             Backwards.DI.Difference <- Backwards.Difference(
af47efdf
                 Vector = Ranges$DI.Data, 
                 sparse = sparse,
625325f1
                 sparsity.idx = SparsityIndex, 
6a22d578
                 sparsity_threshold = sparsity_threshold)
625325f1
             Forwards.DI.Difference <- Forwards.Difference(
af47efdf
                 Vector = Ranges$DI.Data, 
                 sparse = sparse,
625325f1
                 sparsity.idx = SparsityIndex, 
6a22d578
                 sparsity_threshold = sparsity_threshold)
a47def9c
         }else{
625325f1
             Backwards.DI.Difference <- Backwards.Difference(
                 Vector=Ranges$DI.Data)
             Forwards.DI.Difference <- Forwards.Difference(
                 Vector=Ranges$DI.Data)
a47def9c
         }
         Ranges$backward.Differences <- Backwards.DI.Difference
         Ranges$forward.Differences <- Forwards.DI.Difference
d5cc8357
         message("[2] Done\n")
         message("[3] Fetching Outliers ",chr,"\n")
a47def9c
         if(sparse){
625325f1
             Domain.end.candidates <- ComputeOutlierOverIQRInsideWindow(
af47efdf
                 lookup_window = lookup_window,
                 diff.values = Backwards.DI.Difference, 
                 values = Ranges$DI.Data, 
                 sparse = sparse, 
625325f1
                 row.sums = Ranges$row.sums,
af47efdf
                 min_sum = min_sum, 
                 sparsity.idx = SparsityIndex, 
                 sparsity_threshold = sparsity_threshold, 
                 tukeys_constant = tukeys_constant, 
                 tail = "lower.tail",
                 strict = strict)
625325f1
             Domain.start.candidates <- ComputeOutlierOverIQRInsideWindow(
af47efdf
                 lookup_window = lookup_window,
                 diff.values = Forwards.DI.Difference, 
                 values = Ranges$DI.Data, 
                 sparse = sparse, 
                 row.sums = Ranges$row.sums,
                 min_sum = min_sum, 
                 sparsity.idx = SparsityIndex, 
                 sparsity_threshold = sparsity_threshold,
                 tukeys_constant = tukeys_constant, 
                 tail = "upper.tail", 
                 strict = strict)
a47def9c
         }else{
625325f1
             Domain.end.candidates <- ComputeOutlierOverIQRInsideWindow(
6a22d578
                 lookup_window=lookup_window,
625325f1
                 diff.values=Backwards.DI.Difference,
                 values=Ranges$DI.Data, 
                 row.sums = Ranges$row.sums,
6a22d578
                 min_sum = min_sum, 
                 tukeys_constant=tukeys_constant,
625325f1
                 tail="lower.tail",strict=strict)
             Domain.start.candidates <- ComputeOutlierOverIQRInsideWindow(
6a22d578
                 lookup_window=lookup_window,
625325f1
                 diff.values=Forwards.DI.Difference,
                 values=Ranges$DI.Data, 
                 row.sums = Ranges$row.sums,
6a22d578
                 min_sum = min_sum, 
                 tukeys_constant=tukeys_constant,
625325f1
                 tail="upper.tail",
                 strict=strict)
a47def9c
         }
625325f1
         Domain.start.candidates <- Domain.start.candidates[
         Domain.start.candidates != length(Ranges)]
         Domain.end.candidates <- Domain.end.candidates[
         Domain.end.candidates != 1]
d5cc8357
         message("[3] Done\n")
60fdd192
         message("[4] Creating Domain list for ",chr,"\n")
a47def9c
 
         if(!(1 %in% Domain.start.candidates)){
             Domain.start.candidates <- c(1,Domain.start.candidates)
         }
         if(!(length(Ranges) %in% Domain.end.candidates)){
             Domain.end.candidates <- c(Domain.end.candidates,length(Ranges))
         }
         Domain.list <- CreateDomainlist(start.vector=Domain.start.candidates,
af47efdf
             end.vector=Domain.end.candidates, fill_gaps=fill_gaps)
6a22d578
         Domain.Ranges <- Brick_make_ranges(chrom=rep(chr,nrow(Domain.list)),
             start=start(Ranges[Domain.list$startbin]),
             end=end(Ranges[Domain.list$endbin]))
d5cc8357
         message("[4] Done\n")
a47def9c
         Domain.Ranges$gap.fill <- Domain.list$gap.fill
         Domain.Ranges$level <- Domain.list$level
6a22d578
         Domain.Ranges$window.size <- di_window
         Domain.Ranges$lookup_window <- lookup_window
a47def9c
         return(Domain.Ranges)
     })
625325f1
     Chrom.domains.ranges <- do.call(c,unlist(Chrom.domains.ranges.list, 
         use.names = TRUE))
01a05500
     return(Chrom.domains.ranges)
d5cc8357
 }