LBMC » NucleoMiner

count_reads_cur = memoise(function(

marker,

combi,

form,

cur_index,

all_samples, ##<< A table that describe all our samples.

config=NULL ##<< GLOBAL config variable.

) {

	print(paste("Counting reads for ", form, " CUR ", cur_index, " in ", marker, " for [", combi[1], ",", combi[2], "].", sep=""))

	nucs = list()

	for (strain in combi) {

    if (form == "unr") {

  	  nucs[[strain]] = mread.table(paste(config$RESULTS_DIR, "/", strain, "_in_", combi[1], "_", combi[2], "_unr.tab", sep=""), header=TRUE)

    } else {

	    nucs[[strain]] = mread.table(paste(config$RESULTS_DIR, "/", strain, "_wp.tab", sep=""), header=TRUE)

    }

	}

	common_nucs = mread.table(paste(config$RESULTS_DIR, "/", combi[1], "_", combi[2], "_common_", form, ".tab", sep=""), header=TRUE)    

  if (form == "wp") {

    common_nucs = common_nucs[common_nucs$llr_score < 100, ]

  }

	common_nucs = common_nucs[ common_nucs$cur_index == cur_index, ]

	if ( length(common_nucs[,1]) != 0) {

		all_res = apply(t(1:length(common_nucs[,1])), 2, function(i){

    	if (i %% 10 == 1) {

    	  print(paste(i, "/", length(common_nucs[,1]), sep=""))	

      }

    	common_nuc = common_nucs[i,]		

    	res= list()

      # Get nuc infos

      for (strain in combi) {

    	  tmp_nuc_details =  nucs[[strain]][nucs[[strain]]$cur_index == common_nuc$cur_index & nucs[[strain]]$index_nuc == common_nuc[[paste("index_nuc",strain, sep="_")]], ]

        # print(tmp_nuc_details)

    		res[[paste("chr", strain, sep="_")]] = tmp_nuc_details$chr

    		res[[paste("lower_bound", strain, sep="_")]] = tmp_nuc_details$lower_bound

    		res[[paste("upper_bound", strain, sep="_")]] = tmp_nuc_details$upper_bound

    		res[[paste("index_nuc", strain, sep="_")]] = common_nuc[[paste("index_nuc",strain, sep="_")]]

      }	

    	res$cur_index =	common_nuc$cur_index

      # Get normalized reads

    	manip = marker

    	{			

    		for (strain in combi) {				

    			for(sample_id in all_samples[all_samples$marker == manip & all_samples$strain == strain, ]$id) {

            # tmp_filename = paste(config$ALIGN_DIR, "/TF/splited/sample_", sample, "_chr_", res[[paste("chr", strain, sep="_")]], "_splited_sample.tab.gz",sep="")

            if ("tf_input" %in% names(all_samples)) {

              sample_inputs_filename = all_samples$tf_input[all_samples$id==sample_id]

            } else {

              sample_inputs_filename = paste(config$ALIGN_DIR, "/TF/sample_", sample_id, "_TF.txt", sep="")

            }

    				tmp_unnorm_reads = mread.table(sample_inputs_filename, stringsAsFactors=FALSE)					

            names(tmp_unnorm_reads) = c("chr", "pos", "strand", "nb_reads")

    				orig_cur = list(name = "foo", 

                                begin = res[[paste("lower_bound", strain, sep="_")]], 

    														end = res[[paste("upper_bound", strain, sep="_")]], 

    														chr = res[[paste("chr", strain, sep="_")]], 

    														strain_ref = strain)

    				orig_cur$length = orig_cur$end - orig_cur$begin + 1

    				tmp_nuc_reads = filter_tf_inputs(tmp_unnorm_reads, orig_cur$chr, orig_cur$begin, orig_cur$end, orig_cur$length) 

    				res[[paste(strain, manip, sample_id, sep="_")]] = sum(tmp_nuc_reads[,4])

    			}

    		}

    	}

    	res

    })   

		return(all_res)

	} else {

		return(NULL)

	}

})

build_count_table = function(# Build count table for a set of samples.

### This function build a count table for a set of sample.

marker, ##<< The marker that we want to build the count table.

combi, ##<< The combinations of strains that we want to build the count table.

form, ##<< The nucleosome that we want to observe: "wp" for sel;l position and "unr" for UNR.

curs, ##<< The list of CURs

all_samples, ##<< A table that describe all our samples.

config=NULL ##<< GLOBAL config variable.

) {

  dir.create(config$RESULTS_DIR, recursive=TRUE, showWarnings=FALSE)

	all_res = apply(t(1:nrow(curs)), 2, function(cur_index) {

    res = count_reads_cur(marker=marker, combi=combi, form=form, cur_index=cur_index, all_samples, config)

		return(res)

	})

	vec_names = names(all_res[[1]][[1]])

	all_res = data.frame(matrix(unlist(all_res), ncol= length(vec_names), byrow=TRUE))

  names(all_res) = vec_names

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

	write.table(all_res, file=out_filename, row.names=FALSE, quote=FALSE)				

}

extract_maps = memoise(function(# Extract maps from TemplateFilter  outputs

### This function extracts from TemplateFilter outputs for a given CUR./ This is from there that aggregate_intra_strain_nucs and align_inter_strain_nucs fucntions are calles. This is an internal fucntion

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

strains, ##<< The strains for which we want to extract intra strain information.

combis, ##<< All combinations of strains for which we want to extract inter strain information.

all_samples, ##<< A table that describe all our samples.

config=NULL##<< GLOBAL config variable.

) {

  cur = as.list(curs[cur_index,])

  print(paste("cur length", cur$length))

  mnase_replicates = list()

  aggregated_intra_strain_nucs = list()

  for (strain in strains) {

    mnase_replicates[[strain]] = fetch_mnase_replicates(strain, as.list(curs[cur_index,]), all_samples, config)

    aggregated_intra_strain_nucs[[strain]] = aggregate_intra_strain_nucs(mnase_replicates[[strain]])

  }

  aligned_inter_strain_nucs = list()

  for(combi in combis) {

    if (!(length(aggregated_intra_strain_nucs[[combi[1]]][[1]]) == 0 | length(aggregated_intra_strain_nucs[[combi[2]]][[1]]) == 0 )) {

      aligned_inter_strain_nucs[[paste(combi[1], combi[2], sep="_")]] = align_inter_strain_nucs(replicates = list(mnase_replicates[[combi[1]]], mnase_replicates[[combi[2]]]),

                                                                                                wp_nucs_strain_ref1 = aggregated_intra_strain_nucs[[combi[1]]][[1]],

                                                                                                wp_nucs_strain_ref2 = aggregated_intra_strain_nucs[[combi[2]]][[1]], config = config)

    } else {

      print("WARNING! no aggregated_intra_strain_nucs.")

      aligned_inter_strain_nucs[[paste(combi[1], combi[2], sep="_")]] = NULL

    }

  }

  return(list(aggregated_intra_strain_nucs, aligned_inter_strain_nucs))

})

build_maps = function(# Extract maps from TemplateFilter  outputs

### This function extracts from TemplateFilter outputs./ This is from there that aggregate_intra_strain_nucs and align_inter_strain_nucs fucntions are calles. This fucntion write well positionned, fuzzy and both maps in the config$RESULTS_DIR directory.

strains, ##<< The strains for which we want to extract intra strain information.

combis, ##<< The combinations of strains for which we want to extract inter strain information.

all_samples, ##<< A table that describe all our samples.

curs, ##<< The list of CURs

config=NULL##<< GLOBAL config variable.

) {

  dir.create(config$RESULTS_DIR, recursive=TRUE, showWarnings=FALSE)

  ### build wp maps

	glo_results = apply(t(1:1:nrow(curs)), 2, function(cur_index) {

  	dyad_shift = wp_llr = strain_maps = common_nuc_results = intra_llrs = inter_llrs = list()

		print(paste("Collecting maps", cur_index , "/", nrow(curs)))

		extracted_maps = extract_maps(cur_index, strains, combis, all_samples, config)

		aggregated_intra_strain_nucs = extracted_maps[[1]]

		for (strain in strains) {

			print(paste("Collecting maps for", strain, "..."))

			partial_strain_maps = aggregated_intra_strain_nucs[[strain]][[1]]

      print(strain)

      nb_tracks = sum(all_samples$strain == strain & all_samples$marker == "Mnase_Seq")

      tmp_nuc_map = flat_aggregated_intra_strain_nucs(partial_strain_maps, cur_index, nb_tracks)

      wp_indexes = which(tmp_nuc_map$wp == 1)

      tmp_wp_llr = apply(t(wp_indexes), 2, function(wp_index){

        tmp_wp = partial_strain_maps[[wp_index]]

        l = lapply(tmp_wp$nucs, function(tmp_wp_nucs){

          tmp_wp_nucs$original_reads

        })  

        res = llr_score_nvecs(l)

        return(res)

      })

      tmp_dyad_shift = apply(t(wp_indexes), 2, function(wp_index){

        tmp_wp = partial_strain_maps[[wp_index]]  

        c = sapply(tmp_wp$nucs, function(tmp_wp_nucs){

          tmp_wp_nucs$center

        })  

        ds = max(c) - min(c)

        return(ds)

      })

      nb_rep = length(partial_strain_maps[[wp_indexes[1]]]$nucs)

      tmp_nuc_map$wp_llr = rep(NA,nrow(tmp_nuc_map))

      tmp_nuc_map[which(tmp_nuc_map$wp == 1),]$wp_llr = signif(tmp_wp_llr,5)

      tmp_nuc_map$wp_pval = signif(1-pchisq(2*tmp_nuc_map$wp_llr, df=(nb_rep * 2) - 2),5)

      tmp_nuc_map$dyad_shift = rep(NA,nrow(tmp_nuc_map))

      tmp_nuc_map[which(tmp_nuc_map$wp == 1),]$dyad_shift = tmp_dyad_shift

      strain_maps[[strain]] = tmp_nuc_map

      intra_llrs[[strain]] = aggregated_intra_strain_nucs[[strain]][[2]]

		}

    for(combi in combis) {

      print(paste("Collecting results for", combi[1], combi[2], "..."))

      partial_common_nuc_results = extracted_maps[[2]][[paste(combi[1], combi[2], sep="_")]][[1]]

      if (length(partial_common_nuc_results) > 0) {

        partial_common_nuc_results$cur_index = rep(cur_index, length(partial_common_nuc_results[,1]))

        # partial_common_nuc_results$common_wp_pval = signif(1-pchisq(2*partial_common_nuc_results$llr_score, df=2),5)

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

      } else {

        print(paste("No partial_common_nuc_results for cur", cur_index, "ands combi", combi[1], combi[2], "." ))

        common_nuc_results[[paste(combi[1], combi[2], sep="_")]] = data.frame()

      }

      inter_llrs[[paste(combi[1], combi[2], sep="_")]] = extracted_maps[[2]][[paste(combi[1], combi[2], sep="_")]][[2]]

    }

    return(list(strain_maps, common_nuc_results, intra_llrs, inter_llrs))

  })

  glo_results = do.call("rbind", glo_results)

  llrs_intra = do.call("rbind", glo_results[,3])

  llrs_inter = do.call("rbind", glo_results[,4])

  for (strain in strains) {

    outpout_joint_filename = paste(config$RESULTS_DIR, "/", strain, "_llrs.tab", sep="")

    write.table(unlist(llrs_intra[,strain]), file=outpout_joint_filename, row.names=FALSE, quote=FALSE)

  }

  for(combi in combis) {

    outpout_joint_filename = paste(config$RESULTS_DIR, "/", combi[1], "_", combi[2], "_llrs.tab", sep="")

    write.table(unlist(llrs_inter[,paste(combi[1], combi[2], sep="_")]), file=outpout_joint_filename, row.names=FALSE, quote=FALSE)

  }

  glo_strain_maps = do.call("rbind", glo_results[,1])

  for (strain in strains) {

      nuc_map_filename = paste(config$RESULTS_DIR, "/", strain, "_wp.tab", sep="")

      write.table(do.call("rbind", glo_strain_maps[,strain]), file=nuc_map_filename, row.names=FALSE, quote=FALSE)

  }

  glo_common_nuc_results = do.call("rbind", glo_results[,2])

  for(combi in combis) {

    outpout_joint_filename = paste(config$RESULTS_DIR, "/", combi[1], "_", combi[2], "_common_wp.tab", sep="")

    write.table(do.call("rbind", glo_common_nuc_results[,paste(combi[1], combi[2], sep="_")])[, c("cur_index", paste("index_nuc", combi[1], sep="_"), paste("index_nuc", combi[2], sep="_"), "llr_score", "common_wp_pval", "diff")], file=outpout_joint_filename, row.names=FALSE, quote=FALSE)

  }

  ### build fuzzy maps

    wp_maps  = list()

    for (strain in strains) {

      wp_map_filename = paste(config$RESULTS_DIR, "/", strain, "_wp.tab", sep="")

      wp_maps[[strain]] = mread.table(wp_map_filename, header=TRUE)

    }

    glo_fuzzy = list()

    for (strain in strains) {

      tmp_wp = wp_maps[[strain]]

    	glo_fuzzy[[strain]] = apply(t(1:1:nrow(curs)), 2, function(cur_index) {

        print(paste("Building fuzzy for", strain, cur_index, "/", nrow(curs)))

        tmp_wp = tmp_wp[tmp_wp$cur_index == cur_index, ]

        tmp_fuzzy = get_intra_strain_fuzzy(tmp_wp, as.list(curs[cur_index,]), strain, config=config)

        return(tmp_fuzzy)

      })

      fuzzy_map_filename = paste(config$RESULTS_DIR, "/", strain, "_fuzzy.tab", sep="")

      write.table(do.call("rbind", glo_fuzzy[[strain]]), fuzzy_map_filename, row.names=FALSE, quote=FALSE)

    }

  ### build unr maps

  wp_maps  = list()

  for (strain in strains) {

    wp_map_filename = paste(config$RESULTS_DIR, "/", strain, "_wp.tab", sep="")

    wp_maps[[strain]] = mread.table(wp_map_filename, header=TRUE)

  }

  fuzzy_maps = list()

  for (strain in strains) {

    fuzzy_map_filename = paste(config$RESULTS_DIR, "/", strain, "_fuzzy.tab", sep="")

    fuzzy_maps[[strain]] = mread.table(fuzzy_map_filename, header=TRUE)

  }

  common_nuc_results = list()

  for(combi in combis) {

    outpout_joint_filename = paste(config$RESULTS_DIR, "/", combi[1], "_", combi[2], "_common_wp.tab", sep="")

    tmp_common_wp = mread.table(outpout_joint_filename, header=TRUE)

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

  }

  unr = list()

  for(combi in combis) {

    unr[[paste(combi[1], combi[2], sep="_")]] = do.call("rbind", apply(t(1:1:nrow(curs)), 2, function(cur_index) {

      print(paste("Building UNRs", combi[1], combi[2], cur_index, "/", nrow(curs)))

      cur = as.list(curs[cur_index,])

      get_unrs(combi, cur, cur_index, wp_maps, fuzzy_maps, common_nuc_results, config=config)

    }))

    # unr for strain 1

  	unr_map_filename_1 = paste(config$RESULTS_DIR, "/", combi[1], "_in_", combi[1], "_", combi[2], "_unr.tab", sep="")

  	write.table(unr[[paste(combi[1], combi[2], sep="_")]][, c("chr", "lower_bound", "upper_bound", "cur_index", "index_nuc")], file=unr_map_filename_1, row.names=FALSE, quote=FALSE)

    # unr for strain 2

    print(paste(combi[1], combi[2], "translation..."))

    to_tr_unr = unr[[paste(combi[1], combi[2], sep="_")]]

    tr_unr = apply(t(1:length(to_tr_unr[,1])), 2, function(i){

      if (i %% 100 == 1) {print(paste(combi[1], combi[2], "translation... ", i , "/", length(to_tr_unr[,1])))}

      tmp_unr = to_tr_unr[i,]

      tmp_reg = list(name="foo", begin=tmp_unr$lower_bound, end=tmp_unr$upper_bound, chr=tmp_unr$chr, strain_ref = combi[1])

      big_cur = as.list(curs[tmp_unr$cur_index,])

      trs_tmp_reg = translate_cur(tmp_reg, combi[2], config=config, big_cur=big_cur)

      if (!is.null(trs_tmp_reg)) {

        to_tr_unr[i,]$chr = trs_tmp_reg$chr

        to_tr_unr[i,]$lower_bound = min(trs_tmp_reg$begin, trs_tmp_reg$end)

        to_tr_unr[i,]$upper_bound = max(trs_tmp_reg$begin, trs_tmp_reg$end)

        return(to_tr_unr[i,])

        } else {

          return(NULL)

        }

      })

    tr_unr = do.call("rbind", tr_unr)

    # write it.

  	unr_map_filename_2 = paste(config$RESULTS_DIR, "/", combi[2], "_in_", combi[1], "_", combi[2], "_unr.tab", sep="")

  	write.table(tr_unr[, c("chr", "lower_bound", "upper_bound", "cur_index", "index_nuc")], file=unr_map_filename_2, row.names=FALSE, quote=FALSE)

    # joint the two tables

    common_unr = cbind(tr_unr$cur_index, tr_unr$index_nuc, tr_unr$index_nuc)

    common_unr = data.frame(common_unr, stringsAsFactors=FALSE)

    names(common_unr) = c("cur_index", paste("index_nuc", combi[1], sep="_"), paste("index_nuc", combi[2], sep="_"))

    outpout_joint_unr_filename = paste(config$RESULTS_DIR, "/", combi[1], "_", combi[2], "_common_unr.tab", sep="")

    write.table(common_unr, file=outpout_joint_unr_filename, row.names=FALSE, quote=FALSE)

  }

}

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")

})

llr_score_nvecs = structure(function # Likelihood ratio

### Compute the log likelihood ratio of two or more set of value.

(

  xs ##<< list of vectors.

) {

  l = length(xs)

  if (l < 1) {

    return(NA)

  }

  if (l == 1) {

    return(1)

  }

  sumllX = 0

  for (i in 1:l) {

    x = xs[[i]]

  	if (length(x) <= 1) {

  		return(NA)

  	}

    meanX = mean(x)

    sdX = sd(x)

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

    sumllX = sumllX + llX

  }

  meanXglo = mean(unlist(xs))

  sdXglo = sd(unlist(xs))

  llXYZ = sum(log(dnorm(unlist(xs),mean=meanXglo,sd=sdXglo)))

  ratio = sumllX - llXYZ

  return(ratio)

### Returns the log likelihood ratio.

}, ex=function(){

  # LLR 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)

  llr_score_nvecs(list(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 overlapping 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 - tf_outputs$width/2 >= -x_max & tf_outputs$center + tf_outputs$width/2 <=  -x_min,]

	} else {

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

  }

  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.score >= corr_thres,]

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

  i = 1

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

    lb = tf_outputs[i,]$lower_bound

    ub = tf_outputs[i,]$upper_bound

    tf_outputs = tf_outputs[!(tf_outputs$lower_bound <= (ub-ol_bp) & tf_outputs$upper_bound > ub) & !(tf_outputs$upper_bound >= (lb+ol_bp) & tf_outputs$lower_bound < 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.

}

# 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?

# USE_DPLYR = FALSE ##<< Use dplyr lib to filter reads.

# ) {

#   n = names(inputs)

#

#   if (!USE_DPLYR) {

#     if (only_f) {

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

#     } else if (only_r) {

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

#     } else {

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

#     }

#   } else {

#     names(inputs) = c("chr", "pos", "str", "lev")

#     if (only_f) {

#       inputs_out = filter(inputs, chr == chr,  pos >= x_min - nuc_width, str == "F", pos <= x_max + nuc_width)

#     } else if (only_r) {

#       inputs_out = filter(inputs, chr == chr, pos >= x_min - nuc_width, str == "R" & pos <= x_max + nuc_width)

#     } else {

#       inputs_out = filter(inputs, chr == chr, pos >= x_min - nuc_width, pos <= x_max + nuc_width)

#     }

#       # if (!filter_for_coverage) {

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

#       #   inputs = filter(inputs, chr == chr, corrected_inputs_coords >= x_min, corrected_inputs_coords <= x_max)

#       #   inputs$corrected_inputs_coords = NULL

#       # }

#   }

#

#   if (!filter_for_coverage) {

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

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

#   }

#

#   names(inputs_out) = n

#   return(inputs_out)

# ### 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 nucleosomes from 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 (dyad). A chain of nucleosomes is builts across all replicates. Adjacent nucleosomes of the chain are compared two by two. Comparison is based on a log likelihood ratio (LLR1). depending on the LLR1 value nucleosomes are merged (low LLR) or separated (high LLR). Finally the function returns a list of clusters and all computed llr_scores. Each cluster ows an attribute wp for "well positioned". This attribute is set to TRUE if the cluster is composed of exactly one nucleosome 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=...)}.

llr_thres=20, ##<< Log likelihood ratio threshold to decide between merging and separating

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

	llr_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

	llr_score = llr_thres + 1

  # Read nucs from TF outputs

  tf_outs = list()

	i = 1

  for (sample in samples) {

		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

  }

  nb_tracks = length(tf_outs)

	# print(track_readers)

  new_cluster = NULL

  nb_nucs_in_cluster = 0

  nuc_from_track = c()

  for (i in 1:nb_tracks){

    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)) {

			llr_score = llr_score_nvecs(list(current_nuc$original_reads,new_nuc$original_reads))

			llr_scores = c(llr_scores,llr_score)

		}

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

		if (is.na(llr_score)) {

			llr_score = llr_thres + 1

		}

		# Store llr_score

		new_nuc$llr_score = llr_score

	  if (llr_score < llr_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, nb_tracks, min_nuc_center, max_nuc_center)

			}

      # Reinit current cluster composition stuff

      nb_nucs_in_cluster = 0

      nuc_from_track = c()

      for (i in 1:nb_tracks){

        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, nb_tracks, min_nuc_center, max_nuc_center)

  }

	return(list(clusters, llr_scores))

### Returns a list of clusterized nucleosomes, and all computed llr 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

cur_index, ##<< the index of the roi involved

nb_tracks=3 ##<< the number of replicates

) {

	if  (length(partial_strain_maps) == 0 ){

		print(paste("Empty partial_strain_maps for roi", cur_index, "ands current 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[["cur_index"]] = cur_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)

			tmp_nuc_as_list[["nb_nucs"]] = length(tmp_nuc$nucs)

			if (tmp_nuc$wp) {

        for (i in 1:(nb_tracks-1)) {

  				tmp_nuc_as_list[[paste("llr", i, sep="_")]] = signif(tmp_nuc$nucs[[i + 1]]$llr_score,5)          

        }

			} else {

        for (i in 1:(nb_tracks-1)) {

  				tmp_nuc_as_list[[paste("llr", i, sep="_")]] = NA

        }

			}

      return(tmp_nuc_as_list)

    })

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

	}

  return(data.frame(lapply(data.frame(tmp_strain_maps, stringsAsFactors=FALSE), unlist), stringsAsFactors=FALSE))

### 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 nucleosomes 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.

llr_thres=100, ##<< Log likelihood ratio threshold to decide between merging and separating

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

	llr_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_cur = replicates[[1]][[1]]$roi

  orig_big_cur = replicates[[1]][[1]]$orig_roi

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

		tmp_begin = big_cur$begin

		big_cur$begin =  big_cur$end

		big_cur$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]]

  }

	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_cur(roi_strain_ref1, strain_ref2, big_cur=orig_big_cur, 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_cur(nuc_strain_ref2_to_roi, strain_ref1, big_cur=orig_big_cur, 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))

								llr_score = llr_score_nvecs(list(reads_strain_ref1, reads_strain_ref2))

								llr_scores = c(llr_scores, llr_score)

								# Filtering on LLR Score

                if (llr_score < llr_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

									# 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

									# common

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

                  tmp_nuc[["common_wp_pval"]] = signif(1-pchisq(2*tmp_nuc[["llr_score"]], df=2),5)

									tmp_nuc[["diff"]] = abs(abs((roi_strain_ref2$begin - roi_strain_ref2$end)) / 2 - abs((nuc_strain_ref2$lower_bound - nuc_strain_ref2$upper_bound)) / 2)

									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, llr_scores))

	} else {

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

		return(NULL)

	}

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

}, ex=function(){

    # Define new translate_cur function...

    translate_cur = function(roi, strain2, big_cur=NULL, config=NULL) {

      return(roi)

    }

    # Binding it by uncomment follwing lines.

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

    unlockBinding("translate_cur", getNamespace("nucleominer"))

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

    assign("translate_cur", translate_cur, getNamespace("nucleominer"))

    lockBinding("translate_cur", getNamespace("nucleominer"))

    lockBinding("translate_cur", 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_cur(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))

})