-
Notifications
You must be signed in to change notification settings - Fork 4
/
Copy pathsupermap.pl
executable file
·1622 lines (1389 loc) · 70.6 KB
/
supermap.pl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#!/usr/bin/perl
# Supermap: Piecewise monotonic alignment map generator for Shuffle-LAGAN
# Author: Andrey Kislyuk ([email protected])
package Supermap;
require 5.005;
my ($VERSION) = ('$Id: supermap.pl,v 1.50 2005/06/15 22:40:04 kislyuk Exp $' =~ /,v\s+(\d+\S+)/o);
# Default constant values
my $overlap_factor = 0.8; # Aligns will be discarded if another align overlaps them by this factor or more in both seqs and has the same orientation
my $max_asym = 10; # Chains will be formed only if the resulting region's lengths differ by at most this factor
my $min_seq_score; # All aligns for sequences with this total score will be discarded. See getMinSeqScore
my $max_expand_len = 30000; # Aligns will be expanded or contracted on both sides on both strands by this amount up to the total length below
my $expand_factor = 4; # When one of an align's sequences is constrained in its expansion by a neighbor/start/end, the other one will be expanded by this times more than the first one
my $max_chainlen = 1500000; # Aligns will not be joined if the total length on either strand exceeds this. Set 0 to disable (no chain length limit)
my $max_job_size = 50000; # Maximum job size, in blat hits, for chunking when running glocal in parallel
my $erode_align = 15; # Amount by which to erode the coords of each align loaded (to avoid overlap problems when chaining)
my ($c1, $c2, $c3, $c4) = (100, 50, 400, 25); # BLAT->CHAOS score conversion parameters
#my $max_dist_y = 10000; # Join x-monotonic into same single-chain only if at most that apart in y-species.
my $default_lagan_dir = "/home/genome/glocal";
my $glocal_name = (0 ? "SLAGAN" : "glocal");
use Getopt::Long;
use File::Path;
use File::Copy;
use Cwd;
use IPC::Open2;
use IO::Handle;
#use Carp;
use strict;
use warnings;
no warnings "uninitialized";
sub main();
sub init();
sub getSeqSizes($$$);
sub prepareHits();
sub runSLAGAN();
sub reprintInputHits($$$);
sub processResults();
sub removeSLAGANOutput();
sub seqBelowMinScore($);
sub alignHashID($);
sub printChainToTemp($$$$);
sub chainBase1Hits($$);
sub chainBase2Hits($);
sub load2MHashes($);
sub loadBase2Hashes($);
sub postProcessRegions();
sub workerRun($$$$);
sub dequeueClustJobs($);
sub get_all_seqs($$);
sub isBLAT($);
sub useIf($$);
sub writeSizes($$);
sub getMinSeqScore($);
sub checkAlignCoords($);
sub expandSeq1($$);
sub expandSeq2($$);
sub finalExpand($$);
sub expSeq1Reg($$$$$);
sub expSeq2Reg($$$$$);
sub finalExpReg($$$$$);
# array index constants
use constant START1 => 0; use constant END1 => 1;
use constant START2 => 2; use constant END2 => 3;
use constant SEQ1 => 4; use constant SEQ2 => 5;
use constant ORIENT => 6; use constant ORIGIN => 7;
use constant SCORE => 8; use constant TOTSC => 9;
use constant HASHID => 10; use constant FLIPPED=> 11;
use constant CHALO1 => 12; use constant CHAHI1 => 13;
use constant CHALO2 => 14; use constant CHAHI2 => 15;
use constant CHALO1E=> 16; use constant CHAHI1E=> 17;
use constant CHALO2E=> 18; use constant CHAHI2E=> 19;
#use constant PREV1 => 8; use constant NEXT1 => 9;
#use constant PREV2 => 10; use constant NEXT2 => 11;
#use constant OSTART1=> 12; use constant OEND1 => 13;
#use constant OSTART2=> 14; use constant OEND2 => 15;
$SIG{'INT'} = $SIG{'QUIT'} = $SIG{'HUP'} = $SIG{'TRAP'} = $SIG{'ABRT'} = $SIG{'STOP'} = $SIG{'TERM'} = \&dequeueClustJobs;
my ($debug, $quiet, $outfile, $proflip, $skip, $no_pid, $input_glob, $input_dir,
$server, $db, $gen1, $gen2, $gen1sizefile, $gen2sizefile, $write_sizes1, $write_sizes2,
$score_file, $cfg, $cfg_file, $sizes1, $sizes2, $dbh, $tmp_dir, $tmp_prefix, $nodelete,
$clust_run_pid, $print_chains, $no_aligntotals, $no_clust_run, $num_jobs, $input_is_blat,
$force_overwrite, $print_csv, $using_GP, $slagan_params, $tmp_existed, $print_stats, $lagan_dir, $glocal_out_logfile);
my (@input_files);
my (%offsets1, %offsets2, %aligns1, %aligns2, %flipped_aligns);
my $supermapexec = $0; my $mycwd = getcwd(); $supermapexec =~ s/^\./$mycwd/ unless $supermapexec =~ /^\.\./; $supermapexec = $mycwd."/".$supermapexec if $supermapexec =~ /^\.\./;
die("$0: Problem resolving my name, \'$supermapexec\' is not a file") unless -f $supermapexec or $ARGV[0] eq "worker";
$0 = rindex($0, "/") > -1 ? substr($0, rindex($0, "/")+1) : $0;
$lagan_dir = $ENV{"LAGAN_DIR"} if defined $ENV{"LAGAN_DIR"};
$lagan_dir = $ENV{"LAGAN_DIR"} = $default_lagan_dir unless defined $ENV{"LAGAN_DIR"};
$lagan_dir =~ s/^\.\./$mycwd\/\.\./;
$lagan_dir =~ s/^\./$mycwd\//;
$ENV{"LAGAN_DIR"} = $lagan_dir;
print STDERR "$0: Warning: LAGAN_DIR=$lagan_dir is not a valid directory\n" unless -d $lagan_dir;
push @INC, $lagan_dir;
my $SLAGAN = $lagan_dir."/".$glocal_name;
my $error_file = "./$0.$$.error.log";
my $default_score_file = $lagan_dir."/test.score";
my $default_outfile = "$0.out";
my $worker_tmp_dir = "/tmp/$0.$$.worker/"; # The directory where workers store their intermediate files (two workers should not use the same directory)
my $usage = "
-infile=file \t Name of input file containing all hits for the two genomes
-outfile=file \t Output filename (default: $default_outfile)
-gen1=id \t First genome ID (must exist in the GPDB)
-gen2=id \t Second genome ID (must exist in the GPDB)
-sizes1=file \t File with sequence sizes for first genome
-sizes2=file \t File with sequence sizes for second genome
-bacteria \t Rearrange circular DNA to find a better alignment map
-server=hostname GPDB server (default: lemur)
-db=dbname \t GPDB name (default: GP)
-config=file \t GPDB config file (default: ~/.gprc)
-score=file \t Score file for SLAGAN (default: $default_score_file)
-glocal_out=file \t Save intermediate GLOCAL alignment hits to this file
-no_clust_run \t Run CPU/memory intensive jobs locally, not on the GP cluster
-tmp_dir=dir \t Working directory (default: /tmp/$0.pid)
-f \t\t Overwrite output file without prompting if it exists
-v \t\t Verbose mode
-q \t\t Quiet mode
-k \t\t Keep all temporary files
-expand_length=N Maximum length by which to expand alignments (default: $max_expand_len)
-max_length=N \t Maximum length for any alignment chain in either strand
\t\t (default: $max_chainlen)
-min_seq_score=N Sequences with total align score below this threshold will be
\t\t discarded (default: U penalty in SLAGAN score file)
-max_job_size=N Threshold, in hits, for splitting workload into separate jobs
\t\t for clust_run (default: $max_job_size)
-c1, c2, c3, c4=N: Score factors for BLAT->CHAOS conversion
\t\t (default: $c1, $c2, $c3, $c4)
Options may be abbreviated.
Input file format is BLAT or CHAOS. Sequence names should not contain spaces.
Alignments with negative scores are discarded.
Sequence size file format, one sequence per line: seq_name seq_size
";
exit(main());
# ___ Subroutines _______________
sub main() {
if ($ARGV[0] eq "worker") { workerRun($ARGV[1], $ARGV[2], $ARGV[3], $ARGV[4]); exit(0); } # Running SLAGAN in distributed mode
init();
print("$0: Retrieving sequence info...\n") unless $quiet;
$sizes1 = getSeqSizes($dbh, $gen1, $gen1sizefile);
(writeSizes($sizes1, $write_sizes1), exit(0)) if defined $write_sizes1;
$sizes2 = getSeqSizes($dbh, $gen2, $gen2sizefile);
(writeSizes($sizes2, $write_sizes2), exit(0)) if defined $write_sizes2;
die("$0: No sequence size data found. Stopped") if (keys(%$sizes1) < 1 or keys(%$sizes2) < 1);
die("$0: Flip mode is only applicable for two single-sequence organisms. Stopped") if ($proflip and not (keys(%$sizes1) == 1 and keys(%$sizes2) == 1));
# Sort and separate the alignments, run SLAGAN on them
prepareHits();
runSLAGAN();
# Chain SLAGAN alignments into supermonotonic chain and save the intermediate results
my ($dc, $sc1, $sc2) = processResults();
# Load the results back and expand regions, then print them
postProcessRegions();
print "$0: Output written to $outfile\n" unless $quiet;
print "$0: Intermediate files kept in $tmp_dir\n" if $nodelete and not $quiet;
rmdir $tmp_dir unless $tmp_existed or $nodelete;
return 0;
}
# Startup tasks
sub init() {
system('export LC_ALL="C"'); # Things may misbehave if locale is set to UTF-8
# Berkeley Genome Pipeline functionality is used if corresponding Perl modules are found in @INC
foreach my $dir (@INC) {
$using_GP = 1 if -f $dir."/GPDBI.pm" and -f $dir."/GPutils.pm";
}
useIf $using_GP, "GPDBI";
useIf $using_GP, "GPutils";
useIf 1, "Utils";
# useIf 1, "Desoverlap";
die("$0: GetOptions failed to retrieve options. Check the input options. Usage:".$usage) unless
GetOptions(
"server=s" => \$server,
"gen1=s" => \$gen1,
"gen2=s" => \$gen2,
"sizes1=s" => \$gen1sizefile,
"sizes2=s" => \$gen2sizefile,
"blatfile=s" => \$input_glob,
"infile=s" => \$input_glob,
"outfile=s" => \$outfile,
"glocal_out=s" => \$glocal_out_logfile,
"bacteria" => \$proflip,
"server=s" => \$server,
"db=s" => \$db,
"config=s" => \$cfg_file,
"tmp_dir=s" => \$tmp_dir,
"skip" => \$skip,
"no_pid" => \$no_pid,
"no_clust_run" => \$no_clust_run,
"print_chains" => \$print_chains,
"print_stats" => \$print_stats,
"no_aligntotals"=> \$no_aligntotals,
"print_csv" => \$print_csv,
"max_job_size" => \$max_job_size,
"max_length=i" => \$max_chainlen,
"expand_length=i"=>\$max_expand_len,
"min_seq_score=i"=>\$min_seq_score,
"max_asym=i" => \$max_asym,
"overlap_factor"=> \$overlap_factor,
"score=s" => \$score_file,
"c1=i" => \$c1,
"c2=i" => \$c2,
"c3=i" => \$c3,
"c4=i" => \$c4,
"slagan_params" => \$slagan_params,
"write_sizes1=s"=> \$write_sizes1,
"write_sizes2=s"=> \$write_sizes2,
"keep" => \$nodelete,
"f" => \$force_overwrite,
"v" => \$debug,
"q" => \$quiet
);
undef $quiet if $debug;
my @uinfo = getpwuid($>);
print("$0: Version ".$VERSION." started ".localtime()." by ".$uinfo[0]."\n") unless $quiet;
$tmp_prefix = $0.($no_pid ? "" : ".".$$);
unless ($no_clust_run) {
$no_clust_run = `which clust_run 2> /dev/null`; $no_clust_run = not $no_clust_run;
print("$0: clust_run not found - cluster operation disabled\n") if $no_clust_run and not $quiet;
}
if ($tmp_dir) {
$tmp_existed = 1 if -d $tmp_dir;
mkdir $tmp_dir unless -d $tmp_dir;
$tmp_dir .= "/" unless /\/^Z/;
} else {
$tmp_dir = "/tmp/".$tmp_prefix;
mkdir $tmp_dir;
$tmp_dir .= "/";
}
die("$0: No write permissions in working directory $tmp_dir. Stopped") unless -w $tmp_dir;
die("$0: Genome IDs or size files not specified. Usage:".$usage) unless ($gen1 or $gen1sizefile) and ($gen2 or $gen2sizefile);
die("$0: '-gen' options are invalid because GPDB is not available. Use '-sizes'. Stopped") if (($gen1 or $gen2) and not $using_GP);
die("$0: Sequence size file $gen1sizefile not found. Stopped") unless -f $gen1sizefile or $gen1;
die("$0: Sequence size file $gen2sizefile not found. Stopped") unless -f $gen2sizefile or $gen2;
die("$0: Maximum job size too small, must exceed 10000 hits. Stopped") if $max_job_size < 10000;
die("$0: Overlap factor must be between 0 and 1. Stopped") if $overlap_factor < 0 or $overlap_factor > 1;
print("$0: SLAGAN score file not specified, using default $default_score_file\n") unless $score_file or $quiet;
print("$0: Output file not specified, using default $default_outfile\n") unless $outfile or $quiet;
# Check input file or glob
if (defined $input_glob) {
if ($input_glob =~ /\//) { ($input_dir, $input_glob) = ($input_glob =~ /\A(.*\/)([^\/]+)\Z/); }
$input_glob .= "\$" unless $input_glob =~ /\$$/;
$input_glob = "^".$input_glob unless $input_glob =~ /^\^/;
@input_files = Utils::safe_glob($input_glob, $input_dir);
} elsif (@ARGV > 0) {
foreach my $file (@ARGV) {
if ($file =~ /\//) { ($input_dir, $file) = ($file =~ /\A(.*\/)([^\/]+)\Z/); }
push @input_files, $file;
}
} else { # TODO: split stdin for >2GB input
open(FH, "> $tmp_dir$tmp_prefix.in");
print FH while <STDIN>;
close FH;
push @input_files, "$tmp_prefix.in";
$input_dir = $tmp_dir;
}
unless ($input_dir =~ /\A\//) { $input_dir = $mycwd."/".$input_dir; }
die("$0: No input files matching \"$input_dir$input_glob\" found. Stopped") unless @input_files > 0;
print "$0: ".@input_files." input file(s)\n" if $debug;
# Check output file
$outfile = $default_outfile unless $outfile;
if (-f $outfile and not $force_overwrite and -t STDERR) {
print STDERR "$0: $outfile exists. Overwrite? (y/N, '-f' to force) ";
my $overwrite = <STDIN>; chomp $overwrite;
(print("Move \"$outfile\" or use option '-f'.\n"), exit(1)) unless ($overwrite eq "Y" or $overwrite eq "y" or $overwrite eq "yes");
}
open(FH, "> ".$outfile) or die("$0: Cannot open $outfile for writing: $!");
close FH;
# Check SLAGAN score file
$score_file = $default_score_file unless $score_file;
unless ($score_file =~ /\A\//) { $score_file = $mycwd."/".$score_file; }
$max_expand_len += $erode_align;
die("$0: max_length cannot be less than 0. Stopped") if $max_chainlen < 0;
$max_chainlen = 1000000000 if $max_chainlen == 0;
$max_chainlen -= 2*$max_expand_len;
# SLAGAN output for a given sequence will be discarded if the total score for the sequence is below this threshold. Default value is the SLAGAN unrelated gap penalty.
$min_seq_score = getMinSeqScore($score_file) unless defined $min_seq_score;
# Connect to GPDB
if ($using_GP) {
$GPutils::Error = "";
$cfg = read_gp_config(Get_Abs_Path($cfg_file)) or die($GPutils::Error);
$server ||= $cfg->Get_Val("DB", "server");
$db ||= $cfg->Get_Val("DB", "main_db");
$dbh = GPDBI->connect($server, 0, $db, undef, undef, "gp_cgi", undef, {PrintError => 0, RaiseError => 1});
}
}
# Load sequence names and sizes either from GPDB or from file
sub getSeqSizes($$$) {
my ($dbh, $dataset, $gen_size_file) = @_;
if ($dataset) {
return get_all_seqs($dbh, $dataset);
} else {
my %sizes;
open(FH, "< ".$gen_size_file) or die("$0: Could not open file $gen_size_file for reading: ".$!);
while (<FH>) {
chomp;
my ($seq, $size) = split;
die("$0: Invalid format in file $gen_size_file") unless $seq and $size;
$sizes{$seq} = $size;
}
close FH;
return \%sizes;
}
}
# Convert BLAT to CHAOS if necessary
# Flip hits on circular sequence if necessary
sub prepareHits() {
my ($cur_align);
local (*FH, *OUT1);
print "$0: Preparing files...\n" unless $quiet;
$input_is_blat = 1 if isBLAT($input_dir.$input_files[0]);
if ($input_is_blat) {
foreach my $file (@input_files) {
system('awk \'{$13=($13+$15)?$13:1; print $1,$2,$3";",$5,$6,$7"; '.
'score = "' . $c1 . '*$8-' . $c2 . '*$9-' . $c3 . '*($12+$14)-' . $c4 .
'*log($13+$15),"("$4")"}\''.
"< $input_dir$file > $tmp_dir$file.chaos");
}
} else {
foreach my $file (@input_files) {
system('ln -s "'.$input_dir.$file.'" "'.$tmp_dir.$file.'.chaos"');
}
}
if ($proflip) {
open(FH, "< ".$tmp_dir.$input_files[0].".chaos") or die("$0: Could not open file ".$tmp_dir.$input_files[0].".chaos for reading: ".$!);
open(OUT1, "> ".$tmp_dir.$input_files[0].".flipped.chaos") or die("$0: Could not open file ".$tmp_dir.$input_files[0].".flipped.chaos for writing: ".$!);
my (@seq1s, @seq1e, @seq2s, @seq2e, @scores, @orientations, @seqn1, @seqn2);
my ($seq1center, $seq2center, $seq1median, $seq2median);
my $i = 0;
while (<FH>) {
/\A[\s]*.*\s([\d]+)\s([\d]+)\;\s.*\s([\d]+)\s([\d]+)\;\sscore\s\=\s([e\d\.\+\-]+)\s\(([\+\-]+)\)/;
# ($seqn1[$i], $seq1s[$i], $seq1e[$i], $seqn2[$i], $seq2s[$i], $seq2e[$i], $scores[$i], $orientations[$i]) = ($1, $2, $3, $4, $5, $6, $7, $8);
($seq1s[$i], $seq1e[$i], $seq2s[$i], $seq2e[$i], $scores[$i], $orientations[$i]) = ($1, $2, $3, $4, $5, $6);
if ($seq1s[$i] > $seq1e[$i]) { my $j = $seq1s[$i]; $seq1s[$i] = $seq1e[$i]; $seq1e[$i] = $j; }
if ($seq2s[$i] > $seq2e[$i]) { my $j = $seq2s[$i]; $seq2s[$i] = $seq2e[$i]; $seq2e[$i] = $j; }
$i++;
}
# For each interval pair,
# if the seq1 interval median is greater than seq1 median, and the corresponding interval median in seq2 is less than seq2 median,
# OR if the seq1 interval median is less than seq1 median, and the corresponding interval median in seq2 is greater than seq2 median,
# set start of interval in seq1 to 2CoM1 - previous end of interval
# set end of interval in seq1 to 2CoM1 - previous start of interval
# flip the orientation (+/-)
$seq1center = $$sizes1{(keys(%$sizes1))[0]} / 2;
$seq2center = $$sizes2{(keys(%$sizes2))[0]} / 2;
my $flip_counter = 0;
foreach $i (0..@seq1s-1) {
$seq1median = ($seq1s[$i] + $seq1e[$i]) / 2;
$seq2median = ($seq2s[$i] + $seq2e[$i]) / 2;
if (($seq1median > $seq1center and $seq2median < $seq2center)
or ($seq1median < $seq1center and $seq2median > $seq2center)) {
my $j = $seq2s[$i];
$seq2s[$i] = (2 * $seq2center) - $seq2e[$i];
$seq2e[$i] = (2 * $seq2center) - $j;
if ($orientations[$i] eq "+") { $orientations[$i] = "-"; } else { $orientations[$i] = "+"; }
$cur_align = [];
$$cur_align[START1] = $seq1s[$i]; $$cur_align[START2] = $seq2s[$i];
$$cur_align[END1] = $seq1e[$i]; $$cur_align[END2] = $seq2e[$i];
$$cur_align[SCORE] = $scores[$i]; $$cur_align[ORIENT] = $orientations[$i];
$$cur_align[SEQ1] = (keys(%$sizes1))[0]; $$cur_align[SEQ2] = (keys(%$sizes2))[0];
$$cur_align[START1] += $erode_align; $$cur_align[END1] -= $erode_align;
$$cur_align[START2] += $erode_align; $$cur_align[END2] -= $erode_align;
$flipped_aligns{alignHashID($cur_align)} = $cur_align;
$flip_counter++;
}
print OUT1 "seq1 ".$seq1s[$i]." ".$seq1e[$i]."; seq2 ".$seq2s[$i]." ".$seq2e[$i]."; score = ".$scores[$i]." (".$orientations[$i].")\n";
}
close FH; close OUT1;
print "$0: Single-sequence flip mode: ".($flip_counter+0)." hits flipped\n" if $debug;
}
}
# Load all hits into a hash table, then write the hits for each sequence into a file
# Run SLAGAN on each of these files, via worker instances either on the cluster or sequentially
sub runSLAGAN() {
my ($clust_run_invoke, $num_jobs, $sort_pid1, $sort_pid2, $sort_pid3, $one_seq_mode,
$cur_align, $next_align, $curlen1, $curlen2, $nextlen1, $nextlen2, $overlap1, $overlap2, $dump_count);
local (*RH1, *WH1, *RH2, *WH2, *RH3, *WH3, *IN, *DUPES);
# my $filter = Desoverlap->new($overlap_factor, $debug);
print "$0: Sorting input hits...\n" if $debug;
open(DUPES, "> supermap.duplicates") if $debug;
$one_seq_mode = 1 if (keys(%$sizes1) == 1 and keys(%$sizes2) == 1);
$sort_pid1 = open2(\*RH1, \*WH1, "sort -k 1,1 -k 2,2n"); # pre-scan
$sort_pid2 = open2(\*RH2, \*WH2, "sort -k 1,1 -k 2,2n"); # gen1base
$sort_pid3 = open2(\*RH3, \*WH3, "sort -k 4,4 -k 5,5n"); # gen2base
# Sort input on seq1
foreach my $file (@input_files) {
open(IN, "< $tmp_dir$file".($proflip?".flipped":"").".chaos");
print WH1 while <IN>;
close IN;
}
close WH1;
# Scan input, check if start2, end2 are ascending for sorting, erode alignments
while (<RH1>) {
/\A[\s]*(.*)\s([\d]+)\s([\d]+)\;\s(.*)\s([\d]+)\s([\d]+)\;\sscore\s\=\s([e\d\.\+\-]+)\s\(([\+\-]+)\)/o;
$next_align=[];
($$next_align[SEQ1], $$next_align[START1], $$next_align[END1], $$next_align[SEQ2], $$next_align[START2], $$next_align[END2], $$next_align[SCORE], $$next_align[ORIENT])
= ($1, $2, $3, $4, $5, $6, $7, $8);
next if $$next_align[SCORE] <= 0;
if ($one_seq_mode) { $$next_align[SEQ1] = (keys(%$sizes1))[0]; $$next_align[SEQ2] = (keys(%$sizes2))[0]; }
checkAlignCoords($next_align);
unless ($$next_align[END1]-$$next_align[START1] <= $erode_align*2 or $$next_align[END2]-$$next_align[START2] <= $erode_align*2) {
$$next_align[START1] += $erode_align; $$next_align[END1] -= $erode_align;
$$next_align[START2] += $erode_align; $$next_align[END2] -= $erode_align;
}
=head1
# Overlap scan
if ($$next_align[START1] <= $$cur_align[END1] and $$next_align[END1] >= $$cur_align[START1] # overlap in seq1
and $$next_align[START2] <= $$cur_align[END2] and $$next_align[END2] >= $$cur_align[START2] # overlap in seq2
and $$cur_align[SEQ1] eq $$next_align[SEQ1] and $$cur_align[SEQ2] eq $$next_align[SEQ2]
and $$cur_align[ORIENT] eq $$next_align[ORIENT]) {
($curlen1, $curlen2, $nextlen1, $nextlen2)
= ($$cur_align[END1] - $$cur_align[START1] + 1, $$cur_align[END2] - $$cur_align[START2] + 1,
$$next_align[END1] - $$next_align[START1] + 1, $$next_align[END2] - $$next_align[START2] + 1);
if ($$next_align[START1] <= $$cur_align[START1] and $$next_align[END1] >= $$cur_align[END1]) {
$overlap1 = $$cur_align[END1] - $$cur_align[START1] + 1; # next covers cur
} elsif ($$next_align[START1] <= $$cur_align[START1]) {
$overlap1 = $$next_align[END1] - $$cur_align[START1] + 1; # next is to the left
} elsif ($$next_align[END1] >= $$cur_align[END1]) {
$overlap1 = $$cur_align[END1] - $$next_align[START1] + 1; # next is to the right
} else {
$overlap1 = $$next_align[END1] - $$next_align[START1] + 1; # cur covers next
}
if ($$next_align[START2] <= $$cur_align[START2] and $$next_align[END2] >= $$cur_align[END2]) {
$overlap2 = $$cur_align[END2] - $$cur_align[START2] + 1;
} elsif ($$next_align[START2] <= $$cur_align[START2]) {
$overlap2 = $$next_align[END2] - $$cur_align[START2] + 1;
} elsif ($$next_align[END2] >= $$cur_align[END2]) {
$overlap2 = $$cur_align[END2] - $$next_align[START2] + 1;
} else {
$overlap2 = $$next_align[END2] - $$next_align[START2] + 1;
}
die("$0: Bad internal state") if $overlap1 < 0 or $overlap2 < 0;
if (($overlap1 / $curlen1 > $overlap_factor) and ($overlap2 / $curlen2 > $overlap_factor)
and $$cur_align[SCORE] <= $$next_align[SCORE]) {
$dump_count++;
print DUPES "Cur: (".$$cur_align[START1]."-".$$cur_align[END1].")(".$$cur_align[START2]."-".$$cur_align[END2].") ".$$cur_align[SCORE]." over with (".$$next_align[START1]."-".$$next_align[END1].")(".$$next_align[START2]."-".$$next_align[END2].") ".$$next_align[SCORE]."\n" if $debug;
$cur_align = $next_align; next; # discard current align
} elsif (($overlap1 / $nextlen1 > $overlap_factor) and ($overlap2 / $nextlen2 > $overlap_factor)
and $$cur_align[SCORE] >= $$next_align[SCORE]) {
$dump_count++;
print DUPES "Nxt: (".$$next_align[START1]."-".$$next_align[END1].")(".$$next_align[START2]."-".$$next_align[END2].") ".$$next_align[SCORE]." over with (".$$cur_align[START1]."-".$$cur_align[END1].")(".$$cur_align[START2]."-".$$cur_align[END2].") ".$$cur_align[SCORE]."\n" if $debug;
next; # discard next align
}
}
=cut
foreach my $cur_align ($next_align){ # (@{$filter->put($next_align)}) {
print WH2 $$cur_align[SEQ1]." ".$$cur_align[START1]." ".$$cur_align[END1]."; ".$$cur_align[SEQ2]." ".$$cur_align[START2]." ".$$cur_align[END2]."; "."score = ".$$cur_align[SCORE]." (".$$cur_align[ORIENT].")\n";
print WH3 $$cur_align[SEQ1]." ".$$cur_align[START1]." ".$$cur_align[END1]."; ".$$cur_align[SEQ2]." ".$$cur_align[START2]." ".$$cur_align[END2]."; "."score = ".$$cur_align[SCORE]." (".$$cur_align[ORIENT].")\n";
}
# print WH2 $$cur_align[SEQ1]." ".$$cur_align[START1]." ".$$cur_align[END1]."; ".$$cur_align[SEQ2]." ".$$cur_align[START2]." ".$$cur_align[END2]."; "."score = ".$$cur_align[SCORE]." (".$$cur_align[ORIENT].")\n" if @$cur_align;
# print WH3 $$cur_align[SEQ1]." ".$$cur_align[START1]." ".$$cur_align[END1]."; ".$$cur_align[SEQ2]." ".$$cur_align[START2]." ".$$cur_align[END2]."; "."score = ".$$cur_align[SCORE]." (".$$cur_align[ORIENT].")\n" if @$cur_align;
# $cur_align = $next_align;
}
# $filter->printAll();
# Flush alignments remaining in filter buffer
# foreach my $cur_align (@{$filter->getBuffer()}) {
# print WH2 $$cur_align[SEQ1]." ".$$cur_align[START1]." ".$$cur_align[END1]."; ".$$cur_align[SEQ2]." ".$$cur_align[START2]." ".$$cur_align[END2]."; "."score = ".$$cur_align[SCORE]." (".$$cur_align[ORIENT].")\n" if $cur_align != 0;
# print WH3 $$cur_align[SEQ1]." ".$$cur_align[START1]." ".$$cur_align[END1]."; ".$$cur_align[SEQ2]." ".$$cur_align[START2]." ".$$cur_align[END2]."; "."score = ".$$cur_align[SCORE]." (".$$cur_align[ORIENT].")\n" if $cur_align != 0;
# }
close RH1; waitpid $sort_pid1, 0;
close WH2;
$num_jobs = reprintInputHits(1, 1, \*RH2);
close RH2; waitpid $sort_pid2, 0;
close WH3;
$num_jobs = reprintInputHits(2, $num_jobs, \*RH3);
close RH3; waitpid $sort_pid3, 0;
close DUPES if defined fileno DUPES;
# print STDERR "$0: Warning: ".$filter->{dump_count}." near duplicate alignments discarded (overlap factor $overlap_factor)\n" if $filter->{dump_count} and not $quiet;
open(FH, "> ".$tmp_dir."CLUSTER_JOB_PARAMS") or die;
foreach my $i (1..$num_jobs-1) {
print FH "worker JOB".$i.".tar ".$score_file." ".$SLAGAN.($debug ? " -v" : "");
print FH " << JOB$i.tar > CLUSTER_JOB_MESSAGES.$i >> CLUSTER_JOB_ERRMSG.$i" unless $no_clust_run;
print FH "\n";
}
close FH;
if ($no_clust_run) {
open(FH, "< ".$tmp_dir."CLUSTER_JOB_PARAMS") or die;
print "$0: Running ".($num_jobs-1)." SLAGAN jobs locally...\n" unless $quiet;
while (<FH>) {
chomp;
print("Job $.: \"$0 $_\"\n") if $debug;
system("cd $tmp_dir; $supermapexec ".$_);
}
close FH;
} else {
$clust_run_invoke = "clust_run -program=".$supermapexec." -parameters=".$tmp_dir."CLUSTER_JOB_PARAMS -init_dir=$tmp_dir -wait";
print "$0: Running ".($num_jobs-1)." distributed SLAGAN jobs with clust_run...\n" unless $quiet;
print "$0: \"$clust_run_invoke\"\n" if $debug;
if ($clust_run_pid = fork()) { # I am the parent
waitpid($clust_run_pid, 0);
} elsif (not defined $clust_run_pid) {
die("$0: Could not fork");
} else { # I am the child
die("$0: Could not exec \"$clust_run_invoke\"") unless exec($clust_run_invoke);
}
undef $clust_run_pid;
}
foreach my $i (1..$num_jobs-1) {
system("cd $tmp_dir; tar -xf ".$tmp_dir."JOB".$i.".results.tar");
unlink $tmp_dir."JOB".$i.".tar" unless $nodelete;
unlink $tmp_dir."JOB".$i.".results.tar" unless $nodelete;
unlink $tmp_dir."CLUSTER_JOB_MESSAGES.$i" unless $nodelete;
unlink $tmp_dir."CLUSTER_JOB_ERRMSG.$i" unless $nodelete;
}
unlink "$tmp_dir$input_glob.chaos" unless $nodelete;
unlink $tmp_dir."CLUSTER_JOB_PARAMS" unless $nodelete;
foreach my $file (@input_files) {
unlink $tmp_dir.$file.".chaos" unless $nodelete;
}
}
sub reprintInputHit($$$) {
my ($base_gen, $align, $FH) = @_;
if ($base_gen == 1 and $$align[ORIENT] eq "+") {
print $FH $$align[SEQ1]." ".$$align[START1]." ".$$align[END1]."; ".$$align[SEQ2]." ".$$align[START2]." ".$$align[END2]."; "."score = ".$$align[SCORE]." (".$$align[ORIENT].")\n";
} elsif ($base_gen == 1 and $$align[ORIENT] eq "-") {
print $FH $$align[SEQ1]." ".$$align[START1]." ".$$align[END1]."; ".$$align[SEQ2]." ".$$align[END2]." ".$$align[START2]."; "."score = ".$$align[SCORE]." (".$$align[ORIENT].")\n";
} elsif ($base_gen == 2 and $$align[ORIENT] eq "+") {
print $FH $$align[SEQ2]." ".$$align[START2]." ".$$align[END2]."; ".$$align[SEQ1]." ".$$align[START1]." ".$$align[END1]."; "."score = ".$$align[SCORE]." (".$$align[ORIENT].")\n";
} elsif ($base_gen == 2 and $$align[ORIENT] eq "-") {
print $FH $$align[SEQ2]." ".$$align[START2]." ".$$align[END2]."; ".$$align[SEQ1]." ".$$align[END1]." ".$$align[START1]."; "."score = ".$$align[SCORE]." (".$$align[ORIENT].")\n";
} else {
die("$0: Bad internal state from hit ".$$align[SEQ1]." ".$$align[START1]." ".$$align[END1]."; ".$$align[SEQ2]." ".$$align[START2]." ".$$align[END2]."; "."score = ".$$align[SCORE]." (".$$align[ORIENT].")");
}
}
sub writeJobFile($$) {
my ($job_id, $seq_list) = @_;
local *LIST;
open(LIST, "| cd $tmp_dir; xargs tar --append --file=".$tmp_dir."JOB".$job_id.".tar");
foreach my $file (sort alnum keys(%$seq_list)) { $file =~ /\/([^\/]+)$/; print LIST $1." "; }
close LIST;
foreach my $file (sort alnum keys(%$seq_list)) { unlink $file unless $nodelete; }
}
# Separate input into files based on sequence name and reverse order in gen2base hits
sub reprintInputHits($$$) {
my ($base_gen, $job_id, $RH) = @_;
my ($one_seq_mode, $line_count, $prev_seq, $cur_seq, $cur_align);
my (%cur_seq_list, %pruned_sizes);
local (*OUT, *LIST);
$one_seq_mode = 1 if (keys(%$sizes1) == 1 and keys(%$sizes2) == 1);
print "$0: Reprinting hits (base genome $base_gen)..." if $debug;
$line_count = 0;
while (<$RH>) {
/\A[\s]*(.*)\s([\d]+)\s([\d]+)\;\s(.*)\s([\d]+)\s([\d]+)\;\sscore\s\=\s([e\d\.\+\-]+)\s\(([\+\-]+)\)/o;
$cur_align=[];
($$cur_align[SEQ1], $$cur_align[START1], $$cur_align[END1], $$cur_align[SEQ2], $$cur_align[START2], $$cur_align[END2], $$cur_align[SCORE], $$cur_align[ORIENT])
= ($1, $2, $3, $4, $5, $6, $7, $8);
$cur_seq = ($base_gen == 1 ? $$cur_align[SEQ1] : $$cur_align[SEQ2]);
if ($cur_seq ne $prev_seq) {
$pruned_sizes{$cur_seq} = ($base_gen == 1 ? $$sizes1{$cur_seq} : $$sizes2{$cur_seq});
print " ".$cur_seq if $debug;
close OUT if defined fileno OUT;
open(OUT, "> ".$tmp_dir.$input_files[0].".gen".$base_gen."base.".$cur_seq.".chaos") or die("$0: Could not open file ".$tmp_dir.$input_files[0].".gen".$base_gen."base.".$cur_seq.".chaos for writing: ".$!);
if ($line_count > $max_job_size) {
writeJobFile($job_id, \%cur_seq_list);
undef %cur_seq_list; $line_count = 0; $job_id++;
}
$cur_seq_list{$tmp_dir.$input_files[0].".gen".$base_gen."base.".$cur_seq.".chaos"} = 1;
}
reprintInputHit($base_gen, $cur_align, \*OUT) if @$cur_align;
$prev_seq = $cur_seq;
# $cur_align = $next_align;
$line_count++;
}
# reprintInputHit($base_gen, $next_align, \*OUT) if @$next_align;
writeJobFile($job_id, \%cur_seq_list);
$job_id++;
close OUT;
print "\n" if $debug;
$sizes1 = \%pruned_sizes if $base_gen == 1;
$sizes2 = \%pruned_sizes if $base_gen == 2;
return $job_id;
}
sub seqBelowMinScore($) {
my ($line) = @_;
$line =~ /\A[\s]*\([\d]+\s[\d]+\)\=\([\d]+\s[\d]+\)\s([\d\.\-]+)\s[\+\-]+\s\[([\d\.\-]+)\][\s]*s1\:.*[\s]*s2\:.*\n\Z/;
die("$0: Unable to extract score values from SLAGAN output:\n$line") if not defined $2;
return ($2 < $min_seq_score);
}
sub processResults() {
my ($cur_seq, $input_prefix, $dropped_seqs, $sort_pid, $sort_pid2);
local (*RH, *WH, *IN, *OUT, *hashesDM_RH, *hashesDM_WH);
print "$0: Loading SLAGAN output...\n" unless $quiet;
open(GLOCAL_OUT_LOG, "> ".$glocal_out_logfile) if $glocal_out_logfile;
# Sort gen2base aligns on seq1, then seq2, then start2, then print them to separate files, one file per gen1 seq
# These files will be loaded on demand when scanning gen1base aligns (chainBase1Hits())
$sort_pid = open2(\*RH, \*WH, "sort -k 9,9 -k 7,7 -k 1.2,1n"); # input is base 2, key is 9 because a space is expected between s2: and seq2name
$input_prefix = $tmp_dir.$input_files[0].".gen2base";
foreach my $seq (sort alnum keys(%$sizes2)) {
open(IN, "< $input_prefix.$seq.chaos.glocal-out") or (delete($$sizes2{$seq}), next);
my $line = <IN>;
die("$0: Empty SLAGAN output file $input_prefix.$seq.chaos.glocal-out, check corresponding job logs. Stopped") unless $line;
if (seqBelowMinScore($line)) { print "$0: Discarding file $input_prefix.$seq.chaos.glocal-out - score too low ($1<$min_seq_score)\n" if $debug; next; }
seek IN, 0, 0; # back to start
print WH while <IN>;
close IN;
}
close WH or die("$0: Error executing sort");
while (<RH>) {
/\ss2\:[\s]*([^\s]+)[\s]*\n\Z/;
if ($1 ne $cur_seq or not defined $cur_seq) {
next unless $1;
close OUT if defined fileno OUT;
$cur_seq = $1;
open(OUT, "> $input_prefix.sorted-gen1.$cur_seq.chaos.glocal-out") or die("$0: Could not open file $input_prefix.sorted-gen1.$cur_seq.chaos.glocal-out for writing: ".$!);
}
print OUT $_;
}
close RH; close OUT if defined fileno OUT;
waitpid $sort_pid, 0;
# Sort gen1base aligns on seq1, then start1
$sort_pid = open2(\*RH, \*WH, "sort -k 7,7 -k 1.2,1n"); # input is base 1
$input_prefix = $tmp_dir.$input_files[0].".gen1base";
foreach my $seq (sort alnum keys(%$sizes1)) {
open(IN, "< $input_prefix.$seq.chaos.glocal-out") or (delete($$sizes1{$seq}), next);
my $line = <IN>;
if (seqBelowMinScore($line)) { $dropped_seqs++; print "$0: Discarding file $input_prefix.$seq.chaos.glocal-out - score too low ($1<$min_seq_score)\n" if $debug; next; }
seek IN, 0, 0; # back to start
print WH while <IN>;
if ($glocal_out_logfile) { seek IN, 0, 0; print GLOCAL_OUT_LOG while <IN>; }
close IN;
unlink "$input_prefix.$seq.chaos.glocal-out" unless $nodelete;
}
unlink $input_prefix.".chaos" unless $nodelete;
close WH or die("$0: Error executing sort");
# Feed the gen1base aligns to the 2M/1M1 chain scanner (chainBase1Hits())
# The hashesDM handle is used to write 2M aligns' hashes to be sorted in seq2 order
print "$0: Generating supermonotonic map...\n" unless $quiet;
$sort_pid2 = open2(\*hashesDM_RH, \*hashesDM_WH, "sort -k 2,2");
chainBase1Hits(*RH, *hashesDM_WH);
close RH;
waitpid $sort_pid, 0;
close hashesDM_WH or die("$0: Error executing sort");
# Print sorted 2M aligns' hashes, one file per gen2 seq
undef $cur_seq;
while(<hashesDM_RH>) {
my $line = $_;
$line =~ /\A[^\s]+\s([^\s]+)\s[^\s]+\n\Z/;
if ($1 ne $cur_seq or not defined $cur_seq) {
close OUT if defined fileno OUT;
$cur_seq = $1;
open(OUT, "> $tmp_dir".$input_files[0].".hashesDM.gen2.$cur_seq") or die("$0: Could not open file $tmp_dir".$input_files[0].".hashesDM.gen2.$cur_seq for writing: ".$!);
}
print OUT $line;
}
close hashesDM_RH;
waitpid $sort_pid2, 0;
# Sort gen2base aligns on seq2, then start2
$sort_pid = open2(\*RH, \*WH, "sort -k 7,7 -k 1.2,1n"); # input is base 2
$input_prefix = $tmp_dir.$input_files[0].".gen2base";
foreach my $seq (sort alnum keys(%$sizes2)) {
open(IN, "< $input_prefix.$seq.chaos.glocal-out") or next;
my $line = <IN>;
if (seqBelowMinScore($line)) { $dropped_seqs++; print "$0: Discarding file $input_prefix.$seq.chaos.glocal-out - score too low ($1<$min_seq_score)\n" if $debug; next; }
seek IN, 0, 0; # back to start
print WH while <IN>;
close IN;
unlink "$input_prefix.$seq.chaos.glocal-out" unless $nodelete;
}
unlink $input_prefix.".chaos" unless $nodelete;
close WH or die("$0: Error executing sort");
# Feed the gen2base aligns to the 1M2 chain scanner (chainBase2Hits())
chainBase2Hits(*RH);
close RH;
waitpid $sort_pid, 0;
close GLOCAL_OUT_LOG if defined fileno GLOCAL_OUT_LOG;
removeSLAGANOutput();
print STDERR "$0: Warning: Alignments for $dropped_seqs sequences discarded due to total score below cutoff ($min_seq_score)\n" if $dropped_seqs and not $quiet;
}
sub removeSLAGANOutput() {
my $input_prefix = $tmp_dir.$input_files[0].".gen1base";
foreach my $seq (sort alnum keys(%$sizes1)) { unlink "$input_prefix.$seq.chaos.glocal-out" unless $nodelete; }
unlink $input_prefix.".chaos" unless $nodelete;
$input_prefix = $tmp_dir.$input_files[0].".gen2base";
foreach my $seq (sort alnum keys(%$sizes2)) { unlink "$input_prefix.$seq.chaos.glocal-out" unless $nodelete; }
unlink $input_prefix.".chaos" unless $nodelete;
rmdir $tmp_dir;
}
sub alignHashID($) {
my ($align) = @_;
# return 23*$$align[START1] + 41*$$align[START2] + 61*$$align[END1] + 83*$$align[END2];
return $$align[SEQ1].":".$$align[START1]."-".$$align[END1]."=".$$align[SEQ2].":".$$align[START2]."-".$$align[END2];
}
# The chain writer lags the chainer by two chains because the full contents of neighboring chains must be known.
sub printChainToTemp($$$$) {
my ($FH, $prev_chain, $cur_chain, $next_chain) = @_;
return unless defined $cur_chain;
my $type = ${$$cur_chain[0]}[ORIGIN];
my ($first_align, $last_align) = ($$cur_chain[0], $$cur_chain[@$cur_chain-1]);
print $FH ${$$cur_chain[0]}[ORIGIN]." ".@$cur_chain." ".
$$first_align[START1]." ".$$first_align[END1]." ".$$first_align[START2]." ".$$first_align[END2]." ".
$$first_align[SEQ1]." ".$$first_align[SEQ2]." ".$$first_align[ORIENT]." ".$$first_align[SCORE]." ".
$$last_align[START1]." ".$$last_align[END1]." ".$$last_align[START2]." ".$$last_align[END2]." ".
$$last_align[SEQ1]." ".$$last_align[SEQ2]." ".$$last_align[ORIENT]." ".$$last_align[SCORE];
if ($print_chains) {
foreach my $align (@$cur_chain) {
print $FH " ".$$align[START1]." ".$$align[END1]." ".$$align[START2]." ".$$align[END2];
}
}
print $FH "\n";
}
sub chainBase1Hits($$) {
my ($FH, $hashesDM) = @_;
local *OUT;
my ($cur_align, $prev_align, $cur_chain, $prev_chain, $pre_prev_chain, $chain_start_2M, $chain_start_1M1,
$cur_seq, $align_peers, $flip_counter);
my @bad_aligns; my %base2peers;
while (<$FH>) {
/\A[\s]*\(([\d]+)\s([\d]+)\)\=\(([\d]+)\s([\d]+)\)\s([\d\.\-]+)\s([\+\-]+)\s\[([\d\.\-]+)\][\s]*s1\:(.*)[\s]*s2\:(.*)\n\Z/;
next if ($1==$2); # skip null alignments
(push(@bad_aligns, $_), next) unless $1 and $2 and $3 and $4 and $5 and $6;
$cur_align = [];
($$cur_align[START1], $$cur_align[END1], $$cur_align[START2], $$cur_align[END2], $$cur_align[SCORE], $$cur_align[ORIENT], $$cur_align[TOTSC], $$cur_align[SEQ1], $$cur_align[SEQ2])
= ($1, $2, $3, $4, $5, $6, $7, $8, $9);
$$cur_align[SEQ1] =~ s/^\s+//; $$cur_align[SEQ1] =~ s/\s+$//;
$$cur_align[SEQ2] =~ s/^\s+//; $$cur_align[SEQ2] =~ s/\s+$//;
#warn("Seen: ".$_) if $$cur_align[SEQ1] eq "AC002301.1";
checkAlignCoords($cur_align);
if ($proflip and defined $flipped_aligns{alignHashID($cur_align)}) {
my $seq2center = $$sizes2{(keys(%$sizes2))[0]} / 2;
my $j = $$cur_align[START2];
$$cur_align[START2] = (2 * $seq2center) - $$cur_align[END2];
$$cur_align[END2] = (2 * $seq2center) - $j;
if ($$cur_align[ORIENT] eq "+") { $$cur_align[ORIENT] = "-"; } else { $$cur_align[ORIENT] = "+"; }
$$cur_align[FLIPPED]=1;
$flip_counter++;
}
$$cur_align[HASHID] = alignHashID($cur_align);
if ($$cur_align[SEQ1] ne $cur_seq) {
#warn("Handling seq trans") if $prev_align and $$prev_align[SEQ1] eq "AC002301.1";
printChainToTemp(*OUT, $pre_prev_chain, $prev_chain, $cur_chain);# unless defined $cur_seq;
printChainToTemp(*OUT, $prev_chain, $cur_chain, undef);# unless defined $cur_seq;
undef $chain_start_2M; undef $chain_start_1M1; undef $prev_align;
undef $pre_prev_chain; undef $prev_chain; undef $cur_chain;
$cur_seq = $$cur_align[SEQ1];
%base2peers = %{loadBase2Hashes($tmp_dir.$input_files[0].".gen2base.sorted-gen1.$cur_seq.chaos.glocal-out")};
close OUT if defined fileno OUT;
open(OUT, "> ".$tmp_dir.$input_files[0].".2MM1.$cur_seq");
}
$align_peers = $base2peers{$$cur_align[HASHID]};
$$cur_align[ORIGIN] = defined($align_peers) ? 2 : 1;
if ($chain_start_2M and defined $align_peers and defined $prev_align # continue open 2M chain
and (($$cur_align[ORIENT] eq "+" and $$cur_align[START2] > $$prev_align[END2]
and $$prev_align[HASHID] eq $$align_peers[0])
or ($$cur_align[ORIENT] eq "-" and $$cur_align[END2] < $$prev_align[START2]
and $$prev_align[HASHID] eq $$align_peers[1])
or ($$cur_align[FLIPPED] and ($$cur_align[ORIENT] eq "+" and $$cur_align[START2] < $$prev_align[END2]
and $$prev_align[HASHID] eq $$align_peers[0])
or ($$cur_align[ORIENT] eq "-" and $$cur_align[END2] > $$prev_align[START2]
and $$prev_align[HASHID] eq $$align_peers[1])))
and $$cur_align[ORIENT] eq $$prev_align[ORIENT]
and $$cur_align[FLIPPED] eq $$prev_align[FLIPPED]
and $$cur_align[SEQ2] eq $$prev_align[SEQ2]
and ($$cur_align[START1] > $$prev_align[END1] or ($$cur_align[FLIPPED] and $$cur_align[START1] > $$prev_align[END1]))
and abs($$cur_align[END1] - $$chain_start_2M[START1]) < $max_chainlen
and abs($$cur_align[END2] - $$chain_start_2M[START2]) < $max_chainlen
#and abs($$cur_align[END1] - $$chain_start_2M[START1])/abs($$cur_align[END2] - $$chain_start_2M[START2]) < $max_asym
#and abs($$cur_align[END2] - $$chain_start_2M[START2])/abs($$cur_align[END1] - $$chain_start_2M[START1]) < $max_asym
) {
push(@$cur_chain, $cur_align);
print $hashesDM $$cur_align[SEQ1]."\t".$$cur_align[SEQ2]."\t".$$cur_align[HASHID]."\n";
} elsif (defined $align_peers) { # start new 2M chain
printChainToTemp(*OUT, $pre_prev_chain, $prev_chain, $cur_chain);
$chain_start_2M = $cur_align; undef $chain_start_1M1;
$pre_prev_chain = $prev_chain; $prev_chain = $cur_chain;
$cur_chain = [$cur_align];
print $hashesDM $$cur_align[SEQ1]."\t".$$cur_align[SEQ2]."\t".$$cur_align[HASHID]."\n";
} elsif ($chain_start_1M1 and defined $prev_align # continue open 1M1 chain
and ((($$cur_align[ORIENT] eq "+" and $$cur_align[START2] > $$prev_align[END2])
or ($$cur_align[ORIENT] eq "-" and $$cur_align[END2] < $$prev_align[START2]))
or ($$cur_align[FLIPPED] and (($$cur_align[ORIENT] eq "+" and $$cur_align[START2] < $$prev_align[END2])
or ($$cur_align[ORIENT] eq "-" and $$cur_align[END2] > $$prev_align[START2]))))
and $$cur_align[ORIENT] eq $$prev_align[ORIENT]
and $$cur_align[FLIPPED] eq $$prev_align[FLIPPED]
and $$cur_align[SEQ2] eq $$prev_align[SEQ2]
and ($$cur_align[START1] > $$prev_align[END1] or ($$cur_align[FLIPPED] and $$cur_align[START1] > $$prev_align[END1]))
and abs($$cur_align[END1] - $$chain_start_1M1[START1]) < $max_chainlen
and abs($$cur_align[END2] - $$chain_start_1M1[START2]) < $max_chainlen
#and abs($$cur_align[END1] - $$chain_start_1M1[START1])/abs($$cur_align[END2] - $$chain_start_1M1[START2]) < $max_asym
#and abs($$cur_align[END2] - $$chain_start_1M1[START2])/abs($$cur_align[END1] - $$chain_start_1M1[START1]) < $max_asym
) {
push(@$cur_chain, $cur_align);
} else { # start new 1M1 chain
printChainToTemp(*OUT, $pre_prev_chain, $prev_chain, $cur_chain);
$chain_start_1M1 = $cur_align; undef $chain_start_2M;
$pre_prev_chain = $prev_chain; $prev_chain = $cur_chain;
$cur_chain = [$cur_align];
}
$prev_align = $cur_align;
}
printChainToTemp(*OUT, $pre_prev_chain, $prev_chain, $cur_chain);
printChainToTemp(*OUT, $prev_chain, $cur_chain, undef);
print "$0: Single-sequence flip mode: ".($flip_counter+0)." gen1base hits backflipped\n" if $debug and $proflip;
warn "$0: Warning: ".@bad_aligns." bad SLAGAN alignments discarded" if @bad_aligns > 0;
}
# Input is base 2, i.e. (start2 end2)=(start1 end1)...
sub chainBase2Hits($) {
my ($FH) = @_;
local *OUT;
my ($cur_align, $prev_align, $cur_chain, $prev_chain, $pre_prev_chain, $chain_start_2M, $chain_start_1M2,
$cur_seq, $align_is_2M, $flip_counter);
my @bad_aligns; my %aligns2M;
while(<$FH>) {
/\A[\s]*\(([\d]+)\s([\d]+)\)\=\(([\d]+)\s([\d]+)\)\s([\d\.\-]+)\s([\+\-]+)\s\[([\d\.\-]+)\][\s]*s1\:(.*)[\s]*s2\:(.*)\n\Z/;
next if ($1==$2); # skip null alignments
(push(@bad_aligns, $_), next) unless $1 and $2 and $3 and $4 and $5 and $6;
$cur_align = [];
($$cur_align[START2], $$cur_align[END2], $$cur_align[START1], $$cur_align[END1], $$cur_align[SCORE], $$cur_align[ORIENT], $$cur_align[TOTSC], $$cur_align[SEQ2], $$cur_align[SEQ1])
= ($1, $2, $3, $4, $5, $6, $7, $8, $9);
$$cur_align[SEQ1] =~ s/^\s+//; $$cur_align[SEQ1] =~ s/\s+$//;
$$cur_align[SEQ2] =~ s/^\s+//; $$cur_align[SEQ2] =~ s/\s+$//;
checkAlignCoords($cur_align);
if ($proflip and defined $flipped_aligns{alignHashID($cur_align)}) {
my $seq2center = $$sizes2{(keys(%$sizes2))[0]} / 2;
my $j = $$cur_align[START2];
$$cur_align[START2] = (2 * $seq2center) - $$cur_align[END2];
$$cur_align[END2] = (2 * $seq2center) - $j;
if ($$cur_align[ORIENT] eq "+") { $$cur_align[ORIENT] = "-"; } else { $$cur_align[ORIENT] = "+"; }
$$cur_align[FLIPPED] = 1;
$flip_counter++;
}
$$cur_align[HASHID] = alignHashID($cur_align);
if ($$cur_align[SEQ2] ne $cur_seq) {
printChainToTemp(*OUT, $pre_prev_chain, $prev_chain, $cur_chain) if $$prev_chain[0][ORIGIN] == 3;# and not defined $cur_seq;
printChainToTemp(*OUT, $prev_chain, $cur_chain, undef) if $$cur_chain[0][ORIGIN] == 3;# and not defined $cur_seq;
undef $chain_start_1M2; undef $prev_align;
undef $pre_prev_chain; undef $prev_chain; undef $cur_chain;
$cur_seq = $$cur_align[SEQ2];
%aligns2M = %{load2MHashes($tmp_dir.$input_files[0].".hashesDM.gen2.$cur_seq")};
close OUT if defined fileno OUT;
open(OUT, "> ".$tmp_dir.$input_files[0].".M2.$cur_seq");
}
$$cur_align[ORIGIN] = defined($aligns2M{$$cur_align[HASHID]}) ? 2 : 3;
if (defined $aligns2M{$$cur_align[HASHID]}) { # align is 2M
my $prev_ch_last_al = $prev_chain ? $$prev_chain[scalar(@$prev_chain)-1] : [];
printChainToTemp(*OUT, $pre_prev_chain, $prev_chain, $cur_chain) if $$prev_chain[0][ORIGIN] == 3;
undef $chain_start_1M2; # close 1M2 chain
$chain_start_2M = $cur_align;
$pre_prev_chain = $prev_chain; $prev_chain = $cur_chain;
$cur_chain = [$cur_align];
} elsif ($chain_start_1M2 # continue open 1M2 chain
and ((($$cur_align[ORIENT] eq "+" and $$cur_align[START1] > $$prev_align[END1])
or ($$cur_align[ORIENT] eq "-" and $$cur_align[END1] < $$prev_align[START1]))
or ($$cur_align[FLIPPED] and (($$cur_align[ORIENT] eq "+" and $$cur_align[START1] < $$prev_align[END1])
or ($$cur_align[ORIENT] eq "-" and $$cur_align[END1] > $$prev_align[START1]))))
and $$cur_align[ORIENT] eq $$prev_align[ORIENT]
and $$cur_align[SEQ1] eq $$prev_align[SEQ1]
and $$cur_align[FLIPPED] == $$prev_align[FLIPPED]
and ($$cur_align[START2] > $$prev_align[END2] or ($$cur_align[FLIPPED] and $$cur_align[START2] < $$prev_align[END2]))
and abs($$cur_align[END1] - $$chain_start_1M2[START1]) < $max_chainlen
and abs($$cur_align[END2] - $$chain_start_1M2[START2]) < $max_chainlen
#and abs($$cur_align[END1] - $$chain_start_1M2[START1])/abs($$cur_align[END2] - $$chain_start_1M2[START2]) < $max_asym
#and abs($$cur_align[END2] - $$chain_start_1M2[START2])/abs($$cur_align[END1] - $$chain_start_1M2[START1]) < $max_asym
) {
push(@$cur_chain, $cur_align);
} else { # start new 1M2 chain
my $prev_ch_last_al = $prev_chain ? $$prev_chain[scalar(@$prev_chain)-1] : [];
printChainToTemp(*OUT, $pre_prev_chain, $prev_chain, $cur_chain) if $$prev_chain[0][ORIGIN] == 3;
$chain_start_1M2 = $cur_align;
$pre_prev_chain = $prev_chain; $prev_chain = $cur_chain;
$cur_chain = [$cur_align];
}
$prev_align = $cur_align;
}
my $prev_ch_last_al = $prev_chain ? $$prev_chain[scalar(@$prev_chain)-1] : [];
printChainToTemp(*OUT, $pre_prev_chain, $prev_chain, $cur_chain) if $$prev_chain[0][ORIGIN] == 3;
printChainToTemp(*OUT, $prev_chain, $cur_chain, undef) if $$cur_chain[0][ORIGIN] == 3;
print "$0: Single-sequence flip mode: ".($flip_counter+0)." gen2base hits backflipped\n" if $debug and $proflip;
warn "$0: Warning: ".@bad_aligns." bad SLAGAN alignments discarded" if @bad_aligns > 0;
}
# Input: file with lines of the form "seq1 seq2 hash" (seq2 should be the same per file)
# Output: hash(key->align hash ID, value->1). Input file is deleted.
sub load2MHashes($) {
my ($file) = @_;
my %hashes;
local *FH;
open(FH, "< $file") or return {};
while (<FH>) {
/\A[^\s]+\t[^\s]+\t([^\s]+)\n\Z/;
warn("Hash collision in \"$_\" vs. \"".$hashes{$1}."\"") if defined $hashes{$1};