Mercurial > repos > bzeitouni > svdetect

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/svdetect/BAM_preprocessingPairs.pl	Thu Jul 12 12:10:27 2012 -0400
@@ -0,0 +1,340 @@
+#!/usr/bin/perl -w
+
+use strict;
+use warnings;
+use Getopt::Std;
+my $version = '0.4b_galaxy';
+
+my $SAMTOOLS_BIN_DIR="/bioinfo/local/samtools";
+
+my %opts = ( t=>1, p=>1, n=>1000000, f=>3, s=>0, S=>10000, o=>"." );
+
+getopts('dt:p:n:f:s:S:o:b:l:x:N:', \%opts); #GALAXY
+
+my $working_dir=($opts{o} ne ".")? $opts{o}:"working directory";
+
+my $pt_bad_mates_file=$opts{b};  #GALAXY
+my $pt_log_file=$opts{l}; #GALAXY
+my $pt_good_mates_file=$opts{x} if($opts{d}); #GALAXY
+
+
+die(qq/
+
+Description:
+
+    Preprocessing of mates to get anomalously mapped mate-pair\/paired-end reads as input
+    for SVDetect.
+
+    From all pairs mapped onto the reference genome, this script outputs abnormal pairs:
+        - mapped on two different chromosomes
+        - with an incorrect strand orientation and\/or pair order
+        - with an insert size distance +- sigma threshold
+    into a file <prefix.ab.bam\/sam> sorted by read names
+
+    -BAM\/SAM File input format only.
+
+    Version : $version
+    SAMtools required for BAM files
+
+
+Usage:   BAM_preprocessingPairs.pl [options] <all_mate_file.sorted.bam\/sam>
+
+Options: -t BOOLEAN   read type: =1 (Illumina), =0 (SOLiD) [$opts{t}]
+         -p BOOLEAN   pair type: =1 (paired-end), =0 (mate-pair)  [$opts{p}]
+         -n INTEGER   number of pairs for calculating mu and sigma lengths [$opts{n}]
+	 -s INTEGER   minimum value of ISIZE for calculating mu and sigma lengths [$opts{s}]
+	 -S INTEGER   maximum value of ISIZE for calculating mu and sigma lengths [$opts{S}]
+         -f REAL      minimal number of sigma fold for filtering pairs [$opts{f}]
+         -d           dump normal pairs into a file [<prefix.norm.bam\/sam>] (optional)
+	 -o STRING    output directory [$working_dir]
+
+\n/) if (@ARGV == 0 && -t STDIN);
+
+unless (-d $opts{o}){
+	mkdir $opts{o} or die;
+}
+$opts{o}.="/" if($opts{o}!~/\/$/);
+
+my $mates_file=shift(@ARGV);
+
+$mates_file=readlink($mates_file);
+
+my $bad_mates_file=(split(/\//,$mates_file))[$#_];
+
+if($bad_mates_file=~/.(s|b)am$/){
+    $bad_mates_file=~s/.(b|s)am$/.ab.sam/;
+    $bad_mates_file=$opts{o}.$bad_mates_file;
+}
+
+else{
+    die "Error: mate_file with the extension <.bam> or <.sam> needed !\n";
+}
+
+my $good_mates_file;
+if($opts{d}){
+    $good_mates_file=(split(/\//,$mates_file))[$#_];
+    $good_mates_file=~s/.(b|s)am$/.norm.sam/;
+    $good_mates_file=$opts{o}.$good_mates_file;
+}
+
+my $log_file=$opts{o}.$opts{N}.".svdetect_preprocessing.log"; #GALAXY
+
+#------------------------------------------------------------------------------#
+#Calculate mu and sigma
+
+open LOG,">$log_file" or die "$0: can't open ".$log_file.":$!\n";
+
+print LOG "\# Calculating mu and sigma lengths...\n";
+print LOG "-- file=$mates_file\n";
+print LOG "-- n=$opts{n}\n";
+print LOG "-- ISIZE min=$opts{s}, max=$opts{S}\n";
+
+my ($record, $sumX,$sumX2) = (0,0,0);
+my $warn=$opts{n}/10;
+my $prev_pair="FIRST";
+
+my $bam=($mates_file =~ /.bam$/)? 1:0;
+
+if($bam){
+    open(MATES, "${SAMTOOLS_BIN_DIR}/samtools view $mates_file |") or die "$0: can't open ".$mates_file.":$!\n";
+}else{
+    open MATES, "<".$mates_file or die "$0: can't open ".$mates_file.":$!\n";
+}
+
+while(<MATES>){
+
+    my @t=split;
+
+    next if ($t[0]=~/^@/);
+
+    my $current_pair=$t[0];
+    next if($current_pair eq $prev_pair);
+    $prev_pair=$current_pair;
+
+    my ($chr1,$chr2,$length)=($t[2],$t[6],abs($t[8]));
+
+    next if ($chr1 eq "*" || $chr2 eq "*");
+    next if ($length<$opts{s} || $length>$opts{S}) ;
+
+    if($chr2 eq "="){
+
+        $sumX += $length;							#add to sum and sum^2 for mean and variance calculation
+	$sumX2 += $length*$length;
+        $record++;
+    }
+
+    if($record>$warn){
+	print LOG "-- $warn pairs analysed\n";
+        $warn+=$warn;
+    }
+
+    last if ($record>$opts{n});
+
+}
+close (MATES);
+
+$record--;
+my $mu = $sumX/$record;
+my $sigma = sqrt($sumX2/$record - $mu*$mu);
+
+print LOG "-- Total : $record pairs analysed\n";
+print LOG "-- mu length = ".decimal($mu,1).", sigma length = ".decimal($sigma,1)."\n";
+
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+#Preprocessing pairs
+
+$warn=100000;
+
+$record=0;
+my %count=( ab=>0, norm=>0, chr=>0, sense=>0, dist=>0, unmap=>0);
+
+my $read_type=($opts{t})? "Illumina":"SOLiD";
+my $pair_type=($opts{p})? "paired-end":"mate-paired";
+
+print LOG "\# Preprocessing pairs...\n";
+print LOG "-- file= $mates_file\n";
+print LOG "-- type= $read_type $pair_type reads\n";
+print LOG "-- sigma threshold= $opts{f}\n";
+print LOG "-- using ".decimal($mu-$opts{f}*$sigma,4)."-".decimal($mu+$opts{f}*$sigma,4)." as normal range of insert size\n";
+
+my @header;
+
+if($bam){
+    open(HEADER, "${SAMTOOLS_BIN_DIR}/samtools view -H $mates_file |") or die "$0: can't open ".$mates_file.":$!\n";
+    @header=<HEADER>;
+    close HEADER;
+    open(MATES, "${SAMTOOLS_BIN_DIR}/samtools view $mates_file |") or die "$0: can't open ".$mates_file.":$!\n";
+}else{
+    open MATES, "<".$mates_file or die "$0: can't open ".$mates_file.":$!\n";
+}
+
+open AB, ">$bad_mates_file" or die "$0: can't write in the output: $bad_mates_file :$!\n";
+print AB @header if($bam);
+
+if($opts{d}){
+    open NORM, ">$good_mates_file" or die "$0: can't write in the output: $good_mates_file :$!\n";
+    print NORM @header if($bam);
+}
+
+$prev_pair="FIRST";
+my $prev_bad;
+
+while(<MATES>){
+
+    my @t=split;
+    my $bad=0;
+
+    if ($t[0]=~/^@/){
+        print AB;
+        print NORM if ($opts{d});
+        next;
+    }
+
+    my $current_pair=$t[0];
+    if($current_pair eq $prev_pair){
+        next if($prev_bad==-1);
+        if($prev_bad){
+            print AB;
+        }elsif(!$prev_bad){
+            print NORM if($opts{d});
+        }
+        next;
+    }
+
+    $prev_pair=$current_pair;
+
+    my ($chr1,$chr2,$pos1,$pos2,$length)=($t[2],$t[6],$t[3],$t[7], abs($t[8]));
+
+    if ($chr1 eq "*" || $chr2 eq "*"){
+        $prev_bad=-1;
+        $count{unmap}++;
+        $record++;
+        next;
+
+    }
+
+    my $strand1 = (($t[1]&0x0010))? 'R':'F';
+    my $strand2 = (($t[1]&0x0020))? 'R':'F';
+    my $order1  = (($t[1]&0x0040))? '1':'2';
+    my $order2  = (($t[1]&0x0080))? '1':'2';
+
+    if($order1 == 2){
+        ($strand1,$strand2)=($strand2,$strand1);
+        ($chr1,$chr2)=($chr2,$chr1);
+        ($pos1,$pos2)=($pos2,$pos1);
+        ($order1,$order2)=($order2,$order1);
+    }
+
+    my $sense=$strand1.$strand2;
+
+    if($chr1 ne "=" && $chr2 ne "="){
+        $bad=1;
+        $count{chr}++;
+    }
+
+    if($opts{p}){ #paired-end
+        if(!(($sense eq "FR" && $pos1<$pos2) || ($sense eq "RF" && $pos2<$pos1))){
+            $bad=1;
+            $count{sense}++;
+        }
+    }else{ #mate-pair
+        if($opts{t}){ #Illumina
+            if(!(($sense eq "FR" && $pos2<$pos1) || ($sense eq "RF" && $pos1<$pos2))){
+            $bad=1;
+            $count{sense}++;
+            }
+        }else{ #SOLiD
+            if(!(($sense eq "FF" && $pos2<$pos1) || ($sense eq "RR" && $pos1<$pos2))){
+            $bad=1;
+            $count{sense}++;
+            }
+        }
+    }
+
+    if(($chr1 eq "=" || $chr2 eq "=") && ($length <$mu - $opts{f}*$sigma || $length>$mu + $opts{f}*$sigma)){
+        $bad=1;
+        $count{dist}++;
+    }
+
+    if($bad){
+        print AB;
+        $count{ab}++;
+        $prev_bad=$bad;
+    }else{
+        print NORM if ($opts{d});
+        $count{norm}++;
+        $prev_bad=$bad;
+    }
+
+    $record++;
+
+    if($record>$warn){
+        print LOG "-- $warn pairs analysed\n";
+        $warn+=100000;
+    }
+}
+
+close AB;
+close NORM if($opts{d});
+
+print LOG "-- Total : $record pairs analysed\n";
+print LOG "-- $count{unmap} pairs whose one or both reads are unmapped\n";
+print LOG "-- ".($count{ab}+$count{norm})." mapped pairs\n";
+print LOG "---- $count{ab} abnormal mapped pairs\n";
+print LOG "------ $count{chr} pairs mapped on two different chromosomes\n";
+print LOG "------ $count{sense} pairs with incorrect strand orientation and\/or pair order\n";
+print LOG "------ $count{dist} pairs with incorrect insert size distance\n";
+print LOG "--- $count{norm} correct mapped pairs\n";
+
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+#OUTPUT
+
+if($bam){
+
+    my $bam_file=$bad_mates_file;
+    $bam_file=~s/.sam$/.bam/;
+    print LOG "\# Converting sam to bam for abnormal mapped pairs\n";
+    system("${SAMTOOLS_BIN_DIR}/samtools view -bS $bad_mates_file > $bam_file 2>".$opts{o}."samtools.log");
+    unlink($bad_mates_file);
+    print LOG "-- output created: $bam_file\n";
+
+    system "rm $pt_bad_mates_file ; ln -s $bam_file $pt_bad_mates_file"; #GALAXY
+
+    if($opts{d}){
+        $bam_file=$good_mates_file;
+        $bam_file=~s/.sam$/.bam/;
+        print LOG "\# Converting sam to bam for correct mapped pairs\n";
+        system("${SAMTOOLS_BIN_DIR}/samtools view -bS $good_mates_file > $bam_file 2>".$opts{o}."samtools.log");
+        unlink($good_mates_file);
+        print LOG "-- output created: $bam_file\n";
+
+	system "rm $pt_good_mates_file ; ln -s $bam_file $pt_good_mates_file"; #GALAXY
+
+    }
+
+}
+
+else{
+    print LOG "-- output created: $bad_mates_file\n";
+    print LOG "-- output created: $good_mates_file\n" if($opts{d});
+}
+
+close LOG;
+
+system "rm $pt_log_file ; ln -s $log_file $pt_log_file"; #GALAXY
+
+
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub decimal{
+
+  my $num=shift;
+  my $digs_to_cut=shift;
+
+  $num=sprintf("%.".($digs_to_cut-1)."f", $num) if ($num=~/\d+\.(\d){$digs_to_cut,}/);
+
+  return $num;
+}
+#------------------------------------------------------------------------------#
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/svdetect/BAM_preprocessingPairs.xml	Thu Jul 12 12:10:27 2012 -0400
@@ -0,0 +1,77 @@
+<tool id="svdetect_preprocessing" name="BAM preprocessing">
+
+  <description>to get abnormal pairs</description>
+
+  <command interpreter="perl"> BAM_preprocessingPairs.pl -t '$readType' -p '$pairType' -n '$nbrePair' -s '$isizeMin' -S '$isizeMax' -f '$foldPair' -o $__new_file_path__/svdetect -b '$abBAM' -l '$log' -N $sample_name
+	#if $newBam.pairNormal=="yes"
+	-d -x '$normBAM'
+	#end if
+	'$inputBam'
+  </command>
+
+  <inputs>
+    <param name="sample_name" type="text" value="sample" label="Sample Name"/>
+    <param name="inputBam" type="data" format="bam" label="BAM input file"/>
+    <param name="readType" type="select" label="Read type">
+	<option value="1">Illumina</option>
+	<option value="0">SOLiD</option>
+   </param>
+   <param name="pairType" type="select" label="Library type">
+	<option value="1">Paired-end</option>
+	<option value="0">Mate-Pair</option>
+   </param>
+   <conditional name="newBam">
+   <param name="pairNormal" type="select" label="Do you want an additional bam file listing concordant mapped pairs?" help="Dump normal pairs into a file sample_name.norm.bam/sam">
+	<option value="no">No</option>
+	<option value="yes">Yes</option>
+   </param>
+    <when value="yes">
+	  <!-- do nothing here -->
+    </when>
+    <when value="no">
+	 <!-- do nothing here -->
+    </when>
+   </conditional>
+   <param name="nbrePair" value="1000000" type="integer" size="30" label="Number of pairs for calculating mu (µ) and sigma (σ) lengths"/>
+   <param name="isizeMin" value="0" type="integer" size="30" label="Minimum value of ISIZE for calculating mu (µ) and sigma (σ) lengths"/>
+   <param name="isizeMax" value="10000" type="integer" size="30" label="Maximum value of ISIZE for calculating mu (µ)and sigma( σ) lengths"/>
+   <param name="foldPair" value="3" type="float" size="30" label="Minimal number of sigma (σ) fold for filtering pairs"/>
+  </inputs>
+
+  <outputs>
+    <data format="bam" name="abBAM" label="${$sample_name}.ab.bam"/>
+    <data format="txt" name="log" label="${$sample_name}.svdetect_preprocessing.log"/>
+    <data format="bam" name="normBAM" label="${$sample_name}.norm.bam">
+    <filter>newBam['pairNormal'] == 'yes'</filter>
+    </data>
+  </outputs>
+
+  <help>
+
+**What it does**
+
+Bam_preprocessingPairs - Version 0.4b
+
+Preprocessing of mates to get anomalously mapped mate-pair/paired-end reads as input for SVDetect.
+
+From all pairs mapped onto the reference genome, this script outputs abnormal pairs:
+
+	* mapped on two different chromosomes
+	* with an incorrect strand orientation and/or pair order
+	* with an insert size distance +- sigma threshold
+
+into a file prefix.ab.bam/sam sorted by read names
+
+-BAM/SAM File input format only.
+
+SAMtools required for BAM files
+
+-----
+
+.. class:: infomark
+
+Contact Bruno Zeitouni (svdetect@curie.fr) for any questions or concerns about the Galaxy implementation of SVDetect.
+
+  </help>
+
+</tool>
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/svdetect/SVDetect_compare.pl	Thu Jul 12 12:10:27 2012 -0400
@@ -0,0 +1,716 @@
+#!/usr/bin/perl -w
+
+=pod
+
+=head1 NAME
+
+SVDetect Compare for Galaxy
+
+Version: 0.8 for Galaxy
+
+=head1 SYNOPSIS
+
+SVDetect_compare.pl links2compare -conf <configuration_file> [-help] [-man]
+
+=cut
+
+# -------------------------------------------------------------------
+
+use strict;
+use warnings;
+
+use Pod::Usage;
+use Getopt::Long;
+
+use Config::General;
+use Tie::IxHash;
+
+#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::#
+#PARSE THE COMMAND LINE
+my %OPT;
+GetOptions(\%OPT,
+	   'conf=s',
+	   'out1=s', #GALAXY
+	   'out2=s', #GALAXY
+	   'out3=s', #GALAXY
+	   'out4=s', #GALAXY
+	   'out5=s', #GALAXY
+	   'out6=s', #GALAXY
+	   'out7=s', #GALAXY
+	   'out8=s', #GALAXY
+	   'out9=s', #GALAXY
+	   'l=s', #GALAXY
+	   'N=s', #GALAXY
+	   'help',
+           'man'
+	  );
+
+pod2usage() if $OPT{help};
+pod2usage(-verbose=>2) if $OPT{man};
+pod2usage(-message=> "$!", -exitval => 2) if (!defined $OPT{conf});
+
+
+pod2usage() if(@ARGV<1);
+
+tie (my %func, 'Tie::IxHash',links2compare=>\&links2compare);
+
+foreach my $command (@ARGV){
+    pod2usage(-message=> "Unknown command \"$command\"", -exitval => 2) if (!defined($func{$command}));
+}
+#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::#
+
+#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::#
+#READ THE CONFIGURATION FILE
+my $conf=Config::General->new(    -ConfigFile        => $OPT{conf},
+                                  -Tie => "Tie::IxHash",
+                                  -AllowMultiOptions => 1,
+				  -LowerCaseNames    => 1,
+				  -AutoTrue => 1);
+my %CONF= $conf->getall;
+validateconfiguration(\%CONF);							#validation of the configuration parameters
+
+
+my $SAMTOOLS_BIN_DIR="/bioinfo/local/samtools"; #GALAXY
+my $BEDTOOLS_BIN_DIR="/bioinfo/local/BEDTools/bin"; #GALAXY
+
+my $pt_log_file=$OPT{l}; #GALAXY
+my $log_file=$CONF{general}{output_dir}.$OPT{N}.".svdetect_compare.log"; #GALAXY
+open LOG,">$log_file" or die "$0: can't open ".$log_file.":$!\n";#GALAXY
+
+my @pt_sv_file=($OPT{out1},$OPT{out2},$OPT{out3}) if($OPT{out1}); #GALAXY common,sample,reference
+my @pt_circos_file=($OPT{out4},$OPT{out5},$OPT{out6}) if($OPT{out4}); #GALAXY common,sample,reference
+my @pt_bed_file=($OPT{out7},$OPT{out8},$OPT{out9}) if($OPT{out7}); #GALAXY common,sample,reference
+
+$CONF{compare}{sample_link_file}=readlink($CONF{compare}{sample_link_file});#GALAXY
+$CONF{compare}{sample_link_file}=~s/.sv.txt//; #GALAXY
+
+$CONF{compare}{reference_link_file}=readlink($CONF{compare}{reference_link_file});#GALAXY
+$CONF{compare}{reference_link_file}=~s/.sv.txt//; #GALAXY
+
+#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::#
+
+#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::#
+#COMMAND EXECUTION
+foreach my $command (@ARGV){
+    &{$func{$command}}();
+}
+print LOG "-- end\n";
+
+close LOG;#GALAXY
+system "rm $pt_log_file ; ln -s $log_file $pt_log_file"; #GALAXY
+
+exit(0);
+#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::#
+
+#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::#
+#FUNCTIONS
+
+# -----------------------------------------------------------------------------#
+#MAIN FUNCTION number 5:Comparison between samples, common or specific links
+sub links2compare{
+
+    my @compare_files;
+
+    compareSamples($CONF{general}{output_dir},
+		   $CONF{compare}{list_samples},
+		   $CONF{compare}{sample_link_file},
+		   $CONF{compare}{reference_link_file},
+		   $CONF{compare}{min_overlap},
+		   $CONF{compare}{same_sv_type},
+		   \@compare_files);
+
+    my $pt_ind=0;
+
+    for my $input_file (@compare_files){
+
+	$input_file=$CONF{general}{output_dir}.$input_file;
+
+	my $output_file=$input_file;
+	$output_file=~s/unique$/compared/;
+
+	sortLinks($input_file, $output_file,1);
+
+	if($CONF{compare}{circos_output}){
+	    links2segdup($CONF{circos}{organism_id},
+			 $CONF{circos}{colorcode},
+			 $output_file,
+			 $output_file.".segdup.txt");
+	    system "rm $pt_circos_file[$pt_ind]; ln -s $output_file.segdup.txt $pt_circos_file[$pt_ind]" if (defined $pt_circos_file[$pt_ind]); #GALAXY
+	}
+
+	if($CONF{compare}{bed_output}){
+	links2bedfile($CONF{compare}{read_lengths},
+		      $CONF{bed}{colorcode},
+		      $output_file,
+		      $output_file.".bed");
+	system "rm $pt_bed_file[$pt_ind]; ln -s $output_file.bed $pt_bed_file[$pt_ind]" if (defined $pt_bed_file[$pt_ind]); #GALAXY
+	}
+
+	if($CONF{compare}{sv_output}){
+
+	    links2SVfile ($output_file, $output_file.".sv.txt");
+	    system "rm $pt_sv_file[$pt_ind]; ln -s $output_file.sv.txt $pt_sv_file[$pt_ind]" if (defined $pt_sv_file[$pt_ind]); #GALAXY
+	}
+	$pt_ind++;
+
+    }
+    unlink(@compare_files);
+
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub compareSamples{
+
+    my ($dir,$list_samples,$sample_file,$reference_file,$min_overlap,$same_sv_type,$file_names)=@_;
+
+    my @bedpefiles;
+    my @list=split(",",$list_samples);
+    my @list_files=($sample_file,$reference_file);
+
+    print LOG "\# Comparison procedure...\n";
+    print LOG "-- samples=$list_samples\n".
+	 "-- minimum overlap=$min_overlap\n".
+	 "-- same SV type=$same_sv_type\n";
+
+    #conversion of links to bedPE format file
+    print LOG "-- Conversion of links.filtered files to bedPE format\n";
+    for my $s (0..$#list) {
+
+	links2bedPElinksfile($list[$s],$list_files[$s],$list_files[$s].".bedpe.txt");
+	push(@bedpefiles,$list_files[$s].".bedpe.txt");
+
+    }
+
+    #get common links between all samples compared
+    print LOG "-- Getting common links between all samples with BEDTools\n";
+    my $common_name=join(".",@list);
+
+    my $nb=scalar @list;
+    my $command="";
+    my $prompt=">";
+
+    while ($nb>0){
+
+	for my $i (0..$#list_files){
+
+	    $command.="$BEDTOOLS_BIN_DIR/pairToPair -type both -f $min_overlap -a ".$list_files[$i].".bedpe.txt";
+	    my $pipe=0;
+
+	    for my $j ($i+1..$#list_files){
+
+		$command.="| $BEDTOOLS_BIN_DIR/pairToPair -type both -f $min_overlap -a stdin" if($pipe);
+		$command.=" -b ".$list_files[$j].".bedpe.txt";
+		$pipe=1;
+
+	    }
+
+	    $command.=$prompt.$dir.$common_name.".bedpe.tmp;";
+	    $prompt=">>";
+
+	    my $first=shift(@list_files); push(@list_files,$first);
+	    last;
+	}
+	$nb--;
+    }
+
+    system ($command);
+
+    push(@bedpefiles,$dir.$common_name.".bedpe.tmp");
+
+    #Post comparison to get common links if same type only (as an option)
+    open( FILE, "<".$dir.$common_name.".bedpe.tmp") or die "Can't open".$dir.$common_name.".bedpe.tmp : $!";
+    open( OUT, ">".$dir.$common_name.".bedpe.unique") or die "Can't write in ".$dir.$common_name.".bedpe.unique : $!";
+
+    while(<FILE>){
+	my @t=split("\t",$_);
+	my $s=(split("_",$t[6]))[0];
+	my ($sv1,$sv2)=($t[7],$t[18]);
+	splice(@t,11,$#t);
+
+	if($same_sv_type){
+	    print OUT join("\t",@t)."\n" if($sv1 eq $sv2);
+	}else{
+	    print OUT join("\t",@t)."\n";
+	}
+    }
+    close FILE;
+    close OUT;
+
+    bedPElinks2linksfile($dir.$common_name.".bedpe.unique", $dir.$common_name.".unique");
+    push(@bedpefiles,$dir.$common_name.".bedpe.unique");
+    push(@$file_names,$common_name.".unique");
+    print LOG "-- output created: ".$dir.$common_name.".compared\n";
+
+    #get specific links for each sample
+    print LOG "-- Getting specific links for each sample\n";
+    for my $s (0..$#list) {
+	system("grep -Fxv -f ".$dir.$common_name.".bedpe.unique ".$list_files[$s].".bedpe.txt >".$dir.$list[$s].".bedpe.unique");
+	bedPElinks2linksfile($dir.$list[$s].".bedpe.unique",$dir.$list[$s].".unique");
+	push(@bedpefiles,$dir.$list[$s].".bedpe.unique");
+	push(@$file_names,$list[$s].".unique");
+	print LOG "-- output created: ".$dir.$list[$s].".compared\n";
+    }
+
+    unlink(@bedpefiles);
+
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+#convert the links file to the circos format
+sub links2segdup{
+
+    my($id,$color_code,$links_file,$segdup_file)=@_;
+
+    print LOG "\# Converting to the circos format...\n";
+
+    tie (my %hcolor,'Tie::IxHash');						#color-code hash table
+    foreach my $col (keys %{$color_code}){
+	my ($min_links,$max_links)=split(",",$color_code->{$col});
+	$hcolor{$col}=[$min_links,$max_links];
+    }
+
+    open LINKS, "<$links_file" or die "$0: can't open $links_file :$!\n";
+    open SEGDUP, ">$segdup_file" or die "$0: can't write in the output: $segdup_file :$!\n";
+
+    my $index=1;
+    while(<LINKS>){
+
+	my ($chr1,$start1,$end1,$chr2,$start2,$end2,$count)=(split)[0,1,2,3,4,5,6];
+
+	my $color=getColor($count,\%hcolor,"circos");				#get the color-code according the number of links
+
+	print SEGDUP "$index\t$id$chr1\t$start1\t$end1\tcolor=$color\n".	#circos output
+		     "$index\t$id$chr2\t$start2\t$end2\tcolor=$color\n";
+	$index++;
+    }
+
+    close LINKS;
+    close SEGDUP;
+    print LOG "-- output created: $segdup_file\n";
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+#convert the links file to the bedPE format for BEDTools usage
+sub links2bedPElinksfile{
+
+    my ($sample,$links_file,$bedpe_file)=@_;
+
+    open LINKS, "<$links_file" or die "$0: can't open $links_file :$!\n";
+    open BEDPE, ">$bedpe_file" or die "$0: can't write in the output: $bedpe_file :$!\n";
+
+    my $nb_links=1;
+
+    while(<LINKS>){
+
+	chomp;
+	my @t=split("\t",$_);
+	my ($chr1,$start1,$end1,$chr2,$start2,$end2)=splice(@t,0,6);
+	my $type=($chr1 eq $chr2)? "INTRA":"INTER";
+	$type.="_".$t[10];
+
+	$start1--; $start2--;
+
+	print BEDPE "$chr1\t$start1\t$end1\t$chr2\t$start2\t$end2".
+	"\t$sample"."_link$nb_links\t$type\t.\t.".
+	"\t".join("|",@t)."\n";
+
+	$nb_links++;
+    }
+
+    close LINKS;
+    close BEDPE;
+
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub bedPElinks2linksfile{
+
+    my ($bedpe_file,$links_file)=@_;
+
+    open BEDPE, "<$bedpe_file" or die "$0: can't open: $bedpe_file :$!\n";
+    open LINKS, ">$links_file" or die "$0: can't write in the output $links_file :$!\n";
+
+    while(<BEDPE>){
+
+	chomp;
+	my $sample=(split("_",(split("\t",$_))[6]))[0];
+	my @t1=(split("\t",$_))[0,1,2,3,4,5];
+	my @t2=split(/\|/,(split("\t",$_))[10]);
+	push(@t2,$sample);
+
+	print LINKS join("\t",@t1)."\t".join("\t",@t2)."\n";
+
+    }
+    close BEDPE;
+    close LINKS;
+
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+#convert the links file to the bed format
+sub links2bedfile{
+
+    my ($tag_length,$color_code,$links_file,$bed_file)=@_;
+
+    print LOG "\# Converting to the bed format...\n";
+
+    my $compare=1;
+    if($links_file!~/compared$/){
+	$compare=0;
+	$tag_length->{none}->{1}=$tag_length->{1};
+	$tag_length->{none}->{2}=$tag_length->{2};
+    }
+
+    #color-code hash table
+    tie (my %hcolor,'Tie::IxHash');
+    my %color_order;
+    $color_order{"255,255,255"}=0;
+    my $n=1;
+    foreach my $col (keys %{$color_code}){
+	my ($min_links,$max_links)=split(",",$color_code->{$col});
+	$hcolor{$col}=[$min_links,$max_links];
+	$color_order{$col}=$n;
+	$n++;
+    }
+
+    my %pair;
+    my %pt;
+    $n=1;
+    open LINKS, "<$links_file" or die "$0: can't open $links_file:$!\n";
+
+    my %str=( "F"=>"+", "R"=>"-" );
+
+    while(<LINKS>){
+
+	my @t=split;
+	my $sample=($compare)? pop(@t):"none";
+
+	my $chr1=$t[0];
+	my $chr2=$t[3];
+	$chr1 = "chr".$chr1 unless ($chr1 =~ m/chr/i);
+	$chr2 = "chr".$chr2 unless ($chr2 =~ m/chr/i);
+	my $same_chr=($chr1 eq $chr2)? 1:0;
+
+	my $count=$t[6];
+	my $color=getColor($count,\%hcolor,"bed");
+
+	my @pairs=deleteBadOrderSensePairs(split(",",$t[7]));
+	my @strand1=deleteBadOrderSensePairs(split(",",$t[8]));
+	my @strand2=deleteBadOrderSensePairs(split(",",$t[9]));
+	my @ends_order1=deleteBadOrderSensePairs(split(",",$t[10]));
+	my @ends_order2=deleteBadOrderSensePairs(split(",",$t[11]));
+	my @position1=deleteBadOrderSensePairs(split(",",$t[14]));
+	my @position2=deleteBadOrderSensePairs(split(",",$t[15]));
+	my @start1; my @end1; getCoordswithLeftMost(\@start1,\@end1,\@position1,\@strand1,\@ends_order1,$tag_length->{$sample});
+	my @start2; my @end2; getCoordswithLeftMost(\@start2,\@end2,\@position2,\@strand2,\@ends_order2,$tag_length->{$sample});
+
+
+	for my $p (0..$#pairs){
+
+	    if (!exists $pair{$pairs[$p]}){
+
+		if($same_chr){
+
+		    $pair{$pairs[$p]}->{0}=[ $chr1, $start1[$p]-1, $end2[$p], $pairs[$p], 0, $str{$strand1[$p]},
+				    $start1[$p]-1, $end2[$p], $color,
+				    2, $tag_length->{$sample}->{$ends_order1[$p]}.",".$tag_length->{$sample}->{$ends_order2[$p]}, "0,".($start2[$p]-$start1[$p]) ];
+		    $pt{$n}=$pair{$pairs[$p]}->{0};
+		    $n++;
+
+		}else{
+
+		    $pair{$pairs[$p]}->{1}=[ $chr1, $start1[$p]-1, $end1[$p] , $pairs[$p]."/1", 0, $str{$strand1[$p]},
+						$start1[$p]-1, $end1[$p], $color,
+						1, $tag_length->{$sample}->{$ends_order1[$p]}, 0];
+		    $pt{$n}=$pair{$pairs[$p]}->{1};
+		    $n++;
+
+
+		    $pair{$pairs[$p]}->{2}=[ $chr2, $start2[$p]-1, $end2[$p], $pairs[$p]."/2", 0, $str{$strand2[$p]},
+						$start2[$p]-1, $end2[$p], $color,
+						1, $tag_length->{$sample}->{$ends_order2[$p]}, 0];
+		    $pt{$n}=$pair{$pairs[$p]}->{2};
+		    $n++;
+		}
+	    }else{
+
+		if($same_chr){
+		    ${$pair{$pairs[$p]}->{0}}[8]=$color if($color_order{$color}>$color_order{${$pair{$pairs[$p]}->{0}}[8]});
+		}else{
+		    ${$pair{$pairs[$p]}->{1}}[8]=$color if($color_order{$color}>$color_order{${$pair{$pairs[$p]}->{1}}[8]});
+		    ${$pair{$pairs[$p]}->{2}}[8]=$color if($color_order{$color}>$color_order{${$pair{$pairs[$p]}->{2}}[8]});
+		}
+	    }
+	}
+    }
+    close LINKS;
+
+    my $nb_pairs=$n-1;
+
+    open BED, ">$bed_file" or die "$0: can't write in the output: $bed_file :$!\n";
+    print BED "track name=\"$bed_file\" description=\"mate pairs involved in links\" ".
+	      "visibility=2 itemRgb=\"On\"\n";
+
+    for my $i (1..$nb_pairs){
+	print BED join("\t",@{$pt{$i}})."\n";
+    }
+
+    close BED;
+
+    print LOG "-- output created: $bed_file\n";
+
+    undef %pair;
+    undef %pt;
+
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub links2SVfile{
+
+    my($links_file,$sv_file)=@_;
+
+    print LOG "\# Converting to the sv output table...\n";
+    open LINKS, "<$links_file" or die "$0: can't open $links_file:$!\n";
+    open SV, ">$sv_file" or die "$0: can't write in the output: $sv_file :$!\n";
+
+    my @header=qw(chr_type SV_type BAL_type chromosome1 start1-end1 average_dist
+    chromosome2 start2-end2 nb_pairs score_strand_filtering score_order_filtering score_insert_size_filtering
+    final_score breakpoint1_start1-end1 breakpoint2_start2-end2);
+
+    my $nb_links=0;
+
+    while (<LINKS>){
+
+	my @t=split;
+	my @sv=();
+	my $sv_type="-";
+	my $strand_ratio="-";
+	my $eq_ratio="-";
+	my $eq_type="-";
+	my $insert_ratio="-";
+	my $link="-";
+	my ($bk1, $bk2)=("-","-");
+	my $score="-";
+
+	my ($chr1,$start1,$end1)=($t[0],$t[1],$t[2]);
+	my ($chr2,$start2,$end2)=($t[3],$t[4],$t[5]);
+	my $nb_pairs=$t[6];
+	$chr1 = "chr".$chr1 unless ($chr1 =~ m/chr/i);
+	$chr2 = "chr".$chr2 unless ($chr2 =~ m/chr/i);
+	my $chr_type=($chr1 eq $chr2)? "INTRA":"INTER";
+
+	#if strand filtering
+	if (defined $t[16]){
+	    #if inter-chr link
+	    $sv_type=$t[16];
+	    if(defined $t[17] && $t[17]=~/^(\d+)\/(\d+)$/){
+		$strand_ratio=floor($1/$2*100)."%";
+		$score=$t[18];
+	    }
+	    if(defined $t[18] && $t[18]=~/^(\d+)\/(\d+)$/){
+	    #if intra-chr link with insert size filtering
+		$strand_ratio=floor($1/$2*100)."%";
+		$link=floor($t[17]);
+		if($sv_type!~/^INV/){
+		    $insert_ratio=floor($1/$2*100)."%" if($t[19]=~/^(\d+)\/(\d+)$/);
+		    $score=$t[20];
+		}else{
+		    $score=$t[19];
+		}
+	    }
+	}
+
+	if(defined $t[18] && ($t[18] eq "UNBAL" || $t[18] eq "BAL")){
+
+	    #if strand and order filtering only and/or interchr link
+	    $eq_type=$t[18];
+	    $eq_ratio=floor($1/$2*100)."%" if($t[19]=~/^(\d+)\/(\d+)$/);
+	    ($bk1, $bk2)=($t[20],$t[21]);
+	    foreach my $bk ($bk1, $bk2){$bk=~s/\),\(/ /g; $bk=~s/(\(|\))//g; $bk=~s/,/-/g;}
+	    $score=$t[22];
+
+	}elsif(defined $t[19] && ($t[19] eq "UNBAL" || $t[19] eq "BAL")){
+
+	    #if all three filtering
+	    $link=floor($t[17]);
+	    $eq_type=$t[19];
+	    $eq_ratio=floor($1/$2*100)."%" if($t[20]=~/^(\d+)\/(\d+)$/);
+
+	    if(defined $t[21] && $t[21]=~/^(\d+)\/(\d+)$/){
+		$insert_ratio=floor($1/$2*100)."%";
+		($bk1, $bk2)=($t[22],$t[23]);
+		$score=$t[24];
+
+	    }else{
+		($bk1, $bk2)=($t[21],$t[22]);
+		$score=$t[23];
+	    }
+	    foreach my $bk ($bk1, $bk2){$bk=~s/\),\(/ /g; $bk=~s/(\(|\))//g; $bk=~s/,/-/g;}
+
+	}
+
+
+	push(@sv, $chr_type, $sv_type,$eq_type);
+	push(@sv,"$chr1\t$start1-$end1");
+	push(@sv, $link);
+	push(@sv,"$chr2\t$start2-$end2",
+	     $nb_pairs,$strand_ratio,$eq_ratio,$insert_ratio, decimal($score,4), $bk1, $bk2);
+
+
+	print SV join("\t",@sv)."\n";
+    }
+
+    close LINKS;
+    close SV;
+
+    system "sort  -k 9,9nr -k 13,13nr $sv_file > $sv_file.sorted";
+
+    open SV, "<".$sv_file.".sorted" or die "$0: can't open in the output: $sv_file".".sorted :$!\n";
+    my @links=<SV>;
+    close SV;
+
+    open SV, ">$sv_file" or die "$0: can't write in the output: $sv_file :$!\n";
+
+    print SV join("\t",@header)."\n";
+    print SV @links;
+    close SV;
+
+    unlink($sv_file.".sorted");
+
+    print LOG "-- output created: $sv_file\n";
+
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub deleteBadOrderSensePairs{
+
+    my (@tab)=@_;
+    my @tab2;
+
+    foreach my $v (@tab){
+
+	$v=~s/[\(\)]//g;
+	push(@tab2,$v) if($v!~/[\$\*\@]$/);
+    }
+    return @tab2;
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+#gets starts and ends Coords when start=leftmost given positions, directions and orders
+sub getCoordswithLeftMost {
+
+    my ($starts,$ends,$positions,$strand,$end_order,$tag_length) = @_;
+
+    for my $i (0..scalar(@{$positions})-1) {
+
+	if($strand->[$i] eq 'F'){
+	    $starts->[$i]=$positions->[$i];
+	    $ends->[$i]=$positions->[$i]+$tag_length->{$end_order->[$i]}-1;
+	}else{
+	    $starts->[$i]=$positions->[$i]-$tag_length->{$end_order->[$i]}+1;
+	    $ends->[$i]=$positions->[$i];
+	}
+    }
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub floor {
+    my $nb = $_[0];
+    $nb=~ s/\..*//;
+    return $nb;
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub decimal{
+
+  my $num=shift;
+  my $digs_to_cut=shift;
+
+  $num=sprintf("%.".($digs_to_cut-1)."f", $num) if ($num=~/\d+\.(\d){$digs_to_cut,}/);
+
+  return $num;
+}
+
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+#Sort links according the concerned chromosomes and their coordinates
+sub sortLinks{
+
+    my ($links_file,$sortedlinks_file,$unique)=@_;
+
+    print LOG "# Sorting links...\n";
+
+    my $pipe=($unique)? "| sort -u":"";
+    system "sort -k 1,1 -k 4,4 -k 2,2n -k 5,5n -k 8,8n $links_file $pipe > $sortedlinks_file";
+
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub getColor{
+
+    my($count,$hcolor,$format)=@_;
+    for my $col ( keys % { $hcolor} ) {
+       return $col if($count>=$hcolor->{$col}->[0] && $count<=$hcolor->{$col}->[1]);
+    }
+    return "white" if($format eq "circos");
+    return "255,255,255" if($format eq "bed");
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+#check if the configuration file is correct
+sub validateconfiguration{
+
+    my %conf=%{$_[0]};
+    my $list_prgs="@ARGV";
+
+    my @circos_params=qw(organism_id colorcode);
+    my @bed_params=qw(colorcode);
+    my @compare_params=qw(list_samples list_read_lengths sample_link_file reference_link_file);
+
+    unless (defined($conf{general}{output_dir})) {
+	$conf{general}{output_dir} = ".";
+    }
+    unless (-d $conf{general}{output_dir}){
+	mkdir $conf{general}{output_dir} or die;
+    }
+    $conf{general}{output_dir}.="/" if($conf{general}{output_dir}!~/\/$/);
+
+
+    if($list_prgs=~/links2compare/){
+	foreach my $p (@compare_params) {
+	    die("Error Config : The compare parameter \"$p\" is not defined\n") if (!defined $conf{compare}{$p});
+	}
+
+	unless (defined($conf{compare}{same_sv_type})) {
+	    $conf{compare}{same_sv_type} = 0;
+	}
+
+	unless (defined($conf{compare}{min_overlap})) {
+	    $conf{compare}{min_overlap} = 1E-9;
+	}
+
+	if($conf{compare}{circos_output}){
+	    foreach my $p (@circos_params) {
+		next if($list_prgs=~/^ratio/ && $p eq "colorcode");
+		die("Error Config : The circos parameter \"$p\" is not defined\n") if (!defined $conf{circos}{$p});
+	    }
+	}
+	if($conf{compare}{bed_output}){
+	    foreach my $p (@bed_params) {
+		die("Error Config : The bed parameter \"$p\" is not defined\n") if (!defined $conf{bed}{$p});
+	    }
+	    die("Error Config : The compare parameter \"list_read_lengths\" is not defined\n") if (!defined $conf{compare}{list_read_lengths});
+
+	    my @samples=split(",",$conf{compare}{list_samples});
+	    my @read_lengths=split(",",$conf{compare}{list_read_lengths});
+	    for my $i (0..$#samples){
+		my @l=split("-",$read_lengths[$i]);
+		$conf{compare}{read_lengths}{$samples[$i]}={ 1=> $l[0], 2=> $l[1]};
+	    }
+	}
+    }
+
+
+}
+#------------------------------------------------------------------------------#
+#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::#
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/svdetect/SVDetect_compare.xml	Thu Jul 12 12:10:27 2012 -0400
@@ -0,0 +1,218 @@
+<tool id="svdetect_compare" name="Compare">
+
+<description>structural variants between two samples</description>
+
+<command interpreter="perl">SVDetect_compare.pl links2compare -conf '$config_file' -l '$log_file' -N '$sample_name.$reference_name'
+
+#if $links2SV
+-out1 '$common_sv_file'
+-out2 '$sample_sv_file'
+-out3 '$reference_sv_file'
+#end if
+
+#if $file_conversion.file_conversion_select=="convert" and $file_conversion.links2circos
+-out4 '$common_circos_file'
+-out5 '$sample_circos_file'
+-out6 '$reference_circos_file'
+#end if
+
+#if $file_conversion.file_conversion_select=="convert" and $file_conversion.links2bed
+-out7 '$common_bed_file'
+-out8 '$sample_bed_file'
+-out9 '$reference_bed_file'
+#end if
+
+</command>
+
+<inputs>
+	<param name="sample_name" type="text" size="20" value="sample" label="Sample Name"/>
+	<param name="sample_read1_length" type="integer" size="10" value="50" label="Sample read 1 length (bp)"/>
+	<param name="sample_read2_length" type="integer" size="10" value="50" label="Sample read 2 length (bp)"/>
+   	<param name="sample_mates_file" type="data" format="sv" label="Sample input file" help=".sv file"/>
+
+	<param name="reference_name" type="text" size="20" value="reference" label="Reference Name"/>
+	<param name="reference_read1_length" type="integer" size="10" value="50" label="Reference read 1 length (bp)"/>
+	<param name="reference_read2_length" type="integer" size="10" value="50" label="Reference read 2 length (bp)"/>
+   	<param name="reference_mates_file" type="data" format="sv" label="Reference input file" help=".sv file"/>
+
+	<param name="min_overlap" type="float" size="10"  value="0.05" label="Minimum overlap of links required as a fraction"/>
+	<param name="same_sv_type" label="Comparison of SVs with the same type only ?" type="boolean" truevalue="1" falsevalue="0" checked="True"/>
+
+	<param name="links2SV" label="Do you want to have filtered links in a tabulated file format showing significant SVs?" type="boolean" truevalue="1" falsevalue="0" checked="True"/>
+
+	<conditional name="file_conversion">
+		<param name="file_conversion_select" type="select" label="Output file conversion" help="Converts filtered links to Circos/BED files format for graphical view of SVs">
+			<option value="do_not_convert">No</option>
+			<option value="convert">Yes</option>
+		</param>
+		<when value="do_not_convert">
+				<!-- do nothing here -->
+		</when>
+		<when value="convert">
+			<param name="links2circos" label="Converts the link list to the Circos link format" type="boolean" truevalue="1" falsevalue="0" checked="True"/>
+			<param name="links2bed" label="Converts the link list to the UCSC BED format" type="boolean" truevalue="1" falsevalue="0" checked="False"/>
+			<param name="organism_id" type="text" size="10" value="hs" label="Organism ID"/>
+			<repeat name="color_code" title="Color-code" min="1" max="7">
+				<param name="color" type="select" label="Color">
+					<option value="grey">grey</option>
+					<option value="black">black</option>
+					<option value="blue">blue</option>
+					<option value="green">green</option>
+					<option value="purple">purple</option>
+					<option value="orange">orange</option>
+					<option value="red">red</option>
+				</param>
+				<param name="interval" type="text" value="1,3" label="Interval"/>
+			</repeat>
+		</when>
+	</conditional>
+</inputs>
+
+
+
+<outputs>
+	<data format="sv" name="common_sv_file" label="common.compared.sv">
+		<filter>links2SV is True</filter>
+	</data>
+	<data format="sv" name="sample_sv_file" label="${sample_name}.compared.sv">
+		<filter>links2SV is True</filter>
+	</data>
+	<data format="sv" name="reference_sv_file" label="${reference_name}.compared.sv">
+		<filter>links2SV is True</filter>
+	</data>
+
+	<data format="segdup" name="common_circos_file" label="common.compared.segdup">
+		<filter>(
+			file_conversion['file_conversion_select']=="convert" and
+			file_conversion['links2circos'] is True
+			)
+		</filter>
+	</data>
+	<data format="segdup" name="sample_circos_file" label="${sample_name}.compared.segdup">
+		<filter>(
+			file_conversion['file_conversion_select']=="convert" and
+			file_conversion['links2circos'] is True
+			)
+		</filter>
+	</data>
+	<data format="segdup" name="reference_circos_file" label="${reference_name}.compared.segdup">
+		<filter>(
+			file_conversion['file_conversion_select']=="convert" and
+			file_conversion['links2circos'] is True
+			)
+		</filter>
+	</data>
+
+	<data format="bed" name="common_bed_file" label="common.compared.bed">
+		<filter>(
+			file_conversion['file_conversion_select']=="convert" and
+			file_conversion['links2bed'] is True
+			)
+		</filter>
+	</data>
+	<data format="bed" name="sample_bed_file" label="${sample_name}.compared.bed">
+		<filter>(
+			file_conversion['file_conversion_select']=="convert" and
+			file_conversion['links2bed'] is True
+			)
+		</filter>
+	</data>
+	<data format="bed" name="reference_bed_file" label="${reference_name}.compared.bed">
+		<filter>(
+			file_conversion['file_conversion_select']=="convert" and
+			file_conversion['links2bed'] is True
+			)
+		</filter>
+	</data>
+
+	<data format="txt" name="log_file" label="${sample_name}.${reference_name}.svdetect_compare.log"/>
+</outputs>
+
+
+
+<configfiles>
+	<configfile name="config_file">
+&lt;general&gt;
+output_dir=$__new_file_path__/svdetect
+&lt;/general&gt;
+
+#if $file_conversion.file_conversion_select == "convert"
+#if $file_conversion.links2circos
+&lt;circos&gt;
+organism_id=${file_conversion.organism_id}
+&lt;colorcode&gt;
+#for $color_repeat in $file_conversion.color_code
+${color_repeat.color}=${color_repeat.interval}
+#end for
+&lt;/colorcode&gt;
+&lt;/circos&gt;
+#end if
+#if $file_conversion.links2bed
+&lt;bed&gt;
+&lt;colorcode&gt;
+#for $color_repeat in $file_conversion.color_code
+#if str($color_repeat.color)== "grey"
+190,190,190=${color_repeat.interval}
+#end if
+#if str($color_repeat.color)== "black"
+0,0,0=${color_repeat.interval}
+#end if
+#if str($color_repeat.color)== "blue"
+0,0,255=${color_repeat.interval}
+#end if
+#if str($color_repeat.color)== "green"
+0,255,0=${color_repeat.interval}
+#end if
+#if str($color_repeat.color)== "purple"
+153,50,205=${color_repeat.interval}
+#end if
+#if str($color_repeat.color)== "orange"
+255,140,0=${color_repeat.interval}
+#end if
+#if str($color_repeat.color)== "red"
+255,0,0=${color_repeat.interval}
+#end if
+#end for
+&lt;/colorcode&gt;
+&lt;/bed&gt;
+#end if
+#end if
+
+&lt;compare&gt;
+list_samples=${sample_name},${reference_name}
+list_read_lengths=${sample_read1_length}-${sample_read2_length},${reference_read1_length}-${reference_read2_length}
+sample_link_file=${sample_mates_file}
+reference_link_file=${reference_mates_file}
+min_overlap=${min_overlap}
+same_sv_type=${same_sv_type}
+sv_output=${links2SV}
+#if $file_conversion.file_conversion_select == "convert"
+circos_output=${$file_conversion.links2circos}
+bed_output=${$file_conversion.links2bed}
+#end if
+&lt;/compare&gt;
+
+	</configfile>
+</configfiles>
+
+  <help>
+**What it does**
+
+SVDetect - Version : 0.8
+
+Comparison of clusters between two samples to get common or sample-specific SVs
+
+This program is designed to compare filtered links between two anomalously mapped mate-pair/paired-end datasets
+and to identify common and sample-specific SVs (like the usual sample/reference design).
+Overlaps between coordinates of clusters and types of SVs are used as parameters of comparison.
+
+Manual documentation available at the http://svdetect.sourceforge.net/Site/Manual.html
+
+-----
+
+.. class:: infomark
+
+Contact Bruno Zeitouni (svdetect@curie.fr) for any questions or concerns about the Galaxy implementation of SVDetect.
+  </help>
+
+</tool>
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/svdetect/SVDetect_import.sh	Thu Jul 12 12:10:27 2012 -0400
@@ -0,0 +1,15 @@
+#!/bin/bash
+
+
+while getopts "i:o:" optionName; do
+case "$optionName" in
+
+i) INPUT="$OPTARG";;
+o) OUTPUT="$OPTARG";;
+
+esac
+done
+
+rm $OUTPUT
+
+ln -s $INPUT $OUTPUT
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/svdetect/SVDetect_import.xml	Thu Jul 12 12:10:27 2012 -0400
@@ -0,0 +1,85 @@
+<tool id="svdetect_import" name="Import data">
+  <description>BAM, chromosome info or sv files</description>
+  <command interpreter="bash">SVDetect_import.sh -i $file_path
+  #if str($type.file_type)=="bam"
+  -o $outbamfile
+  #elif str($type.file_type)=="len"
+  -o $outlenfile
+  #elif str($type.file_type)=="sv"
+  -o $outsvfile
+  #end if
+  </command>
+  <inputs>
+    <param name="file_name" type="text" value="file1" label="File Name"/>
+    <conditional name="type">
+      <param name="file_type" type="select" label="Select the file type to import" help="BAM file (BAM) or text file (SAM, chromosome list or a SV tabulated text file)">
+          <option value="bam">BAM file (.bam)</option>
+          <option value="len">Chromosome info file (.len)</option>
+          <option value="sv">SVDetect output file (.sv)</option>
+      </param>
+    <when value="bam">
+      <!-- do nothing here -->
+    </when>
+    <when value="len">
+      <!-- do nothing here -->
+    </when>
+     <when value="sv">
+      <!-- do nothing here -->
+    </when>
+    </conditional>
+     <param name="file_path" type="text" size="150" label="Path to file"/>
+  </inputs>
+  <outputs>
+    <data format="bam" name="outbamfile" label="${file_name}.bam">
+      <filter>type['file_type']=="bam"</filter>
+    </data>
+    <data format="len" name="outlenfile" label="${file_name}.len">
+      <filter>type['file_type']=="len"</filter>
+    </data>
+    <data format="sv" name="outsvfile" label="${file_name}.sv">
+      <filter>type['file_type']=="sv"</filter>
+    </data>
+  </outputs>
+  <help>
+**What it does**
+
+This tool allows you to import quickly a BAM file, a chromosome info file or a SVDetect output file from you computer as inputs for SVDetect.
+
+
+**Example of chromosome file**
+
+Input len file::
+
+    1  chr1  247249719
+    2  chr2  242951149
+    3  chr3  199501827
+    4  chr4  191273063
+    5  chr5  180857866
+    6  chr6  170899992
+    7  chr7  158821424
+    8  chr8  146274826
+    9  chr9  140273252
+    10  chr10  135374737
+    11  chr11  134452384
+    12  chr12  132349534
+    13  chr13  114142980
+    14  chr14  106368585
+    15  chr15  100338915
+    16  chr16  88827254
+    17  chr17  78774742
+    18  chr18  76117153
+    19  chr19  63811651
+    20  chr20  62435964
+    21  chr21  46944323
+    22  chr22  49691432
+    23  chrX  154913754
+    24  chrY  57772954
+
+-----
+
+.. class:: infomark
+
+Contact Bruno Zeitouni (svdetect@curie.fr) for any questions or concerns about the Galaxy implementation of SVDetect.
+  </help>
+
+</tool>
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/svdetect/SVDetect_run_parallel.pl	Thu Jul 12 12:10:27 2012 -0400
@@ -0,0 +1,3537 @@
+#!/usr/bin/perl -w
+
+=pod
+
+=head1 NAME
+
+SVDetect - Program designed to the detection of structural variations
+from paired-end/mate-pair sequencing data, compatible with SOLiD and Illumina (>=1.3) reads
+
+Version: 0.8 for Galaxy
+
+=head1 SYNOPSIS
+
+SVDetect <command> -conf <configuration_file> [-help] [-man]
+
+    Command:
+
+    	linking		detection and isolation of links
+        filtering	filtering of links according different parameters
+        links2circos	links conversion to circos format
+	links2bed 	paired-ends of links converted to bed format (UCSC)
+	links2SV	formatted output to show most significant SVs
+	cnv		calculate copy-number profiles
+	ratio2circos	ratio conversion to circos density format
+	ratio2bedgraph	ratio conversion to bedGraph density format (UCSC)
+
+=head1 DESCRIPTION
+
+This is a command-line interface to SVDetect.
+
+
+=head1 AUTHORS
+
+Bruno Zeitouni E<lt>bruno.zeitouni@curie.frE<gt>,
+Valentina Boeva E<lt>valentina.boeva@curie.frE<gt>
+
+=cut
+
+# -------------------------------------------------------------------
+
+use strict;
+use warnings;
+
+use Pod::Usage;
+use Getopt::Long;
+
+use Config::General;
+use Tie::IxHash;
+use FileHandle;
+use Parallel::ForkManager;
+
+#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::#
+#PARSE THE COMMAND LINE
+my %OPT;
+GetOptions(\%OPT,
+	   'conf=s',
+	   'out1=s', #GALAXY
+	   'out2=s', #GALAXY
+	   'out3=s', #GALAXY
+	   'out4=s', #GALAXY
+	   'out5=s', #GALAXY
+	   'l=s', #GALAXY
+	   'N=s',#GALAXY
+	   'help',#GALAXY
+           'man'
+	  );
+
+pod2usage() if $OPT{help};
+pod2usage(-verbose=>2) if $OPT{man};
+pod2usage(-message=> "$!", -exitval => 2) if (!defined $OPT{conf});
+
+pod2usage() if(@ARGV<1);
+
+tie (my %func, 'Tie::IxHash',linking=>\&createlinks,
+			     filtering=>\&filterlinks,
+			     links2circos=>\&links2circos,
+			     links2bed=>\&links2bed,
+			     links2compare=>\&links2compare,
+			     links2SV=>\&links2SV,
+			     cnv=>\&cnv,
+			     ratio2circos=>\&ratio2circos,
+			     ratio2bedgraph=>\&ratio2bedgraph);
+
+foreach my $command (@ARGV){
+    pod2usage(-message=> "Unknown command \"$command\"", -exitval => 2) if (!defined($func{$command}));
+}
+#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::#
+
+
+#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::#
+#READ THE CONFIGURATION FILE
+my $conf=Config::General->new(    -ConfigFile        => $OPT{conf},
+                                  -Tie => "Tie::IxHash",
+                                  -AllowMultiOptions => 1,
+				  -LowerCaseNames    => 1,
+				  -AutoTrue => 1);
+my %CONF= $conf->getall;
+validateconfiguration(\%CONF);							#validation of the configuration parameters
+
+my $SAMTOOLS_BIN_DIR="/bioinfo/local/samtools"; #GALAXY
+
+my $pt_log_file=$OPT{l}; #GALAXY
+my $pt_links_file=$OPT{out1} if($OPT{out1}); #GALAXY
+my $pt_flinks_file=$OPT{out2} if($OPT{out2}); #GALAXY
+my $pt_sv_file=$OPT{out3} if($OPT{out3}); #GALAXY
+my $pt_circos_file=$OPT{out4} if($OPT{out4}); #GALAXY
+my $pt_bed_file=$OPT{out5} if($OPT{out5}); #GALAXY
+
+$CONF{general}{mates_file}=readlink($CONF{general}{mates_file});#GALAXY
+$CONF{general}{cmap_file}=readlink($CONF{general}{cmap_file});#GALAXY
+
+my $log_file=$CONF{general}{output_dir}.$OPT{N}.".svdetect_run.log"; #GALAXY
+open LOG,">$log_file" or die "$0: can't open ".$log_file.":$!\n";#GALAXY
+#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::#
+
+#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::#
+#COMMAND EXECUTION
+foreach my $command (@ARGV){
+    &{$func{$command}}();
+}
+print LOG "-- end\n";#GALAXY
+
+close LOG;#GALAXY
+system "rm $pt_log_file ; ln -s $log_file $pt_log_file"; #GALAXY
+exit(0);
+#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::#
+
+
+#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::#
+#FUNCTIONS
+#------------------------------------------------------------------------------#
+#MAIN FUNCTION number 1: Detection of links from mate-pairs data
+sub createlinks{
+
+    my %CHR;									#main hash table 1: fragments, links
+    my %CHRID;
+    my @MATEFILES;
+
+    my $output_prefix=$CONF{general}{mates_file}.".".$CONF{general}{sv_type};
+    my @path=split(/\//,$output_prefix);
+    $output_prefix=$CONF{general}{output_dir}.$path[$#path];
+    my $tmp_mates_prefix=$CONF{general}{tmp_dir}."mates/".$path[$#path];
+    my $tmp_links_prefix=$CONF{general}{tmp_dir}."links/".$path[$#path];
+
+    shearingChromosome(\%CHR, \%CHRID,						#making the genomic fragment library with the detection parameters
+		     $CONF{detection}{window_size},
+		     $CONF{detection}{step_length},
+		     $CONF{general}{cmap_file});
+
+    if($CONF{detection}{split_mate_file}){
+
+	splitMateFile(\%CHR, \%CHRID, \@MATEFILES, $tmp_mates_prefix,
+		      $CONF{general}{sv_type},
+		      $CONF{general}{mates_file},
+		      $CONF{general}{input_format},
+		      $CONF{general}{read_lengths}
+		      );
+    }else{
+
+	@MATEFILES=qx{ls $tmp_mates_prefix*} or die "# Error: No splitted mate files already created at $CONF{general}{tmp_dir} :$!";
+	chomp(@MATEFILES);
+	print LOG "# Splitted mate files already created.\n";
+    }
+
+
+    #Parallelization of the linking per chromosome for intra +  interchrs
+    my $pm = new Parallel::ForkManager($CONF{general}{num_threads});
+
+    foreach my $matefile (@MATEFILES){
+
+	my $pid = $pm->start and next;
+	getlinks(\%CHR, \%CHRID, $matefile);
+	$pm->finish;
+
+    }
+    $pm->wait_all_children;
+
+    #Merge the chromosome links file into only one
+    my @LINKFILES= qx{ls $tmp_links_prefix*links} or die "# Error: No links files created at $CONF{general}{tmp_dir} :$!";
+    chomp(@LINKFILES);
+    catFiles( \@LINKFILES => "$output_prefix.links" );
+
+    system "rm $pt_links_file; ln -s $output_prefix.links $pt_links_file" if (defined $pt_links_file); #GALAXY
+    print LOG "# Linking end procedure : output created: $output_prefix.links\n";
+    #unlink(@LINKFILES);
+    #unlink(@MATEFILES);
+
+    undef %CHR;
+    undef %CHRID;
+
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub getlinks {
+
+    my ($chr,$chrID,$tmp_mates_prefix)=@_;
+
+    my $tmp_links_prefix=$tmp_mates_prefix;
+    $tmp_links_prefix=~s/\/mates\//\/links\//;
+
+    my %PAIR;									#main hash table 2: pairs
+
+    linking($chr,$chrID, \%PAIR,						#creation of all links from chromosome coordinates of pairs
+	  $CONF{general}{read_lengths},
+	  $CONF{detection}{window_size},
+	  $CONF{detection}{step_length},
+	  $tmp_mates_prefix,
+	  $CONF{general}{input_format},
+	  $CONF{general}{sv_type},
+	  "$tmp_links_prefix.links.mapped"
+	  );
+
+    getUniqueLinks("$tmp_links_prefix.links.mapped",					#remove the doublons
+		   "$tmp_links_prefix.links.unique");
+
+    defineCoordsLinks($chr,$chrID, \%PAIR,					#definition of the precise coordinates of links
+		      $CONF{general}{input_format},
+		      $CONF{general}{sv_type},
+		      $CONF{general}{read_lengths},
+		      "$tmp_links_prefix.links.unique",
+		      "$tmp_links_prefix.links.unique_defined");
+
+    sortLinks("$tmp_links_prefix.links.unique_defined",				#sorting links from coordinates
+	      "$tmp_links_prefix.links.sorted");
+
+    removeFullyOverlappedLinks("$tmp_links_prefix.links.sorted",			#remove redundant links
+			       "$tmp_links_prefix.links",1);				#file output
+
+
+    undef %PAIR;
+
+    unlink("$tmp_links_prefix.links.mapped",
+	   "$tmp_links_prefix.links.unique",
+	   "$tmp_links_prefix.links.unique_defined",
+	   "$tmp_links_prefix.links.sorted");
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub splitMateFile{
+
+    my ($chr,$chrID,$files_list,$output_prefix,$sv_type,$mates_file,$input_format,$tag_length)=@_;
+
+    print LOG "# Splitting the mate file \"$mates_file\" for parallel processing...\n";
+
+    my %filesHandle;
+
+    #fichier matefile inter
+    if($sv_type=~/^(all|inter)$/){
+	my $newFileName="$output_prefix.interchrs";
+	push(@{$files_list},$newFileName);
+	my $fh = new FileHandle;
+	$fh->open(">$newFileName");
+	$filesHandle{inter}=$fh;
+    }
+
+    #fichiers matefiles intra
+    if($sv_type=~/^(all|intra)$/){
+	foreach my $k (1..$chr->{nb_chrs}){
+	    my $newFileName=$output_prefix.".".$chr->{$k}->{name};
+	    push(@{$files_list},$newFileName);
+	    my $fh = new FileHandle;
+	    $fh->open(">$newFileName");
+	    $filesHandle{$k}=$fh;
+	}
+    }
+
+    if ($mates_file =~ /.gz$/) {
+	open(MATES, "gunzip -c $mates_file |") or die "$0: can't open ".$mates_file.":$!\n"; #gzcat
+    }elsif($mates_file =~ /.bam$/){
+	open(MATES, "$SAMTOOLS_BIN_DIR/samtools view $mates_file |") or die "$0: can't open ".$mates_file.":$!\n";#GALAXY
+    }else{
+	open MATES, "<".$mates_file or die "$0: can't open ".$mates_file.":$!\n";
+    }
+
+
+    while(<MATES>){
+
+	my @t=split;
+	my ($chr_read1, $chr_read2, $firstbase_read1, $firstbase_read2, $end_order_read1,$end_order_read2);
+
+	next if (!readMateFile(\$chr_read1, \$chr_read2, \$firstbase_read1, \$firstbase_read2, \$end_order_read1, \$end_order_read2, \@t, $input_format,$tag_length));
+
+	next unless (exists $chrID->{$chr_read1} && exists $chrID->{$chr_read2});
+
+	($chr_read1, $chr_read2)= ($chrID->{$chr_read1},$chrID->{$chr_read2});
+
+	if( ($sv_type=~/^(all|inter)$/) && ($chr_read1 ne $chr_read2) ){
+	    my $fh2print=$filesHandle{inter};
+	    print $fh2print join("\t",@t)."\n";
+	}
+
+	if( ($sv_type=~/^(all|intra)$/) && ($chr_read1 eq $chr_read2) ){
+	    my $fh2print=$filesHandle{$chr_read1};
+	    print $fh2print join("\t",@t)."\n";
+
+	}
+    }
+
+    foreach my $name (keys %filesHandle){
+	my $fh=$filesHandle{$name};
+	$fh->close;
+    }
+
+    print LOG "# Splitted mate files of \"$mates_file\" created.\n";
+}
+
+
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub filterlinks{
+
+    my %CHR;
+    my %CHRID;
+    my @LINKFILES;
+    my @FLINKFILES;
+
+    my $output_prefix=$CONF{general}{mates_file}.".".$CONF{general}{sv_type};
+    my @path=split(/\//,$output_prefix);
+    $output_prefix=$CONF{general}{output_dir}.$path[$#path];
+    my $tmp_links_prefix=$CONF{general}{tmp_dir}."links/".$path[$#path];
+
+    createChrHashTables(\%CHR,\%CHRID,
+			$CONF{general}{cmap_file});
+
+    if($CONF{filtering}{split_link_file}){
+
+	splitLinkFile(\%CHR, \%CHRID, \@LINKFILES,
+		      $tmp_links_prefix,
+		      $CONF{general}{sv_type},
+		      "$output_prefix.links",
+		      );
+    }else{
+
+	@LINKFILES=qx{ls $tmp_links_prefix*links} or die "# Error: No splitted link files already created\n";
+	chomp(@LINKFILES);
+	print LOG "# Splitted link files already created.\n";
+    }
+
+    #Parallelization of the filtering per chromosome for intra +  interchrs
+    my $pm = new Parallel::ForkManager($CONF{general}{num_threads});
+
+    foreach my $linkfile (@LINKFILES){
+
+	my $pid = $pm->start and next;
+	getFilteredlinks(\%CHR, \%CHRID, $linkfile);
+	$pm->finish;
+
+    }
+    $pm->wait_all_children;
+
+    #Merge the chromosome links file into only one
+    @FLINKFILES= qx{ls $tmp_links_prefix*filtered} or die "# Error: No links files created\n";
+    chomp(@FLINKFILES);
+    catFiles( \@FLINKFILES => "$output_prefix.links.filtered" );
+
+    system "rm $pt_flinks_file; ln -s $output_prefix.links.filtered $pt_flinks_file" if (defined $pt_flinks_file); #GALAXY
+    print LOG"# Filtering end procedure : output created: $output_prefix.links.filtered\n";
+
+    undef %CHR;
+    undef %CHRID;
+
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub splitLinkFile{
+
+    my ($chr,$chrID,$files_list,$input_prefix,$sv_type,$link_file)=@_;
+
+    print LOG "# Splitting the link file for parallel processing...\n";
+
+    my %filesHandle;
+
+    #fichier matefile inter
+    if($sv_type=~/^(all|inter)$/){
+	my $newFileName="$input_prefix.interchrs.links";
+	push(@{$files_list},$newFileName);
+	my $fh = new FileHandle;
+	$fh->open(">$newFileName");
+	$filesHandle{inter}=$fh;
+    }
+
+    #fichiers matefiles intra
+    if($sv_type=~/^(all|intra)$/){
+	foreach my $k (1..$chr->{nb_chrs}){
+	    my $newFileName=$input_prefix.".".$chr->{$k}->{name}.".links";
+	    push(@{$files_list},$newFileName);
+	    my $fh = new FileHandle;
+	    $fh->open(">$newFileName");
+	    $filesHandle{$k}=$fh;
+	}
+    }
+
+    open LINKS, "<".$link_file or die "$0: can't open ".$link_file.":$!\n";
+    while(<LINKS>){
+
+	my @t=split;
+	my ($chr_read1,$chr_read2)=($t[0],$t[3]);
+
+	next unless (exists $chrID->{$chr_read1} && exists $chrID->{$chr_read2});
+
+	($chr_read1, $chr_read2)= ($chrID->{$chr_read1},$chrID->{$chr_read2});
+
+	if( ($sv_type=~/^(all|inter)$/) && ($chr_read1 ne $chr_read2) ){
+	    my $fh2print=$filesHandle{inter};
+	    print $fh2print join("\t",@t)."\n";
+	}
+
+	if( ($sv_type=~/^(all|intra)$/) && ($chr_read1 eq $chr_read2) ){
+	    my $fh2print=$filesHandle{$chr_read1};
+	    print $fh2print join("\t",@t)."\n";
+
+	}
+    }
+
+    foreach my $name (keys %filesHandle){
+	my $fh=$filesHandle{$name};
+	$fh->close;
+    }
+
+    print LOG "# Splitted link files created.\n";
+}
+
+
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+#MAIN FUNCTION number 2: Filtering processing
+sub getFilteredlinks {
+
+    my ($chr,$chrID,$tmp_links_prefix)=@_;
+    my %PAIR;
+
+    strandFiltering($chr,$chrID,
+		     $CONF{filtering}{nb_pairs_threshold},			#filtering of links
+		     $CONF{filtering}{strand_filtering},
+		     $CONF{filtering}{chromosomes},
+		     $CONF{general}{input_format},
+		     $CONF{general}{cmap_file},
+		     $CONF{general}{mates_orientation},
+		     $CONF{general}{read_lengths},
+		    $tmp_links_prefix,
+		    "$tmp_links_prefix.filtered",
+		    );
+
+    if($CONF{filtering}{strand_filtering}){					#re-definition of links coordinates with strand filtering
+
+	my @tmpfiles;
+
+	rename("$tmp_links_prefix.filtered","$tmp_links_prefix.filtered_unique");
+
+	getUniqueLinks("$tmp_links_prefix.filtered_unique",
+		       "$tmp_links_prefix.filtered");
+
+	push(@tmpfiles,"$tmp_links_prefix.filtered_unique");
+
+	if($CONF{filtering}{order_filtering}){					#filtering using the order
+
+	    rename("$tmp_links_prefix.filtered","$tmp_links_prefix.filtered_ordered");
+
+	    orderFiltering($chr,$chrID,
+			      $CONF{filtering}{nb_pairs_threshold},
+			      $CONF{filtering}{nb_pairs_order_threshold},
+			      $CONF{filtering}{mu_length},
+			      $CONF{filtering}{sigma_length},
+			      $CONF{general}{mates_orientation},
+			      $CONF{general}{read_lengths},
+			      "$tmp_links_prefix.filtered_ordered",
+			      "$tmp_links_prefix.filtered",
+			     );
+
+	    push(@tmpfiles,"$tmp_links_prefix.filtered_ordered");
+	}
+
+	if (($CONF{filtering}{insert_size_filtering})&&
+	    ($CONF{general}{sv_type} ne 'inter')){
+
+	    rename("$tmp_links_prefix.filtered","$tmp_links_prefix.filtered_withoutIndelSize");
+
+	    addInsertionInfo($chr,$chrID,
+			     $CONF{filtering}{nb_pairs_threshold},
+			     $CONF{filtering}{order_filtering},
+			     $CONF{filtering}{indel_sigma_threshold},
+			     $CONF{filtering}{dup_sigma_threshold},
+			     $CONF{filtering}{singleton_sigma_threshold},
+			     $CONF{filtering}{mu_length},
+			     $CONF{filtering}{sigma_length},
+			     $CONF{general}{mates_orientation},
+			     $CONF{general}{read_lengths},
+			     "$tmp_links_prefix.filtered_withoutIndelSize",
+			     "$tmp_links_prefix.filtered"
+			     );
+
+	    push(@tmpfiles,"$tmp_links_prefix.filtered_withoutIndelSize");
+	}
+
+	sortLinks("$tmp_links_prefix.filtered",
+		  "$tmp_links_prefix.filtered_sorted");
+
+	removeFullyOverlappedLinks("$tmp_links_prefix.filtered_sorted",
+			           "$tmp_links_prefix.filtered_nodup",
+				   );
+
+	postFiltering("$tmp_links_prefix.filtered_nodup",
+		      "$tmp_links_prefix.filtered",
+		       $CONF{filtering}{final_score_threshold});
+
+	push(@tmpfiles,"$tmp_links_prefix.filtered_sorted","$tmp_links_prefix.filtered_nodup");
+
+	unlink(@tmpfiles);
+
+
+    }
+    undef %PAIR;
+
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+#MAIN FUNCTION number 3: Circos format conversion for links
+sub links2circos{
+
+    my $input_file=$CONF{general}{mates_file}.".".$CONF{general}{sv_type}.".links.filtered";
+    my @path=split(/\//,$input_file);
+    $input_file=$CONF{general}{output_dir}.$path[$#path];
+
+    my $output_file.=$input_file.".segdup.txt";
+
+    links2segdup($CONF{circos}{organism_id},
+		 $CONF{circos}{colorcode},
+		 $input_file,
+		 $output_file);							#circos file output
+
+    system "rm $pt_circos_file; ln -s $output_file $pt_circos_file" if (defined $pt_circos_file); #GALAXY
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+#MAIN FUNCTION number 4: Bed format conversion for links
+sub links2bed{
+
+    my $input_file=$CONF{general}{mates_file}.".".$CONF{general}{sv_type}.".links.filtered";
+    my @path=split(/\//,$input_file);
+    $input_file=$CONF{general}{output_dir}.$path[$#path];
+
+    my $output_file.=$input_file.".bed.txt";
+
+    links2bedfile($CONF{general}{read_lengths},
+		  $CONF{bed}{colorcode},
+		  $input_file,
+		  $output_file);						#bed file output
+
+    system "rm $pt_bed_file; ln -s $output_file $pt_bed_file" if (defined $pt_bed_file); #GALAXY
+
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+#MAIN FUNCTION number 6: Bed format conversion for links
+sub links2SV{
+
+    my $input_file=$CONF{general}{mates_file}.".".$CONF{general}{sv_type}.".links.filtered";
+
+    my @path=split(/\//,$input_file);
+    $input_file=$CONF{general}{output_dir}.$path[$#path];
+
+    my $output_file.=$input_file.".sv.txt";
+
+
+    links2SVfile( $input_file,
+		  $output_file);
+
+     system "rm $pt_sv_file; ln -s $output_file $pt_sv_file" if (defined $pt_sv_file); #GALAXY
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+#MAIN FUNCTION number 7: copy number variations, coverage ratio calculation
+sub cnv{
+
+    my %CHR;
+    my %CHRID;
+    my @MATEFILES;
+    my @MATEFILES_REF;
+
+    my $output_prefix=$CONF{general}{mates_file};
+    my $output_prefix_ref=$CONF{detection}{mates_file_ref};
+    my @path=split(/\//,$output_prefix);
+    my @path_ref=split(/\//,$output_prefix_ref);
+    $output_prefix=$CONF{general}{output_dir}.$path[$#path];
+    $output_prefix_ref=$CONF{general}{output_dir}.$path_ref[$#path_ref];
+    my $tmp_mates_prefix=$CONF{general}{tmp_dir}."mates/".$path[$#path];
+    my $tmp_mates_prefix_ref=$CONF{general}{tmp_dir}."mates/".$path_ref[$#path_ref];
+    my $tmp_density_prefix=$CONF{general}{tmp_dir}."density/".$path[$#path];
+
+    shearingChromosome(\%CHR, \%CHRID,						#making the genomic fragment library with the detection parameters
+		     $CONF{detection}{window_size},
+		     $CONF{detection}{step_length},
+		     $CONF{general}{cmap_file});
+
+    if($CONF{detection}{split_mate_file}){
+
+	splitMateFile(\%CHR, \%CHRID, \@MATEFILES, $tmp_mates_prefix,
+		      "intra",
+		      $CONF{general}{mates_file},
+		      $CONF{general}{input_format},
+		      $CONF{general}{read_lengths}
+		      );
+
+	splitMateFile(\%CHR, \%CHRID, \@MATEFILES_REF, $tmp_mates_prefix_ref,
+		      "intra",
+		      $CONF{detection}{mates_file_ref},
+		      $CONF{general}{input_format},
+		      $CONF{general}{read_lengths}
+		      );
+
+
+    }else{
+
+	@MATEFILES=qx{ls $tmp_mates_prefix*} or die "# Error: No splitted sample mate files of \"$CONF{general}{mates_file}\" already created at $CONF{general}{tmp_dir} :$!";
+	chomp(@MATEFILES);
+	@MATEFILES_REF=qx{ls $tmp_mates_prefix_ref*} or die "# Error: No splitted reference mate files of \"$CONF{detection}{mates_file_ref}\" already created at $CONF{general}{tmp_dir} :$!";
+	chomp(@MATEFILES_REF);
+	print LOG "# Splitted sample and reference mate files already created.\n";
+    }
+
+    #Parallelization of the cnv per chromosome
+    my $pm = new Parallel::ForkManager($CONF{general}{num_threads});
+
+    foreach my $file (0..$#MATEFILES){
+
+	my $pid = $pm->start and next;
+
+	densityCalculation(\%CHR, \%CHRID, $file,
+	  $CONF{general}{read_lengths},
+	  $CONF{detection}{window_size},
+	  $CONF{detection}{step_length},
+	  \@MATEFILES,
+	  \@MATEFILES_REF,
+	  $MATEFILES[$file].".density",
+	  $CONF{general}{input_format});
+
+	$pm->finish;
+
+    }
+    $pm->wait_all_children;
+
+    #Merge the chromosome links file into only one
+    my @DENSITYFILES= qx{ls $tmp_density_prefix*density} or die "# Error: No density files created at $CONF{general}{tmp_dir} :$!";
+    chomp(@DENSITYFILES);
+    catFiles( \@DENSITYFILES => "$output_prefix.density" );
+
+    print LOG "# cnv end procedure : output created: $output_prefix.density\n";
+
+
+    undef %CHR;
+    undef %CHRID;
+
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+#MAIN FUNCTION number 8: Circos format conversion for cnv ratios
+sub ratio2circos{
+
+    my $input_file=$CONF{general}{mates_file}.".density";
+
+    my @path=split(/\//,$input_file);
+    $input_file=$CONF{general}{output_dir}.$path[$#path];
+
+    my $output_file.=$input_file.".segdup.txt";
+
+    ratio2segdup($CONF{circos}{organism_id},
+		 $input_file,
+		 $output_file);
+}
+#------------------------------------------------------------------------------#
+#MAIN FUNCTION number 9: BedGraph format conversion for cnv ratios
+sub ratio2bedgraph{
+
+    my $input_file=$CONF{general}{mates_file}.".density";
+
+    my @path=split(/\//,$input_file);
+    $input_file=$CONF{general}{output_dir}.$path[$#path];
+
+    my $output_file.=$input_file.".bedgraph.txt";
+
+    ratio2bedfile($input_file,
+		  $output_file);						#bed file output
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+#Creation of the fragment library
+sub shearingChromosome{
+
+    print LOG "# Making the fragments library...\n";
+
+    my ($chr,$chrID,$window,$step,$cmap_file)=@_;				#window and step sizes parameters
+
+    createChrHashTables($chr,$chrID,$cmap_file);				#hash tables: chromosome ID <=> chromsomes Name
+
+    foreach my $k (1..$chr->{nb_chrs}){
+
+	print LOG"-- $chr->{$k}->{name}\n";
+
+	my $frag=1;
+	for (my $start=0; $start<$chr->{$k}->{length}; $start+=$step){
+
+	    my $end=($start<($chr->{$k}->{length})-$window)? $start+$window-1:($chr->{$k}->{length})-1;
+	    $chr->{$k}->{$frag}=[$start,$end];					#creation of fragments, coordinates storage
+
+	    if($end==($chr->{$k}->{length})-1){
+		$chr->{$k}->{nb_frag}=$frag;					#nb of fragments per chromosome
+		last;
+	    }
+	    $frag++;
+	}
+    }
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+#Creation of chromosome hash tables from the cmap file
+sub createChrHashTables{
+
+    my ($chr,$chrID,$cmap_file)=@_;
+    $chr->{nb_chrs}=0;
+
+    open CMAP, "<".$cmap_file or die "$0: can't open ".$cmap_file.":$!\n";
+    while(<CMAP>){
+
+	if(/^\s+$/){ next;}
+	my ($k,$name,$length) = split;
+	$chr->{$k}->{name}=$name;
+	$chr->{$k}->{length}=$length;
+	$chrID->{$name}=$k;
+	$chr->{nb_chrs}++;
+
+    }
+    close CMAP;
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+#Read the mate file according the input format file (solid, eland or sam)
+sub readMateFile{
+
+    my ($chr1,$chr2,$pos1,$pos2,$order1,$order2,$t,$file_type,$tag_length)=@_;
+    my ($strand1,$strand2);
+
+    if($file_type eq "solid"){
+
+	($$chr1,$$chr2,$$pos1,$$pos2,$$order1,$$order2)=($$t[6],$$t[7],$$t[8]+1,$$t[9]+1,1,2); #0-based
+
+    }else{
+	my ($tag_length1,$tag_length2);
+	($$chr1,$$chr2,$$pos1,$strand1,$$pos2,$strand2,$$order1,$$order2,$tag_length1,$tag_length2)=($$t[11],$$t[12],$$t[7],$$t[8],$$t[9],$$t[10],1,2,length($$t[1]),length($$t[2])) #1-based
+	if($file_type eq "eland");
+
+	if($file_type eq "sam"){
+
+	    return 0 if ($$t[0]=~/^@/);						#header sam filtered out
+
+	    ($$chr1,$$chr2,$$pos1,$$pos2)=($$t[2],$$t[6],$$t[3],$$t[7]);
+
+	    return 0 if ($$chr1 eq "*" || $$chr2 eq "*");
+
+	    $$chr2=$$chr1 if($$chr2 eq "=");
+
+	    $strand1 = (($$t[1]&0x0010))? 'R':'F';
+	    $strand2 = (($$t[1]&0x0020))? 'R':'F';
+
+	    $$order1=  (($$t[1]&0x0040))? '1':'2';
+	    $$order2=  (($$t[1]&0x0080))? '1':'2';
+	    $tag_length1 = $tag_length->{$$order1};
+	    $tag_length2 = $tag_length->{$$order2};
+	}
+
+	$$pos1 = -($$pos1+$tag_length1) if ($strand1 eq "R");		#get sequencing starts
+	$$pos2 = -($$pos2+$tag_length2) if ($strand2 eq "R");
+    }
+    return 1;
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+#Parsing of the mates files and creation of links between 2 chromosomal fragments
+sub linking{
+
+    my ($chr,$chrID,$pair,$tag_length,$window_dist,$step,$mates_file,$input_format,$sv_type,$links_file)=@_;
+    my %link;
+
+    my $record=0;
+    my $nb_links=0;
+    my $warn=10000;
+
+    my @sfile=split(/\./,$mates_file);
+    my $fchr=$sfile[$#sfile];
+
+    my $fh = new FileHandle;
+
+    print LOG "# $fchr : Linking procedure...\n";
+    print LOG "-- file=$mates_file\n".
+	 "-- chromosome=$fchr\n".
+	 "-- input format=$input_format\n".
+	 "-- type=$sv_type\n".
+	 "-- read1 length=$tag_length->{1}, read2 length=$tag_length->{2}\n".
+	 "-- window size=$window_dist, step length=$step\n";
+
+    if ($mates_file =~ /.gz$/) {
+	$fh->open("gunzip -c $mates_file |") or die "$0: can't open ".$mates_file.":$!\n"; #gzcat
+    }elsif($mates_file =~ /.bam$/){
+	$fh->open("$SAMTOOLS_BIN_DIR/samtools view $mates_file |") or die "$0: can't open ".$mates_file.":$!\n";#GALAXY
+    }else{
+	$fh->open("<".$mates_file) or die "$0: can't open ".$mates_file.":$!\n";
+    }
+
+
+    while(<$fh>){
+
+	my @t=split;								#for each mate-pair
+	my $mate=$t[0];
+	my ($chr_read1, $chr_read2, $firstbase_read1, $firstbase_read2, $end_order_read1,$end_order_read2);
+
+	next if(exists $$pair{$mate});
+
+	next if (!readMateFile(\$chr_read1, \$chr_read2, \$firstbase_read1, \$firstbase_read2, \$end_order_read1, \$end_order_read2, \@t, $input_format,$tag_length));
+
+	next unless (exists $chrID->{$chr_read1} && exists $chrID->{$chr_read2});
+
+	($chr_read1, $chr_read2)= ($chrID->{$chr_read1},$chrID->{$chr_read2});
+
+	if($sv_type ne "all"){
+	    if( ($sv_type eq "inter") && ($chr_read1 ne $chr_read2) ||
+	        ($sv_type eq "intra") && ($chr_read1 eq $chr_read2) ){
+	    }else{
+		next;
+	    }
+	}
+
+	$$pair{$mate}=[$chr_read1, $chr_read2, $firstbase_read1, $firstbase_read2, $end_order_read1, $end_order_read2 ]; 	#fill out the hash pair table (ready for the defineCoordsLinks function)
+
+	$record++;
+
+	my ($coord_start_read1,$coord_end_read1,$coord_start_read2,$coord_end_read2);	#get the coordinates of each read
+
+	recupCoords($firstbase_read1,\$coord_start_read1,\$coord_end_read1,$tag_length->{$end_order_read1},$input_format);
+	recupCoords($firstbase_read2,\$coord_start_read2,\$coord_end_read2,$tag_length->{$end_order_read2},$input_format);
+
+	for(my $i=1;$i<=$chr->{$chr_read1}->{'nb_frag'};$i++){			#fast genome parsing for link creation
+
+	    if (abs ($coord_start_read1-${$chr->{$chr_read1}->{$i}}[0]) <= $window_dist){
+
+		if(overlap($coord_start_read1,$coord_end_read1,${$chr->{$chr_read1}->{$i}}[0],${$chr->{$chr_read1}->{$i}}[1])){
+
+		    for(my $j=1;$j<=$chr->{$chr_read2}->{'nb_frag'};$j++){
+
+			if (abs ($coord_start_read2-${$chr->{$chr_read2}->{$j}}[0]) <= $window_dist) {
+
+			    if(overlap($coord_start_read2,$coord_end_read2,${$chr->{$chr_read2}->{$j}}[0],${$chr->{$chr_read2}->{$j}}[1])){
+
+				makeLink(\%link,$chr_read1,$i,$chr_read2,$j,$mate,\$nb_links);	#make the link
+			    }
+
+			}else{
+
+			    $j=getNextFrag($coord_start_read2,$j,${$chr->{$chr_read2}->{$j}}[0],$chr->{$chr_read2}->{nb_frag},$window_dist,$step);
+			}
+		    }
+		}
+
+	    }else{
+
+		$i=getNextFrag($coord_start_read1,$i,${$chr->{$chr_read1}->{$i}}[0],$chr->{$chr_read1}->{nb_frag},$window_dist,$step);
+	    }
+	}
+
+	if($record>=$warn){
+	    print LOG "-- $fchr : $warn mate-pairs analysed - $nb_links links done\n";
+	    $warn+=10000;
+	}
+    }
+    $fh->close;
+
+    if(!$nb_links){
+	print LOG "-- $fchr : No mate-pairs !\n".
+	 "-- $fchr : No links have been found with the selected type of structural variations \($sv_type\)\n";
+    }
+
+    print LOG "-- $fchr : Total : $record mate-pairs analysed - $nb_links links done\n";
+
+    print LOG "-- $fchr : writing...\n";
+
+    $fh = new FileHandle;
+
+    $fh->open(">".$links_file) or die "$0: can't write in the output ".$links_file." :$!\n";
+
+    foreach my $chr1 ( sort { $a <=> $b} keys %link){				#Sorted links output
+
+	foreach my $chr2 ( sort { $a <=> $b} keys %{$link{$chr1}}){
+
+	    foreach my $frag1 ( sort { $a <=> $b} keys %{$link{$chr1}{$chr2}}){
+
+		foreach my $frag2 ( sort { $a <=> $b} keys %{$link{$chr1}{$chr2}{$frag1}}){
+
+		    my @count=split(",",$link{$chr1}{$chr2}{$frag1}{$frag2});
+		    print $fh "$chr->{$chr1}->{name}\t".(${$chr->{$chr1}->{$frag1}}[0]+1)."\t".(${$chr->{$chr1}->{$frag1}}[1]+1)."\t".
+				"$chr->{$chr2}->{name}\t".(${$chr->{$chr2}->{$frag2}}[0]+1)."\t".(${$chr->{$chr2}->{$frag2}}[1]+1)."\t".
+				scalar @count."\t".				#nb of read
+				$link{$chr1}{$chr2}{$frag1}{$frag2}."\n";	#mate list
+		}
+	    }
+	}
+    }
+
+    $fh->close;
+
+    undef %link;
+
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+#remove exact links doublons according to the mate list
+sub getUniqueLinks{
+
+    my ($links_file,$nrlinks_file)=@_;
+    my %links;
+    my %pt;
+    my $nb_links;
+    my $n=1;
+
+    my $record=0;
+    my $warn=300000;
+
+    my @sfile=split(/\./,$links_file);
+    my $fchr=$sfile[$#sfile-2];
+
+    my $fh = new FileHandle;
+
+    print LOG "# $fchr : Getting unique links...\n";
+    $fh->open("<$links_file") or die "$0: can't open $links_file :$!\n";
+
+    while(<$fh>){
+
+	my @t=split;
+	my $mates=$t[7];
+	$record++;
+
+	if(!exists $links{$mates}){						#Unique links selection
+
+	    $links{$mates}=[@t];
+	    $pt{$n}=$links{$mates};
+	    $n++;
+
+
+	}else{									#get the link coordinates from the mate-pairs list
+
+	    for my $i (1,2,4,5){						#get the shortest regions
+
+		$links{$mates}->[$i]=($t[$i]>$links{$mates}->[$i])? $t[$i]:$links{$mates}->[$i]	#maximum start
+		if($i==1 || $i==4);
+		$links{$mates}->[$i]=($t[$i]<$links{$mates}->[$i])? $t[$i]:$links{$mates}->[$i]	#minimum end
+		if($i==2 || $i==5);
+	    }
+	}
+	if($record>=$warn){
+	    print LOG "-- $fchr : $warn links analysed - ".($n-1)." unique links done\n";
+	    $warn+=300000;
+	}
+    }
+    $fh->close;
+
+    $nb_links=$n-1;
+    print LOG "-- $fchr : Total : $record links analysed - $nb_links unique links done\n";
+
+    $fh = new FileHandle;
+    $fh->open(">$nrlinks_file") or die "$0: can't write in the output: $nrlinks_file :$!\n";
+    print LOG "-- $fchr : writing...\n";
+    for my $i (1..$nb_links){
+
+	print $fh join("\t",@{$pt{$i}})."\n";				#all links output
+    }
+
+    $fh->close;
+
+    undef %links;
+    undef %pt;
+
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+#get the new coordinates of each link from the mate list
+sub defineCoordsLinks{
+
+    my ($chr,$chrID,$pair,$input_format,$sv_type,$tag_length,$links_file,$clinks_file)=@_;
+
+    my @sfile=split(/\./,$links_file);
+    my $fchr=$sfile[$#sfile-2];
+
+    my $fh = new FileHandle;
+    my $fh2 = new FileHandle;
+
+    $fh->open("<$links_file") or die "$0: can't open $links_file :$!\n";
+    $fh2->open(">$clinks_file") or die "$0: can't write in the output: $clinks_file :$!\n";
+
+    print LOG "# $fchr : Defining precise link coordinates...\n";
+
+    my $record=0;
+    my $warn=100000;
+
+    my %coords;
+    my %strands;
+    my %order;
+    my %ends_order;
+
+    while(<$fh>){
+
+
+	my ($col1,$col2)=(1,2);							#for an intrachromosomal link
+	my $diffchr=0;								#difference between chr1 and chr2
+	my ($chr1,$chr2,$mates_list,$npairs)=(split)[0,3,7,8];
+	($chr1,$chr2) = ($chrID->{$chr1},$chrID->{$chr2});
+	if ($chr1 != $chr2){							#for an interchromosomal link
+	    $col1=$col2=0;							#no distinction
+	    $diffchr=1;
+	}
+
+	my @pairs=split(",",$mates_list);
+
+	$coords{$col1}{$chr1}->{start}=undef;
+	$coords{$col1}{$chr1}->{end}=undef;
+	$coords{$col2}{$chr2}->{start}=undef;
+	$coords{$col2}{$chr2}->{end}=undef;
+	$strands{$col1}{$chr1}=undef;
+	$strands{$col2}{$chr2}=undef;
+	$ends_order{$col1}{$chr1}=undef;
+	$ends_order{$col2}{$chr2}=undef;
+
+
+	$order{$col1}{$chr1}->{index}->{1}=undef;
+	$order{$col1}{$chr1}->{index}->{2}=undef;
+	$order{$col2}{$chr2}->{index}->{1}=undef;
+	$order{$col2}{$chr2}->{index}->{2}=undef;
+	$order{$col1}{$chr1}->{order}=undef;
+	$order{$col2}{$chr2}->{order}=undef;
+
+	$record++;
+
+	for my $p (0..$#pairs){							#for each pair
+
+	    my ($coord_start_read1,$coord_end_read1);
+	    my ($coord_start_read2,$coord_end_read2);
+	    my $strand_read1=recupCoords(${$$pair{$pairs[$p]}}[2],\$coord_start_read1,\$coord_end_read1,$tag_length->{${$$pair{$pairs[$p]}}[4]},$input_format);
+	    my $strand_read2=recupCoords(${$$pair{$pairs[$p]}}[3],\$coord_start_read2,\$coord_end_read2,$tag_length->{${$$pair{$pairs[$p]}}[5]},$input_format);
+
+	    if(!$diffchr){							#for a intrachromosomal link
+		if($coord_start_read2<$coord_start_read1){				#get the closer start coordinate for each column
+		    ($col1,$col2)=(2,1);
+		}else{
+		    ($col1,$col2)=(1,2);
+		}
+	    }
+
+	    push(@{$coords{$col1}{${$$pair{$pairs[$p]}}[0]}->{start}},$coord_start_read1);	#get coords and strands of f3 and r3 reads
+	    push(@{$coords{$col1}{${$$pair{$pairs[$p]}}[0]}->{end}},$coord_end_read1);
+	    push(@{$coords{$col2}{${$$pair{$pairs[$p]}}[1]}->{start}},$coord_start_read2);
+	    push(@{$coords{$col2}{${$$pair{$pairs[$p]}}[1]}->{end}},$coord_end_read2);
+	    push(@{$strands{$col1}{${$$pair{$pairs[$p]}}[0]}},$strand_read1);
+	    push(@{$strands{$col2}{${$$pair{$pairs[$p]}}[1]}},$strand_read2);
+	    push(@{$ends_order{$col1}{${$$pair{$pairs[$p]}}[0]}},${$$pair{$pairs[$p]}}[4]);
+	    push(@{$ends_order{$col2}{${$$pair{$pairs[$p]}}[1]}},${$$pair{$pairs[$p]}}[5]);
+	}
+
+	($col1,$col2)=(1,2) if(!$diffchr);
+
+	my $coord_start_chr1=min(min(@{$coords{$col1}{$chr1}->{start}}),min(@{$coords{$col1}{$chr1}->{end}}));		#get the biggest region
+	my $coord_end_chr1=max(max(@{$coords{$col1}{$chr1}->{start}}),max(@{$coords{$col1}{$chr1}->{end}}));
+	my $coord_start_chr2=min(min(@{$coords{$col2}{$chr2}->{start}}),min(@{$coords{$col2}{$chr2}->{end}}));
+	my $coord_end_chr2=max(max(@{$coords{$col2}{$chr2}->{start}}),max(@{$coords{$col2}{$chr2}->{end}}));
+
+	@{$order{$col1}{$chr1}->{index}->{1}}= sort {${$coords{$col1}{$chr1}->{start}}[$a] <=>  ${$coords{$col1}{$chr1}->{start}}[$b]} 0 .. $#{$coords{$col1}{$chr1}->{start}};
+	@{$order{$col2}{$chr2}->{index}->{1}}= sort {${$coords{$col2}{$chr2}->{start}}[$a] <=>  ${$coords{$col2}{$chr2}->{start}}[$b]} 0 .. $#{$coords{$col2}{$chr2}->{start}};
+
+	foreach my $i (@{$order{$col1}{$chr1}->{index}->{1}}){			#get the rank of the chr2 reads according to the sorted chr1 reads (start coordinate sorting)
+	    foreach my $j (@{$order{$col2}{$chr2}->{index}->{1}}){
+
+		if(${$order{$col1}{$chr1}->{index}->{1}}[$i] == ${$order{$col2}{$chr2}->{index}->{1}}[$j]){
+		    ${$order{$col1}{$chr1}->{index}->{2}}[$i]=$i;
+		    ${$order{$col2}{$chr2}->{index}->{2}}[$i]=$j;
+		    last;
+		}
+	    }
+	}
+
+	foreach my $i (@{$order{$col1}{$chr1}->{index}->{2}}){			#use rank chr1 as an ID
+	    foreach my $j (@{$order{$col2}{$chr2}->{index}->{2}}){
+
+		if(${$order{$col1}{$chr1}->{index}->{2}}[$i] == ${$order{$col2}{$chr2}->{index}->{2}}[$j]){
+		    ${$order{$col1}{$chr1}->{order}}[$i]=$i+1;
+		    ${$order{$col2}{$chr2}->{order}}[$i]=$j+1;
+		    last;
+		}
+	    }
+	}
+
+	@pairs=sortTablebyIndex(\@{$order{$col1}{$chr1}->{index}->{1}},\@pairs);#sorting of the pairs, strands, and start coords from the sorted chr2 reads
+	@{$strands{$col1}{$chr1}}=sortTablebyIndex(\@{$order{$col1}{$chr1}->{index}->{1}},$strands{$col1}{$chr1});
+	@{$strands{$col2}{$chr2}}=sortTablebyIndex(\@{$order{$col1}{$chr1}->{index}->{1}},$strands{$col2}{$chr2});
+	@{$ends_order{$col1}{$chr1}}=sortTablebyIndex(\@{$order{$col1}{$chr1}->{index}->{1}},$ends_order{$col1}{$chr1});
+	@{$ends_order{$col2}{$chr2}}=sortTablebyIndex(\@{$order{$col1}{$chr1}->{index}->{1}},$ends_order{$col2}{$chr2});
+	@{$coords{$col1}{$chr1}->{start}}=sortTablebyIndex(\@{$order{$col1}{$chr1}->{index}->{1}},$coords{$col1}{$chr1}->{start});
+	@{$coords{$col2}{$chr2}->{start}}=sortTablebyIndex(\@{$order{$col1}{$chr1}->{index}->{1}},$coords{$col2}{$chr2}->{start});
+
+
+	my @link=($chr->{$chr1}->{name}, $coord_start_chr1 , $coord_end_chr1,	#all information output
+		  $chr->{$chr2}->{name}, $coord_start_chr2 , $coord_end_chr2,
+		  scalar @pairs,
+		  join(",",@pairs),
+		  join(",",@{$strands{$col1}{$chr1}}),
+		  join(",",@{$strands{$col2}{$chr2}}),
+		  join(",",@{$ends_order{$col1}{$chr1}}),
+		  join(",",@{$ends_order{$col2}{$chr2}}),
+		  join(",",@{$order{$col1}{$chr1}->{order}}),
+		  join(",",@{$order{$col2}{$chr2}->{order}}),
+		  join(",",@{$coords{$col1}{$chr1}->{start}}),
+		  join(",",@{$coords{$col2}{$chr2}->{start}}));
+
+	print $fh2 join("\t",@link)."\n";
+
+	if($record>=$warn){
+	    print LOG "-- $fchr : $warn links processed\n";
+	    $warn+=100000;
+	}
+    }
+    $fh->close;
+    $fh2->close;
+
+    print LOG "-- $fchr : Total : $record links processed\n";
+
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+#Sort links according the concerned chromosomes and their coordinates
+sub sortLinks{
+
+    my ($links_file,$sortedlinks_file,$unique)=@_;
+
+    my @sfile=split(/\./,$links_file);
+    my $fchr=$sfile[$#sfile-2];
+
+
+    print LOG "# $fchr : Sorting links...\n";
+
+    my $pipe=($unique)? "| sort -u":"";
+    system "sort -k 1,1 -k 4,4 -k 2,2n -k 5,5n -k 8,8n $links_file $pipe > $sortedlinks_file";
+
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+#removal of fully overlapped links
+sub removeFullyOverlappedLinks{
+
+    my ($links_file,$nrlinks_file,$warn_out)=@_;
+
+    my %pt;
+    my $n=1;
+
+    my @sfile=split(/\./,$links_file);
+    my $fchr=$sfile[$#sfile-2];
+
+    my $fh = new FileHandle;
+
+    $fh->open("<$links_file") or die "$0: can't open $links_file :$!\n";
+    while(<$fh>){
+
+	my @t=split("\t",$_);
+	$pt{$n}=[@t];
+	$n++;
+    }
+    $fh->close;
+
+    my $nb_links=$n-1;
+    my $nb=$nb_links;
+
+    my %pt2;
+    my $nb2=1;
+    my $record=0;
+    my $warn=10000;
+
+    print LOG "# $fchr : Removing fully overlapped links...\n";
+
+    LINK:
+
+    for my $i (1..$nb){
+
+	my @link=();
+	my @next_link=();
+	my $ind1=$i;
+
+	$record++;
+	if($record>=$warn){
+	    print LOG "-- $fchr : $warn unique links analysed - ".($nb2-1)." non-overlapped links done\n";
+	    $warn+=10000;
+	}
+
+	if(exists $pt{$ind1}){
+	    @link=@{$pt{$ind1}};                                                #link1
+	}else{
+	    next LINK;
+	}
+
+	my ($chr1,$start1,$end1,$chr2,$start2,$end2)=($link[0],$link[1],$link[2],$link[3],$link[4],$link[5]);	#get info of link1
+	my @mates=deleteBadOrderSensePairs(split(",",$link[7]));
+
+	my $ind2=$ind1+1;
+	$ind2++ while (!exists $pt{$ind2}&& $ind2<=$nb);			#get the next found link
+
+	if($ind2<=$nb){
+
+	    @next_link=@{$pt{$ind2}};					        #link2
+	    my ($chr3,$start3,$end3,$chr4,$start4,$end4)=($next_link[0],$next_link[1],$next_link[2],$next_link[3],$next_link[4],$next_link[5]); #get info of link2
+	    my @next_mates=deleteBadOrderSensePairs(split(",",$next_link[7]));
+
+	    while(($chr1 eq $chr3 && $chr2 eq $chr4) && overlap($start1,$end1,$start3,$end3)){	#loop here according to the chr1 coordinates, need an overlap between links to enter
+
+		if(!overlap($start2,$end2,$start4,$end4)){			#if no overlap with chr2 coordinates ->next link2
+
+		    $ind2++;
+		    $ind2++ while (!exists $pt{$ind2}&& $ind2<=$nb);
+
+		    if($ind2>$nb){						#if no more link in the file -> save link1
+
+			$pt2{$nb2}=\@link;
+			$nb2++;
+			next LINK;
+		    }
+
+		    @next_link=@{$pt{$ind2}};
+		    ($chr3,$start3,$end3,$chr4,$start4,$end4)=($next_link[0],$next_link[1],$next_link[2],$next_link[3],$next_link[4],$next_link[5]);
+		    @next_mates=deleteBadOrderSensePairs(split(",",$next_link[7]));
+		    next;
+		}
+
+		my %mates=map{$_ =>1} @mates;					#get the equal number of mates
+		my @same_mates = grep( $mates{$_}, @next_mates );
+		my $nb_mates= scalar @same_mates;
+
+		if($nb_mates == scalar @mates){
+
+		    delete $pt{$ind1};						#if pairs of link 1 are all included in link 2 -> delete link1
+		    next LINK;							#go to link2, link2 becomes link1
+
+		}else{
+			delete $pt{$ind2} if($nb_mates == scalar @next_mates);	#if pairs of link2 are all included in link 1 -> delete link2
+			$ind2++;							#we continue by checking the next link2
+			$ind2++ while (!exists $pt{$ind2}&& $ind2<=$nb);
+
+			if($ind2>$nb){						#if no more link in the file -> save link1
+
+			    $pt2{$nb2}=\@link;
+			    $nb2++;
+			    next LINK;
+			}
+
+			@next_link=@{$pt{$ind2}};					#get info of link2
+			($chr3,$start3,$end3,$chr4,$start4,$end4)=($next_link[0],$next_link[1],$next_link[2],$next_link[3],$next_link[4],$next_link[5]);
+			@next_mates=deleteBadOrderSensePairs(split(",",$next_link[7]));
+
+		}
+	    }
+	}
+	$pt2{$nb2}=\@link;							#if no (more) link with chr1 coordinates overlap -> save link1
+	$nb2++;
+    }
+
+    print LOG "-- $fchr : Total : $nb_links unique links analysed - ".($nb2-1)." non-overlapped links done\n";
+
+    #OUTPUT
+
+    $fh = new FileHandle;
+    $fh->open(">$nrlinks_file") or die "$0: can't write in the output: $nrlinks_file :$!\n";
+    print LOG "-- $fchr :  writing...\n";
+    for my $i (1..$nb2-1){
+
+	print $fh join("\t",@{$pt2{$i}});				#all links output
+    }
+
+    close $fh;
+
+    print LOG "-- $fchr : output created: $nrlinks_file\n" if($warn_out);
+
+    undef %pt;
+    undef %pt2;
+}
+#------------------------------------------------------------------------------#
+sub postFiltering {
+
+    my ($links_file,$pflinks_file, $finalScore_thres)=@_;
+
+    my @sfile=split(/\./,$links_file);
+    my $fchr=$sfile[$#sfile-2];
+
+
+    my ($nb,$nb2)=(0,0);
+
+    print LOG "# $fchr : Post-filtering links...\n";
+    print LOG "-- $fchr : final score threshold = $finalScore_thres\n";
+
+    my $fh = new FileHandle;
+    my $fh2 = new FileHandle;
+
+    $fh->open("<$links_file") or die "$0: can't open $links_file :$!\n";
+    $fh2->open(">$pflinks_file") or die "$0: can't write in the output: $pflinks_file :$!\n";
+
+
+    while(<$fh>){
+
+	my @t=split("\t",$_);
+	my $score=$t[$#t-1];
+
+	if($score >= $finalScore_thres){
+	    print $fh2 join("\t", @t);
+	    $nb2++;
+	}
+	$nb++;
+    }
+    $fh->close;
+    $fh2->close;
+
+    print LOG "-- $fchr : Total : $nb unique links analysed - $nb2 links kept\n";
+    print LOG "-- $fchr : output created: $pflinks_file\n";
+}
+
+
+
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+#Filtering of the links
+sub strandFiltering{
+
+    my($chr,$chrID,$pairs_threshold,$strand_filtering,$chromosomes,
+       $input_format,$cmap_file,$mate_sense, $tag_length,$links_file,$flinks_file)=@_;
+
+    my @sfile=split(/\./,$links_file);
+    my $fchr=$sfile[$#sfile-1];
+
+
+    my %chrs;
+    my %chrs1;
+    my %chrs2;
+    my $nb_chrs;
+    my $exclude;
+
+    if($chromosomes){
+	my @chrs=split(",",$chromosomes);
+	$nb_chrs=scalar @chrs;
+	$exclude=($chrs[0]=~/^\-/)? 1:0;
+	for my $chrName (@chrs){
+	    $chrName=~s/^(\-)//;
+	    my $col=($chrName=~s/_(1|2)$//);
+
+	    if(!$col){
+		$chrs{$chrID->{$chrName}}=undef
+	    }else{
+		$chrs1{$chrID->{$chrName}}=undef if($1==1);
+		$chrs2{$chrID->{$chrName}}=undef if($1==2);
+	    }
+	}
+    }
+
+    my $record=0;
+    my $nb_links=0;
+    my $warn=10000;
+
+    my $sens_ratio_threshold=0.6;
+
+    print LOG "\# Filtering procedure...\n";
+    print LOG "\# Number of pairs and strand filtering...\n";
+    print LOG "-- file=$links_file\n";
+    print LOG "-- nb_pairs_threshold=$pairs_threshold, strand_filtering=".(($strand_filtering)? "yes":"no").
+	 ", chromosomes=".(($chromosomes)? "$chromosomes":"all")."\n";
+
+
+
+    my $fh = new FileHandle;
+    my $fh2 = new FileHandle;
+
+    $fh->open("<$links_file") or die "$0: can't open $links_file :$!\n";
+    $fh2->open(">$flinks_file") or die "$0: can't write in the output: $flinks_file :$!\n";
+
+    while(<$fh>){
+
+	my @t=split;	#for each link
+	my $is_good=1;
+	$record++;
+
+
+	if($chromosomes){
+
+	    my ($chr1,$chr2)=($chrID->{$t[0]},$chrID->{$t[3]});
+
+	    if(!$exclude){
+		$is_good=(exists $chrs{$chr1} && exists $chrs{$chr2})? 1:0;
+		$is_good=(exists $chrs1{$chr1} && exists $chrs2{$chr2})? 1:0 if(!$is_good);
+		$is_good=($nb_chrs==1 && (exists $chrs1{$chr1} || exists $chrs2{$chr2}))? 1:0 if(!$is_good);
+	    }else{
+		$is_good=(exists $chrs{$chr1} || exists $chrs{$chr2})? 0:1;
+		$is_good=(exists $chrs1{$chr1} || exists $chrs2{$chr2})? 0:1 if($is_good);
+	    }
+	}
+
+	$is_good = ($is_good && $t[6] >= $pairs_threshold)? 1 :0;		#filtering according the number of pairs
+	if($is_good && $strand_filtering){					#if filtering according the strand sense
+
+	    my @mates=split(/,/,$t[7]);						#get the concordant pairs in the strand sense
+	    my @strands1=split(/,/,$t[8]);
+	    my @strands2=split(/,/,$t[9]);
+
+	    my %mate_class=( 'FF' => 0, 'RR' => 0, 'FR' => 0,  'RF' => 0);
+
+	    my %mate_reverse=( 'FF' => 'RR', 'RR' => 'FF',			#group1: FF,RR
+		       'FR' => 'RF', 'RF' => 'FR');				#group2: FR,RF
+
+	    my %mate_class2=( $mate_sense=>"NORMAL_SENSE", inverseSense($mate_sense)=>"NORMAL_SENSE",
+		    substr($mate_sense,0,1).inverseSense(substr($mate_sense,1,1))=>"REVERSE_SENSE",
+		    inverseSense(substr($mate_sense,0,1)).substr($mate_sense,1,1)=>"REVERSE_SENSE");
+
+	    if($t[6] == 1){
+
+		push(@t,$mate_class2{$strands1[0].$strands2[0]},"1/1",1,1);
+
+	    }else{
+
+		tie (my %class,'Tie::IxHash');
+		my $split;
+
+		foreach my $i (0..$#mates){
+		    $mate_class{$strands1[$i].$strands2[$i]}++;		#get the over-represented group
+		}
+
+		my $nb_same_sens_class=$mate_class{FF}+$mate_class{RR};
+		my $nb_diff_sens_class=$mate_class{FR}+$mate_class{RF};
+		my $sens_ratio=max($nb_same_sens_class,$nb_diff_sens_class)/($nb_same_sens_class+$nb_diff_sens_class);
+
+		if($sens_ratio < $sens_ratio_threshold){
+		    %class=(1=>'FF', 2=>'FR');
+		    $split=1;
+		}else{
+		    $class{1}=($nb_same_sens_class >  $nb_diff_sens_class)? 'FF':'FR';	#if yes get the concerned class
+		    $split=0;
+		}
+
+		$is_good=getConsistentSenseLinks(\@t,\@mates,\@strands1,\@strands2,$tag_length,$mate_sense,\%mate_reverse,\%mate_class2,\%class,$split,$pairs_threshold);
+	    }
+	}
+
+	if($is_good){								#PRINT
+
+	    my $nb=scalar @t;
+	    if($nb > 20){
+		my @t2=splice(@t,0,20);
+		print $fh2 join("\t",@t2)."\n";
+		$nb_links++;
+	    }
+	    $nb_links++;
+	    print $fh2 join("\t",@t)."\n";
+	}
+
+	if($record>=$warn){
+	   print LOG "-- $fchr : $warn links analysed - $nb_links links kept\n";
+	    $warn+=10000;
+	}
+    }
+    $fh->close;
+    $fh2->close;
+
+    print LOG "-- $fchr : No links have been found with the selected filtering parameters\n" if(!$nb_links);
+
+    print LOG "-- $fchr : Total : $record links analysed - $nb_links links kept\n";
+
+
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub getConsistentSenseLinks{
+
+    my ($t,$mates,$strands1,$strands2,$tag_length,$mate_sense, $mate_reverse,$mate_class2, $class, $split,$thres)=@_;
+
+    my $npairs=scalar @$mates;
+
+    my @ends_order1 = split (/,/,$$t[10]);
+    my @ends_order2 = split (/,/,$$t[11]);
+    my @order1 = split (/,/,$$t[12]);
+    my @order2 = split (/,/,$$t[13]);
+    my @positions1 = split (/,/,$$t[14]);
+    my @positions2 = split (/,/,$$t[15]);
+
+    my @newlink;
+
+    foreach my $ind (keys %{$class} ){
+
+	tie (my %flink,'Tie::IxHash');
+	my @orders2remove=();
+
+	foreach my $i (0..$#{$mates}){					#get the pairs belonging the over-represented group
+
+	    if((($$strands1[$i].$$strands2[$i]) eq $$class{$ind}) || (($$strands1[$i].$$strands2[$i]) eq $$mate_reverse{$$class{$ind}})){
+		push(@{$flink{mates}},$$mates[$i]);
+		push(@{$flink{strands1}},$$strands1[$i]);
+		push(@{$flink{strands2}},$$strands2[$i]);
+		push(@{$flink{ends_order1}},$ends_order1[$i]);
+		push(@{$flink{ends_order2}},$ends_order2[$i]);
+		push(@{$flink{positions1}},$positions1[$i]);
+		push(@{$flink{positions2}},$positions2[$i]);
+
+	    }else{
+		push(@orders2remove,$order1[$i]);
+	    }
+	}
+
+	@{$flink{order1}}=();
+	@{$flink{order2}}=();
+	if(scalar @orders2remove > 0){
+	    getNewOrders(\@order1,\@order2,\@orders2remove,$flink{order1},$flink{order2})
+	}else{
+	    @{$flink{order1}}=@order1;
+	    @{$flink{order2}}=@order2;
+	}
+
+	my @ends1; getEnds(\@ends1,$flink{positions1},$flink{strands1},$flink{ends_order1},$tag_length);
+	my @ends2; getEnds(\@ends2,$flink{positions2},$flink{strands2},$flink{ends_order2},$tag_length);
+
+	my $fnpairs=scalar @{$flink{mates}};
+	my $strand_filtering_ratio=$fnpairs."/".$npairs;
+	my $real_ratio=$fnpairs/$npairs;
+
+	if($fnpairs>=$thres){							#filtering according the number of pairs
+
+	    push(@newlink,
+		$$t[0],
+		min(min(@{$flink{positions1}}),min(@ends1)),
+		max(max(@{$flink{positions1}}),max(@ends1)),
+		$$t[3],
+		min(min(@{$flink{positions2}}),min(@ends2)),
+		max(max(@{$flink{positions2}}),max(@ends2)),
+		$fnpairs,
+		join(",",@{$flink{mates}}),
+		join(",",@{$flink{strands1}}),
+		join(",",@{$flink{strands2}}),
+		join(",",@{$flink{ends_order1}}),
+		join(",",@{$flink{ends_order2}}),
+		join(",",@{$flink{order1}}),
+		join(",",@{$flink{order2}}),
+		join(",",@{$flink{positions1}}),
+		join(",",@{$flink{positions2}}),
+		$$mate_class2{${$flink{strands1}}[0].${$flink{strands2}}[0]},
+		$strand_filtering_ratio,
+		$real_ratio,
+		$npairs
+	    );
+	}
+    }
+
+    if (grep {defined($_)} @newlink) {
+	@$t=@newlink;
+	return 1
+    }
+    return 0;
+
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub getNewOrders{
+
+    my($tab1,$tab2,$list,$newtab1,$newtab2)=@_;
+    my $j=1;
+    my $k=1;
+    for my $i (0..$#{$tab2}){
+	my $c=0;
+	for my $j (0..$#{$list}){
+	    $c++ if(${$list}[$j] < ${$tab2}[$i]);
+	    if(${$list}[$j] == ${$tab2}[$i]){
+		$c=-1; last;
+	    }
+	}
+	if($c!=-1){
+	    push(@{$newtab2}, ${$tab2}[$i]-$c);
+	    push(@{$newtab1}, $k);
+	    $k++;
+	}
+    }
+}
+
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+#Filtering of the links using their order
+sub orderFiltering {
+
+    my ($chr,$chrID,$nb_pairs_threshold,$nb_pairs_order_threshold,$mu,$sigma,$mate_sense,$tag_length,$links_file,$flinks_file)=@_;
+
+    my @sfile=split(/\./,$links_file);
+    my $fchr=$sfile[$#sfile-2];
+
+
+    my $diff_sense_ends=(($mate_sense eq "FR") || ($mate_sense eq "RF"))? 1:0;
+
+    my $record=0;
+    my $warn=10000;
+    my $nb_links=0;
+
+    my $quant05 = 1.644854;
+    my $quant001 = 3.090232;
+    my $alphaDist = $quant05 * 2 * $sigma;
+    my $maxFragmentLength = &floor($quant001 * $sigma + $mu);
+
+    print LOG "\# Filtering by order...\n";
+    print LOG "-- mu length=$mu, sigma length=$sigma, nb pairs order threshold=$nb_pairs_order_threshold\n";
+    print LOG "-- distance between comparable pairs was set to $alphaDist\n";
+    print LOG "-- maximal fragment length was set to $maxFragmentLength\n";
+
+
+    my $fh = new FileHandle;
+    my $fh2 = new FileHandle;
+
+    $fh->open("<$links_file") or die "$0: can't open $links_file :$!\n";
+    $fh2->open(">$flinks_file") or die "$0: can't write in the output: $flinks_file :$!\n";
+
+    while(<$fh>){
+
+	$record++;
+	my @t = split;
+	my ($chr1,$chr2,$mates_list)=@t[0,3,7];
+	my @pairs=split(",",$mates_list);
+	($chr1,$chr2) = ($chrID->{$chr1},$chrID->{$chr2});
+	my ($coord_start_chr1,$coord_end_chr1,$coord_start_chr2,$coord_end_chr2) = @t[1,2,4,5];
+	my $numberOfPairs = $t[6];
+	my @strand1 = split (/,/,$t[8]);
+	my @strand2 = split (/,/,$t[9]);
+	my @ends_order1 = split (/,/,$t[10]);
+	my @ends_order2 = split (/,/,$t[11]);
+	my @order1 = split (/,/,$t[12]);
+	my @order2 = split (/,/,$t[13]);
+	my @positions1 = split (/,/,$t[14]);
+	my @positions2 = split (/,/,$t[15]);
+	my @ends1; getEnds(\@ends1,\@positions1,\@strand1,\@ends_order1,$tag_length);
+	my @ends2; getEnds(\@ends2,\@positions2,\@strand2,\@ends_order2,$tag_length);
+	my $clusterCoordinates_chr1;
+	my $clusterCoordinates_chr2;
+	my $reads_left = 0;
+
+	my $ifRenv = $t[16];
+	my $strand_ratio_filtering=$t[17];
+
+	#kind of strand filtering. For example, will keep only FFF-RRR from a link FFRF-RRRF if <F-R> orientation is correct
+	my ($singleBreakpoint, %badInFRSense) = findBadInFRSenseSOLiDSolexa(\@strand1,\@strand2,\@ends_order1,\@ends_order2,\@order1,\@order2,$mate_sense);
+										#find pairs type F-RRRR or FFFF-R in the case if <R-F> orientation is correct
+										#These pairs are annotated as BED pairs forever! They won't be recycled!
+	my $table;
+	for my $i (0..$numberOfPairs-1) {					#fill the table with non adequate pairs: pairID	numberOfNonAdPairs nonAdPairIDs
+	    my $nonAdeq = 0;
+	    for my $j (0..$i-1) {
+		if (exists($table->{$j}->{$i})) {
+		    $nonAdeq++;
+		    $table->{$i}->{$j} = 1;
+		}
+	    }
+	    for my $j ($i+1..$numberOfPairs-1) {
+		if ($positions1[$j]-$positions1[$i]>$alphaDist) {
+		    if (&reversed ($i,$j,$ifRenv,\@positions2)) {
+			$nonAdeq++;
+			$table->{$i}->{$j} = 1;
+		    }
+		}
+	    }
+	    $table->{$i}->{nonAdeq} = $nonAdeq;
+	}
+
+	for my $bad (keys %badInFRSense) { #remove pairs type F-RRRR or FFFF-R in the case of <R-F> orientation
+	    &remove($bad,$table);
+	}
+
+	my @falseReads;
+							#RRRR-F -> RRRR or R-FFFF -> FFFF
+	@falseReads = findBadInRFSenseSOLiDSolexa(\@strand1,\@ends_order1,$mate_sense, keys %{$table});
+										#these pairs will be recycled later as $secondTable
+	for my $bad (@falseReads) {
+	    &remove($bad,$table);
+	}
+
+	my $bad = &check($table);
+	while ($bad ne "OK") {							#clear the table to reject non adequate pairs in the sense of ORDER
+	   # push (@falseReads, $bad);  remove completely!!!
+	    &remove($bad,$table);
+	    $bad = &check($table);
+	}
+
+	$reads_left = scalar keys %{$table};
+	my $coord_start_chr1_cluster1 = min(min(@positions1[sort {$a<=>$b} keys %{$table}]),min(@ends1[sort {$a<=>$b} keys %{$table}]));
+	my $coord_end_chr1_cluster1 = max(max(@positions1[sort {$a<=>$b} keys %{$table}]),max(@ends1[sort {$a<=>$b} keys %{$table}]));
+	my $coord_start_chr2_cluster1 = min(min(@positions2[sort {$a<=>$b} keys %{$table}]),min(@ends2[sort {$a<=>$b} keys %{$table}]));
+	my $coord_end_chr2_cluster1 = max(max(@positions2[sort {$a<=>$b} keys %{$table}]),max(@ends2[sort {$a<=>$b} keys %{$table}]));
+
+	$clusterCoordinates_chr1 = '('.$coord_start_chr1_cluster1.','.$coord_end_chr1_cluster1.')';
+	$clusterCoordinates_chr2 = '('.$coord_start_chr2_cluster1.','.$coord_end_chr2_cluster1.')';
+
+	my $ifBalanced = 'UNBAL';
+	my $secondTable;
+	my $clusterCoordinates;
+	my ($break_pont_chr1,$break_pont_chr2);
+
+	my $signatureType="";
+
+	my $maxCoord1 =$chr->{$chr1}->{length};
+	my $maxCoord2 =$chr->{$chr2}->{length};
+
+	if (scalar @falseReads) {
+	    @falseReads = sort @falseReads;
+	    #now delete FRFR choosing the majority
+	    my @newfalseReads;	   							#find and remove pairs type RRRR-F or R-FFFF
+	    @newfalseReads = findBadInRFSenseSOLiDSolexa(\@strand1,\@ends_order1,$mate_sense,@falseReads);	     #these @newfalseReads  won't be recycled
+	    my %hashTmp;
+	    for my $count1 (0..scalar(@falseReads)-1) {
+		my $i = $falseReads[$count1];
+		$hashTmp{$i} = 1;
+		for my $bad (@newfalseReads) {
+		    if ($bad == $i) {
+			delete $hashTmp{$i};
+			next;
+		    }
+		}
+	    }
+	    @falseReads = sort keys %hashTmp;  #what is left
+	    for my $count1 (0..scalar(@falseReads)-1) {				#fill the table for reads which were previously rejected
+		my $nonAdeq = 0;
+		my $i = $falseReads[$count1];
+
+		for my $count2 (0..$count1-1) {
+		    my $j = $falseReads[$count2];
+		    if (exists($secondTable->{$j}->{$i})) {
+			$nonAdeq++;
+			$secondTable->{$i}->{$j} = 1;
+		    }
+		}
+		for my $count2 ($count1+1..scalar(@falseReads)-1) {
+		    my $j = $falseReads[$count2];
+		    if ($positions1[$j]-$positions1[$i]>$alphaDist) {
+			if (&reversed ($i,$j,$ifRenv,\@positions2)) {
+			    $nonAdeq++;
+			    $secondTable->{$i}->{$j} = 1;
+			}
+		    }
+		}
+		$secondTable->{$i}->{nonAdeq} = $nonAdeq;
+	    }
+
+	    my @falseReads2;
+	    my $bad = &check($secondTable);
+	    while ($bad ne "OK") {						#clear the table to reject non adequate pairs
+		push (@falseReads2, $bad);
+		&remove($bad,$secondTable);
+		$bad = &check($secondTable);
+	    }
+	    if (scalar keys %{$secondTable} >= $nb_pairs_order_threshold) {
+		my $coord_start_chr1_cluster2 = min(min(@positions1[sort {$a<=>$b} keys %{$secondTable}]),min(@ends1[sort {$a<=>$b} keys %{$secondTable}]));
+		my $coord_end_chr1_cluster2 = max(max(@positions1[sort {$a<=>$b} keys %{$secondTable}]),max(@ends1[sort {$a<=>$b} keys %{$secondTable}]));
+		my $coord_start_chr2_cluster2 = min(min(@positions2[sort {$a<=>$b} keys %{$secondTable}]),min(@ends2[sort {$a<=>$b} keys %{$secondTable}]));
+		my $coord_end_chr2_cluster2 = max(max(@positions2[sort {$a<=>$b} keys %{$secondTable}]),max(@ends2[sort {$a<=>$b} keys %{$secondTable}]));
+
+		$ifBalanced = 'BAL';
+
+		if ($ifBalanced eq 'BAL') {
+
+		    if (scalar keys %{$table} < $nb_pairs_order_threshold) {
+			$ifBalanced = 'UNBAL';					#kill cluster 1!
+			($table,$secondTable)=($secondTable,$table);		#this means that one needs to exchange cluster1 with cluster2
+			$reads_left = scalar keys %{$table};
+			$coord_start_chr1_cluster1 = $coord_start_chr1_cluster2;
+			$coord_end_chr1_cluster1 = $coord_end_chr1_cluster2;
+			$coord_start_chr2_cluster1 = $coord_start_chr2_cluster2;
+			$coord_end_chr2_cluster1 = $coord_end_chr2_cluster2;
+			$clusterCoordinates_chr1 = '('.$coord_start_chr1_cluster1.','.$coord_end_chr1_cluster1.')';
+			$clusterCoordinates_chr2 = '('.$coord_start_chr2_cluster1.','.$coord_end_chr2_cluster1.')';
+
+		    } else {
+
+			$reads_left += scalar keys %{$secondTable};
+			next if ($reads_left < $nb_pairs_threshold);
+
+			if ($coord_end_chr1_cluster2 < $coord_start_chr1_cluster1) {
+			    ($table,$secondTable)=($secondTable,$table);	 #this means that one needs to exchange cluster1 with cluster2
+
+			    ($coord_start_chr1_cluster1,$coord_start_chr1_cluster2) = ($coord_start_chr1_cluster2,$coord_start_chr1_cluster1);
+			    ($coord_end_chr1_cluster1,$coord_end_chr1_cluster2)=($coord_end_chr1_cluster2,$coord_end_chr1_cluster1);
+			    ($coord_start_chr2_cluster1,$coord_start_chr2_cluster2)=($coord_start_chr2_cluster2,$coord_start_chr2_cluster1);
+			    ($coord_end_chr2_cluster1 , $coord_end_chr2_cluster2)=($coord_end_chr2_cluster2 , $coord_end_chr2_cluster1);
+
+			    $clusterCoordinates_chr1 = '('.$coord_start_chr1_cluster1.','.$coord_end_chr1_cluster1.'),'.$clusterCoordinates_chr1;
+			    $clusterCoordinates_chr2 = '('.$coord_start_chr2_cluster1.','.$coord_end_chr2_cluster1.'),'.$clusterCoordinates_chr2;
+			}else {
+			    $clusterCoordinates_chr1 .= ',('.$coord_start_chr1_cluster2.','.$coord_end_chr1_cluster2.')';
+			    $clusterCoordinates_chr2 .= ',('.$coord_start_chr2_cluster2.','.$coord_end_chr2_cluster2.')';
+			}
+			$coord_start_chr1 = min($coord_start_chr1_cluster1,$coord_start_chr1_cluster2);
+			$coord_end_chr1 = max($coord_end_chr1_cluster1,$coord_end_chr1_cluster2);
+			$coord_start_chr2 =  min($coord_start_chr2_cluster1,$coord_start_chr2_cluster2);
+			$coord_end_chr2 = max($coord_end_chr2_cluster1,$coord_end_chr2_cluster2);
+			#to calculate breakpoints one need to take into account read orientation in claster..
+			my $leftLetterOk = substr($mate_sense, 0, 1);  #R
+			my $rightLetterOk = substr($mate_sense, 1, 1); #F
+
+
+			my @index1 = keys %{$table};
+			my @index2 = keys %{$secondTable};
+
+			my (@generalStrand1,@generalStrand2) = 0;
+
+			if ($leftLetterOk eq $rightLetterOk) { #SOLID mate-pairs
+			    $leftLetterOk = 'R';
+			    $rightLetterOk = 'F';
+			    @generalStrand1 = translateSolidToRF(\@strand1,\@ends_order1);
+			    @generalStrand2 = translateSolidToRF(\@strand2,\@ends_order2);
+			} else {
+			    @generalStrand1 = @strand1;
+			    @generalStrand2 = @strand2; # TODO check if it is correct
+			}
+			if ($generalStrand1[$index1[0]] eq $leftLetterOk && $generalStrand1[$index2[0]] eq $rightLetterOk) { #(R,F)
+			    $break_pont_chr1 = '('.$coord_end_chr1_cluster1.','.$coord_start_chr1_cluster2.')';
+
+			    if ($generalStrand2[$index1[0]] eq $rightLetterOk && $generalStrand2[$index2[0]] eq $leftLetterOk) {
+				if ($coord_end_chr2_cluster1 >= $coord_end_chr2_cluster2) {
+				    $break_pont_chr2 = '('.$coord_end_chr2_cluster2.','.$coord_start_chr2_cluster1.')';
+				    $signatureType = "TRANSLOC";
+				} else {
+				    $break_pont_chr2 = '('.max(($coord_end_chr2_cluster1-$maxFragmentLength),1).','.$coord_start_chr2_cluster1.')';
+				    $break_pont_chr2 .= ',('.$coord_end_chr2_cluster2.','.min(($coord_start_chr2_cluster2+$maxFragmentLength),$maxCoord2).')';
+				    $signatureType = "INS_FRAGMT";
+				}
+
+			    } elsif ($generalStrand2[$index1[0]] eq $leftLetterOk && $generalStrand2[$index2[0]] eq $rightLetterOk) {
+				if ($coord_end_chr2_cluster1 >= $coord_end_chr2_cluster2) {
+				    $break_pont_chr2 = '('.max(($coord_end_chr2_cluster2-$maxFragmentLength),1).','.$coord_start_chr2_cluster2.')';
+				    $break_pont_chr2 .= ',('.$coord_end_chr2_cluster1.','.min(($coord_start_chr2_cluster1+$maxFragmentLength),$maxCoord2).')';
+				    $signatureType = "INV_INS_FRAGMT";
+				} else {
+				    $break_pont_chr2 = '('.$coord_end_chr2_cluster1.','.$coord_start_chr2_cluster2.')';
+				    $signatureType = "INV_TRANSLOC";
+				}
+			    } else {
+				#should not occur
+				print STDERR "\nError in orderFiltering\n\n";
+			    }
+			}
+
+			elsif ($generalStrand1[$index1[0]] eq $rightLetterOk && $generalStrand1[$index2[0]] eq $leftLetterOk) { #(F,R)
+			    $break_pont_chr1 = '('.max(($coord_end_chr1_cluster1-$maxFragmentLength),1).','.$coord_start_chr1_cluster1.')';
+			    $break_pont_chr1 .= ',('.$coord_end_chr1_cluster2.','.min(($coord_start_chr1_cluster2+$maxFragmentLength),$maxCoord1).')';
+			    if ($generalStrand2[$index1[0]] eq $rightLetterOk && $generalStrand2[$index2[0]] eq $leftLetterOk) {
+				if ($coord_end_chr2_cluster1 >= $coord_end_chr2_cluster2) {
+				    $break_pont_chr2 = '('.$coord_end_chr2_cluster2.','.$coord_start_chr2_cluster1.')';
+				    $signatureType = "INV_INS_FRAGMT";
+				} else {
+				    $break_pont_chr2 = '('.max(($coord_end_chr2_cluster1-$maxFragmentLength),1).','.$coord_start_chr2_cluster1.')';
+				    $break_pont_chr2 .= ',('.$coord_end_chr2_cluster2.','.min(($coord_start_chr2_cluster2+$maxFragmentLength),$maxCoord2).')';
+				    $signatureType = "INV_COAMPLICON";
+				}
+
+			    } elsif ($generalStrand2[$index1[0]] eq $leftLetterOk && $generalStrand2[$index2[0]] eq $rightLetterOk) {
+				if ($coord_end_chr2_cluster1 >= $coord_end_chr2_cluster2) {
+				    $break_pont_chr2 = '('.max(($coord_end_chr2_cluster2-$maxFragmentLength),1).','.$coord_start_chr2_cluster2.')';
+				    $break_pont_chr2 .= ',('.$coord_end_chr2_cluster1.','.min(($coord_start_chr2_cluster1+$maxFragmentLength),$maxCoord2).')';
+				    $signatureType = "COAMPLICON";
+				} else {
+				    $break_pont_chr2 = '('.$coord_end_chr2_cluster1.','.$coord_start_chr2_cluster2.')';
+				    $signatureType = "INS_FRAGMT";
+				}
+			    } else {
+				#should not occur
+				$signatureType = "UNDEFINED";
+			    }
+			}
+			else { # (F,F) or (R,R)  something strange. We will discard the smallest cluster
+			    $ifBalanced = 'UNBAL';
+			    if (scalar keys %{$secondTable} > scalar keys %{$table}) {
+				($table,$secondTable)=($secondTable,$table);		#this means that one needs to exchange cluster1 with cluster2
+
+				$coord_start_chr1_cluster1 = $coord_start_chr1_cluster2;
+				$coord_end_chr1_cluster1 = $coord_end_chr1_cluster2;
+				$coord_start_chr2_cluster1 = $coord_start_chr2_cluster2;
+				$coord_end_chr2_cluster1 = $coord_end_chr2_cluster2;
+				$clusterCoordinates_chr1 = '('.$coord_start_chr1_cluster1.','.$coord_end_chr1_cluster1.')';
+				$clusterCoordinates_chr2 = '('.$coord_start_chr2_cluster1.','.$coord_end_chr2_cluster1.')';
+			    }
+			    $reads_left = scalar keys %{$table};
+			}
+			if ($ifBalanced eq 'BAL') {
+			    $ifRenv = $signatureType;
+			}
+		    }
+		}
+	    }
+	}
+	if ($ifBalanced ne 'BAL') {
+	    #define possible break point
+	    $coord_start_chr1 = $coord_start_chr1_cluster1;
+	    $coord_end_chr1 = $coord_end_chr1_cluster1;
+	    $coord_start_chr2 = $coord_start_chr2_cluster1;
+	    $coord_end_chr2 = $coord_end_chr2_cluster1;
+
+	    my $region_length_chr1 = $coord_end_chr1-$coord_start_chr1;
+	    my $region_length_chr2 = $coord_end_chr2-$coord_start_chr2;
+
+	    my $leftLetterOk = substr($mate_sense, 0, 1);  #R
+	    my $rightLetterOk = substr($mate_sense, 1, 1); #F
+
+	    my @index = keys %{$table};
+	    unless ($diff_sense_ends) {
+		my $firstEndOrder1 = $ends_order1[$index[0]];
+		my $firstEndOrder2 = $ends_order2[$index[0]];
+	        $break_pont_chr1 = (($strand1[$index[0]] eq 'R' && $firstEndOrder1 == 2) || ($strand1[$index[0]] eq 'F' && $firstEndOrder1 == 1))?'('.$coord_end_chr1.','.min(($coord_start_chr1+$maxFragmentLength),$maxCoord1).')':'('.max(($coord_end_chr1-$maxFragmentLength),1).','.$coord_start_chr1.')';
+	        $break_pont_chr2 = (($strand2[$index[0]] eq 'R' && $firstEndOrder2 == 2) || ($strand2[$index[0]] eq 'F' && $firstEndOrder2 == 1))?'('.$coord_end_chr2.','.min(($coord_start_chr2+$maxFragmentLength),$maxCoord2).')':'('.max(($coord_end_chr2-$maxFragmentLength),1).','.$coord_start_chr2.')';
+	    } else {
+	        $break_pont_chr1 = ($strand1[$index[0]] eq $leftLetterOk )?'('.$coord_end_chr1.','.min(($coord_start_chr1+$maxFragmentLength),$maxCoord1).')':'('.max(($coord_end_chr1-$maxFragmentLength),1).','.$coord_start_chr1.')';
+	        $break_pont_chr2 = ($strand2[$index[0]] eq $leftLetterOk )?'('.$coord_end_chr2.','.min(($coord_start_chr2+$maxFragmentLength),$maxCoord2).')':'('.max(($coord_end_chr2-$maxFragmentLength),1).','.$coord_start_chr2.')';
+	    }
+
+	    if ($chr1 ne $chr2){
+		$ifRenv="INV_TRANSLOC" if($ifRenv eq "REVERSE_SENSE");
+		$ifRenv="TRANSLOC" if($ifRenv eq "NORMAL_SENSE");
+	    }
+	}
+
+	if (($ifBalanced eq 'BAL')&&( (scalar keys %{$table}) + (scalar keys %{$secondTable}) < $nb_pairs_threshold)) {
+	    next; 								#discard the link
+	}
+ 	if (($ifBalanced eq 'UNBAL')&&(scalar keys %{$table} < $nb_pairs_threshold)) {
+	    next; 								#discard the link
+	}
+	my $ratioTxt = "$reads_left/".(scalar @pairs);
+	my ($n1,$nTot) = split ("/",$strand_ratio_filtering);
+	my $ratioReal = $reads_left/$nTot;
+
+	if ($coord_start_chr1<=0) {
+	    $coord_start_chr1=1;
+	}
+	if ($coord_start_chr2<=0) {
+	    $coord_start_chr2=1;
+	}
+	#create output
+	my @link=($chr->{$chr1}->{name}, $coord_start_chr1 , $coord_end_chr1,	#all information output
+		  $chr->{$chr2}->{name}, $coord_start_chr2 , $coord_end_chr2,
+		  $reads_left,
+		  &redraw(1,$table,$secondTable,\%badInFRSense,$ifBalanced,\@pairs),
+		  &redraw(1,$table,$secondTable,\%badInFRSense,$ifBalanced,\@strand1),
+		  &redraw(1,$table,$secondTable,\%badInFRSense,$ifBalanced,\@strand2),
+		  &redraw(1,$table,$secondTable,\%badInFRSense,$ifBalanced,\@ends_order1),
+		  &redraw(1,$table,$secondTable,\%badInFRSense,$ifBalanced,\@ends_order2),
+		  &redraw(2,$table,$secondTable,\%badInFRSense,$ifBalanced,\@order1),
+		  &redraw(2,$table,$secondTable,\%badInFRSense,$ifBalanced,\@order2),
+		  &redraw(1,$table,$secondTable,\%badInFRSense,$ifBalanced,\@positions1),
+		  &redraw(1,$table,$secondTable,\%badInFRSense,$ifBalanced,\@positions2),
+		  $ifRenv,
+		  $strand_ratio_filtering,
+		  $ifBalanced, $ratioTxt, $break_pont_chr1, $break_pont_chr2,
+		  $ratioReal, $nTot);
+
+	$nb_links++;
+	print $fh2 join("\t",@link)."\n";
+
+	if($record>=$warn){
+	    print LOG "-- $fchr : $warn links analysed - $nb_links links kept\n";
+	    $warn+=10000;
+	}
+
+    }
+    $fh->close;
+    $fh2->close;
+
+    print LOG "-- $fchr : Total : $record links analysed - $nb_links links kept\n";
+
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+#gets information about ends positions given start, direction and order
+sub getEnds {
+    my ($ends,$starts,$strand,$end_order,$tag_length) = @_;
+    for my $i (0..scalar(@{$starts})-1) {
+	$ends->[$i] = getEnd($starts->[$i],$strand->[$i],$end_order->[$i],$tag_length);
+    }
+}
+sub getEnd {
+    my ($start,$strand, $end_order,$tag_length) = @_;
+    return ($strand eq 'F')? $start+$tag_length->{$end_order}-1:$start-$tag_length->{$end_order}+1;
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+#gets starts and ends Coords when start=leftmost given positions, directions and orders
+sub getCoordswithLeftMost {
+
+    my ($starts,$ends,$positions,$strand,$end_order,$tag_length) = @_;
+
+    for my $i (0..scalar(@{$positions})-1) {
+
+	if($strand->[$i] eq 'F'){
+	    $starts->[$i]=$positions->[$i];
+	    $ends->[$i]=$positions->[$i]+$tag_length->{$end_order->[$i]}-1;
+	}else{
+	    $starts->[$i]=$positions->[$i]-$tag_length->{$end_order->[$i]}+1;
+	    $ends->[$i]=$positions->[$i];
+	}
+    }
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub addInsertionInfo { 								#add field with INS,DEL,NA and distance between clusters and performs filtering
+
+    my ($chr,$chrID,$nb_pairs_threshold,$order_filtering,$indel_sigma_threshold,$dup_sigma_threshold,$singleton_sigma_threshold,$mu,$sigma,$mate_sense,$tag_length,$links_file,$flinks_file)=@_;
+
+    my @sfile=split(/\./,$links_file);
+    my $fchr=$sfile[$#sfile-2];
+
+
+    my $diff_sense_ends=(($mate_sense eq "FR") || ($mate_sense eq "RF"))? 1:0;
+
+    my $record=0;
+    my $nb_links=0;
+    my $warn=10000;
+
+    print LOG "\# Filtering out normal pairs using insert size...\n";
+    print LOG "-- mu length=$mu, sigma length=$sigma, indel sigma threshold=$indel_sigma_threshold, dup sigma threshold=$dup_sigma_threshold\n";
+    print LOG "-- using ".($mu-$indel_sigma_threshold*$sigma)."-".
+		      ($mu+$indel_sigma_threshold*$sigma)." as normal range of insert size for indels\n";
+    print LOG "-- using ".($mu-$dup_sigma_threshold*$sigma)."-".
+		      ($mu+$dup_sigma_threshold*$sigma)." as normal range of insert size for duplications\n";
+    print LOG "-- using ".($mu-$singleton_sigma_threshold*$sigma)." as the upper limit of insert size for singletons\n" if($mate_sense eq "RF");
+
+    my $fh = new FileHandle;
+    my $fh2 = new FileHandle;
+
+    $fh->open("<$links_file") or die "$0: can't open $links_file :$!\n";
+    $fh2->open(">$flinks_file") or die "$0: can't write in the output: $flinks_file :$!\n";
+
+    while(<$fh>){
+
+	$record++;
+	my @t = split;
+	my ($chr1,$chr2,$mates_list)=@t[0,3,7];
+
+	if($chrID->{$chr1} ne $chrID->{$chr2}) {				#if inter-chromosomal link here (because sv_type=all),
+	    $nb_links++;
+
+	    $t[16]="INV_TRANSLOC" if($t[16] eq "REVERSE_SENSE");
+	    $t[16]="TRANSLOC" if($t[16] eq "NORMAL_SENSE");
+
+	    $t[16].= "\t";
+	    $t[19].= "\t";
+
+	    print $fh2 join("\t",@t)."\n";
+
+	    if($record>=$warn){
+		print LOG "-- $fchr : $warn links processed - $nb_links links kept\n";
+		$warn+=10000;
+	    }
+	    next;
+	}
+
+	my $ifRenv = $t[16];
+	my $ifBalanced = "UNBAL";
+	$ifBalanced = $t[18] if ($order_filtering);
+
+	my $numberOfPairs = $t[6];
+	my @positions1 = deleteBadOrderSensePairs(split (/,/,$t[14]));
+	my @positions2 = deleteBadOrderSensePairs(split (/,/,$t[15]));
+
+	if ($ifBalanced eq "BAL") {
+
+	    if ($ifRenv eq "INV_TRANSLOC") {
+		$ifRenv = "INV_FRAGMT"; #for intrachromosomal inverted translocation is the same as inverted fragment
+	    }
+	    if ($ifRenv eq "NORMAL_SENSE") {
+		$ifRenv = "TRANSLOC";
+	    }
+	    if ($ifRenv eq "REVERSE_SENSE") {
+		$ifRenv = "INV_FRAGMT"; #for intrachromosomal inverted translocation is the same as inverted fragment
+	    }
+	    $t[19].= "\t";
+
+	    my $meanDistance = 0;
+
+	    for my $i (0..$numberOfPairs-1) {
+		$meanDistance += $positions2[$i]-$positions1[$i];
+	    }
+	    $meanDistance /= $numberOfPairs;
+
+	    $t[16] = $ifRenv."\t".$meanDistance;
+	    #dont touch the annotation. It should be already OK.
+
+	} else {
+	    #only for unbalanced
+
+	    my $ifoverlap=overlap($t[1],$t[2],$t[4],$t[5]);
+
+	    my $ends_sense_class = (deleteBadOrderSensePairs(split (/,/,$t[8])))[0].
+				   (deleteBadOrderSensePairs(split (/,/,$t[9])))[0];
+	    my $ends_order_class = (deleteBadOrderSensePairs(split (/,/,$t[10])))[0].
+				   (deleteBadOrderSensePairs(split (/,/,$t[11])))[0];
+
+	    my $indel_type = $ifRenv;
+
+	    my $meanDistance = "N/A";
+
+	    ($meanDistance, $indel_type) = checkIndel ($numberOfPairs,		#identify insertion type for rearrangments without inversion, calculates distance between cluster
+						       \@positions1,		#assign N/A to $indel_type if unknown
+						       \@positions2,
+						       $ifRenv,
+						       $ifoverlap,
+						       $indel_sigma_threshold,
+						       $dup_sigma_threshold,
+						       $singleton_sigma_threshold,
+						       $mu,
+						       $sigma,
+						       $ifBalanced,
+						       $ends_sense_class,
+						       $ends_order_class,
+						       $mate_sense,
+						       $diff_sense_ends,
+						       );
+
+	    #filtering of pairs with distance inconsistant with the SV
+	    if ($ifRenv ne "REVERSE_SENSE") {
+		my $maxCoord1 =$chr->{$chrID->{$chr1}}->{length};
+		my $maxCoord2 =$chr->{$chrID->{$chr2}}->{length};
+		$meanDistance = recalc_t_usingInsertSizeInfo(\@t,$mu,$sigma,$meanDistance,$tag_length,$diff_sense_ends,$mate_sense,
+							     $maxCoord1,$maxCoord2,$ends_sense_class,$ends_order_class,$nb_pairs_threshold,$order_filtering);
+		next if ($t[6] < $nb_pairs_threshold);
+	    }else{
+		 $t[19].= "\t";
+	    }
+	    $t[16] = $indel_type."\t".$meanDistance;
+	}
+
+	$nb_links++;
+
+	print $fh2 join("\t",@t)."\n";
+	if($record>=$warn){
+	     print LOG "-- $fchr : $warn links processed - $nb_links links kept\n";
+	    $warn+=10000;
+	}
+    }
+    $fh->close;
+    $fh2->close;
+
+    print LOG "-- $fchr : Total : $record links analysed - $nb_links links kept\n";
+
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub checkIndel  {
+
+    my ($numberOfPairs,	$positions1, $positions2, $ifRenv, $ifoverlap, $indel_sigma_threshold, $dup_sigma_threshold, $singleton_sigma_threshold,
+	$mu, $sigma, $ifBalanced,$ends_sense_class,$ends_order_class,$mate_sense,$diff_sense_ends) = @_;
+
+    my $meanDistance = 0;
+
+    for my $i (0..$numberOfPairs-1) {
+	$meanDistance += $positions2->[$i]-$positions1->[$i];
+    }
+    $meanDistance /= $numberOfPairs;
+
+    return ($meanDistance,"INV_DUPLI") if (($ifRenv eq "REVERSE_SENSE") && ($meanDistance<$mu+$dup_sigma_threshold*$sigma) );
+
+    return ($meanDistance,"INVERSION") if ($ifRenv eq "REVERSE_SENSE");
+
+    if($diff_sense_ends){
+	return ($meanDistance, "LARGE_DUPLI") if ($ends_sense_class ne $mate_sense) && ($meanDistance>$mu+$dup_sigma_threshold*$sigma) ;
+	return ($meanDistance, "SINGLETON") if (($meanDistance<$mu-$singleton_sigma_threshold*$sigma) && $mate_sense eq "RF" && ($ends_sense_class eq inverseSense($mate_sense)));
+    }else{
+	return ($meanDistance, "LARGE_DUPLI") if (($ends_sense_class eq $mate_sense) && ($ends_order_class eq "12") || ($ends_sense_class eq inverseSense($mate_sense)) && ($ends_order_class eq "21")) &&
+	($meanDistance>$mu+$dup_sigma_threshold*$sigma) ;
+    }
+
+    return ($meanDistance, "SMALL_DUPLI") if (($meanDistance<$mu-$dup_sigma_threshold*$sigma) && $ifoverlap);
+
+    return ($meanDistance, "DUPLICATION") if ($diff_sense_ends && ($ends_sense_class ne $mate_sense) && ($meanDistance<$mu-$dup_sigma_threshold*$sigma) ) ;
+
+    return ($meanDistance, "INSERTION") if ($meanDistance<$mu -$indel_sigma_threshold*$sigma);
+    return ($meanDistance, "DELETION") if ($meanDistance>$mu+$indel_sigma_threshold*$sigma);
+
+    return ($meanDistance, "UNDEFINED");
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+#sub reacalulate @t so that get rid of unconsistent pairs (unconsistent insert size )
+sub recalc_t_usingInsertSizeInfo {
+    my($t,$mu,$sigma,$meanDistance,$tag_length,$diff_sense_ends,$mate_sense,$maxCoord1,$maxCoord2,$ends_sense_class,$ends_order_class,$nb_pairs_threshold,$order_filtering) = @_;
+
+    my @badPairs;
+
+    my @positions1 = getAllEntries($t->[14]);
+    my @positions2 = getAllEntries($t->[15]);
+
+    if ($meanDistance < $mu) {
+	for my $i (0..scalar(@positions1)-1) {
+	    if (substr($positions2[$i],-1,1) ne '$' && substr($positions2[$i],-1,1) ne '*' && $positions2[$i]-$positions1[$i]>=$mu) {
+	    	push(@badPairs,$i);
+	    }
+    	}
+    } else {
+	for my $i (0..scalar(@positions1)-1) {
+	    if (substr($positions2[$i],-1,1) ne '$' && substr($positions2[$i],-1,1) ne '*' && $positions2[$i]-$positions1[$i]<=$mu) {
+		push(@badPairs,$i);
+	    }
+	}
+    }
+
+    if (scalar (@badPairs)>0) {
+	#print join("\t",@badPairs).": ".join("\t",@t)."\n";
+	#remove these inconsistant links
+	$t->[6] -= scalar(@badPairs); #numberOfPairs
+	return if ($t->[6] < $nb_pairs_threshold);
+
+	$t->[7] = mark_values(\@badPairs, $t->[7]);
+	$t->[8] = mark_values(\@badPairs, $t->[8]);
+	$t->[9] = mark_values(\@badPairs, $t->[9]);
+	$t->[10] = mark_values(\@badPairs, $t->[10]);
+	$t->[11] = mark_values(\@badPairs, $t->[11]);
+
+	$t->[12] = mark_indexes(\@badPairs, $t->[12]);
+	$t->[13] = mark_indexes(\@badPairs, $t->[13]);
+
+	$t->[14] = mark_values(\@badPairs, $t->[14]);
+	$t->[15] = mark_values(\@badPairs, $t->[15]);
+	$t->[19] = recalculate_ratio($t->[6],$t->[19]) if ($order_filtering); #add the second ratio
+	$t->[17] = recalculate_ratio($t->[6],$t->[17]) unless ($order_filtering);
+	($t->[1],$t->[2]) = recalculate_boundaries($t->[14],$t->[8],$t->[10],$tag_length);
+	($t->[4],$t->[5]) = recalculate_boundaries($t->[15],$t->[9],$t->[11],$tag_length);
+
+	#recalc breakpoints:
+	my $quant001 = 3.090232;
+	my $maxFragmentLength = &floor($quant001 * $sigma + $mu);
+	$t->[20] = recalc_breakpoints($mate_sense,$maxCoord1,$t->[14],substr($ends_sense_class,0,1),substr($ends_order_class,0,1),$t->[1],$t->[2],$maxFragmentLength,$diff_sense_ends ) if ($order_filtering);
+	$t->[21] = recalc_breakpoints($mate_sense,$maxCoord2,$t->[15],substr($ends_sense_class,1,1),substr($ends_order_class,1,1),$t->[4],$t->[5],$maxFragmentLength,$diff_sense_ends ) if ($order_filtering);
+	#recalc total ratio
+	$t->[22] = $t->[6] / $t->[23] if ($order_filtering);
+	$t->[18] = $t->[6] / $t->[19] unless ($order_filtering);
+
+    	@positions1 = deleteBadOrderSensePairs(split (/,/,$t->[14]));
+    	@positions2 = deleteBadOrderSensePairs(split (/,/,$t->[15]));
+
+	$meanDistance = 0;
+
+	for my $i (0..scalar(@positions1)-1) {
+	    $meanDistance += $positions2[$i]-$positions1[$i];
+	}
+	$meanDistance /= scalar(@positions1);
+
+    } else {
+	$t->[17] = recalculate_ratio((split(/\//,$t->[17]))[0],$t->[17]) unless ($order_filtering);
+	$t->[19] = recalculate_ratio((split(/\//,$t->[19]))[0],$t->[19]) if ($order_filtering);
+
+    } #nothing has been filtered
+    return $meanDistance;
+}
+
+sub recalculate_ratio {
+    my ($left, $ratio) = @_;
+    my @elements = split (/\//,$ratio);
+    $elements[1]= $elements[0];
+    $elements[0]=$left;
+    return $ratio."\t".join("/",@elements);
+}
+
+sub recalc_breakpoints {
+    my ($mate_sense,$maxCoord,$startString,$strand,$firstEndOrder,$coord_start_chr,$coord_end_chr,$maxFragmentLength,$diff_sense_ends ) = @_;
+    my $break_pont_chr;
+
+    my $leftLetterOk = substr($mate_sense, 0, 1);  #R
+    my $rightLetterOk = substr($mate_sense, 1, 1); #F
+
+
+    my @positions = deleteBadOrderSensePairs(split (/,/,$startString));
+
+    unless ($diff_sense_ends) {
+	$break_pont_chr = (($strand eq 'R' && $firstEndOrder == 2) || ($strand eq 'F' && $firstEndOrder == 1))?'('.$coord_end_chr.','.min(($coord_start_chr+$maxFragmentLength),$maxCoord).')':'('.max(($coord_end_chr-$maxFragmentLength),1).','.$coord_start_chr.')';
+    } else {
+	$break_pont_chr = ($strand eq $leftLetterOk)?'('.$coord_end_chr.','.min(($coord_start_chr+$maxFragmentLength),$maxCoord).')':'('.max(($coord_end_chr-$maxFragmentLength),1).','.$coord_start_chr.')';
+    }
+    return $break_pont_chr;
+}
+sub recalculate_boundaries {
+    my ($startString,$senseString,$endsOrderString,$tag_length) = @_;
+    my @positions = deleteBadOrderSensePairs(split (/,/,$startString));
+    my @strands = deleteBadOrderSensePairs(split (/,/,$senseString));
+    my @ends_orders = deleteBadOrderSensePairs(split (/,/,$endsOrderString));
+    my @ends; getEnds(\@ends,\@positions,\@strands,\@ends_orders,$tag_length);
+    my $coord_start_cluster = min(min(@positions),min(@ends));
+    my $coord_end_cluster = max(max(@positions),max(@ends));
+    return ($coord_start_cluster,$coord_end_cluster);
+}
+
+sub remove_indexes {
+    my ($bads, $string) = @_;
+    my @elements = deleteBadOrderSensePairs(split (/,/,$string));
+    for my $i (reverse sort %{$bads}) {
+	delete $elements[$i];
+    }
+    return "(".join(",",@elements).")";
+}
+##add @ to to elements
+sub mark_values {
+    my ($bads, $string) = @_;
+    my @elements = getAllEntries($string);
+    for my $i (@{$bads}) {
+	$elements[$i] .= "@";
+    }
+    return "(".join(",",@elements).")";
+}
+##add @ to to indexes
+sub mark_indexes {
+    my ($bads, $string) = @_;
+    my @elements = getAllEntries($string);
+    for my $i ((0..scalar(@elements)-1)) {
+	for my $j (@{$bads}) {
+	    $elements[$i] .= "@" if ($elements[$i] eq ($j+1));
+	}
+    }
+
+    return "(".join(",",@elements).")";
+}
+
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub redraw {
+
+    my ($type,$table,$secondTable,$badInFRSense,$ifBalanced,$arr) = @_;
+
+    my $out;
+    my @first_arr;
+    if ($ifBalanced eq 'BAL') {
+	my @second_arr;
+	my $lastPushed = 1;
+	if ($type == 1) {
+	    for my $i (0 .. scalar(@{$arr})-1) {
+		if (exists ($table->{$i})) {
+	    	    push(@first_arr,$arr->[$i]);
+		    $lastPushed = 1;
+		}elsif (exists ($secondTable->{$i})) {
+		    push(@second_arr,$arr->[$i]);
+		    $lastPushed = 2;
+		} elsif ($lastPushed == 1) {
+		    if (exists ($badInFRSense->{$i})) {
+		        push(@first_arr,$arr->[$i]."\$");
+		    }else {
+		        push(@first_arr,$arr->[$i]."*");
+		    }
+		} elsif ($lastPushed == 2) {
+		    if (exists ($badInFRSense->{$i})) {
+		        push(@second_arr,$arr->[$i]."\$");
+		    }else {
+			push(@second_arr,$arr->[$i]."*");
+		    }
+		} else {print "Error!";exit;}
+	    }
+	} else {
+	     for my $i (@{$arr}) {
+		if (exists ($table->{$i-1})) {
+	    	    push(@first_arr,$i);
+		    $lastPushed = 1;
+		}elsif (exists ($secondTable->{$i-1})) {
+		    push(@second_arr,$i);
+		    $lastPushed = 2;
+		} elsif ($lastPushed == 1) {
+		    if (exists ($badInFRSense->{$i-1})) {
+		        push(@first_arr,$i."\$");
+		    }else {
+			push(@first_arr,$i."*");
+		    }
+		} elsif ($lastPushed == 2) {
+		    if (exists ($badInFRSense->{$i-1})) {
+		        push(@second_arr,$i."\$");
+		    }else {
+			push(@second_arr,$i."*");
+		    }
+		} else {print "Error!";exit;}
+	    }
+	}
+	$out = '('.join(",",@first_arr).'),('.join(",",@second_arr).')';
+    }
+    else {
+	if ($type == 1) {
+	    for my $i (0 .. scalar(@{$arr})-1) {
+		if (exists ($table->{$i})) {
+	    	    push(@first_arr,$arr->[$i]);
+		} else {
+		    if (exists ($badInFRSense->{$i})) {
+			push(@first_arr,$arr->[$i]."\$");
+		    }else {
+		        push(@first_arr,$arr->[$i]."*");
+		    }
+		}
+	    }
+	} else {
+	    for my $i (@{$arr}) {
+		if (exists ($table->{$i-1})) {
+	    	    push(@first_arr,$i);
+		} else {
+		    if (exists ($badInFRSense->{$i-1})) {
+			push(@first_arr,$i."\$");
+		    }else {
+			push(@first_arr,$i."*");
+		    }
+		}
+	    }
+	}
+	$out = '('.join(",",@first_arr).')';
+    }
+    return $out;
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub check {
+
+    my $table = $_[0];
+    my $bad = 'OK';
+    my $max = 0;
+    for my $i (sort {$a<=>$b} keys %{$table}) {
+	unless ($table->{$i}->{nonAdeq} == 0) {
+	    if ($max<$table->{$i}->{nonAdeq}) {
+		$max=$table->{$i}->{nonAdeq};
+		$bad = $i;
+	    }
+	}
+    }
+    return $bad;
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub reversed {
+
+    my ($i,$j,$ifRenv,$positions) = @_;
+    if (($ifRenv eq 'REVERSE_SENSE' && $positions->[$i]<$positions->[$j]) || ($ifRenv ne 'REVERSE_SENSE' && $positions->[$i]>$positions->[$j])){
+	return 1;
+    }
+    return 0;
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub remove {
+
+    my ($bad,$table) = @_;
+    for my $i (sort {$a<=>$b} keys %{$table}) {
+	if ($bad == $i) {
+	    delete($table->{$i});;
+	} else {
+	    if (exists($table->{$i}->{$bad})) {
+		delete($table->{$i}->{$bad});
+		$table->{$i}->{nonAdeq}--;
+	    }
+	}
+    }
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub findBadInRFSenseSOLiDSolexa { #choose maximum: FFFFs or RRRRs
+
+    my ($strand,$ends_order,$mate_sense,@keysLeft) = @_;
+
+    my $leftLetterOk = substr($mate_sense, 0, 1);  #R
+    my $rightLetterOk = substr($mate_sense, 1, 1); #F
+
+    my (@standardArray);
+    if ($leftLetterOk eq $rightLetterOk) { #SOLID mate-pairs
+	$leftLetterOk = 'R';
+	$rightLetterOk = 'F';
+	@standardArray = translateSolidToRF($strand,$ends_order);
+    } else {
+	@standardArray = @{$strand};
+    }
+
+    my $ifR = 0;
+    my @Rs;
+
+    for my $i (@keysLeft) {
+	if ($standardArray[$i] eq $leftLetterOk) {
+	    $ifR++;
+	    push(@Rs,$i);
+	}
+    }
+
+
+    my $ifF = 0;
+    my @Fs;
+
+    for my $i (@keysLeft) {
+	if ($standardArray[$i] eq $rightLetterOk) {
+	    $ifF++;
+	    push(@Fs,$i);
+	}
+    }
+
+    if($ifR>=$ifF) {
+	return @Fs;
+    }
+    return @Rs;
+}
+
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub findBadInFRSenseSOLiDSolexa { #should work both for SOLiD and Solexa
+
+    my ($strand1,$strand2,$ends_order1,$ends_order2,$order1,$order2) = ($_[0],$_[1],$_[2],$_[3],$_[4],$_[5]);
+    my $mate_sense =  $_[6];
+
+    my $leftLetterOk = substr($mate_sense, 0, 1);  #R
+    my $rightLetterOk = substr($mate_sense, 1, 1); #F
+
+    my (@standardArray1,@standardArray2);
+
+    if ($leftLetterOk eq $rightLetterOk) { #SOLID mate-pairs
+	$leftLetterOk = 'R';
+	$rightLetterOk = 'F';
+	@standardArray1 = translateSolidToRF($strand1,$ends_order1);
+	my @arr = getOrderedStrands($strand2,$order2);
+	my @ends2 = getOrderedStrands($ends_order2,$order2);
+	@standardArray2 = translateSolidToRF(\@arr,\@ends2);
+
+    } else {
+	@standardArray1 = @{$strand1};
+	@standardArray2 = getOrderedStrands($strand2,$order2);
+    }
+
+    #we will try 4 possibilities, 2 for each end of the link: RFRR-FFF->RFFFF , RFRR-FFF->RRRFFF
+
+    #for the first end:
+
+    my @array = @standardArray1;
+    my %badInFRSense1;
+    for my $i (1..scalar (@array)-1){ # FRFRFFFF -> FFFFFF and RRFRFRFFFF -> RRFFFFFF
+	if ($array[$i-1] eq $rightLetterOk && $array[$i] eq $leftLetterOk) {
+	    $badInFRSense1{$i}=1;
+	    $array[$i] = $rightLetterOk;
+	}
+    }
+    my $numberRRRFFF_or_FFF_1 = scalar(@array)-scalar(keys %badInFRSense1);
+    @array = @standardArray1;
+    my %badInFRSense0;
+    for my $i (reverse(1..scalar (@array)-1)){ # FRFRFFFFRR -> FFFFFFRR
+	if ($array[$i-1] eq $rightLetterOk && $array[$i] eq $leftLetterOk) {
+	    $badInFRSense0{$i-1}=1;
+	    $array[$i-1] = $leftLetterOk;
+
+	}
+    }
+    my $numberRRF1 = scalar(@array)-scalar(keys %badInFRSense0);
+
+    #for the second end:
+    @array = @standardArray2;
+
+    my %badInFRSense3;
+    for my $i (1..scalar(@array)-1){
+	if ($array[$i-1] eq $rightLetterOk && $array[$i] eq $leftLetterOk) {
+	    $badInFRSense3{$order2->[$i]}=1;
+	    $array[$i] = $rightLetterOk;
+	}
+    }
+    my $numberRRRFFF_or_FFF_2 = scalar(@array)-scalar(keys %badInFRSense3);
+
+    @array = @standardArray2;
+    my %badInFRSense5;
+    for my $i (reverse(1..scalar (@array)-1)){ # FRFRFFFF -> FFFFFF
+	if ($array[$i-1] eq $rightLetterOk && $array[$i] eq $leftLetterOk) {
+	    $badInFRSense5{$i-1}=1;
+	    $array[$i-1] = $leftLetterOk;
+	}
+    }
+    my $numberRRF2 = scalar(@array)-scalar(keys %badInFRSense5);
+
+    if ($numberRRF1>=$numberRRRFFF_or_FFF_1 && $numberRRF1 >= $numberRRRFFF_or_FFF_2 && $numberRRF1 >=$numberRRF2) {
+	return (1,%badInFRSense0);
+    }
+
+    if ($numberRRRFFF_or_FFF_1 >=$numberRRF1 && $numberRRRFFF_or_FFF_1 >= $numberRRRFFF_or_FFF_2 && $numberRRRFFF_or_FFF_1 >= $numberRRF2) {
+	return (1,%badInFRSense1);
+    }
+
+    if ($numberRRRFFF_or_FFF_2 >= $numberRRF1 && $numberRRRFFF_or_FFF_2 >= $numberRRRFFF_or_FFF_1 && $numberRRRFFF_or_FFF_2 >=$numberRRF2) {
+	return (2,%badInFRSense3);
+    }
+
+    if ($numberRRF2 >=  $numberRRF1 && $numberRRF2 >= $numberRRRFFF_or_FFF_1 && $numberRRF2 >= $numberRRRFFF_or_FFF_2 ) {
+	return (2,%badInFRSense5);
+    }
+
+    #should not get here:
+    print STDERR "Error in findBadInFRSenseSOLiDSolexa()!\n";
+    return (1,%badInFRSense1);
+}
+
+sub getOrderedStrands {
+    my ($strand,$order) = ($_[0],$_[1]);
+    my @arr;
+    for my $i (0..scalar(@{$strand})-1) {
+	push(@arr,$strand->[$order->[$i]-1]);
+    }
+    return @arr;
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub checkClusters {
+
+    my ($ifRenv,$coord_start_chr1_cluster1,$coord_start_chr1_cluster2,$coord_start_chr2_cluster1,$coord_start_chr2_cluster2) = @_;
+    if ($ifRenv eq 'REVERSE_SENSE') {
+	if ($coord_start_chr1_cluster1 <= $coord_start_chr1_cluster2) {
+	    return ($coord_start_chr2_cluster1 <= $coord_start_chr2_cluster2)?1:0;
+	}
+	return ($coord_start_chr2_cluster1 >= $coord_start_chr2_cluster2)?1:0;
+    }
+    #if NORM
+    if ($coord_start_chr1_cluster1 <= $coord_start_chr1_cluster2) {
+        return ($coord_start_chr2_cluster1 >= $coord_start_chr2_cluster2)?1:0;
+    }
+    return ($coord_start_chr2_cluster1 <= $coord_start_chr2_cluster2)?1:0;
+}
+
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub translateSolidToRF {
+    my ($strandArr,$ends_orderArr)=@_;
+    my @array;
+    for my $i (0..scalar(@{$strandArr})-1) {
+	if ($ends_orderArr->[$i]==1 && $strandArr->[$i] eq 'F') {
+	    push(@array,'F');
+	}
+	if ($ends_orderArr->[$i]==2 && $strandArr->[$i] eq 'F') {
+	    push(@array,'R');
+	}
+	if ($ends_orderArr->[$i]==1 && $strandArr->[$i] eq 'R') {
+	    push(@array,'R');
+	}
+	if ($ends_orderArr->[$i]==2 && $strandArr->[$i] eq 'R') {
+	    push(@array,'F');
+	}
+    }
+    return @array;
+}
+
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+#convert the links file to the circos format
+sub links2segdup{
+
+    my($id,$color_code,$links_file,$segdup_file)=@_;
+
+    print LOG "# Converting to the circos format...\n";
+
+    tie (my %hcolor,'Tie::IxHash');						#color-code hash table
+    foreach my $col (keys %{$color_code}){
+	my ($min_links,$max_links)=split(",",$color_code->{$col});
+	$hcolor{$col}=[$min_links,$max_links];
+    }
+
+    open LINKS, "<$links_file" or die "$0: can't open $links_file :$!\n";
+    open SEGDUP, ">$segdup_file" or die "$0: can't write in the output: $segdup_file :$!\n";
+
+    my $index=1;
+    while(<LINKS>){
+
+	my ($chr1,$start1,$end1,$chr2,$start2,$end2,$count)=(split)[0,1,2,3,4,5,6];
+
+	my $color=getColor($count,\%hcolor,"circos");				#get the color-code according the number of links
+
+	print SEGDUP "$index\t$id$chr1\t$start1\t$end1\tcolor=$color\n".	#circos output
+		     "$index\t$id$chr2\t$start2\t$end2\tcolor=$color\n";
+	$index++;
+    }
+
+    close LINKS;
+    close SEGDUP;
+    print LOG "-- output created: $segdup_file\n";
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+#convert the links file to the bedPE format for BEDTools usage
+sub links2bedPElinksfile{
+
+    my ($sample,$links_file,$bedpe_file)=@_;
+
+    open LINKS, "<$links_file" or die "$0: can't open $links_file :$!\n";
+    open BEDPE, ">$bedpe_file" or die "$0: can't write in the output: $bedpe_file :$!\n";
+
+    my $nb_links=1;
+
+    while(<LINKS>){
+
+	chomp;
+	my @t=split("\t",$_);
+	my ($chr1,$start1,$end1,$chr2,$start2,$end2)=splice(@t,0,6);
+	my $type=($chr1 eq $chr2)? "INTRA":"INTER";
+	$type.="_".$t[10];
+
+	$start1--; $start2--;
+
+	print BEDPE "$chr1\t$start1\t$end1\t$chr2\t$start2\t$end2".
+	"\t$sample"."_link$nb_links\t$type\t.\t.".
+	"\t".join("|",@t)."\n";
+
+	$nb_links++;
+    }
+
+    close LINKS;
+    close BEDPE;
+
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub bedPElinks2linksfile{
+
+    my ($bedpe_file,$links_file)=@_;
+
+    open BEDPE, "<$bedpe_file" or die "$0: can't open: $bedpe_file :$!\n";
+    open LINKS, ">$links_file" or die "$0: can't write in the output $links_file :$!\n";
+
+    while(<BEDPE>){
+
+	chomp;
+	my $sample=(split("_",(split("\t",$_))[6]))[0];
+	my @t1=(split("\t",$_))[0,1,2,3,4,5];
+	my @t2=split(/\|/,(split("\t",$_))[10]);
+	push(@t2,$sample);
+
+	print LINKS join("\t",@t1)."\t".join("\t",@t2)."\n";
+
+    }
+    close BEDPE;
+    close LINKS;
+
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+#convert the links file to the bed format
+sub links2bedfile{
+
+    my ($tag_length,$color_code,$links_file,$bed_file)=@_;
+
+    print LOG "# Converting to the bed format...\n";
+
+    my $compare=1;
+    if($links_file!~/compared$/){
+	$compare=0;
+	$tag_length->{none}->{1}=$tag_length->{1};
+	$tag_length->{none}->{2}=$tag_length->{2};
+    }
+
+    #color-code hash table
+    tie (my %hcolor,'Tie::IxHash');
+    my %color_order;
+    my $n=1;
+    foreach my $col (keys %{$color_code}){
+	my ($min_links,$max_links)=split(",",$color_code->{$col});
+	$hcolor{$col}=[$min_links,$max_links];
+	$color_order{$col}=$n;
+	$n++;
+    }
+
+    my %pair;
+    my %pt;
+    $n=1;
+    open LINKS, "<$links_file" or die "$0: can't open $links_file:$!\n";
+
+    my %str=( "F"=>"+", "R"=>"-" );
+
+    while(<LINKS>){
+
+	my @t=split;
+	my $sample=($compare)? pop(@t):"none";
+
+	my $chr1=$t[0];
+	my $chr2=$t[3];
+	$chr1 = "chr".$chr1 unless ($chr1 =~ m/chr/i);
+	$chr2 = "chr".$chr2 unless ($chr2 =~ m/chr/i);
+	my $same_chr=($chr1 eq $chr2)? 1:0;
+
+	my $count=$t[6];
+	my $color=getColor($count,\%hcolor,"bed");
+
+	my @pairs=deleteBadOrderSensePairs(split(",",$t[7]));
+	my @strand1=deleteBadOrderSensePairs(split(",",$t[8]));
+	my @strand2=deleteBadOrderSensePairs(split(",",$t[9]));
+	my @ends_order1=deleteBadOrderSensePairs(split(",",$t[10]));
+	my @ends_order2=deleteBadOrderSensePairs(split(",",$t[11]));
+	my @position1=deleteBadOrderSensePairs(split(",",$t[14]));
+	my @position2=deleteBadOrderSensePairs(split(",",$t[15]));
+	my @start1; my @end1; getCoordswithLeftMost(\@start1,\@end1,\@position1,\@strand1,\@ends_order1,$tag_length->{$sample});
+	my @start2; my @end2; getCoordswithLeftMost(\@start2,\@end2,\@position2,\@strand2,\@ends_order2,$tag_length->{$sample});
+
+
+	for my $p (0..$#pairs){
+
+	    if (!exists $pair{$pairs[$p]}){
+
+		if($same_chr){
+
+		    $pair{$pairs[$p]}->{0}=[ $chr1, $start1[$p]-1, $end2[$p], $pairs[$p], 0, $str{$strand1[$p]},
+				    $start1[$p]-1, $end2[$p], $color,
+				    2, $tag_length->{$sample}->{$ends_order1[$p]}.",".$tag_length->{$sample}->{$ends_order2[$p]}, "0,".($start2[$p]-$start1[$p]) ];
+		    $pt{$n}=$pair{$pairs[$p]}->{0};
+		    $n++;
+
+		}else{
+
+		    $pair{$pairs[$p]}->{1}=[ $chr1, $start1[$p]-1, $end1[$p] , $pairs[$p]."/1", 0, $str{$strand1[$p]},
+						$start1[$p]-1, $end1[$p], $color,
+						1, $tag_length->{$sample}->{$ends_order1[$p]}, 0];
+		    $pt{$n}=$pair{$pairs[$p]}->{1};
+		    $n++;
+
+
+		    $pair{$pairs[$p]}->{2}=[ $chr2, $start2[$p]-1, $end2[$p], $pairs[$p]."/2", 0, $str{$strand2[$p]},
+						$start2[$p]-1, $end2[$p], $color,
+						1, $tag_length->{$sample}->{$ends_order2[$p]}, 0];
+		    $pt{$n}=$pair{$pairs[$p]}->{2};
+		    $n++;
+		}
+	    }else{
+
+		if($same_chr){
+		    ${$pair{$pairs[$p]}->{0}}[8]=$color if($color_order{$color}>$color_order{${$pair{$pairs[$p]}->{0}}[8]});
+		}else{
+		    ${$pair{$pairs[$p]}->{1}}[8]=$color if($color_order{$color}>$color_order{${$pair{$pairs[$p]}->{1}}[8]});
+		    ${$pair{$pairs[$p]}->{2}}[8]=$color if($color_order{$color}>$color_order{${$pair{$pairs[$p]}->{2}}[8]});
+		}
+	    }
+	}
+    }
+    close LINKS;
+
+    my $nb_pairs=$n-1;
+
+    open BED, ">$bed_file" or die "$0: can't write in the output: $bed_file :$!\n";
+    print BED "track name=\"$bed_file\" description=\"mate pairs involved in links\" ".
+	      "visibility=2 itemRgb=\"On\"\n";
+
+    for my $i (1..$nb_pairs){
+	print BED join("\t",@{$pt{$i}})."\n";
+    }
+
+    close BED;
+
+    print LOG "-- output created: $bed_file\n";
+
+    undef %pair;
+    undef %pt;
+
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub deleteBadOrderSensePairs{
+
+    my (@tab)=@_;
+    my @tab2;
+
+    foreach my $v (@tab){
+
+	$v=~s/[\(\)]//g;
+	push(@tab2,$v) if($v!~/[\$\*\@]$/);
+    }
+    return @tab2;
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub getAllEntries{
+    my (@tab)=split (/,/,$_[0]);
+    my @tab2;
+
+    foreach my $v (@tab){
+
+	$v=~s/[\(\)]//g;
+	push(@tab2,$v);
+    }
+    return @tab2;
+}#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub getAllEntriesWOspecialChar{
+    my (@tab)=split (/,/,$_[0]);
+    my @tab2;
+
+    foreach my $v (@tab){
+
+	$v=~s/[\(\)\$\*\@]//g;
+	push(@tab2,$v);
+    }
+    return @tab2;
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub links2SVfile{
+
+    my($links_file,$sv_file)=@_;
+
+    print LOG "# Converting to the sv output table...\n";
+    open LINKS, "<$links_file" or die "$0: can't open $links_file:$!\n";
+    open SV, ">$sv_file" or die "$0: can't write in the output: $sv_file :$!\n";
+
+    my @header=qw(chr_type SV_type BAL_type chromosome1 start1-end1 average_dist
+    chromosome2 start2-end2 nb_pairs score_strand_filtering score_order_filtering score_insert_size_filtering
+    final_score breakpoint1_start1-end1 breakpoint2_start2-end2);
+
+    my $nb_links=0;
+
+    while (<LINKS>){
+
+	my @t=split;
+	my @sv=();
+	my $sv_type="-";
+	my $strand_ratio="-";
+	my $eq_ratio="-";
+	my $eq_type="-";
+	my $insert_ratio="-";
+	my $link="-";
+	my ($bk1, $bk2)=("-","-");
+	my $score="-";
+
+	my ($chr1,$start1,$end1)=($t[0],$t[1],$t[2]);
+	my ($chr2,$start2,$end2)=($t[3],$t[4],$t[5]);
+	my $nb_pairs=$t[6];
+	$chr1 = "chr".$chr1 unless ($chr1 =~ m/chr/i);
+	$chr2 = "chr".$chr2 unless ($chr2 =~ m/chr/i);
+	my $chr_type=($chr1 eq $chr2)? "INTRA":"INTER";
+
+	#if strand filtering
+	if (defined $t[16]){
+	    #if inter-chr link
+	    $sv_type=$t[16];
+	    if(defined $t[17] && $t[17]=~/^(\d+)\/(\d+)$/){
+		$strand_ratio=floor($1/$2*100)."%";
+		$score=$t[18];
+	    }
+	    if(defined $t[18] && $t[18]=~/^(\d+)\/(\d+)$/){
+	    #if intra-chr link with insert size filtering
+		$strand_ratio=floor($1/$2*100)."%";
+		$link=floor($t[17]);
+		if($sv_type!~/^INV/){
+		    $insert_ratio=floor($1/$2*100)."%" if($t[19]=~/^(\d+)\/(\d+)$/);
+		    $score=$t[20];
+		}else{
+		    $score=$t[19];
+		}
+	    }
+	}
+
+	if(defined $t[18] && ($t[18] eq "UNBAL" || $t[18] eq "BAL")){
+
+	    #if strand and order filtering only and/or interchr link
+	    $eq_type=$t[18];
+	    $eq_ratio=floor($1/$2*100)."%" if($t[19]=~/^(\d+)\/(\d+)$/);
+	    ($bk1, $bk2)=($t[20],$t[21]);
+	    foreach my $bk ($bk1, $bk2){$bk=~s/\),\(/ /g; $bk=~s/(\(|\))//g; $bk=~s/,/-/g;}
+	    $score=$t[22];
+
+	}elsif(defined $t[19] && ($t[19] eq "UNBAL" || $t[19] eq "BAL")){
+
+	    #if all three filtering
+	    $link=floor($t[17]);
+	    $eq_type=$t[19];
+	    $eq_ratio=floor($1/$2*100)."%" if($t[20]=~/^(\d+)\/(\d+)$/);
+
+	    if(defined $t[21] && $t[21]=~/^(\d+)\/(\d+)$/){
+		$insert_ratio=floor($1/$2*100)."%";
+		($bk1, $bk2)=($t[22],$t[23]);
+		$score=$t[24];
+
+	    }else{
+		($bk1, $bk2)=($t[21],$t[22]);
+		$score=$t[23];
+	    }
+	    foreach my $bk ($bk1, $bk2){$bk=~s/\),\(/ /g; $bk=~s/(\(|\))//g; $bk=~s/,/-/g;}
+
+	}
+
+
+	push(@sv, $chr_type, $sv_type,$eq_type);
+	push(@sv,"$chr1\t$start1-$end1");
+	push(@sv, $link);
+	push(@sv,"$chr2\t$start2-$end2",
+	     $nb_pairs,$strand_ratio,$eq_ratio,$insert_ratio, decimal($score,4), $bk1, $bk2);
+
+
+	print SV join("\t",@sv)."\n";
+    }
+
+    close LINKS;
+    close SV;
+
+    system "sort  -k 9,9nr -k 13,13nr $sv_file > $sv_file.sorted";
+
+    open SV, "<".$sv_file.".sorted" or die "$0: can't open in the output: $sv_file".".sorted :$!\n";
+    my @links=<SV>;
+    close SV;
+
+    open SV, ">$sv_file" or die "$0: can't write in the output: $sv_file :$!\n";
+
+    print SV join("\t",@header)."\n";
+    print SV @links;
+    close SV;
+
+    unlink($sv_file.".sorted");
+
+    print LOG "-- output created: $sv_file\n";
+
+}
+#------------------------------------------------------------------------------#
+sub densityCalculation{
+
+    my ($chr,$chrID,$file,$tag_length,$window_dist,$step,$mates_file,$mates_file_ref,$density_file,$input_format)=@_;
+
+    my @sfile=split(/\./,$$mates_file[$file]);
+    my $fchr=$sfile[$#sfile];
+
+    my $fh = new FileHandle;
+
+    my %density;
+    my %density_ref;
+    my @ratio;
+    my ($cov,$cov_ref);
+
+    #FREQUENCY CALCULATION PROCEDURE
+    print LOG "# $fchr : Frequency calculation procedure...\n";
+    &FreqCalculation(\%density,$chr,$chrID,$tag_length,$window_dist,$step,$$mates_file[$file],$input_format);
+    &FreqCalculation(\%density_ref,$chr,$chrID,$tag_length,$window_dist,$step,$$mates_file_ref[$file],$input_format);
+
+    #MAKING RATIO AND OUTPUT
+    print LOG "\# Ratio calculation procedure...\n";
+    $density_file=~s/\/mates\//\/density\//;
+    $fh->open(">".$density_file) or die "$0: can't write in the output ".$density_file." :$!\n";
+
+    foreach my $k (1..$chr->{nb_chrs}){
+	foreach my $frag (1..$chr->{$k}->{nb_frag}){
+
+	    @ratio= ($chr->{$k}->{name},
+		    (${$chr->{$k}->{$frag}}[0]+1),
+		    (${$chr->{$k}->{$frag}}[1]+1));
+
+	    $cov=(exists $density{$k}{$frag}->{count})? $density{$k}{$frag}->{count}:0;
+	    $cov_ref=(exists $density_ref{$k}{$frag}->{count})? $density_ref{$k}{$frag}->{count}:0;
+
+	    push(@ratio,$cov,$cov_ref);
+	    push(@ratio,log($cov/$cov_ref)) if($cov && $cov_ref);
+	    push(@ratio,-log($cov_ref+1)) if(!$cov && $cov_ref);
+	    push(@ratio,log($cov+1)) if($cov && !$cov_ref);
+	    next if(!$cov && !$cov_ref);
+
+	    print $fh join("\t",@ratio)."\n";
+	}
+    }
+
+    $fh->close;
+    print LOG "-- output created: $density_file\n";
+
+    undef %density;
+    undef %density_ref;
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub FreqCalculation{
+
+    my ($density,$chr,$chrID,$tag_length,$window_dist,$step,$mates_file,$input_format) = @_;
+
+    my @sfile=split(/\./,$mates_file);
+    my $fchr=$sfile[$#sfile];
+    my $fh = new FileHandle;
+
+    my $nb_windows=0;
+    my $warn=100000;
+    my $record=0;
+    my %pair;
+
+    my ($sumX,$sumX2) = (0,0);
+
+    print LOG "\# Frequency calculation for $mates_file...\n";
+
+     if ($mates_file =~ /.gz$/) {
+	$fh->open("gunzip -c $mates_file |") or die "$0: can't open ".$mates_file.":$!\n";
+    }elsif($mates_file =~ /.bam$/){
+	o$fh->open("$SAMTOOLS_BIN_DIR/samtools view $mates_file |") or die "$0: can't open ".$mates_file.":$!\n";#GALAXY
+    }else{
+	$fh->open("<".$mates_file) or die "$0: can't open ".$mates_file.":$!\n";
+    }
+
+    while(<$fh>){
+
+	my @t=split;
+	my $mate=$t[0];
+
+	my ($chr_read1, $chr_read2, $firstbase_read1, $firstbase_read2, $end_order_read1, $end_order_read2,);
+
+	next if(exists $pair{$mate});
+
+	next if (!readMateFile(\$chr_read1, \$chr_read2, \$firstbase_read1, \$firstbase_read2,\$end_order_read1, \$end_order_read2, \@t, $input_format,$tag_length));
+
+	next unless (exists $chrID->{$chr_read1} || exists $chrID->{$chr_read2});
+	($chr_read1, $chr_read2)= ($chrID->{$chr_read1},$chrID->{$chr_read2});
+
+	$pair{$mate}=undef;
+	$record++;
+
+	my ($coord_start_read1,$coord_end_read1, $coord_start_read2,$coord_end_read2);
+
+	recupCoords($firstbase_read1,\$coord_start_read1,\$coord_end_read1,$tag_length->{$end_order_read1},$input_format);
+	recupCoords($firstbase_read2,\$coord_start_read2,\$coord_end_read2,$tag_length->{$end_order_read2},$input_format);
+
+	my $length = abs($coord_start_read1-$coord_start_read2);
+	$sumX += $length;							#add to sum and sum^2 for mean and variance calculation
+	$sumX2 += $length*$length;
+
+	for(my $i=1;$i<=$chr->{$chr_read1}->{'nb_frag'};$i++){
+
+	    if (abs ($coord_start_read1-${$chr->{$chr_read1}->{$i}}[0]) <= $window_dist){
+
+		&addToDensity($density,$chr_read1,$i,\$nb_windows)
+		if(overlap($coord_start_read1,$coord_end_read2,${$chr->{$chr_read1}->{$i}}[0],${$chr->{$chr_read1}->{$i}}[1]));
+
+	    }else{
+
+		$i=getNextFrag($coord_start_read1,$i,${$chr->{$chr_read1}->{$i}}[0],$chr->{$chr_read1}->{nb_frag},$window_dist,$step);
+	    }
+	}
+
+	if($record>=$warn){
+	    print LOG "-- $warn mate-pairs analysed - $nb_windows points created\n";
+	    $warn+=100000;
+	}
+    }
+    $fh->close;
+
+    print LOG "-- $fchr : Total : $record mate-pairs analysed - $nb_windows points created\n";
+
+    if($record>0){
+
+	my $mu = $sumX/$record;
+	my $sigma = sqrt($sumX2/$record - $mu*$mu);
+	print LOG "-- $fchr : mu length = $mu, sigma length = $sigma\n";
+    }
+
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub ratio2segdup{
+
+    my($id,$density_file,$segdup_file)=@_;
+
+    print LOG "# Converting to circos format...\n";
+
+    open RATIO, "<$density_file" or die "$0: can't open $density_file :$!\n";
+    open SEGDUP, ">$segdup_file" or die "$0: can't write in the output: $segdup_file :$!\n";
+
+    while(<RATIO>){
+	chomp;
+	my ($chr1,$start1,$end1,$ratio)=(split /\t/)[0,1,2,5];
+	print SEGDUP "$id$chr1\t$start1\t$end1\t$ratio\n";
+    }
+
+    close RATIO;
+    close SEGDUP;
+    print LOG "-- output created: $segdup_file\n";
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub ratio2bedfile{
+
+    my($density_file,$bed_file)=@_;
+
+    print LOG "# Converting to bedGraph format...\n";
+
+    open RATIO, "<$density_file" or die "$0: can't open $density_file :$!\n";
+    open BED, ">$bed_file" or die "$0: can't write in the output: $bed_file :$!\n";
+    print BED "track type=bedGraph name=\"$bed_file\" description=\"log ratios for cnv detection\" ".
+	      "visibility=2 color=255,0,0 alwaysZero=\"On\"\n";
+
+    while(<RATIO>){
+	chomp;
+	my ($chr1,$start1,$end1,$ratio)=(split /\t/)[0,1,2,5];
+	$chr1 = "chr".$chr1 unless ($chr1 =~ m/chr/);
+	print BED "$chr1\t".($start1-1)."\t$end1\t$ratio\n";
+    }
+
+    close RATIO;
+    close BED;
+    print LOG "-- output created: $bed_file\n";
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub inverseSense{
+
+    my $mate_sense=$_[0];
+    my %reverse=( 'F'  => 'R' , 'R'  => 'F' ,
+		  'FF' => 'RR', 'RR' => 'FF',
+		  'FR' => 'RF', 'RF' => 'FR');
+    return $reverse{$mate_sense};
+}
+
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub getNextFrag{
+
+    my ($read_start,$frag_num,$frag_start,$frag_last,$window_dist,$step)=@_;
+
+    my $how_far = $read_start-$frag_start;
+    my $nb_windows_toskip;
+
+    if($how_far>0){
+
+	$nb_windows_toskip=($how_far/$step)-($window_dist/$step);
+	$nb_windows_toskip=~ s/\..*//;
+	$nb_windows_toskip=0 if($nb_windows_toskip<0);
+	return ($frag_num + $nb_windows_toskip);
+    }
+    return $frag_last;
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub getColor{
+
+    my($count,$hcolor,$format)=@_;
+    for my $col ( keys % { $hcolor} ) {
+       return $col if($count>=$hcolor->{$col}->[0] && $count<=$hcolor->{$col}->[1]);
+    }
+    return "white" if($format eq "circos");
+    return "255,255,255" if($format eq "bed");
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub recupCoords{
+
+    my($c_hit,$cs_hit,$ce_hit,$tag_length,$input_format)=@_;
+    my $strand = 'F';
+
+    if ($c_hit=~s/^\-//) {
+	$strand='R';
+	    $$cs_hit=$c_hit;
+	    $$ce_hit=$c_hit-($tag_length-1);
+    }else{
+	$$cs_hit=$c_hit;
+	$$ce_hit=$c_hit+($tag_length-1);
+    }
+    return $strand;
+
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub overlap {
+    my($cs_hit,$ce_hit,$cs_region,$ce_region)=@_;
+    if( (($cs_hit < $cs_region) && ($ce_hit < $cs_region )) || (($cs_hit > $ce_region) && ($ce_hit > $ce_region )) ) {
+        return 0;
+    }
+    return 1;
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub makeLink {
+
+    my ($link,$chr1,$frag1,$chr2,$frag2,$mt,$nb)=@_;
+
+    if($chr1>$chr2){
+        ($chr1,$chr2)= ($chr2,$chr1);
+        ($frag1,$frag2)= ($frag2,$frag1);
+    }
+
+     if($chr1 == $chr2){
+        if($frag1>$frag2){
+	    ($frag1,$frag2)= ($frag2,$frag1);
+	}
+    }
+
+    if(!exists $link->{$chr1}->{$chr2}->{$frag1}->{$frag2}){
+        $link->{$chr1}->{$chr2}->{$frag1}->{$frag2}=$mt;
+	$$nb++;
+    }elsif($link->{$chr1}->{$chr2}->{$frag1}->{$frag2}!~/(^|,)$mt(,|$)/){
+        $link->{$chr1}->{$chr2}->{$frag1}->{$frag2}.=",$mt";
+    }
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+#fonction of adding the read to the density profile
+sub addToDensity {
+
+    my ($density,$chr1,$frag1,$nb)=@_;
+
+    if(!exists $density->{$chr1}->{$frag1}->{count}){
+            $density->{$chr1}->{$frag1}->{count}=1;
+	    $$nb++;
+    }else{
+            $density->{$chr1}->{$frag1}->{count}++;
+    }
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub floor {
+    my $nb = $_[0];
+    $nb=~ s/\..*//;
+    return $nb;
+}
+#------------------------------------------------------------------------------#
+sub decimal{
+
+  my $num=shift;
+  my $digs_to_cut=shift;
+
+  $num=sprintf("%.".($digs_to_cut-1)."f", $num) if ($num=~/\d+\.(\d){$digs_to_cut,}/);
+
+  return $num;
+}
+
+#------------------------------------------------------------------------------#
+sub max {
+
+    my($max) = shift(@_);
+    foreach my $temp (@_) {
+        $max = $temp if $temp > $max;
+    }
+    return($max);
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub min {
+
+    my($min) = shift(@_);
+    foreach my $temp (@_) {
+        $min = $temp if $temp < $min;
+    }
+    return($min);
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub sortTablebyIndex{
+    my ($tab1,$tab2)=@_;
+    my @tab3;
+
+    foreach my $i (@$tab1){
+	$tab3[$i]=$$tab2[$$tab1[$i]];
+    }
+    return @tab3;
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub round {
+  my $number = shift || 0;
+  my $dec = 10 ** (shift || 0);
+  return int( $dec * $number + .5 * ($number <=> 0)) / $dec;
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub getUniqueTable{
+
+    my (@tab)=@_;
+    my (%saw,@out)=();
+    undef %saw;
+    return sort(grep(!$saw{$_}++, @tab));
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+sub catFiles {
+
+    unlink("$_[1]") if(exists $_[1]);
+    system qq( cat "$_" >> "$_[1]" ) for @{$_[0]};
+}
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+#check if the configuration file is correct
+sub validateconfiguration{
+
+    my %conf=%{$_[0]};
+    my $list_prgs="@ARGV";
+
+    my @general_params=qw(input_format mates_orientation read1_length read2_length mates_file cmap_file);
+    my @detection_params=qw(split_mate_file window_size step_length split_mate_file);
+    my @filtering_params=qw(split_link_file nb_pairs_threshold strand_filtering split_link_file);
+    my @circos_params=qw(organism_id colorcode);
+    my @bed_params=qw(colorcode);
+    my @compare_params=qw(list_samples file_suffix);
+
+    foreach my $dir ($conf{general}{output_dir},$conf{general}{tmp_dir}){
+
+	unless (defined($dir)) {
+	    $dir = ".";
+	}
+	unless (-d $dir){
+	    mkdir $dir or die;
+	}
+	$dir.="/" if($dir!~/\/$/);
+    }
+
+    unless (defined($conf{general}{num_threads})) {
+	    $conf{general}{num_threads} = 1;
+	}
+    $conf{general}{num_threads}=24 if($conf{general}{num_threads}>24);
+
+    if($list_prgs!~/links2compare/){
+
+	foreach my $p (@general_params){
+	    die("Error Config : The parameter \"$p\" is not defined\n") if (!defined $conf{general}{$p});
+	}
+
+	$conf{general}{input_format}="sam" if($conf{general}{input_format} eq "bam");
+
+	unless (defined($conf{general}{sv_type})) {
+	    $conf{general}{sv_type} = "all";
+	}
+
+	$conf{general}{read_lengths}={ 1=> $conf{general}{read1_length}, 2=> $conf{general}{read2_length}};
+    }
+
+    if($list_prgs=~/(linking|cnv)/){
+
+	foreach my $p (@detection_params){
+	    die("Error Config : The parameter \"$p\" is not defined\n") if (!defined $conf{detection}{$p});
+	}
+
+	die("Error Config : The parameter \"mates_file_ref\" is not defined\n") if($list_prgs=~/cnv/ && !defined $conf{detection}{mates_file_ref});
+
+	if($conf{detection}{step_length}>$conf{detection}{window_size}){
+	    die("Error Config : Parameter \"step_length\" should not exceed \"window size\"\n");
+	}
+
+	unless (-d $conf{general}{tmp_dir}."/mates"){
+	    mkdir $conf{general}{tmp_dir}."/mates" or die;
+	}
+
+	if($list_prgs=~/linking/){
+	    unless (-d $conf{general}{tmp_dir}."/links"){
+		mkdir $conf{general}{tmp_dir}."/links" or die;
+	    }
+	}
+	if($list_prgs=~/cnv/){
+	    unless (-d $conf{general}{tmp_dir}."/density"){
+		mkdir $conf{general}{tmp_dir}."/density" or die;
+	    }
+	}
+
+    }
+
+    if($list_prgs=~/filtering/){
+
+	foreach my $p (@filtering_params) {
+	    die("Error Config : The filtering parameter \"$p\" is not defined\n") if (!defined $conf{filtering}{$p});
+
+	}
+
+	if(defined($conf{filtering}{chromosomes})) {
+	    my @chrs=split(",",$conf{filtering}{chromosomes});
+	    my $exclude=($chrs[0]=~/^\-/)? 1:0;
+	    for my $chrName (@chrs){
+
+		die("Error Config : The filtering parameter \"chromosomes\" is not valid\n")
+		if(($chrName!~/^\-/ && $exclude) || ($chrName=~/^\-/ && !$exclude));
+
+	    }
+	}
+
+	if (( $conf{filtering}{order_filtering} )&& !$conf{filtering}{strand_filtering}) {
+	    die("Error Config : The parameter strand_filtering is set to \"0\" while order_filtering is selected".
+		"\nChange strand_filtering to \"1\" if you want to use the order filtering\n");
+	}
+	if (( !defined($conf{filtering}{mu_length}) || !defined($conf{filtering}{sigma_length}) )&& $conf{filtering}{order_filtering}) {
+	    die("Error Config : You should set parameters \"mu_length\" and \"sigma_length\" to use order filtering\n");
+	}
+	if (( $conf{filtering}{insert_size_filtering} )&& !$conf{filtering}{strand_filtering}) {
+	    die("Error Config : The parameter strand_filtering is set to \"0\" while insert_size_filtering is selected".
+		"\nChange strand_filtering to \"1\" if you want to use the insert size filtering\n");
+	}
+	if (( !defined($conf{filtering}{mu_length}) || !defined($conf{filtering}{sigma_length}) )&& $conf{filtering}{insert_size_filtering}) {
+	    die("Error Config : You should set parameters \"mu_length\" and \"sigma_length\" to use discriminate insertions from deletions\n");
+	}
+
+	if (!defined($conf{filtering}{indel_sigma_threshold})) {
+	    $conf{filtering}{indel_sigma_threshold} = 2;
+	}
+	if (!defined($conf{filtering}{dup_sigma_threshold})) {
+	    $conf{filtering}{dup_sigma_threshold} = 2;
+	}
+	if (!defined($conf{filtering}{singleton_sigma_threshold})) {
+	    $conf{filtering}{singleton_sigma_threshold} = 4;
+	}
+
+	if (!defined($conf{filtering}{nb_pairs_order_threshold})) {
+	    $conf{filtering}{nb_pairs_order_threshold} = 1;
+	}
+
+	if (!defined($conf{filtering}{final_score_threshold})) {
+	    $conf{filtering}{final_score_threshold} = 0.8;
+	}
+
+	if ($conf{filtering}{nb_pairs_order_threshold}>$conf{filtering}{nb_pairs_threshold}) {
+	    die("Error Config : Parameter \"nb_pairs_order_threshold\" should not exceed \"nb_pairs_threshold\"\n");
+	}
+
+    }
+
+    if($list_prgs=~/2circos$/){
+	foreach my $p (@circos_params) {
+	    next if($list_prgs=~/^ratio/ && $p eq "colorcode");
+	    die("Error Config : The circos parameter \"$p\" is not defined\n") if (!defined $conf{circos}{$p});
+	}
+    }
+
+    if($list_prgs=~/2bed$/){
+	foreach my $p (@bed_params) {
+	    die("Error Config : The bed parameter \"$p\" is not defined\n") if (!defined $conf{bed}{$p});
+	}
+    }
+
+    if($list_prgs=~/links2compare/){
+	foreach my $p (@compare_params) {
+	    die("Error Config : The compare parameter \"$p\" is not defined\n") if (!defined $conf{compare}{$p});
+	}
+
+	unless (defined($conf{compare}{same_sv_type})) {
+	    $conf{compare}{same_sv_type} = 0;
+	}
+
+	unless (defined($conf{compare}{min_overlap})) {
+	    $conf{compare}{min_overlap} = 1E-9;
+	}
+
+	if($conf{compare}{circos_output}){
+	    foreach my $p (@circos_params) {
+		next if($list_prgs=~/^ratio/ && $p eq "colorcode");
+		die("Error Config : The circos parameter \"$p\" is not defined\n") if (!defined $conf{circos}{$p});
+	    }
+	}
+	if($conf{compare}{bed_output}){
+	    foreach my $p (@bed_params) {
+		die("Error Config : The bed parameter \"$p\" is not defined\n") if (!defined $conf{bed}{$p});
+	    }
+	    die("Error Config : The compare parameter \"list_read_lengths\" is not defined\n") if (!defined $conf{compare}{list_read_lengths});
+
+	    my @samples=split(",",$conf{compare}{list_samples});
+	    my @read_lengths=split(",",$conf{compare}{list_read_lengths});
+	    for my $i (0..$#samples){
+		my @l=split("-",$read_lengths[$i]);
+		$conf{compare}{read_lengths}{$samples[$i]}={ 1=> $l[0], 2=> $l[1]};
+	    }
+	}
+    }
+
+
+}
+#------------------------------------------------------------------------------#
+#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::#
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/svdetect/SVDetect_run_parallel.xml	Thu Jul 12 12:10:27 2012 -0400
@@ -0,0 +1,324 @@
+<tool id="svdetect_run_parallel" name="Detect clusters of anomalously mapped pairs">
+
+<description>and identify structural variants</description>
+
+<command interpreter="perl">SVDetect_run_parallel.pl
+
+#if $getLinks.linking == "linking"
+linking
+<!-- -out1 '$links_file' -->
+#end if
+#if $getFilteredLinks.filtering == "filtering"
+filtering
+<!--- out2 '$flinks_file' -->
+#if str($getFilteredLinks.links2SV) == "create"
+links2SV
+-out3 '$sv_file'
+#end if
+#if  $getFilteredLinks.file_conversion.file_conversion_select=="convert" and str($getFilteredLinks.file_conversion.links2circos) == "create"
+links2circos
+-out4 '$circos_file'
+#end if
+#if  $getFilteredLinks.file_conversion.file_conversion_select=="convert" and str($getFilteredLinks.file_conversion.links2bed) == "create"
+links2bed
+-out5 '$bed_file'
+#end if
+#end if
+-conf '$config_file'
+-l '$log_file'
+-N '$sample_name'
+
+</command>
+
+<inputs>
+	<param name="sample_name" type="text" value="sample" label="Sample Name"/>
+	<param name="mates_file" format="bam" type="data" label="Input BAM file (.ab.bam)"/>
+   	<param name="cmap_file" format="len" type="data" label="Chromosomes list file (.len)" help="Tabulated file format with Chromosome ID (integer from 1), name and length"/>
+  	<param name="mates_orientation" type="select" format="txt" label="Type of sequencing technology and libraries">
+		<option value="FR">Illumina paired-ends</option>
+		<option value="RF">Illumina mate-pairs</option>
+		<option value="FR">SOLiD paired-ends</option>
+		<option value="RR">SOLiD mate-pairs</option>
+   	</param>
+	<param name="read1_length" type="integer" size="10" value="50" label="Read 1 length (bp)" help="Length of the first read in a pair (left read)"/>
+	<param name="read2_length" type="integer" size="10" value="50" label="Read 2 length (bp)" help="Length of the second read in a pair (right read)"/>
+	<param name="sv_type" type="select" format="txt" label="Type of SV to detect">
+		<option value="all">all types of SVs</option>
+		<option value="intra">intrachromosomal SVs only</option>
+		<option value="inter">interchromosomal SVs only</option>
+  	</param>
+
+   	<conditional name="getLinks">
+   		<param name="linking" type="select" label="Linking procedure" help="Detection and isolation of links">
+			<option value="linking">Yes</option>
+			<option value="">No, already done</option>
+  		</param>
+		<when value="">
+    			<!-- do nothing here -->
+   		</when>
+    		<when value="linking">
+			<param name="splitmate" label="Do you want to split the original mate file per chromosome for parallel computing?" type="boolean" truevalue="split" falsevalue="do_not_split" checked="True" help="Untick it if already done"/>
+    			<param name="window_size" type="integer" size="20" value="3000" label="Window size (bp)" help="Equal to at least “2µ+2√2σ"/>
+			<param name="step_length" type="integer" size="20" value="250" label="Step length size (bp)" help="Equal to 1/2 or 1/4 of the window size"/>
+    		</when>
+   	</conditional>
+
+   	<conditional name="getFilteredLinks">
+	   	<param name="filtering" type="select" label="Filtering procedure" help="Filtering of links according different parameters and thresholds">
+			<option value="filtering">Yes</option>
+                        <option value="">No</option>
+	  	</param>
+		<when value="">
+	    		<!-- do nothing here -->
+	   	</when>
+	    	<when value="filtering">
+
+			<param name="splitlink" label="Do you want to split the original link file per chromosome for parallel computing?" type="boolean" truevalue="split" falsevalue="do_not_split" checked="False" help="Untick it if (the linking is) already done"/>
+			<param name="chromosomes" type="text" size="20" label="List of chromosome names to keep or exclude"/>
+			<param name="nb_pairs_threshold" type="integer" size="20" value="5" label="Minimum number of pairs in a cluster"/>
+
+			<conditional name="filter1">
+	   			<param name="strand_filtering" type="select" label="Strand filtering procedure">
+					<option value="strand">Yes</option>
+					<option value="">No</option>
+	  			</param>
+				<when value="">
+	    				<!-- do nothing here -->
+	   			</when>
+	    			<when value="strand">
+
+					<conditional name="filter2">
+			   			<param name="order_filtering" type="select" label="Order filtering procedure">
+							<option value="order">Yes</option>
+							<option value="">No</option>
+			  			</param>
+						<when value="">
+			    				<!-- do nothing here -->
+			   			</when>
+			    			<when value="order">
+
+							<conditional name="filter3">
+					   			<param name="insert_size_filtering" type="select" label="Insert-size filtering procedure">
+									<option value="insert">Yes</option>
+									<option value="">No</option>
+					  			</param>
+								<when value="">
+					    				<!-- do nothing here -->
+					   			</when>
+					    			<when value="insert">
+									<param name="indel_sigma_threshold" type="float" size="20" value="3" label="Minimal number of sigma fold for the insert size filtering and to call insertions and deletions"/>
+									<param name="dup_sigma_threshold" type="float" size="20" value="3" label="minimal number of sigma fold for the insert size filtering to call tandem duplications"/>
+									<param name="singleton_sigma_threshold" type="float" size="20" value="4" label="Minimal number of sigma fold for the insert size filtering to call singletons" help="for Illumina mate-pairs only"/>
+			    					</when>
+			   				</conditional>
+
+							<param name="mu_length" type="integer" size="20" value="3000" label="Mean insert size value (µ) of normally mapped mate-pairs, in bp"/>
+							<param name="sigma_length" type="integer" size="20" value="250" label="Calculated sd value (σ) from the distribution of normally mapped  mate-pairs, in bp"/>
+			    				<param name="nb_pairs_order_threshold" type="integer" size="20" value="2" label="Minimal number of pairs in a subgroup of paired-end reads for balanced events"/>
+			    			</when>
+			   		</conditional>
+
+					<param name="final_score_threshold" type="float" size="20"  value="1.0" label="Minimal final filtering score for calling SVs" help="A value of 1 means all the pairs in a cluster were consistent between each other after applying filters"/>
+	    			</when>
+	   		</conditional>
+
+			<param name="links2SV" label="Do you want to have filtered links in a tabulated file format showing significant SVs?" type="boolean" truevalue="create" falsevalue="do_not_create" checked="True"/>
+
+			<conditional name="file_conversion">
+				<param name="file_conversion_select" type="select" label="Output file conversion" help="Converts filtered links to Circos/BED files format for graphical view of SVs">
+					<option value="do_not_convert">No</option>
+					<option value="convert">Yes</option>
+				</param>
+				<when value="do_not_convert">
+					    <!-- do nothing here -->
+				</when>
+				<when value="convert">
+					<param name="links2circos" label="Converts the link list to the Circos link format" type="boolean" truevalue="create" falsevalue="do_not_create" checked="True"/>
+					<param name="links2bed" label="Converts the link list to the UCSC BED format" type="boolean" truevalue="create" falsevalue="do_not_create" checked="False"/>
+					<param name="organism_id" type="text" size="10" value="hs" label="Organism ID"/>
+					<repeat name="color_code" title="Color-code" min="1" max="7">
+						<param name="color" type="select" label="Color">
+							<option value="grey">grey</option>
+							<option value="black">black</option>
+							<option value="blue">blue</option>
+							<option value="green">green</option>
+							<option value="purple">purple</option>
+							<option value="orange">orange</option>
+							<option value="red">red</option>
+						</param>
+						<param name="interval" type="text" value="1,3" label="Interval"/>
+					</repeat>
+				</when>
+			</conditional>
+    		</when>
+	</conditional>
+</inputs>
+
+
+<outputs>
+	<!--<data format="txt" name="links_file" label="svdetect.links">
+		<filter>getLinks['linking']=="linking"</filter>
+	</data>
+	<data format="txt" name="flinks_file" label="svdetect.links.filtered">
+		<filter>getFilteredLinks['filtering']=="filtering"</filter>
+	</data>-->
+	<data format="sv" name="sv_file" label="${sample_name}.sv">
+		<filter>(
+			getFilteredLinks['filtering']=="filtering" and
+			getFilteredLinks['links2SV'] is True
+			)
+	         </filter>
+	</data>
+	<data format="segdup" name="circos_file" label="${sample_name}.segdup">
+		<filter>(
+			getFilteredLinks['filtering']=="filtering" and
+			getFilteredLinks['file_conversion']['file_conversion_select']=="convert" and
+			getFilteredLinks['file_conversion']['links2circos'] is True
+			)
+	         </filter>
+	</data>
+	<data format="bed" name="bed_file" label="${sample_name}.bed">
+		<filter>(
+			getFilteredLinks['filtering']=="filtering" and
+			getFilteredLinks['file_conversion']['file_conversion_select']=="convert" and
+			getFilteredLinks['file_conversion']['links2bed'] is True
+			)
+	         </filter>
+	</data>
+	<data format="txt" name="log_file" label="${sample_name}.svdetect_run.log"/>
+</outputs>
+
+
+
+<configfiles>
+	<configfile name="config_file">
+&lt;general&gt;
+input_format = bam
+sv_type = ${sv_type}
+mates_orientation=${mates_orientation}
+read1_length=${read1_length}
+read2_length=${read2_length}
+mates_file=${mates_file}
+cmap_file=${cmap_file}
+tmp_dir=$__new_file_path__/svdetect/tmp
+output_dir=$__new_file_path__/svdetect
+num_threads=8
+&lt;/general&gt;
+
+#if $getLinks.linking == "linking"
+&lt;detection&gt;
+#if str($getLinks.splitmate) == "split"
+split_mate_file=1
+#else
+split_mate_file=0
+#end if
+window_size=${getLinks.window_size}
+step_length=${getLinks.step_length}
+&lt;/detection&gt;
+#end if
+
+#if $getFilteredLinks.filtering == "filtering"
+&lt;filtering&gt;
+#if str($getFilteredLinks.splitlink) == "split"
+split_link_file=1
+#else
+split_link_file=0
+#end if
+#if str($getFilteredLinks.chromosomes) != ""
+chromosomes=${getFilteredLinks.chromosomes}
+#end if
+nb_pairs_threshold=${getFilteredLinks.nb_pairs_threshold}
+#if $getFilteredLinks.filter1.strand_filtering == "strand"
+strand_filtering=1
+final_score_threshold=${getFilteredLinks.filter1.final_score_threshold}
+#if $getFilteredLinks.filter1.filter2.order_filtering == "order"
+order_filtering=1
+mu_length=${getFilteredLinks.filter1.filter2.mu_length}
+sigma_length=${getFilteredLinks.filter1.filter2.sigma_length}
+nb_pairs_order_threshold=${getFilteredLinks.filter1.filter2.nb_pairs_order_threshold}
+#if $getFilteredLinks.filter1.filter2.filter3.insert_size_filtering == "insert"
+insert_size_filtering=1
+indel_sigma_threshold=${getFilteredLinks.filter1.filter2.filter3.indel_sigma_threshold}
+dup_sigma_threshold=${getFilteredLinks.filter1.filter2.filter3.dup_sigma_threshold}
+singleton_sigma_threshold=${getFilteredLinks.filter1.filter2.filter3.singleton_sigma_threshold}
+#else
+insert_size_filtering=0
+#end if
+#else
+order_filtering=0
+#end if
+#else
+strand_filtering=0
+#end if
+&lt;/filtering&gt;
+#end if
+
+#if $getFilteredLinks.filtering == "filtering"
+#if $getFilteredLinks.file_conversion.file_conversion_select == "convert"
+#if str($getFilteredLinks.file_conversion.links2circos) == "create"
+&lt;circos&gt;
+organism_id=${getFilteredLinks.file_conversion.organism_id}
+&lt;colorcode&gt;
+#for $color_repeat in $getFilteredLinks.file_conversion.color_code
+${color_repeat.color}=${color_repeat.interval}
+#end for
+&lt;/colorcode&gt;
+&lt;/circos&gt;
+#end if
+#if str($getFilteredLinks.file_conversion.links2bed) == "create"
+&lt;bed&gt;
+&lt;colorcode&gt;
+#for $color_repeat in $getFilteredLinks.file_conversion.color_code
+#if str($color_repeat.color)== "grey"
+190,190,190=${color_repeat.interval}
+#end if
+#if str($color_repeat.color)== "black"
+0,0,0=${color_repeat.interval}
+#end if
+#if str($color_repeat.color)== "blue"
+0,0,255=${color_repeat.interval}
+#end if
+#if str($color_repeat.color)== "green"
+0,255,0=${color_repeat.interval}
+#end if
+#if str($color_repeat.color)== "purple"
+153,50,205=${color_repeat.interval}
+#end if
+#if str($color_repeat.color)== "orange"
+255,140,0=${color_repeat.interval}
+#end if
+#if str($color_repeat.color)== "red"
+255,0,0=${color_repeat.interval}
+#end if
+#end for
+&lt;/colorcode&gt;
+&lt;/bed&gt;
+#end if
+#end if
+#end if
+	</configfile>
+</configfiles>
+
+  <help>
+**What it does**
+
+SVDetect - Version : 0.8
+
+Parallel version (nCPU=8)
+
+SVDetect is a application for the isolation and the type prediction of intra- and inter-chromosomal rearrangements from paired-end/mate-pair sequencing data provided by the high-throughput sequencing technologies
+
+This tool aims to identifying structural variations (SVs) with both clustering and sliding-window strategies, and helping in their visualization at the genome scale.
+SVDetect is compatible with SOLiD and Illumina (>=1.3) reads.
+
+Manual documentation available at the http://svdetect.sourceforge.net/Site/Manual.html
+
+-----
+
+.. class:: infomark
+
+Contact Bruno Zeitouni (svdetect@curie.fr) for any questions or concerns about the Galaxy implementation of SVDetect.
+
+  </help>
+
+</tool>
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/svdetect/circos_graph.sh	Thu Jul 12 12:10:27 2012 -0400
@@ -0,0 +1,4 @@
+#Copy your circos directory installation in the galaxy svdetect folder
+#>> ls SVDetect/circos
+#bin/  Build.PL  CHANGES  data/  etc/  fonts/  install-unix*  lib/  MANIFEST.SKIP  README  svg/  TODO  tools/  tutorials/
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/svdetect/circos_graph.xml	Thu Jul 12 12:10:27 2012 -0400
@@ -0,0 +1,289 @@
+<tool id="circos_graph" name="Circos">
+
+<description>plots</description>
+
+<command interpreter="perl">circos
+
+-conf '$circos_config_file'
+-outputdir $__new_file_path__/svdetect
+-outputfile '$outputfile'
+-png
+-out1 '$out1_png'
+
+> '$log_file'
+
+</command>
+
+
+<inputs>
+	<param name="graph_name" type="text" size="20" value="graph1" label="Graph name"/>
+
+	<param name="karyotype" type="select" format="txt" label="Type of model organism">
+		<option value="data/karyotype.human_hg19.txt">Human (homo sapiens, hs) -hg19-</option>
+		<option value="data/karyotype.human.txt">Human (homo sapiens, hs) -hg18-</option>
+		<option value="data/2/karyotype.mouse.txt">Mouse (Mus Musculus, mm)</option>
+		<option value="data/2/karyotype.dog.txt">Dog (Canis familiaris, cf)</option>
+		<option value="data/2/karyotype.rt.txt">Rat (Rattus norvegicus, rn)</option>
+		<option value="data/karyotype.yeast.txt">Yeast (Saccharomyces Cerevisiae, sc) -SGD-</option>
+
+  	</param>
+	<param name="chromosomes_units" type="integer" size="50" value="1000000" label="Chromosomes units"/>
+	<param name="chromosomes" type="text" size="100" value="" label="List of chromosome names to keep or exclude" help="ex: hs2;hs3 or -hsX;-hsY">
+		<sanitizer>
+			<valid initial="string.printable">
+				<add value=";"/>
+			</valid>
+		</sanitizer>
+	</param>
+	<param name="link_file" format="segdup" type="data" label="Input link file (.segdup)"/>
+	<param name="outputfile" type="text" size="50" value="graph1.circos.png" label="Output image name file (png)"/>
+</inputs>
+
+<outputs>
+	<data format="txt" name="log_file" label="${graph_name}.circos.log"/>
+	<data format="png" name="out1_png" label="${outputfile}"/>
+</outputs>
+
+
+
+<configfiles>
+	<configfile name="ideogram_config_file">
+
+&lt;ideogram&gt;
+
+&lt;spacing&gt;
+
+default = 5u
+break   = 1u
+
+axis_break_at_edge = yes
+axis_break         = yes
+axis_break_style   = 2
+
+&lt;break_style 1&gt;
+stroke_color = black
+fill_color   = blue
+thickness    = 0.25r
+stroke_thickness = 2
+&lt;/break&gt;
+
+&lt;break_style 2&gt;
+stroke_color     = black
+stroke_thickness = 3p
+thickness        = 1.5r
+&lt;/break&gt;
+
+&lt;/spacing&gt;
+
+## thickness (px) of chromosome ideogram
+thickness        = 100p
+stroke_thickness = 2
+## ideogram border color
+stroke_color     = black
+fill             = yes
+## the default chromosome color is set here and any value
+## defined in the karyotype file overrides it
+fill_color       = black
+
+## fractional radius position of chromosome ideogram within image
+radius         = 0.85r
+show_label     = yes
+label_with_tag = yes
+label_font     = condensedbold
+label_radius   = dims(ideogram,radius) + 0.075r
+label_size     = 60p
+
+## cytogenetic bands
+band_stroke_thickness = 2
+
+## show_bands determines whether the outline of cytogenetic bands
+## will be seen
+show_bands            = yes
+## in order to fill the bands with the color defined in the karyotype
+## file you must set fill_bands
+fill_bands            = yes
+
+&lt;/ideogram&gt;
+
+	</configfile>
+
+	<configfile name="ticks_config_file">
+
+show_ticks          = yes
+show_tick_labels    = yes
+
+&lt;ticks&gt;
+radius               = dims(ideogram,radius_outer)
+multiplier           = 1e-6
+
+&lt;tick&gt;
+spacing        = 0.5u
+size           = 2p
+thickness      = 2p
+color          = grey
+show_label     = no
+label_size     = 12p
+label_offset   = 0p
+format         = %.2f
+&lt;/tick&gt;
+
+&lt;tick&gt;
+spacing        = 1u
+size           = 3p
+thickness      = 2p
+color          = dgrey
+show_label     = no
+label_size     = 12p
+label_offset   = 0p
+format         = %.2f
+&lt;/tick&gt;
+
+&lt;tick&gt;
+spacing        = 5u
+size           = 5p
+thickness      = 2p
+color          = black
+show_label     = yes
+label_size     = 16p
+label_offset   = 0p
+format         = %d
+&lt;/tick&gt;
+
+&lt;tick&gt;
+spacing        = 10u
+size           = 8p
+thickness      = 2p
+color          = black
+show_label     = yes
+label_size     = 20p
+label_offset   = 5p
+format         = %d
+&lt;/tick&gt;
+&lt;/ticks&gt;
+	</configfile>
+
+
+	<configfile name="circos_config_file">
+&lt;colors&gt;
+&lt;&lt;include etc/colors.conf&gt;&gt;
+&lt;/colors&gt;
+
+&lt;fonts&gt;
+&lt;&lt;include etc/fonts.conf&gt;&gt;
+&lt;/fonts&gt;
+
+&lt;&lt;include $ideogram_config_file&gt;&gt;
+&lt;&lt;include $ticks_config_file&gt;&gt;
+
+karyotype = $karyotype
+
+&lt;image&gt;
+24bit = yes
+##png = yes
+##svg = no
+## radius of inscribed circle in image
+radius         = 1500p
+background     = white
+## by default angle=0 is at 3 o'clock position
+angle_offset   = -90
+#angle_orientation = counterclockwise
+
+auto_alpha_colors = yes
+auto_alpha_steps  = 5
+&lt;/image&gt;
+
+chromosomes_units= $chromosomes_units
+
+#if str($chromosomes)==""
+chromosomes_display_default = yes
+#else
+chromosomes_display_default = no
+chromosomes = $chromosomes
+#end if
+
+&lt;links&gt;
+
+z      = 0
+radius = 0.95r
+bezier_radius = 0.2r
+
+&lt;link segdup&gt;
+show         = yes
+color        = dgrey_a5
+thickness    = 2
+file         = $link_file
+record_limit = 1000
+&lt;/link&gt;
+
+&lt;/links&gt;
+
+
+anglestep       = 0.5
+minslicestep    = 10
+beziersamples   = 40
+debug           = no
+warnings        = no
+imagemap        = no
+
+units_ok = bupr
+units_nounit = n
+	</configfile>
+</configfiles>
+
+  <help>
+**What it does**
+
+Circos
+
+Manual documentation available at the http://circos.ca/
+
+
+**Example of link segdup file**
+
+segdup file::
+
+    1	hs1	1077096	1078746	color=red
+    1	hs1	1080923	1082805	color=red
+    2	hs1	1137684	1137961	color=red
+    2	hs3	1138138	1138423	color=red
+    3	hs11	1169417	1170000	color=red
+    3	hs11	1170025	1170975	color=red
+    4	hs11	1222480	1224271	color=green
+    4	hs11	1223328	1225675	color=green
+    5	hs12	1223336	1225812	color=grey
+    5	hs13	1224709	1227633	color=grey
+    6	hs11	1223621	1226460	color=red
+    6	hs11	1224918	1227633	color=red
+    7	hs11	1399510	1401513	color=white
+    7	hs11	1401628	1403697	color=white
+    8	hs15	1652045	1653746	color=red
+    8	hs15	1657167	1658940	color=red
+    9	hs11	165333	165887	color=white
+    9	hs11	165981	168016	color=white
+    10	hs11	1702700	1702841	color=red
+    10	hs11	1702903	1703057	color=red
+    11	hs11	1912272	1915186	color=white
+    11	hs11	1937111	1939824	color=white
+    12	hs11	1983211	1983355	color=red
+    12	hs11	1983591	1983748	color=red
+    13	hs11	2913657	2913898	color=white
+    13	hs11	2914048	2914341	color=white
+    14	hs11	3090593	3090749	color=purple
+    14	hs11	3090709	3090864	color=purple
+    15	hs21	3466365	3466434	color=red
+    15	hs21	3466554	3466620	color=red
+    16	hsX	3603073	3603321	color=white
+    16	hsX	3603295	3603520	color=white
+
+
+
+-----
+
+.. class:: infomark
+
+Contact Bruno Zeitouni (svdetect@curie.fr) for any questions or concerns about the Galaxy implementation of Circos.
+
+
+  </help>
+
+</tool>