LBMC » NucleoMiner

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?

USE_DPLYR = TRUE ##<< 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))

})

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_cur(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[,1:4], region2[,1:4])

  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 themselves.

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

) {

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

  if (nrow(regions) == 0) {return(regions)}

  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)] <= 1, 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]])

      cur_index = regions$cur_index[vec_beg_1[i]]

      data.frame(list(chr=chr, lower_bound=lower_bound, upper_bound=upper_bound, cur_index=cur_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.

# cur_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$cur_index == cur_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 elaborated call to translate_cur.

regions, ##<< Regions to be translated.

combi, ##<< Combination of strains.

cur_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_cur =  roi

    trs_tmp_regions_ref2 = translate_cur(tmp_regions_ref2, combi[1], config = config, big_cur = big_cur)

    if (is.null(trs_tmp_regions_ref2)) {

      return(NULL)

    } else {

      return(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), cur_index=cur_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

) {

  if (is.null(regions)) {

    return(NULL)

  } else {

    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)

  }

}

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

### The fucntion is no more necessary since we remove "big_cur" bug in translate_cur 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_cur(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_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 == "unr") {

		  out_filename = paste(config$RESULTS_DIR, "/",combi[1],"_",combi[2],"_",manip,"_unr_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 == "wpunr") {

		 	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,"_unr_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 = "unr"

      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 normalizes data, fit data to GLM model with and without interaction term and compares the two models.

snep_design, ##<< The design to consider.

reads ##<< The data to consider.

) {

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

}