|
13
|
1 #!/usr/bin/perl -w
|
|
|
2
|
|
|
3 =pod
|
|
|
4
|
|
|
5 =head1 NAME
|
|
|
6
|
|
|
7 SVDetect - Program designed to the detection of structural variations
|
|
|
8 from paired-end/mate-pair sequencing data, compatible with SOLiD and Illumina (>=1.3) reads
|
|
|
9
|
|
|
10 Version: 0.8 for Galaxy
|
|
|
11
|
|
|
12 =head1 SYNOPSIS
|
|
|
13
|
|
|
14 SVDetect <command> -conf <configuration_file> [-help] [-man]
|
|
|
15
|
|
|
16 Command:
|
|
|
17
|
|
|
18 linking detection and isolation of links
|
|
|
19 filtering filtering of links according different parameters
|
|
|
20 links2circos links conversion to circos format
|
|
|
21 links2bed paired-ends of links converted to bed format (UCSC)
|
|
|
22 links2SV formatted output to show most significant SVs
|
|
|
23 cnv calculate copy-number profiles
|
|
|
24 ratio2circos ratio conversion to circos density format
|
|
|
25 ratio2bedgraph ratio conversion to bedGraph density format (UCSC)
|
|
|
26
|
|
|
27 =head1 DESCRIPTION
|
|
|
28
|
|
|
29 This is a command-line interface to SVDetect.
|
|
|
30
|
|
|
31
|
|
|
32 =head1 AUTHORS
|
|
|
33
|
|
|
34 Bruno Zeitouni E<lt>bruno.zeitouni@curie.frE<gt>,
|
|
|
35 Valentina Boeva E<lt>valentina.boeva@curie.frE<gt>
|
|
|
36
|
|
|
37 =cut
|
|
|
38
|
|
|
39 # -------------------------------------------------------------------
|
|
|
40
|
|
|
41 use strict;
|
|
|
42 use warnings;
|
|
|
43
|
|
|
44 use Pod::Usage;
|
|
|
45 use Getopt::Long;
|
|
|
46
|
|
|
47 use Config::General;
|
|
|
48 use Tie::IxHash;
|
|
|
49 use FileHandle;
|
|
|
50 use Parallel::ForkManager;
|
|
|
51
|
|
|
52 #::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::#
|
|
|
53 #PARSE THE COMMAND LINE
|
|
|
54 my %OPT;
|
|
|
55 GetOptions(\%OPT,
|
|
|
56 'conf=s',
|
|
|
57 'out1=s', #GALAXY
|
|
|
58 'out2=s', #GALAXY
|
|
|
59 'out3=s', #GALAXY
|
|
|
60 'out4=s', #GALAXY
|
|
|
61 'out5=s', #GALAXY
|
|
|
62 'l=s', #GALAXY
|
|
|
63 'N=s',#GALAXY
|
|
|
64 'help',#GALAXY
|
|
|
65 'man'
|
|
|
66 );
|
|
|
67
|
|
|
68 pod2usage() if $OPT{help};
|
|
|
69 pod2usage(-verbose=>2) if $OPT{man};
|
|
|
70 pod2usage(-message=> "$!", -exitval => 2) if (!defined $OPT{conf});
|
|
|
71
|
|
|
72 pod2usage() if(@ARGV<1);
|
|
|
73
|
|
|
74 tie (my %func, 'Tie::IxHash',linking=>\&createlinks,
|
|
|
75 filtering=>\&filterlinks,
|
|
|
76 links2circos=>\&links2circos,
|
|
|
77 links2bed=>\&links2bed,
|
|
|
78 links2compare=>\&links2compare,
|
|
|
79 links2SV=>\&links2SV,
|
|
|
80 cnv=>\&cnv,
|
|
|
81 ratio2circos=>\&ratio2circos,
|
|
|
82 ratio2bedgraph=>\&ratio2bedgraph);
|
|
|
83
|
|
|
84 foreach my $command (@ARGV){
|
|
|
85 pod2usage(-message=> "Unknown command \"$command\"", -exitval => 2) if (!defined($func{$command}));
|
|
|
86 }
|
|
|
87 #::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::#
|
|
|
88
|
|
|
89
|
|
|
90 #::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::#
|
|
|
91 #READ THE CONFIGURATION FILE
|
|
|
92 my $conf=Config::General->new( -ConfigFile => $OPT{conf},
|
|
|
93 -Tie => "Tie::IxHash",
|
|
|
94 -AllowMultiOptions => 1,
|
|
|
95 -LowerCaseNames => 1,
|
|
|
96 -AutoTrue => 1);
|
|
|
97 my %CONF= $conf->getall;
|
|
|
98 validateconfiguration(\%CONF); #validation of the configuration parameters
|
|
|
99
|
|
|
100 my $SAMTOOLS_BIN_DIR="/bioinfo/local/samtools"; #GALAXY
|
|
|
101
|
|
|
102 my $pt_log_file=$OPT{l}; #GALAXY
|
|
|
103 my $pt_links_file=$OPT{out1} if($OPT{out1}); #GALAXY
|
|
|
104 my $pt_flinks_file=$OPT{out2} if($OPT{out2}); #GALAXY
|
|
|
105 my $pt_sv_file=$OPT{out3} if($OPT{out3}); #GALAXY
|
|
|
106 my $pt_circos_file=$OPT{out4} if($OPT{out4}); #GALAXY
|
|
|
107 my $pt_bed_file=$OPT{out5} if($OPT{out5}); #GALAXY
|
|
|
108
|
|
|
109 $CONF{general}{mates_file}=readlink($CONF{general}{mates_file});#GALAXY
|
|
|
110 $CONF{general}{cmap_file}=readlink($CONF{general}{cmap_file});#GALAXY
|
|
|
111
|
|
|
112 my $log_file=$CONF{general}{output_dir}.$OPT{N}.".svdetect_run.log"; #GALAXY
|
|
|
113 open LOG,">$log_file" or die "$0: can't open ".$log_file.":$!\n";#GALAXY
|
|
|
114 #::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::#
|
|
|
115
|
|
|
116 #::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::#
|
|
|
117 #COMMAND EXECUTION
|
|
|
118 foreach my $command (@ARGV){
|
|
|
119 &{$func{$command}}();
|
|
|
120 }
|
|
|
121 print LOG "-- end\n";#GALAXY
|
|
|
122
|
|
|
123 close LOG;#GALAXY
|
|
|
124 system "rm $pt_log_file ; ln -s $log_file $pt_log_file"; #GALAXY
|
|
|
125 exit(0);
|
|
|
126 #::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::#
|
|
|
127
|
|
|
128
|
|
|
129 #::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::#
|
|
|
130 #FUNCTIONS
|
|
|
131 #------------------------------------------------------------------------------#
|
|
|
132 #MAIN FUNCTION number 1: Detection of links from mate-pairs data
|
|
|
133 sub createlinks{
|
|
|
134
|
|
|
135 my %CHR; #main hash table 1: fragments, links
|
|
|
136 my %CHRID;
|
|
|
137 my @MATEFILES;
|
|
|
138
|
|
|
139 my $output_prefix=$CONF{general}{mates_file}.".".$CONF{general}{sv_type};
|
|
|
140 my @path=split(/\//,$output_prefix);
|
|
|
141 $output_prefix=$CONF{general}{output_dir}.$path[$#path];
|
|
|
142 my $tmp_mates_prefix=$CONF{general}{tmp_dir}."mates/".$path[$#path];
|
|
|
143 my $tmp_links_prefix=$CONF{general}{tmp_dir}."links/".$path[$#path];
|
|
|
144
|
|
|
145 shearingChromosome(\%CHR, \%CHRID, #making the genomic fragment library with the detection parameters
|
|
|
146 $CONF{detection}{window_size},
|
|
|
147 $CONF{detection}{step_length},
|
|
|
148 $CONF{general}{cmap_file});
|
|
|
149
|
|
|
150 if($CONF{detection}{split_mate_file}){
|
|
|
151
|
|
|
152 splitMateFile(\%CHR, \%CHRID, \@MATEFILES, $tmp_mates_prefix,
|
|
|
153 $CONF{general}{sv_type},
|
|
|
154 $CONF{general}{mates_file},
|
|
|
155 $CONF{general}{input_format},
|
|
|
156 $CONF{general}{read_lengths}
|
|
|
157 );
|
|
|
158 }else{
|
|
|
159
|
|
|
160 @MATEFILES=qx{ls $tmp_mates_prefix*} or die "# Error: No splitted mate files already created at $CONF{general}{tmp_dir} :$!";
|
|
|
161 chomp(@MATEFILES);
|
|
|
162 print LOG "# Splitted mate files already created.\n";
|
|
|
163 }
|
|
|
164
|
|
|
165
|
|
|
166 #Parallelization of the linking per chromosome for intra + interchrs
|
|
|
167 my $pm = new Parallel::ForkManager($CONF{general}{num_threads});
|
|
|
168
|
|
|
169 foreach my $matefile (@MATEFILES){
|
|
|
170
|
|
|
171 my $pid = $pm->start and next;
|
|
|
172 getlinks(\%CHR, \%CHRID, $matefile);
|
|
|
173 $pm->finish;
|
|
|
174
|
|
|
175 }
|
|
|
176 $pm->wait_all_children;
|
|
|
177
|
|
|
178 #Merge the chromosome links file into only one
|
|
|
179 my @LINKFILES= qx{ls $tmp_links_prefix*links} or die "# Error: No links files created at $CONF{general}{tmp_dir} :$!";
|
|
|
180 chomp(@LINKFILES);
|
|
|
181 catFiles( \@LINKFILES => "$output_prefix.links" );
|
|
|
182
|
|
|
183 system "rm $pt_links_file; ln -s $output_prefix.links $pt_links_file" if (defined $pt_links_file); #GALAXY
|
|
|
184 print LOG "# Linking end procedure : output created: $output_prefix.links\n";
|
|
|
185 #unlink(@LINKFILES);
|
|
|
186 #unlink(@MATEFILES);
|
|
|
187
|
|
|
188 undef %CHR;
|
|
|
189 undef %CHRID;
|
|
|
190
|
|
|
191 }
|
|
|
192 #------------------------------------------------------------------------------#
|
|
|
193 #------------------------------------------------------------------------------#
|
|
|
194 sub getlinks {
|
|
|
195
|
|
|
196 my ($chr,$chrID,$tmp_mates_prefix)=@_;
|
|
|
197
|
|
|
198 my $tmp_links_prefix=$tmp_mates_prefix;
|
|
|
199 $tmp_links_prefix=~s/\/mates\//\/links\//;
|
|
|
200
|
|
|
201 my %PAIR; #main hash table 2: pairs
|
|
|
202
|
|
|
203 linking($chr,$chrID, \%PAIR, #creation of all links from chromosome coordinates of pairs
|
|
|
204 $CONF{general}{read_lengths},
|
|
|
205 $CONF{detection}{window_size},
|
|
|
206 $CONF{detection}{step_length},
|
|
|
207 $tmp_mates_prefix,
|
|
|
208 $CONF{general}{input_format},
|
|
|
209 $CONF{general}{sv_type},
|
|
|
210 "$tmp_links_prefix.links.mapped"
|
|
|
211 );
|
|
|
212
|
|
|
213 getUniqueLinks("$tmp_links_prefix.links.mapped", #remove the doublons
|
|
|
214 "$tmp_links_prefix.links.unique");
|
|
|
215
|
|
|
216 defineCoordsLinks($chr,$chrID, \%PAIR, #definition of the precise coordinates of links
|
|
|
217 $CONF{general}{input_format},
|
|
|
218 $CONF{general}{sv_type},
|
|
|
219 $CONF{general}{read_lengths},
|
|
|
220 "$tmp_links_prefix.links.unique",
|
|
|
221 "$tmp_links_prefix.links.unique_defined");
|
|
|
222
|
|
|
223 sortLinks("$tmp_links_prefix.links.unique_defined", #sorting links from coordinates
|
|
|
224 "$tmp_links_prefix.links.sorted");
|
|
|
225
|
|
|
226 removeFullyOverlappedLinks("$tmp_links_prefix.links.sorted", #remove redundant links
|
|
|
227 "$tmp_links_prefix.links",1); #file output
|
|
|
228
|
|
|
229
|
|
|
230 undef %PAIR;
|
|
|
231
|
|
|
232 unlink("$tmp_links_prefix.links.mapped",
|
|
|
233 "$tmp_links_prefix.links.unique",
|
|
|
234 "$tmp_links_prefix.links.unique_defined",
|
|
|
235 "$tmp_links_prefix.links.sorted");
|
|
|
236 }
|
|
|
237 #------------------------------------------------------------------------------#
|
|
|
238 #------------------------------------------------------------------------------#
|
|
|
239 sub splitMateFile{
|
|
|
240
|
|
|
241 my ($chr,$chrID,$files_list,$output_prefix,$sv_type,$mates_file,$input_format,$tag_length)=@_;
|
|
|
242
|
|
|
243 print LOG "# Splitting the mate file \"$mates_file\" for parallel processing...\n";
|
|
|
244
|
|
|
245 my %filesHandle;
|
|
|
246
|
|
|
247 #fichier matefile inter
|
|
|
248 if($sv_type=~/^(all|inter)$/){
|
|
|
249 my $newFileName="$output_prefix.interchrs";
|
|
|
250 push(@{$files_list},$newFileName);
|
|
|
251 my $fh = new FileHandle;
|
|
|
252 $fh->open(">$newFileName");
|
|
|
253 $filesHandle{inter}=$fh;
|
|
|
254 }
|
|
|
255
|
|
|
256 #fichiers matefiles intra
|
|
|
257 if($sv_type=~/^(all|intra)$/){
|
|
|
258 foreach my $k (1..$chr->{nb_chrs}){
|
|
|
259 my $newFileName=$output_prefix.".".$chr->{$k}->{name};
|
|
|
260 push(@{$files_list},$newFileName);
|
|
|
261 my $fh = new FileHandle;
|
|
|
262 $fh->open(">$newFileName");
|
|
|
263 $filesHandle{$k}=$fh;
|
|
|
264 }
|
|
|
265 }
|
|
|
266
|
|
|
267 if ($mates_file =~ /.gz$/) {
|
|
|
268 open(MATES, "gunzip -c $mates_file |") or die "$0: can't open ".$mates_file.":$!\n"; #gzcat
|
|
|
269 }elsif($mates_file =~ /.bam$/){
|
|
|
270 open(MATES, "$SAMTOOLS_BIN_DIR/samtools view $mates_file |") or die "$0: can't open ".$mates_file.":$!\n";#GALAXY
|
|
|
271 }else{
|
|
|
272 open MATES, "<".$mates_file or die "$0: can't open ".$mates_file.":$!\n";
|
|
|
273 }
|
|
|
274
|
|
|
275
|
|
|
276 while(<MATES>){
|
|
|
277
|
|
|
278 my @t=split;
|
|
|
279 my ($chr_read1, $chr_read2, $firstbase_read1, $firstbase_read2, $end_order_read1,$end_order_read2);
|
|
|
280
|
|
|
281 next if (!readMateFile(\$chr_read1, \$chr_read2, \$firstbase_read1, \$firstbase_read2, \$end_order_read1, \$end_order_read2, \@t, $input_format,$tag_length));
|
|
|
282
|
|
|
283 next unless (exists $chrID->{$chr_read1} && exists $chrID->{$chr_read2});
|
|
|
284
|
|
|
285 ($chr_read1, $chr_read2)= ($chrID->{$chr_read1},$chrID->{$chr_read2});
|
|
|
286
|
|
|
287 if( ($sv_type=~/^(all|inter)$/) && ($chr_read1 ne $chr_read2) ){
|
|
|
288 my $fh2print=$filesHandle{inter};
|
|
|
289 print $fh2print join("\t",@t)."\n";
|
|
|
290 }
|
|
|
291
|
|
|
292 if( ($sv_type=~/^(all|intra)$/) && ($chr_read1 eq $chr_read2) ){
|
|
|
293 my $fh2print=$filesHandle{$chr_read1};
|
|
|
294 print $fh2print join("\t",@t)."\n";
|
|
|
295
|
|
|
296 }
|
|
|
297 }
|
|
|
298
|
|
|
299 foreach my $name (keys %filesHandle){
|
|
|
300 my $fh=$filesHandle{$name};
|
|
|
301 $fh->close;
|
|
|
302 }
|
|
|
303
|
|
|
304 print LOG "# Splitted mate files of \"$mates_file\" created.\n";
|
|
|
305 }
|
|
|
306
|
|
|
307
|
|
|
308 #------------------------------------------------------------------------------#
|
|
|
309 #------------------------------------------------------------------------------#
|
|
|
310 sub filterlinks{
|
|
|
311
|
|
|
312 my %CHR;
|
|
|
313 my %CHRID;
|
|
|
314 my @LINKFILES;
|
|
|
315 my @FLINKFILES;
|
|
|
316
|
|
|
317 my $output_prefix=$CONF{general}{mates_file}.".".$CONF{general}{sv_type};
|
|
|
318 my @path=split(/\//,$output_prefix);
|
|
|
319 $output_prefix=$CONF{general}{output_dir}.$path[$#path];
|
|
|
320 my $tmp_links_prefix=$CONF{general}{tmp_dir}."links/".$path[$#path];
|
|
|
321
|
|
|
322 createChrHashTables(\%CHR,\%CHRID,
|
|
|
323 $CONF{general}{cmap_file});
|
|
|
324
|
|
|
325 if($CONF{filtering}{split_link_file}){
|
|
|
326
|
|
|
327 splitLinkFile(\%CHR, \%CHRID, \@LINKFILES,
|
|
|
328 $tmp_links_prefix,
|
|
|
329 $CONF{general}{sv_type},
|
|
|
330 "$output_prefix.links",
|
|
|
331 );
|
|
|
332 }else{
|
|
|
333
|
|
|
334 @LINKFILES=qx{ls $tmp_links_prefix*links} or die "# Error: No splitted link files already created\n";
|
|
|
335 chomp(@LINKFILES);
|
|
|
336 print LOG "# Splitted link files already created.\n";
|
|
|
337 }
|
|
|
338
|
|
|
339 #Parallelization of the filtering per chromosome for intra + interchrs
|
|
|
340 my $pm = new Parallel::ForkManager($CONF{general}{num_threads});
|
|
|
341
|
|
|
342 foreach my $linkfile (@LINKFILES){
|
|
|
343
|
|
|
344 my $pid = $pm->start and next;
|
|
|
345 getFilteredlinks(\%CHR, \%CHRID, $linkfile);
|
|
|
346 $pm->finish;
|
|
|
347
|
|
|
348 }
|
|
|
349 $pm->wait_all_children;
|
|
|
350
|
|
|
351 #Merge the chromosome links file into only one
|
|
|
352 @FLINKFILES= qx{ls $tmp_links_prefix*filtered} or die "# Error: No links files created\n";
|
|
|
353 chomp(@FLINKFILES);
|
|
|
354 catFiles( \@FLINKFILES => "$output_prefix.links.filtered" );
|
|
|
355
|
|
|
356 system "rm $pt_flinks_file; ln -s $output_prefix.links.filtered $pt_flinks_file" if (defined $pt_flinks_file); #GALAXY
|
|
|
357 print LOG"# Filtering end procedure : output created: $output_prefix.links.filtered\n";
|
|
|
358
|
|
|
359 undef %CHR;
|
|
|
360 undef %CHRID;
|
|
|
361
|
|
|
362 }
|
|
|
363 #------------------------------------------------------------------------------#
|
|
|
364 #------------------------------------------------------------------------------#
|
|
|
365 sub splitLinkFile{
|
|
|
366
|
|
|
367 my ($chr,$chrID,$files_list,$input_prefix,$sv_type,$link_file)=@_;
|
|
|
368
|
|
|
369 print LOG "# Splitting the link file for parallel processing...\n";
|
|
|
370
|
|
|
371 my %filesHandle;
|
|
|
372
|
|
|
373 #fichier matefile inter
|
|
|
374 if($sv_type=~/^(all|inter)$/){
|
|
|
375 my $newFileName="$input_prefix.interchrs.links";
|
|
|
376 push(@{$files_list},$newFileName);
|
|
|
377 my $fh = new FileHandle;
|
|
|
378 $fh->open(">$newFileName");
|
|
|
379 $filesHandle{inter}=$fh;
|
|
|
380 }
|
|
|
381
|
|
|
382 #fichiers matefiles intra
|
|
|
383 if($sv_type=~/^(all|intra)$/){
|
|
|
384 foreach my $k (1..$chr->{nb_chrs}){
|
|
|
385 my $newFileName=$input_prefix.".".$chr->{$k}->{name}.".links";
|
|
|
386 push(@{$files_list},$newFileName);
|
|
|
387 my $fh = new FileHandle;
|
|
|
388 $fh->open(">$newFileName");
|
|
|
389 $filesHandle{$k}=$fh;
|
|
|
390 }
|
|
|
391 }
|
|
|
392
|
|
|
393 open LINKS, "<".$link_file or die "$0: can't open ".$link_file.":$!\n";
|
|
|
394 while(<LINKS>){
|
|
|
395
|
|
|
396 my @t=split;
|
|
|
397 my ($chr_read1,$chr_read2)=($t[0],$t[3]);
|
|
|
398
|
|
|
399 next unless (exists $chrID->{$chr_read1} && exists $chrID->{$chr_read2});
|
|
|
400
|
|
|
401 ($chr_read1, $chr_read2)= ($chrID->{$chr_read1},$chrID->{$chr_read2});
|
|
|
402
|
|
|
403 if( ($sv_type=~/^(all|inter)$/) && ($chr_read1 ne $chr_read2) ){
|
|
|
404 my $fh2print=$filesHandle{inter};
|
|
|
405 print $fh2print join("\t",@t)."\n";
|
|
|
406 }
|
|
|
407
|
|
|
408 if( ($sv_type=~/^(all|intra)$/) && ($chr_read1 eq $chr_read2) ){
|
|
|
409 my $fh2print=$filesHandle{$chr_read1};
|
|
|
410 print $fh2print join("\t",@t)."\n";
|
|
|
411
|
|
|
412 }
|
|
|
413 }
|
|
|
414
|
|
|
415 foreach my $name (keys %filesHandle){
|
|
|
416 my $fh=$filesHandle{$name};
|
|
|
417 $fh->close;
|
|
|
418 }
|
|
|
419
|
|
|
420 print LOG "# Splitted link files created.\n";
|
|
|
421 }
|
|
|
422
|
|
|
423
|
|
|
424 #------------------------------------------------------------------------------#
|
|
|
425 #------------------------------------------------------------------------------#
|
|
|
426 #MAIN FUNCTION number 2: Filtering processing
|
|
|
427 sub getFilteredlinks {
|
|
|
428
|
|
|
429 my ($chr,$chrID,$tmp_links_prefix)=@_;
|
|
|
430 my %PAIR;
|
|
|
431
|
|
|
432 strandFiltering($chr,$chrID,
|
|
|
433 $CONF{filtering}{nb_pairs_threshold}, #filtering of links
|
|
|
434 $CONF{filtering}{strand_filtering},
|
|
|
435 $CONF{filtering}{chromosomes},
|
|
|
436 $CONF{general}{input_format},
|
|
|
437 $CONF{general}{cmap_file},
|
|
|
438 $CONF{general}{mates_orientation},
|
|
|
439 $CONF{general}{read_lengths},
|
|
|
440 $tmp_links_prefix,
|
|
|
441 "$tmp_links_prefix.filtered",
|
|
|
442 );
|
|
|
443
|
|
|
444 if($CONF{filtering}{strand_filtering}){ #re-definition of links coordinates with strand filtering
|
|
|
445
|
|
|
446 my @tmpfiles;
|
|
|
447
|
|
|
448 rename("$tmp_links_prefix.filtered","$tmp_links_prefix.filtered_unique");
|
|
|
449
|
|
|
450 getUniqueLinks("$tmp_links_prefix.filtered_unique",
|
|
|
451 "$tmp_links_prefix.filtered");
|
|
|
452
|
|
|
453 push(@tmpfiles,"$tmp_links_prefix.filtered_unique");
|
|
|
454
|
|
|
455 if($CONF{filtering}{order_filtering}){ #filtering using the order
|
|
|
456
|
|
|
457 rename("$tmp_links_prefix.filtered","$tmp_links_prefix.filtered_ordered");
|
|
|
458
|
|
|
459 orderFiltering($chr,$chrID,
|
|
|
460 $CONF{filtering}{nb_pairs_threshold},
|
|
|
461 $CONF{filtering}{nb_pairs_order_threshold},
|
|
|
462 $CONF{filtering}{mu_length},
|
|
|
463 $CONF{filtering}{sigma_length},
|
|
|
464 $CONF{general}{mates_orientation},
|
|
|
465 $CONF{general}{read_lengths},
|
|
|
466 "$tmp_links_prefix.filtered_ordered",
|
|
|
467 "$tmp_links_prefix.filtered",
|
|
|
468 );
|
|
|
469
|
|
|
470 push(@tmpfiles,"$tmp_links_prefix.filtered_ordered");
|
|
|
471 }
|
|
|
472
|
|
|
473 if (($CONF{filtering}{insert_size_filtering})&&
|
|
|
474 ($CONF{general}{sv_type} ne 'inter')){
|
|
|
475
|
|
|
476 rename("$tmp_links_prefix.filtered","$tmp_links_prefix.filtered_withoutIndelSize");
|
|
|
477
|
|
|
478 addInsertionInfo($chr,$chrID,
|
|
|
479 $CONF{filtering}{nb_pairs_threshold},
|
|
|
480 $CONF{filtering}{order_filtering},
|
|
|
481 $CONF{filtering}{indel_sigma_threshold},
|
|
|
482 $CONF{filtering}{dup_sigma_threshold},
|
|
|
483 $CONF{filtering}{singleton_sigma_threshold},
|
|
|
484 $CONF{filtering}{mu_length},
|
|
|
485 $CONF{filtering}{sigma_length},
|
|
|
486 $CONF{general}{mates_orientation},
|
|
|
487 $CONF{general}{read_lengths},
|
|
|
488 "$tmp_links_prefix.filtered_withoutIndelSize",
|
|
|
489 "$tmp_links_prefix.filtered"
|
|
|
490 );
|
|
|
491
|
|
|
492 push(@tmpfiles,"$tmp_links_prefix.filtered_withoutIndelSize");
|
|
|
493 }
|
|
|
494
|
|
|
495 sortLinks("$tmp_links_prefix.filtered",
|
|
|
496 "$tmp_links_prefix.filtered_sorted");
|
|
|
497
|
|
|
498 removeFullyOverlappedLinks("$tmp_links_prefix.filtered_sorted",
|
|
|
499 "$tmp_links_prefix.filtered_nodup",
|
|
|
500 );
|
|
|
501
|
|
|
502 postFiltering("$tmp_links_prefix.filtered_nodup",
|
|
|
503 "$tmp_links_prefix.filtered",
|
|
|
504 $CONF{filtering}{final_score_threshold});
|
|
|
505
|
|
|
506 push(@tmpfiles,"$tmp_links_prefix.filtered_sorted","$tmp_links_prefix.filtered_nodup");
|
|
|
507
|
|
|
508 unlink(@tmpfiles);
|
|
|
509
|
|
|
510
|
|
|
511 }
|
|
|
512 undef %PAIR;
|
|
|
513
|
|
|
514 }
|
|
|
515 #------------------------------------------------------------------------------#
|
|
|
516 #------------------------------------------------------------------------------#
|
|
|
517 #MAIN FUNCTION number 3: Circos format conversion for links
|
|
|
518 sub links2circos{
|
|
|
519
|
|
|
520 my $input_file=$CONF{general}{mates_file}.".".$CONF{general}{sv_type}.".links.filtered";
|
|
|
521 my @path=split(/\//,$input_file);
|
|
|
522 $input_file=$CONF{general}{output_dir}.$path[$#path];
|
|
|
523
|
|
|
524 my $output_file.=$input_file.".segdup.txt";
|
|
|
525
|
|
|
526 links2segdup($CONF{circos}{organism_id},
|
|
|
527 $CONF{circos}{colorcode},
|
|
|
528 $input_file,
|
|
|
529 $output_file); #circos file output
|
|
|
530
|
|
|
531 system "rm $pt_circos_file; ln -s $output_file $pt_circos_file" if (defined $pt_circos_file); #GALAXY
|
|
|
532 }
|
|
|
533 #------------------------------------------------------------------------------#
|
|
|
534 #------------------------------------------------------------------------------#
|
|
|
535 #MAIN FUNCTION number 4: Bed format conversion for links
|
|
|
536 sub links2bed{
|
|
|
537
|
|
|
538 my $input_file=$CONF{general}{mates_file}.".".$CONF{general}{sv_type}.".links.filtered";
|
|
|
539 my @path=split(/\//,$input_file);
|
|
|
540 $input_file=$CONF{general}{output_dir}.$path[$#path];
|
|
|
541
|
|
|
542 my $output_file.=$input_file.".bed.txt";
|
|
|
543
|
|
|
544 links2bedfile($CONF{general}{read_lengths},
|
|
|
545 $CONF{bed}{colorcode},
|
|
|
546 $input_file,
|
|
|
547 $output_file); #bed file output
|
|
|
548
|
|
|
549 system "rm $pt_bed_file; ln -s $output_file $pt_bed_file" if (defined $pt_bed_file); #GALAXY
|
|
|
550
|
|
|
551 }
|
|
|
552 #------------------------------------------------------------------------------#
|
|
|
553 #------------------------------------------------------------------------------#
|
|
|
554 #MAIN FUNCTION number 6: Bed format conversion for links
|
|
|
555 sub links2SV{
|
|
|
556
|
|
|
557 my $input_file=$CONF{general}{mates_file}.".".$CONF{general}{sv_type}.".links.filtered";
|
|
|
558
|
|
|
559 my @path=split(/\//,$input_file);
|
|
|
560 $input_file=$CONF{general}{output_dir}.$path[$#path];
|
|
|
561
|
|
|
562 my $output_file.=$input_file.".sv.txt";
|
|
|
563
|
|
|
564
|
|
|
565 links2SVfile( $input_file,
|
|
|
566 $output_file);
|
|
|
567
|
|
|
568 system "rm $pt_sv_file; ln -s $output_file $pt_sv_file" if (defined $pt_sv_file); #GALAXY
|
|
|
569 }
|
|
|
570 #------------------------------------------------------------------------------#
|
|
|
571 #------------------------------------------------------------------------------#
|
|
|
572 #MAIN FUNCTION number 7: copy number variations, coverage ratio calculation
|
|
|
573 sub cnv{
|
|
|
574
|
|
|
575 my %CHR;
|
|
|
576 my %CHRID;
|
|
|
577 my @MATEFILES;
|
|
|
578 my @MATEFILES_REF;
|
|
|
579
|
|
|
580 my $output_prefix=$CONF{general}{mates_file};
|
|
|
581 my $output_prefix_ref=$CONF{detection}{mates_file_ref};
|
|
|
582 my @path=split(/\//,$output_prefix);
|
|
|
583 my @path_ref=split(/\//,$output_prefix_ref);
|
|
|
584 $output_prefix=$CONF{general}{output_dir}.$path[$#path];
|
|
|
585 $output_prefix_ref=$CONF{general}{output_dir}.$path_ref[$#path_ref];
|
|
|
586 my $tmp_mates_prefix=$CONF{general}{tmp_dir}."mates/".$path[$#path];
|
|
|
587 my $tmp_mates_prefix_ref=$CONF{general}{tmp_dir}."mates/".$path_ref[$#path_ref];
|
|
|
588 my $tmp_density_prefix=$CONF{general}{tmp_dir}."density/".$path[$#path];
|
|
|
589
|
|
|
590 shearingChromosome(\%CHR, \%CHRID, #making the genomic fragment library with the detection parameters
|
|
|
591 $CONF{detection}{window_size},
|
|
|
592 $CONF{detection}{step_length},
|
|
|
593 $CONF{general}{cmap_file});
|
|
|
594
|
|
|
595 if($CONF{detection}{split_mate_file}){
|
|
|
596
|
|
|
597 splitMateFile(\%CHR, \%CHRID, \@MATEFILES, $tmp_mates_prefix,
|
|
|
598 "intra",
|
|
|
599 $CONF{general}{mates_file},
|
|
|
600 $CONF{general}{input_format},
|
|
|
601 $CONF{general}{read_lengths}
|
|
|
602 );
|
|
|
603
|
|
|
604 splitMateFile(\%CHR, \%CHRID, \@MATEFILES_REF, $tmp_mates_prefix_ref,
|
|
|
605 "intra",
|
|
|
606 $CONF{detection}{mates_file_ref},
|
|
|
607 $CONF{general}{input_format},
|
|
|
608 $CONF{general}{read_lengths}
|
|
|
609 );
|
|
|
610
|
|
|
611
|
|
|
612 }else{
|
|
|
613
|
|
|
614 @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} :$!";
|
|
|
615 chomp(@MATEFILES);
|
|
|
616 @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} :$!";
|
|
|
617 chomp(@MATEFILES_REF);
|
|
|
618 print LOG "# Splitted sample and reference mate files already created.\n";
|
|
|
619 }
|
|
|
620
|
|
|
621 #Parallelization of the cnv per chromosome
|
|
|
622 my $pm = new Parallel::ForkManager($CONF{general}{num_threads});
|
|
|
623
|
|
|
624 foreach my $file (0..$#MATEFILES){
|
|
|
625
|
|
|
626 my $pid = $pm->start and next;
|
|
|
627
|
|
|
628 densityCalculation(\%CHR, \%CHRID, $file,
|
|
|
629 $CONF{general}{read_lengths},
|
|
|
630 $CONF{detection}{window_size},
|
|
|
631 $CONF{detection}{step_length},
|
|
|
632 \@MATEFILES,
|
|
|
633 \@MATEFILES_REF,
|
|
|
634 $MATEFILES[$file].".density",
|
|
|
635 $CONF{general}{input_format});
|
|
|
636
|
|
|
637 $pm->finish;
|
|
|
638
|
|
|
639 }
|
|
|
640 $pm->wait_all_children;
|
|
|
641
|
|
|
642 #Merge the chromosome links file into only one
|
|
|
643 my @DENSITYFILES= qx{ls $tmp_density_prefix*density} or die "# Error: No density files created at $CONF{general}{tmp_dir} :$!";
|
|
|
644 chomp(@DENSITYFILES);
|
|
|
645 catFiles( \@DENSITYFILES => "$output_prefix.density" );
|
|
|
646
|
|
|
647 print LOG "# cnv end procedure : output created: $output_prefix.density\n";
|
|
|
648
|
|
|
649
|
|
|
650 undef %CHR;
|
|
|
651 undef %CHRID;
|
|
|
652
|
|
|
653 }
|
|
|
654 #------------------------------------------------------------------------------#
|
|
|
655 #------------------------------------------------------------------------------#
|
|
|
656 #MAIN FUNCTION number 8: Circos format conversion for cnv ratios
|
|
|
657 sub ratio2circos{
|
|
|
658
|
|
|
659 my $input_file=$CONF{general}{mates_file}.".density";
|
|
|
660
|
|
|
661 my @path=split(/\//,$input_file);
|
|
|
662 $input_file=$CONF{general}{output_dir}.$path[$#path];
|
|
|
663
|
|
|
664 my $output_file.=$input_file.".segdup.txt";
|
|
|
665
|
|
|
666 ratio2segdup($CONF{circos}{organism_id},
|
|
|
667 $input_file,
|
|
|
668 $output_file);
|
|
|
669 }
|
|
|
670 #------------------------------------------------------------------------------#
|
|
|
671 #MAIN FUNCTION number 9: BedGraph format conversion for cnv ratios
|
|
|
672 sub ratio2bedgraph{
|
|
|
673
|
|
|
674 my $input_file=$CONF{general}{mates_file}.".density";
|
|
|
675
|
|
|
676 my @path=split(/\//,$input_file);
|
|
|
677 $input_file=$CONF{general}{output_dir}.$path[$#path];
|
|
|
678
|
|
|
679 my $output_file.=$input_file.".bedgraph.txt";
|
|
|
680
|
|
|
681 ratio2bedfile($input_file,
|
|
|
682 $output_file); #bed file output
|
|
|
683 }
|
|
|
684 #------------------------------------------------------------------------------#
|
|
|
685 #------------------------------------------------------------------------------#
|
|
|
686 #Creation of the fragment library
|
|
|
687 sub shearingChromosome{
|
|
|
688
|
|
|
689 print LOG "# Making the fragments library...\n";
|
|
|
690
|
|
|
691 my ($chr,$chrID,$window,$step,$cmap_file)=@_; #window and step sizes parameters
|
|
|
692
|
|
|
693 createChrHashTables($chr,$chrID,$cmap_file); #hash tables: chromosome ID <=> chromsomes Name
|
|
|
694
|
|
|
695 foreach my $k (1..$chr->{nb_chrs}){
|
|
|
696
|
|
|
697 print LOG"-- $chr->{$k}->{name}\n";
|
|
|
698
|
|
|
699 my $frag=1;
|
|
|
700 for (my $start=0; $start<$chr->{$k}->{length}; $start+=$step){
|
|
|
701
|
|
|
702 my $end=($start<($chr->{$k}->{length})-$window)? $start+$window-1:($chr->{$k}->{length})-1;
|
|
|
703 $chr->{$k}->{$frag}=[$start,$end]; #creation of fragments, coordinates storage
|
|
|
704
|
|
|
705 if($end==($chr->{$k}->{length})-1){
|
|
|
706 $chr->{$k}->{nb_frag}=$frag; #nb of fragments per chromosome
|
|
|
707 last;
|
|
|
708 }
|
|
|
709 $frag++;
|
|
|
710 }
|
|
|
711 }
|
|
|
712 }
|
|
|
713 #------------------------------------------------------------------------------#
|
|
|
714 #------------------------------------------------------------------------------#
|
|
|
715 #Creation of chromosome hash tables from the cmap file
|
|
|
716 sub createChrHashTables{
|
|
|
717
|
|
|
718 my ($chr,$chrID,$cmap_file)=@_;
|
|
|
719 $chr->{nb_chrs}=0;
|
|
|
720
|
|
|
721 open CMAP, "<".$cmap_file or die "$0: can't open ".$cmap_file.":$!\n";
|
|
|
722 while(<CMAP>){
|
|
|
723
|
|
|
724 if(/^\s+$/){ next;}
|
|
|
725 my ($k,$name,$length) = split;
|
|
|
726 $chr->{$k}->{name}=$name;
|
|
|
727 $chr->{$k}->{length}=$length;
|
|
|
728 $chrID->{$name}=$k;
|
|
|
729 $chr->{nb_chrs}++;
|
|
|
730
|
|
|
731 }
|
|
|
732 close CMAP;
|
|
|
733 }
|
|
|
734 #------------------------------------------------------------------------------#
|
|
|
735 #------------------------------------------------------------------------------#
|
|
|
736 #Read the mate file according the input format file (solid, eland or sam)
|
|
|
737 sub readMateFile{
|
|
|
738
|
|
|
739 my ($chr1,$chr2,$pos1,$pos2,$order1,$order2,$t,$file_type,$tag_length)=@_;
|
|
|
740 my ($strand1,$strand2);
|
|
|
741
|
|
|
742 if($file_type eq "solid"){
|
|
|
743
|
|
|
744 ($$chr1,$$chr2,$$pos1,$$pos2,$$order1,$$order2)=($$t[6],$$t[7],$$t[8]+1,$$t[9]+1,1,2); #0-based
|
|
|
745
|
|
|
746 }else{
|
|
|
747 my ($tag_length1,$tag_length2);
|
|
|
748 ($$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
|
|
|
749 if($file_type eq "eland");
|
|
|
750
|
|
|
751 if($file_type eq "sam"){
|
|
|
752
|
|
|
753 return 0 if ($$t[0]=~/^@/); #header sam filtered out
|
|
|
754
|
|
|
755 ($$chr1,$$chr2,$$pos1,$$pos2)=($$t[2],$$t[6],$$t[3],$$t[7]);
|
|
|
756
|
|
|
757 return 0 if ($$chr1 eq "*" || $$chr2 eq "*");
|
|
|
758
|
|
|
759 $$chr2=$$chr1 if($$chr2 eq "=");
|
|
|
760
|
|
|
761 $strand1 = (($$t[1]&0x0010))? 'R':'F';
|
|
|
762 $strand2 = (($$t[1]&0x0020))? 'R':'F';
|
|
|
763
|
|
|
764 $$order1= (($$t[1]&0x0040))? '1':'2';
|
|
|
765 $$order2= (($$t[1]&0x0080))? '1':'2';
|
|
|
766 $tag_length1 = $tag_length->{$$order1};
|
|
|
767 $tag_length2 = $tag_length->{$$order2};
|
|
|
768 }
|
|
|
769
|
|
|
770 $$pos1 = -($$pos1+$tag_length1) if ($strand1 eq "R"); #get sequencing starts
|
|
|
771 $$pos2 = -($$pos2+$tag_length2) if ($strand2 eq "R");
|
|
|
772 }
|
|
|
773 return 1;
|
|
|
774 }
|
|
|
775 #------------------------------------------------------------------------------#
|
|
|
776 #------------------------------------------------------------------------------#
|
|
|
777 #Parsing of the mates files and creation of links between 2 chromosomal fragments
|
|
|
778 sub linking{
|
|
|
779
|
|
|
780 my ($chr,$chrID,$pair,$tag_length,$window_dist,$step,$mates_file,$input_format,$sv_type,$links_file)=@_;
|
|
|
781 my %link;
|
|
|
782
|
|
|
783 my $record=0;
|
|
|
784 my $nb_links=0;
|
|
|
785 my $warn=10000;
|
|
|
786
|
|
|
787 my @sfile=split(/\./,$mates_file);
|
|
|
788 my $fchr=$sfile[$#sfile];
|
|
|
789
|
|
|
790 my $fh = new FileHandle;
|
|
|
791
|
|
|
792 print LOG "# $fchr : Linking procedure...\n";
|
|
|
793 print LOG "-- file=$mates_file\n".
|
|
|
794 "-- chromosome=$fchr\n".
|
|
|
795 "-- input format=$input_format\n".
|
|
|
796 "-- type=$sv_type\n".
|
|
|
797 "-- read1 length=$tag_length->{1}, read2 length=$tag_length->{2}\n".
|
|
|
798 "-- window size=$window_dist, step length=$step\n";
|
|
|
799
|
|
|
800 if ($mates_file =~ /.gz$/) {
|
|
|
801 $fh->open("gunzip -c $mates_file |") or die "$0: can't open ".$mates_file.":$!\n"; #gzcat
|
|
|
802 }elsif($mates_file =~ /.bam$/){
|
|
|
803 $fh->open("$SAMTOOLS_BIN_DIR/samtools view $mates_file |") or die "$0: can't open ".$mates_file.":$!\n";#GALAXY
|
|
|
804 }else{
|
|
|
805 $fh->open("<".$mates_file) or die "$0: can't open ".$mates_file.":$!\n";
|
|
|
806 }
|
|
|
807
|
|
|
808
|
|
|
809 while(<$fh>){
|
|
|
810
|
|
|
811 my @t=split; #for each mate-pair
|
|
|
812 my $mate=$t[0];
|
|
|
813 my ($chr_read1, $chr_read2, $firstbase_read1, $firstbase_read2, $end_order_read1,$end_order_read2);
|
|
|
814
|
|
|
815 next if(exists $$pair{$mate});
|
|
|
816
|
|
|
817 next if (!readMateFile(\$chr_read1, \$chr_read2, \$firstbase_read1, \$firstbase_read2, \$end_order_read1, \$end_order_read2, \@t, $input_format,$tag_length));
|
|
|
818
|
|
|
819 next unless (exists $chrID->{$chr_read1} && exists $chrID->{$chr_read2});
|
|
|
820
|
|
|
821 ($chr_read1, $chr_read2)= ($chrID->{$chr_read1},$chrID->{$chr_read2});
|
|
|
822
|
|
|
823 if($sv_type ne "all"){
|
|
|
824 if( ($sv_type eq "inter") && ($chr_read1 ne $chr_read2) ||
|
|
|
825 ($sv_type eq "intra") && ($chr_read1 eq $chr_read2) ){
|
|
|
826 }else{
|
|
|
827 next;
|
|
|
828 }
|
|
|
829 }
|
|
|
830
|
|
|
831 $$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)
|
|
|
832
|
|
|
833 $record++;
|
|
|
834
|
|
|
835 my ($coord_start_read1,$coord_end_read1,$coord_start_read2,$coord_end_read2); #get the coordinates of each read
|
|
|
836
|
|
|
837 recupCoords($firstbase_read1,\$coord_start_read1,\$coord_end_read1,$tag_length->{$end_order_read1},$input_format);
|
|
|
838 recupCoords($firstbase_read2,\$coord_start_read2,\$coord_end_read2,$tag_length->{$end_order_read2},$input_format);
|
|
|
839
|
|
|
840 for(my $i=1;$i<=$chr->{$chr_read1}->{'nb_frag'};$i++){ #fast genome parsing for link creation
|
|
|
841
|
|
|
842 if (abs ($coord_start_read1-${$chr->{$chr_read1}->{$i}}[0]) <= $window_dist){
|
|
|
843
|
|
|
844 if(overlap($coord_start_read1,$coord_end_read1,${$chr->{$chr_read1}->{$i}}[0],${$chr->{$chr_read1}->{$i}}[1])){
|
|
|
845
|
|
|
846 for(my $j=1;$j<=$chr->{$chr_read2}->{'nb_frag'};$j++){
|
|
|
847
|
|
|
848 if (abs ($coord_start_read2-${$chr->{$chr_read2}->{$j}}[0]) <= $window_dist) {
|
|
|
849
|
|
|
850 if(overlap($coord_start_read2,$coord_end_read2,${$chr->{$chr_read2}->{$j}}[0],${$chr->{$chr_read2}->{$j}}[1])){
|
|
|
851
|
|
|
852 makeLink(\%link,$chr_read1,$i,$chr_read2,$j,$mate,\$nb_links); #make the link
|
|
|
853 }
|
|
|
854
|
|
|
855 }else{
|
|
|
856
|
|
|
857 $j=getNextFrag($coord_start_read2,$j,${$chr->{$chr_read2}->{$j}}[0],$chr->{$chr_read2}->{nb_frag},$window_dist,$step);
|
|
|
858 }
|
|
|
859 }
|
|
|
860 }
|
|
|
861
|
|
|
862 }else{
|
|
|
863
|
|
|
864 $i=getNextFrag($coord_start_read1,$i,${$chr->{$chr_read1}->{$i}}[0],$chr->{$chr_read1}->{nb_frag},$window_dist,$step);
|
|
|
865 }
|
|
|
866 }
|
|
|
867
|
|
|
868 if($record>=$warn){
|
|
|
869 print LOG "-- $fchr : $warn mate-pairs analysed - $nb_links links done\n";
|
|
|
870 $warn+=10000;
|
|
|
871 }
|
|
|
872 }
|
|
|
873 $fh->close;
|
|
|
874
|
|
|
875 if(!$nb_links){
|
|
|
876 print LOG "-- $fchr : No mate-pairs !\n".
|
|
|
877 "-- $fchr : No links have been found with the selected type of structural variations \($sv_type\)\n";
|
|
|
878 }
|
|
|
879
|
|
|
880 print LOG "-- $fchr : Total : $record mate-pairs analysed - $nb_links links done\n";
|
|
|
881
|
|
|
882 print LOG "-- $fchr : writing...\n";
|
|
|
883
|
|
|
884 $fh = new FileHandle;
|
|
|
885
|
|
|
886 $fh->open(">".$links_file) or die "$0: can't write in the output ".$links_file." :$!\n";
|
|
|
887
|
|
|
888 foreach my $chr1 ( sort { $a <=> $b} keys %link){ #Sorted links output
|
|
|
889
|
|
|
890 foreach my $chr2 ( sort { $a <=> $b} keys %{$link{$chr1}}){
|
|
|
891
|
|
|
892 foreach my $frag1 ( sort { $a <=> $b} keys %{$link{$chr1}{$chr2}}){
|
|
|
893
|
|
|
894 foreach my $frag2 ( sort { $a <=> $b} keys %{$link{$chr1}{$chr2}{$frag1}}){
|
|
|
895
|
|
|
896 my @count=split(",",$link{$chr1}{$chr2}{$frag1}{$frag2});
|
|
|
897 print $fh "$chr->{$chr1}->{name}\t".(${$chr->{$chr1}->{$frag1}}[0]+1)."\t".(${$chr->{$chr1}->{$frag1}}[1]+1)."\t".
|
|
|
898 "$chr->{$chr2}->{name}\t".(${$chr->{$chr2}->{$frag2}}[0]+1)."\t".(${$chr->{$chr2}->{$frag2}}[1]+1)."\t".
|
|
|
899 scalar @count."\t". #nb of read
|
|
|
900 $link{$chr1}{$chr2}{$frag1}{$frag2}."\n"; #mate list
|
|
|
901 }
|
|
|
902 }
|
|
|
903 }
|
|
|
904 }
|
|
|
905
|
|
|
906 $fh->close;
|
|
|
907
|
|
|
908 undef %link;
|
|
|
909
|
|
|
910 }
|
|
|
911 #------------------------------------------------------------------------------#
|
|
|
912 #------------------------------------------------------------------------------#
|
|
|
913 #remove exact links doublons according to the mate list
|
|
|
914 sub getUniqueLinks{
|
|
|
915
|
|
|
916 my ($links_file,$nrlinks_file)=@_;
|
|
|
917 my %links;
|
|
|
918 my %pt;
|
|
|
919 my $nb_links;
|
|
|
920 my $n=1;
|
|
|
921
|
|
|
922 my $record=0;
|
|
|
923 my $warn=300000;
|
|
|
924
|
|
|
925 my @sfile=split(/\./,$links_file);
|
|
|
926 my $fchr=$sfile[$#sfile-2];
|
|
|
927
|
|
|
928 my $fh = new FileHandle;
|
|
|
929
|
|
|
930 print LOG "# $fchr : Getting unique links...\n";
|
|
|
931 $fh->open("<$links_file") or die "$0: can't open $links_file :$!\n";
|
|
|
932
|
|
|
933 while(<$fh>){
|
|
|
934
|
|
|
935 my @t=split;
|
|
|
936 my $mates=$t[7];
|
|
|
937 $record++;
|
|
|
938
|
|
|
939 if(!exists $links{$mates}){ #Unique links selection
|
|
|
940
|
|
|
941 $links{$mates}=[@t];
|
|
|
942 $pt{$n}=$links{$mates};
|
|
|
943 $n++;
|
|
|
944
|
|
|
945
|
|
|
946 }else{ #get the link coordinates from the mate-pairs list
|
|
|
947
|
|
|
948 for my $i (1,2,4,5){ #get the shortest regions
|
|
|
949
|
|
|
950 $links{$mates}->[$i]=($t[$i]>$links{$mates}->[$i])? $t[$i]:$links{$mates}->[$i] #maximum start
|
|
|
951 if($i==1 || $i==4);
|
|
|
952 $links{$mates}->[$i]=($t[$i]<$links{$mates}->[$i])? $t[$i]:$links{$mates}->[$i] #minimum end
|
|
|
953 if($i==2 || $i==5);
|
|
|
954 }
|
|
|
955 }
|
|
|
956 if($record>=$warn){
|
|
|
957 print LOG "-- $fchr : $warn links analysed - ".($n-1)." unique links done\n";
|
|
|
958 $warn+=300000;
|
|
|
959 }
|
|
|
960 }
|
|
|
961 $fh->close;
|
|
|
962
|
|
|
963 $nb_links=$n-1;
|
|
|
964 print LOG "-- $fchr : Total : $record links analysed - $nb_links unique links done\n";
|
|
|
965
|
|
|
966 $fh = new FileHandle;
|
|
|
967 $fh->open(">$nrlinks_file") or die "$0: can't write in the output: $nrlinks_file :$!\n";
|
|
|
968 print LOG "-- $fchr : writing...\n";
|
|
|
969 for my $i (1..$nb_links){
|
|
|
970
|
|
|
971 print $fh join("\t",@{$pt{$i}})."\n"; #all links output
|
|
|
972 }
|
|
|
973
|
|
|
974 $fh->close;
|
|
|
975
|
|
|
976 undef %links;
|
|
|
977 undef %pt;
|
|
|
978
|
|
|
979 }
|
|
|
980 #------------------------------------------------------------------------------#
|
|
|
981 #------------------------------------------------------------------------------#
|
|
|
982 #get the new coordinates of each link from the mate list
|
|
|
983 sub defineCoordsLinks{
|
|
|
984
|
|
|
985 my ($chr,$chrID,$pair,$input_format,$sv_type,$tag_length,$links_file,$clinks_file)=@_;
|
|
|
986
|
|
|
987 my @sfile=split(/\./,$links_file);
|
|
|
988 my $fchr=$sfile[$#sfile-2];
|
|
|
989
|
|
|
990 my $fh = new FileHandle;
|
|
|
991 my $fh2 = new FileHandle;
|
|
|
992
|
|
|
993 $fh->open("<$links_file") or die "$0: can't open $links_file :$!\n";
|
|
|
994 $fh2->open(">$clinks_file") or die "$0: can't write in the output: $clinks_file :$!\n";
|
|
|
995
|
|
|
996 print LOG "# $fchr : Defining precise link coordinates...\n";
|
|
|
997
|
|
|
998 my $record=0;
|
|
|
999 my $warn=100000;
|
|
|
1000
|
|
|
1001 my %coords;
|
|
|
1002 my %strands;
|
|
|
1003 my %order;
|
|
|
1004 my %ends_order;
|
|
|
1005
|
|
|
1006 while(<$fh>){
|
|
|
1007
|
|
|
1008
|
|
|
1009 my ($col1,$col2)=(1,2); #for an intrachromosomal link
|
|
|
1010 my $diffchr=0; #difference between chr1 and chr2
|
|
|
1011 my ($chr1,$chr2,$mates_list,$npairs)=(split)[0,3,7,8];
|
|
|
1012 ($chr1,$chr2) = ($chrID->{$chr1},$chrID->{$chr2});
|
|
|
1013 if ($chr1 != $chr2){ #for an interchromosomal link
|
|
|
1014 $col1=$col2=0; #no distinction
|
|
|
1015 $diffchr=1;
|
|
|
1016 }
|
|
|
1017
|
|
|
1018 my @pairs=split(",",$mates_list);
|
|
|
1019
|
|
|
1020 $coords{$col1}{$chr1}->{start}=undef;
|
|
|
1021 $coords{$col1}{$chr1}->{end}=undef;
|
|
|
1022 $coords{$col2}{$chr2}->{start}=undef;
|
|
|
1023 $coords{$col2}{$chr2}->{end}=undef;
|
|
|
1024 $strands{$col1}{$chr1}=undef;
|
|
|
1025 $strands{$col2}{$chr2}=undef;
|
|
|
1026 $ends_order{$col1}{$chr1}=undef;
|
|
|
1027 $ends_order{$col2}{$chr2}=undef;
|
|
|
1028
|
|
|
1029
|
|
|
1030 $order{$col1}{$chr1}->{index}->{1}=undef;
|
|
|
1031 $order{$col1}{$chr1}->{index}->{2}=undef;
|
|
|
1032 $order{$col2}{$chr2}->{index}->{1}=undef;
|
|
|
1033 $order{$col2}{$chr2}->{index}->{2}=undef;
|
|
|
1034 $order{$col1}{$chr1}->{order}=undef;
|
|
|
1035 $order{$col2}{$chr2}->{order}=undef;
|
|
|
1036
|
|
|
1037 $record++;
|
|
|
1038
|
|
|
1039 for my $p (0..$#pairs){ #for each pair
|
|
|
1040
|
|
|
1041 my ($coord_start_read1,$coord_end_read1);
|
|
|
1042 my ($coord_start_read2,$coord_end_read2);
|
|
|
1043 my $strand_read1=recupCoords(${$$pair{$pairs[$p]}}[2],\$coord_start_read1,\$coord_end_read1,$tag_length->{${$$pair{$pairs[$p]}}[4]},$input_format);
|
|
|
1044 my $strand_read2=recupCoords(${$$pair{$pairs[$p]}}[3],\$coord_start_read2,\$coord_end_read2,$tag_length->{${$$pair{$pairs[$p]}}[5]},$input_format);
|
|
|
1045
|
|
|
1046 if(!$diffchr){ #for a intrachromosomal link
|
|
|
1047 if($coord_start_read2<$coord_start_read1){ #get the closer start coordinate for each column
|
|
|
1048 ($col1,$col2)=(2,1);
|
|
|
1049 }else{
|
|
|
1050 ($col1,$col2)=(1,2);
|
|
|
1051 }
|
|
|
1052 }
|
|
|
1053
|
|
|
1054 push(@{$coords{$col1}{${$$pair{$pairs[$p]}}[0]}->{start}},$coord_start_read1); #get coords and strands of f3 and r3 reads
|
|
|
1055 push(@{$coords{$col1}{${$$pair{$pairs[$p]}}[0]}->{end}},$coord_end_read1);
|
|
|
1056 push(@{$coords{$col2}{${$$pair{$pairs[$p]}}[1]}->{start}},$coord_start_read2);
|
|
|
1057 push(@{$coords{$col2}{${$$pair{$pairs[$p]}}[1]}->{end}},$coord_end_read2);
|
|
|
1058 push(@{$strands{$col1}{${$$pair{$pairs[$p]}}[0]}},$strand_read1);
|
|
|
1059 push(@{$strands{$col2}{${$$pair{$pairs[$p]}}[1]}},$strand_read2);
|
|
|
1060 push(@{$ends_order{$col1}{${$$pair{$pairs[$p]}}[0]}},${$$pair{$pairs[$p]}}[4]);
|
|
|
1061 push(@{$ends_order{$col2}{${$$pair{$pairs[$p]}}[1]}},${$$pair{$pairs[$p]}}[5]);
|
|
|
1062 }
|
|
|
1063
|
|
|
1064 ($col1,$col2)=(1,2) if(!$diffchr);
|
|
|
1065
|
|
|
1066 my $coord_start_chr1=min(min(@{$coords{$col1}{$chr1}->{start}}),min(@{$coords{$col1}{$chr1}->{end}})); #get the biggest region
|
|
|
1067 my $coord_end_chr1=max(max(@{$coords{$col1}{$chr1}->{start}}),max(@{$coords{$col1}{$chr1}->{end}}));
|
|
|
1068 my $coord_start_chr2=min(min(@{$coords{$col2}{$chr2}->{start}}),min(@{$coords{$col2}{$chr2}->{end}}));
|
|
|
1069 my $coord_end_chr2=max(max(@{$coords{$col2}{$chr2}->{start}}),max(@{$coords{$col2}{$chr2}->{end}}));
|
|
|
1070
|
|
|
1071 @{$order{$col1}{$chr1}->{index}->{1}}= sort {${$coords{$col1}{$chr1}->{start}}[$a] <=> ${$coords{$col1}{$chr1}->{start}}[$b]} 0 .. $#{$coords{$col1}{$chr1}->{start}};
|
|
|
1072 @{$order{$col2}{$chr2}->{index}->{1}}= sort {${$coords{$col2}{$chr2}->{start}}[$a] <=> ${$coords{$col2}{$chr2}->{start}}[$b]} 0 .. $#{$coords{$col2}{$chr2}->{start}};
|
|
|
1073
|
|
|
1074 foreach my $i (@{$order{$col1}{$chr1}->{index}->{1}}){ #get the rank of the chr2 reads according to the sorted chr1 reads (start coordinate sorting)
|
|
|
1075 foreach my $j (@{$order{$col2}{$chr2}->{index}->{1}}){
|
|
|
1076
|
|
|
1077 if(${$order{$col1}{$chr1}->{index}->{1}}[$i] == ${$order{$col2}{$chr2}->{index}->{1}}[$j]){
|
|
|
1078 ${$order{$col1}{$chr1}->{index}->{2}}[$i]=$i;
|
|
|
1079 ${$order{$col2}{$chr2}->{index}->{2}}[$i]=$j;
|
|
|
1080 last;
|
|
|
1081 }
|
|
|
1082 }
|
|
|
1083 }
|
|
|
1084
|
|
|
1085 foreach my $i (@{$order{$col1}{$chr1}->{index}->{2}}){ #use rank chr1 as an ID
|
|
|
1086 foreach my $j (@{$order{$col2}{$chr2}->{index}->{2}}){
|
|
|
1087
|
|
|
1088 if(${$order{$col1}{$chr1}->{index}->{2}}[$i] == ${$order{$col2}{$chr2}->{index}->{2}}[$j]){
|
|
|
1089 ${$order{$col1}{$chr1}->{order}}[$i]=$i+1;
|
|
|
1090 ${$order{$col2}{$chr2}->{order}}[$i]=$j+1;
|
|
|
1091 last;
|
|
|
1092 }
|
|
|
1093 }
|
|
|
1094 }
|
|
|
1095
|
|
|
1096 @pairs=sortTablebyIndex(\@{$order{$col1}{$chr1}->{index}->{1}},\@pairs);#sorting of the pairs, strands, and start coords from the sorted chr2 reads
|
|
|
1097 @{$strands{$col1}{$chr1}}=sortTablebyIndex(\@{$order{$col1}{$chr1}->{index}->{1}},$strands{$col1}{$chr1});
|
|
|
1098 @{$strands{$col2}{$chr2}}=sortTablebyIndex(\@{$order{$col1}{$chr1}->{index}->{1}},$strands{$col2}{$chr2});
|
|
|
1099 @{$ends_order{$col1}{$chr1}}=sortTablebyIndex(\@{$order{$col1}{$chr1}->{index}->{1}},$ends_order{$col1}{$chr1});
|
|
|
1100 @{$ends_order{$col2}{$chr2}}=sortTablebyIndex(\@{$order{$col1}{$chr1}->{index}->{1}},$ends_order{$col2}{$chr2});
|
|
|
1101 @{$coords{$col1}{$chr1}->{start}}=sortTablebyIndex(\@{$order{$col1}{$chr1}->{index}->{1}},$coords{$col1}{$chr1}->{start});
|
|
|
1102 @{$coords{$col2}{$chr2}->{start}}=sortTablebyIndex(\@{$order{$col1}{$chr1}->{index}->{1}},$coords{$col2}{$chr2}->{start});
|
|
|
1103
|
|
|
1104
|
|
|
1105 my @link=($chr->{$chr1}->{name}, $coord_start_chr1 , $coord_end_chr1, #all information output
|
|
|
1106 $chr->{$chr2}->{name}, $coord_start_chr2 , $coord_end_chr2,
|
|
|
1107 scalar @pairs,
|
|
|
1108 join(",",@pairs),
|
|
|
1109 join(",",@{$strands{$col1}{$chr1}}),
|
|
|
1110 join(",",@{$strands{$col2}{$chr2}}),
|
|
|
1111 join(",",@{$ends_order{$col1}{$chr1}}),
|
|
|
1112 join(",",@{$ends_order{$col2}{$chr2}}),
|
|
|
1113 join(",",@{$order{$col1}{$chr1}->{order}}),
|
|
|
1114 join(",",@{$order{$col2}{$chr2}->{order}}),
|
|
|
1115 join(",",@{$coords{$col1}{$chr1}->{start}}),
|
|
|
1116 join(",",@{$coords{$col2}{$chr2}->{start}}));
|
|
|
1117
|
|
|
1118 print $fh2 join("\t",@link)."\n";
|
|
|
1119
|
|
|
1120 if($record>=$warn){
|
|
|
1121 print LOG "-- $fchr : $warn links processed\n";
|
|
|
1122 $warn+=100000;
|
|
|
1123 }
|
|
|
1124 }
|
|
|
1125 $fh->close;
|
|
|
1126 $fh2->close;
|
|
|
1127
|
|
|
1128 print LOG "-- $fchr : Total : $record links processed\n";
|
|
|
1129
|
|
|
1130 }
|
|
|
1131 #------------------------------------------------------------------------------#
|
|
|
1132 #------------------------------------------------------------------------------#
|
|
|
1133 #Sort links according the concerned chromosomes and their coordinates
|
|
|
1134 sub sortLinks{
|
|
|
1135
|
|
|
1136 my ($links_file,$sortedlinks_file,$unique)=@_;
|
|
|
1137
|
|
|
1138 my @sfile=split(/\./,$links_file);
|
|
|
1139 my $fchr=$sfile[$#sfile-2];
|
|
|
1140
|
|
|
1141
|
|
|
1142 print LOG "# $fchr : Sorting links...\n";
|
|
|
1143
|
|
|
1144 my $pipe=($unique)? "| sort -u":"";
|
|
|
1145 system "sort -k 1,1 -k 4,4 -k 2,2n -k 5,5n -k 8,8n $links_file $pipe > $sortedlinks_file";
|
|
|
1146
|
|
|
1147 }
|
|
|
1148 #------------------------------------------------------------------------------#
|
|
|
1149 #------------------------------------------------------------------------------#
|
|
|
1150 #removal of fully overlapped links
|
|
|
1151 sub removeFullyOverlappedLinks{
|
|
|
1152
|
|
|
1153 my ($links_file,$nrlinks_file,$warn_out)=@_;
|
|
|
1154
|
|
|
1155 my %pt;
|
|
|
1156 my $n=1;
|
|
|
1157
|
|
|
1158 my @sfile=split(/\./,$links_file);
|
|
|
1159 my $fchr=$sfile[$#sfile-2];
|
|
|
1160
|
|
|
1161 my $fh = new FileHandle;
|
|
|
1162
|
|
|
1163 $fh->open("<$links_file") or die "$0: can't open $links_file :$!\n";
|
|
|
1164 while(<$fh>){
|
|
|
1165
|
|
|
1166 my @t=split("\t",$_);
|
|
|
1167 $pt{$n}=[@t];
|
|
|
1168 $n++;
|
|
|
1169 }
|
|
|
1170 $fh->close;
|
|
|
1171
|
|
|
1172 my $nb_links=$n-1;
|
|
|
1173 my $nb=$nb_links;
|
|
|
1174
|
|
|
1175 my %pt2;
|
|
|
1176 my $nb2=1;
|
|
|
1177 my $record=0;
|
|
|
1178 my $warn=10000;
|
|
|
1179
|
|
|
1180 print LOG "# $fchr : Removing fully overlapped links...\n";
|
|
|
1181
|
|
|
1182 LINK:
|
|
|
1183
|
|
|
1184 for my $i (1..$nb){
|
|
|
1185
|
|
|
1186 my @link=();
|
|
|
1187 my @next_link=();
|
|
|
1188 my $ind1=$i;
|
|
|
1189
|
|
|
1190 $record++;
|
|
|
1191 if($record>=$warn){
|
|
|
1192 print LOG "-- $fchr : $warn unique links analysed - ".($nb2-1)." non-overlapped links done\n";
|
|
|
1193 $warn+=10000;
|
|
|
1194 }
|
|
|
1195
|
|
|
1196 if(exists $pt{$ind1}){
|
|
|
1197 @link=@{$pt{$ind1}}; #link1
|
|
|
1198 }else{
|
|
|
1199 next LINK;
|
|
|
1200 }
|
|
|
1201
|
|
|
1202 my ($chr1,$start1,$end1,$chr2,$start2,$end2)=($link[0],$link[1],$link[2],$link[3],$link[4],$link[5]); #get info of link1
|
|
|
1203 my @mates=deleteBadOrderSensePairs(split(",",$link[7]));
|
|
|
1204
|
|
|
1205 my $ind2=$ind1+1;
|
|
|
1206 $ind2++ while (!exists $pt{$ind2}&& $ind2<=$nb); #get the next found link
|
|
|
1207
|
|
|
1208 if($ind2<=$nb){
|
|
|
1209
|
|
|
1210 @next_link=@{$pt{$ind2}}; #link2
|
|
|
1211 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
|
|
|
1212 my @next_mates=deleteBadOrderSensePairs(split(",",$next_link[7]));
|
|
|
1213
|
|
|
1214 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
|
|
|
1215
|
|
|
1216 if(!overlap($start2,$end2,$start4,$end4)){ #if no overlap with chr2 coordinates ->next link2
|
|
|
1217
|
|
|
1218 $ind2++;
|
|
|
1219 $ind2++ while (!exists $pt{$ind2}&& $ind2<=$nb);
|
|
|
1220
|
|
|
1221 if($ind2>$nb){ #if no more link in the file -> save link1
|
|
|
1222
|
|
|
1223 $pt2{$nb2}=\@link;
|
|
|
1224 $nb2++;
|
|
|
1225 next LINK;
|
|
|
1226 }
|
|
|
1227
|
|
|
1228 @next_link=@{$pt{$ind2}};
|
|
|
1229 ($chr3,$start3,$end3,$chr4,$start4,$end4)=($next_link[0],$next_link[1],$next_link[2],$next_link[3],$next_link[4],$next_link[5]);
|
|
|
1230 @next_mates=deleteBadOrderSensePairs(split(",",$next_link[7]));
|
|
|
1231 next;
|
|
|
1232 }
|
|
|
1233
|
|
|
1234 my %mates=map{$_ =>1} @mates; #get the equal number of mates
|
|
|
1235 my @same_mates = grep( $mates{$_}, @next_mates );
|
|
|
1236 my $nb_mates= scalar @same_mates;
|
|
|
1237
|
|
|
1238 if($nb_mates == scalar @mates){
|
|
|
1239
|
|
|
1240 delete $pt{$ind1}; #if pairs of link 1 are all included in link 2 -> delete link1
|
|
|
1241 next LINK; #go to link2, link2 becomes link1
|
|
|
1242
|
|
|
1243 }else{
|
|
|
1244 delete $pt{$ind2} if($nb_mates == scalar @next_mates); #if pairs of link2 are all included in link 1 -> delete link2
|
|
|
1245 $ind2++; #we continue by checking the next link2
|
|
|
1246 $ind2++ while (!exists $pt{$ind2}&& $ind2<=$nb);
|
|
|
1247
|
|
|
1248 if($ind2>$nb){ #if no more link in the file -> save link1
|
|
|
1249
|
|
|
1250 $pt2{$nb2}=\@link;
|
|
|
1251 $nb2++;
|
|
|
1252 next LINK;
|
|
|
1253 }
|
|
|
1254
|
|
|
1255 @next_link=@{$pt{$ind2}}; #get info of link2
|
|
|
1256 ($chr3,$start3,$end3,$chr4,$start4,$end4)=($next_link[0],$next_link[1],$next_link[2],$next_link[3],$next_link[4],$next_link[5]);
|
|
|
1257 @next_mates=deleteBadOrderSensePairs(split(",",$next_link[7]));
|
|
|
1258
|
|
|
1259 }
|
|
|
1260 }
|
|
|
1261 }
|
|
|
1262 $pt2{$nb2}=\@link; #if no (more) link with chr1 coordinates overlap -> save link1
|
|
|
1263 $nb2++;
|
|
|
1264 }
|
|
|
1265
|
|
|
1266 print LOG "-- $fchr : Total : $nb_links unique links analysed - ".($nb2-1)." non-overlapped links done\n";
|
|
|
1267
|
|
|
1268 #OUTPUT
|
|
|
1269
|
|
|
1270 $fh = new FileHandle;
|
|
|
1271 $fh->open(">$nrlinks_file") or die "$0: can't write in the output: $nrlinks_file :$!\n";
|
|
|
1272 print LOG "-- $fchr : writing...\n";
|
|
|
1273 for my $i (1..$nb2-1){
|
|
|
1274
|
|
|
1275 print $fh join("\t",@{$pt2{$i}}); #all links output
|
|
|
1276 }
|
|
|
1277
|
|
|
1278 close $fh;
|
|
|
1279
|
|
|
1280 print LOG "-- $fchr : output created: $nrlinks_file\n" if($warn_out);
|
|
|
1281
|
|
|
1282 undef %pt;
|
|
|
1283 undef %pt2;
|
|
|
1284 }
|
|
|
1285 #------------------------------------------------------------------------------#
|
|
|
1286 sub postFiltering {
|
|
|
1287
|
|
|
1288 my ($links_file,$pflinks_file, $finalScore_thres)=@_;
|
|
|
1289
|
|
|
1290 my @sfile=split(/\./,$links_file);
|
|
|
1291 my $fchr=$sfile[$#sfile-2];
|
|
|
1292
|
|
|
1293
|
|
|
1294 my ($nb,$nb2)=(0,0);
|
|
|
1295
|
|
|
1296 print LOG "# $fchr : Post-filtering links...\n";
|
|
|
1297 print LOG "-- $fchr : final score threshold = $finalScore_thres\n";
|
|
|
1298
|
|
|
1299 my $fh = new FileHandle;
|
|
|
1300 my $fh2 = new FileHandle;
|
|
|
1301
|
|
|
1302 $fh->open("<$links_file") or die "$0: can't open $links_file :$!\n";
|
|
|
1303 $fh2->open(">$pflinks_file") or die "$0: can't write in the output: $pflinks_file :$!\n";
|
|
|
1304
|
|
|
1305
|
|
|
1306 while(<$fh>){
|
|
|
1307
|
|
|
1308 my @t=split("\t",$_);
|
|
|
1309 my $score=$t[$#t-1];
|
|
|
1310
|
|
|
1311 if($score >= $finalScore_thres){
|
|
|
1312 print $fh2 join("\t", @t);
|
|
|
1313 $nb2++;
|
|
|
1314 }
|
|
|
1315 $nb++;
|
|
|
1316 }
|
|
|
1317 $fh->close;
|
|
|
1318 $fh2->close;
|
|
|
1319
|
|
|
1320 print LOG "-- $fchr : Total : $nb unique links analysed - $nb2 links kept\n";
|
|
|
1321 print LOG "-- $fchr : output created: $pflinks_file\n";
|
|
|
1322 }
|
|
|
1323
|
|
|
1324
|
|
|
1325
|
|
|
1326 #------------------------------------------------------------------------------#
|
|
|
1327 #------------------------------------------------------------------------------#
|
|
|
1328 #Filtering of the links
|
|
|
1329 sub strandFiltering{
|
|
|
1330
|
|
|
1331 my($chr,$chrID,$pairs_threshold,$strand_filtering,$chromosomes,
|
|
|
1332 $input_format,$cmap_file,$mate_sense, $tag_length,$links_file,$flinks_file)=@_;
|
|
|
1333
|
|
|
1334 my @sfile=split(/\./,$links_file);
|
|
|
1335 my $fchr=$sfile[$#sfile-1];
|
|
|
1336
|
|
|
1337
|
|
|
1338 my %chrs;
|
|
|
1339 my %chrs1;
|
|
|
1340 my %chrs2;
|
|
|
1341 my $nb_chrs;
|
|
|
1342 my $exclude;
|
|
|
1343
|
|
|
1344 if($chromosomes){
|
|
|
1345 my @chrs=split(",",$chromosomes);
|
|
|
1346 $nb_chrs=scalar @chrs;
|
|
|
1347 $exclude=($chrs[0]=~/^\-/)? 1:0;
|
|
|
1348 for my $chrName (@chrs){
|
|
|
1349 $chrName=~s/^(\-)//;
|
|
|
1350 my $col=($chrName=~s/_(1|2)$//);
|
|
|
1351
|
|
|
1352 if(!$col){
|
|
|
1353 $chrs{$chrID->{$chrName}}=undef
|
|
|
1354 }else{
|
|
|
1355 $chrs1{$chrID->{$chrName}}=undef if($1==1);
|
|
|
1356 $chrs2{$chrID->{$chrName}}=undef if($1==2);
|
|
|
1357 }
|
|
|
1358 }
|
|
|
1359 }
|
|
|
1360
|
|
|
1361 my $record=0;
|
|
|
1362 my $nb_links=0;
|
|
|
1363 my $warn=10000;
|
|
|
1364
|
|
|
1365 my $sens_ratio_threshold=0.6;
|
|
|
1366
|
|
|
1367 print LOG "\# Filtering procedure...\n";
|
|
|
1368 print LOG "\# Number of pairs and strand filtering...\n";
|
|
|
1369 print LOG "-- file=$links_file\n";
|
|
|
1370 print LOG "-- nb_pairs_threshold=$pairs_threshold, strand_filtering=".(($strand_filtering)? "yes":"no").
|
|
|
1371 ", chromosomes=".(($chromosomes)? "$chromosomes":"all")."\n";
|
|
|
1372
|
|
|
1373
|
|
|
1374
|
|
|
1375 my $fh = new FileHandle;
|
|
|
1376 my $fh2 = new FileHandle;
|
|
|
1377
|
|
|
1378 $fh->open("<$links_file") or die "$0: can't open $links_file :$!\n";
|
|
|
1379 $fh2->open(">$flinks_file") or die "$0: can't write in the output: $flinks_file :$!\n";
|
|
|
1380
|
|
|
1381 while(<$fh>){
|
|
|
1382
|
|
|
1383 my @t=split; #for each link
|
|
|
1384 my $is_good=1;
|
|
|
1385 $record++;
|
|
|
1386
|
|
|
1387
|
|
|
1388 if($chromosomes){
|
|
|
1389
|
|
|
1390 my ($chr1,$chr2)=($chrID->{$t[0]},$chrID->{$t[3]});
|
|
|
1391
|
|
|
1392 if(!$exclude){
|
|
|
1393 $is_good=(exists $chrs{$chr1} && exists $chrs{$chr2})? 1:0;
|
|
|
1394 $is_good=(exists $chrs1{$chr1} && exists $chrs2{$chr2})? 1:0 if(!$is_good);
|
|
|
1395 $is_good=($nb_chrs==1 && (exists $chrs1{$chr1} || exists $chrs2{$chr2}))? 1:0 if(!$is_good);
|
|
|
1396 }else{
|
|
|
1397 $is_good=(exists $chrs{$chr1} || exists $chrs{$chr2})? 0:1;
|
|
|
1398 $is_good=(exists $chrs1{$chr1} || exists $chrs2{$chr2})? 0:1 if($is_good);
|
|
|
1399 }
|
|
|
1400 }
|
|
|
1401
|
|
|
1402 $is_good = ($is_good && $t[6] >= $pairs_threshold)? 1 :0; #filtering according the number of pairs
|
|
|
1403 if($is_good && $strand_filtering){ #if filtering according the strand sense
|
|
|
1404
|
|
|
1405 my @mates=split(/,/,$t[7]); #get the concordant pairs in the strand sense
|
|
|
1406 my @strands1=split(/,/,$t[8]);
|
|
|
1407 my @strands2=split(/,/,$t[9]);
|
|
|
1408
|
|
|
1409 my %mate_class=( 'FF' => 0, 'RR' => 0, 'FR' => 0, 'RF' => 0);
|
|
|
1410
|
|
|
1411 my %mate_reverse=( 'FF' => 'RR', 'RR' => 'FF', #group1: FF,RR
|
|
|
1412 'FR' => 'RF', 'RF' => 'FR'); #group2: FR,RF
|
|
|
1413
|
|
|
1414 my %mate_class2=( $mate_sense=>"NORMAL_SENSE", inverseSense($mate_sense)=>"NORMAL_SENSE",
|
|
|
1415 substr($mate_sense,0,1).inverseSense(substr($mate_sense,1,1))=>"REVERSE_SENSE",
|
|
|
1416 inverseSense(substr($mate_sense,0,1)).substr($mate_sense,1,1)=>"REVERSE_SENSE");
|
|
|
1417
|
|
|
1418 if($t[6] == 1){
|
|
|
1419
|
|
|
1420 push(@t,$mate_class2{$strands1[0].$strands2[0]},"1/1",1,1);
|
|
|
1421
|
|
|
1422 }else{
|
|
|
1423
|
|
|
1424 tie (my %class,'Tie::IxHash');
|
|
|
1425 my $split;
|
|
|
1426
|
|
|
1427 foreach my $i (0..$#mates){
|
|
|
1428 $mate_class{$strands1[$i].$strands2[$i]}++; #get the over-represented group
|
|
|
1429 }
|
|
|
1430
|
|
|
1431 my $nb_same_sens_class=$mate_class{FF}+$mate_class{RR};
|
|
|
1432 my $nb_diff_sens_class=$mate_class{FR}+$mate_class{RF};
|
|
|
1433 my $sens_ratio=max($nb_same_sens_class,$nb_diff_sens_class)/($nb_same_sens_class+$nb_diff_sens_class);
|
|
|
1434
|
|
|
1435 if($sens_ratio < $sens_ratio_threshold){
|
|
|
1436 %class=(1=>'FF', 2=>'FR');
|
|
|
1437 $split=1;
|
|
|
1438 }else{
|
|
|
1439 $class{1}=($nb_same_sens_class > $nb_diff_sens_class)? 'FF':'FR'; #if yes get the concerned class
|
|
|
1440 $split=0;
|
|
|
1441 }
|
|
|
1442
|
|
|
1443 $is_good=getConsistentSenseLinks(\@t,\@mates,\@strands1,\@strands2,$tag_length,$mate_sense,\%mate_reverse,\%mate_class2,\%class,$split,$pairs_threshold);
|
|
|
1444 }
|
|
|
1445 }
|
|
|
1446
|
|
|
1447 if($is_good){ #PRINT
|
|
|
1448
|
|
|
1449 my $nb=scalar @t;
|
|
|
1450 if($nb > 20){
|
|
|
1451 my @t2=splice(@t,0,20);
|
|
|
1452 print $fh2 join("\t",@t2)."\n";
|
|
|
1453 $nb_links++;
|
|
|
1454 }
|
|
|
1455 $nb_links++;
|
|
|
1456 print $fh2 join("\t",@t)."\n";
|
|
|
1457 }
|
|
|
1458
|
|
|
1459 if($record>=$warn){
|
|
|
1460 print LOG "-- $fchr : $warn links analysed - $nb_links links kept\n";
|
|
|
1461 $warn+=10000;
|
|
|
1462 }
|
|
|
1463 }
|
|
|
1464 $fh->close;
|
|
|
1465 $fh2->close;
|
|
|
1466
|
|
|
1467 print LOG "-- $fchr : No links have been found with the selected filtering parameters\n" if(!$nb_links);
|
|
|
1468
|
|
|
1469 print LOG "-- $fchr : Total : $record links analysed - $nb_links links kept\n";
|
|
|
1470
|
|
|
1471
|
|
|
1472 }
|
|
|
1473 #------------------------------------------------------------------------------#
|
|
|
1474 #------------------------------------------------------------------------------#
|
|
|
1475 sub getConsistentSenseLinks{
|
|
|
1476
|
|
|
1477 my ($t,$mates,$strands1,$strands2,$tag_length,$mate_sense, $mate_reverse,$mate_class2, $class, $split,$thres)=@_;
|
|
|
1478
|
|
|
1479 my $npairs=scalar @$mates;
|
|
|
1480
|
|
|
1481 my @ends_order1 = split (/,/,$$t[10]);
|
|
|
1482 my @ends_order2 = split (/,/,$$t[11]);
|
|
|
1483 my @order1 = split (/,/,$$t[12]);
|
|
|
1484 my @order2 = split (/,/,$$t[13]);
|
|
|
1485 my @positions1 = split (/,/,$$t[14]);
|
|
|
1486 my @positions2 = split (/,/,$$t[15]);
|
|
|
1487
|
|
|
1488 my @newlink;
|
|
|
1489
|
|
|
1490 foreach my $ind (keys %{$class} ){
|
|
|
1491
|
|
|
1492 tie (my %flink,'Tie::IxHash');
|
|
|
1493 my @orders2remove=();
|
|
|
1494
|
|
|
1495 foreach my $i (0..$#{$mates}){ #get the pairs belonging the over-represented group
|
|
|
1496
|
|
|
1497 if((($$strands1[$i].$$strands2[$i]) eq $$class{$ind}) || (($$strands1[$i].$$strands2[$i]) eq $$mate_reverse{$$class{$ind}})){
|
|
|
1498 push(@{$flink{mates}},$$mates[$i]);
|
|
|
1499 push(@{$flink{strands1}},$$strands1[$i]);
|
|
|
1500 push(@{$flink{strands2}},$$strands2[$i]);
|
|
|
1501 push(@{$flink{ends_order1}},$ends_order1[$i]);
|
|
|
1502 push(@{$flink{ends_order2}},$ends_order2[$i]);
|
|
|
1503 push(@{$flink{positions1}},$positions1[$i]);
|
|
|
1504 push(@{$flink{positions2}},$positions2[$i]);
|
|
|
1505
|
|
|
1506 }else{
|
|
|
1507 push(@orders2remove,$order1[$i]);
|
|
|
1508 }
|
|
|
1509 }
|
|
|
1510
|
|
|
1511 @{$flink{order1}}=();
|
|
|
1512 @{$flink{order2}}=();
|
|
|
1513 if(scalar @orders2remove > 0){
|
|
|
1514 getNewOrders(\@order1,\@order2,\@orders2remove,$flink{order1},$flink{order2})
|
|
|
1515 }else{
|
|
|
1516 @{$flink{order1}}=@order1;
|
|
|
1517 @{$flink{order2}}=@order2;
|
|
|
1518 }
|
|
|
1519
|
|
|
1520 my @ends1; getEnds(\@ends1,$flink{positions1},$flink{strands1},$flink{ends_order1},$tag_length);
|
|
|
1521 my @ends2; getEnds(\@ends2,$flink{positions2},$flink{strands2},$flink{ends_order2},$tag_length);
|
|
|
1522
|
|
|
1523 my $fnpairs=scalar @{$flink{mates}};
|
|
|
1524 my $strand_filtering_ratio=$fnpairs."/".$npairs;
|
|
|
1525 my $real_ratio=$fnpairs/$npairs;
|
|
|
1526
|
|
|
1527 if($fnpairs>=$thres){ #filtering according the number of pairs
|
|
|
1528
|
|
|
1529 push(@newlink,
|
|
|
1530 $$t[0],
|
|
|
1531 min(min(@{$flink{positions1}}),min(@ends1)),
|
|
|
1532 max(max(@{$flink{positions1}}),max(@ends1)),
|
|
|
1533 $$t[3],
|
|
|
1534 min(min(@{$flink{positions2}}),min(@ends2)),
|
|
|
1535 max(max(@{$flink{positions2}}),max(@ends2)),
|
|
|
1536 $fnpairs,
|
|
|
1537 join(",",@{$flink{mates}}),
|
|
|
1538 join(",",@{$flink{strands1}}),
|
|
|
1539 join(",",@{$flink{strands2}}),
|
|
|
1540 join(",",@{$flink{ends_order1}}),
|
|
|
1541 join(",",@{$flink{ends_order2}}),
|
|
|
1542 join(",",@{$flink{order1}}),
|
|
|
1543 join(",",@{$flink{order2}}),
|
|
|
1544 join(",",@{$flink{positions1}}),
|
|
|
1545 join(",",@{$flink{positions2}}),
|
|
|
1546 $$mate_class2{${$flink{strands1}}[0].${$flink{strands2}}[0]},
|
|
|
1547 $strand_filtering_ratio,
|
|
|
1548 $real_ratio,
|
|
|
1549 $npairs
|
|
|
1550 );
|
|
|
1551 }
|
|
|
1552 }
|
|
|
1553
|
|
|
1554 if (grep {defined($_)} @newlink) {
|
|
|
1555 @$t=@newlink;
|
|
|
1556 return 1
|
|
|
1557 }
|
|
|
1558 return 0;
|
|
|
1559
|
|
|
1560 }
|
|
|
1561 #------------------------------------------------------------------------------#
|
|
|
1562 #------------------------------------------------------------------------------#
|
|
|
1563 sub getNewOrders{
|
|
|
1564
|
|
|
1565 my($tab1,$tab2,$list,$newtab1,$newtab2)=@_;
|
|
|
1566 my $j=1;
|
|
|
1567 my $k=1;
|
|
|
1568 for my $i (0..$#{$tab2}){
|
|
|
1569 my $c=0;
|
|
|
1570 for my $j (0..$#{$list}){
|
|
|
1571 $c++ if(${$list}[$j] < ${$tab2}[$i]);
|
|
|
1572 if(${$list}[$j] == ${$tab2}[$i]){
|
|
|
1573 $c=-1; last;
|
|
|
1574 }
|
|
|
1575 }
|
|
|
1576 if($c!=-1){
|
|
|
1577 push(@{$newtab2}, ${$tab2}[$i]-$c);
|
|
|
1578 push(@{$newtab1}, $k);
|
|
|
1579 $k++;
|
|
|
1580 }
|
|
|
1581 }
|
|
|
1582 }
|
|
|
1583
|
|
|
1584 #------------------------------------------------------------------------------#
|
|
|
1585 #------------------------------------------------------------------------------#
|
|
|
1586 #Filtering of the links using their order
|
|
|
1587 sub orderFiltering {
|
|
|
1588
|
|
|
1589 my ($chr,$chrID,$nb_pairs_threshold,$nb_pairs_order_threshold,$mu,$sigma,$mate_sense,$tag_length,$links_file,$flinks_file)=@_;
|
|
|
1590
|
|
|
1591 my @sfile=split(/\./,$links_file);
|
|
|
1592 my $fchr=$sfile[$#sfile-2];
|
|
|
1593
|
|
|
1594
|
|
|
1595 my $diff_sense_ends=(($mate_sense eq "FR") || ($mate_sense eq "RF"))? 1:0;
|
|
|
1596
|
|
|
1597 my $record=0;
|
|
|
1598 my $warn=10000;
|
|
|
1599 my $nb_links=0;
|
|
|
1600
|
|
|
1601 my $quant05 = 1.644854;
|
|
|
1602 my $quant001 = 3.090232;
|
|
|
1603 my $alphaDist = $quant05 * 2 * $sigma;
|
|
|
1604 my $maxFragmentLength = &floor($quant001 * $sigma + $mu);
|
|
|
1605
|
|
|
1606 print LOG "\# Filtering by order...\n";
|
|
|
1607 print LOG "-- mu length=$mu, sigma length=$sigma, nb pairs order threshold=$nb_pairs_order_threshold\n";
|
|
|
1608 print LOG "-- distance between comparable pairs was set to $alphaDist\n";
|
|
|
1609 print LOG "-- maximal fragment length was set to $maxFragmentLength\n";
|
|
|
1610
|
|
|
1611
|
|
|
1612 my $fh = new FileHandle;
|
|
|
1613 my $fh2 = new FileHandle;
|
|
|
1614
|
|
|
1615 $fh->open("<$links_file") or die "$0: can't open $links_file :$!\n";
|
|
|
1616 $fh2->open(">$flinks_file") or die "$0: can't write in the output: $flinks_file :$!\n";
|
|
|
1617
|
|
|
1618 while(<$fh>){
|
|
|
1619
|
|
|
1620 $record++;
|
|
|
1621 my @t = split;
|
|
|
1622 my ($chr1,$chr2,$mates_list)=@t[0,3,7];
|
|
|
1623 my @pairs=split(",",$mates_list);
|
|
|
1624 ($chr1,$chr2) = ($chrID->{$chr1},$chrID->{$chr2});
|
|
|
1625 my ($coord_start_chr1,$coord_end_chr1,$coord_start_chr2,$coord_end_chr2) = @t[1,2,4,5];
|
|
|
1626 my $numberOfPairs = $t[6];
|
|
|
1627 my @strand1 = split (/,/,$t[8]);
|
|
|
1628 my @strand2 = split (/,/,$t[9]);
|
|
|
1629 my @ends_order1 = split (/,/,$t[10]);
|
|
|
1630 my @ends_order2 = split (/,/,$t[11]);
|
|
|
1631 my @order1 = split (/,/,$t[12]);
|
|
|
1632 my @order2 = split (/,/,$t[13]);
|
|
|
1633 my @positions1 = split (/,/,$t[14]);
|
|
|
1634 my @positions2 = split (/,/,$t[15]);
|
|
|
1635 my @ends1; getEnds(\@ends1,\@positions1,\@strand1,\@ends_order1,$tag_length);
|
|
|
1636 my @ends2; getEnds(\@ends2,\@positions2,\@strand2,\@ends_order2,$tag_length);
|
|
|
1637 my $clusterCoordinates_chr1;
|
|
|
1638 my $clusterCoordinates_chr2;
|
|
|
1639 my $reads_left = 0;
|
|
|
1640
|
|
|
1641 my $ifRenv = $t[16];
|
|
|
1642 my $strand_ratio_filtering=$t[17];
|
|
|
1643
|
|
|
1644 #kind of strand filtering. For example, will keep only FFF-RRR from a link FFRF-RRRF if <F-R> orientation is correct
|
|
|
1645 my ($singleBreakpoint, %badInFRSense) = findBadInFRSenseSOLiDSolexa(\@strand1,\@strand2,\@ends_order1,\@ends_order2,\@order1,\@order2,$mate_sense);
|
|
|
1646 #find pairs type F-RRRR or FFFF-R in the case if <R-F> orientation is correct
|
|
|
1647 #These pairs are annotated as BED pairs forever! They won't be recycled!
|
|
|
1648 my $table;
|
|
|
1649 for my $i (0..$numberOfPairs-1) { #fill the table with non adequate pairs: pairID numberOfNonAdPairs nonAdPairIDs
|
|
|
1650 my $nonAdeq = 0;
|
|
|
1651 for my $j (0..$i-1) {
|
|
|
1652 if (exists($table->{$j}->{$i})) {
|
|
|
1653 $nonAdeq++;
|
|
|
1654 $table->{$i}->{$j} = 1;
|
|
|
1655 }
|
|
|
1656 }
|
|
|
1657 for my $j ($i+1..$numberOfPairs-1) {
|
|
|
1658 if ($positions1[$j]-$positions1[$i]>$alphaDist) {
|
|
|
1659 if (&reversed ($i,$j,$ifRenv,\@positions2)) {
|
|
|
1660 $nonAdeq++;
|
|
|
1661 $table->{$i}->{$j} = 1;
|
|
|
1662 }
|
|
|
1663 }
|
|
|
1664 }
|
|
|
1665 $table->{$i}->{nonAdeq} = $nonAdeq;
|
|
|
1666 }
|
|
|
1667
|
|
|
1668 for my $bad (keys %badInFRSense) { #remove pairs type F-RRRR or FFFF-R in the case of <R-F> orientation
|
|
|
1669 &remove($bad,$table);
|
|
|
1670 }
|
|
|
1671
|
|
|
1672 my @falseReads;
|
|
|
1673 #RRRR-F -> RRRR or R-FFFF -> FFFF
|
|
|
1674 @falseReads = findBadInRFSenseSOLiDSolexa(\@strand1,\@ends_order1,$mate_sense, keys %{$table});
|
|
|
1675 #these pairs will be recycled later as $secondTable
|
|
|
1676 for my $bad (@falseReads) {
|
|
|
1677 &remove($bad,$table);
|
|
|
1678 }
|
|
|
1679
|
|
|
1680 my $bad = &check($table);
|
|
|
1681 while ($bad ne "OK") { #clear the table to reject non adequate pairs in the sense of ORDER
|
|
|
1682 # push (@falseReads, $bad); remove completely!!!
|
|
|
1683 &remove($bad,$table);
|
|
|
1684 $bad = &check($table);
|
|
|
1685 }
|
|
|
1686
|
|
|
1687 $reads_left = scalar keys %{$table};
|
|
|
1688 my $coord_start_chr1_cluster1 = min(min(@positions1[sort {$a<=>$b} keys %{$table}]),min(@ends1[sort {$a<=>$b} keys %{$table}]));
|
|
|
1689 my $coord_end_chr1_cluster1 = max(max(@positions1[sort {$a<=>$b} keys %{$table}]),max(@ends1[sort {$a<=>$b} keys %{$table}]));
|
|
|
1690 my $coord_start_chr2_cluster1 = min(min(@positions2[sort {$a<=>$b} keys %{$table}]),min(@ends2[sort {$a<=>$b} keys %{$table}]));
|
|
|
1691 my $coord_end_chr2_cluster1 = max(max(@positions2[sort {$a<=>$b} keys %{$table}]),max(@ends2[sort {$a<=>$b} keys %{$table}]));
|
|
|
1692
|
|
|
1693 $clusterCoordinates_chr1 = '('.$coord_start_chr1_cluster1.','.$coord_end_chr1_cluster1.')';
|
|
|
1694 $clusterCoordinates_chr2 = '('.$coord_start_chr2_cluster1.','.$coord_end_chr2_cluster1.')';
|
|
|
1695
|
|
|
1696 my $ifBalanced = 'UNBAL';
|
|
|
1697 my $secondTable;
|
|
|
1698 my $clusterCoordinates;
|
|
|
1699 my ($break_pont_chr1,$break_pont_chr2);
|
|
|
1700
|
|
|
1701 my $signatureType="";
|
|
|
1702
|
|
|
1703 my $maxCoord1 =$chr->{$chr1}->{length};
|
|
|
1704 my $maxCoord2 =$chr->{$chr2}->{length};
|
|
|
1705
|
|
|
1706 if (scalar @falseReads) {
|
|
|
1707 @falseReads = sort @falseReads;
|
|
|
1708 #now delete FRFR choosing the majority
|
|
|
1709 my @newfalseReads; #find and remove pairs type RRRR-F or R-FFFF
|
|
|
1710 @newfalseReads = findBadInRFSenseSOLiDSolexa(\@strand1,\@ends_order1,$mate_sense,@falseReads); #these @newfalseReads won't be recycled
|
|
|
1711 my %hashTmp;
|
|
|
1712 for my $count1 (0..scalar(@falseReads)-1) {
|
|
|
1713 my $i = $falseReads[$count1];
|
|
|
1714 $hashTmp{$i} = 1;
|
|
|
1715 for my $bad (@newfalseReads) {
|
|
|
1716 if ($bad == $i) {
|
|
|
1717 delete $hashTmp{$i};
|
|
|
1718 next;
|
|
|
1719 }
|
|
|
1720 }
|
|
|
1721 }
|
|
|
1722 @falseReads = sort keys %hashTmp; #what is left
|
|
|
1723 for my $count1 (0..scalar(@falseReads)-1) { #fill the table for reads which were previously rejected
|
|
|
1724 my $nonAdeq = 0;
|
|
|
1725 my $i = $falseReads[$count1];
|
|
|
1726
|
|
|
1727 for my $count2 (0..$count1-1) {
|
|
|
1728 my $j = $falseReads[$count2];
|
|
|
1729 if (exists($secondTable->{$j}->{$i})) {
|
|
|
1730 $nonAdeq++;
|
|
|
1731 $secondTable->{$i}->{$j} = 1;
|
|
|
1732 }
|
|
|
1733 }
|
|
|
1734 for my $count2 ($count1+1..scalar(@falseReads)-1) {
|
|
|
1735 my $j = $falseReads[$count2];
|
|
|
1736 if ($positions1[$j]-$positions1[$i]>$alphaDist) {
|
|
|
1737 if (&reversed ($i,$j,$ifRenv,\@positions2)) {
|
|
|
1738 $nonAdeq++;
|
|
|
1739 $secondTable->{$i}->{$j} = 1;
|
|
|
1740 }
|
|
|
1741 }
|
|
|
1742 }
|
|
|
1743 $secondTable->{$i}->{nonAdeq} = $nonAdeq;
|
|
|
1744 }
|
|
|
1745
|
|
|
1746 my @falseReads2;
|
|
|
1747 my $bad = &check($secondTable);
|
|
|
1748 while ($bad ne "OK") { #clear the table to reject non adequate pairs
|
|
|
1749 push (@falseReads2, $bad);
|
|
|
1750 &remove($bad,$secondTable);
|
|
|
1751 $bad = &check($secondTable);
|
|
|
1752 }
|
|
|
1753 if (scalar keys %{$secondTable} >= $nb_pairs_order_threshold) {
|
|
|
1754 my $coord_start_chr1_cluster2 = min(min(@positions1[sort {$a<=>$b} keys %{$secondTable}]),min(@ends1[sort {$a<=>$b} keys %{$secondTable}]));
|
|
|
1755 my $coord_end_chr1_cluster2 = max(max(@positions1[sort {$a<=>$b} keys %{$secondTable}]),max(@ends1[sort {$a<=>$b} keys %{$secondTable}]));
|
|
|
1756 my $coord_start_chr2_cluster2 = min(min(@positions2[sort {$a<=>$b} keys %{$secondTable}]),min(@ends2[sort {$a<=>$b} keys %{$secondTable}]));
|
|
|
1757 my $coord_end_chr2_cluster2 = max(max(@positions2[sort {$a<=>$b} keys %{$secondTable}]),max(@ends2[sort {$a<=>$b} keys %{$secondTable}]));
|
|
|
1758
|
|
|
1759 $ifBalanced = 'BAL';
|
|
|
1760
|
|
|
1761 if ($ifBalanced eq 'BAL') {
|
|
|
1762
|
|
|
1763 if (scalar keys %{$table} < $nb_pairs_order_threshold) {
|
|
|
1764 $ifBalanced = 'UNBAL'; #kill cluster 1!
|
|
|
1765 ($table,$secondTable)=($secondTable,$table); #this means that one needs to exchange cluster1 with cluster2
|
|
|
1766 $reads_left = scalar keys %{$table};
|
|
|
1767 $coord_start_chr1_cluster1 = $coord_start_chr1_cluster2;
|
|
|
1768 $coord_end_chr1_cluster1 = $coord_end_chr1_cluster2;
|
|
|
1769 $coord_start_chr2_cluster1 = $coord_start_chr2_cluster2;
|
|
|
1770 $coord_end_chr2_cluster1 = $coord_end_chr2_cluster2;
|
|
|
1771 $clusterCoordinates_chr1 = '('.$coord_start_chr1_cluster1.','.$coord_end_chr1_cluster1.')';
|
|
|
1772 $clusterCoordinates_chr2 = '('.$coord_start_chr2_cluster1.','.$coord_end_chr2_cluster1.')';
|
|
|
1773
|
|
|
1774 } else {
|
|
|
1775
|
|
|
1776 $reads_left += scalar keys %{$secondTable};
|
|
|
1777 next if ($reads_left < $nb_pairs_threshold);
|
|
|
1778
|
|
|
1779 if ($coord_end_chr1_cluster2 < $coord_start_chr1_cluster1) {
|
|
|
1780 ($table,$secondTable)=($secondTable,$table); #this means that one needs to exchange cluster1 with cluster2
|
|
|
1781
|
|
|
1782 ($coord_start_chr1_cluster1,$coord_start_chr1_cluster2) = ($coord_start_chr1_cluster2,$coord_start_chr1_cluster1);
|
|
|
1783 ($coord_end_chr1_cluster1,$coord_end_chr1_cluster2)=($coord_end_chr1_cluster2,$coord_end_chr1_cluster1);
|
|
|
1784 ($coord_start_chr2_cluster1,$coord_start_chr2_cluster2)=($coord_start_chr2_cluster2,$coord_start_chr2_cluster1);
|
|
|
1785 ($coord_end_chr2_cluster1 , $coord_end_chr2_cluster2)=($coord_end_chr2_cluster2 , $coord_end_chr2_cluster1);
|
|
|
1786
|
|
|
1787 $clusterCoordinates_chr1 = '('.$coord_start_chr1_cluster1.','.$coord_end_chr1_cluster1.'),'.$clusterCoordinates_chr1;
|
|
|
1788 $clusterCoordinates_chr2 = '('.$coord_start_chr2_cluster1.','.$coord_end_chr2_cluster1.'),'.$clusterCoordinates_chr2;
|
|
|
1789 }else {
|
|
|
1790 $clusterCoordinates_chr1 .= ',('.$coord_start_chr1_cluster2.','.$coord_end_chr1_cluster2.')';
|
|
|
1791 $clusterCoordinates_chr2 .= ',('.$coord_start_chr2_cluster2.','.$coord_end_chr2_cluster2.')';
|
|
|
1792 }
|
|
|
1793 $coord_start_chr1 = min($coord_start_chr1_cluster1,$coord_start_chr1_cluster2);
|
|
|
1794 $coord_end_chr1 = max($coord_end_chr1_cluster1,$coord_end_chr1_cluster2);
|
|
|
1795 $coord_start_chr2 = min($coord_start_chr2_cluster1,$coord_start_chr2_cluster2);
|
|
|
1796 $coord_end_chr2 = max($coord_end_chr2_cluster1,$coord_end_chr2_cluster2);
|
|
|
1797 #to calculate breakpoints one need to take into account read orientation in claster..
|
|
|
1798 my $leftLetterOk = substr($mate_sense, 0, 1); #R
|
|
|
1799 my $rightLetterOk = substr($mate_sense, 1, 1); #F
|
|
|
1800
|
|
|
1801
|
|
|
1802 my @index1 = keys %{$table};
|
|
|
1803 my @index2 = keys %{$secondTable};
|
|
|
1804
|
|
|
1805 my (@generalStrand1,@generalStrand2) = 0;
|
|
|
1806
|
|
|
1807 if ($leftLetterOk eq $rightLetterOk) { #SOLID mate-pairs
|
|
|
1808 $leftLetterOk = 'R';
|
|
|
1809 $rightLetterOk = 'F';
|
|
|
1810 @generalStrand1 = translateSolidToRF(\@strand1,\@ends_order1);
|
|
|
1811 @generalStrand2 = translateSolidToRF(\@strand2,\@ends_order2);
|
|
|
1812 } else {
|
|
|
1813 @generalStrand1 = @strand1;
|
|
|
1814 @generalStrand2 = @strand2; # TODO check if it is correct
|
|
|
1815 }
|
|
|
1816 if ($generalStrand1[$index1[0]] eq $leftLetterOk && $generalStrand1[$index2[0]] eq $rightLetterOk) { #(R,F)
|
|
|
1817 $break_pont_chr1 = '('.$coord_end_chr1_cluster1.','.$coord_start_chr1_cluster2.')';
|
|
|
1818
|
|
|
1819 if ($generalStrand2[$index1[0]] eq $rightLetterOk && $generalStrand2[$index2[0]] eq $leftLetterOk) {
|
|
|
1820 if ($coord_end_chr2_cluster1 >= $coord_end_chr2_cluster2) {
|
|
|
1821 $break_pont_chr2 = '('.$coord_end_chr2_cluster2.','.$coord_start_chr2_cluster1.')';
|
|
|
1822 $signatureType = "TRANSLOC";
|
|
|
1823 } else {
|
|
|
1824 $break_pont_chr2 = '('.max(($coord_end_chr2_cluster1-$maxFragmentLength),1).','.$coord_start_chr2_cluster1.')';
|
|
|
1825 $break_pont_chr2 .= ',('.$coord_end_chr2_cluster2.','.min(($coord_start_chr2_cluster2+$maxFragmentLength),$maxCoord2).')';
|
|
|
1826 $signatureType = "INS_FRAGMT";
|
|
|
1827 }
|
|
|
1828
|
|
|
1829 } elsif ($generalStrand2[$index1[0]] eq $leftLetterOk && $generalStrand2[$index2[0]] eq $rightLetterOk) {
|
|
|
1830 if ($coord_end_chr2_cluster1 >= $coord_end_chr2_cluster2) {
|
|
|
1831 $break_pont_chr2 = '('.max(($coord_end_chr2_cluster2-$maxFragmentLength),1).','.$coord_start_chr2_cluster2.')';
|
|
|
1832 $break_pont_chr2 .= ',('.$coord_end_chr2_cluster1.','.min(($coord_start_chr2_cluster1+$maxFragmentLength),$maxCoord2).')';
|
|
|
1833 $signatureType = "INV_INS_FRAGMT";
|
|
|
1834 } else {
|
|
|
1835 $break_pont_chr2 = '('.$coord_end_chr2_cluster1.','.$coord_start_chr2_cluster2.')';
|
|
|
1836 $signatureType = "INV_TRANSLOC";
|
|
|
1837 }
|
|
|
1838 } else {
|
|
|
1839 #should not occur
|
|
|
1840 print STDERR "\nError in orderFiltering\n\n";
|
|
|
1841 }
|
|
|
1842 }
|
|
|
1843
|
|
|
1844 elsif ($generalStrand1[$index1[0]] eq $rightLetterOk && $generalStrand1[$index2[0]] eq $leftLetterOk) { #(F,R)
|
|
|
1845 $break_pont_chr1 = '('.max(($coord_end_chr1_cluster1-$maxFragmentLength),1).','.$coord_start_chr1_cluster1.')';
|
|
|
1846 $break_pont_chr1 .= ',('.$coord_end_chr1_cluster2.','.min(($coord_start_chr1_cluster2+$maxFragmentLength),$maxCoord1).')';
|
|
|
1847 if ($generalStrand2[$index1[0]] eq $rightLetterOk && $generalStrand2[$index2[0]] eq $leftLetterOk) {
|
|
|
1848 if ($coord_end_chr2_cluster1 >= $coord_end_chr2_cluster2) {
|
|
|
1849 $break_pont_chr2 = '('.$coord_end_chr2_cluster2.','.$coord_start_chr2_cluster1.')';
|
|
|
1850 $signatureType = "INV_INS_FRAGMT";
|
|
|
1851 } else {
|
|
|
1852 $break_pont_chr2 = '('.max(($coord_end_chr2_cluster1-$maxFragmentLength),1).','.$coord_start_chr2_cluster1.')';
|
|
|
1853 $break_pont_chr2 .= ',('.$coord_end_chr2_cluster2.','.min(($coord_start_chr2_cluster2+$maxFragmentLength),$maxCoord2).')';
|
|
|
1854 $signatureType = "INV_COAMPLICON";
|
|
|
1855 }
|
|
|
1856
|
|
|
1857 } elsif ($generalStrand2[$index1[0]] eq $leftLetterOk && $generalStrand2[$index2[0]] eq $rightLetterOk) {
|
|
|
1858 if ($coord_end_chr2_cluster1 >= $coord_end_chr2_cluster2) {
|
|
|
1859 $break_pont_chr2 = '('.max(($coord_end_chr2_cluster2-$maxFragmentLength),1).','.$coord_start_chr2_cluster2.')';
|
|
|
1860 $break_pont_chr2 .= ',('.$coord_end_chr2_cluster1.','.min(($coord_start_chr2_cluster1+$maxFragmentLength),$maxCoord2).')';
|
|
|
1861 $signatureType = "COAMPLICON";
|
|
|
1862 } else {
|
|
|
1863 $break_pont_chr2 = '('.$coord_end_chr2_cluster1.','.$coord_start_chr2_cluster2.')';
|
|
|
1864 $signatureType = "INS_FRAGMT";
|
|
|
1865 }
|
|
|
1866 } else {
|
|
|
1867 #should not occur
|
|
|
1868 $signatureType = "UNDEFINED";
|
|
|
1869 }
|
|
|
1870 }
|
|
|
1871 else { # (F,F) or (R,R) something strange. We will discard the smallest cluster
|
|
|
1872 $ifBalanced = 'UNBAL';
|
|
|
1873 if (scalar keys %{$secondTable} > scalar keys %{$table}) {
|
|
|
1874 ($table,$secondTable)=($secondTable,$table); #this means that one needs to exchange cluster1 with cluster2
|
|
|
1875
|
|
|
1876 $coord_start_chr1_cluster1 = $coord_start_chr1_cluster2;
|
|
|
1877 $coord_end_chr1_cluster1 = $coord_end_chr1_cluster2;
|
|
|
1878 $coord_start_chr2_cluster1 = $coord_start_chr2_cluster2;
|
|
|
1879 $coord_end_chr2_cluster1 = $coord_end_chr2_cluster2;
|
|
|
1880 $clusterCoordinates_chr1 = '('.$coord_start_chr1_cluster1.','.$coord_end_chr1_cluster1.')';
|
|
|
1881 $clusterCoordinates_chr2 = '('.$coord_start_chr2_cluster1.','.$coord_end_chr2_cluster1.')';
|
|
|
1882 }
|
|
|
1883 $reads_left = scalar keys %{$table};
|
|
|
1884 }
|
|
|
1885 if ($ifBalanced eq 'BAL') {
|
|
|
1886 $ifRenv = $signatureType;
|
|
|
1887 }
|
|
|
1888 }
|
|
|
1889 }
|
|
|
1890 }
|
|
|
1891 }
|
|
|
1892 if ($ifBalanced ne 'BAL') {
|
|
|
1893 #define possible break point
|
|
|
1894 $coord_start_chr1 = $coord_start_chr1_cluster1;
|
|
|
1895 $coord_end_chr1 = $coord_end_chr1_cluster1;
|
|
|
1896 $coord_start_chr2 = $coord_start_chr2_cluster1;
|
|
|
1897 $coord_end_chr2 = $coord_end_chr2_cluster1;
|
|
|
1898
|
|
|
1899 my $region_length_chr1 = $coord_end_chr1-$coord_start_chr1;
|
|
|
1900 my $region_length_chr2 = $coord_end_chr2-$coord_start_chr2;
|
|
|
1901
|
|
|
1902 my $leftLetterOk = substr($mate_sense, 0, 1); #R
|
|
|
1903 my $rightLetterOk = substr($mate_sense, 1, 1); #F
|
|
|
1904
|
|
|
1905 my @index = keys %{$table};
|
|
|
1906 unless ($diff_sense_ends) {
|
|
|
1907 my $firstEndOrder1 = $ends_order1[$index[0]];
|
|
|
1908 my $firstEndOrder2 = $ends_order2[$index[0]];
|
|
|
1909 $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.')';
|
|
|
1910 $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.')';
|
|
|
1911 } else {
|
|
|
1912 $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.')';
|
|
|
1913 $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.')';
|
|
|
1914 }
|
|
|
1915
|
|
|
1916 if ($chr1 ne $chr2){
|
|
|
1917 $ifRenv="INV_TRANSLOC" if($ifRenv eq "REVERSE_SENSE");
|
|
|
1918 $ifRenv="TRANSLOC" if($ifRenv eq "NORMAL_SENSE");
|
|
|
1919 }
|
|
|
1920 }
|
|
|
1921
|
|
|
1922 if (($ifBalanced eq 'BAL')&&( (scalar keys %{$table}) + (scalar keys %{$secondTable}) < $nb_pairs_threshold)) {
|
|
|
1923 next; #discard the link
|
|
|
1924 }
|
|
|
1925 if (($ifBalanced eq 'UNBAL')&&(scalar keys %{$table} < $nb_pairs_threshold)) {
|
|
|
1926 next; #discard the link
|
|
|
1927 }
|
|
|
1928 my $ratioTxt = "$reads_left/".(scalar @pairs);
|
|
|
1929 my ($n1,$nTot) = split ("/",$strand_ratio_filtering);
|
|
|
1930 my $ratioReal = $reads_left/$nTot;
|
|
|
1931
|
|
|
1932 if ($coord_start_chr1<=0) {
|
|
|
1933 $coord_start_chr1=1;
|
|
|
1934 }
|
|
|
1935 if ($coord_start_chr2<=0) {
|
|
|
1936 $coord_start_chr2=1;
|
|
|
1937 }
|
|
|
1938 #create output
|
|
|
1939 my @link=($chr->{$chr1}->{name}, $coord_start_chr1 , $coord_end_chr1, #all information output
|
|
|
1940 $chr->{$chr2}->{name}, $coord_start_chr2 , $coord_end_chr2,
|
|
|
1941 $reads_left,
|
|
|
1942 &redraw(1,$table,$secondTable,\%badInFRSense,$ifBalanced,\@pairs),
|
|
|
1943 &redraw(1,$table,$secondTable,\%badInFRSense,$ifBalanced,\@strand1),
|
|
|
1944 &redraw(1,$table,$secondTable,\%badInFRSense,$ifBalanced,\@strand2),
|
|
|
1945 &redraw(1,$table,$secondTable,\%badInFRSense,$ifBalanced,\@ends_order1),
|
|
|
1946 &redraw(1,$table,$secondTable,\%badInFRSense,$ifBalanced,\@ends_order2),
|
|
|
1947 &redraw(2,$table,$secondTable,\%badInFRSense,$ifBalanced,\@order1),
|
|
|
1948 &redraw(2,$table,$secondTable,\%badInFRSense,$ifBalanced,\@order2),
|
|
|
1949 &redraw(1,$table,$secondTable,\%badInFRSense,$ifBalanced,\@positions1),
|
|
|
1950 &redraw(1,$table,$secondTable,\%badInFRSense,$ifBalanced,\@positions2),
|
|
|
1951 $ifRenv,
|
|
|
1952 $strand_ratio_filtering,
|
|
|
1953 $ifBalanced, $ratioTxt, $break_pont_chr1, $break_pont_chr2,
|
|
|
1954 $ratioReal, $nTot);
|
|
|
1955
|
|
|
1956 $nb_links++;
|
|
|
1957 print $fh2 join("\t",@link)."\n";
|
|
|
1958
|
|
|
1959 if($record>=$warn){
|
|
|
1960 print LOG "-- $fchr : $warn links analysed - $nb_links links kept\n";
|
|
|
1961 $warn+=10000;
|
|
|
1962 }
|
|
|
1963
|
|
|
1964 }
|
|
|
1965 $fh->close;
|
|
|
1966 $fh2->close;
|
|
|
1967
|
|
|
1968 print LOG "-- $fchr : Total : $record links analysed - $nb_links links kept\n";
|
|
|
1969
|
|
|
1970 }
|
|
|
1971 #------------------------------------------------------------------------------#
|
|
|
1972 #------------------------------------------------------------------------------#
|
|
|
1973 #gets information about ends positions given start, direction and order
|
|
|
1974 sub getEnds {
|
|
|
1975 my ($ends,$starts,$strand,$end_order,$tag_length) = @_;
|
|
|
1976 for my $i (0..scalar(@{$starts})-1) {
|
|
|
1977 $ends->[$i] = getEnd($starts->[$i],$strand->[$i],$end_order->[$i],$tag_length);
|
|
|
1978 }
|
|
|
1979 }
|
|
|
1980 sub getEnd {
|
|
|
1981 my ($start,$strand, $end_order,$tag_length) = @_;
|
|
|
1982 return ($strand eq 'F')? $start+$tag_length->{$end_order}-1:$start-$tag_length->{$end_order}+1;
|
|
|
1983 }
|
|
|
1984 #------------------------------------------------------------------------------#
|
|
|
1985 #------------------------------------------------------------------------------#
|
|
|
1986 #gets starts and ends Coords when start=leftmost given positions, directions and orders
|
|
|
1987 sub getCoordswithLeftMost {
|
|
|
1988
|
|
|
1989 my ($starts,$ends,$positions,$strand,$end_order,$tag_length) = @_;
|
|
|
1990
|
|
|
1991 for my $i (0..scalar(@{$positions})-1) {
|
|
|
1992
|
|
|
1993 if($strand->[$i] eq 'F'){
|
|
|
1994 $starts->[$i]=$positions->[$i];
|
|
|
1995 $ends->[$i]=$positions->[$i]+$tag_length->{$end_order->[$i]}-1;
|
|
|
1996 }else{
|
|
|
1997 $starts->[$i]=$positions->[$i]-$tag_length->{$end_order->[$i]}+1;
|
|
|
1998 $ends->[$i]=$positions->[$i];
|
|
|
1999 }
|
|
|
2000 }
|
|
|
2001 }
|
|
|
2002 #------------------------------------------------------------------------------#
|
|
|
2003 #------------------------------------------------------------------------------#
|
|
|
2004 sub addInsertionInfo { #add field with INS,DEL,NA and distance between clusters and performs filtering
|
|
|
2005
|
|
|
2006 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)=@_;
|
|
|
2007
|
|
|
2008 my @sfile=split(/\./,$links_file);
|
|
|
2009 my $fchr=$sfile[$#sfile-2];
|
|
|
2010
|
|
|
2011
|
|
|
2012 my $diff_sense_ends=(($mate_sense eq "FR") || ($mate_sense eq "RF"))? 1:0;
|
|
|
2013
|
|
|
2014 my $record=0;
|
|
|
2015 my $nb_links=0;
|
|
|
2016 my $warn=10000;
|
|
|
2017
|
|
|
2018 print LOG "\# Filtering out normal pairs using insert size...\n";
|
|
|
2019 print LOG "-- mu length=$mu, sigma length=$sigma, indel sigma threshold=$indel_sigma_threshold, dup sigma threshold=$dup_sigma_threshold\n";
|
|
|
2020 print LOG "-- using ".($mu-$indel_sigma_threshold*$sigma)."-".
|
|
|
2021 ($mu+$indel_sigma_threshold*$sigma)." as normal range of insert size for indels\n";
|
|
|
2022 print LOG "-- using ".($mu-$dup_sigma_threshold*$sigma)."-".
|
|
|
2023 ($mu+$dup_sigma_threshold*$sigma)." as normal range of insert size for duplications\n";
|
|
|
2024 print LOG "-- using ".($mu-$singleton_sigma_threshold*$sigma)." as the upper limit of insert size for singletons\n" if($mate_sense eq "RF");
|
|
|
2025
|
|
|
2026 my $fh = new FileHandle;
|
|
|
2027 my $fh2 = new FileHandle;
|
|
|
2028
|
|
|
2029 $fh->open("<$links_file") or die "$0: can't open $links_file :$!\n";
|
|
|
2030 $fh2->open(">$flinks_file") or die "$0: can't write in the output: $flinks_file :$!\n";
|
|
|
2031
|
|
|
2032 while(<$fh>){
|
|
|
2033
|
|
|
2034 $record++;
|
|
|
2035 my @t = split;
|
|
|
2036 my ($chr1,$chr2,$mates_list)=@t[0,3,7];
|
|
|
2037
|
|
|
2038 if($chrID->{$chr1} ne $chrID->{$chr2}) { #if inter-chromosomal link here (because sv_type=all),
|
|
|
2039 $nb_links++;
|
|
|
2040
|
|
|
2041 $t[16]="INV_TRANSLOC" if($t[16] eq "REVERSE_SENSE");
|
|
|
2042 $t[16]="TRANSLOC" if($t[16] eq "NORMAL_SENSE");
|
|
|
2043
|
|
|
2044 $t[16].= "\t";
|
|
|
2045 $t[19].= "\t";
|
|
|
2046
|
|
|
2047 print $fh2 join("\t",@t)."\n";
|
|
|
2048
|
|
|
2049 if($record>=$warn){
|
|
|
2050 print LOG "-- $fchr : $warn links processed - $nb_links links kept\n";
|
|
|
2051 $warn+=10000;
|
|
|
2052 }
|
|
|
2053 next;
|
|
|
2054 }
|
|
|
2055
|
|
|
2056 my $ifRenv = $t[16];
|
|
|
2057 my $ifBalanced = "UNBAL";
|
|
|
2058 $ifBalanced = $t[18] if ($order_filtering);
|
|
|
2059
|
|
|
2060 my $numberOfPairs = $t[6];
|
|
|
2061 my @positions1 = deleteBadOrderSensePairs(split (/,/,$t[14]));
|
|
|
2062 my @positions2 = deleteBadOrderSensePairs(split (/,/,$t[15]));
|
|
|
2063
|
|
|
2064 if ($ifBalanced eq "BAL") {
|
|
|
2065
|
|
|
2066 if ($ifRenv eq "INV_TRANSLOC") {
|
|
|
2067 $ifRenv = "INV_FRAGMT"; #for intrachromosomal inverted translocation is the same as inverted fragment
|
|
|
2068 }
|
|
|
2069 if ($ifRenv eq "NORMAL_SENSE") {
|
|
|
2070 $ifRenv = "TRANSLOC";
|
|
|
2071 }
|
|
|
2072 if ($ifRenv eq "REVERSE_SENSE") {
|
|
|
2073 $ifRenv = "INV_FRAGMT"; #for intrachromosomal inverted translocation is the same as inverted fragment
|
|
|
2074 }
|
|
|
2075 $t[19].= "\t";
|
|
|
2076
|
|
|
2077 my $meanDistance = 0;
|
|
|
2078
|
|
|
2079 for my $i (0..$numberOfPairs-1) {
|
|
|
2080 $meanDistance += $positions2[$i]-$positions1[$i];
|
|
|
2081 }
|
|
|
2082 $meanDistance /= $numberOfPairs;
|
|
|
2083
|
|
|
2084 $t[16] = $ifRenv."\t".$meanDistance;
|
|
|
2085 #dont touch the annotation. It should be already OK.
|
|
|
2086
|
|
|
2087 } else {
|
|
|
2088 #only for unbalanced
|
|
|
2089
|
|
|
2090 my $ifoverlap=overlap($t[1],$t[2],$t[4],$t[5]);
|
|
|
2091
|
|
|
2092 my $ends_sense_class = (deleteBadOrderSensePairs(split (/,/,$t[8])))[0].
|
|
|
2093 (deleteBadOrderSensePairs(split (/,/,$t[9])))[0];
|
|
|
2094 my $ends_order_class = (deleteBadOrderSensePairs(split (/,/,$t[10])))[0].
|
|
|
2095 (deleteBadOrderSensePairs(split (/,/,$t[11])))[0];
|
|
|
2096
|
|
|
2097 my $indel_type = $ifRenv;
|
|
|
2098
|
|
|
2099 my $meanDistance = "N/A";
|
|
|
2100
|
|
|
2101 ($meanDistance, $indel_type) = checkIndel ($numberOfPairs, #identify insertion type for rearrangments without inversion, calculates distance between cluster
|
|
|
2102 \@positions1, #assign N/A to $indel_type if unknown
|
|
|
2103 \@positions2,
|
|
|
2104 $ifRenv,
|
|
|
2105 $ifoverlap,
|
|
|
2106 $indel_sigma_threshold,
|
|
|
2107 $dup_sigma_threshold,
|
|
|
2108 $singleton_sigma_threshold,
|
|
|
2109 $mu,
|
|
|
2110 $sigma,
|
|
|
2111 $ifBalanced,
|
|
|
2112 $ends_sense_class,
|
|
|
2113 $ends_order_class,
|
|
|
2114 $mate_sense,
|
|
|
2115 $diff_sense_ends,
|
|
|
2116 );
|
|
|
2117
|
|
|
2118 #filtering of pairs with distance inconsistant with the SV
|
|
|
2119 if ($ifRenv ne "REVERSE_SENSE") {
|
|
|
2120 my $maxCoord1 =$chr->{$chrID->{$chr1}}->{length};
|
|
|
2121 my $maxCoord2 =$chr->{$chrID->{$chr2}}->{length};
|
|
|
2122 $meanDistance = recalc_t_usingInsertSizeInfo(\@t,$mu,$sigma,$meanDistance,$tag_length,$diff_sense_ends,$mate_sense,
|
|
|
2123 $maxCoord1,$maxCoord2,$ends_sense_class,$ends_order_class,$nb_pairs_threshold,$order_filtering);
|
|
|
2124 next if ($t[6] < $nb_pairs_threshold);
|
|
|
2125 }else{
|
|
|
2126 $t[19].= "\t";
|
|
|
2127 }
|
|
|
2128 $t[16] = $indel_type."\t".$meanDistance;
|
|
|
2129 }
|
|
|
2130
|
|
|
2131 $nb_links++;
|
|
|
2132
|
|
|
2133 print $fh2 join("\t",@t)."\n";
|
|
|
2134 if($record>=$warn){
|
|
|
2135 print LOG "-- $fchr : $warn links processed - $nb_links links kept\n";
|
|
|
2136 $warn+=10000;
|
|
|
2137 }
|
|
|
2138 }
|
|
|
2139 $fh->close;
|
|
|
2140 $fh2->close;
|
|
|
2141
|
|
|
2142 print LOG "-- $fchr : Total : $record links analysed - $nb_links links kept\n";
|
|
|
2143
|
|
|
2144 }
|
|
|
2145 #------------------------------------------------------------------------------#
|
|
|
2146 #------------------------------------------------------------------------------#
|
|
|
2147 sub checkIndel {
|
|
|
2148
|
|
|
2149 my ($numberOfPairs, $positions1, $positions2, $ifRenv, $ifoverlap, $indel_sigma_threshold, $dup_sigma_threshold, $singleton_sigma_threshold,
|
|
|
2150 $mu, $sigma, $ifBalanced,$ends_sense_class,$ends_order_class,$mate_sense,$diff_sense_ends) = @_;
|
|
|
2151
|
|
|
2152 my $meanDistance = 0;
|
|
|
2153
|
|
|
2154 for my $i (0..$numberOfPairs-1) {
|
|
|
2155 $meanDistance += $positions2->[$i]-$positions1->[$i];
|
|
|
2156 }
|
|
|
2157 $meanDistance /= $numberOfPairs;
|
|
|
2158
|
|
|
2159 return ($meanDistance,"INV_DUPLI") if (($ifRenv eq "REVERSE_SENSE") && ($meanDistance<$mu+$dup_sigma_threshold*$sigma) );
|
|
|
2160
|
|
|
2161 return ($meanDistance,"INVERSION") if ($ifRenv eq "REVERSE_SENSE");
|
|
|
2162
|
|
|
2163 if($diff_sense_ends){
|
|
|
2164 return ($meanDistance, "LARGE_DUPLI") if ($ends_sense_class ne $mate_sense) && ($meanDistance>$mu+$dup_sigma_threshold*$sigma) ;
|
|
|
2165 return ($meanDistance, "SINGLETON") if (($meanDistance<$mu-$singleton_sigma_threshold*$sigma) && $mate_sense eq "RF" && ($ends_sense_class eq inverseSense($mate_sense)));
|
|
|
2166 }else{
|
|
|
2167 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")) &&
|
|
|
2168 ($meanDistance>$mu+$dup_sigma_threshold*$sigma) ;
|
|
|
2169 }
|
|
|
2170
|
|
|
2171 return ($meanDistance, "SMALL_DUPLI") if (($meanDistance<$mu-$dup_sigma_threshold*$sigma) && $ifoverlap);
|
|
|
2172
|
|
|
2173 return ($meanDistance, "DUPLICATION") if ($diff_sense_ends && ($ends_sense_class ne $mate_sense) && ($meanDistance<$mu-$dup_sigma_threshold*$sigma) ) ;
|
|
|
2174
|
|
|
2175 return ($meanDistance, "INSERTION") if ($meanDistance<$mu -$indel_sigma_threshold*$sigma);
|
|
|
2176 return ($meanDistance, "DELETION") if ($meanDistance>$mu+$indel_sigma_threshold*$sigma);
|
|
|
2177
|
|
|
2178 return ($meanDistance, "UNDEFINED");
|
|
|
2179 }
|
|
|
2180 #------------------------------------------------------------------------------#
|
|
|
2181 #------------------------------------------------------------------------------#
|
|
|
2182 #sub reacalulate @t so that get rid of unconsistent pairs (unconsistent insert size )
|
|
|
2183 sub recalc_t_usingInsertSizeInfo {
|
|
|
2184 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) = @_;
|
|
|
2185
|
|
|
2186 my @badPairs;
|
|
|
2187
|
|
|
2188 my @positions1 = getAllEntries($t->[14]);
|
|
|
2189 my @positions2 = getAllEntries($t->[15]);
|
|
|
2190
|
|
|
2191 if ($meanDistance < $mu) {
|
|
|
2192 for my $i (0..scalar(@positions1)-1) {
|
|
|
2193 if (substr($positions2[$i],-1,1) ne '$' && substr($positions2[$i],-1,1) ne '*' && $positions2[$i]-$positions1[$i]>=$mu) {
|
|
|
2194 push(@badPairs,$i);
|
|
|
2195 }
|
|
|
2196 }
|
|
|
2197 } else {
|
|
|
2198 for my $i (0..scalar(@positions1)-1) {
|
|
|
2199 if (substr($positions2[$i],-1,1) ne '$' && substr($positions2[$i],-1,1) ne '*' && $positions2[$i]-$positions1[$i]<=$mu) {
|
|
|
2200 push(@badPairs,$i);
|
|
|
2201 }
|
|
|
2202 }
|
|
|
2203 }
|
|
|
2204
|
|
|
2205 if (scalar (@badPairs)>0) {
|
|
|
2206 #print join("\t",@badPairs).": ".join("\t",@t)."\n";
|
|
|
2207 #remove these inconsistant links
|
|
|
2208 $t->[6] -= scalar(@badPairs); #numberOfPairs
|
|
|
2209 return if ($t->[6] < $nb_pairs_threshold);
|
|
|
2210
|
|
|
2211 $t->[7] = mark_values(\@badPairs, $t->[7]);
|
|
|
2212 $t->[8] = mark_values(\@badPairs, $t->[8]);
|
|
|
2213 $t->[9] = mark_values(\@badPairs, $t->[9]);
|
|
|
2214 $t->[10] = mark_values(\@badPairs, $t->[10]);
|
|
|
2215 $t->[11] = mark_values(\@badPairs, $t->[11]);
|
|
|
2216
|
|
|
2217 $t->[12] = mark_indexes(\@badPairs, $t->[12]);
|
|
|
2218 $t->[13] = mark_indexes(\@badPairs, $t->[13]);
|
|
|
2219
|
|
|
2220 $t->[14] = mark_values(\@badPairs, $t->[14]);
|
|
|
2221 $t->[15] = mark_values(\@badPairs, $t->[15]);
|
|
|
2222 $t->[19] = recalculate_ratio($t->[6],$t->[19]) if ($order_filtering); #add the second ratio
|
|
|
2223 $t->[17] = recalculate_ratio($t->[6],$t->[17]) unless ($order_filtering);
|
|
|
2224 ($t->[1],$t->[2]) = recalculate_boundaries($t->[14],$t->[8],$t->[10],$tag_length);
|
|
|
2225 ($t->[4],$t->[5]) = recalculate_boundaries($t->[15],$t->[9],$t->[11],$tag_length);
|
|
|
2226
|
|
|
2227 #recalc breakpoints:
|
|
|
2228 my $quant001 = 3.090232;
|
|
|
2229 my $maxFragmentLength = &floor($quant001 * $sigma + $mu);
|
|
|
2230 $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);
|
|
|
2231 $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);
|
|
|
2232 #recalc total ratio
|
|
|
2233 $t->[22] = $t->[6] / $t->[23] if ($order_filtering);
|
|
|
2234 $t->[18] = $t->[6] / $t->[19] unless ($order_filtering);
|
|
|
2235
|
|
|
2236 @positions1 = deleteBadOrderSensePairs(split (/,/,$t->[14]));
|
|
|
2237 @positions2 = deleteBadOrderSensePairs(split (/,/,$t->[15]));
|
|
|
2238
|
|
|
2239 $meanDistance = 0;
|
|
|
2240
|
|
|
2241 for my $i (0..scalar(@positions1)-1) {
|
|
|
2242 $meanDistance += $positions2[$i]-$positions1[$i];
|
|
|
2243 }
|
|
|
2244 $meanDistance /= scalar(@positions1);
|
|
|
2245
|
|
|
2246 } else {
|
|
|
2247 $t->[17] = recalculate_ratio((split(/\//,$t->[17]))[0],$t->[17]) unless ($order_filtering);
|
|
|
2248 $t->[19] = recalculate_ratio((split(/\//,$t->[19]))[0],$t->[19]) if ($order_filtering);
|
|
|
2249
|
|
|
2250 } #nothing has been filtered
|
|
|
2251 return $meanDistance;
|
|
|
2252 }
|
|
|
2253
|
|
|
2254 sub recalculate_ratio {
|
|
|
2255 my ($left, $ratio) = @_;
|
|
|
2256 my @elements = split (/\//,$ratio);
|
|
|
2257 $elements[1]= $elements[0];
|
|
|
2258 $elements[0]=$left;
|
|
|
2259 return $ratio."\t".join("/",@elements);
|
|
|
2260 }
|
|
|
2261
|
|
|
2262 sub recalc_breakpoints {
|
|
|
2263 my ($mate_sense,$maxCoord,$startString,$strand,$firstEndOrder,$coord_start_chr,$coord_end_chr,$maxFragmentLength,$diff_sense_ends ) = @_;
|
|
|
2264 my $break_pont_chr;
|
|
|
2265
|
|
|
2266 my $leftLetterOk = substr($mate_sense, 0, 1); #R
|
|
|
2267 my $rightLetterOk = substr($mate_sense, 1, 1); #F
|
|
|
2268
|
|
|
2269
|
|
|
2270 my @positions = deleteBadOrderSensePairs(split (/,/,$startString));
|
|
|
2271
|
|
|
2272 unless ($diff_sense_ends) {
|
|
|
2273 $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.')';
|
|
|
2274 } else {
|
|
|
2275 $break_pont_chr = ($strand eq $leftLetterOk)?'('.$coord_end_chr.','.min(($coord_start_chr+$maxFragmentLength),$maxCoord).')':'('.max(($coord_end_chr-$maxFragmentLength),1).','.$coord_start_chr.')';
|
|
|
2276 }
|
|
|
2277 return $break_pont_chr;
|
|
|
2278 }
|
|
|
2279 sub recalculate_boundaries {
|
|
|
2280 my ($startString,$senseString,$endsOrderString,$tag_length) = @_;
|
|
|
2281 my @positions = deleteBadOrderSensePairs(split (/,/,$startString));
|
|
|
2282 my @strands = deleteBadOrderSensePairs(split (/,/,$senseString));
|
|
|
2283 my @ends_orders = deleteBadOrderSensePairs(split (/,/,$endsOrderString));
|
|
|
2284 my @ends; getEnds(\@ends,\@positions,\@strands,\@ends_orders,$tag_length);
|
|
|
2285 my $coord_start_cluster = min(min(@positions),min(@ends));
|
|
|
2286 my $coord_end_cluster = max(max(@positions),max(@ends));
|
|
|
2287 return ($coord_start_cluster,$coord_end_cluster);
|
|
|
2288 }
|
|
|
2289
|
|
|
2290 sub remove_indexes {
|
|
|
2291 my ($bads, $string) = @_;
|
|
|
2292 my @elements = deleteBadOrderSensePairs(split (/,/,$string));
|
|
|
2293 for my $i (reverse sort %{$bads}) {
|
|
|
2294 delete $elements[$i];
|
|
|
2295 }
|
|
|
2296 return "(".join(",",@elements).")";
|
|
|
2297 }
|
|
|
2298 ##add @ to to elements
|
|
|
2299 sub mark_values {
|
|
|
2300 my ($bads, $string) = @_;
|
|
|
2301 my @elements = getAllEntries($string);
|
|
|
2302 for my $i (@{$bads}) {
|
|
|
2303 $elements[$i] .= "@";
|
|
|
2304 }
|
|
|
2305 return "(".join(",",@elements).")";
|
|
|
2306 }
|
|
|
2307 ##add @ to to indexes
|
|
|
2308 sub mark_indexes {
|
|
|
2309 my ($bads, $string) = @_;
|
|
|
2310 my @elements = getAllEntries($string);
|
|
|
2311 for my $i ((0..scalar(@elements)-1)) {
|
|
|
2312 for my $j (@{$bads}) {
|
|
|
2313 $elements[$i] .= "@" if ($elements[$i] eq ($j+1));
|
|
|
2314 }
|
|
|
2315 }
|
|
|
2316
|
|
|
2317 return "(".join(",",@elements).")";
|
|
|
2318 }
|
|
|
2319
|
|
|
2320 #------------------------------------------------------------------------------#
|
|
|
2321 #------------------------------------------------------------------------------#
|
|
|
2322 sub redraw {
|
|
|
2323
|
|
|
2324 my ($type,$table,$secondTable,$badInFRSense,$ifBalanced,$arr) = @_;
|
|
|
2325
|
|
|
2326 my $out;
|
|
|
2327 my @first_arr;
|
|
|
2328 if ($ifBalanced eq 'BAL') {
|
|
|
2329 my @second_arr;
|
|
|
2330 my $lastPushed = 1;
|
|
|
2331 if ($type == 1) {
|
|
|
2332 for my $i (0 .. scalar(@{$arr})-1) {
|
|
|
2333 if (exists ($table->{$i})) {
|
|
|
2334 push(@first_arr,$arr->[$i]);
|
|
|
2335 $lastPushed = 1;
|
|
|
2336 }elsif (exists ($secondTable->{$i})) {
|
|
|
2337 push(@second_arr,$arr->[$i]);
|
|
|
2338 $lastPushed = 2;
|
|
|
2339 } elsif ($lastPushed == 1) {
|
|
|
2340 if (exists ($badInFRSense->{$i})) {
|
|
|
2341 push(@first_arr,$arr->[$i]."\$");
|
|
|
2342 }else {
|
|
|
2343 push(@first_arr,$arr->[$i]."*");
|
|
|
2344 }
|
|
|
2345 } elsif ($lastPushed == 2) {
|
|
|
2346 if (exists ($badInFRSense->{$i})) {
|
|
|
2347 push(@second_arr,$arr->[$i]."\$");
|
|
|
2348 }else {
|
|
|
2349 push(@second_arr,$arr->[$i]."*");
|
|
|
2350 }
|
|
|
2351 } else {print "Error!";exit;}
|
|
|
2352 }
|
|
|
2353 } else {
|
|
|
2354 for my $i (@{$arr}) {
|
|
|
2355 if (exists ($table->{$i-1})) {
|
|
|
2356 push(@first_arr,$i);
|
|
|
2357 $lastPushed = 1;
|
|
|
2358 }elsif (exists ($secondTable->{$i-1})) {
|
|
|
2359 push(@second_arr,$i);
|
|
|
2360 $lastPushed = 2;
|
|
|
2361 } elsif ($lastPushed == 1) {
|
|
|
2362 if (exists ($badInFRSense->{$i-1})) {
|
|
|
2363 push(@first_arr,$i."\$");
|
|
|
2364 }else {
|
|
|
2365 push(@first_arr,$i."*");
|
|
|
2366 }
|
|
|
2367 } elsif ($lastPushed == 2) {
|
|
|
2368 if (exists ($badInFRSense->{$i-1})) {
|
|
|
2369 push(@second_arr,$i."\$");
|
|
|
2370 }else {
|
|
|
2371 push(@second_arr,$i."*");
|
|
|
2372 }
|
|
|
2373 } else {print "Error!";exit;}
|
|
|
2374 }
|
|
|
2375 }
|
|
|
2376 $out = '('.join(",",@first_arr).'),('.join(",",@second_arr).')';
|
|
|
2377 }
|
|
|
2378 else {
|
|
|
2379 if ($type == 1) {
|
|
|
2380 for my $i (0 .. scalar(@{$arr})-1) {
|
|
|
2381 if (exists ($table->{$i})) {
|
|
|
2382 push(@first_arr,$arr->[$i]);
|
|
|
2383 } else {
|
|
|
2384 if (exists ($badInFRSense->{$i})) {
|
|
|
2385 push(@first_arr,$arr->[$i]."\$");
|
|
|
2386 }else {
|
|
|
2387 push(@first_arr,$arr->[$i]."*");
|
|
|
2388 }
|
|
|
2389 }
|
|
|
2390 }
|
|
|
2391 } else {
|
|
|
2392 for my $i (@{$arr}) {
|
|
|
2393 if (exists ($table->{$i-1})) {
|
|
|
2394 push(@first_arr,$i);
|
|
|
2395 } else {
|
|
|
2396 if (exists ($badInFRSense->{$i-1})) {
|
|
|
2397 push(@first_arr,$i."\$");
|
|
|
2398 }else {
|
|
|
2399 push(@first_arr,$i."*");
|
|
|
2400 }
|
|
|
2401 }
|
|
|
2402 }
|
|
|
2403 }
|
|
|
2404 $out = '('.join(",",@first_arr).')';
|
|
|
2405 }
|
|
|
2406 return $out;
|
|
|
2407 }
|
|
|
2408 #------------------------------------------------------------------------------#
|
|
|
2409 #------------------------------------------------------------------------------#
|
|
|
2410 sub check {
|
|
|
2411
|
|
|
2412 my $table = $_[0];
|
|
|
2413 my $bad = 'OK';
|
|
|
2414 my $max = 0;
|
|
|
2415 for my $i (sort {$a<=>$b} keys %{$table}) {
|
|
|
2416 unless ($table->{$i}->{nonAdeq} == 0) {
|
|
|
2417 if ($max<$table->{$i}->{nonAdeq}) {
|
|
|
2418 $max=$table->{$i}->{nonAdeq};
|
|
|
2419 $bad = $i;
|
|
|
2420 }
|
|
|
2421 }
|
|
|
2422 }
|
|
|
2423 return $bad;
|
|
|
2424 }
|
|
|
2425 #------------------------------------------------------------------------------#
|
|
|
2426 #------------------------------------------------------------------------------#
|
|
|
2427 sub reversed {
|
|
|
2428
|
|
|
2429 my ($i,$j,$ifRenv,$positions) = @_;
|
|
|
2430 if (($ifRenv eq 'REVERSE_SENSE' && $positions->[$i]<$positions->[$j]) || ($ifRenv ne 'REVERSE_SENSE' && $positions->[$i]>$positions->[$j])){
|
|
|
2431 return 1;
|
|
|
2432 }
|
|
|
2433 return 0;
|
|
|
2434 }
|
|
|
2435 #------------------------------------------------------------------------------#
|
|
|
2436 #------------------------------------------------------------------------------#
|
|
|
2437 sub remove {
|
|
|
2438
|
|
|
2439 my ($bad,$table) = @_;
|
|
|
2440 for my $i (sort {$a<=>$b} keys %{$table}) {
|
|
|
2441 if ($bad == $i) {
|
|
|
2442 delete($table->{$i});;
|
|
|
2443 } else {
|
|
|
2444 if (exists($table->{$i}->{$bad})) {
|
|
|
2445 delete($table->{$i}->{$bad});
|
|
|
2446 $table->{$i}->{nonAdeq}--;
|
|
|
2447 }
|
|
|
2448 }
|
|
|
2449 }
|
|
|
2450 }
|
|
|
2451 #------------------------------------------------------------------------------#
|
|
|
2452 #------------------------------------------------------------------------------#
|
|
|
2453 sub findBadInRFSenseSOLiDSolexa { #choose maximum: FFFFs or RRRRs
|
|
|
2454
|
|
|
2455 my ($strand,$ends_order,$mate_sense,@keysLeft) = @_;
|
|
|
2456
|
|
|
2457 my $leftLetterOk = substr($mate_sense, 0, 1); #R
|
|
|
2458 my $rightLetterOk = substr($mate_sense, 1, 1); #F
|
|
|
2459
|
|
|
2460 my (@standardArray);
|
|
|
2461 if ($leftLetterOk eq $rightLetterOk) { #SOLID mate-pairs
|
|
|
2462 $leftLetterOk = 'R';
|
|
|
2463 $rightLetterOk = 'F';
|
|
|
2464 @standardArray = translateSolidToRF($strand,$ends_order);
|
|
|
2465 } else {
|
|
|
2466 @standardArray = @{$strand};
|
|
|
2467 }
|
|
|
2468
|
|
|
2469 my $ifR = 0;
|
|
|
2470 my @Rs;
|
|
|
2471
|
|
|
2472 for my $i (@keysLeft) {
|
|
|
2473 if ($standardArray[$i] eq $leftLetterOk) {
|
|
|
2474 $ifR++;
|
|
|
2475 push(@Rs,$i);
|
|
|
2476 }
|
|
|
2477 }
|
|
|
2478
|
|
|
2479
|
|
|
2480 my $ifF = 0;
|
|
|
2481 my @Fs;
|
|
|
2482
|
|
|
2483 for my $i (@keysLeft) {
|
|
|
2484 if ($standardArray[$i] eq $rightLetterOk) {
|
|
|
2485 $ifF++;
|
|
|
2486 push(@Fs,$i);
|
|
|
2487 }
|
|
|
2488 }
|
|
|
2489
|
|
|
2490 if($ifR>=$ifF) {
|
|
|
2491 return @Fs;
|
|
|
2492 }
|
|
|
2493 return @Rs;
|
|
|
2494 }
|
|
|
2495
|
|
|
2496 #------------------------------------------------------------------------------#
|
|
|
2497 #------------------------------------------------------------------------------#
|
|
|
2498 sub findBadInFRSenseSOLiDSolexa { #should work both for SOLiD and Solexa
|
|
|
2499
|
|
|
2500 my ($strand1,$strand2,$ends_order1,$ends_order2,$order1,$order2) = ($_[0],$_[1],$_[2],$_[3],$_[4],$_[5]);
|
|
|
2501 my $mate_sense = $_[6];
|
|
|
2502
|
|
|
2503 my $leftLetterOk = substr($mate_sense, 0, 1); #R
|
|
|
2504 my $rightLetterOk = substr($mate_sense, 1, 1); #F
|
|
|
2505
|
|
|
2506 my (@standardArray1,@standardArray2);
|
|
|
2507
|
|
|
2508 if ($leftLetterOk eq $rightLetterOk) { #SOLID mate-pairs
|
|
|
2509 $leftLetterOk = 'R';
|
|
|
2510 $rightLetterOk = 'F';
|
|
|
2511 @standardArray1 = translateSolidToRF($strand1,$ends_order1);
|
|
|
2512 my @arr = getOrderedStrands($strand2,$order2);
|
|
|
2513 my @ends2 = getOrderedStrands($ends_order2,$order2);
|
|
|
2514 @standardArray2 = translateSolidToRF(\@arr,\@ends2);
|
|
|
2515
|
|
|
2516 } else {
|
|
|
2517 @standardArray1 = @{$strand1};
|
|
|
2518 @standardArray2 = getOrderedStrands($strand2,$order2);
|
|
|
2519 }
|
|
|
2520
|
|
|
2521 #we will try 4 possibilities, 2 for each end of the link: RFRR-FFF->RFFFF , RFRR-FFF->RRRFFF
|
|
|
2522
|
|
|
2523 #for the first end:
|
|
|
2524
|
|
|
2525 my @array = @standardArray1;
|
|
|
2526 my %badInFRSense1;
|
|
|
2527 for my $i (1..scalar (@array)-1){ # FRFRFFFF -> FFFFFF and RRFRFRFFFF -> RRFFFFFF
|
|
|
2528 if ($array[$i-1] eq $rightLetterOk && $array[$i] eq $leftLetterOk) {
|
|
|
2529 $badInFRSense1{$i}=1;
|
|
|
2530 $array[$i] = $rightLetterOk;
|
|
|
2531 }
|
|
|
2532 }
|
|
|
2533 my $numberRRRFFF_or_FFF_1 = scalar(@array)-scalar(keys %badInFRSense1);
|
|
|
2534 @array = @standardArray1;
|
|
|
2535 my %badInFRSense0;
|
|
|
2536 for my $i (reverse(1..scalar (@array)-1)){ # FRFRFFFFRR -> FFFFFFRR
|
|
|
2537 if ($array[$i-1] eq $rightLetterOk && $array[$i] eq $leftLetterOk) {
|
|
|
2538 $badInFRSense0{$i-1}=1;
|
|
|
2539 $array[$i-1] = $leftLetterOk;
|
|
|
2540
|
|
|
2541 }
|
|
|
2542 }
|
|
|
2543 my $numberRRF1 = scalar(@array)-scalar(keys %badInFRSense0);
|
|
|
2544
|
|
|
2545 #for the second end:
|
|
|
2546 @array = @standardArray2;
|
|
|
2547
|
|
|
2548 my %badInFRSense3;
|
|
|
2549 for my $i (1..scalar(@array)-1){
|
|
|
2550 if ($array[$i-1] eq $rightLetterOk && $array[$i] eq $leftLetterOk) {
|
|
|
2551 $badInFRSense3{$order2->[$i]}=1;
|
|
|
2552 $array[$i] = $rightLetterOk;
|
|
|
2553 }
|
|
|
2554 }
|
|
|
2555 my $numberRRRFFF_or_FFF_2 = scalar(@array)-scalar(keys %badInFRSense3);
|
|
|
2556
|
|
|
2557 @array = @standardArray2;
|
|
|
2558 my %badInFRSense5;
|
|
|
2559 for my $i (reverse(1..scalar (@array)-1)){ # FRFRFFFF -> FFFFFF
|
|
|
2560 if ($array[$i-1] eq $rightLetterOk && $array[$i] eq $leftLetterOk) {
|
|
|
2561 $badInFRSense5{$i-1}=1;
|
|
|
2562 $array[$i-1] = $leftLetterOk;
|
|
|
2563 }
|
|
|
2564 }
|
|
|
2565 my $numberRRF2 = scalar(@array)-scalar(keys %badInFRSense5);
|
|
|
2566
|
|
|
2567 if ($numberRRF1>=$numberRRRFFF_or_FFF_1 && $numberRRF1 >= $numberRRRFFF_or_FFF_2 && $numberRRF1 >=$numberRRF2) {
|
|
|
2568 return (1,%badInFRSense0);
|
|
|
2569 }
|
|
|
2570
|
|
|
2571 if ($numberRRRFFF_or_FFF_1 >=$numberRRF1 && $numberRRRFFF_or_FFF_1 >= $numberRRRFFF_or_FFF_2 && $numberRRRFFF_or_FFF_1 >= $numberRRF2) {
|
|
|
2572 return (1,%badInFRSense1);
|
|
|
2573 }
|
|
|
2574
|
|
|
2575 if ($numberRRRFFF_or_FFF_2 >= $numberRRF1 && $numberRRRFFF_or_FFF_2 >= $numberRRRFFF_or_FFF_1 && $numberRRRFFF_or_FFF_2 >=$numberRRF2) {
|
|
|
2576 return (2,%badInFRSense3);
|
|
|
2577 }
|
|
|
2578
|
|
|
2579 if ($numberRRF2 >= $numberRRF1 && $numberRRF2 >= $numberRRRFFF_or_FFF_1 && $numberRRF2 >= $numberRRRFFF_or_FFF_2 ) {
|
|
|
2580 return (2,%badInFRSense5);
|
|
|
2581 }
|
|
|
2582
|
|
|
2583 #should not get here:
|
|
|
2584 print STDERR "Error in findBadInFRSenseSOLiDSolexa()!\n";
|
|
|
2585 return (1,%badInFRSense1);
|
|
|
2586 }
|
|
|
2587
|
|
|
2588 sub getOrderedStrands {
|
|
|
2589 my ($strand,$order) = ($_[0],$_[1]);
|
|
|
2590 my @arr;
|
|
|
2591 for my $i (0..scalar(@{$strand})-1) {
|
|
|
2592 push(@arr,$strand->[$order->[$i]-1]);
|
|
|
2593 }
|
|
|
2594 return @arr;
|
|
|
2595 }
|
|
|
2596 #------------------------------------------------------------------------------#
|
|
|
2597 #------------------------------------------------------------------------------#
|
|
|
2598 sub checkClusters {
|
|
|
2599
|
|
|
2600 my ($ifRenv,$coord_start_chr1_cluster1,$coord_start_chr1_cluster2,$coord_start_chr2_cluster1,$coord_start_chr2_cluster2) = @_;
|
|
|
2601 if ($ifRenv eq 'REVERSE_SENSE') {
|
|
|
2602 if ($coord_start_chr1_cluster1 <= $coord_start_chr1_cluster2) {
|
|
|
2603 return ($coord_start_chr2_cluster1 <= $coord_start_chr2_cluster2)?1:0;
|
|
|
2604 }
|
|
|
2605 return ($coord_start_chr2_cluster1 >= $coord_start_chr2_cluster2)?1:0;
|
|
|
2606 }
|
|
|
2607 #if NORM
|
|
|
2608 if ($coord_start_chr1_cluster1 <= $coord_start_chr1_cluster2) {
|
|
|
2609 return ($coord_start_chr2_cluster1 >= $coord_start_chr2_cluster2)?1:0;
|
|
|
2610 }
|
|
|
2611 return ($coord_start_chr2_cluster1 <= $coord_start_chr2_cluster2)?1:0;
|
|
|
2612 }
|
|
|
2613
|
|
|
2614 #------------------------------------------------------------------------------#
|
|
|
2615 #------------------------------------------------------------------------------#
|
|
|
2616 sub translateSolidToRF {
|
|
|
2617 my ($strandArr,$ends_orderArr)=@_;
|
|
|
2618 my @array;
|
|
|
2619 for my $i (0..scalar(@{$strandArr})-1) {
|
|
|
2620 if ($ends_orderArr->[$i]==1 && $strandArr->[$i] eq 'F') {
|
|
|
2621 push(@array,'F');
|
|
|
2622 }
|
|
|
2623 if ($ends_orderArr->[$i]==2 && $strandArr->[$i] eq 'F') {
|
|
|
2624 push(@array,'R');
|
|
|
2625 }
|
|
|
2626 if ($ends_orderArr->[$i]==1 && $strandArr->[$i] eq 'R') {
|
|
|
2627 push(@array,'R');
|
|
|
2628 }
|
|
|
2629 if ($ends_orderArr->[$i]==2 && $strandArr->[$i] eq 'R') {
|
|
|
2630 push(@array,'F');
|
|
|
2631 }
|
|
|
2632 }
|
|
|
2633 return @array;
|
|
|
2634 }
|
|
|
2635
|
|
|
2636 #------------------------------------------------------------------------------#
|
|
|
2637 #------------------------------------------------------------------------------#
|
|
|
2638 #convert the links file to the circos format
|
|
|
2639 sub links2segdup{
|
|
|
2640
|
|
|
2641 my($id,$color_code,$links_file,$segdup_file)=@_;
|
|
|
2642
|
|
|
2643 print LOG "# Converting to the circos format...\n";
|
|
|
2644
|
|
|
2645 tie (my %hcolor,'Tie::IxHash'); #color-code hash table
|
|
|
2646 foreach my $col (keys %{$color_code}){
|
|
|
2647 my ($min_links,$max_links)=split(",",$color_code->{$col});
|
|
|
2648 $hcolor{$col}=[$min_links,$max_links];
|
|
|
2649 }
|
|
|
2650
|
|
|
2651 open LINKS, "<$links_file" or die "$0: can't open $links_file :$!\n";
|
|
|
2652 open SEGDUP, ">$segdup_file" or die "$0: can't write in the output: $segdup_file :$!\n";
|
|
|
2653
|
|
|
2654 my $index=1;
|
|
|
2655 while(<LINKS>){
|
|
|
2656
|
|
|
2657 my ($chr1,$start1,$end1,$chr2,$start2,$end2,$count)=(split)[0,1,2,3,4,5,6];
|
|
|
2658
|
|
|
2659 my $color=getColor($count,\%hcolor,"circos"); #get the color-code according the number of links
|
|
|
2660
|
|
|
2661 print SEGDUP "$index\t$id$chr1\t$start1\t$end1\tcolor=$color\n". #circos output
|
|
|
2662 "$index\t$id$chr2\t$start2\t$end2\tcolor=$color\n";
|
|
|
2663 $index++;
|
|
|
2664 }
|
|
|
2665
|
|
|
2666 close LINKS;
|
|
|
2667 close SEGDUP;
|
|
|
2668 print LOG "-- output created: $segdup_file\n";
|
|
|
2669 }
|
|
|
2670 #------------------------------------------------------------------------------#
|
|
|
2671 #------------------------------------------------------------------------------#
|
|
|
2672 #convert the links file to the bedPE format for BEDTools usage
|
|
|
2673 sub links2bedPElinksfile{
|
|
|
2674
|
|
|
2675 my ($sample,$links_file,$bedpe_file)=@_;
|
|
|
2676
|
|
|
2677 open LINKS, "<$links_file" or die "$0: can't open $links_file :$!\n";
|
|
|
2678 open BEDPE, ">$bedpe_file" or die "$0: can't write in the output: $bedpe_file :$!\n";
|
|
|
2679
|
|
|
2680 my $nb_links=1;
|
|
|
2681
|
|
|
2682 while(<LINKS>){
|
|
|
2683
|
|
|
2684 chomp;
|
|
|
2685 my @t=split("\t",$_);
|
|
|
2686 my ($chr1,$start1,$end1,$chr2,$start2,$end2)=splice(@t,0,6);
|
|
|
2687 my $type=($chr1 eq $chr2)? "INTRA":"INTER";
|
|
|
2688 $type.="_".$t[10];
|
|
|
2689
|
|
|
2690 $start1--; $start2--;
|
|
|
2691
|
|
|
2692 print BEDPE "$chr1\t$start1\t$end1\t$chr2\t$start2\t$end2".
|
|
|
2693 "\t$sample"."_link$nb_links\t$type\t.\t.".
|
|
|
2694 "\t".join("|",@t)."\n";
|
|
|
2695
|
|
|
2696 $nb_links++;
|
|
|
2697 }
|
|
|
2698
|
|
|
2699 close LINKS;
|
|
|
2700 close BEDPE;
|
|
|
2701
|
|
|
2702 }
|
|
|
2703 #------------------------------------------------------------------------------#
|
|
|
2704 #------------------------------------------------------------------------------#
|
|
|
2705 sub bedPElinks2linksfile{
|
|
|
2706
|
|
|
2707 my ($bedpe_file,$links_file)=@_;
|
|
|
2708
|
|
|
2709 open BEDPE, "<$bedpe_file" or die "$0: can't open: $bedpe_file :$!\n";
|
|
|
2710 open LINKS, ">$links_file" or die "$0: can't write in the output $links_file :$!\n";
|
|
|
2711
|
|
|
2712 while(<BEDPE>){
|
|
|
2713
|
|
|
2714 chomp;
|
|
|
2715 my $sample=(split("_",(split("\t",$_))[6]))[0];
|
|
|
2716 my @t1=(split("\t",$_))[0,1,2,3,4,5];
|
|
|
2717 my @t2=split(/\|/,(split("\t",$_))[10]);
|
|
|
2718 push(@t2,$sample);
|
|
|
2719
|
|
|
2720 print LINKS join("\t",@t1)."\t".join("\t",@t2)."\n";
|
|
|
2721
|
|
|
2722 }
|
|
|
2723 close BEDPE;
|
|
|
2724 close LINKS;
|
|
|
2725
|
|
|
2726 }
|
|
|
2727 #------------------------------------------------------------------------------#
|
|
|
2728 #------------------------------------------------------------------------------#
|
|
|
2729 #convert the links file to the bed format
|
|
|
2730 sub links2bedfile{
|
|
|
2731
|
|
|
2732 my ($tag_length,$color_code,$links_file,$bed_file)=@_;
|
|
|
2733
|
|
|
2734 print LOG "# Converting to the bed format...\n";
|
|
|
2735
|
|
|
2736 my $compare=1;
|
|
|
2737 if($links_file!~/compared$/){
|
|
|
2738 $compare=0;
|
|
|
2739 $tag_length->{none}->{1}=$tag_length->{1};
|
|
|
2740 $tag_length->{none}->{2}=$tag_length->{2};
|
|
|
2741 }
|
|
|
2742
|
|
|
2743 #color-code hash table
|
|
|
2744 tie (my %hcolor,'Tie::IxHash');
|
|
|
2745 my %color_order;
|
|
|
2746 my $n=1;
|
|
|
2747 foreach my $col (keys %{$color_code}){
|
|
|
2748 my ($min_links,$max_links)=split(",",$color_code->{$col});
|
|
|
2749 $hcolor{$col}=[$min_links,$max_links];
|
|
|
2750 $color_order{$col}=$n;
|
|
|
2751 $n++;
|
|
|
2752 }
|
|
|
2753
|
|
|
2754 my %pair;
|
|
|
2755 my %pt;
|
|
|
2756 $n=1;
|
|
|
2757 open LINKS, "<$links_file" or die "$0: can't open $links_file:$!\n";
|
|
|
2758
|
|
|
2759 my %str=( "F"=>"+", "R"=>"-" );
|
|
|
2760
|
|
|
2761 while(<LINKS>){
|
|
|
2762
|
|
|
2763 my @t=split;
|
|
|
2764 my $sample=($compare)? pop(@t):"none";
|
|
|
2765
|
|
|
2766 my $chr1=$t[0];
|
|
|
2767 my $chr2=$t[3];
|
|
|
2768 $chr1 = "chr".$chr1 unless ($chr1 =~ m/chr/i);
|
|
|
2769 $chr2 = "chr".$chr2 unless ($chr2 =~ m/chr/i);
|
|
|
2770 my $same_chr=($chr1 eq $chr2)? 1:0;
|
|
|
2771
|
|
|
2772 my $count=$t[6];
|
|
|
2773 my $color=getColor($count,\%hcolor,"bed");
|
|
|
2774
|
|
|
2775 my @pairs=deleteBadOrderSensePairs(split(",",$t[7]));
|
|
|
2776 my @strand1=deleteBadOrderSensePairs(split(",",$t[8]));
|
|
|
2777 my @strand2=deleteBadOrderSensePairs(split(",",$t[9]));
|
|
|
2778 my @ends_order1=deleteBadOrderSensePairs(split(",",$t[10]));
|
|
|
2779 my @ends_order2=deleteBadOrderSensePairs(split(",",$t[11]));
|
|
|
2780 my @position1=deleteBadOrderSensePairs(split(",",$t[14]));
|
|
|
2781 my @position2=deleteBadOrderSensePairs(split(",",$t[15]));
|
|
|
2782 my @start1; my @end1; getCoordswithLeftMost(\@start1,\@end1,\@position1,\@strand1,\@ends_order1,$tag_length->{$sample});
|
|
|
2783 my @start2; my @end2; getCoordswithLeftMost(\@start2,\@end2,\@position2,\@strand2,\@ends_order2,$tag_length->{$sample});
|
|
|
2784
|
|
|
2785
|
|
|
2786 for my $p (0..$#pairs){
|
|
|
2787
|
|
|
2788 if (!exists $pair{$pairs[$p]}){
|
|
|
2789
|
|
|
2790 if($same_chr){
|
|
|
2791
|
|
|
2792 $pair{$pairs[$p]}->{0}=[ $chr1, $start1[$p]-1, $end2[$p], $pairs[$p], 0, $str{$strand1[$p]},
|
|
|
2793 $start1[$p]-1, $end2[$p], $color,
|
|
|
2794 2, $tag_length->{$sample}->{$ends_order1[$p]}.",".$tag_length->{$sample}->{$ends_order2[$p]}, "0,".($start2[$p]-$start1[$p]) ];
|
|
|
2795 $pt{$n}=$pair{$pairs[$p]}->{0};
|
|
|
2796 $n++;
|
|
|
2797
|
|
|
2798 }else{
|
|
|
2799
|
|
|
2800 $pair{$pairs[$p]}->{1}=[ $chr1, $start1[$p]-1, $end1[$p] , $pairs[$p]."/1", 0, $str{$strand1[$p]},
|
|
|
2801 $start1[$p]-1, $end1[$p], $color,
|
|
|
2802 1, $tag_length->{$sample}->{$ends_order1[$p]}, 0];
|
|
|
2803 $pt{$n}=$pair{$pairs[$p]}->{1};
|
|
|
2804 $n++;
|
|
|
2805
|
|
|
2806
|
|
|
2807 $pair{$pairs[$p]}->{2}=[ $chr2, $start2[$p]-1, $end2[$p], $pairs[$p]."/2", 0, $str{$strand2[$p]},
|
|
|
2808 $start2[$p]-1, $end2[$p], $color,
|
|
|
2809 1, $tag_length->{$sample}->{$ends_order2[$p]}, 0];
|
|
|
2810 $pt{$n}=$pair{$pairs[$p]}->{2};
|
|
|
2811 $n++;
|
|
|
2812 }
|
|
|
2813 }else{
|
|
|
2814
|
|
|
2815 if($same_chr){
|
|
|
2816 ${$pair{$pairs[$p]}->{0}}[8]=$color if($color_order{$color}>$color_order{${$pair{$pairs[$p]}->{0}}[8]});
|
|
|
2817 }else{
|
|
|
2818 ${$pair{$pairs[$p]}->{1}}[8]=$color if($color_order{$color}>$color_order{${$pair{$pairs[$p]}->{1}}[8]});
|
|
|
2819 ${$pair{$pairs[$p]}->{2}}[8]=$color if($color_order{$color}>$color_order{${$pair{$pairs[$p]}->{2}}[8]});
|
|
|
2820 }
|
|
|
2821 }
|
|
|
2822 }
|
|
|
2823 }
|
|
|
2824 close LINKS;
|
|
|
2825
|
|
|
2826 my $nb_pairs=$n-1;
|
|
|
2827
|
|
|
2828 open BED, ">$bed_file" or die "$0: can't write in the output: $bed_file :$!\n";
|
|
|
2829 print BED "track name=\"$bed_file\" description=\"mate pairs involved in links\" ".
|
|
|
2830 "visibility=2 itemRgb=\"On\"\n";
|
|
|
2831
|
|
|
2832 for my $i (1..$nb_pairs){
|
|
|
2833 print BED join("\t",@{$pt{$i}})."\n";
|
|
|
2834 }
|
|
|
2835
|
|
|
2836 close BED;
|
|
|
2837
|
|
|
2838 print LOG "-- output created: $bed_file\n";
|
|
|
2839
|
|
|
2840 undef %pair;
|
|
|
2841 undef %pt;
|
|
|
2842
|
|
|
2843 }
|
|
|
2844 #------------------------------------------------------------------------------#
|
|
|
2845 #------------------------------------------------------------------------------#
|
|
|
2846 sub deleteBadOrderSensePairs{
|
|
|
2847
|
|
|
2848 my (@tab)=@_;
|
|
|
2849 my @tab2;
|
|
|
2850
|
|
|
2851 foreach my $v (@tab){
|
|
|
2852
|
|
|
2853 $v=~s/[\(\)]//g;
|
|
|
2854 push(@tab2,$v) if($v!~/[\$\*\@]$/);
|
|
|
2855 }
|
|
|
2856 return @tab2;
|
|
|
2857 }
|
|
|
2858 #------------------------------------------------------------------------------#
|
|
|
2859 #------------------------------------------------------------------------------#
|
|
|
2860 sub getAllEntries{
|
|
|
2861 my (@tab)=split (/,/,$_[0]);
|
|
|
2862 my @tab2;
|
|
|
2863
|
|
|
2864 foreach my $v (@tab){
|
|
|
2865
|
|
|
2866 $v=~s/[\(\)]//g;
|
|
|
2867 push(@tab2,$v);
|
|
|
2868 }
|
|
|
2869 return @tab2;
|
|
|
2870 }#------------------------------------------------------------------------------#
|
|
|
2871 #------------------------------------------------------------------------------#
|
|
|
2872 sub getAllEntriesWOspecialChar{
|
|
|
2873 my (@tab)=split (/,/,$_[0]);
|
|
|
2874 my @tab2;
|
|
|
2875
|
|
|
2876 foreach my $v (@tab){
|
|
|
2877
|
|
|
2878 $v=~s/[\(\)\$\*\@]//g;
|
|
|
2879 push(@tab2,$v);
|
|
|
2880 }
|
|
|
2881 return @tab2;
|
|
|
2882 }
|
|
|
2883 #------------------------------------------------------------------------------#
|
|
|
2884 #------------------------------------------------------------------------------#
|
|
|
2885 sub links2SVfile{
|
|
|
2886
|
|
|
2887 my($links_file,$sv_file)=@_;
|
|
|
2888
|
|
|
2889 print LOG "# Converting to the sv output table...\n";
|
|
|
2890 open LINKS, "<$links_file" or die "$0: can't open $links_file:$!\n";
|
|
|
2891 open SV, ">$sv_file" or die "$0: can't write in the output: $sv_file :$!\n";
|
|
|
2892
|
|
|
2893 my @header=qw(chr_type SV_type BAL_type chromosome1 start1-end1 average_dist
|
|
|
2894 chromosome2 start2-end2 nb_pairs score_strand_filtering score_order_filtering score_insert_size_filtering
|
|
|
2895 final_score breakpoint1_start1-end1 breakpoint2_start2-end2);
|
|
|
2896
|
|
|
2897 my $nb_links=0;
|
|
|
2898
|
|
|
2899 while (<LINKS>){
|
|
|
2900
|
|
|
2901 my @t=split;
|
|
|
2902 my @sv=();
|
|
|
2903 my $sv_type="-";
|
|
|
2904 my $strand_ratio="-";
|
|
|
2905 my $eq_ratio="-";
|
|
|
2906 my $eq_type="-";
|
|
|
2907 my $insert_ratio="-";
|
|
|
2908 my $link="-";
|
|
|
2909 my ($bk1, $bk2)=("-","-");
|
|
|
2910 my $score="-";
|
|
|
2911
|
|
|
2912 my ($chr1,$start1,$end1)=($t[0],$t[1],$t[2]);
|
|
|
2913 my ($chr2,$start2,$end2)=($t[3],$t[4],$t[5]);
|
|
|
2914 my $nb_pairs=$t[6];
|
|
|
2915 $chr1 = "chr".$chr1 unless ($chr1 =~ m/chr/i);
|
|
|
2916 $chr2 = "chr".$chr2 unless ($chr2 =~ m/chr/i);
|
|
|
2917 my $chr_type=($chr1 eq $chr2)? "INTRA":"INTER";
|
|
|
2918
|
|
|
2919 #if strand filtering
|
|
|
2920 if (defined $t[16]){
|
|
|
2921 #if inter-chr link
|
|
|
2922 $sv_type=$t[16];
|
|
|
2923 if(defined $t[17] && $t[17]=~/^(\d+)\/(\d+)$/){
|
|
|
2924 $strand_ratio=floor($1/$2*100)."%";
|
|
|
2925 $score=$t[18];
|
|
|
2926 }
|
|
|
2927 if(defined $t[18] && $t[18]=~/^(\d+)\/(\d+)$/){
|
|
|
2928 #if intra-chr link with insert size filtering
|
|
|
2929 $strand_ratio=floor($1/$2*100)."%";
|
|
|
2930 $link=floor($t[17]);
|
|
|
2931 if($sv_type!~/^INV/){
|
|
|
2932 $insert_ratio=floor($1/$2*100)."%" if($t[19]=~/^(\d+)\/(\d+)$/);
|
|
|
2933 $score=$t[20];
|
|
|
2934 }else{
|
|
|
2935 $score=$t[19];
|
|
|
2936 }
|
|
|
2937 }
|
|
|
2938 }
|
|
|
2939
|
|
|
2940 if(defined $t[18] && ($t[18] eq "UNBAL" || $t[18] eq "BAL")){
|
|
|
2941
|
|
|
2942 #if strand and order filtering only and/or interchr link
|
|
|
2943 $eq_type=$t[18];
|
|
|
2944 $eq_ratio=floor($1/$2*100)."%" if($t[19]=~/^(\d+)\/(\d+)$/);
|
|
|
2945 ($bk1, $bk2)=($t[20],$t[21]);
|
|
|
2946 foreach my $bk ($bk1, $bk2){$bk=~s/\),\(/ /g; $bk=~s/(\(|\))//g; $bk=~s/,/-/g;}
|
|
|
2947 $score=$t[22];
|
|
|
2948
|
|
|
2949 }elsif(defined $t[19] && ($t[19] eq "UNBAL" || $t[19] eq "BAL")){
|
|
|
2950
|
|
|
2951 #if all three filtering
|
|
|
2952 $link=floor($t[17]);
|
|
|
2953 $eq_type=$t[19];
|
|
|
2954 $eq_ratio=floor($1/$2*100)."%" if($t[20]=~/^(\d+)\/(\d+)$/);
|
|
|
2955
|
|
|
2956 if(defined $t[21] && $t[21]=~/^(\d+)\/(\d+)$/){
|
|
|
2957 $insert_ratio=floor($1/$2*100)."%";
|
|
|
2958 ($bk1, $bk2)=($t[22],$t[23]);
|
|
|
2959 $score=$t[24];
|
|
|
2960
|
|
|
2961 }else{
|
|
|
2962 ($bk1, $bk2)=($t[21],$t[22]);
|
|
|
2963 $score=$t[23];
|
|
|
2964 }
|
|
|
2965 foreach my $bk ($bk1, $bk2){$bk=~s/\),\(/ /g; $bk=~s/(\(|\))//g; $bk=~s/,/-/g;}
|
|
|
2966
|
|
|
2967 }
|
|
|
2968
|
|
|
2969
|
|
|
2970 push(@sv, $chr_type, $sv_type,$eq_type);
|
|
|
2971 push(@sv,"$chr1\t$start1-$end1");
|
|
|
2972 push(@sv, $link);
|
|
|
2973 push(@sv,"$chr2\t$start2-$end2",
|
|
|
2974 $nb_pairs,$strand_ratio,$eq_ratio,$insert_ratio, decimal($score,4), $bk1, $bk2);
|
|
|
2975
|
|
|
2976
|
|
|
2977 print SV join("\t",@sv)."\n";
|
|
|
2978 }
|
|
|
2979
|
|
|
2980 close LINKS;
|
|
|
2981 close SV;
|
|
|
2982
|
|
|
2983 system "sort -k 9,9nr -k 13,13nr $sv_file > $sv_file.sorted";
|
|
|
2984
|
|
|
2985 open SV, "<".$sv_file.".sorted" or die "$0: can't open in the output: $sv_file".".sorted :$!\n";
|
|
|
2986 my @links=<SV>;
|
|
|
2987 close SV;
|
|
|
2988
|
|
|
2989 open SV, ">$sv_file" or die "$0: can't write in the output: $sv_file :$!\n";
|
|
|
2990
|
|
|
2991 print SV join("\t",@header)."\n";
|
|
|
2992 print SV @links;
|
|
|
2993 close SV;
|
|
|
2994
|
|
|
2995 unlink($sv_file.".sorted");
|
|
|
2996
|
|
|
2997 print LOG "-- output created: $sv_file\n";
|
|
|
2998
|
|
|
2999 }
|
|
|
3000 #------------------------------------------------------------------------------#
|
|
|
3001 sub densityCalculation{
|
|
|
3002
|
|
|
3003 my ($chr,$chrID,$file,$tag_length,$window_dist,$step,$mates_file,$mates_file_ref,$density_file,$input_format)=@_;
|
|
|
3004
|
|
|
3005 my @sfile=split(/\./,$$mates_file[$file]);
|
|
|
3006 my $fchr=$sfile[$#sfile];
|
|
|
3007
|
|
|
3008 my $fh = new FileHandle;
|
|
|
3009
|
|
|
3010 my %density;
|
|
|
3011 my %density_ref;
|
|
|
3012 my @ratio;
|
|
|
3013 my ($cov,$cov_ref);
|
|
|
3014
|
|
|
3015 #FREQUENCY CALCULATION PROCEDURE
|
|
|
3016 print LOG "# $fchr : Frequency calculation procedure...\n";
|
|
|
3017 &FreqCalculation(\%density,$chr,$chrID,$tag_length,$window_dist,$step,$$mates_file[$file],$input_format);
|
|
|
3018 &FreqCalculation(\%density_ref,$chr,$chrID,$tag_length,$window_dist,$step,$$mates_file_ref[$file],$input_format);
|
|
|
3019
|
|
|
3020 #MAKING RATIO AND OUTPUT
|
|
|
3021 print LOG "\# Ratio calculation procedure...\n";
|
|
|
3022 $density_file=~s/\/mates\//\/density\//;
|
|
|
3023 $fh->open(">".$density_file) or die "$0: can't write in the output ".$density_file." :$!\n";
|
|
|
3024
|
|
|
3025 foreach my $k (1..$chr->{nb_chrs}){
|
|
|
3026 foreach my $frag (1..$chr->{$k}->{nb_frag}){
|
|
|
3027
|
|
|
3028 @ratio= ($chr->{$k}->{name},
|
|
|
3029 (${$chr->{$k}->{$frag}}[0]+1),
|
|
|
3030 (${$chr->{$k}->{$frag}}[1]+1));
|
|
|
3031
|
|
|
3032 $cov=(exists $density{$k}{$frag}->{count})? $density{$k}{$frag}->{count}:0;
|
|
|
3033 $cov_ref=(exists $density_ref{$k}{$frag}->{count})? $density_ref{$k}{$frag}->{count}:0;
|
|
|
3034
|
|
|
3035 push(@ratio,$cov,$cov_ref);
|
|
|
3036 push(@ratio,log($cov/$cov_ref)) if($cov && $cov_ref);
|
|
|
3037 push(@ratio,-log($cov_ref+1)) if(!$cov && $cov_ref);
|
|
|
3038 push(@ratio,log($cov+1)) if($cov && !$cov_ref);
|
|
|
3039 next if(!$cov && !$cov_ref);
|
|
|
3040
|
|
|
3041 print $fh join("\t",@ratio)."\n";
|
|
|
3042 }
|
|
|
3043 }
|
|
|
3044
|
|
|
3045 $fh->close;
|
|
|
3046 print LOG "-- output created: $density_file\n";
|
|
|
3047
|
|
|
3048 undef %density;
|
|
|
3049 undef %density_ref;
|
|
|
3050 }
|
|
|
3051 #------------------------------------------------------------------------------#
|
|
|
3052 #------------------------------------------------------------------------------#
|
|
|
3053 sub FreqCalculation{
|
|
|
3054
|
|
|
3055 my ($density,$chr,$chrID,$tag_length,$window_dist,$step,$mates_file,$input_format) = @_;
|
|
|
3056
|
|
|
3057 my @sfile=split(/\./,$mates_file);
|
|
|
3058 my $fchr=$sfile[$#sfile];
|
|
|
3059 my $fh = new FileHandle;
|
|
|
3060
|
|
|
3061 my $nb_windows=0;
|
|
|
3062 my $warn=100000;
|
|
|
3063 my $record=0;
|
|
|
3064 my %pair;
|
|
|
3065
|
|
|
3066 my ($sumX,$sumX2) = (0,0);
|
|
|
3067
|
|
|
3068 print LOG "\# Frequency calculation for $mates_file...\n";
|
|
|
3069
|
|
|
3070 if ($mates_file =~ /.gz$/) {
|
|
|
3071 $fh->open("gunzip -c $mates_file |") or die "$0: can't open ".$mates_file.":$!\n";
|
|
|
3072 }elsif($mates_file =~ /.bam$/){
|
|
|
3073 o$fh->open("$SAMTOOLS_BIN_DIR/samtools view $mates_file |") or die "$0: can't open ".$mates_file.":$!\n";#GALAXY
|
|
|
3074 }else{
|
|
|
3075 $fh->open("<".$mates_file) or die "$0: can't open ".$mates_file.":$!\n";
|
|
|
3076 }
|
|
|
3077
|
|
|
3078 while(<$fh>){
|
|
|
3079
|
|
|
3080 my @t=split;
|
|
|
3081 my $mate=$t[0];
|
|
|
3082
|
|
|
3083 my ($chr_read1, $chr_read2, $firstbase_read1, $firstbase_read2, $end_order_read1, $end_order_read2,);
|
|
|
3084
|
|
|
3085 next if(exists $pair{$mate});
|
|
|
3086
|
|
|
3087 next if (!readMateFile(\$chr_read1, \$chr_read2, \$firstbase_read1, \$firstbase_read2,\$end_order_read1, \$end_order_read2, \@t, $input_format,$tag_length));
|
|
|
3088
|
|
|
3089 next unless (exists $chrID->{$chr_read1} || exists $chrID->{$chr_read2});
|
|
|
3090 ($chr_read1, $chr_read2)= ($chrID->{$chr_read1},$chrID->{$chr_read2});
|
|
|
3091
|
|
|
3092 $pair{$mate}=undef;
|
|
|
3093 $record++;
|
|
|
3094
|
|
|
3095 my ($coord_start_read1,$coord_end_read1, $coord_start_read2,$coord_end_read2);
|
|
|
3096
|
|
|
3097 recupCoords($firstbase_read1,\$coord_start_read1,\$coord_end_read1,$tag_length->{$end_order_read1},$input_format);
|
|
|
3098 recupCoords($firstbase_read2,\$coord_start_read2,\$coord_end_read2,$tag_length->{$end_order_read2},$input_format);
|
|
|
3099
|
|
|
3100 my $length = abs($coord_start_read1-$coord_start_read2);
|
|
|
3101 $sumX += $length; #add to sum and sum^2 for mean and variance calculation
|
|
|
3102 $sumX2 += $length*$length;
|
|
|
3103
|
|
|
3104 for(my $i=1;$i<=$chr->{$chr_read1}->{'nb_frag'};$i++){
|
|
|
3105
|
|
|
3106 if (abs ($coord_start_read1-${$chr->{$chr_read1}->{$i}}[0]) <= $window_dist){
|
|
|
3107
|
|
|
3108 &addToDensity($density,$chr_read1,$i,\$nb_windows)
|
|
|
3109 if(overlap($coord_start_read1,$coord_end_read2,${$chr->{$chr_read1}->{$i}}[0],${$chr->{$chr_read1}->{$i}}[1]));
|
|
|
3110
|
|
|
3111 }else{
|
|
|
3112
|
|
|
3113 $i=getNextFrag($coord_start_read1,$i,${$chr->{$chr_read1}->{$i}}[0],$chr->{$chr_read1}->{nb_frag},$window_dist,$step);
|
|
|
3114 }
|
|
|
3115 }
|
|
|
3116
|
|
|
3117 if($record>=$warn){
|
|
|
3118 print LOG "-- $warn mate-pairs analysed - $nb_windows points created\n";
|
|
|
3119 $warn+=100000;
|
|
|
3120 }
|
|
|
3121 }
|
|
|
3122 $fh->close;
|
|
|
3123
|
|
|
3124 print LOG "-- $fchr : Total : $record mate-pairs analysed - $nb_windows points created\n";
|
|
|
3125
|
|
|
3126 if($record>0){
|
|
|
3127
|
|
|
3128 my $mu = $sumX/$record;
|
|
|
3129 my $sigma = sqrt($sumX2/$record - $mu*$mu);
|
|
|
3130 print LOG "-- $fchr : mu length = $mu, sigma length = $sigma\n";
|
|
|
3131 }
|
|
|
3132
|
|
|
3133 }
|
|
|
3134 #------------------------------------------------------------------------------#
|
|
|
3135 #------------------------------------------------------------------------------#
|
|
|
3136 sub ratio2segdup{
|
|
|
3137
|
|
|
3138 my($id,$density_file,$segdup_file)=@_;
|
|
|
3139
|
|
|
3140 print LOG "# Converting to circos format...\n";
|
|
|
3141
|
|
|
3142 open RATIO, "<$density_file" or die "$0: can't open $density_file :$!\n";
|
|
|
3143 open SEGDUP, ">$segdup_file" or die "$0: can't write in the output: $segdup_file :$!\n";
|
|
|
3144
|
|
|
3145 while(<RATIO>){
|
|
|
3146 chomp;
|
|
|
3147 my ($chr1,$start1,$end1,$ratio)=(split /\t/)[0,1,2,5];
|
|
|
3148 print SEGDUP "$id$chr1\t$start1\t$end1\t$ratio\n";
|
|
|
3149 }
|
|
|
3150
|
|
|
3151 close RATIO;
|
|
|
3152 close SEGDUP;
|
|
|
3153 print LOG "-- output created: $segdup_file\n";
|
|
|
3154 }
|
|
|
3155 #------------------------------------------------------------------------------#
|
|
|
3156 #------------------------------------------------------------------------------#
|
|
|
3157 sub ratio2bedfile{
|
|
|
3158
|
|
|
3159 my($density_file,$bed_file)=@_;
|
|
|
3160
|
|
|
3161 print LOG "# Converting to bedGraph format...\n";
|
|
|
3162
|
|
|
3163 open RATIO, "<$density_file" or die "$0: can't open $density_file :$!\n";
|
|
|
3164 open BED, ">$bed_file" or die "$0: can't write in the output: $bed_file :$!\n";
|
|
|
3165 print BED "track type=bedGraph name=\"$bed_file\" description=\"log ratios for cnv detection\" ".
|
|
|
3166 "visibility=2 color=255,0,0 alwaysZero=\"On\"\n";
|
|
|
3167
|
|
|
3168 while(<RATIO>){
|
|
|
3169 chomp;
|
|
|
3170 my ($chr1,$start1,$end1,$ratio)=(split /\t/)[0,1,2,5];
|
|
|
3171 $chr1 = "chr".$chr1 unless ($chr1 =~ m/chr/);
|
|
|
3172 print BED "$chr1\t".($start1-1)."\t$end1\t$ratio\n";
|
|
|
3173 }
|
|
|
3174
|
|
|
3175 close RATIO;
|
|
|
3176 close BED;
|
|
|
3177 print LOG "-- output created: $bed_file\n";
|
|
|
3178 }
|
|
|
3179 #------------------------------------------------------------------------------#
|
|
|
3180 #------------------------------------------------------------------------------#
|
|
|
3181 sub inverseSense{
|
|
|
3182
|
|
|
3183 my $mate_sense=$_[0];
|
|
|
3184 my %reverse=( 'F' => 'R' , 'R' => 'F' ,
|
|
|
3185 'FF' => 'RR', 'RR' => 'FF',
|
|
|
3186 'FR' => 'RF', 'RF' => 'FR');
|
|
|
3187 return $reverse{$mate_sense};
|
|
|
3188 }
|
|
|
3189
|
|
|
3190 #------------------------------------------------------------------------------#
|
|
|
3191 #------------------------------------------------------------------------------#
|
|
|
3192 sub getNextFrag{
|
|
|
3193
|
|
|
3194 my ($read_start,$frag_num,$frag_start,$frag_last,$window_dist,$step)=@_;
|
|
|
3195
|
|
|
3196 my $how_far = $read_start-$frag_start;
|
|
|
3197 my $nb_windows_toskip;
|
|
|
3198
|
|
|
3199 if($how_far>0){
|
|
|
3200
|
|
|
3201 $nb_windows_toskip=($how_far/$step)-($window_dist/$step);
|
|
|
3202 $nb_windows_toskip=~ s/\..*//;
|
|
|
3203 $nb_windows_toskip=0 if($nb_windows_toskip<0);
|
|
|
3204 return ($frag_num + $nb_windows_toskip);
|
|
|
3205 }
|
|
|
3206 return $frag_last;
|
|
|
3207 }
|
|
|
3208 #------------------------------------------------------------------------------#
|
|
|
3209 #------------------------------------------------------------------------------#
|
|
|
3210 sub getColor{
|
|
|
3211
|
|
|
3212 my($count,$hcolor,$format)=@_;
|
|
|
3213 for my $col ( keys % { $hcolor} ) {
|
|
|
3214 return $col if($count>=$hcolor->{$col}->[0] && $count<=$hcolor->{$col}->[1]);
|
|
|
3215 }
|
|
|
3216 return "white" if($format eq "circos");
|
|
|
3217 return "255,255,255" if($format eq "bed");
|
|
|
3218 }
|
|
|
3219 #------------------------------------------------------------------------------#
|
|
|
3220 #------------------------------------------------------------------------------#
|
|
|
3221 sub recupCoords{
|
|
|
3222
|
|
|
3223 my($c_hit,$cs_hit,$ce_hit,$tag_length,$input_format)=@_;
|
|
|
3224 my $strand = 'F';
|
|
|
3225
|
|
|
3226 if ($c_hit=~s/^\-//) {
|
|
|
3227 $strand='R';
|
|
|
3228 $$cs_hit=$c_hit;
|
|
|
3229 $$ce_hit=$c_hit-($tag_length-1);
|
|
|
3230 }else{
|
|
|
3231 $$cs_hit=$c_hit;
|
|
|
3232 $$ce_hit=$c_hit+($tag_length-1);
|
|
|
3233 }
|
|
|
3234 return $strand;
|
|
|
3235
|
|
|
3236 }
|
|
|
3237 #------------------------------------------------------------------------------#
|
|
|
3238 #------------------------------------------------------------------------------#
|
|
|
3239 sub overlap {
|
|
|
3240 my($cs_hit,$ce_hit,$cs_region,$ce_region)=@_;
|
|
|
3241 if( (($cs_hit < $cs_region) && ($ce_hit < $cs_region )) || (($cs_hit > $ce_region) && ($ce_hit > $ce_region )) ) {
|
|
|
3242 return 0;
|
|
|
3243 }
|
|
|
3244 return 1;
|
|
|
3245 }
|
|
|
3246 #------------------------------------------------------------------------------#
|
|
|
3247 #------------------------------------------------------------------------------#
|
|
|
3248 sub makeLink {
|
|
|
3249
|
|
|
3250 my ($link,$chr1,$frag1,$chr2,$frag2,$mt,$nb)=@_;
|
|
|
3251
|
|
|
3252 if($chr1>$chr2){
|
|
|
3253 ($chr1,$chr2)= ($chr2,$chr1);
|
|
|
3254 ($frag1,$frag2)= ($frag2,$frag1);
|
|
|
3255 }
|
|
|
3256
|
|
|
3257 if($chr1 == $chr2){
|
|
|
3258 if($frag1>$frag2){
|
|
|
3259 ($frag1,$frag2)= ($frag2,$frag1);
|
|
|
3260 }
|
|
|
3261 }
|
|
|
3262
|
|
|
3263 if(!exists $link->{$chr1}->{$chr2}->{$frag1}->{$frag2}){
|
|
|
3264 $link->{$chr1}->{$chr2}->{$frag1}->{$frag2}=$mt;
|
|
|
3265 $$nb++;
|
|
|
3266 }elsif($link->{$chr1}->{$chr2}->{$frag1}->{$frag2}!~/(^|,)$mt(,|$)/){
|
|
|
3267 $link->{$chr1}->{$chr2}->{$frag1}->{$frag2}.=",$mt";
|
|
|
3268 }
|
|
|
3269 }
|
|
|
3270 #------------------------------------------------------------------------------#
|
|
|
3271 #------------------------------------------------------------------------------#
|
|
|
3272 #fonction of adding the read to the density profile
|
|
|
3273 sub addToDensity {
|
|
|
3274
|
|
|
3275 my ($density,$chr1,$frag1,$nb)=@_;
|
|
|
3276
|
|
|
3277 if(!exists $density->{$chr1}->{$frag1}->{count}){
|
|
|
3278 $density->{$chr1}->{$frag1}->{count}=1;
|
|
|
3279 $$nb++;
|
|
|
3280 }else{
|
|
|
3281 $density->{$chr1}->{$frag1}->{count}++;
|
|
|
3282 }
|
|
|
3283 }
|
|
|
3284 #------------------------------------------------------------------------------#
|
|
|
3285 #------------------------------------------------------------------------------#
|
|
|
3286 sub floor {
|
|
|
3287 my $nb = $_[0];
|
|
|
3288 $nb=~ s/\..*//;
|
|
|
3289 return $nb;
|
|
|
3290 }
|
|
|
3291 #------------------------------------------------------------------------------#
|
|
|
3292 sub decimal{
|
|
|
3293
|
|
|
3294 my $num=shift;
|
|
|
3295 my $digs_to_cut=shift;
|
|
|
3296
|
|
|
3297 $num=sprintf("%.".($digs_to_cut-1)."f", $num) if ($num=~/\d+\.(\d){$digs_to_cut,}/);
|
|
|
3298
|
|
|
3299 return $num;
|
|
|
3300 }
|
|
|
3301
|
|
|
3302 #------------------------------------------------------------------------------#
|
|
|
3303 sub max {
|
|
|
3304
|
|
|
3305 my($max) = shift(@_);
|
|
|
3306 foreach my $temp (@_) {
|
|
|
3307 $max = $temp if $temp > $max;
|
|
|
3308 }
|
|
|
3309 return($max);
|
|
|
3310 }
|
|
|
3311 #------------------------------------------------------------------------------#
|
|
|
3312 #------------------------------------------------------------------------------#
|
|
|
3313 sub min {
|
|
|
3314
|
|
|
3315 my($min) = shift(@_);
|
|
|
3316 foreach my $temp (@_) {
|
|
|
3317 $min = $temp if $temp < $min;
|
|
|
3318 }
|
|
|
3319 return($min);
|
|
|
3320 }
|
|
|
3321 #------------------------------------------------------------------------------#
|
|
|
3322 #------------------------------------------------------------------------------#
|
|
|
3323 sub sortTablebyIndex{
|
|
|
3324 my ($tab1,$tab2)=@_;
|
|
|
3325 my @tab3;
|
|
|
3326
|
|
|
3327 foreach my $i (@$tab1){
|
|
|
3328 $tab3[$i]=$$tab2[$$tab1[$i]];
|
|
|
3329 }
|
|
|
3330 return @tab3;
|
|
|
3331 }
|
|
|
3332 #------------------------------------------------------------------------------#
|
|
|
3333 #------------------------------------------------------------------------------#
|
|
|
3334 sub round {
|
|
|
3335 my $number = shift || 0;
|
|
|
3336 my $dec = 10 ** (shift || 0);
|
|
|
3337 return int( $dec * $number + .5 * ($number <=> 0)) / $dec;
|
|
|
3338 }
|
|
|
3339 #------------------------------------------------------------------------------#
|
|
|
3340 #------------------------------------------------------------------------------#
|
|
|
3341 sub getUniqueTable{
|
|
|
3342
|
|
|
3343 my (@tab)=@_;
|
|
|
3344 my (%saw,@out)=();
|
|
|
3345 undef %saw;
|
|
|
3346 return sort(grep(!$saw{$_}++, @tab));
|
|
|
3347 }
|
|
|
3348 #------------------------------------------------------------------------------#
|
|
|
3349 #------------------------------------------------------------------------------#
|
|
|
3350 sub catFiles {
|
|
|
3351
|
|
|
3352 unlink("$_[1]") if(exists $_[1]);
|
|
|
3353 system qq( cat "$_" >> "$_[1]" ) for @{$_[0]};
|
|
|
3354 }
|
|
|
3355 #------------------------------------------------------------------------------#
|
|
|
3356 #------------------------------------------------------------------------------#
|
|
|
3357 #check if the configuration file is correct
|
|
|
3358 sub validateconfiguration{
|
|
|
3359
|
|
|
3360 my %conf=%{$_[0]};
|
|
|
3361 my $list_prgs="@ARGV";
|
|
|
3362
|
|
|
3363 my @general_params=qw(input_format mates_orientation read1_length read2_length mates_file cmap_file);
|
|
|
3364 my @detection_params=qw(split_mate_file window_size step_length split_mate_file);
|
|
|
3365 my @filtering_params=qw(split_link_file nb_pairs_threshold strand_filtering split_link_file);
|
|
|
3366 my @circos_params=qw(organism_id colorcode);
|
|
|
3367 my @bed_params=qw(colorcode);
|
|
|
3368 my @compare_params=qw(list_samples file_suffix);
|
|
|
3369
|
|
|
3370 foreach my $dir ($conf{general}{output_dir},$conf{general}{tmp_dir}){
|
|
|
3371
|
|
|
3372 unless (defined($dir)) {
|
|
|
3373 $dir = ".";
|
|
|
3374 }
|
|
|
3375 unless (-d $dir){
|
|
|
3376 mkdir $dir or die;
|
|
|
3377 }
|
|
|
3378 $dir.="/" if($dir!~/\/$/);
|
|
|
3379 }
|
|
|
3380
|
|
|
3381 unless (defined($conf{general}{num_threads})) {
|
|
|
3382 $conf{general}{num_threads} = 1;
|
|
|
3383 }
|
|
|
3384 $conf{general}{num_threads}=24 if($conf{general}{num_threads}>24);
|
|
|
3385
|
|
|
3386 if($list_prgs!~/links2compare/){
|
|
|
3387
|
|
|
3388 foreach my $p (@general_params){
|
|
|
3389 die("Error Config : The parameter \"$p\" is not defined\n") if (!defined $conf{general}{$p});
|
|
|
3390 }
|
|
|
3391
|
|
|
3392 $conf{general}{input_format}="sam" if($conf{general}{input_format} eq "bam");
|
|
|
3393
|
|
|
3394 unless (defined($conf{general}{sv_type})) {
|
|
|
3395 $conf{general}{sv_type} = "all";
|
|
|
3396 }
|
|
|
3397
|
|
|
3398 $conf{general}{read_lengths}={ 1=> $conf{general}{read1_length}, 2=> $conf{general}{read2_length}};
|
|
|
3399 }
|
|
|
3400
|
|
|
3401 if($list_prgs=~/(linking|cnv)/){
|
|
|
3402
|
|
|
3403 foreach my $p (@detection_params){
|
|
|
3404 die("Error Config : The parameter \"$p\" is not defined\n") if (!defined $conf{detection}{$p});
|
|
|
3405 }
|
|
|
3406
|
|
|
3407 die("Error Config : The parameter \"mates_file_ref\" is not defined\n") if($list_prgs=~/cnv/ && !defined $conf{detection}{mates_file_ref});
|
|
|
3408
|
|
|
3409 if($conf{detection}{step_length}>$conf{detection}{window_size}){
|
|
|
3410 die("Error Config : Parameter \"step_length\" should not exceed \"window size\"\n");
|
|
|
3411 }
|
|
|
3412
|
|
|
3413 unless (-d $conf{general}{tmp_dir}."/mates"){
|
|
|
3414 mkdir $conf{general}{tmp_dir}."/mates" or die;
|
|
|
3415 }
|
|
|
3416
|
|
|
3417 if($list_prgs=~/linking/){
|
|
|
3418 unless (-d $conf{general}{tmp_dir}."/links"){
|
|
|
3419 mkdir $conf{general}{tmp_dir}."/links" or die;
|
|
|
3420 }
|
|
|
3421 }
|
|
|
3422 if($list_prgs=~/cnv/){
|
|
|
3423 unless (-d $conf{general}{tmp_dir}."/density"){
|
|
|
3424 mkdir $conf{general}{tmp_dir}."/density" or die;
|
|
|
3425 }
|
|
|
3426 }
|
|
|
3427
|
|
|
3428 }
|
|
|
3429
|
|
|
3430 if($list_prgs=~/filtering/){
|
|
|
3431
|
|
|
3432 foreach my $p (@filtering_params) {
|
|
|
3433 die("Error Config : The filtering parameter \"$p\" is not defined\n") if (!defined $conf{filtering}{$p});
|
|
|
3434
|
|
|
3435 }
|
|
|
3436
|
|
|
3437 if(defined($conf{filtering}{chromosomes})) {
|
|
|
3438 my @chrs=split(",",$conf{filtering}{chromosomes});
|
|
|
3439 my $exclude=($chrs[0]=~/^\-/)? 1:0;
|
|
|
3440 for my $chrName (@chrs){
|
|
|
3441
|
|
|
3442 die("Error Config : The filtering parameter \"chromosomes\" is not valid\n")
|
|
|
3443 if(($chrName!~/^\-/ && $exclude) || ($chrName=~/^\-/ && !$exclude));
|
|
|
3444
|
|
|
3445 }
|
|
|
3446 }
|
|
|
3447
|
|
|
3448 if (( $conf{filtering}{order_filtering} )&& !$conf{filtering}{strand_filtering}) {
|
|
|
3449 die("Error Config : The parameter strand_filtering is set to \"0\" while order_filtering is selected".
|
|
|
3450 "\nChange strand_filtering to \"1\" if you want to use the order filtering\n");
|
|
|
3451 }
|
|
|
3452 if (( !defined($conf{filtering}{mu_length}) || !defined($conf{filtering}{sigma_length}) )&& $conf{filtering}{order_filtering}) {
|
|
|
3453 die("Error Config : You should set parameters \"mu_length\" and \"sigma_length\" to use order filtering\n");
|
|
|
3454 }
|
|
|
3455 if (( $conf{filtering}{insert_size_filtering} )&& !$conf{filtering}{strand_filtering}) {
|
|
|
3456 die("Error Config : The parameter strand_filtering is set to \"0\" while insert_size_filtering is selected".
|
|
|
3457 "\nChange strand_filtering to \"1\" if you want to use the insert size filtering\n");
|
|
|
3458 }
|
|
|
3459 if (( !defined($conf{filtering}{mu_length}) || !defined($conf{filtering}{sigma_length}) )&& $conf{filtering}{insert_size_filtering}) {
|
|
|
3460 die("Error Config : You should set parameters \"mu_length\" and \"sigma_length\" to use discriminate insertions from deletions\n");
|
|
|
3461 }
|
|
|
3462
|
|
|
3463 if (!defined($conf{filtering}{indel_sigma_threshold})) {
|
|
|
3464 $conf{filtering}{indel_sigma_threshold} = 2;
|
|
|
3465 }
|
|
|
3466 if (!defined($conf{filtering}{dup_sigma_threshold})) {
|
|
|
3467 $conf{filtering}{dup_sigma_threshold} = 2;
|
|
|
3468 }
|
|
|
3469 if (!defined($conf{filtering}{singleton_sigma_threshold})) {
|
|
|
3470 $conf{filtering}{singleton_sigma_threshold} = 4;
|
|
|
3471 }
|
|
|
3472
|
|
|
3473 if (!defined($conf{filtering}{nb_pairs_order_threshold})) {
|
|
|
3474 $conf{filtering}{nb_pairs_order_threshold} = 1;
|
|
|
3475 }
|
|
|
3476
|
|
|
3477 if (!defined($conf{filtering}{final_score_threshold})) {
|
|
|
3478 $conf{filtering}{final_score_threshold} = 0.8;
|
|
|
3479 }
|
|
|
3480
|
|
|
3481 if ($conf{filtering}{nb_pairs_order_threshold}>$conf{filtering}{nb_pairs_threshold}) {
|
|
|
3482 die("Error Config : Parameter \"nb_pairs_order_threshold\" should not exceed \"nb_pairs_threshold\"\n");
|
|
|
3483 }
|
|
|
3484
|
|
|
3485 }
|
|
|
3486
|
|
|
3487 if($list_prgs=~/2circos$/){
|
|
|
3488 foreach my $p (@circos_params) {
|
|
|
3489 next if($list_prgs=~/^ratio/ && $p eq "colorcode");
|
|
|
3490 die("Error Config : The circos parameter \"$p\" is not defined\n") if (!defined $conf{circos}{$p});
|
|
|
3491 }
|
|
|
3492 }
|
|
|
3493
|
|
|
3494 if($list_prgs=~/2bed$/){
|
|
|
3495 foreach my $p (@bed_params) {
|
|
|
3496 die("Error Config : The bed parameter \"$p\" is not defined\n") if (!defined $conf{bed}{$p});
|
|
|
3497 }
|
|
|
3498 }
|
|
|
3499
|
|
|
3500 if($list_prgs=~/links2compare/){
|
|
|
3501 foreach my $p (@compare_params) {
|
|
|
3502 die("Error Config : The compare parameter \"$p\" is not defined\n") if (!defined $conf{compare}{$p});
|
|
|
3503 }
|
|
|
3504
|
|
|
3505 unless (defined($conf{compare}{same_sv_type})) {
|
|
|
3506 $conf{compare}{same_sv_type} = 0;
|
|
|
3507 }
|
|
|
3508
|
|
|
3509 unless (defined($conf{compare}{min_overlap})) {
|
|
|
3510 $conf{compare}{min_overlap} = 1E-9;
|
|
|
3511 }
|
|
|
3512
|
|
|
3513 if($conf{compare}{circos_output}){
|
|
|
3514 foreach my $p (@circos_params) {
|
|
|
3515 next if($list_prgs=~/^ratio/ && $p eq "colorcode");
|
|
|
3516 die("Error Config : The circos parameter \"$p\" is not defined\n") if (!defined $conf{circos}{$p});
|
|
|
3517 }
|
|
|
3518 }
|
|
|
3519 if($conf{compare}{bed_output}){
|
|
|
3520 foreach my $p (@bed_params) {
|
|
|
3521 die("Error Config : The bed parameter \"$p\" is not defined\n") if (!defined $conf{bed}{$p});
|
|
|
3522 }
|
|
|
3523 die("Error Config : The compare parameter \"list_read_lengths\" is not defined\n") if (!defined $conf{compare}{list_read_lengths});
|
|
|
3524
|
|
|
3525 my @samples=split(",",$conf{compare}{list_samples});
|
|
|
3526 my @read_lengths=split(",",$conf{compare}{list_read_lengths});
|
|
|
3527 for my $i (0..$#samples){
|
|
|
3528 my @l=split("-",$read_lengths[$i]);
|
|
|
3529 $conf{compare}{read_lengths}{$samples[$i]}={ 1=> $l[0], 2=> $l[1]};
|
|
|
3530 }
|
|
|
3531 }
|
|
|
3532 }
|
|
|
3533
|
|
|
3534
|
|
|
3535 }
|
|
|
3536 #------------------------------------------------------------------------------#
|
|
|
3537 #::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::#
|
