commit 729c934eba59d0044626c5da6d67bfbeb993fc01
Author: Florent Chuffart <florent.chuffart@ens-lyon.fr>
Date:   Thu Feb 27 16:51:50 2014 +0100

    Adding chain visualiZation feature in watch_sample function

diff --git a/src/R/nucleominer.R b/src/R/nucleominer.R
index be0a2d8..283d032 100644
--- a/src/R/nucleominer.R
+++ b/src/R/nucleominer.R
@@ -124,7 +124,7 @@ y ##<< Second vector.
   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)
+  ratio = llX + llY - llXY
   return(ratio)
 ### Returns the likelihood ratio.
 }, ex=function(){
@@ -267,7 +267,7 @@ strand ##<< The original 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.
+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) {
@@ -311,7 +311,7 @@ coord_max=20000000 ##<< A too big value to be a coord for a nucleosome lower bou
   indexes = c()
   track_readers = c()
   current_nuc = NULL
-	lod_score = lod_thres - 1 
+	lod_score = lod_thres + 1 
   # Read nucs from TF outputs  
   tf_outs = list()
 	i = 1
@@ -361,11 +361,11 @@ coord_max=20000000 ##<< A too big value to be a coord for a nucleosome lower bou
 		}
 		# 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 
+			lod_score = lod_thres + 1 
 		}
-	  if (lod_score > lod_thres) {
-			# Store lod_score
-			new_nuc$lod_score = lod_score 			
+		# 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) {
@@ -456,7 +456,7 @@ 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.
+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.	
 ) {
@@ -544,7 +544,7 @@ config=NULL, ##<< GLOBAL config variable
 								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) {
+								if (lod_score < lod_thres) {
 									tmp_nuc = list()
 									# strain_ref1
 									tmp_nuc[[paste("chr_", strain_ref1, sep="")]] = chr
@@ -1730,6 +1730,7 @@ plot = FALSE ##<< Plot CURs or not
 		}
 	}
 
+
 	print(length(rois_3rd_round[,1]))
 	print(sum(rois_3rd_round$length))
 	rois = rois_3rd_round
@@ -2061,6 +2062,7 @@ plot_common_nucs = TRUE,  ##<< Plot (or not) aligned reads.
 plot_anovas = FALSE,  ##<< Plot (or not) scatter for each nuc.
 plot_anova_boxes =  FALSE,  ##<< Plot (or not) boxplot for each nuc.
 plot_wp_nucs_4_nonmnase = FALSE,  ##<< Plot (or not) clusters for non inputs samples.
+plot_chain = FALSE,  ##<< Plot (or not) clusterised nuceosomes between mnase samples.
 aggregated_intra_strain_nucs = NULL, ##<< list of aggregated intra strain nucs. If NULL, it will be computed. 
 aligned_inter_strain_nucs = NULL, ##<< list of aligned inter strain nucs. If NULL, it will be computed.
 height = 10, ##<< Number of reads in per million read for each sample, graphical parametre for the y axis.
@@ -2196,34 +2198,65 @@ config=NULL ##<< GLOBAL config variable
 		    }
 			}
 	  }
+    
 	  # Plot wp nucs
-		if ((plot_wp_nucs_4_nonmnase | sample$marker == "Mnase_Seq") & (plot_wp_nucs | plot_common_nucs)) {
+		if ((plot_wp_nucs_4_nonmnase | sample$marker == "Mnase_Seq") & (plot_wp_nucs | plot_common_nucs | plot_chain)) {
 			if (samples[[1]]$marker == "Mnase_Seq") {
 				if (is.null(aggregated_intra_strain_nucs)) {
 	  			wp_nucs = aggregate_intra_strain_nucs(samples)[[1]]			
+          foo <<- wp_nucs
 				} else {
 					wp_nucs = aggregated_intra_strain_nucs[[replicate_rank]]
 				}
 		  } else {
   			wp_nucs = replicates_wp_nucs[[replicate_rank-2]]
 		  }
-			replicates_wp_nucs[[replicate_rank]] = wp_nucs
-			for (wp_nuc in wp_nucs) {
-				if (wp_nuc$wp){
-					rect(sign(x_min) * wp_nuc$lower_bound + shift, y_min, sign(x_min) * wp_nuc$upper_bound + shift, y_max, col=adjustcolor(2, alpha.f = 0.1), border=1)
-					all_original_reads = c()
-					for(initial_nuc in wp_nuc$nucs) {
-						all_original_reads = c(all_original_reads, initial_nuc$original_reads)						
-					}
-					delta_x = wp_nuc$lower_bound:wp_nuc$upper_bound
-					if (FALSE) {
-					  rect(sign(x_min) * wp_nuc$lower_bound + shift, y_min, sign(x_min) * wp_nuc$upper_bound + shift, y_max, col="#EEEEEE", border=1)
-				  }
-					if (plot_wp_nuc_model) {
-					  lines(sign(x_min) * delta_x + shift, dnorm(delta_x, mean(all_original_reads), sd(all_original_reads)) * length(all_original_reads) * height/5 + y_min, col=1)				
-				  }
-				}
-			}
+      
+      if (plot_chain) {
+        tf_nucs = lapply(wp_nucs, function(nuc) {      
+          bar = apply(t(nuc$nucs), 2, function(tmp_nuc){
+            tmp_nuc = tmp_nuc[[1]]
+            tmp_nuc$inputs = NULL
+            tmp_nuc$original_reads = NULL
+            tmp_nuc$wp = nuc$wp
+            # print(tmp_nuc)
+            return(tmp_nuc)
+            }) 
+            return(do.call(rbind, bar))
+          })
+        tf_nucs = data.frame(do.call(rbind, tf_nucs))
+
+        points((unlist(tf_nucs$lower_bound) + unlist(tf_nucs$upper_bound)) / 2 , y_min + (unlist(tf_nucs$track) - 0.5) * delta_y/nb_rank, cex=4, pch=16, col="white")
+        points((unlist(tf_nucs$lower_bound) + unlist(tf_nucs$upper_bound)) / 2 , y_min + (unlist(tf_nucs$track) - 0.5) * delta_y/nb_rank, cex=4)
+        text((unlist(tf_nucs$lower_bound) + unlist(tf_nucs$upper_bound)) / 2 , y_min + (unlist(tf_nucs$track) - 0.5) * delta_y/nb_rank, 1:nrow(tf_nucs))
+
+        xs = (1:nrow(tf_nucs) - 1) * (x_max_glo - x_min_glo) / (nrow(tf_nucs)) + lb
+        points(xs, rep(0, nrow(tf_nucs)), cex=2, pch=16, col="white")
+        points(xs, rep(0, nrow(tf_nucs)), cex=2, col = unlist(tf_nucs$wp)+2)
+        text(xs, 0, 1:nrow(tf_nucs), cex=0.5)
+        xs = (1:nrow(tf_nucs) - 1.5) * (x_max_glo - x_min_glo) / (nrow(tf_nucs)) + lb
+        text(xs, 0, signif(unlist(tf_nucs$lod_score)), cex=0.5, col=(unlist(tf_nucs$lod_score) < 150) + 2)  
+      }
+  
+      if (plot_wp_nucs | plot_common_nucs ) {
+    		replicates_wp_nucs[[replicate_rank]] = wp_nucs
+  			for (wp_nuc in wp_nucs) {
+  				if (wp_nuc$wp){
+  					rect(sign(x_min) * wp_nuc$lower_bound + shift, y_min, sign(x_min) * wp_nuc$upper_bound + shift, y_max, col=adjustcolor(2, alpha.f = 0.1), border=1)
+  					all_original_reads = c()
+  					for(initial_nuc in wp_nuc$nucs) {
+  						all_original_reads = c(all_original_reads, initial_nuc$original_reads)						
+  					}
+  					delta_x = wp_nuc$lower_bound:wp_nuc$upper_bound
+  					if (FALSE) {
+  					  rect(sign(x_min) * wp_nuc$lower_bound + shift, y_min, sign(x_min) * wp_nuc$upper_bound + shift, y_max, col="#EEEEEE", border=1)
+  				  }
+  					if (plot_wp_nuc_model) {
+  					  lines(sign(x_min) * delta_x + shift, dnorm(delta_x, mean(all_original_reads), sd(all_original_reads)) * length(all_original_reads) * height/5 + y_min, col=1)				
+  				  }
+  				}
+  			}
+      }
 		}
 	}	
   
