LBMC » NucleoMiner

# get_content = function(# Get content from cached 

# ### Acces to the cached content of a file via the global variable NM_CACHE. Load content in NM_CACHE if needed.

# obj, ##<< The object that we want its content.

# ... ##<< Parameters that will be passed to my_read

# ) {

#   UseMethod("get_content", obj)

#   ### Returns the cached content of the current object. 

# }

# 

# get_content.default = structure( function(# Get content from cached 

# ### Acces to the cached content of a file via the global variable NM_CACHE. Load content in NM_CACHE if needed.

# obj, ##<< The object that we want its content.

# ... ##<< Parameters that will be passed to my_read

# ) {

#   if (inherits(try(NM_CACHE,TRUE), "try-error") || is.null(NM_CACHE)) {

#     NM_CACHE <<- list()

#   }

#   if (is.null(NM_CACHE[[obj$filename]])) {

#     print(paste("Loading file ",obj$filename, sep=""))

#     tmp_content = my_read(obj, ...)

#     print("affect it...")

#     NM_CACHE[[obj$filename]] <<- tmp_content

#     print("done.")

#   }  

#   return(NM_CACHE[[obj$filename]])

# ### Returns the cached content of the current object. 

# }, ex=function(){

#   # Create a dataframe

#   df = NULL

#   df = dfadd(df, list(key1 = "value1", key2 = "value2"))

#   df = dfadd(df, list(key1 = "value1'", key2 = "value3'"))

#   # Dump it into tmp file

#   write.table(df, file="/tmp/tmp.dump.table.tmp")

#   # Load it into cache using feature.

#   cached_content_object  = list(filename="/tmp/tmp.dump.table.tmp")

#   class(cached_content_object) = "table"  

#   # First time it will be load into cache

#   print(get_content(cached_content_object))

#   # Second time not

#   print(get_content(cached_content_object))

#   

# })

#   

# my_read = function(

# ### Abstract my_read function.  

#   obj, ...) {

#   UseMethod("my_read", obj)

# }

#   

# my_read.default = function(

# ### Default my_read function.  

#   obj, ...){

#   stop(paste("ERROR, my_read is not defined for any Objects (file ", obj$filename," )", sep=""))      

# }

# 

# my_read.fasta = function(

# ### my_read function for fasta files.  

#   obj, ...){

#   # require(seqinr)

#   return(read.fasta(obj$filename, ...))

# }

# 

# my_read.table = function(

# ### my_read function for table files.  

#   obj, ...){

#   if (rev(unlist(strsplit(obj$filename, ".", fixed=TRUE)))[1] == "gz") {

#     return(read.table(file=gzfile(obj$filename), ...))      

#   } else { 

#     return(read.table(file=obj$filename, ...))

#   }

# }  

# 

# my_read.cvs = function(

# ### my_read function for cvs files.  

#   obj, ...){

#   return(read.csv(file=obj$filename, ...))

# }  

FDR = structure(function#  False Discovery Rate

### From a vector x of independent p-values, extract the cutoff corresponding to the specified FDR. See Benjamini & Hochberg 1995 paper

##author<< Gael Yvert, 

(

x, ##<< A vector x of independent p-values.

FDR ##<< The specified FDR.

) {

  x <- sort(na.omit(x))

  N = length(x)

  i = 1;

  while(N*x[i]/i < FDR & i <= N) i = i + 1; # we search for the highest i where Nrandom / Nobserved < FDR

  if (i == 1)

    return (NA)

  else

    return( x[i-1] )

### Return the the corresponding cutoff.

}, ex=function(){

  print("example")

})

lod_score_vecs = structure(function # Likelihood ratio

### Compute the likelihood log of two set of value from two models Vs. a unique model.

(

x ,##<< First vector.

y ##<< Second vector.

) {

	if (length(x) <=1 | length(y) <= 1) {

		return(NA)

	}

  meanX = mean(x)

  sdX = sd(x)

  meanY = mean(y)

  sdY = sd(y)

  meanXY = mean(c(x,y))

  sdXY = sd(c(x,y))

  llX = sum(log(dnorm(x,mean=meanX,sd=sdX)))

  llY = sum(log(dnorm(y,mean=meanY,sd=sdY)))

  llXY = sum(log(dnorm(c(x,y),mean=meanXY,sd=sdXY)))

  ratio = llX + llY - llXY

  return(ratio)

### Returns the likelihood ratio.

}, ex=function(){

  # LOD score for 2 set of values

  mean1=5; sd1=2; card2 = 250

  mean2=6; sd2=3; card1 = 200

  x1 = rnorm(card1, mean1, sd1)

  x2 = rnorm(card2, mean2, sd2)  

  min = floor(min(c(x1,x2)))

  max = ceiling(max(c(x1,x2)))

  hist(c(x1,x2), xlim=c(min, max), breaks=min:max)

  lines(min:max,dnorm(min:max,mean1,sd1)*card1,col=2)

  lines(min:max,dnorm(min:max,mean2,sd2)*card2,col=3)

  lines(min:max,dnorm(min:max,mean(c(x1,x2)),sd(c(x1,x2)))*card2,col=4)

  lod_score_vecs(x1,x2)

 })

dfadd = structure(function# Adding list to a dataframe.

### Add a list \emph{l} to a dataframe \emph{df}. Create it if \emph{df} is \emph{NULL}. Return the dataframe \emph{df}.

	(df, ##<<  A dataframe

		l ##<<  A list

	) {

  if (is.null(df)) {

    df = data.frame(l,stringsAsFactors=FALSE)

  } else {

    df = rbind(df, data.frame(l,stringsAsFactors=FALSE))

  }

  return(df)

### Return the dataframe \emph{df}.

}, ex=function(){

		## Here dataframe is NULL

		print(df)

		df = NULL

		# Initialize df

		df = dfadd(df, list(key1 = "value1", key2 = "value2"))

		print(df)

		# Adding elements to df

		df = dfadd(df, list(key1 = "value1'", key2 = "value2'"))

		print(df)

})

sign_from_strand = function(

### Get the sign of strand 

strands) {

	apply(t(strands), 2, function(strand) {	if (strand == "F") return(1) else return(-1)})

### If strand in forward then returns 1 else returns -1

}

flat_reads = function(

### Extract reads coordinates from TempleteFilter input sequence

reads, ##<< TemplateFilter input reads 

nuc_width ##<< Width used to shift F and R reads.

) {

	F_flatted_reads = unlist(apply(t(reads[reads$V3=="F",]),2,function(r){rep(as.integer(r[2]), r[4])}))

	R_flatted_reads = unlist(apply(t(reads[reads$V3=="R",]),2,function(r){rep(as.integer(r[2]), r[4])}))

	flatted_reads = c(F_flatted_reads + rep(nuc_width/2, length(F_flatted_reads)), R_flatted_reads - rep(nuc_width/2, length(R_flatted_reads))  )

	return(list(F_flatted_reads, R_flatted_reads, flatted_reads))

### Returns a list of F reads, R reads and joint/shifted F and R reads.

}

filter_tf_outputs = function(# Filter TemplateFilter outputs

### This function filters TemplateFilter outputs according, not only genome area observerved properties, but also correlation and overlap threshold.

tf_outputs, ##<< TemplateFilter outputs.

chr, ##<< Chromosome observed, here chr is an integer.

x_min, ##<< Coordinate of the first bp observed.

x_max, ##<< Coordinate of the last bp observed.

nuc_width = 160, ##<< Nucleosome width.

ol_bp = 59, ##<< Overlap Threshold.

corr_thres = 0.5 ##<< Correlation threshold.

) {

  if (x_min < 0) {

    tf_outputs = tf_outputs[tf_outputs$chr == paste("chr", chr, sep="") & tf_outputs$center > (-x_max - nuc_width/2) & tf_outputs$center <  (-x_min + nuc_width/2),]

	} else {

    tf_outputs = tf_outputs[tf_outputs$chr == paste("chr", chr, sep="") & tf_outputs$center > (x_min - nuc_width/2) & tf_outputs$center < (x_max + nuc_width/2),]

  }

  tf_outputs$lower_bound = tf_outputs$center - tf_outputs$width/2

  tf_outputs$upper_bound = tf_outputs$center + tf_outputs$width/2

  tf_outputs = tf_outputs[tf_outputs$correlation >= corr_thres,]  

  tf_outputs = tf_outputs[order(tf_outputs$correlation,decreasing=TRUE),]  

  i = 1

  while (i <= length(tf_outputs[,1])) {    

    lb = tf_outputs[i,]$low 

    ub = tf_outputs[i,]$up

    tf_outputs = tf_outputs[!(tf_outputs$low <= (ub-ol_bp) & tf_outputs$up > ub) & !(tf_outputs$up >= (lb+ol_bp) & tf_outputs$low < lb),]

    i = i+1

  }

  return(tf_outputs)

### Returns filtered TemplateFilter Outputs

}

filter_tf_inputs = function(# Filter TemplateFilter inputs

### This function filters TemplateFilter inputs according genome area observed properties. It takes into account reads that are at the frontier of this area and the strand of these reads.

inputs, ##<< TF inputs to be filtered.

chr, ##<< Chromosome observed, here chr is an integer.

x_min, ##<< Coordinate of the first bp observed.

x_max, ##<< Coordinate of the last bp observed.

nuc_width = 160, ##<< Nucleosome width.

only_f = FALSE, ##<< Filter only F reads.

only_r = FALSE, ##<< Filter only R reads.

filter_for_coverage = FALSE ##<< Does it filter for plot coverage?

) {

	if (only_f) {

		inputs = inputs[inputs[,1]==chr & inputs[,2] >= x_min - nuc_width & inputs[,3] == "F" & inputs[,2] <= x_max + nuc_width,]

	} else if (only_r) {

		inputs = inputs[inputs[,1]==chr & inputs[,2] >= x_min - nuc_width & inputs[,3] == "R" & inputs[,2] <= x_max + nuc_width,]			

	} else {

		inputs = inputs[inputs[,1]==chr & inputs[,2] >= x_min - nuc_width & inputs[,2] <= x_max + nuc_width,]

	}

  if (!filter_for_coverage) {

    corrected_inputs_coords = inputs[,2] + nuc_width/2 * sign_from_strand(inputs[,3])

    inputs = inputs[inputs[,1]==chr & corrected_inputs_coords >= x_min & corrected_inputs_coords <= x_max,]    

  }

	return(inputs)

### Returns filtred inputs.

}

get_comp_strand = function(

### Compute the complementatry strand.

strand ##<< The original strand.

) {

	apply(t(strand),2, function(n){

	  if (n=="a") {return("t")} 

		if (n=="t") {return("a")}

		if (n=="c") {return("g")} 

		if (n=="g") {return("c")} 

	})

### Returns the complementatry strand.

}

aggregate_intra_strain_nucs = structure(function(# Aggregate replicated sample's nucleosomes. 

### This function aggregates nucleosome for replicated samples. It uses TemplateFilter ouput of each sample as replicate. Each sample owns a set of nucleosomes computed using TemplateFilter and ordered by the position of their center. Adajacent nucleosomes are compared two by two. Comparison is based on a log likelihood ratio score. The issue of comparison is adjacents nucleosomes merge or separation. Finally the function returns a list of clusters and all computed \emph{lod_scores}. Each cluster ows an attribute \emph{wp} for "well positionned". This attribute is set as \emph{TRUE} if the cluster is composed of exactly one nucleosomes of each sample.

samples, ##<< A list of samples. Each sample is a list like \emph{sample = list(id=..., marker=..., strain=..., roi=..., inputs=..., outputs=...)} with \emph{roi = list(name=..., begin=...,  end=..., chr=..., genome=...)}.

lod_thres=20, ##<< Log likelihood ration threshold.

coord_max=20000000 ##<< A too big value to be a coord for a nucleosome lower bound.

){

	end_of_tracks = function(tracks) {

		if (length(tracks) == 0) {

			return(TRUE)

		}

	  for (lower_bound in tracks) {

			if(!is.na(lower_bound)) {

	      if (lower_bound < coord_max) {

	        return(FALSE)

	      }

	  	}

	  }

	  return(TRUE)

	}

	store_cluster = function(clusters, new_cluster, nb_nucs_in_cluster,nuc_from_track,nb_tracks, min_nuc_center, max_nuc_center) {

		if ( nb_nucs_in_cluster==nb_tracks & sum(nuc_from_track)==nb_tracks) {

			new_cluster$wp = TRUE

			center = (new_cluster$lower_bound + new_cluster$upper_bound) / 2

			if (is.null(min_nuc_center) | ((min_nuc_center <= center) & (center < max_nuc_center))) {

		  	clusters[[length(clusters) + 1]] = new_cluster

				# print(new_cluster)

		  }

		} else {

			new_cluster$wp = FALSE

			center = (new_cluster$lower_bound + new_cluster$upper_bound) / 2

			if (is.null(min_nuc_center) | ((min_nuc_center <= center) & (center < max_nuc_center))) {

			  clusters[[length(clusters) + 1]] = new_cluster

			}

		}

		return(clusters)

	}

	strain = samples[[1]]$strain

	lod_scores = c()

  min_nuc_center = min(samples[[1]]$roi$begin, samples[[1]]$roi$end)

	max_nuc_center = max(samples[[1]]$roi$begin, samples[[1]]$roi$end)

  # compute clusters

  clusters = list()

  cluster_contents = list()

  # Init reader

  indexes = c()

  track_readers = c()

  current_nuc = NULL

	lod_score = lod_thres + 1 

  # Read nucs from TF outputs  

  tf_outs = list()

	i = 1

  for (sample in samples) {

		# print(sample$roi$chr)

		# print(min_nuc_center)

		# print(max_nuc_center)

		# print(sample$outputs)

		# tf_outs[[i]] = filter_tf_outputs(sample$outputs, sample$roi$chr, min_nuc_center, max_nuc_center)

		# print(tf_outs[[i]])

		tf_outs[[i]] = sample$outputs

		tf_outs[[i]] = tf_outs[[i]][order(tf_outs[[i]]$center),]  

    indexes[i] = 1

		if (is.na(tf_outs[[i]][indexes[i],]$center)) {

      track_readers[i] = coord_max				

	  } else {

      track_readers[i] = tf_outs[[i]][indexes[i],]$center				

		}

		i = i+1		

  }

	# print(track_readers)

  new_cluster = NULL

  nb_nucs_in_cluster = 0

  nuc_from_track = c()

  for (i in 1:length(tf_outs)){

    nuc_from_track[i] = FALSE

  }

  # Start clustering

  while (!end_of_tracks(track_readers)) {

    new_center = min(track_readers)

		current_track = which(track_readers == new_center)[1]			

    new_nuc = as.list(tf_outs[[current_track]][indexes[current_track],])

		new_nuc$chr = substr(new_nuc$chr,4,1000000L)

		new_nuc$inputs = samples[[current_track]]$inputs

		new_nuc$chr = samples[[current_track]]$roi$chr

		new_nuc$track = current_track

		new_nuc$inputs = filter_tf_inputs(samples[[current_track]]$inputs, new_nuc$chr, new_nuc$lower_bound, new_nuc$upper_bound, new_nuc$width) 

		flatted_reads = flat_reads(new_nuc$inputs, new_nuc$width)

		new_nuc$original_reads = flatted_reads[[3]]

    new_upper_bound = new_nuc$upper_bound

    if (!is.null(current_nuc)) {

			lod_score = lod_score_vecs(current_nuc$original_reads,new_nuc$original_reads)

			lod_scores = c(lod_scores,lod_score)

		}

		# print(paste(lod_score, length(current_nuc$original_reads), length(new_nuc$original_reads), sep=" "))

		if (is.na(lod_score)) {

			lod_score = lod_thres + 1 

		}

		# Store lod_score

		new_nuc$lod_score = lod_score 			

	  if (lod_score < lod_thres) {

      # aggregate to current cluster

      #   update bound

      if (new_nuc$upper_bound > new_cluster$upper_bound) {

        new_cluster$upper_bound = new_nuc$upper_bound      

      }

      if (new_nuc$lower_bound < new_cluster$lower_bound) {

        new_cluster$lower_bound = new_nuc$lower_bound         

      }

      #   add nucleosome to current cluster

      nuc_from_track[current_track] = TRUE

      nb_nucs_in_cluster = nb_nucs_in_cluster + 1 

			new_cluster$nucs[[length(new_cluster$nucs)+1]] = new_nuc

    } else {

			if (!is.null(new_cluster)) {

        # store old cluster

	      clusters = store_cluster(clusters, new_cluster, nb_nucs_in_cluster,nuc_from_track,length(tf_outs),min_nuc_center, max_nuc_center)

			}

      # Reinit current cluster composition stuff

      nb_nucs_in_cluster = 0

      nuc_from_track = c()

      for (i in 1:length(tf_outs)){

        nuc_from_track[i] = FALSE

      }

      # create new cluster        

      new_cluster = list(lower_bound=new_nuc$low, upper_bound=new_nuc$up, chr=new_nuc$chr, strain_ref=strain , nucs=list())

      # update upper bound

      current_upper_bound = new_upper_bound

      # add nucleosome to current cluster

      nb_nucs_in_cluster = nb_nucs_in_cluster + 1 

      nuc_from_track[current_track] = TRUE

			new_cluster$nucs[[length(new_cluster$nucs)+1]] = new_nuc

		}

		current_nuc = new_nuc

    # update indexes

    if (indexes[current_track] < length(tf_outs[[current_track]]$center)) {      

      indexes[current_track] = indexes[current_track] + 1

      # update track

      track_readers[current_track] = tf_outs[[current_track]][indexes[current_track],]$center

    } else {

      # update track

      track_readers[current_track] = coord_max

    }

  }

  # store last cluster

  if (!is.null(new_cluster)) {

    # store old cluster

    clusters = store_cluster(clusters, new_cluster, nb_nucs_in_cluster,nuc_from_track,length(tf_outs),min_nuc_center, max_nuc_center)

  }

	return(list(clusters, lod_scores))

### Returns a list of clusterized nucleosomes, and all computed lod scores.

}, ex=function(){	

	# Dealing with a region of interest

	roi =list(name="example", begin=1000,  end=1300, chr="1", genome=rep("A",301))

	samples = list()

	for (i in 1:3) {

		# Create TF output

		tf_nuc = list("chr"=paste("chr", roi$chr, sep=""), "center"=(roi$end + roi$begin)/2, "width"= 150, "correlation.score"= 0.9)

		outputs = dfadd(NULL,tf_nuc)

		outputs = filter_tf_outputs(outputs, roi$chr, roi$begin, roi$end)

		# Generate corresponding reads

		nb_reads = round(runif(1,170,230))

		reads = round(rnorm(nb_reads, tf_nuc$center,20))

		u_reads = sort(unique(reads))

		strands = sample(c(rep("R",ceiling(length(u_reads)/2)),rep("F",floor(length(u_reads)/2))))

		counts = apply(t(u_reads), 2, function(r) { sum(reads == r)})

		shifts = apply(t(strands), 2, function(s) { if (s == "F") return(-tf_nuc$width/2) else return(tf_nuc$width/2)})

		u_reads = u_reads + shifts

		inputs = data.frame(list("V1" = rep(roi$chr, length(u_reads)), 

		                         "V2" = u_reads, 

														 "V3" = strands, 

														 "V4" = counts), stringsAsFactors=FALSE)

		samples[[length(samples) + 1]] = list(id=1, marker="Mnase_Seq", strain="strain_ex", total_reads = 10000000, roi=roi, inputs=inputs, outputs=outputs)

	}

	print(aggregate_intra_strain_nucs(samples))

})

flat_aggregated_intra_strain_nucs = function(# to flat aggregate_intra_strain_nucs function output

### This function builds a dataframe of all clusters obtain from aggregate_intra_strain_nucs function.

partial_strain_maps, ##<< the output of aggregate_intra_strain_nucs function

roi_index ##<< the index of the roi involved

) {

	if  (length(partial_strain_maps) == 0 ){

		print(paste("Empty partial_strain_maps for roi", roi_index, "ands strain", strain, "." ))

    tmp_strain_maps = list()  

	} else {

		tmp_strain_map = apply(t(1:length(partial_strain_maps)), 2, function(i){

			tmp_nuc = partial_strain_maps[[i]]

			tmp_nuc_as_list = list()

			tmp_nuc_as_list[["chr"]] = tmp_nuc[["chr"]]

			tmp_nuc_as_list[["lower_bound"]] = ceiling(tmp_nuc[["lower_bound"]])

			tmp_nuc_as_list[["upper_bound"]] = floor(tmp_nuc[["upper_bound"]])

			tmp_nuc_as_list[["roi_index"]] = roi_index

			tmp_nuc_as_list[["index_nuc"]] = i

			tmp_nuc_as_list[["wp"]] = as.integer(tmp_nuc$wp)

			all_original_reads = c()

			for (j in 1:length(tmp_nuc$nucs)) {

				all_original_reads = c(all_original_reads, tmp_nuc$nucs[[j]]$original_reads)

			}

			tmp_nuc_as_list[["nb_reads"]] = length(all_original_reads)

			if (tmp_nuc$wp) {

				tmp_nuc_as_list[["lod_1"]] = signif(tmp_nuc$nucs[[2]]$lod_score,5)

				tmp_nuc_as_list[["lod_2"]] = signif(tmp_nuc$nucs[[3]]$lod_score,5)

			} else {

				tmp_nuc_as_list[["lod_1"]] = NA

				tmp_nuc_as_list[["lod_2"]] = NA							

			}

      return(tmp_nuc_as_list)

    })

    tmp_strain_maps = do.call("rbind", tmp_strain_map)

	}

  return(data.frame(tmp_strain_maps))

### Returns a dataframe of all clusters obtain from aggregate_intra_strain_nucs function.

}

align_inter_strain_nucs = structure(function(# Aligns nucleosomes between 2 strains.

### This function aligns nucs between two strains for a given genome region.

replicates, ##<< Set of replicates, ideally 3 per strain. 

wp_nucs_strain_ref1=NULL, ##<< List of aggregates nucleosome for strain 1. If it's null this list will be computed.

wp_nucs_strain_ref2=NULL, ##<< List of aggregates nucleosome for strain 2. If it's null this list will be computed.

corr_thres=0.5, ##<< Correlation threshold.

lod_thres=100, ##<< LOD cut off.

config=NULL, ##<< GLOBAL config variable

... ##<< A list of parameters that will be passed to \emph{aggregate_intra_strain_nucs} if needed.	

) {

	if (length(replicates) < 2) {

		stop("ERROR, align_inter_strain_nucs needs 2 replicate sets.")

	} else if (length(replicates) > 2) {

		print("WARNING, align_inter_strain_nucs will use 2 first sets of replicates as inputs.")			

	}

	common_nuc = NULL

	lod_scores = c()

	chr = replicates[[1]][[1]]$roi$chr

  min_nuc_center = min(replicates[[1]][[1]]$roi$begin, replicates[[1]][[1]]$roi$end)

	max_nuc_center = max(replicates[[1]][[1]]$roi$begin, replicates[[1]][[1]]$roi$end)

	strain_ref1 = replicates[[1]][[1]]$strain

	strain_ref2 = replicates[[2]][[1]]$strain

	big_roi = replicates[[1]][[1]]$roi

  orig_big_roi = replicates[[1]][[1]]$orig_roi

	if(big_roi$end - big_roi$begin < 0) {

		tmp_begin = big_roi$begin

		big_roi$begin =  big_roi$end

		big_roi$end =  tmp_begin

	}

	# GO!

	if (is.null(wp_nucs_strain_ref1)) {

		wp_nucs_strain_ref1 = aggregate_intra_strain_nucs(replicates[[1]], ...)[[1]]

	}

	if (is.null(wp_nucs_strain_ref2)) {

	  wp_nucs_strain_ref2 = aggregate_intra_strain_nucs(replicates[[2]], ...)[[1]]

  }

  # foo <<- wp_nucs_strain_ref1

  # print(apply(t(wp_nucs_strain_ref1), 2, function(l){c(l[[1]]$lower_bound, l[[1]]$upper_bound, l[[1]]$wp)}))

  # print(apply(t(wp_nucs_strain_ref2), 2, function(l){c(l[[1]]$lower_bound, l[[1]]$upper_bound, l[[1]]$wp)}))

	# dealing with matching_nas

	lws = c()

	ups = c()

	for (na in wp_nucs_strain_ref2) {

		lws = c(lws, na$lower_bound)

		ups = c(ups, na$upper_bound)

	}

	print(paste("Exploring chr" , chr , ", " , length(wp_nucs_strain_ref1) , ", [" , min_nuc_center , ", " , max_nuc_center , "] nucs...", sep=""))			

	roi_strain_ref1 = NULL

	roi_strain_ref2 = NULL

	if (length(wp_nucs_strain_ref1) > 0) {

		for(index_nuc_strain_ref1 in 1:length(wp_nucs_strain_ref1)){

			# print(paste("" , index_nuc_strain_ref1 , "/" , length(wp_nucs_strain_ref1), sep=""))

			nuc_strain_ref1 = wp_nucs_strain_ref1[[index_nuc_strain_ref1]]

			# Filtering on Well Positionned

			if (nuc_strain_ref1$wp) {

				roi_strain_ref1 = list(name=paste("strain_chr_id_" , strain_ref1 , "_" , chr , "_" , "i" , "_", sep=""), begin=nuc_strain_ref1$lower_bound, end=nuc_strain_ref1$upper_bound, chr=chr, strain_ref = strain_ref1)

				roi_strain_ref2 = translate_roi(roi_strain_ref1, strain_ref2, big_roi=orig_big_roi, config=config) 

        if (!is.null(roi_strain_ref2)){

					# LOADING INTRA_STRAIN_NUCS_FILENAME_STRAIN_REF2 FILE(S) TO COMPUTE MATCHING_NAS (FILTER)

					lower_bound_roi_strain_ref2 = min(roi_strain_ref2$end,roi_strain_ref2$begin)

					upper_bound_roi_strain_ref2 = max(roi_strain_ref2$end,roi_strain_ref2$begin)

					matching_nas = which( lower_bound_roi_strain_ref2 <= ups & lws <= upper_bound_roi_strain_ref2)

					for (index_nuc_strain_ref2 in matching_nas) {

						nuc_strain_ref2 = wp_nucs_strain_ref2[[index_nuc_strain_ref2]]

						# Filtering on Well Positionned

    				nuc_strain_ref2_to_roi = list(begin=nuc_strain_ref2$lower_bound, end=nuc_strain_ref2$upper_bound, chr=nuc_strain_ref2$chr, strain_ref = strain_ref2)             

						if (!is.null(translate_roi(nuc_strain_ref2_to_roi, strain_ref1, big_roi=orig_big_roi, config=config)) &  

                nuc_strain_ref2$wp) {

							# Filtering on correlation Score and collecting reads					

							SKIP = FALSE

							# TODO: This for loop could be done before working on strain_ref2. Isn't it?

							reads_strain_ref1 = c()

							for (nuc in nuc_strain_ref1$nucs){

								reads_strain_ref1 = c(reads_strain_ref1, nuc$original_reads)							

								if (nuc$corr < corr_thres) {

									SKIP = TRUE

								}

							}

							reads_strain_ref2 = c()

							for (nuc in nuc_strain_ref2$nucs){

								reads_strain_ref2 = c(reads_strain_ref2, nuc$original_reads)							

								if (nuc$corr < corr_thres) {

									SKIP = TRUE

								}

							}

							# Filtering on correlation Score						

							if (!SKIP) {

								# tranlation of reads into strain 2 coords

								diff = ((roi_strain_ref1$begin + roi_strain_ref1$end) - (roi_strain_ref2$begin + roi_strain_ref2$end)) / 2

								reads_strain_ref1 = reads_strain_ref1 - rep(diff, length(reads_strain_ref1))

								lod_score = lod_score_vecs(reads_strain_ref1, reads_strain_ref2)

								lod_scores = c(lod_scores, lod_score)

								# Filtering on LOD Score						

								if (lod_score < lod_thres) {

									tmp_nuc = list()

									# strain_ref1

									tmp_nuc[[paste("chr_", strain_ref1, sep="")]] = chr

									tmp_nuc[[paste("lower_bound_", strain_ref1, sep="")]] = nuc_strain_ref1$lower_bound

									tmp_nuc[[paste("upper_bound_", strain_ref1, sep="")]] = nuc_strain_ref1$upper_bound

									tmp_nuc[[paste("mean_", strain_ref1, sep="")]] = signif(mean(reads_strain_ref1),5)

									tmp_nuc[[paste("sd_", strain_ref1, sep="")]] = signif(sd(reads_strain_ref1),5)

									tmp_nuc[[paste("nb_reads_", strain_ref1, sep="")]] = length(reads_strain_ref1)

									tmp_nuc[[paste("index_nuc_", strain_ref1, sep="")]] = index_nuc_strain_ref1

									# tmp_nuc[[paste("corr1_", strain_ref1, sep="")]] = signif(nuc_strain_ref1$nucs[[1]]$corr,5)

									# tmp_nuc[[paste("corr2_", strain_ref1, sep="")]] = signif(nuc_strain_ref1$nucs[[2]]$corr,5)

									# tmp_nuc[[paste("corr3_", strain_ref1, sep="")]] = signif(nuc_strain_ref1$nucs[[3]]$corr,5)

									# strain_ref2

									tmp_nuc[[paste("chr_", strain_ref2, sep="")]] = roi_strain_ref2$chr

									tmp_nuc[[paste("lower_bound_", strain_ref2, sep="")]] = nuc_strain_ref2$lower_bound

									tmp_nuc[[paste("upper_bound_", strain_ref2, sep="")]] = nuc_strain_ref2$upper_bound

									tmp_nuc[[paste("means_", strain_ref2, sep="")]] = signif(mean(reads_strain_ref2),5)

									tmp_nuc[[paste("sd_", strain_ref2, sep="")]] = signif(sd(reads_strain_ref2),5)

									tmp_nuc[[paste("nb_reads_", strain_ref2, sep="")]] = length(reads_strain_ref2)

									tmp_nuc[[paste("index_nuc_", strain_ref2, sep="")]] = index_nuc_strain_ref2

									# tmp_nuc[[paste("corr1_", strain_ref2, sep="")]] = signif(nuc_strain_ref2$nucs[[1]]$corr,5)

									# tmp_nuc[[paste("corr2_", strain_ref2, sep="")]] = signif(nuc_strain_ref2$nucs[[2]]$corr,5)

									# tmp_nuc[[paste("corr3_", strain_ref2, sep="")]] = signif(nuc_strain_ref2$nucs[[3]]$corr,5)

									# common

									tmp_nuc[["lod_score"]] = signif(lod_score,5)

									# print(tmp_nuc)

									common_nuc = dfadd(common_nuc, tmp_nuc)

								}

							}

						}

					}				

				} else {

		      print("WARNING! No roi for strain ref 2.")

			  }

		  }

		}

		if(length(unique(common_nuc[,1:3])[,1]) != length((common_nuc[,1:3])[,1])) {

			index_redundant = which(apply(common_nuc[,1:3][-length(common_nuc[,1]),] ==  common_nuc[,1:3][-1,] ,1,sum) == 3)

			to_remove_list = c()			

			for (i in 1:length(index_redundant)) {	

				if (common_nuc[index_redundant[i],15] < common_nuc[index_redundant[i]+1,15]) {

				  to_remove = index_redundant[i]

				}	 else {

					to_remove = index_redundant[i] + 1

			  } 		

				to_remove_list = c(to_remove_list, to_remove)

			}

			common_nuc = common_nuc[-to_remove_list,]

		}

		if(length(unique(common_nuc[,8:10])[,1]) != length((common_nuc[,8:10])[,1])) {

			index_redundant = which(apply(common_nuc[,8:10][-length(common_nuc[,1]),] == common_nuc[,8:10][-1,] ,1,sum) == 3)

			to_remove_list = c()			

			for (i in 1:length(index_redundant)) {	

				if (common_nuc[index_redundant[i],15] < common_nuc[index_redundant[i]+1,15]) {

				  to_remove = index_redundant[i]

				}	 else {

					to_remove = index_redundant[i] + 1

			  } 		

				to_remove_list = c(to_remove_list, to_remove)

			}

			common_nuc = common_nuc[-to_remove_list,]

		}

		return(list(common_nuc, lod_scores))

	} else {

		print("WARNING, no nucs for strain_ref1.")

		return(NULL)

	}

### Returns a list of clusterized nucleosomes, and all computed lod scores.

}, ex=function(){	

    # Define new translate_roi function...

    translate_roi = function(roi, strain2, big_roi=NULL, config=NULL) {

      return(roi)

    }

    # Binding it by uncomment follwing lines.

    unlockBinding("translate_roi", as.environment("package:nucleominer"))

    unlockBinding("translate_roi", getNamespace("nucleominer"))

    assign("translate_roi", translate_roi, "package:nucleominer")

    assign("translate_roi", translate_roi, getNamespace("nucleominer"))

    lockBinding("translate_roi", getNamespace("nucleominer"))

    lockBinding("translate_roi", as.environment("package:nucleominer"))  

	# Dealing with a region of interest

	roi =list(name="example", begin=1000,  end=1300, chr="1", genome=rep("A",301), strain_ref1 = "STRAINREF1")

	roi2 = translate_roi(roi, roi$strain_ref1)

	replicates = list()

	for (j in 1:2) {

		samples = list()

		for (i in 1:3) {

			# Create TF output

			tf_nuc = list("chr"=paste("chr", roi$chr, sep=""), "center"=(roi$end + roi$begin)/2, "width"= 150, "correlation.score"= 0.9)

			outputs = dfadd(NULL,tf_nuc)

			outputs = filter_tf_outputs(outputs, roi$chr, roi$begin, roi$end)

			# Generate corresponding reads

			nb_reads = round(runif(1,170,230))

			reads = round(rnorm(nb_reads, tf_nuc$center,20))

			u_reads = sort(unique(reads))

			strands = sample(c(rep("R",ceiling(length(u_reads)/2)),rep("F",floor(length(u_reads)/2))))

			counts = apply(t(u_reads), 2, function(r) { sum(reads == r)})

			shifts = apply(t(strands), 2, function(s) { if (s == "F") return(-tf_nuc$width/2) else return(tf_nuc$width/2)})

			u_reads = u_reads + shifts

			inputs = data.frame(list("V1" = rep(roi$chr, length(u_reads)), 

			                         "V2" = u_reads, 

															 "V3" = strands, 

															 "V4" = counts), stringsAsFactors=FALSE)

			samples[[length(samples) + 1]] = list(id=1, marker="Mnase_Seq", strain=paste("strain_ex",j,sep=""), total_reads = 10000000, roi=roi, inputs=inputs, outputs=outputs)

		}

		replicates[[length(replicates) + 1]] = samples

	}

	print(align_inter_strain_nucs(replicates))

})

fetch_mnase_replicates = function(# Prefetch data

### Fetch and filter inputs and outpouts per region of interest. Organize it per replicates.

strain, ##<< The strain we want mnase replicatesList of replicates. Each replicates is a vector of sample ids.

roi, ##<< Region of interest.

all_samples, ##<< Global list of samples.

config=NULL, ##<< GLOBAL config variable

only_fetch=FALSE, ##<< If TRUE, only fetch and not filtering. It is used tio load sample files into memory before forking.

get_genome=FALSE, ##<< If TRUE, load corresponding genome sequence.

get_ouputs=TRUE##<< If TRUE, get also ouput corresponding TF output files.

) {

	samples=list()

  samples_ids = unique(all_samples[all_samples$marker == "Mnase_Seq" & all_samples$strain == strain,]$id)

	for (i in samples_ids) {

		sample = as.list(all_samples[all_samples$id==i,])

    sample$orig_roi = roi

    sample$roi = translate_roi(roi, sample$strain, config = config)

		if (get_genome) {

			# Get Genome

      sample$roi$genome = get_content(config$FASTA_REFERENCE_GENOME_FILES[[sample$strain]], "fasta")[[switch_pairlist(config$FASTA_INDEXES[[sample$strain]])[[sample$roi$chr]]]][sample$roi$begin:sample$roi$end]		

		}

		# Get inputs

		sample$inputs = get_content(paste(config$ALIGN_DIR, "/TF/sample_", i, "_TF.txt", sep=""), "table", stringsAsFactors=FALSE) 

		sample$total_reads = sum(sample$inputs[,4]) 	

		if (!only_fetch) {

		  sample$inputs = filter_tf_inputs(sample$inputs, sample$roi$chr, min(sample$roi$begin, sample$roi$end), max(sample$roi$begin, sample$roi$end), 300)

	  }

	  # Get TF outputs for Mnase_Seq samples

		if (sample$marker == "Mnase_Seq" & get_ouputs) {	

			sample$outputs = get_content(paste(config$ALIGN_DIR, "/TF/sample_", i, "_all_nucs.tab", sep=""), "table", header=TRUE, sep="\t")

			if (!only_fetch) {

	  		sample$outputs = filter_tf_outputs(sample$outputs, sample$roi$chr,  min(sample$roi$begin, sample$roi$end), max(sample$roi$begin, sample$roi$end), 300)

  		}

		}

		samples[[length(samples) + 1]] = sample		

	}

  return(samples)

}

substract_region = function(# Substract to a list of regions an other list of regions that intersect it.

### This fucntion embed a recursive part. It occurs when a substracted region split an original region on two.

region1, ##<< Original regions. 

region2 ##<< Regions to substract.

) {

  rec_substract_region = function(region1, region2) {

  non_inter_fuzzy = apply(t(1:length(region1[,1])), 2, function(i) {

    cur_fuzzy = region1[i,]

    inter_wp = region2[region2$lower_bound <= cur_fuzzy$upper_bound & region2$upper_bound >= cur_fuzzy$lower_bound,]

    if (length(inter_wp[,1]) > 0) {

      ret = c()

      for (j in 1:length(inter_wp[,1])) {

        cur_wp = inter_wp[j,]

        if (cur_wp$lower_bound <= cur_fuzzy$lower_bound & cur_fuzzy$upper_bound <= cur_wp$upper_bound) {

          # remove cur_fuzzy

          ret = c()

          break

        } else if (cur_wp$lower_bound <= cur_fuzzy$lower_bound & cur_wp$upper_bound < cur_fuzzy$upper_bound) {

          # crop fuzzy 

          cur_fuzzy$lower_bound = cur_wp$upper_bound + 1 

          ret = cur_fuzzy

        } else if (cur_fuzzy$lower_bound < cur_wp$lower_bound & cur_fuzzy$upper_bound <= cur_wp$upper_bound) {

          # crop fuzzy 

          cur_fuzzy$upper_bound = cur_wp$lower_bound - 1 

          ret = cur_fuzzy

        } else if (cur_wp$lower_bound > cur_fuzzy$lower_bound & cur_wp$upper_bound < cur_fuzzy$upper_bound) {

          # split fuzzy

          tmp_ret_fuzzy_1 = cur_fuzzy

          tmp_ret_fuzzy_1$upper_bound = cur_wp$lower_bound - 1 

          tmp_ret_fuzzy_2 = cur_fuzzy

          tmp_ret_fuzzy_2$lower_bound = cur_wp$upper_bound + 1 

          ret = rec_substract_region(rbind(tmp_ret_fuzzy_1, tmp_ret_fuzzy_2), inter_wp)

          # print(ret)

          # ret = cur_fuzzy

          break

        } else {

          stop("WARNING NO ADAPTED CASE!")

        }          

      }

      return(ret)

    } else {

      return(cur_fuzzy)

    }

  })

  }

  non_inter_fuzzy = rec_substract_region(region1, region2)

  if (is.null(non_inter_fuzzy)) {return(non_inter_fuzzy)}

  tmp_ulist = unlist(non_inter_fuzzy)

  tmp_names = names(tmp_ulist)[1:4]

  non_inter_fuzzy = data.frame(matrix(tmp_ulist, ncol=4, byrow=TRUE), stringsAsFactors=FALSE) 

  names(non_inter_fuzzy) = tmp_names

  non_inter_fuzzy$chr = as.character(non_inter_fuzzy$chr)     

  non_inter_fuzzy$chr = as.numeric(non_inter_fuzzy$chr)     

  non_inter_fuzzy$lower_bound = as.numeric(non_inter_fuzzy$lower_bound)

  non_inter_fuzzy$upper_bound = as.numeric(non_inter_fuzzy$upper_bound)

  non_inter_fuzzy = non_inter_fuzzy[order(non_inter_fuzzy$lower_bound),]      

  return(non_inter_fuzzy)

}

union_regions = function(# Aggregate regions that intersect themnselves.

### This function is based on sort of lower bounds to detect regions that intersect. We compare lower bound and upper bound of the porevious item. This function embed a while loop and break break regions list become stable.  

regions ##<< The Regions to be aggregated

) {

  old_length = length(regions[,1])

  new_length = 0

  while (old_length != new_length) {

    regions = regions[order(regions$lower_bound), ]

    regions$stop = !c(regions$lower_bound[-1] - regions$upper_bound[-length(regions$lower_bound)] <= 0, TRUE)      

    vec_end_1 = which(regions$stop)

    if (length(vec_end_1) == 0) {

      vec_end_1 = c(length(regions$stop))      

    } 

    if (vec_end_1[length(vec_end_1)] != length(regions$stop)) {

      vec_end_1 = c(vec_end_1, length(regions$stop))

    }

    vec_beg_1 = c(1, vec_end_1[-length(vec_end_1)] + 1)

    union = apply(t(1:length(vec_beg_1)), 2, function(i) {

      chr = regions$chr[vec_beg_1[i]]

      lower_bound = min(regions$lower_bound[vec_beg_1[i]:vec_end_1[i]])

      upper_bound = max(regions$upper_bound[vec_beg_1[i]:vec_end_1[i]])

      roi_index = regions$roi_index[vec_beg_1[i]] 

      data.frame(list(chr=chr, lower_bound=lower_bound, upper_bound=upper_bound, roi_index=roi_index))

      })

    union = collapse_regions(union)

    old_length = length(regions[,1])

    new_length = length(union[,1])

    regions = union

  }

  return(union)

}

remove_aligned_wp = function(# Remove wp nucs from common nucs list.

### It is based on common wp nucs index on nucs and region.  

strain_maps, ##<< Nuc maps.

roi_index, ##<< The region of interest index.

tmp_common_nucs, ##<< the list of wp nucs.

strain##<< The strain to consider.

){

  fuzzy_nucs = strain_maps[[strain]]

  fuzzy_nucs = fuzzy_nucs[fuzzy_nucs$roi_index == roi_index,]

  fuzzy_nucs = fuzzy_nucs[order(fuzzy_nucs$index_nuc),]

  if (length(fuzzy_nucs[,1]) == 0) {return(fuzzy_nucs)}

  if (sum(fuzzy_nucs$index_nuc == min(fuzzy_nucs$index_nuc):max(fuzzy_nucs$index_nuc)) != max(fuzzy_nucs$index_nuc)) {"Warning in index!"}

  anti_index_1 = tmp_common_nucs[[paste("index_nuc", strain, sep="_")]]

  fuzzy_nucs = fuzzy_nucs[-anti_index_1,]

  return(fuzzy_nucs)

}

translate_regions = function(# Translate a list of regions from a strain ref to another.

### This function is an eloborated call to translate_roi.

regions, ##<< Regions to be translated. 

combi, ##<< Combination of strains.

roi_index, ##<< The region of interest index.

config=NULL, ##<< GLOBAL config variable

roi ##<< The region of interest.

) {

  tr_regions = apply(t(1:length(regions[,1])), 2, function(i) {

    tmp_regions_ref2 = list(name="foo", begin=regions[i,]$lower_bound, end=regions[i,]$upper_bound, chr=as.character(regions[i,]$chr), strain_ref = combi[2])

    big_roi =  roi

    trs_tmp_regions_ref2 = translate_roi(tmp_regions_ref2, combi[1], config = config, big_roi = big_roi) 

    data.frame(list(chr=trs_tmp_regions_ref2$chr, lower_bound=min(trs_tmp_regions_ref2$begin, trs_tmp_regions_ref2$end), upper_bound=max(trs_tmp_regions_ref2$begin, trs_tmp_regions_ref2$end), roi_index=roi_index))

    })

  return(collapse_regions(tr_regions))

}

collapse_regions = function(# reformat an "apply  manipulated" list of regions

### Utils to reformat an "apply  manipulated" list of regions

regions ##< a list of regions

) {

  regions = do.call(rbind, regions)      

  regions$chr = as.character(regions$chr)     

  regions$chr = as.numeric(regions$chr)     

  regions$lower_bound = as.numeric(regions$lower_bound)

  regions$upper_bound = as.numeric(regions$upper_bound)

  regions = regions[order(regions$lower_bound),]      

  return(regions)  

}

extract_wp = function(# Extract wp nucs from nuc map. 

### Function based on common wp nuc index and roi_index.

strain_maps, ##<< Nuc maps.

roi_index, ##<< The region of interest index.

strain, ##<< The strain to consider.

tmp_common_nucs ##<< the list of wp nucs.

) {

  wp_nucs = apply(t(tmp_common_nucs[[paste("index_nuc", strain, sep="_")]]), 2, function(i) {

    tmp_wp_nucs = strain_maps[[strain]]

    tmp_wp_nucs = tmp_wp_nucs[tmp_wp_nucs$roi_index == roi_index & tmp_wp_nucs$index_nuc == i,]

    return(tmp_wp_nucs)

    })

  return(collapse_regions(wp_nucs))

}

crop_fuzzy = function(# Crop bound of regions according to region of interest bound

### The fucntion is no more necessary since we remove "big_roi" bug in translate_roi function.

tmp_fuzzy_nucs, ##<< the regiuons to be croped.

roi, ##<< The region of interest.

strain, ##<< The strain to consider.

config=NULL ##<< GLOBAL config variable

) {

  tr_roi = translate_roi(roi, strain, config = config)

  tr_roi_begin = min(tr_roi$begin, tr_roi$end)

  tr_roi_end = max(tr_roi$begin, tr_roi$end)

  if (length(tmp_fuzzy_nucs[tmp_fuzzy_nucs$lower_bound < tr_roi_begin,1]) > 0) {

    tmp_fuzzy_nucs[tmp_fuzzy_nucs$lower_bound < tr_roi_begin,]$lower_bound = tr_roi_begin

  }

  if (length(tmp_fuzzy_nucs[tmp_fuzzy_nucs$upper_bound < tr_roi_begin,1]) > 0) {

    tmp_fuzzy_nucs[tmp_fuzzy_nucs$upper_bound < tr_roi_begin,]$upper_bound = tr_roi_begin

  }

  if (length(tmp_fuzzy_nucs[tmp_fuzzy_nucs$lower_bound > tr_roi_end,1]) > 0) {

    tmp_fuzzy_nucs[tmp_fuzzy_nucs$lower_bound > tr_roi_end,]$lower_bound = tr_roi_end

  }

  if (length(tmp_fuzzy_nucs[tmp_fuzzy_nucs$upper_bound > tr_roi_end,1]) > 0) {

    tmp_fuzzy_nucs[tmp_fuzzy_nucs$upper_bound > tr_roi_end,]$upper_bound = tr_roi_end

  }

  tmp_fuzzy_nucs = tmp_fuzzy_nucs[tmp_fuzzy_nucs$upper_bound != tmp_fuzzy_nucs$lower_bound,]

  return(tmp_fuzzy_nucs)

}

get_fuzzy = function(# Compute the fuzzy nucs.

### This function aggregate non common wp nucs for each strain and substract common wp nucs. It does not take care about the size of the resulting fuzzy regions. It will be take into account in the count read part og the pipeline.

combi, ##<< The strain combination to consider.

roi, ##<< The region of interest.

roi_index, ##<< The region of interest index.

strain_maps, ##<< Nuc maps.

common_nuc_results, ##<< Common wp nuc maps

config=NULL ##<< GLOBAL config variable

) {

  print(roi_index)

  PLOT = FALSE

  tmp_common_nucs = common_nuc_results[[paste(combi[1], combi[2], sep="_")]]

  tmp_common_nucs = tmp_common_nucs[tmp_common_nucs$roi_index == roi_index, ]

  print(paste("Dealing with fuzzy from", combi[1]))

  tmp_fuzzy_nucs_1 = remove_aligned_wp(strain_maps, roi_index, tmp_common_nucs, combi[1])

  tmp_fuzzy_nucs_1 = crop_fuzzy(tmp_fuzzy_nucs_1, roi, combi[1], config)

  if (length(tmp_fuzzy_nucs_1[,1]) == 0) {return(NULL)}

  agg_fuzzy_1 = union_regions(tmp_fuzzy_nucs_1)    

  if (PLOT) for (i in 1:length(agg_fuzzy_1[,1])) {

    lines(c(agg_fuzzy_1[i,]$lower_bound, agg_fuzzy_1[i,]$upper_bound), c(+3.1,+3.1), col=2)

  }

  print(paste("Dealing with fuzzy from ", combi[2]))

  tmp_fuzzy_nucs_2 = remove_aligned_wp(strain_maps, roi_index, tmp_common_nucs, combi[2])

  if (length(tmp_fuzzy_nucs_2[,1]) == 0) {return(NULL)}

  agg_fuzzy_2 = union_regions(tmp_fuzzy_nucs_2)      

  agg_fuzzy_2 = crop_fuzzy(agg_fuzzy_2, roi, combi[2], config)

  tr_agg_fuzzy_2 = translate_regions(agg_fuzzy_2, combi, roi_index, roi=roi, config=config)

  tr_agg_fuzzy_2 = crop_fuzzy(tr_agg_fuzzy_2, roi, combi[2], config)

  # tr_agg_fuzzy_2 = union_regions(tr_agg_fuzzy_2)

  if (PLOT) for (i in 1:length(tr_agg_fuzzy_2[,1])) {

    lines(c(tr_agg_fuzzy_2[i,]$lower_bound, tr_agg_fuzzy_2[i,]$upper_bound), c(+3.3,+3.3), col=2)

  }

  print("Dealing with fuzzy from both...")

  all_fuzzy = union_regions(rbind(agg_fuzzy_1, tr_agg_fuzzy_2))

  if (PLOT) for (i in 1:length(all_fuzzy[,1])) {

    lines(c(all_fuzzy[i,]$lower_bound, all_fuzzy[i,]$upper_bound), c(+3.2, +3.2), col=1)

  }

  print(paste("Dealing with wp from", combi[1]))

  wp_nucs_1 = extract_wp(strain_maps, roi_index, combi[1], tmp_common_nucs)

  if (PLOT) for (i in 1:length(wp_nucs_1[,1])) {

    lines(c(wp_nucs_1[i,]$lower_bound, wp_nucs_1[i,]$upper_bound), c(+3.5,+3.5), col=3)

  }

  print(paste("Dealing with wp from", combi[2]))

  wp_nucs_2 = extract_wp(strain_maps, roi_index, combi[2], tmp_common_nucs)

  tr_wp_nucs_2 = translate_regions(wp_nucs_2, combi, roi_index, roi=roi, config=config)

  if (PLOT) for (i in 1:length(tr_wp_nucs_2[,1])) {

    lines(c(tr_wp_nucs_2[i,]$lower_bound, tr_wp_nucs_2[i,]$upper_bound), c(+3.7,+3.7), col=3)

  }

  print("Dealing with wp from both...")

  all_wp = union_regions(rbind(wp_nucs_1[,1:4], tr_wp_nucs_2))

  if (PLOT) for (i in 1:length(all_wp[,1])) {

    lines(c(all_wp[i,]$lower_bound, all_wp[i,]$upper_bound), c(+3.6, +3.6), col=1)

  }  

  print("Dealing with fuzzy and wp...")

  non_inter_fuzzy = substract_region(all_fuzzy, all_wp)

  if (is.null(non_inter_fuzzy)) { return(NULL) }

  non_inter_fuzzy$len = non_inter_fuzzy$upper_bound - non_inter_fuzzy$lower_bound

  # non_inter_fuzzy = non_inter_fuzzy[non_inter_fuzzy$len >= min_fuzz_width,]

  if (PLOT) for (i in 1:length(non_inter_fuzzy[,1])) {

    lines(c(non_inter_fuzzy[i,]$lower_bound, non_inter_fuzzy[i,]$upper_bound), c(+3.9, +3.9), col=1)

  }

  non_inter_fuzzy$index_nuc = 1:length(non_inter_fuzzy[,1])

  return (non_inter_fuzzy)

}

get_all_reads = function(# Retrieve Reads 

### Retrieve reads for a given marker, combi, form.

marker, ##<< The marker to considere.

combi, ##<< The starin combination to considere.

form="wp", ##<< The nuc form to considere.

config=NULL ##<< GLOBAL config variable

) {

	all_reads = NULL

  for (manip in c("Mnase_Seq", marker)) {

    if (form == "fuzzy") {

		  out_filename = paste(config$RESULTS_DIR, "/",combi[1],"_",combi[2],"_",manip,"_fuzzy_and_nbreads.tab",sep="")

  		tmp_res = read.table(file=out_filename, header=TRUE)				

			tmp_res = tmp_res[tmp_res[,3] - tmp_res[,2] > 75,]

      tmp_res$form = form

    } else if (form == "wp") {

		 	out_filename = paste(config$RESULTS_DIR, "/",combi[1],"_",combi[2],"_",manip,"_wp_and_nbreads.tab",sep="")

  		tmp_res = read.table(file=out_filename, header=TRUE)				

      tmp_res$form = form

    } else if (form == "wpfuzzy") {

		 	out_filename = paste(config$RESULTS_DIR, "/",combi[1],"_",combi[2],"_",manip,"_wp_and_nbreads.tab",sep="")

  		tmp_res = read.table(file=out_filename, header=TRUE)				

      tmp_res$form = "wp"

		  out_filename = paste(config$RESULTS_DIR, "/",combi[1],"_",combi[2],"_",manip,"_fuzzy_and_nbreads.tab",sep="")

  		tmp_res2 = read.table(file=out_filename, header=TRUE)				

			tmp_res2 = tmp_res2[tmp_res2[,3] - tmp_res2[,2] > 75,]

      tmp_res2$form = "fuzzy"

      tmp_res = rbind(tmp_res, tmp_res2)

    }

		if (is.null(all_reads)) {

			all_reads = tmp_res[,c(1:9,length(tmp_res))]

		}     

		tmp_res = tmp_res[,-c(1:9,length(tmp_res))]

		all_reads = cbind(all_reads, tmp_res)

  }

  return(all_reads)	

}

get_design = function(# Build the design for deseq

### This function build the design according sample properties.

marker, ##<< The marker to considere.

combi, ##<< The starin combination to considere.

all_samples ##<< Global list of samples.

) {

  off1 = 0

  off2 = 0

	manips = c("Mnase_Seq", marker)

	design_rownames = c()

	design_manip = c()

	design_strain = c()

  off2index = function(off) {	

  	switch(toString(off), 

  		"1"=c(0,1,1),

  	  "2"=c(1,0,1),

    	"3"=c(1,1,0), 

  		c(1,1,1)

  		)	

  }

	for (manip in manips) {

		tmp_samples = all_samples[ ((all_samples$strain == combi[1] | all_samples$strain == combi[2]) &  all_samples$marker == manip), ]

		tmp_samples = tmp_samples[order(tmp_samples$strain), ]

		if (manip == "H3K4me1" & (off1 != 0 & off2 ==0 )) {

			tmp_samples = tmp_samples[c(off2index(off1), c(1,1)) == 1,]

		} else {

			if (manip != "Mnase_Seq" & (off1 != 0 | off2 !=0)) {

				tmp_samples = tmp_samples[c(off2index(off1), off2index(off2)) == 1,]

			}

		}

		design_manip = c(design_manip, rep(manip, length(tmp_samples$id)))				

		for (strain in combi) {			

			cols = apply(t(tmp_samples[ (tmp_samples$strain == strain &  tmp_samples$marker == manip), ]$id), 2, function(i){paste(strain, manip, i, sep="_")})

			design_strain = c(design_strain, rep(strain, length(cols)))

			design_rownames = c(design_rownames, cols)

		}

	}

	snep_design = data.frame( row.names=design_rownames, manip=design_manip, strain=design_strain)

	return(snep_design)	

}

plot_dist_samples = function(# Plot the distribution of reads.

### This fuxntion use the deseq nomalization feature to compare qualitatively the distribution.

strain, ##<< The strain to considere.

marker, ##<< The marker to considere.

res, ##<< Data

all_samples, ##<< Global list of samples.

NEWPLOT = TRUE ##<< If FALSE the curve will be add to the current plot.

) {			

	cols = apply(t(all_samples[ (all_samples$strain == strain &  all_samples$marker == marker), ]$id), 2, function(i){paste(strain, marker, i, sep="_")})

	snepCountTable = res[,cols]

	snepDesign = data.frame(

		row.names = cols,

		manip = rep(marker, length(cols)),

		strain = rep(strain, length(cols))

		)

	cdsFull = newCountDataSet(snepCountTable, snepDesign)

	sizeFactors = estimateSizeFactors(cdsFull)

	# print(sizeFactors[[1]])

	sample_ids = all_samples[ (all_samples$strain == strain &  all_samples$marker == marker), ]$id

	if (NEWPLOT) {

		plot(density(res[,paste(strain, marker, sample_ids[1], sep="_")] / sizeFactors[[1]][1]), col=0, main=paste(strain, marker))				

		NEWPLOT = FALSE

	}

	for (it in 1:length(sample_ids)) {

		sample_id = sample_ids[it]

		lines(density(res[,paste(strain, marker, sample_id, sep="_")] / sizeFactors[[1]][it]), col = it + 1, lty = it)

	}

  legend("topright", col=(1:length(sample_ids))+1, lty=1:length(sample_ids), legend=cols)

}

analyse_design = function(# Launch deseq methods.

### This function is based on deseq example. It mormalizes data, fit data to GLM model with and without interaction term and compare the two l;=models.

snep_design, ##<< The design to considere.

reads ##<< The data to considere. 

) {

	snep_count_table = reads[, rownames(snep_design)]

	cdsFull = newCountDataSet(snep_count_table, snep_design)

	cdsFull1 = estimateDispersions(estimateSizeFactors(cdsFull), fitType="local", method="pooled", sharingMode="maximum")

	fit1 = fitNbinomGLMs(cdsFull1, count ~ manip * strain)

	cdsFull0 = estimateDispersions(estimateSizeFactors(cdsFull), fitType="local", method="pooled", sharingMode="maximum")

	fit0 = fitNbinomGLMs(cdsFull0, count ~ manip + strain)

	pvalsGLM = nbinomGLMTest( fit1, fit0 )

	return(list(fit1, fit0, snep_design, pvalsGLM))	

}

get_sneps = structure(function(# Compute the list of SNEPs for a given set of marker, strain combination and nuc form.

### This function uses 

marker, ##<< The marker involved.

combi, ##<< The strain combination involved.

form, ##<< the nuc form involved. 

all_samples, ##<< Global list of samples.

config=NULL ##<< GLOBAL config variable

) {

  # PRETREAT

  d = get_design(marker, combi, all_samples)

  reads = get_all_reads(marker, combi, form, config=config)

  # RUN ANALYSE

  tmp_analyse = analyse_design(d, reads)

  # RESULTS

	fit1 = tmp_analyse[[1]]

	fit0 = tmp_analyse[[2]]

  k = names(fit1)

  reads[[k[2]]] = signif(fit1[[k[2]]], 5)

  reads[[k[3]]] = signif(fit1[[k[3]]], 5)

  reads[[k[4]]] = signif(fit1[[k[4]]], 5)

	reads$pvalsGLM = signif(tmp_analyse[[4]], 5)

	snep_design = tmp_analyse[[3]]

  print(snep_design)

	fdr = 0.0001

	thres = FDR(reads$pvalsGLM, fdr) 

	reads$snep_index = reads$pvalsGLM < thres

	print(paste(sum(reads$snep_index), " SNEPs found for ", length(reads[,1])," nucs and ", fdr*100,"% of FDR.", sep = ""))

  return(reads)

  },  ex=function(){	

    marker = "H3K4me1"

    combi = c("BY", "YJM") 

    form = "wpfuzzy" # "wp" | "fuzzy" | "wpfuzzy"

    # foo = get_sneps(marker, combi, form)

    # foo = get_sneps("H4K12ac", c("BY", "RM"), "wp")

})

ROM2ARAB = function(# Roman to Arabic pair list.

### Util to convert Roman to Arabic

){list(

  "I" = 1,

  "II" = 2,

  "III" = 3,

  "IV" = 4,

  "V" = 5,

  "VI" = 6,

  "VII" = 7,

  "VIII" = 8,

  "IX" = 9,

  "X" = 10,

  "XI" = 11,

  "XII" = 12,

  "XIII" = 13,

  "XIV" = 14,

  "XV" = 15,

  "XVI" = 16,

  "XVII" = 17,

  "XVIII" = 18,

  "XIX" = 19,

  "XX" = 20

)}

switch_pairlist = structure(function(# Switch a pairlist

### Take a pairlist key:value and return the switched pairlist value:key.

l ##<< The pairlist to switch. 

) {

	ret = list()

	for (name in names(l)) {

		ret[[as.character(l[[name]])]] = name

	}

	ret

### The switched pairlist.

}, ex=function(){

	l = list(key1 = "value1", key2 = "value2")

	print(switch_pairlist(l))

})

ARAB2ROM = function(# Arabic to Roman pair list.

### Util to convert Arabicto Roman

){switch_pairlist(ROM2ARAB())}

translate_roi = structure(function(# Translate coords of a genome region.

### This function is used in the examples, usualy you have to define your own translation function and overwrite this one using \emph{unlockBinding} features. Please, refer to the example.

roi, ##<< Original genome region of interest.

strain2, ##<< The strain in wich you want the genome region of interest.

config=NULL, ##<< GLOBAL config variable

big_roi=NULL ##<< A largest region than roi use to filter c2c if it is needed.

) {

	strain1 = roi$strain_ref

  reverse = (strain1=="RM" & strain2=="BY") | strain1=="YJM"

	if (strain1 == strain2) {

		roi$length = roi$end - roi$begin + sign(roi$end - roi$begin) * 1	

		return(roi)

	}

	# Launch c2c file

	if (reverse) {

		c2c_file = list(filename=config$C2C_FILES[[paste(strain2, "-", strain1, sep="")]])

	} else {

		c2c_file = list(filename=config$C2C_FILES[[paste(strain1, "-", strain2, sep="")]])

	}

	c2c = get_content(c2c_file$filename, "table", stringsAsFactors=FALSE)	

	# filtering it

  c2c = c2c[c2c$V6=="-",]

	# Reverse

	if (reverse) {

		tmp_col = c2c$V1

		c2c$V1 = c2c$V7

		c2c$V7 = tmp_col

		tmp_col = c2c$V2

		c2c$V2 = c2c$V9

		c2c$V9 = tmp_col

		tmp_col = c2c$V3

		c2c$V3 = c2c$V10

		c2c$V10 = tmp_col

	}

	# Restrict c2c to big_roi

	# if (FALSE) {

	if (!is.null(big_roi)) {

		if (roi$strain_ref != big_roi$strain_ref) {

      # print("WARNING, big_roi and roi not in the same strain_ref... translating big_roi.")

      big_roi = translate_roi(big_roi, roi$strain_ref, config=config)

    }

    if (big_roi$end < big_roi$begin) {

      tmp_var = big_roi$begin

      big_roi$begin = big_roi$end 

      big_roi$end = tmp_var        

      big_roi$length = big_roi$length        

    }    

    if (big_roi$chr!=roi$chr | roi$end > big_roi$end | roi$end < big_roi$begin | roi$begin > big_roi$end | roi$begin < big_roi$begin) {

      print("WARNING! Trying to translate a roi not included in a big_roi.")

      return(NULL)

    }

		if (strain1 == "BY") {

			big_chro_1 = paste("chr", ARAB2ROM()[[big_roi$chr]], sep="")

		} else if (strain1 == "RM") {			

		  big_chro_1 = paste("supercontig_1.",big_roi$chr,sep="")

		} else if (strain1 == "YJM") {			

		  big_chro_1 = switch_pairlist(config$FASTA_INDEXES$YJM)[[big_roi$chr]]

		}

		big_begin_1 = big_roi$begin

	  big_end_1 = big_roi$end

		c2c = c2c[c2c$V1==big_chro_1,]

    if (length(c2c[c2c$V3 < big_begin_1 & c2c$V2 < c2c$V3, 1] > 0)) {c2c[c2c$V3 < big_begin_1 & c2c$V2 < c2c$V3,c("V2", "V3") ] = big_begin_1}

    if (length(c2c[c2c$V2 > big_end_1   & c2c$V2 < c2c$V3, 1] > 0)) {c2c[c2c$V2 > big_end_1   & c2c$V2 < c2c$V3, c("V2", "V3")] = big_end_1}

    if (length(c2c[c2c$V2 < big_begin_1 & c2c$V3 < c2c$V2, 1] > 0)) {c2c[c2c$V2 < big_begin_1 & c2c$V3 < c2c$V2,c("V2", "V3") ] = big_begin_1}

    if (length(c2c[c2c$V3 > big_end_1   & c2c$V3 < c2c$V2, 1] > 0)) {c2c[c2c$V3 > big_end_1   & c2c$V3 < c2c$V2, c("V2", "V3")] = big_end_1}

    c2c = c2c[c2c$V2 - c2c$V3 != 0,]

	}

  #	Convert initial roi$chr into c2c format

	if (strain1 == "BY") {

		chro_1 = paste("chr", ARAB2ROM()[[roi$chr]], sep="")

	} else if (strain1 == "RM") {			

	  chro_1 = paste("supercontig_1.",roi$chr,sep="")

	} else if (strain1 == "YJM") {			

	  chro_1 = switch_pairlist(config$FASTA_INDEXES$YJM)[[roi$chr]]

	}

	begin_1 = roi$begin

  end_1 = roi$end

  # Computing equivalent strain_2 alignment coordinates	

	if (reverse) {

  	tmptransfostart = c2c[c2c$V1==chro_1 & ((c2c$V3>=begin_1 & c2c$V2<=begin_1 & c2c$V8==1) | (c2c$V2>=begin_1 & c2c$V3<=begin_1 & c2c$V8==-1)),]

    tmptransfostop = c2c[c2c$V1==chro_1 &  ((c2c$V3>=end_1   & c2c$V2<=end_1   & c2c$V8==1) | (c2c$V2>=end_1   & c2c$V3<=end_1   & c2c$V8==-1)),]

	} else {

		tmptransfostart = c2c[c2c$V1==chro_1 & c2c$V3>=begin_1 & c2c$V2<=begin_1,]

	  tmptransfostop = c2c[c2c$V1==chro_1 & c2c$V3>=end_1 & c2c$V2<=end_1,]

	}

	# Never happend conditions ...	

	{

		if (length(tmptransfostart$V8) == 0) {

			# begin_1 is between to lines: shift begin_1 to the start of 2nd line.

      if (sum(c2c$V3 >= c2c$V2) != 0) {

  			tmp_c2c = c2c[c2c$V1==chro_1 & c2c$V2>=begin_1,]

  			begin_1 = min(tmp_c2c$V2)

      } else {

  			tmp_c2c = c2c[c2c$V1==chro_1 & c2c$V3>=begin_1,]

  			begin_1 = min(tmp_c2c$V3)

      }

			if (reverse) {

		  	tmptransfostart = c2c[c2c$V1==chro_1 & ((c2c$V3>=begin_1 & c2c$V2<=begin_1 & c2c$V8==1) | (c2c$V2>=begin_1 & c2c$V3<=begin_1 & c2c$V8==-1)),]

			} else {

				tmptransfostart = c2c[c2c$V1==chro_1 & c2c$V3>=begin_1 & c2c$V2<=begin_1,]

			}

			if (length(tmptransfostart$V8) == 0) {

				if (!is.null(big_roi)) {

					return(NULL)

					tmptransfostart = c2c[c2c$V1==chro_1 & c2c$V3>=big_roi$begin & c2c$V2<=big_roi$begin,]

				} else { 

					# return(NULL)

					# print(c2c[c2c$V1==chro_1 & c2c$V2<=end_1 & c2c$V3>=begin_1,])

					# print(c2c[c2c$V1==chro_1,])

					print(tmptransfostart)

					print(tmptransfostop)

					stop("Never happend condition 1.")					

				}

			}

		} 

		if (length(tmptransfostop$V8) == 0) {

			# end_1 is between to lines: shift end_1 to the end of 2nd line.

      if (sum(c2c$V3 >= c2c$V2) != 0) {

  			tmp_c2c = c2c[c2c$V1==chro_1 & c2c$V3<=end_1,]

  			end_1 = max(tmp_c2c$V3)

      } else {

  			tmp_c2c = c2c[c2c$V1==chro_1 & c2c$V2<=end_1,]

  			end_1 = max(tmp_c2c$V2)

      }

			if (reverse) {

		    tmptransfostop = c2c[c2c$V1==chro_1 &  ((c2c$V3>=end_1   & c2c$V2<=end_1   & c2c$V8==1) | (c2c$V2>=end_1   & c2c$V3<=end_1   & c2c$V8==-1)),]

			} else {

			  tmptransfostop = c2c[c2c$V1==chro_1 & c2c$V3>=end_1 & c2c$V2<=end_1,]

			}

			if (length(tmptransfostop$V8) == 0) {

				if (!is.null(big_roi)) {

					return(NULL)

				  tmptransfostop = c2c[c2c$V1==chro_1 & c2c$V3>=big_roi$end & c2c$V2<=big_roi$end,]

				} else {

					# return(NULL)

					print(c2c[c2c$V1==chro_1,])

					print(tmptransfostart)

					print(tmptransfostop)

					stop("Never happend condition 2.")

				} 

			} 

		}

		if (length(tmptransfostart$V8) != 1) {

			# tmptransfostart = tmptransfostart[1,]

			# print("many start")

			# print(c2c[c2c$V1==chro_1,])

			tmptransfostart = tmptransfostart[tmptransfostart$V1==chro_1 & tmptransfostart$V3>=begin_1 & tmptransfostart$V2==begin_1,]

			if (length(tmptransfostart$V8) != 1) {

				# return(NULL)

				print(tmptransfostart)

				print(tmptransfostop)

  			stop("Never happend condition 3.")

			}

		}

		if (length(tmptransfostop$V8) != 1) {

			# tmptransfostop = tmptransfostop[length(tmptransfostop$V8),]

			# print("many stop")

			# print(tmptransfostop)

			# print(roi)

		  tmptransfostop = tmptransfostop[tmptransfostop$V1==chro_1 & tmptransfostop$V3==end_1 & tmptransfostop$V2<=end_1,]

			if (length(tmptransfostop$V8) != 1) {

				# return(NULL)

				print(tmptransfostart)

				print(tmptransfostop)

  			stop("Never happend condition 4.")

			}

		}		

		if (tmptransfostart$V7 != tmptransfostop$V7) {

			print(tmptransfostart)

			print(tmptransfostop)		

 			stop("Problem with genome region of interest of strain 1. \nIt is translated over many contigs into strain 2 ref. \nSorry, but you have to redefine your region of interest.")

		}

	} 

  # Deal with strand

  if (tmptransfostart$V8 == 1) {

    begin_2 = tmptransfostart$V9 + (begin_1 - tmptransfostart$V2)

    end_2 = tmptransfostop$V9 + (end_1 - tmptransfostop$V2)

  } else {

    begin_2 = tmptransfostart$V9 - (begin_1 - tmptransfostart$V2)

    end_2 = tmptransfostop$V9 - (end_1 - tmptransfostop$V2)

  }

	# Build returned roi

	roi$strain_ref = strain2

	if (roi$strain_ref == "BY") {

		roi$chr = ROM2ARAB()[[substr(tmptransfostart$V7, 4, 12)]]

	} else {

		roi$chr = config$FASTA_INDEXES[[strain2]][[tmptransfostart$V7]]

	}

  roi$begin = begin_2

  roi$end = end_2

	if (sign(roi$end - roi$begin) == 0) {

		roi$length = 1			

	} else {		

		roi$length = roi$end - roi$begin + sign(roi$end - roi$begin) * 1	

	}

  return(roi)  

}, ex=function(){

	# Define new translate_roi function...

	translate_roi = function(roi, strain2, config) {

		strain1 = roi$strain_ref

		if (strain1 == strain2) {

			return(roi)

		} else {

		  stop("Here is my new translate_roi function...")		

		}	

	}

	# Binding it by uncomment follwing lines.

	# unlockBinding("translate_roi", as.environment("package:nm"))

	# unlockBinding("translate_roi", getNamespace("nm"))

	# assign("translate_roi", translate_roi, "package:nm")

	# assign("translate_roi", translate_roi, getNamespace("nm"))

	# lockBinding("translate_roi", getNamespace("nm"))

	# lockBinding("translate_roi", as.environment("package:nm"))	

})

compute_inter_all_strain_curs = function (# Compute Common Uninterrupted Regions (CUR)

### CURs are regions that can be aligned between the genomes 

diff_allowed = 10, ##<< the maximum indel width allowe din a CUR

min_cur_width = 200, ##<< The minimum width of a CUR

config = NULL, ##<< GLOBAL config variable

plot = FALSE ##<< Plot CURs or not

) {

  get_inter_strain_rois = function(strain1, strain2, diff_allowed = 10, min_cur_width = 200, plot=FALSE) {

  	c2c = get_content(config$C2C_FILES[[paste(strain1, "-", strain2, sep="")]], "table", stringsAsFactors=FALSE)

    # Filtering unagapped

    c2c = c2c[c2c$V6=="-",]

    # filtering some things (chr...)

    # c2c = c2c[c2c$V1 == "chrIV",]

    diff = c2c$V2[-1] - c2c$V3[-length(c2c$V2)] 

    diff2 = c2c$V9[-1] - c2c$V10[-length(c2c$V2)] 

  	# Plot diffs to define a threshold (diff_allowed)

  	# hist(abs(c(diff2, diff)),breaks=c(0:2000, 200000000000), xlim=c(0,100))

    # Filtering

  	indexes_stop = which(abs(diff) > diff_allowed | abs(diff2) > diff_allowed)

  	indexes_start = c(1, indexes_stop[-length(indexes_stop)] + rep(1, length(indexes_stop) -1))

    rois = NULL	

  	for(i in 1:length(indexes_start)) {

  		start = indexes_start[i]

  		stop = indexes_stop[i]

  		sub_c2c = c2c[start:stop,]

  		if (strain1 == "BY") {

  			chr = ROM2ARAB()[[substr(sub_c2c[1,]$V1,4,10)]]

  		} else {			

  			chr = config$FASTA_INDEXES[[strain1]][[sub_c2c[1,]$V1]]

  		}

  		roi = list(chr=chr, begin=sub_c2c[1,]$V2, end=sub_c2c[length(sub_c2c$V1),]$V3, strain_ref=strain1)

  		roi[["length"]] = roi$end - roi$begin

  		if (roi$length >= min_cur_width) {

        rois = dfadd(rois,roi)

  	  }

  		if (length(unique(sub_c2c[,c(1,7,8)])[,2]) != 1) {

  			print("*************** ERROR, non homogenous region! ********************")

  		}

  		# print(i)

  		# print(roi)

  		# print(sub_c2c)

  		# print("________________________________________________________________")

  	}

  	if (plot) {

  		print(paste(length(indexes_stop), "area of interest."))

  	  # Plot rois

  	  genome = get_content(config$FASTA_REFERENCE_GENOME_FILES[[strain1]], "fasta")   

  		plot(0,0, ylim=(c(1,length(genome))), xlim = c(0, max(apply(t(genome), 2, function(chr){length(unlist(chr))}))))

  		for (name in names(genome)) {

  			if (strain1 == "BY") {

  				chr_ref = paste("chr", ARAB2ROM()[[config$FASTA_INDEXES[[strain1]][[name]]]], sep="")

  			} else {			

  				chr_ref = name

  			}

  			y_lev = as.integer(config$FASTA_INDEXES[[strain1]][[name]])

  			lines(c(0,length(unlist(genome[[name]]))), c(y_lev,y_lev))

  			text( length(unlist(genome[[name]]))/2, y_lev, labels = chr_ref)	

  		}

  	  col=1

  	  for (roi_index in 1:length(rois$chr)) {

  			roi = rois[roi_index,]	

  			y_lev = as.integer(roi$chr) + 0.3

  			lines(c(roi$begin,roi$end), c(y_lev,y_lev), col=col)

  			text( mean(c(roi$begin,roi$end)), y_lev, labels = roi_index)

  	  	col = col + 1

  	  }	

  	}

  	return(rois)

  }

	rois = NULL	

	rois_BY_RM = get_inter_strain_rois("BY", "RM", min_cur_width = min_cur_width, diff_allowed = diff_allowed)

	rois_BY_YJM = get_inter_strain_rois("BY", "YJM", min_cur_width = min_cur_width, diff_allowed = diff_allowed)

	for (roi_1_index in 1:length(rois_BY_RM[,1])) {

		roi_1 = rois_BY_RM[roi_1_index,]

		roi_2_candidates = rois_BY_YJM[rois_BY_YJM$chr== roi_1$chr & rois_BY_YJM$begin <= roi_1$end & rois_BY_YJM$end >= roi_1$begin , ] ;  

		# print(length(roi_2_candidates[,1]))

		if (length(roi_2_candidates[,1]) > 0) {

			for(roi_2_index in 1:length(roi_2_candidates[,1])) {

				roi_2 = roi_2_candidates[roi_2_index,]

				roi = list(chr=roi_1$chr, begin=max(roi_1$begin, roi_2$begin), end=min(roi_1$end, roi_2$end), strain_ref="BY")

				roi[["length"]] = roi$end - roi$begin + 1

				if (roi$length >= min_cur_width) {

					# if (length(rois[,1]) == 153) {

					# 	print(paste(length(rois[,1]), roi_1_index, roi_2_index ))

					# 	print(roi_1)

					# 	print(roi_2)

					# 	print(roi)

					# }

			    rois = dfadd(rois,roi)

			  }			

			}

		}

	}

	print(length(rois[,1]))

	print(sum(rois$length))

	rois_1st_round = rois

	rois_2nd_round = NULL	

	rois_RM_YJM = get_inter_strain_rois("RM", "YJM", min_cur_width = min_cur_width, diff_allowed = diff_allowed)

	for (roi_1_index in 1:length(rois_1st_round[,1])) {

		roi_1 = rois_1st_round[roi_1_index,]

		translated_roi_1 = translate_roi(roi_1, "RM", config = config)

		t_b = min(translated_roi_1$begin, translated_roi_1$end)

		t_e = max(translated_roi_1$begin, translated_roi_1$end)

		roi_2_candidates = rois_RM_YJM[rois_RM_YJM$chr== translated_roi_1$chr & rois_RM_YJM$begin <= t_e & rois_RM_YJM$end >= t_b , ] ;  

		if (length(roi_2_candidates[,1]) > 0) {

			for(roi_2_index in 1:length(roi_2_candidates[,1])) {

				roi_2 = roi_2_candidates[roi_2_index,]

				roi = list(chr=translated_roi_1$chr, begin=max(t_b, roi_2$begin), end=min(t_e, roi_2$end), strain_ref="RM")

				roi[["length"]] = roi$end - roi$begin + 1