1#!@PERL@
2
3################################################################################
4#   Programmer: Adam M Phillippy, The Institute for Genomic Research
5#         File: mummerplot
6#         Date: 01 / 08 / 03
7#               01 / 06 / 05 rewritten (v3.0)
8#
9#        Usage:
10#    mummerplot  [options]  <match file>
11#
12#                Try 'mummerplot -h' for more information.
13#
14#      Purpose: To generate a gnuplot plot for the display of mummer, nucmer,
15#               promer, and show-tiling alignments.
16#
17################################################################################
18
19use lib "@LIB_DIR@";
20use Foundation;
21use strict;
22use IO::Socket;
23
24my $BIN_DIR     = "@BIN_DIR@";
25my $LIB_DIR     = "@LIB_DIR@";
26my $GNUPLOT_EXE = "@GNUPLOT_EXE@";
27
28
29#================================================================= Globals ====#
30#-- terminal types
31my $X11    = "x11";
32my $PS     = "postscript";
33my $PNG    = "png";
34
35#-- terminal sizes
36my $SMALL  = "small";
37my $MEDIUM = "medium";
38my $LARGE  = "large";
39
40my %TERMSIZE =
41    (
42     $X11 => { $SMALL => 500, $MEDIUM => 700,  $LARGE => 900  }, # screen pix
43     $PS  => { $SMALL => 1,   $MEDIUM => 2,    $LARGE => 3    }, # pages
44     $PNG => { $SMALL => 800, $MEDIUM => 1024, $LARGE => 1400 }  # image pix
45     );
46
47#-- terminal format
48my $FFACE    = "Courier";
49my $FSIZE    = "8";
50my $TFORMAT  = "%.0f";
51my $MFORMAT  = "[%.0f, %.0f]";
52
53#-- output suffixes
54my $FILTER  = "filter";
55my $FWDPLOT = "fplot";
56my $REVPLOT = "rplot";
57my $HLTPLOT = "hplot";
58my $GNUPLOT = "gnuplot";
59
60my %SUFFIX =
61    (
62     $FILTER  => ".filter",
63     $FWDPLOT => ".fplot",
64     $REVPLOT => ".rplot",
65     $HLTPLOT => ".hplot",
66     $GNUPLOT => ".gp",
67     $PS      => ".ps",
68     $PNG     => ".png"
69     );
70
71
72#================================================================= Options ====#
73my $OPT_breaklen;                  # -b option
74my $OPT_color;                     # --[no]color option
75my $OPT_coverage;                  # --[no]coverage option
76my $OPT_filter;                    # -f option
77my $OPT_layout;                    # -l option
78my $OPT_prefix    = "out";         # -p option
79my $OPT_rv;                        # --rv option
80my $OPT_terminal  = $X11;          # -t option
81my $OPT_IdR;                       # -r option
82my $OPT_IdQ;                       # -q option
83my $OPT_IDRfile;                   # -R option
84my $OPT_IDQfile;                   # -Q option
85my $OPT_rport;                     # -rport option
86my $OPT_qport;                     # -qport option
87my $OPT_size      = $SMALL;        # -small, -medium, -large
88my $OPT_SNP;                       # -S option
89my $OPT_xrange;                    # -x option
90my $OPT_yrange;                    # -y option
91my $OPT_title;                     # -title option
92
93my $OPT_Mfile;                     # match file
94my $OPT_Dfile;                     # delta filter file
95my $OPT_Ffile;                     # .fplot output
96my $OPT_Rfile;                     # .rplot output
97my $OPT_Hfile;                     # .hplot output
98my $OPT_Gfile;                     # .gp output
99my $OPT_Pfile;                     # .ps .png output
100
101my $OPT_gpstatus;                  # gnuplot status
102
103my $OPT_ONLY_USE_FATTEST;          # Only use fattest alignment for layout
104
105
106#============================================================== Foundation ====#
107my $VERSION = '3.5';
108
109my $USAGE = qq~
110  USAGE: mummerplot  [options]  <match file>
111    ~;
112
113my $HELP = qq~
114  USAGE: mummerplot  [options]  <match file>
115
116  DESCRIPTION:
117    mummerplot generates plots of alignment data produced by mummer, nucmer,
118    promer or show-tiling by using the GNU gnuplot utility. After generating
119    the appropriate scripts and datafiles, mummerplot will attempt to run
120    gnuplot to generate the plot. If this attempt fails, a warning will be
121    output and the resulting .gp and .[frh]plot files will remain so that the
122    user may run gnuplot independently. If the attempt succeeds, either an x11
123    window will be spawned or an additional output file will be generated
124    (.ps or .png depending on the selected terminal). Feel free to edit the
125    resulting gnuplot script (.gp) and rerun gnuplot to change line thinkness,
126    labels, colors, plot size etc.
127
128  MANDATORY:
129    match file      Set the alignment input to 'match file'
130                    Valid inputs are from mummer, nucmer, promer and
131                    show-tiling (.out, .cluster, .delta and .tiling)
132
133  OPTIONS:
134    -b|breaklen     Highlight alignments with breakpoints further than
135                    breaklen nucleotides from the nearest sequence end
136    --[no]color     Color plot lines with a percent similarity gradient or
137                    turn off all plot color (default color by match dir)
138                    If the plot is very sparse, edit the .gp script to plot
139                    with 'linespoints' instead of 'lines'
140    -c
141    --[no]coverage  Generate a reference coverage plot (default for .tiling)
142    --depend        Print the dependency information and exit
143    -f
144    --filter        Only display .delta alignments which represent the "best"
145                    hit to any particular spot on either sequence, i.e. a
146                    one-to-one mapping of reference and query subsequences
147    -h
148    --help          Display help information and exit
149    -l
150    --layout        Layout a .delta multiplot in an intelligible fashion,
151                    this option requires the -R -Q options
152    --fat           Layout sequences using fattest alignment only
153    -p|prefix       Set the prefix of the output files (default '$OPT_prefix')
154    -rv             Reverse video for x11 plots
155    -r|IdR          Plot a particular reference sequence ID on the X-axis
156    -q|IdQ          Plot a particular query sequence ID on the Y-axis
157    -R|Rfile        Plot an ordered set of reference sequences from Rfile
158    -Q|Qfile        Plot an ordered set of query sequences from Qfile
159                    Rfile/Qfile Can either be the original DNA multi-FastA
160                    files or lists of sequence IDs, lens and dirs [ /+/-]
161    -r|rport        Specify the port to send reference ID and position on
162                    mouse double click in X11 plot window
163    -q|qport        Specify the port to send query IDs and position on mouse
164                    double click in X11 plot window
165    -s|size         Set the output size to small, medium or large
166                    --small --medium --large (default '$OPT_size')
167    -S
168    --SNP           Highlight SNP locations in each alignment
169    -t|terminal     Set the output terminal to x11, postscript or png
170                    --x11 --postscript --png (default '$OPT_terminal')
171    -t|title        Specify the gnuplot plot title (default none)
172    -x|xrange       Set the xrange for the plot '[min:max]'
173    -y|yrange       Set the yrange for the plot '[min:max]'
174    -V
175    --version       Display the version information and exit
176    ~;
177
178my @DEPEND =
179    (
180     "$LIB_DIR/Foundation.pm",
181     "$BIN_DIR/delta-filter",
182     "$BIN_DIR/show-coords",
183     "$BIN_DIR/show-snps",
184     "gnuplot"
185     );
186
187my $tigr = new TIGR::Foundation
188    or die "ERROR: TIGR::Foundation could not be initialized\n";
189
190$tigr -> setVersionInfo ($VERSION);
191$tigr -> setUsageInfo ($USAGE);
192$tigr -> setHelpInfo ($HELP);
193$tigr -> addDependInfo (@DEPEND);
194
195
196#=========================================================== Function Decs ====#
197sub GetParseFunc( );
198
199sub ParseIDs($$);
200
201sub ParseDelta($);
202sub ParseCluster($);
203sub ParseMummer($);
204sub ParseTiling($);
205
206sub LayoutIDs($$);
207sub SpanXwY ($$$$$);
208
209sub PlotData($$$);
210sub WriteGP($$);
211sub RunGP( );
212sub ListenGP($$);
213
214sub ParseOptions( );
215
216
217#=========================================================== Function Defs ====#
218MAIN:
219{
220    my @aligns;                # (sR eR sQ eQ sim lenR lenQ idR idQ)
221    my %refs;                  # (id => (off, len, [1/-1]))
222    my %qrys;                  # (id => (off, len, [1/-1]))
223
224    #-- Get the command line options (sets OPT_ global vars)
225    ParseOptions( );
226
227
228    #-- Get the alignment type
229    my $parsefunc = GetParseFunc( );
230
231    if ( $parsefunc != \&ParseDelta &&
232         ($OPT_filter || $OPT_layout || $OPT_SNP) ) {
233        print STDERR "WARNING: -f -l -S only work with delta input\n";
234        undef $OPT_filter;
235        undef $OPT_layout;
236        undef $OPT_SNP;
237    }
238
239    #-- Parse the reference and query IDs
240    if    ( defined $OPT_IdR ) { $refs{$OPT_IdR} = [ 0, 0, 1 ]; }
241    elsif ( defined $OPT_IDRfile ) {
242        ParseIDs ($OPT_IDRfile, \%refs);
243    }
244
245    if    ( defined $OPT_IdQ ) { $qrys{$OPT_IdQ} = [ 0, 0, 1 ]; }
246    elsif ( defined $OPT_IDQfile ) {
247        ParseIDs ($OPT_IDQfile, \%qrys);
248    }
249
250
251    #-- Filter the alignments
252    if ( $OPT_filter || $OPT_layout ) {
253        print STDERR "Writing filtered delta file $OPT_Dfile\n";
254        system ("$BIN_DIR/delta-filter -r -q $OPT_Mfile > $OPT_Dfile")
255            and die "ERROR: Could not run delta-filter, $!\n";
256        if ( $OPT_filter ) { $OPT_Mfile = $OPT_Dfile; }
257    }
258
259
260    #-- Parse the alignment data
261    $parsefunc->(\@aligns);
262
263
264    #-- Layout the alignment data if requested
265    if ( $OPT_layout ) {
266        if ( scalar (keys %refs) || scalar (keys %qrys) ) {
267            LayoutIDs (\%refs, \%qrys);
268        }
269        else {
270            print STDERR "WARNING: --layout option only works with -R or -Q\n";
271            undef $OPT_layout;
272        }
273    }
274
275
276    #-- Plot the alignment data
277    PlotData (\@aligns, \%refs, \%qrys);
278
279
280    #-- Write the gnuplot script
281    WriteGP (\%refs, \%qrys);
282
283
284    #-- Run gnuplot script and fork a clipboard listener
285    unless ( $OPT_gpstatus == -1 ) {
286
287        my $child = 1;
288        if ( $OPT_gpstatus == 0 && $OPT_terminal eq $X11 ) {
289            print STDERR "Forking mouse listener\n";
290            $child = fork;
291        }
292
293        #-- parent runs gnuplot
294        if ( $child ) {
295            RunGP( );
296            kill 1, $child;
297        }
298        #-- child listens to clipboard
299        elsif ( defined $child ) {
300            ListenGP(\%refs, \%qrys);
301        }
302        else {
303            print STDERR "WARNING: Could not fork mouse listener\n";
304        }
305    }
306
307    exit (0);
308}
309
310
311#------------------------------------------------------------ GetParseFunc ----#
312sub GetParseFunc ( )
313{
314    my $fref;
315
316    open (MFILE, "<$OPT_Mfile")
317        or die "ERROR: Could not open $OPT_Mfile, $!\n";
318
319    $_ = <MFILE>;
320    if ( !defined ) { die "ERROR: Could not read $OPT_Mfile, File is empty\n" }
321
322  SWITCH: {
323      #-- tiling
324      if ( /^>\S+ \d+ bases/ ) {
325          $fref = \&ParseTiling;
326          last SWITCH;
327      }
328
329      #-- mummer
330      if ( /^> \S+/ ) {
331          $fref = \&ParseMummer;
332          last SWITCH;
333      }
334
335      #-- nucmer/promer
336      if ( /^(\S+) (\S+)/ ) {
337          if ( ! defined $OPT_IDRfile ) {
338              $OPT_IDRfile = $1;
339          }
340          if ( ! defined $OPT_IDQfile ) {
341              $OPT_IDQfile = $2;
342          }
343
344          $_ = <MFILE>;
345          if ( (defined)  &&  (/^NUCMER$/  ||  /^PROMER$/) ) {
346              $_ = <MFILE>;   # sequence header
347              $_ = <MFILE>;   # alignment header
348              if ( !defined ) {
349                  $fref = \&ParseDelta;
350                  last SWITCH;
351              }
352              elsif ( /^\d+ \d+ \d+ \d+ \d+ \d+ \d+$/ ) {
353                  $fref = \&ParseDelta;
354                  last SWITCH;
355              }
356              elsif ( /^[ \-][1-3] [ \-][1-3]$/ ) {
357                  $fref = \&ParseCluster;
358                  last SWITCH;
359              }
360          }
361      }
362
363      #-- default
364      die "ERROR: Could not read $OPT_Mfile, Unrecognized file type\n";
365  }
366
367    close (MFILE)
368        or print STDERR "WARNING: Trouble closing $OPT_Mfile, $!\n";
369
370    return $fref;
371}
372
373
374#---------------------------------------------------------------- ParseIDs ----#
375sub ParseIDs ($$)
376{
377    my $file = shift;
378    my $href = shift;
379
380    open (IDFILE, "<$file")
381        or print STDERR "WARNING: Could not open $file, $!\n";
382
383    my $dir;
384    my $aref;
385    my $isfasta;
386    my $offset = 0;
387    while ( <IDFILE> ) {
388        #-- Ignore blank lines
389        if ( /^\s*$/ ) { next; }
390
391        #-- FastA header
392        if ( /^>(\S+)/ ) {
393            if ( exists $href->{$1} ) {
394                print STDERR "WARNING: Duplicate sequence '$1' ignored\n";
395                undef $aref;
396                next;
397            }
398
399            if ( !$isfasta ) { $isfasta = 1; }
400            if ( defined $aref ) { $offset += $aref->[1] - 1; }
401
402            $aref = [ $offset, 0, 1 ];
403            $href->{$1} = $aref;
404            next;
405        }
406
407        #-- FastA sequence
408        if ( $isfasta  &&  /^\S+$/ ) {
409            if ( defined $aref ) { $aref->[1] += (length) - 1; }
410            next;
411        }
412
413        #-- ID len dir
414        if ( !$isfasta  &&  /^(\S+)\s+(\d+)\s+([+-]?)$/ ) {
415            if ( exists $href->{$1} ) {
416                print STDERR "WARNING: Duplicate sequence '$1' ignored\n";
417                undef $aref;
418                next;
419            }
420
421            $dir = (defined $3 && $3 eq "-") ? -1 : 1;
422            $aref = [ $offset, $2, $dir ];
423            $offset += $2 - 1;
424            $href->{$1} = $aref;
425            next;
426        }
427
428        #-- default
429        print STDERR "WARNING: Could not parse $file\n$_";
430        undef %$href;
431        last;
432    }
433
434    close (IDFILE)
435        or print STDERR "WARNING: Trouble closing $file, $!\n";
436}
437
438
439#-------------------------------------------------------------- ParseDelta ----#
440sub ParseDelta ($)
441{
442    my $aref = shift;
443
444    print STDERR "Reading delta file $OPT_Mfile\n";
445
446    open (MFILE, "<$OPT_Mfile")
447        or die "ERROR: Could not open $OPT_Mfile, $!\n";
448
449    my @align;
450    my $ispromer;
451    my ($sim, $tot);
452    my ($lenR, $lenQ, $idR, $idQ);
453
454    $_ = <MFILE>;
455    $_ = <MFILE>;
456    $ispromer = /^PROMER/;
457
458    while ( <MFILE> ) {
459        #-- delta int
460        if ( /^([-]?\d+)$/ ) {
461            if ( $1 < 0 ) {
462                $tot ++;
463            }
464            elsif ( $1 == 0 ) {
465                $align[4] = ($tot - $sim) / $tot * 100.0;
466                push @$aref, [ @align ];
467                $tot = 0;
468            }
469            next;
470        }
471
472        #-- alignment header
473        if ( /^(\d+) (\d+) (\d+) (\d+) \d+ (\d+) \d+$/ ) {
474            if ( $tot == 0 ) {
475                @align = ($1, $2, $3, $4, 0, $lenR, $lenQ, $idR, $idQ);
476                $tot = abs($1 - $2) + 1;
477                $sim = $5;
478                if ( $ispromer ) { $tot /= 3.0; }
479                next;
480            }
481            #-- drop to default
482        }
483
484        #-- sequence header
485        if ( /^>(\S+) (\S+) (\d+) (\d+)$/ ) {
486            ($idR, $idQ, $lenR, $lenQ) = ($1, $2, $3, $4);
487            $tot = 0;
488            next;
489        }
490
491        #-- default
492        die "ERROR: Could not parse $OPT_Mfile\n$_";
493    }
494
495    close (MFILE)
496        or print STDERR "WARNING: Trouble closing $OPT_Mfile, $!\n";
497}
498
499
500#------------------------------------------------------------ ParseCluster ----#
501sub ParseCluster ($)
502{
503    my $aref = shift;
504
505    print STDERR "Reading cluster file $OPT_Mfile\n";
506
507    open (MFILE, "<$OPT_Mfile")
508        or die "ERROR: Could not open $OPT_Mfile, $!\n";
509
510    my @align;
511    my ($dR, $dQ, $len);
512    my ($lenR, $lenQ, $idR, $idQ);
513
514    $_ = <MFILE>;
515    $_ = <MFILE>;
516
517    while ( <MFILE> ) {
518        #-- match
519        if ( /^\s+(\d+)\s+(\d+)\s+(\d+)\s+\S+\s+\S+$/ ) {
520            @align = ($1, $1, $2, $2, 100, $lenR, $lenQ, $idR, $idQ);
521            $len = $3 - 1;
522            $align[1] += $dR == 1 ? $len : -$len;
523            $align[3] += $dQ == 1 ? $len : -$len;
524            push @$aref, [ @align ];
525            next;
526        }
527
528        #-- cluster header
529        if ( /^[ \-][1-3] [ \-][1-3]$/ ) {
530            $dR = /^-/ ? -1 : 1;
531            $dQ = /-[1-3]$/ ? -1 : 1;
532            next;
533        }
534
535        #-- sequence header
536        if ( /^>(\S+) (\S+) (\d+) (\d+)$/ ) {
537            ($idR, $idQ, $lenR, $lenQ) = ($1, $2, $3, $4);
538            next;
539        }
540
541        #-- default
542        die "ERROR: Could not parse $OPT_Mfile\n$_";
543    }
544
545    close (MFILE)
546        or print STDERR "WARNING: Trouble closing $OPT_Mfile, $!\n";
547}
548
549
550#------------------------------------------------------------- ParseMummer ----#
551sub ParseMummer ($)
552{
553    my $aref = shift;
554
555    print STDERR "Reading mummer file $OPT_Mfile (use mummer -c)\n";
556
557    open (MFILE, "<$OPT_Mfile")
558        or die "ERROR: Could not open $OPT_Mfile, $!\n";
559
560    my @align;
561    my ($dQ, $len);
562    my ($lenQ, $idQ);
563
564    while ( <MFILE> ) {
565        #-- 3 column match
566        if ( /^\s+(\d+)\s+(\d+)\s+(\d+)$/ ) {
567            @align = ($1, $1, $2, $2, 100, 0, $lenQ, "REF", $idQ);
568            $len = $3 - 1;
569            $align[1] += $len;
570            $align[3] += $dQ == 1 ? $len : -$len;
571            push @$aref, [ @align ];
572            next;
573        }
574
575        #-- 4 column match
576        if ( /^\s+(\S+)\s+(\d+)\s+(\d+)\s+(\d+)$/ ) {
577            @align = ($2, $2, $3, $3, 100, 0, $lenQ, $1, $idQ);
578            $len = $4 - 1;
579            $align[1] += $len;
580            $align[3] += $dQ == 1 ? $len : -$len;
581            push @$aref, [ @align ];
582            next;
583        }
584
585        #-- sequence header
586        if ( /^> (\S+)/ ) {
587            $idQ = $1;
588            $dQ = /^> \S+ Reverse/ ? -1 : 1;
589            $lenQ = /Len = (\d+)/ ? $1 : 0;
590            next;
591        }
592
593        #-- default
594        die "ERROR: Could not parse $OPT_Mfile\n$_";
595    }
596
597    close (MFILE)
598        or print STDERR "WARNING: Trouble closing $OPT_Mfile, $!\n";
599}
600
601
602#------------------------------------------------------------- ParseTiling ----#
603sub ParseTiling ($)
604{
605    my $aref = shift;
606
607    print STDERR "Reading tiling file $OPT_Mfile\n";
608
609    open (MFILE, "<$OPT_Mfile")
610        or die "ERROR: Could not open $OPT_Mfile, $!\n";
611
612    my @align;
613    my ($dR, $dQ, $len);
614    my ($lenR, $lenQ, $idR, $idQ);
615
616    while ( <MFILE> ) {
617        #-- tile
618        if ( /^(\S+)\s+\S+\s+\S+\s+(\d+)\s+\S+\s+(\S+)\s+([+-])\s+(\S+)$/ ) {
619            @align = ($1, $1, 1, 1, $3, $lenR, $2, $idR, $5);
620            $len = $2 - 1;
621            $align[1] += $len;
622            $align[($4 eq "-" ? 2 : 3)] += $len;
623            push @$aref, [ @align ];
624            next;
625        }
626
627        #-- sequence header
628        if ( /^>(\S+) (\d+) bases$/ ) {
629            ($idR, $lenR) = ($1, $2);
630            next;
631        }
632
633        #-- default
634        die "ERROR: Could not parse $OPT_Mfile\n$_";
635    }
636
637    close (MFILE)
638        or print STDERR "WARNING: Trouble closing $OPT_Mfile, $!\n";
639
640    if ( ! defined $OPT_coverage ) { $OPT_coverage = 1; }
641}
642
643
644#--------------------------------------------------------------- LayoutIDs ----#
645# For each reference and query sequence, find the set of alignments that
646# produce the heaviest (both in non-redundant coverage and percent
647# identity) alignment subset of each sequence using a modified version
648# of the longest increasing subset algorithm. Let R be the union of all
649# reference LIS subsets, and Q be the union of all query LIS
650# subsets. Let S be the intersection of R and Q. Using this LIS subset,
651# recursively span reference and query sequences by their smaller
652# counterparts until all spanning sequences have been placed. The goal
653# is to cluster all the "major" alignment information along the main
654# diagonal for easy viewing and interpretation.
655sub LayoutIDs ($$)
656{
657    my $rref = shift;
658    my $qref = shift;
659
660    my %rc;          # chains of qry seqs needed to span each ref
661    my %qc;          # chains of ref seqs needed to span each qry
662    #  {idR} -> [ placed, len, {idQ} -> [ \slope, \loR, \hiR, \loQ, \hiQ ] ]
663    #  {idQ} -> [ placed, len, {idR} -> [ \slope, \loQ, \hiQ, \loR, \hiR ] ]
664
665    my @rl;          # oo of ref seqs
666    my @ql;          # oo of qry seqs
667    #  [ [idR, slope] ]
668    #  [ [idQ, slope] ]
669
670    #-- get the filtered alignments
671    open (BTAB, "$BIN_DIR/show-coords -B $OPT_Dfile |")
672        or die "ERROR: Could not open show-coords pipe, $!\n";
673
674    my @align;
675    my ($sR, $eR, $sQ, $eQ, $lenR, $lenQ, $idR, $idQ);
676    my ($loR, $hiR, $loQ, $hiQ);
677    my ($dR, $dQ, $slope);
678    while ( <BTAB> ) {
679        chomp;
680        @align = split "\t";
681        if ( scalar @align != 21 ) {
682            die "ERROR: Could not read show-coords pipe, invalid btab format\n";
683        }
684
685        $sR   = $align[8];   $eR   = $align[9];
686        $sQ   = $align[6];   $eQ   = $align[7];
687        $lenR = $align[18];  $lenQ = $align[2];
688        $idR  = $align[5];   $idQ  = $align[0];
689
690        #-- skip it if not on include list
691        if ( !exists $rref->{$idR} || !exists $qref->{$idQ} ) { next; }
692
693        #-- get orientation of both alignments and alignment slope
694        $dR = $sR < $eR ? 1 : -1;
695        $dQ = $sQ < $eQ ? 1 : -1;
696        $slope = $dR == $dQ ? 1 : -1;
697
698        #-- get lo's and hi's
699        $loR = $dR == 1 ? $sR : $eR;
700        $hiR = $dR == 1 ? $eR : $sR;
701
702        $loQ = $dQ == 1 ? $sQ : $eQ;
703        $hiQ = $dQ == 1 ? $eQ : $sQ;
704
705        if ($OPT_ONLY_USE_FATTEST)
706        {
707          #-- Check to see if there is another better alignment
708          if (exists $qc{$idQ})
709          {
710            my ($oldR) = keys %{$qc{$idQ}[2]};
711            my $val = $qc{$idQ}[2]{$oldR};
712
713            if (${$val->[4]} - ${$val->[3]} > $hiR - $loR)
714            {
715              #-- Old alignment is better, skip this one
716              next;
717            }
718            else
719            {
720              #-- This alignment is better, prune old alignment
721              delete $rc{$oldR}[2]{$idQ};
722              delete $qc{$idQ};
723            }
724          }
725        }
726
727        #-- initialize
728        if ( !exists $rc{$idR} ) { $rc{$idR} = [ 0, $lenR, { } ]; }
729        if ( !exists $qc{$idQ} ) { $qc{$idQ} = [ 0, $lenQ, { } ]; }
730
731        #-- if no alignments for these two exist OR
732        #-- this alignment is bigger than the current
733        if ( !exists $rc{$idR}[2]{$idQ} || !exists $qc{$idQ}[2]{$idR} ||
734             $hiR - $loR >
735             ${$rc{$idR}[2]{$idQ}[2]} - ${$rc{$idR}[2]{$idQ}[1]} ) {
736
737            #-- rc and qc reference these anonymous values
738            my $aref = [ $slope, $loR, $hiR, $loQ, $hiQ ];
739
740            #-- rc is ordered [ slope, loR, hiR, loQ, hiQ ]
741            #-- qc is ordered [ slope, loQ, hiQ, loR, hiR ]
742            $rc{$idR}[2]{$idQ}[0] = $qc{$idQ}[2]{$idR}[0] = \$aref->[0];
743            $rc{$idR}[2]{$idQ}[1] = $qc{$idQ}[2]{$idR}[3] = \$aref->[1];
744            $rc{$idR}[2]{$idQ}[2] = $qc{$idQ}[2]{$idR}[4] = \$aref->[2];
745            $rc{$idR}[2]{$idQ}[3] = $qc{$idQ}[2]{$idR}[1] = \$aref->[3];
746            $rc{$idR}[2]{$idQ}[4] = $qc{$idQ}[2]{$idR}[2] = \$aref->[4];
747        }
748    }
749
750    close (BTAB)
751        or print STDERR "WARNING: Trouble closing show-coords pipe, $!\n";
752
753    #-- recursively span sequences to generate the layout
754    foreach $idR ( sort { $rc{$b}[1] <=> $rc{$a}[1] } keys %rc ) {
755        SpanXwY ($idR, \%rc, \@rl, \%qc, \@ql);
756    }
757
758    #-- undefine the current offsets
759    foreach $idR ( keys %{$rref} ) { undef $rref->{$idR}[0]; }
760    foreach $idQ ( keys %{$qref} ) { undef $qref->{$idQ}[0]; }
761
762    #-- redefine the offsets according to the new layout
763    my $roff = 0;
764    foreach my $r ( @rl ) {
765        $idR = $r->[0];
766        $rref->{$idR}[0] = $roff;
767        $rref->{$idR}[2] = $r->[1];
768        $roff += $rref->{$idR}[1] - 1;
769    }
770    #-- append the guys left out of the layout
771    foreach $idR ( keys %{$rref} ) {
772        if ( !defined $rref->{$idR}[0] ) {
773            $rref->{$idR}[0] = $roff;
774            $roff += $rref->{$idR}[1] - 1;
775        }
776    }
777
778    #-- redefine the offsets according to the new layout
779    my $qoff = 0;
780    foreach my $q ( @ql ) {
781        $idQ = $q->[0];
782        $qref->{$idQ}[0] = $qoff;
783        $qref->{$idQ}[2] = $q->[1];
784        $qoff += $qref->{$idQ}[1] - 1;
785    }
786    #-- append the guys left out of the layout
787    foreach $idQ ( keys %{$qref} ) {
788        if ( !defined $qref->{$idQ}[0] ) {
789            $qref->{$idQ}[0] = $qoff;
790            $qoff += $qref->{$idQ}[1] - 1;
791        }
792    }
793}
794
795
796#----------------------------------------------------------------- SpanXwY ----#
797sub SpanXwY ($$$$$) {
798    my $x   = shift;   # idX
799    my $xcr = shift;   # xc ref
800    my $xlr = shift;   # xl ref
801    my $ycr = shift;   # yc ref
802    my $ylr = shift;   # yl ref
803
804    my @post;
805    foreach my $y ( sort { ${$xcr->{$x}[2]{$a}[1]} <=> ${$xcr->{$x}[2]{$b}[1]} }
806                    keys %{$xcr->{$x}[2]} ) {
807
808        #-- skip if already placed (RECURSION BASE)
809        if ( $ycr->{$y}[0] ) { next; }
810        else { $ycr->{$y}[0] = 1; }
811
812        #-- get len and slope info for y
813        my $len = $ycr->{$y}[1];
814        my $slope = ${$xcr->{$x}[2]{$y}[0]};
815
816        #-- if we need to flip, reverse complement all y records
817        if ( $slope == -1 ) {
818            foreach my $xx ( keys %{$ycr->{$y}[2]} ) {
819                ${$ycr->{$y}[2]{$xx}[0]} *= -1;
820
821                my $loy = ${$ycr->{$y}[2]{$xx}[1]};
822                my $hiy = ${$ycr->{$y}[2]{$xx}[2]};
823                ${$ycr->{$y}[2]{$xx}[1]} = $len - $hiy + 1;
824                ${$ycr->{$y}[2]{$xx}[2]} = $len - $loy + 1;
825            }
826        }
827
828        #-- place y
829        push @{$ylr}, [ $y, $slope ];
830
831        #-- RECURSE if y > x, else save for later
832        if ( $len > $xcr->{$x}[1] ) { SpanXwY ($y, $ycr, $ylr, $xcr, $xlr); }
833        else { push @post, $y; }
834    }
835
836    #-- RECURSE for all y < x
837    foreach my $y ( @post ) { SpanXwY ($y, $ycr, $ylr, $xcr, $xlr); }
838}
839
840
841#---------------------------------------------------------------- PlotData ----#
842sub PlotData ($$$)
843{
844    my $aref = shift;
845    my $rref = shift;
846    my $qref = shift;
847
848    print STDERR "Writing plot files $OPT_Ffile, $OPT_Rfile",
849    (defined $OPT_Hfile ? ", $OPT_Hfile\n" : "\n");
850
851    open (FFILE, ">$OPT_Ffile")
852        or die "ERROR: Could not open $OPT_Ffile, $!\n";
853    print FFILE "#-- forward hits sorted by %sim\n0 0 0\n0 0 0\n\n\n";
854
855    open (RFILE, ">$OPT_Rfile")
856        or die "ERROR: Could not open $OPT_Rfile, $!\n";
857    print RFILE "#-- reverse hits sorted by %sim\n0 0 0\n0 0 0\n\n\n";
858
859    if ( defined $OPT_Hfile ) {
860        open (HFILE, ">$OPT_Hfile")
861            or die "ERROR: Could not open $OPT_Hfile, $!\n";
862        print HFILE "#-- highlighted hits sorted by %sim\n0 0 0\n0 0 0\n\n\n";
863    }
864
865    my $fh;
866    my $align;
867    my $isplotted;
868    my $ismultiref;
869    my $ismultiqry;
870    my ($plenR, $plenQ, $pidR, $pidQ);
871
872    #-- for each alignment sorted by ascending identity
873    foreach $align ( sort { $a->[4] <=> $b->[4] } @$aref ) {
874
875        my ($sR, $eR, $sQ, $eQ, $sim, $lenR, $lenQ, $idR, $idQ) = @$align;
876
877        if ( ! defined $pidR ) {
878            ($plenR, $plenQ, $pidR, $pidQ) = ($lenR, $lenQ, $idR, $idQ);
879        }
880
881        #-- set the sequence offset, length, direction, etc...
882        my ($refoff, $reflen, $refdir);
883        my ($qryoff, $qrylen, $qrydir);
884
885        if ( %$rref ) {
886            #-- skip reference sequence or set atts from hash
887            if ( !exists ($rref->{$idR}) ) { next; }
888            else { ($refoff, $reflen, $refdir) = @{$rref->{$idR}}; }
889        }
890        else {
891            #-- no reference hash, so default atts
892            ($refoff, $reflen, $refdir) = (0, $lenR, 1);
893        }
894
895        if ( %$qref ) {
896            #-- skip query sequence or set atts from hash
897            if ( !exists ($qref->{$idQ}) ) { next; }
898            else { ($qryoff, $qrylen, $qrydir) = @{$qref->{$idQ}}; }
899        }
900        else {
901            #-- no query hash, so default atts
902            ($qryoff, $qrylen, $qrydir) = (0, $lenQ, 1);
903        }
904
905        #-- get the orientation right
906        if ( $refdir == -1 ) {
907            $sR = $reflen - $sR + 1;
908            $eR = $reflen - $eR + 1;
909        }
910        if ( $qrydir == -1 ) {
911            $sQ = $qrylen - $sQ + 1;
912            $eQ = $qrylen - $eQ + 1;
913        }
914
915        #-- forward file, reverse file, highlight file
916        my @fha;
917
918        if ( defined $OPT_breaklen &&
919             ( ($sR - 1 > $OPT_breaklen &&
920                $sQ - 1 > $OPT_breaklen &&
921                $reflen - $sR > $OPT_breaklen &&
922                $qrylen - $sQ > $OPT_breaklen)
923               ||
924               ($eR - 1 > $OPT_breaklen &&
925                $eQ - 1 > $OPT_breaklen &&
926                $reflen - $eR > $OPT_breaklen &&
927                $qrylen - $eQ > $OPT_breaklen) ) ) {
928            push @fha, \*HFILE;
929        }
930
931        push @fha, (($sR < $eR) == ($sQ < $eQ) ? \*FFILE : \*RFILE);
932
933        #-- plot it
934        $sR += $refoff; $eR += $refoff;
935        $sQ += $qryoff; $eQ += $qryoff;
936
937        if ( $OPT_coverage ) {
938            foreach $fh ( @fha ) {
939                print $fh
940                    "$sR 10 $sim\n", "$eR 10 $sim\n\n\n",
941                    "$sR $sim 0\n", "$eR $sim 0\n\n\n";
942            }
943        }
944        else {
945            foreach $fh ( @fha ) {
946                print $fh "$sR $sQ $sim\n", "$eR $eQ $sim\n\n\n";
947            }
948        }
949
950        #-- set some flags
951        if ( !$ismultiref && $idR ne $pidR ) { $ismultiref = 1; }
952        if ( !$ismultiqry && $idQ ne $pidQ ) { $ismultiqry = 1; }
953        if ( !$isplotted ) { $isplotted = 1; }
954    }
955
956
957    #-- highlight the SNPs
958    if ( defined $OPT_SNP ) {
959
960        print STDERR "Determining SNPs from sequence and alignment data\n";
961
962        open (SNPS, "$BIN_DIR/show-snps -H -T -l $OPT_Mfile |")
963            or die "ERROR: Could not open show-snps pipe, $!\n";
964
965        my @snps;
966        my ($pR, $pQ, $lenR, $lenQ, $idR, $idQ);
967        while ( <SNPS> ) {
968            chomp;
969            @snps = split "\t";
970            if ( scalar @snps != 14 ) {
971                die "ERROR: Could not read show-snps pipe, invalid format\n";
972            }
973
974            $pR   = $snps[0];   $pQ   = $snps[3];
975            $lenR = $snps[8];   $lenQ = $snps[9];
976            $idR  = $snps[12];  $idQ  = $snps[13];
977
978            #-- set the sequence offset, length, direction, etc...
979            my ($refoff, $reflen, $refdir);
980            my ($qryoff, $qrylen, $qrydir);
981
982            if ( %$rref ) {
983                #-- skip reference sequence or set atts from hash
984                if ( !exists ($rref->{$idR}) ) { next; }
985                else { ($refoff, $reflen, $refdir) = @{$rref->{$idR}}; }
986            }
987            else {
988                #-- no reference hash, so default atts
989                ($refoff, $reflen, $refdir) = (0, $lenR, 1);
990            }
991
992            if ( %$qref ) {
993                #-- skip query sequence or set atts from hash
994                if ( !exists ($qref->{$idQ}) ) { next; }
995                else { ($qryoff, $qrylen, $qrydir) = @{$qref->{$idQ}}; }
996            }
997            else {
998                #-- no query hash, so default atts
999                ($qryoff, $qrylen, $qrydir) = (0, $lenQ, 1);
1000            }
1001
1002            #-- get the orientation right
1003            if ( $refdir == -1 ) { $pR = $reflen - $pR + 1; }
1004            if ( $qrydir == -1 ) { $pQ = $qrylen - $pQ + 1; }
1005
1006            #-- plot it
1007            $pR += $refoff;
1008            $pQ += $qryoff;
1009
1010            if ( $OPT_coverage ) {
1011                print HFILE "$pR 10 0\n", "$pR 10 0\n\n\n",
1012            }
1013            else {
1014                print HFILE "$pR $pQ 0\n", "$pR $pQ 0\n\n\n";
1015            }
1016        }
1017
1018        close (SNPS)
1019            or print STDERR "WARNING: Trouble closing show-snps pipe, $!\n";
1020    }
1021
1022
1023    close (FFILE)
1024        or print STDERR "WARNING: Trouble closing $OPT_Ffile, $!\n";
1025
1026    close (RFILE)
1027        or print STDERR "WARNING: Trouble closing $OPT_Rfile, $!\n";
1028
1029    if ( defined $OPT_Hfile ) {
1030        close (HFILE)
1031            or print STDERR "WARNING: Trouble closing $OPT_Hfile, $!\n";
1032    }
1033
1034
1035    if ( !%$rref ) {
1036        if ( $ismultiref ) {
1037            print STDERR
1038                "WARNING: Multiple ref sequences overlaid, try -R or -r\n";
1039        }
1040        elsif ( defined $pidR ) {
1041            $rref->{$pidR} = [ 0, $plenR, 1 ];
1042        }
1043    }
1044
1045    if ( !%$qref ) {
1046        if ( $ismultiqry && !$OPT_coverage ) {
1047            print STDERR
1048                "WARNING: Multiple qry sequences overlaid, try -Q, -q or -c\n";
1049        }
1050        elsif ( defined $pidQ ) {
1051            $qref->{$pidQ} = [ 0, $plenQ, 1 ];
1052        }
1053    }
1054
1055    if ( !$isplotted ) {
1056        die "ERROR: No alignment data to plot\n";
1057    }
1058}
1059
1060
1061#----------------------------------------------------------------- WriteGP ----#
1062sub WriteGP ($$)
1063{
1064    my $rref = shift;
1065    my $qref = shift;
1066
1067    print STDERR "Writing gnuplot script $OPT_Gfile\n";
1068
1069    open (GFILE, ">$OPT_Gfile")
1070        or die "ERROR: Could not open $OPT_Gfile, $!\n";
1071
1072    my ($FWD, $REV, $HLT) = (1, 2, 3);
1073    my $SIZE = $TERMSIZE{$OPT_terminal}{$OPT_size};
1074
1075    #-- terminal specific stuff
1076    my ($P_TERM, $P_SIZE, %P_PS, %P_LW);
1077    foreach ( $OPT_terminal ) {
1078        /^$X11/    and do {
1079            $P_TERM = $OPT_gpstatus == 0 ?
1080                "$X11 font \"$FFACE,$FSIZE\"" : "$X11";
1081
1082            %P_PS = ( $FWD => 1.0, $REV => 1.0, $HLT => 1.0 );
1083
1084            %P_LW = $OPT_coverage || $OPT_color ?
1085                ( $FWD => 3.0, $REV => 3.0, $HLT => 3.0 ) :
1086                ( $FWD => 2.0, $REV => 2.0, $HLT => 2.0 );
1087
1088            $P_SIZE = $OPT_coverage ?
1089                "set size 1,1" :
1090                "set size 1,1";
1091
1092            last;
1093        };
1094
1095        /^$PS/     and do {
1096            $P_TERM = defined $OPT_color && $OPT_color == 0 ?
1097                "$PS monochrome" : "$PS color";
1098            $P_TERM .= $OPT_gpstatus == 0 ?
1099                " solid \"$FFACE\" $FSIZE" : " solid \"$FFACE\" $FSIZE";
1100
1101            %P_PS = ( $FWD => 0.5, $REV => 0.5, $HLT => 0.5 );
1102
1103            %P_LW = $OPT_coverage || $OPT_color ?
1104                ( $FWD => 4.0, $REV => 4.0, $HLT => 4.0 ) :
1105                ( $FWD => 2.0, $REV => 2.0, $HLT => 2.0 );
1106
1107            $P_SIZE = $OPT_coverage ?
1108                "set size ".(1.0 * $SIZE).",".(0.5 * $SIZE) :
1109                "set size ".(1.0 * $SIZE).",".(1.0 * $SIZE);
1110
1111            last;
1112        };
1113
1114        /^$PNG/    and do {
1115            $P_TERM = $OPT_gpstatus == 0 ?
1116                "$PNG tiny size $SIZE,$SIZE" : "$PNG small";
1117            if ( defined $OPT_color && $OPT_color == 0 ) {
1118                $P_TERM .= " xffffff x000000 x000000";
1119                $P_TERM .= " x000000 x000000 x000000";
1120                $P_TERM .= " x000000 x000000 x000000";
1121            }
1122
1123            %P_PS = ( $FWD => 1.0, $REV => 1.0, $HLT => 1.0 );
1124
1125            %P_LW = $OPT_coverage || $OPT_color ?
1126                ( $FWD => 3.0, $REV => 3.0, $HLT => 3.0 ) :
1127                ( $FWD => 3.0, $REV => 3.0, $HLT => 3.0 );
1128
1129            $P_SIZE = $OPT_coverage ?
1130                "set size 1,.375" :
1131                "set size 1,1";
1132
1133            last;
1134        };
1135
1136        die "ERROR: Don't know how to initialize terminal, $OPT_terminal\n";
1137    }
1138
1139    #-- plot commands
1140    my ($P_WITH, $P_FORMAT, $P_LS, $P_KEY, %P_PT, %P_LT);
1141
1142    %P_PT = ( $FWD => 6, $REV => 6, $HLT => 6 );
1143    %P_LT = defined $OPT_Hfile ?
1144        ( $FWD => 2, $REV => 2, $HLT => 1 ) :
1145        ( $FWD => 1, $REV => 3, $HLT => 2 );
1146
1147    $P_WITH = $OPT_coverage || $OPT_color ? "w l" : "w lp";
1148
1149    $P_FORMAT = "set format \"$TFORMAT\"";
1150    if ( $OPT_gpstatus == 0 ) {
1151        $P_LS = "set style line";
1152        $P_KEY = "unset key";
1153        $P_FORMAT .= "\nset mouse format \"$TFORMAT\"";
1154        $P_FORMAT .= "\nset mouse mouseformat \"$MFORMAT\"";
1155        $P_FORMAT .= "\nif(GPVAL_VERSION < 5) { set mouse clipboardformat \"$MFORMAT\" }";
1156    }
1157    else {
1158        $P_LS = "set linestyle";
1159        $P_KEY = "set nokey";
1160    }
1161
1162
1163    my @refk = keys (%$rref);
1164    my @qryk = keys (%$qref);
1165    my ($xrange, $yrange);
1166    my ($xlabel, $ylabel);
1167    my ($tic, $dir);
1168    my $border = 0;
1169
1170    #-- terminal header and output
1171    print GFILE "set terminal $P_TERM\n";
1172
1173    if ( defined $OPT_Pfile ) {
1174        print GFILE "set output \"$OPT_Pfile\"\n";
1175    }
1176
1177    if ( defined $OPT_title ) {
1178        print GFILE "set title \"$OPT_title\"\n";
1179    }
1180
1181    #-- set tics, determine labels, ranges (ref)
1182    if ( scalar (@refk) == 1 ) {
1183        $xlabel = $refk[0];
1184        $xrange = $rref->{$xlabel}[1];
1185    }
1186    else {
1187        $xrange = 0;
1188        print GFILE "set xtics rotate \( \\\n";
1189        foreach $xlabel ( sort { $rref->{$a}[0] <=> $rref->{$b}[0] } @refk ) {
1190            $xrange += $rref->{$xlabel}[1];
1191            $tic = $rref->{$xlabel}[0] + 1;
1192            $dir = ($rref->{$xlabel}[2] == 1) ? "" : "*";
1193            print GFILE " \"$dir$xlabel\" $tic.0, \\\n";
1194        }
1195        print GFILE " \"\" $xrange \\\n\)\n";
1196        $xlabel = "REF";
1197    }
1198    if ( $xrange == 0 ) { $xrange = "*"; }
1199
1200    #-- set tics, determine labels, ranges (qry)
1201    if ( $OPT_coverage ) {
1202        $ylabel = "%SIM";
1203        $yrange = 110;
1204    }
1205    elsif ( scalar (@qryk) == 1 ) {
1206        $ylabel = $qryk[0];
1207        $yrange = $qref->{$ylabel}[1];
1208    }
1209    else {
1210        $yrange = 0;
1211        print GFILE "set ytics \( \\\n";
1212        foreach $ylabel ( sort { $qref->{$a}[0] <=> $qref->{$b}[0] } @qryk ) {
1213            $yrange += $qref->{$ylabel}[1];
1214            $tic = $qref->{$ylabel}[0] + 1;
1215            $dir = ($qref->{$ylabel}[2] == 1) ? "" : "*";
1216            print GFILE " \"$dir$ylabel\" $tic.0, \\\n";
1217        }
1218        print GFILE " \"\" $yrange \\\n\)\n";
1219        $ylabel = "QRY";
1220    }
1221    if ( $yrange == 0 ) { $yrange = "*"; }
1222
1223    #-- determine borders
1224    if ( $xrange ne "*" && scalar (@refk) == 1 ) { $border |= 10; }
1225    if ( $yrange ne "*" && scalar (@qryk) == 1 ) { $border |= 5; }
1226    if ( $OPT_coverage ) { $border |= 5; }
1227
1228    #-- grid, labels, border
1229    print GFILE
1230        "$P_SIZE\n",
1231        "set grid\n",
1232        "$P_KEY\n",
1233        "set border $border\n",
1234        "set tics scale 0\n",
1235        "set xlabel \"$xlabel\"\n",
1236        "set ylabel \"$ylabel\"\n",
1237        "$P_FORMAT\n";
1238
1239    #-- ranges
1240    if ( defined $OPT_xrange ) { print GFILE "set xrange $OPT_xrange\n"; }
1241    else                       { print GFILE "set xrange [1:$xrange]\n"; }
1242
1243    if ( defined $OPT_yrange ) { print GFILE "set yrange $OPT_yrange\n"; }
1244    else                       { print GFILE "set yrange [1:$yrange]\n"; }
1245
1246    #-- if %sim plot
1247    if ( $OPT_color ) {
1248        print GFILE
1249            "set zrange [0:100]\n",
1250            "set colorbox default\n",
1251            "set cblabel \"%similarity\"\n",
1252            "set cbrange [0:100]\n",
1253            "set cbtics 20\n",
1254            "set pm3d map\n",
1255            "set palette model RGB defined ( \\\n",
1256            "  0 \"#000000\", \\\n",
1257            "  4 \"#DD00DD\", \\\n",
1258            "  6 \"#0000DD\", \\\n",
1259            "  7 \"#00DDDD\", \\\n",
1260            "  8 \"#00DD00\", \\\n",
1261            "  9 \"#DDDD00\", \\\n",
1262            " 10 \"#DD0000\"  \\\n)\n";
1263    }
1264
1265    foreach my $s ( ($FWD, $REV, $HLT) ) {
1266        my $ss = "$P_LS $s ";
1267        $ss .= $OPT_color ? " palette" : " lt $P_LT{$s}";
1268        $ss .= " lw $P_LW{$s}";
1269        if ( ! $OPT_coverage || $s == $HLT ) {
1270            $ss .= " pt $P_PT{$s} ps $P_PS{$s}";
1271        }
1272        print GFILE "$ss\n";
1273    }
1274
1275    #-- plot it
1276    print GFILE
1277        ($OPT_color ? "splot \\\n" : "plot \\\n");
1278    print GFILE
1279        " \"$OPT_Ffile\" title \"FWD\" $P_WITH ls $FWD, \\\n",
1280        " \"$OPT_Rfile\" title \"REV\" $P_WITH ls $REV",
1281        (! defined $OPT_Hfile ? "\n" :
1282         ", \\\n \"$OPT_Hfile\" title \"HLT\" w lp ls $HLT");
1283
1284    #-- interactive mode
1285    if ( $OPT_terminal eq $X11 ) {
1286        print GFILE "\n",
1287        "print \"-- INTERACTIVE MODE --\"\n",
1288        "print \"consult gnuplot docs for command list\"\n",
1289        "print \"mouse 1: coords to clipboard\"\n",
1290        "print \"mouse 2: mark on plot\"\n",
1291        "print \"mouse 3: zoom box\"\n",
1292        "print \"'h' for help in plot window\"\n",
1293        "print \"enter to exit\"\n",
1294        "pause -1\n";
1295    }
1296
1297    close (GFILE)
1298        or print STDERR "WARNING: Trouble closing $OPT_Gfile, $!\n";
1299}
1300
1301
1302#------------------------------------------------------------------- RunGP ----#
1303sub RunGP ( )
1304{
1305    if ( defined $OPT_Pfile ) {
1306        print STDERR "Rendering plot $OPT_Pfile\n";
1307    }
1308    else {
1309        print STDERR "Rendering plot to screen\n";
1310    }
1311
1312    my $cmd = $GNUPLOT_EXE;
1313
1314    #-- x11 specifics
1315    if ( $OPT_terminal eq $X11 ) {
1316        my $size = $TERMSIZE{$OPT_terminal}{$OPT_size};
1317        $cmd .= " -geometry ${size}x";
1318        if ( $OPT_coverage ) { $size = sprintf ("%.0f", $size * .375); }
1319        $cmd .= "${size}+0+0 -title mummerplot";
1320
1321        if ( defined $OPT_color && $OPT_color == 0 ) {
1322            $cmd .= " -mono";
1323            $cmd .= " -xrm 'gnuplot*line1Dashes: 0'";
1324            $cmd .= " -xrm 'gnuplot*line2Dashes: 0'";
1325            $cmd .= " -xrm 'gnuplot*line3Dashes: 0'";
1326        }
1327
1328        if ( $OPT_rv ) {
1329            $cmd .= " -rv";
1330            $cmd .= " -xrm 'gnuplot*background: black'";
1331            $cmd .= " -xrm 'gnuplot*textColor: white'";
1332            $cmd .= " -xrm 'gnuplot*borderColor: white'";
1333            $cmd .= " -xrm 'gnuplot*axisColor: white'";
1334        }
1335    }
1336
1337    $cmd .= " $OPT_Gfile";
1338
1339    system ($cmd)
1340        and print STDERR "WARNING: Unable to run '$cmd', $!\n";
1341}
1342
1343
1344#---------------------------------------------------------------- ListenGP ----#
1345sub ListenGP($$)
1346{
1347    my $rref = shift;
1348    my $qref = shift;
1349
1350    my ($refc, $qryc);
1351    my ($refid, $qryid);
1352    my ($rsock, $qsock);
1353    my $oldclip = "";
1354
1355    #-- get IDs sorted by offset
1356    my @refo = sort { $rref->{$a}[0] <=> $rref->{$b}[0] } keys %$rref;
1357    my @qryo = sort { $qref->{$a}[0] <=> $qref->{$b}[0] } keys %$qref;
1358
1359    #-- attempt to connect sockets
1360    if ( $OPT_rport ) {
1361        $rsock = IO::Socket::INET->new("localhost:$OPT_rport")
1362            or print STDERR "WARNING: Could not connect to rport $OPT_rport\n";
1363    }
1364
1365    if ( $OPT_qport ) {
1366        $qsock = IO::Socket::INET->new("localhost:$OPT_qport")
1367            or print STDERR "WARNING: Could not connect to qport $OPT_qport\n";
1368    }
1369
1370    #-- while parent still exists
1371    while ( getppid != 1 ) {
1372
1373        #-- query the clipboard
1374        $_ = `xclip -o -silent -selection primary`;
1375        if ( $? >> 8 ) {
1376            die "WARNING: Unable to query clipboard with xclip\n";
1377        }
1378
1379        #-- if cliboard has changed and contains a coordinate
1380        if ( $_ ne $oldclip && (($refc, $qryc) = /^\[(\d+), (\d+)\]/) ) {
1381
1382            $oldclip = $_;
1383
1384            #-- translate the reference position
1385            $refid = "NULL";
1386            for ( my $i = 0; $i < (scalar @refo); ++ $i ) {
1387                my $aref = $rref->{$refo[$i]};
1388                if ( $i == $#refo || $aref->[0] + $aref->[1] > $refc ) {
1389                    $refid = $refo[$i];
1390                    $refc -= $aref->[0];
1391                    if ( $aref->[2] == -1 ) {
1392                        $refc = $aref->[1] - $refc + 1;
1393                    }
1394                    last;
1395                }
1396            }
1397
1398            #-- translate the query position
1399            $qryid = "NULL";
1400            for ( my $i = 0; $i < (scalar @qryo); ++ $i ) {
1401                my $aref = $qref->{$qryo[$i]};
1402                if ( $i == $#qryo || $aref->[0] + $aref->[1] > $qryc ) {
1403                    $qryid = $qryo[$i];
1404                    $qryc -= $aref->[0];
1405                    if ( $aref->[2] == -1 ) {
1406                        $qryc = $aref->[1] - $qryc + 1;
1407                    }
1408                    last;
1409                }
1410            }
1411
1412            #-- print the info to stdout and socket
1413            print "$refid\t$qryid\t$refc\t$qryc\n";
1414
1415            if ( $rsock ) {
1416                print $rsock "contig I$refid $refc\n";
1417                print "sent \"contig I$refid $refc\" to $OPT_rport\n";
1418            }
1419            if ( $qsock ) {
1420                print $qsock "contig I$qryid $qryc\n";
1421                print "sent \"contig I$qryid $qryc\" to $OPT_qport\n";
1422            }
1423        }
1424
1425        #-- sleep for half second
1426        select undef, undef, undef, .5;
1427    }
1428
1429    exit (0);
1430}
1431
1432
1433#------------------------------------------------------------ ParseOptions ----#
1434sub ParseOptions ( )
1435{
1436    my ($opt_small, $opt_medium, $opt_large);
1437    my ($opt_ps, $opt_x11, $opt_png);
1438    my $cnt;
1439
1440    #-- Get options
1441    my $err = $tigr -> TIGR_GetOptions
1442        (
1443         "b|breaklen:i" => \$OPT_breaklen,
1444         "color!"       => \$OPT_color,
1445         "c|coverage!"  => \$OPT_coverage,
1446         "f|filter!"    => \$OPT_filter,
1447         "l|layout!"    => \$OPT_layout,
1448         "p|prefix=s"   => \$OPT_prefix,
1449         "rv"           => \$OPT_rv,
1450         "r|IdR=s"      => \$OPT_IdR,
1451         "q|IdQ=s"      => \$OPT_IdQ,
1452         "R|Rfile=s"    => \$OPT_IDRfile,
1453         "Q|Qfile=s"    => \$OPT_IDQfile,
1454         "rport=i"      => \$OPT_rport,
1455         "qport=i"      => \$OPT_qport,
1456         "s|size=s"     => \$OPT_size,
1457         "S|SNP"        => \$OPT_SNP,
1458         "t|terminal=s" => \$OPT_terminal,
1459         "title=s"      => \$OPT_title,
1460         "x|xrange=s"   => \$OPT_xrange,
1461         "y|yrange=s"   => \$OPT_yrange,
1462         "x11"          => \$opt_x11,
1463         "postscript"   => \$opt_ps,
1464         "png"          => \$opt_png,
1465         "small"        => \$opt_small,
1466         "medium"       => \$opt_medium,
1467         "large"        => \$opt_large,
1468         "fat"          => \$OPT_ONLY_USE_FATTEST,
1469         );
1470
1471    if ( !$err  ||  scalar (@ARGV) != 1 ) {
1472        $tigr -> printUsageInfo( );
1473        die "Try '$0 -h' for more information.\n";
1474    }
1475
1476    $cnt = 0;
1477    if ( $opt_png ) { $OPT_terminal = $PNG; $cnt ++; }
1478    if ( $opt_ps  ) { $OPT_terminal = $PS;  $cnt ++; }
1479    if ( $opt_x11 ) { $OPT_terminal = $X11; $cnt ++; }
1480    if ( $cnt > 1 ) {
1481        print STDERR
1482            "WARNING: Multiple terminals not allowed, using '$OPT_terminal'\n";
1483    }
1484
1485    $cnt = 0;
1486    if ( $opt_large  ) { $OPT_size = $LARGE;  $cnt ++; }
1487    if ( $opt_medium ) { $OPT_size = $MEDIUM; $cnt ++; }
1488    if ( $opt_small  ) { $OPT_size = $SMALL;  $cnt ++; }
1489    if ( $cnt > 1 ) {
1490        print STDERR
1491            "WARNING: Multiple sizes now allowed, using '$OPT_size'\n";
1492    }
1493
1494    #-- Check that status of gnuplot
1495    $OPT_gpstatus = system ("gnuplot --version");
1496
1497    if ( $OPT_gpstatus == -1 ) {
1498        print STDERR
1499            "WARNING: Could not find gnuplot, plot will not be rendered\n";
1500    }
1501    elsif ( $OPT_gpstatus ) {
1502        print STDERR
1503            "WARNING: Using outdated gnuplot, use v4 or later for best results\n";
1504
1505        if ( $OPT_color ) {
1506            print STDERR
1507                "WARNING: Turning off --color option for compatibility\n";
1508            undef $OPT_color;
1509        }
1510
1511        if ( $OPT_terminal eq $PNG  &&  $OPT_size ne $SMALL ) {
1512            print STDERR
1513                "WARNING: Turning off --size option for compatibility\n";
1514            $OPT_size = $SMALL;
1515        }
1516    }
1517
1518    #-- Check options
1519    if ( !exists $TERMSIZE{$OPT_terminal} ) {
1520        die "ERROR: Invalid terminal type, $OPT_terminal\n";
1521    }
1522
1523    if ( !exists $TERMSIZE{$OPT_terminal}{$OPT_size} ) {
1524        die "ERROR: Invalid terminal size, $OPT_size\n";
1525    }
1526
1527    if ( $OPT_xrange ) {
1528        $OPT_xrange =~ tr/,/:/;
1529        $OPT_xrange =~ /^\[\d+:\d+\]$/
1530            or die "ERROR: Invalid xrange format, $OPT_xrange\n";
1531    }
1532
1533    if ( $OPT_yrange ) {
1534        $OPT_yrange =~ tr/,/:/;
1535        $OPT_yrange =~ /^\[\d+:\d+\]$/
1536            or die "ERROR: Invalid yrange format, $OPT_yrange\n";
1537    }
1538
1539    #-- Set file names
1540    $OPT_Mfile = $ARGV[0];
1541    $tigr->isReadableFile ($OPT_Mfile)
1542        or die "ERROR: Could not read $OPT_Mfile, $!\n";
1543
1544    $OPT_Ffile = $OPT_prefix . $SUFFIX{$FWDPLOT};
1545    $tigr->isWritableFile ($OPT_Ffile) or $tigr->isCreatableFile ($OPT_Ffile)
1546        or die "ERROR: Could not write $OPT_Ffile, $!\n";
1547
1548    $OPT_Rfile = $OPT_prefix . $SUFFIX{$REVPLOT};
1549    $tigr->isWritableFile ($OPT_Rfile) or $tigr->isCreatableFile ($OPT_Rfile)
1550        or die "ERROR: Could not write $OPT_Rfile, $!\n";
1551
1552    if ( defined $OPT_breaklen || defined $OPT_SNP ) {
1553        $OPT_Hfile = $OPT_prefix . $SUFFIX{$HLTPLOT};
1554        $tigr->isWritableFile($OPT_Hfile) or $tigr->isCreatableFile($OPT_Hfile)
1555            or die "ERROR: Could not write $OPT_Hfile, $!\n";
1556    }
1557
1558    if ($OPT_ONLY_USE_FATTEST)
1559    {
1560      $OPT_layout = 1;
1561    }
1562
1563    if ( $OPT_filter || $OPT_layout ) {
1564        $OPT_Dfile = $OPT_prefix . $SUFFIX{$FILTER};
1565        $tigr->isWritableFile($OPT_Dfile) or $tigr->isCreatableFile($OPT_Dfile)
1566            or die "ERROR: Could not write $OPT_Dfile, $!\n";
1567    }
1568
1569    $OPT_Gfile = $OPT_prefix . $SUFFIX{$GNUPLOT};
1570    $tigr->isWritableFile ($OPT_Gfile) or $tigr->isCreatableFile ($OPT_Gfile)
1571        or die "ERROR: Could not write $OPT_Gfile, $!\n";
1572
1573    if ( exists $SUFFIX{$OPT_terminal} ) {
1574        $OPT_Pfile = $OPT_prefix . $SUFFIX{$OPT_terminal};
1575        $tigr->isWritableFile($OPT_Pfile) or $tigr->isCreatableFile($OPT_Pfile)
1576            or die "ERROR: Could not write $OPT_Pfile, $!\n";
1577    }
1578
1579    if ( defined $OPT_IDRfile ) {
1580        $tigr->isReadableFile ($OPT_IDRfile)
1581            or die "ERROR: Could not read $OPT_IDRfile, $!\n";
1582    }
1583
1584    if ( defined $OPT_IDQfile ) {
1585        $tigr->isReadableFile ($OPT_IDQfile)
1586            or die "ERROR: Could not read $OPT_IDQfile, $!\n";
1587    }
1588
1589    if ( (defined $OPT_rport || defined $OPT_qport) &&
1590         ($OPT_terminal ne $X11 || $OPT_gpstatus ) ) {
1591        print STDERR
1592            "WARNING: Port options available only for v4.0 X11 plots\n";
1593        undef $OPT_rport;
1594        undef $OPT_qport;
1595    }
1596
1597
1598    if ( defined $OPT_color && defined $OPT_Hfile ) {
1599        print STDERR
1600            "WARNING: Turning off --color option so highlighting is visible\n";
1601        undef $OPT_color;
1602    }
1603}
1604