jpayne@69: #!/usr/bin/env perl
jpayne@69: 
jpayne@69: ################################################################################
jpayne@69: #   Programmer: Adam M Phillippy, The Institute for Genomic Research
jpayne@69: #         File: mummerplot
jpayne@69: #         Date: 01 / 08 / 03
jpayne@69: #               01 / 06 / 05 rewritten (v3.0)
jpayne@69: #  
jpayne@69: #        Usage:
jpayne@69: #    mummerplot  [options]  <match file>
jpayne@69: # 
jpayne@69: #                Try 'mummerplot -h' for more information.
jpayne@69: # 
jpayne@69: #      Purpose: To generate a gnuplot plot for the display of mummer, nucmer,
jpayne@69: #               promer, and show-tiling alignments.
jpayne@69: # 
jpayne@69: ################################################################################
jpayne@69: 
jpayne@69: use lib "/mnt/c/Users/crash/Documents/BobLiterman/CSP2_Galaxy/CSP2/CSP2_env/env-d9b9114564458d9d-741b3de822f2aaca6c6caa4325c4afce/opt/mummer-3.23/scripts";
jpayne@69: use Foundation;
jpayne@69: use strict;
jpayne@69: use IO::Socket;
jpayne@69: 
jpayne@69: my $BIN_DIR = "/mnt/c/Users/crash/Documents/BobLiterman/CSP2_Galaxy/CSP2/CSP2_env/env-d9b9114564458d9d-741b3de822f2aaca6c6caa4325c4afce/opt/mummer-3.23";
jpayne@69: my $SCRIPT_DIR = "/mnt/c/Users/crash/Documents/BobLiterman/CSP2_Galaxy/CSP2/CSP2_env/env-d9b9114564458d9d-741b3de822f2aaca6c6caa4325c4afce/opt/mummer-3.23/scripts";
jpayne@69: 
jpayne@69: 
jpayne@69: #================================================================= Globals ====#
jpayne@69: #-- terminal types
jpayne@69: my $X11    = "x11";
jpayne@69: my $PS     = "postscript";
jpayne@69: my $PNG    = "png";
jpayne@69: 
jpayne@69: #-- terminal sizes
jpayne@69: my $SMALL  = "small";
jpayne@69: my $MEDIUM = "medium";
jpayne@69: my $LARGE  = "large";
jpayne@69: 
jpayne@69: my %TERMSIZE =
jpayne@69:     (
jpayne@69:      $X11 => { $SMALL => 500, $MEDIUM => 700,  $LARGE => 900  }, # screen pix
jpayne@69:      $PS  => { $SMALL => 1,   $MEDIUM => 2,    $LARGE => 3    }, # pages
jpayne@69:      $PNG => { $SMALL => 800, $MEDIUM => 1024, $LARGE => 1400 }  # image pix
jpayne@69:      );
jpayne@69: 
jpayne@69: #-- terminal format
jpayne@69: my $FFACE    = "Courier";
jpayne@69: my $FSIZE    = "8";
jpayne@69: my $TFORMAT  = "%.0f";
jpayne@69: my $MFORMAT  = "[%.0f, %.0f]";
jpayne@69: 
jpayne@69: #-- output suffixes
jpayne@69: my $FILTER  = "filter";
jpayne@69: my $FWDPLOT = "fplot";
jpayne@69: my $REVPLOT = "rplot";
jpayne@69: my $HLTPLOT = "hplot";
jpayne@69: my $GNUPLOT = "gnuplot";
jpayne@69: 
jpayne@69: my %SUFFIX =
jpayne@69:     (
jpayne@69:      $FILTER  => ".filter",
jpayne@69:      $FWDPLOT => ".fplot",
jpayne@69:      $REVPLOT => ".rplot",
jpayne@69:      $HLTPLOT => ".hplot",
jpayne@69:      $GNUPLOT => ".gp",
jpayne@69:      $PS      => ".ps",
jpayne@69:      $PNG     => ".png"
jpayne@69:      );
jpayne@69: 
jpayne@69: 
jpayne@69: #================================================================= Options ====#
jpayne@69: my $OPT_breaklen;                  # -b option
jpayne@69: my $OPT_color;                     # --[no]color option
jpayne@69: my $OPT_coverage;                  # --[no]coverage option
jpayne@69: my $OPT_filter;                    # -f option
jpayne@69: my $OPT_layout;                    # -l option
jpayne@69: my $OPT_prefix    = "out";         # -p option
jpayne@69: my $OPT_rv;                        # --rv option
jpayne@69: my $OPT_terminal  = $X11;          # -t option
jpayne@69: my $OPT_IdR;                       # -r option
jpayne@69: my $OPT_IdQ;                       # -q option
jpayne@69: my $OPT_IDRfile;                   # -R option
jpayne@69: my $OPT_IDQfile;                   # -Q option
jpayne@69: my $OPT_rport;                     # -rport option
jpayne@69: my $OPT_qport;                     # -qport option
jpayne@69: my $OPT_size      = $SMALL;        # -small, -medium, -large
jpayne@69: my $OPT_SNP;                       # -S option
jpayne@69: my $OPT_xrange;                    # -x option
jpayne@69: my $OPT_yrange;                    # -y option
jpayne@69: my $OPT_title;                     # -title option
jpayne@69: 
jpayne@69: my $OPT_Mfile;                     # match file
jpayne@69: my $OPT_Dfile;                     # delta filter file
jpayne@69: my $OPT_Ffile;                     # .fplot output
jpayne@69: my $OPT_Rfile;                     # .rplot output
jpayne@69: my $OPT_Hfile;                     # .hplot output
jpayne@69: my $OPT_Gfile;                     # .gp output
jpayne@69: my $OPT_Pfile;                     # .ps .png output
jpayne@69: 
jpayne@69: my $OPT_gpstatus;                  # gnuplot status
jpayne@69: 
jpayne@69: my $OPT_ONLY_USE_FATTEST;          # Only use fattest alignment for layout
jpayne@69: 
jpayne@69: 
jpayne@69: #============================================================== Foundation ====#
jpayne@69: my $VERSION = '3.5';
jpayne@69: 
jpayne@69: my $USAGE = qq~
jpayne@69:   USAGE: mummerplot  [options]  <match file>
jpayne@69:     ~;
jpayne@69: 
jpayne@69: my $HELP = qq~
jpayne@69:   USAGE: mummerplot  [options]  <match file>
jpayne@69: 
jpayne@69:   DESCRIPTION:
jpayne@69:     mummerplot generates plots of alignment data produced by mummer, nucmer,
jpayne@69:     promer or show-tiling by using the GNU gnuplot utility. After generating
jpayne@69:     the appropriate scripts and datafiles, mummerplot will attempt to run
jpayne@69:     gnuplot to generate the plot. If this attempt fails, a warning will be
jpayne@69:     output and the resulting .gp and .[frh]plot files will remain so that the
jpayne@69:     user may run gnuplot independently. If the attempt succeeds, either an x11
jpayne@69:     window will be spawned or an additional output file will be generated
jpayne@69:     (.ps or .png depending on the selected terminal). Feel free to edit the
jpayne@69:     resulting gnuplot script (.gp) and rerun gnuplot to change line thinkness,
jpayne@69:     labels, colors, plot size etc.
jpayne@69: 
jpayne@69:   MANDATORY:
jpayne@69:     match file      Set the alignment input to 'match file'
jpayne@69:                     Valid inputs are from mummer, nucmer, promer and
jpayne@69:                     show-tiling (.out, .cluster, .delta and .tiling)
jpayne@69: 
jpayne@69:   OPTIONS:
jpayne@69:     -b|breaklen     Highlight alignments with breakpoints further than
jpayne@69:                     breaklen nucleotides from the nearest sequence end
jpayne@69:     --[no]color     Color plot lines with a percent similarity gradient or
jpayne@69:                     turn off all plot color (default color by match dir)
jpayne@69:                     If the plot is very sparse, edit the .gp script to plot
jpayne@69:                     with 'linespoints' instead of 'lines'
jpayne@69:     -c
jpayne@69:     --[no]coverage  Generate a reference coverage plot (default for .tiling)
jpayne@69:     --depend        Print the dependency information and exit
jpayne@69:     -f
jpayne@69:     --filter        Only display .delta alignments which represent the "best"
jpayne@69:                     hit to any particular spot on either sequence, i.e. a
jpayne@69:                     one-to-one mapping of reference and query subsequences
jpayne@69:     -h
jpayne@69:     --help          Display help information and exit
jpayne@69:     -l
jpayne@69:     --layout        Layout a .delta multiplot in an intelligible fashion,
jpayne@69:                     this option requires the -R -Q options
jpayne@69:     --fat           Layout sequences using fattest alignment only
jpayne@69:     -p|prefix       Set the prefix of the output files (default '$OPT_prefix')
jpayne@69:     -rv             Reverse video for x11 plots
jpayne@69:     -r|IdR          Plot a particular reference sequence ID on the X-axis
jpayne@69:     -q|IdQ          Plot a particular query sequence ID on the Y-axis
jpayne@69:     -R|Rfile        Plot an ordered set of reference sequences from Rfile
jpayne@69:     -Q|Qfile        Plot an ordered set of query sequences from Qfile
jpayne@69:                     Rfile/Qfile Can either be the original DNA multi-FastA
jpayne@69:                     files or lists of sequence IDs, lens and dirs [ /+/-]
jpayne@69:     -r|rport        Specify the port to send reference ID and position on
jpayne@69:                     mouse double click in X11 plot window
jpayne@69:     -q|qport        Specify the port to send query IDs and position on mouse
jpayne@69:                     double click in X11 plot window
jpayne@69:     -s|size         Set the output size to small, medium or large
jpayne@69:                     --small --medium --large (default '$OPT_size')
jpayne@69:     -S
jpayne@69:     --SNP           Highlight SNP locations in each alignment
jpayne@69:     -t|terminal     Set the output terminal to x11, postscript or png
jpayne@69:                     --x11 --postscript --png (default '$OPT_terminal')
jpayne@69:     -t|title        Specify the gnuplot plot title (default none)
jpayne@69:     -x|xrange       Set the xrange for the plot '[min:max]'
jpayne@69:     -y|yrange       Set the yrange for the plot '[min:max]'
jpayne@69:     -V
jpayne@69:     --version       Display the version information and exit
jpayne@69:     ~;
jpayne@69: 
jpayne@69: my @DEPEND =
jpayne@69:     (
jpayne@69:      "$SCRIPT_DIR/Foundation.pm",
jpayne@69:      "$BIN_DIR/delta-filter",
jpayne@69:      "$BIN_DIR/show-coords",
jpayne@69:      "$BIN_DIR/show-snps",
jpayne@69:      "gnuplot"
jpayne@69:      );
jpayne@69: 
jpayne@69: my $tigr = new TIGR::Foundation
jpayne@69:     or die "ERROR: TIGR::Foundation could not be initialized\n";
jpayne@69: 
jpayne@69: $tigr -> setVersionInfo ($VERSION);
jpayne@69: $tigr -> setUsageInfo ($USAGE);
jpayne@69: $tigr -> setHelpInfo ($HELP);
jpayne@69: $tigr -> addDependInfo (@DEPEND);
jpayne@69: 
jpayne@69: 
jpayne@69: #=========================================================== Function Decs ====#
jpayne@69: sub GetParseFunc( );
jpayne@69: 
jpayne@69: sub ParseIDs($$);
jpayne@69: 
jpayne@69: sub ParseDelta($);
jpayne@69: sub ParseCluster($);
jpayne@69: sub ParseMummer($);
jpayne@69: sub ParseTiling($);
jpayne@69: 
jpayne@69: sub LayoutIDs($$);
jpayne@69: sub SpanXwY ($$$$$);
jpayne@69: 
jpayne@69: sub PlotData($$$);
jpayne@69: sub WriteGP($$);
jpayne@69: sub RunGP( );
jpayne@69: sub ListenGP($$);
jpayne@69: 
jpayne@69: sub ParseOptions( );
jpayne@69: 
jpayne@69: 
jpayne@69: #=========================================================== Function Defs ====#
jpayne@69: MAIN:
jpayne@69: {
jpayne@69:     my @aligns;                # (sR eR sQ eQ sim lenR lenQ idR idQ)
jpayne@69:     my %refs;                  # (id => (off, len, [1/-1]))
jpayne@69:     my %qrys;                  # (id => (off, len, [1/-1]))
jpayne@69: 
jpayne@69:     #-- Get the command line options (sets OPT_ global vars)
jpayne@69:     ParseOptions( );
jpayne@69: 
jpayne@69: 
jpayne@69:     #-- Get the alignment type
jpayne@69:     my $parsefunc = GetParseFunc( );
jpayne@69: 
jpayne@69:     if ( $parsefunc != \&ParseDelta &&
jpayne@69:          ($OPT_filter || $OPT_layout || $OPT_SNP) ) {
jpayne@69:         print STDERR "WARNING: -f -l -S only work with delta input\n";
jpayne@69:         undef $OPT_filter;
jpayne@69:         undef $OPT_layout;
jpayne@69:         undef $OPT_SNP;
jpayne@69:     }
jpayne@69: 
jpayne@69:     #-- Parse the reference and query IDs
jpayne@69:     if    ( defined $OPT_IdR ) { $refs{$OPT_IdR} = [ 0, 0, 1 ]; }
jpayne@69:     elsif ( defined $OPT_IDRfile ) {
jpayne@69:         ParseIDs ($OPT_IDRfile, \%refs);
jpayne@69:     }
jpayne@69: 
jpayne@69:     if    ( defined $OPT_IdQ ) { $qrys{$OPT_IdQ} = [ 0, 0, 1 ]; }
jpayne@69:     elsif ( defined $OPT_IDQfile ) {
jpayne@69:         ParseIDs ($OPT_IDQfile, \%qrys);
jpayne@69:     }
jpayne@69: 
jpayne@69: 
jpayne@69:     #-- Filter the alignments
jpayne@69:     if ( $OPT_filter || $OPT_layout ) {
jpayne@69:         print STDERR "Writing filtered delta file $OPT_Dfile\n";
jpayne@69:         system ("$BIN_DIR/delta-filter -r -q $OPT_Mfile > $OPT_Dfile")
jpayne@69:             and die "ERROR: Could not run delta-filter, $!\n";
jpayne@69:         if ( $OPT_filter ) { $OPT_Mfile = $OPT_Dfile; }
jpayne@69:     }
jpayne@69: 
jpayne@69: 
jpayne@69:     #-- Parse the alignment data
jpayne@69:     $parsefunc->(\@aligns);
jpayne@69: 
jpayne@69:     
jpayne@69:     #-- Layout the alignment data if requested
jpayne@69:     if ( $OPT_layout ) {
jpayne@69:         if ( scalar (keys %refs) || scalar (keys %qrys) ) {
jpayne@69:             LayoutIDs (\%refs, \%qrys);
jpayne@69:         }
jpayne@69:         else {
jpayne@69:             print STDERR "WARNING: --layout option only works with -R or -Q\n";
jpayne@69:             undef $OPT_layout;
jpayne@69:         }
jpayne@69:     }
jpayne@69: 
jpayne@69: 
jpayne@69:     #-- Plot the alignment data
jpayne@69:     PlotData (\@aligns, \%refs, \%qrys);
jpayne@69: 
jpayne@69: 
jpayne@69:     #-- Write the gnuplot script
jpayne@69:     WriteGP (\%refs, \%qrys);
jpayne@69: 
jpayne@69: 
jpayne@69:     #-- Run gnuplot script and fork a clipboard listener
jpayne@69:     unless ( $OPT_gpstatus == -1 ) {
jpayne@69: 
jpayne@69:         my $child = 1;
jpayne@69:         if ( $OPT_gpstatus == 0 && $OPT_terminal eq $X11 ) {
jpayne@69:             print STDERR "Forking mouse listener\n";
jpayne@69:             $child = fork;
jpayne@69:         }
jpayne@69: 
jpayne@69:         #-- parent runs gnuplot
jpayne@69:         if ( $child ) {
jpayne@69:             RunGP( );
jpayne@69:             kill 1, $child;
jpayne@69:         }
jpayne@69:         #-- child listens to clipboard
jpayne@69:         elsif ( defined $child ) {
jpayne@69:             ListenGP(\%refs, \%qrys);
jpayne@69:         }
jpayne@69:         else {
jpayne@69:             print STDERR "WARNING: Could not fork mouse listener\n";
jpayne@69:         }
jpayne@69:     }
jpayne@69: 
jpayne@69:     exit (0);
jpayne@69: }
jpayne@69: 
jpayne@69: 
jpayne@69: #------------------------------------------------------------ GetParseFunc ----#
jpayne@69: sub GetParseFunc ( )
jpayne@69: {
jpayne@69:     my $fref;
jpayne@69: 
jpayne@69:     open (MFILE, "<$OPT_Mfile")
jpayne@69:         or die "ERROR: Could not open $OPT_Mfile, $!\n";
jpayne@69: 
jpayne@69:     $_ = <MFILE>;
jpayne@69:     if ( !defined ) { die "ERROR: Could not read $OPT_Mfile, File is empty\n" }
jpayne@69: 
jpayne@69:   SWITCH: {
jpayne@69:       #-- tiling
jpayne@69:       if ( /^>\S+ \d+ bases/ ) {
jpayne@69:           $fref = \&ParseTiling;
jpayne@69:           last SWITCH;
jpayne@69:       }
jpayne@69: 
jpayne@69:       #-- mummer
jpayne@69:       if ( /^> \S+/ ) {
jpayne@69:           $fref = \&ParseMummer;
jpayne@69:           last SWITCH;
jpayne@69:       }
jpayne@69: 
jpayne@69:       #-- nucmer/promer
jpayne@69:       if ( /^(\S+) (\S+)/ ) {
jpayne@69:           if ( ! defined $OPT_IDRfile ) {
jpayne@69:               $OPT_IDRfile = $1;
jpayne@69:           }
jpayne@69:           if ( ! defined $OPT_IDQfile ) {
jpayne@69:               $OPT_IDQfile = $2;
jpayne@69:           }
jpayne@69: 
jpayne@69:           $_ = <MFILE>;
jpayne@69:           if ( (defined)  &&  (/^NUCMER$/  ||  /^PROMER$/) ) {
jpayne@69:               $_ = <MFILE>;   # sequence header
jpayne@69:               $_ = <MFILE>;   # alignment header
jpayne@69:               if ( !defined ) {
jpayne@69:                   $fref = \&ParseDelta;
jpayne@69:                   last SWITCH;
jpayne@69:               }
jpayne@69:               elsif ( /^\d+ \d+ \d+ \d+ \d+ \d+ \d+$/ ) {
jpayne@69:                   $fref = \&ParseDelta;
jpayne@69:                   last SWITCH;
jpayne@69:               }
jpayne@69:               elsif ( /^[ \-][1-3] [ \-][1-3]$/ ) {
jpayne@69:                   $fref = \&ParseCluster;
jpayne@69:                   last SWITCH;
jpayne@69:               }
jpayne@69:           }
jpayne@69:       }
jpayne@69:       
jpayne@69:       #-- default
jpayne@69:       die "ERROR: Could not read $OPT_Mfile, Unrecognized file type\n";
jpayne@69:   }
jpayne@69: 
jpayne@69:     close (MFILE)
jpayne@69:         or print STDERR "WARNING: Trouble closing $OPT_Mfile, $!\n";
jpayne@69: 
jpayne@69:     return $fref;
jpayne@69: }
jpayne@69: 
jpayne@69: 
jpayne@69: #---------------------------------------------------------------- ParseIDs ----#
jpayne@69: sub ParseIDs ($$)
jpayne@69: {
jpayne@69:     my $file = shift;
jpayne@69:     my $href = shift;
jpayne@69: 
jpayne@69:     open (IDFILE, "<$file")
jpayne@69:         or print STDERR "WARNING: Could not open $file, $!\n";
jpayne@69: 
jpayne@69:     my $dir;
jpayne@69:     my $aref;
jpayne@69:     my $isfasta;
jpayne@69:     my $offset = 0;
jpayne@69:     while ( <IDFILE> ) {
jpayne@69:         #-- Ignore blank lines
jpayne@69:         if ( /^\s*$/ ) { next; }
jpayne@69: 
jpayne@69:         #-- FastA header
jpayne@69:         if ( /^>(\S+)/ ) {
jpayne@69:             if ( exists $href->{$1} ) {
jpayne@69:                 print STDERR "WARNING: Duplicate sequence '$1' ignored\n";
jpayne@69:                 undef $aref;
jpayne@69:                 next;
jpayne@69:             }
jpayne@69: 
jpayne@69:             if ( !$isfasta ) { $isfasta = 1; }
jpayne@69:             if ( defined $aref ) { $offset += $aref->[1] - 1; }
jpayne@69: 
jpayne@69:             $aref = [ $offset, 0, 1 ];
jpayne@69:             $href->{$1} = $aref;
jpayne@69:             next;
jpayne@69:         }
jpayne@69:         
jpayne@69:         #-- FastA sequence
jpayne@69:         if ( $isfasta  &&  /^\S+$/ ) {
jpayne@69:             if ( defined $aref ) { $aref->[1] += (length) - 1; }
jpayne@69:             next;
jpayne@69:         }
jpayne@69: 
jpayne@69:         #-- ID len dir
jpayne@69:         if ( !$isfasta  &&  /^(\S+)\s+(\d+)\s+([+-]?)$/ ) {
jpayne@69:             if ( exists $href->{$1} ) {
jpayne@69:                 print STDERR "WARNING: Duplicate sequence '$1' ignored\n";
jpayne@69:                 undef $aref;
jpayne@69:                 next;
jpayne@69:             }
jpayne@69: 
jpayne@69:             $dir = (defined $3 && $3 eq "-") ? -1 : 1;
jpayne@69:             $aref = [ $offset, $2, $dir ];
jpayne@69:             $offset += $2 - 1;
jpayne@69:             $href->{$1} = $aref;
jpayne@69:             next;
jpayne@69:         }
jpayne@69: 
jpayne@69:         #-- default
jpayne@69:         print STDERR "WARNING: Could not parse $file\n$_";
jpayne@69:         undef %$href;
jpayne@69:         last;
jpayne@69:     }
jpayne@69: 
jpayne@69:     close (IDFILE)
jpayne@69:         or print STDERR "WARNING: Trouble closing $file, $!\n";
jpayne@69: }
jpayne@69: 
jpayne@69: 
jpayne@69: #-------------------------------------------------------------- ParseDelta ----#
jpayne@69: sub ParseDelta ($)
jpayne@69: {
jpayne@69:     my $aref = shift;
jpayne@69: 
jpayne@69:     print STDERR "Reading delta file $OPT_Mfile\n";
jpayne@69: 
jpayne@69:     open (MFILE, "<$OPT_Mfile")
jpayne@69:         or die "ERROR: Could not open $OPT_Mfile, $!\n";
jpayne@69: 
jpayne@69:     my @align;
jpayne@69:     my $ispromer;
jpayne@69:     my ($sim, $tot);
jpayne@69:     my ($lenR, $lenQ, $idR, $idQ);
jpayne@69: 
jpayne@69:     $_ = <MFILE>;
jpayne@69:     $_ = <MFILE>;
jpayne@69:     $ispromer = /^PROMER/;
jpayne@69: 
jpayne@69:     while ( <MFILE> ) {
jpayne@69:         #-- delta int
jpayne@69:         if ( /^([-]?\d+)$/ ) {
jpayne@69:             if ( $1 < 0 ) {
jpayne@69:                 $tot ++;
jpayne@69:             }
jpayne@69:             elsif ( $1 == 0 ) {
jpayne@69:                 $align[4] = ($tot - $sim) / $tot * 100.0;
jpayne@69:                 push @$aref, [ @align ];
jpayne@69:                 $tot = 0;
jpayne@69:             }
jpayne@69:             next;
jpayne@69:         }
jpayne@69: 
jpayne@69:         #-- alignment header
jpayne@69:         if ( /^(\d+) (\d+) (\d+) (\d+) \d+ (\d+) \d+$/ ) {
jpayne@69:             if ( $tot == 0 ) {
jpayne@69:                 @align = ($1, $2, $3, $4, 0, $lenR, $lenQ, $idR, $idQ);
jpayne@69:                 $tot = abs($1 - $2) + 1;
jpayne@69:                 $sim = $5;
jpayne@69:                 if ( $ispromer ) { $tot /= 3.0; }
jpayne@69:                 next;
jpayne@69:             }
jpayne@69:             #-- drop to default
jpayne@69:         }
jpayne@69: 
jpayne@69:         #-- sequence header
jpayne@69:         if ( /^>(\S+) (\S+) (\d+) (\d+)$/ ) {
jpayne@69:             ($idR, $idQ, $lenR, $lenQ) = ($1, $2, $3, $4);
jpayne@69:             $tot = 0;
jpayne@69:             next;
jpayne@69:         }
jpayne@69: 
jpayne@69:         #-- default
jpayne@69:         die "ERROR: Could not parse $OPT_Mfile\n$_";
jpayne@69:     }
jpayne@69: 
jpayne@69:     close (MFILE)
jpayne@69:         or print STDERR "WARNING: Trouble closing $OPT_Mfile, $!\n";
jpayne@69: }
jpayne@69: 
jpayne@69: 
jpayne@69: #------------------------------------------------------------ ParseCluster ----#
jpayne@69: sub ParseCluster ($)
jpayne@69: {
jpayne@69:     my $aref = shift;
jpayne@69: 
jpayne@69:     print STDERR "Reading cluster file $OPT_Mfile\n";
jpayne@69: 
jpayne@69:     open (MFILE, "<$OPT_Mfile")
jpayne@69:         or die "ERROR: Could not open $OPT_Mfile, $!\n";
jpayne@69: 
jpayne@69:     my @align;
jpayne@69:     my ($dR, $dQ, $len);
jpayne@69:     my ($lenR, $lenQ, $idR, $idQ);
jpayne@69: 
jpayne@69:     $_ = <MFILE>;
jpayne@69:     $_ = <MFILE>;
jpayne@69: 
jpayne@69:     while ( <MFILE> ) {
jpayne@69:         #-- match
jpayne@69:         if ( /^\s+(\d+)\s+(\d+)\s+(\d+)\s+\S+\s+\S+$/ ) {
jpayne@69:             @align = ($1, $1, $2, $2, 100, $lenR, $lenQ, $idR, $idQ);
jpayne@69:             $len = $3 - 1;
jpayne@69:             $align[1] += $dR == 1 ? $len : -$len;
jpayne@69:             $align[3] += $dQ == 1 ? $len : -$len;
jpayne@69:             push @$aref, [ @align ];
jpayne@69:             next;
jpayne@69:         }
jpayne@69: 
jpayne@69:         #-- cluster header
jpayne@69:         if ( /^[ \-][1-3] [ \-][1-3]$/ ) {
jpayne@69:             $dR = /^-/ ? -1 : 1;
jpayne@69:             $dQ = /-[1-3]$/ ? -1 : 1;
jpayne@69:             next;
jpayne@69:         }
jpayne@69: 
jpayne@69:         #-- sequence header
jpayne@69:         if ( /^>(\S+) (\S+) (\d+) (\d+)$/ ) {
jpayne@69:             ($idR, $idQ, $lenR, $lenQ) = ($1, $2, $3, $4);
jpayne@69:             next;
jpayne@69:         }
jpayne@69: 
jpayne@69:         #-- default
jpayne@69:         die "ERROR: Could not parse $OPT_Mfile\n$_";
jpayne@69:     }
jpayne@69: 
jpayne@69:     close (MFILE)
jpayne@69:         or print STDERR "WARNING: Trouble closing $OPT_Mfile, $!\n";
jpayne@69: }
jpayne@69: 
jpayne@69: 
jpayne@69: #------------------------------------------------------------- ParseMummer ----#
jpayne@69: sub ParseMummer ($)
jpayne@69: {
jpayne@69:     my $aref = shift;
jpayne@69: 
jpayne@69:     print STDERR "Reading mummer file $OPT_Mfile (use mummer -c)\n";
jpayne@69: 
jpayne@69:     open (MFILE, "<$OPT_Mfile")
jpayne@69:         or die "ERROR: Could not open $OPT_Mfile, $!\n";
jpayne@69: 
jpayne@69:     my @align;
jpayne@69:     my ($dQ, $len);
jpayne@69:     my ($lenQ, $idQ);
jpayne@69: 
jpayne@69:     while ( <MFILE> ) {
jpayne@69:         #-- 3 column match
jpayne@69:         if ( /^\s+(\d+)\s+(\d+)\s+(\d+)$/ ) {
jpayne@69:             @align = ($1, $1, $2, $2, 100, 0, $lenQ, "REF", $idQ);
jpayne@69:             $len = $3 - 1;
jpayne@69:             $align[1] += $len;
jpayne@69:             $align[3] += $dQ == 1 ? $len : -$len;
jpayne@69:             push @$aref, [ @align ];
jpayne@69:             next;
jpayne@69:         }
jpayne@69: 
jpayne@69:         #-- 4 column match
jpayne@69:         if ( /^\s+(\S+)\s+(\d+)\s+(\d+)\s+(\d+)$/ ) {
jpayne@69:             @align = ($2, $2, $3, $3, 100, 0, $lenQ, $1, $idQ);
jpayne@69:             $len = $4 - 1;
jpayne@69:             $align[1] += $len;
jpayne@69:             $align[3] += $dQ == 1 ? $len : -$len;
jpayne@69:             push @$aref, [ @align ];
jpayne@69:             next;
jpayne@69:         }
jpayne@69: 
jpayne@69:         #-- sequence header
jpayne@69:         if ( /^> (\S+)/ ) {
jpayne@69:             $idQ = $1;
jpayne@69:             $dQ = /^> \S+ Reverse/ ? -1 : 1;
jpayne@69:             $lenQ = /Len = (\d+)/ ? $1 : 0;
jpayne@69:             next;
jpayne@69:         }
jpayne@69: 
jpayne@69:         #-- default
jpayne@69:         die "ERROR: Could not parse $OPT_Mfile\n$_";
jpayne@69:     }
jpayne@69: 
jpayne@69:     close (MFILE)
jpayne@69:         or print STDERR "WARNING: Trouble closing $OPT_Mfile, $!\n";
jpayne@69: }
jpayne@69: 
jpayne@69: 
jpayne@69: #------------------------------------------------------------- ParseTiling ----#
jpayne@69: sub ParseTiling ($)
jpayne@69: {
jpayne@69:     my $aref = shift;
jpayne@69: 
jpayne@69:     print STDERR "Reading tiling file $OPT_Mfile\n";
jpayne@69: 
jpayne@69:     open (MFILE, "<$OPT_Mfile")
jpayne@69:         or die "ERROR: Could not open $OPT_Mfile, $!\n";
jpayne@69: 
jpayne@69:     my @align;
jpayne@69:     my ($dR, $dQ, $len);
jpayne@69:     my ($lenR, $lenQ, $idR, $idQ);
jpayne@69: 
jpayne@69:     while ( <MFILE> ) {
jpayne@69:         #-- tile
jpayne@69:         if ( /^(\S+)\s+\S+\s+\S+\s+(\d+)\s+\S+\s+(\S+)\s+([+-])\s+(\S+)$/ ) {
jpayne@69:             @align = ($1, $1, 1, 1, $3, $lenR, $2, $idR, $5);
jpayne@69:             $len = $2 - 1;
jpayne@69:             $align[1] += $len;
jpayne@69:             $align[($4 eq "-" ? 2 : 3)] += $len;
jpayne@69:             push @$aref, [ @align ];
jpayne@69:             next;
jpayne@69:         }
jpayne@69: 
jpayne@69:         #-- sequence header
jpayne@69:         if ( /^>(\S+) (\d+) bases$/ ) {
jpayne@69:             ($idR, $lenR) = ($1, $2);
jpayne@69:             next;
jpayne@69:         }
jpayne@69: 
jpayne@69:         #-- default
jpayne@69:         die "ERROR: Could not parse $OPT_Mfile\n$_";
jpayne@69:     }
jpayne@69: 
jpayne@69:     close (MFILE)
jpayne@69:         or print STDERR "WARNING: Trouble closing $OPT_Mfile, $!\n";
jpayne@69: 
jpayne@69:     if ( ! defined $OPT_coverage ) { $OPT_coverage = 1; }
jpayne@69: }
jpayne@69: 
jpayne@69: 
jpayne@69: #--------------------------------------------------------------- LayoutIDs ----#
jpayne@69: # For each reference and query sequence, find the set of alignments that
jpayne@69: # produce the heaviest (both in non-redundant coverage and percent
jpayne@69: # identity) alignment subset of each sequence using a modified version
jpayne@69: # of the longest increasing subset algorithm. Let R be the union of all
jpayne@69: # reference LIS subsets, and Q be the union of all query LIS
jpayne@69: # subsets. Let S be the intersection of R and Q. Using this LIS subset,
jpayne@69: # recursively span reference and query sequences by their smaller
jpayne@69: # counterparts until all spanning sequences have been placed. The goal
jpayne@69: # is to cluster all the "major" alignment information along the main
jpayne@69: # diagonal for easy viewing and interpretation.
jpayne@69: sub LayoutIDs ($$)
jpayne@69: {
jpayne@69:     my $rref = shift;
jpayne@69:     my $qref = shift;
jpayne@69: 
jpayne@69:     my %rc;          # chains of qry seqs needed to span each ref
jpayne@69:     my %qc;          # chains of ref seqs needed to span each qry
jpayne@69:     #  {idR} -> [ placed, len, {idQ} -> [ \slope, \loR, \hiR, \loQ, \hiQ ] ]
jpayne@69:     #  {idQ} -> [ placed, len, {idR} -> [ \slope, \loQ, \hiQ, \loR, \hiR ] ]
jpayne@69: 
jpayne@69:     my @rl;          # oo of ref seqs
jpayne@69:     my @ql;          # oo of qry seqs
jpayne@69:     #  [ [idR, slope] ]
jpayne@69:     #  [ [idQ, slope] ]
jpayne@69: 
jpayne@69:     #-- get the filtered alignments
jpayne@69:     open (BTAB, "$BIN_DIR/show-coords -B $OPT_Dfile |")
jpayne@69:         or die "ERROR: Could not open show-coords pipe, $!\n";
jpayne@69: 
jpayne@69:     my @align;
jpayne@69:     my ($sR, $eR, $sQ, $eQ, $lenR, $lenQ, $idR, $idQ);
jpayne@69:     my ($loR, $hiR, $loQ, $hiQ);
jpayne@69:     my ($dR, $dQ, $slope);
jpayne@69:     while ( <BTAB> ) {
jpayne@69:         chomp;
jpayne@69:         @align = split "\t";
jpayne@69:         if ( scalar @align != 21 ) {
jpayne@69:             die "ERROR: Could not read show-coords pipe, invalid btab format\n";
jpayne@69:         }
jpayne@69: 
jpayne@69:         $sR   = $align[8];   $eR   = $align[9];
jpayne@69:         $sQ   = $align[6];   $eQ   = $align[7];
jpayne@69:         $lenR = $align[18];  $lenQ = $align[2];
jpayne@69:         $idR  = $align[5];   $idQ  = $align[0];
jpayne@69: 
jpayne@69:         #-- skip it if not on include list
jpayne@69:         if ( !exists $rref->{$idR} || !exists $qref->{$idQ} ) { next; }
jpayne@69: 
jpayne@69:         #-- get orientation of both alignments and alignment slope
jpayne@69:         $dR = $sR < $eR ? 1 : -1;
jpayne@69:         $dQ = $sQ < $eQ ? 1 : -1;
jpayne@69:         $slope = $dR == $dQ ? 1 : -1;
jpayne@69: 
jpayne@69:         #-- get lo's and hi's
jpayne@69:         $loR = $dR == 1 ? $sR : $eR;
jpayne@69:         $hiR = $dR == 1 ? $eR : $sR;
jpayne@69: 
jpayne@69:         $loQ = $dQ == 1 ? $sQ : $eQ;
jpayne@69:         $hiQ = $dQ == 1 ? $eQ : $sQ;
jpayne@69: 
jpayne@69:         if ($OPT_ONLY_USE_FATTEST)
jpayne@69:         {
jpayne@69:           #-- Check to see if there is another better alignment
jpayne@69:           if (exists $qc{$idQ})
jpayne@69:           {
jpayne@69:             my ($oldR) = keys %{$qc{$idQ}[2]};
jpayne@69:             my $val = $qc{$idQ}[2]{$oldR};
jpayne@69: 
jpayne@69:             if (${$val->[4]} - ${$val->[3]} > $hiR - $loR)
jpayne@69:             {
jpayne@69:               #-- Old alignment is better, skip this one
jpayne@69:               next;
jpayne@69:             }
jpayne@69:             else
jpayne@69:             {
jpayne@69:               #-- This alignment is better, prune old alignment
jpayne@69:               delete $rc{$oldR}[2]{$idQ};
jpayne@69:               delete $qc{$idQ};
jpayne@69:             }
jpayne@69:           }
jpayne@69:         }
jpayne@69: 
jpayne@69:         #-- initialize
jpayne@69:         if ( !exists $rc{$idR} ) { $rc{$idR} = [ 0, $lenR, { } ]; }
jpayne@69:         if ( !exists $qc{$idQ} ) { $qc{$idQ} = [ 0, $lenQ, { } ]; }
jpayne@69: 
jpayne@69:         #-- if no alignments for these two exist OR
jpayne@69:         #-- this alignment is bigger than the current
jpayne@69:         if ( !exists $rc{$idR}[2]{$idQ} || !exists $qc{$idQ}[2]{$idR} ||
jpayne@69:              $hiR - $loR >
jpayne@69:              ${$rc{$idR}[2]{$idQ}[2]} - ${$rc{$idR}[2]{$idQ}[1]} ) {
jpayne@69: 
jpayne@69:             #-- rc and qc reference these anonymous values
jpayne@69:             my $aref = [ $slope, $loR, $hiR, $loQ, $hiQ ];
jpayne@69: 
jpayne@69:             #-- rc is ordered [ slope, loR, hiR, loQ, hiQ ]
jpayne@69:             #-- qc is ordered [ slope, loQ, hiQ, loR, hiR ]
jpayne@69:             $rc{$idR}[2]{$idQ}[0] = $qc{$idQ}[2]{$idR}[0] = \$aref->[0];
jpayne@69:             $rc{$idR}[2]{$idQ}[1] = $qc{$idQ}[2]{$idR}[3] = \$aref->[1];
jpayne@69:             $rc{$idR}[2]{$idQ}[2] = $qc{$idQ}[2]{$idR}[4] = \$aref->[2];
jpayne@69:             $rc{$idR}[2]{$idQ}[3] = $qc{$idQ}[2]{$idR}[1] = \$aref->[3];
jpayne@69:             $rc{$idR}[2]{$idQ}[4] = $qc{$idQ}[2]{$idR}[2] = \$aref->[4];
jpayne@69:         }
jpayne@69:     }
jpayne@69: 
jpayne@69:     close (BTAB)
jpayne@69:         or print STDERR "WARNING: Trouble closing show-coords pipe, $!\n";
jpayne@69: 
jpayne@69:     #-- recursively span sequences to generate the layout
jpayne@69:     foreach $idR ( sort { $rc{$b}[1] <=> $rc{$a}[1] } keys %rc ) {
jpayne@69:         SpanXwY ($idR, \%rc, \@rl, \%qc, \@ql);
jpayne@69:     }
jpayne@69: 
jpayne@69:     #-- undefine the current offsets
jpayne@69:     foreach $idR ( keys %{$rref} ) { undef $rref->{$idR}[0]; }
jpayne@69:     foreach $idQ ( keys %{$qref} ) { undef $qref->{$idQ}[0]; }
jpayne@69: 
jpayne@69:     #-- redefine the offsets according to the new layout
jpayne@69:     my $roff = 0;
jpayne@69:     foreach my $r ( @rl ) {
jpayne@69:         $idR = $r->[0];
jpayne@69:         $rref->{$idR}[0] = $roff;
jpayne@69:         $rref->{$idR}[2] = $r->[1];
jpayne@69:         $roff += $rref->{$idR}[1] - 1;
jpayne@69:     }
jpayne@69:     #-- append the guys left out of the layout
jpayne@69:     foreach $idR ( keys %{$rref} ) {
jpayne@69:         if ( !defined $rref->{$idR}[0] ) {
jpayne@69:             $rref->{$idR}[0] = $roff;
jpayne@69:             $roff += $rref->{$idR}[1] - 1;
jpayne@69:         }
jpayne@69:     }
jpayne@69: 
jpayne@69:     #-- redefine the offsets according to the new layout
jpayne@69:     my $qoff = 0;
jpayne@69:     foreach my $q ( @ql ) {
jpayne@69:         $idQ = $q->[0];
jpayne@69:         $qref->{$idQ}[0] = $qoff;
jpayne@69:         $qref->{$idQ}[2] = $q->[1];
jpayne@69:         $qoff += $qref->{$idQ}[1] - 1;
jpayne@69:     }
jpayne@69:     #-- append the guys left out of the layout
jpayne@69:     foreach $idQ ( keys %{$qref} ) {
jpayne@69:         if ( !defined $qref->{$idQ}[0] ) {
jpayne@69:             $qref->{$idQ}[0] = $qoff;
jpayne@69:             $qoff += $qref->{$idQ}[1] - 1;
jpayne@69:         }
jpayne@69:     }
jpayne@69: }
jpayne@69: 
jpayne@69: 
jpayne@69: #----------------------------------------------------------------- SpanXwY ----#
jpayne@69: sub SpanXwY ($$$$$) {
jpayne@69:     my $x   = shift;   # idX
jpayne@69:     my $xcr = shift;   # xc ref
jpayne@69:     my $xlr = shift;   # xl ref
jpayne@69:     my $ycr = shift;   # yc ref
jpayne@69:     my $ylr = shift;   # yl ref
jpayne@69: 
jpayne@69:     my @post;
jpayne@69:     foreach my $y ( sort { ${$xcr->{$x}[2]{$a}[1]} <=> ${$xcr->{$x}[2]{$b}[1]} }
jpayne@69:                     keys %{$xcr->{$x}[2]} ) {
jpayne@69: 
jpayne@69:         #-- skip if already placed (RECURSION BASE)
jpayne@69:         if ( $ycr->{$y}[0] ) { next; }
jpayne@69:         else { $ycr->{$y}[0] = 1; }
jpayne@69: 
jpayne@69:         #-- get len and slope info for y
jpayne@69:         my $len = $ycr->{$y}[1];
jpayne@69:         my $slope = ${$xcr->{$x}[2]{$y}[0]};
jpayne@69: 
jpayne@69:         #-- if we need to flip, reverse complement all y records
jpayne@69:         if ( $slope == -1 ) {
jpayne@69:             foreach my $xx ( keys %{$ycr->{$y}[2]} ) {
jpayne@69:                 ${$ycr->{$y}[2]{$xx}[0]} *= -1;
jpayne@69: 
jpayne@69:                 my $loy = ${$ycr->{$y}[2]{$xx}[1]};
jpayne@69:                 my $hiy = ${$ycr->{$y}[2]{$xx}[2]};
jpayne@69:                 ${$ycr->{$y}[2]{$xx}[1]} = $len - $hiy + 1;
jpayne@69:                 ${$ycr->{$y}[2]{$xx}[2]} = $len - $loy + 1;
jpayne@69:             }
jpayne@69:         }
jpayne@69: 
jpayne@69:         #-- place y
jpayne@69:         push @{$ylr}, [ $y, $slope ];
jpayne@69: 
jpayne@69:         #-- RECURSE if y > x, else save for later
jpayne@69:         if ( $len > $xcr->{$x}[1] ) { SpanXwY ($y, $ycr, $ylr, $xcr, $xlr); }
jpayne@69:         else { push @post, $y; }
jpayne@69:     }
jpayne@69: 
jpayne@69:     #-- RECURSE for all y < x
jpayne@69:     foreach my $y ( @post ) { SpanXwY ($y, $ycr, $ylr, $xcr, $xlr); }
jpayne@69: }
jpayne@69: 
jpayne@69: 
jpayne@69: #---------------------------------------------------------------- PlotData ----#
jpayne@69: sub PlotData ($$$)
jpayne@69: {
jpayne@69:     my $aref = shift;
jpayne@69:     my $rref = shift;
jpayne@69:     my $qref = shift;
jpayne@69: 
jpayne@69:     print STDERR "Writing plot files $OPT_Ffile, $OPT_Rfile",
jpayne@69:     (defined $OPT_Hfile ? ", $OPT_Hfile\n" : "\n");
jpayne@69: 
jpayne@69:     open (FFILE, ">$OPT_Ffile")
jpayne@69:         or die "ERROR: Could not open $OPT_Ffile, $!\n";
jpayne@69:     print FFILE "#-- forward hits sorted by %sim\n0 0 0\n0 0 0\n\n\n";
jpayne@69: 
jpayne@69:     open (RFILE, ">$OPT_Rfile")
jpayne@69:         or die "ERROR: Could not open $OPT_Rfile, $!\n";
jpayne@69:     print RFILE "#-- reverse hits sorted by %sim\n0 0 0\n0 0 0\n\n\n";
jpayne@69: 
jpayne@69:     if ( defined $OPT_Hfile ) {
jpayne@69:         open (HFILE, ">$OPT_Hfile")
jpayne@69:             or die "ERROR: Could not open $OPT_Hfile, $!\n";
jpayne@69:         print HFILE "#-- highlighted hits sorted by %sim\n0 0 0\n0 0 0\n\n\n";
jpayne@69:     }
jpayne@69: 
jpayne@69:     my $fh;
jpayne@69:     my $align;
jpayne@69:     my $isplotted;
jpayne@69:     my $ismultiref;
jpayne@69:     my $ismultiqry;
jpayne@69:     my ($plenR, $plenQ, $pidR, $pidQ);
jpayne@69: 
jpayne@69:     #-- for each alignment sorted by ascending identity
jpayne@69:     foreach $align ( sort { $a->[4] <=> $b->[4] } @$aref ) {
jpayne@69: 
jpayne@69:         my ($sR, $eR, $sQ, $eQ, $sim, $lenR, $lenQ, $idR, $idQ) = @$align;
jpayne@69: 
jpayne@69:         if ( ! defined $pidR ) {
jpayne@69:             ($plenR, $plenQ, $pidR, $pidQ) = ($lenR, $lenQ, $idR, $idQ);
jpayne@69:         }
jpayne@69: 
jpayne@69:         #-- set the sequence offset, length, direction, etc...
jpayne@69:         my ($refoff, $reflen, $refdir);
jpayne@69:         my ($qryoff, $qrylen, $qrydir);
jpayne@69: 
jpayne@69:         if ( (%$rref) ) {
jpayne@69:             #-- skip reference sequence or set atts from hash
jpayne@69:             if ( !exists ($rref->{$idR}) ) { next; }
jpayne@69:             else { ($refoff, $reflen, $refdir) = @{$rref->{$idR}}; }
jpayne@69:         }
jpayne@69:         else {
jpayne@69:             #-- no reference hash, so default atts
jpayne@69:             ($refoff, $reflen, $refdir) = (0, $lenR, 1);
jpayne@69:         }
jpayne@69: 
jpayne@69:         if ( (%$qref) ) {
jpayne@69:             #-- skip query sequence or set atts from hash
jpayne@69:             if ( !exists ($qref->{$idQ}) ) { next; }
jpayne@69:             else { ($qryoff, $qrylen, $qrydir) = @{$qref->{$idQ}}; }
jpayne@69:         }
jpayne@69:         else {
jpayne@69:             #-- no query hash, so default atts
jpayne@69:             ($qryoff, $qrylen, $qrydir) = (0, $lenQ, 1);
jpayne@69:         }
jpayne@69: 
jpayne@69:         #-- get the orientation right
jpayne@69:         if ( $refdir == -1 ) {
jpayne@69:             $sR = $reflen - $sR + 1;
jpayne@69:             $eR = $reflen - $eR + 1;
jpayne@69:         }
jpayne@69:         if ( $qrydir == -1 ) {
jpayne@69:             $sQ = $qrylen - $sQ + 1;
jpayne@69:             $eQ = $qrylen - $eQ + 1;
jpayne@69:         }
jpayne@69: 
jpayne@69:         #-- forward file, reverse file, highlight file
jpayne@69:         my @fha;
jpayne@69: 
jpayne@69:         if ( defined $OPT_breaklen &&
jpayne@69:              ( ($sR - 1 > $OPT_breaklen &&
jpayne@69:                 $sQ - 1 > $OPT_breaklen &&
jpayne@69:                 $reflen - $sR > $OPT_breaklen &&
jpayne@69:                 $qrylen - $sQ > $OPT_breaklen)
jpayne@69:                ||
jpayne@69:                ($eR - 1 > $OPT_breaklen &&
jpayne@69:                 $eQ - 1 > $OPT_breaklen &&
jpayne@69:                 $reflen - $eR > $OPT_breaklen &&
jpayne@69:                 $qrylen - $eQ > $OPT_breaklen) ) ) {
jpayne@69:             push @fha, \*HFILE;
jpayne@69:         }
jpayne@69: 
jpayne@69:         push @fha, (($sR < $eR) == ($sQ < $eQ) ? \*FFILE : \*RFILE);
jpayne@69: 
jpayne@69:         #-- plot it
jpayne@69:         $sR += $refoff; $eR += $refoff;
jpayne@69:         $sQ += $qryoff; $eQ += $qryoff;
jpayne@69: 
jpayne@69:         if ( $OPT_coverage ) {
jpayne@69:             foreach $fh ( @fha ) {
jpayne@69:                 print $fh
jpayne@69:                     "$sR 10 $sim\n", "$eR 10 $sim\n\n\n",
jpayne@69:                     "$sR $sim 0\n", "$eR $sim 0\n\n\n";
jpayne@69:             }
jpayne@69:         }
jpayne@69:         else {
jpayne@69:             foreach $fh ( @fha ) {
jpayne@69:                 print $fh "$sR $sQ $sim\n", "$eR $eQ $sim\n\n\n";
jpayne@69:             }
jpayne@69:         }            
jpayne@69: 
jpayne@69:         #-- set some flags
jpayne@69:         if ( !$ismultiref && $idR ne $pidR ) { $ismultiref = 1; }
jpayne@69:         if ( !$ismultiqry && $idQ ne $pidQ ) { $ismultiqry = 1; }
jpayne@69:         if ( !$isplotted ) { $isplotted = 1; }
jpayne@69:     }
jpayne@69: 
jpayne@69: 
jpayne@69:     #-- highlight the SNPs
jpayne@69:     if ( defined $OPT_SNP ) {
jpayne@69: 
jpayne@69:         print STDERR "Determining SNPs from sequence and alignment data\n";
jpayne@69: 
jpayne@69:         open (SNPS, "$BIN_DIR/show-snps -H -T -l $OPT_Mfile |")
jpayne@69:             or die "ERROR: Could not open show-snps pipe, $!\n";
jpayne@69: 
jpayne@69:         my @snps;
jpayne@69:         my ($pR, $pQ, $lenR, $lenQ, $idR, $idQ);
jpayne@69:         while ( <SNPS> ) {
jpayne@69:             chomp;
jpayne@69:             @snps = split "\t";
jpayne@69:             if ( scalar @snps != 14 ) {
jpayne@69:                 die "ERROR: Could not read show-snps pipe, invalid format\n";
jpayne@69:             }
jpayne@69: 
jpayne@69:             $pR   = $snps[0];   $pQ   = $snps[3];
jpayne@69:             $lenR = $snps[8];   $lenQ = $snps[9];
jpayne@69:             $idR  = $snps[12];  $idQ  = $snps[13];
jpayne@69: 
jpayne@69:             #-- set the sequence offset, length, direction, etc...
jpayne@69:             my ($refoff, $reflen, $refdir);
jpayne@69:             my ($qryoff, $qrylen, $qrydir);
jpayne@69:             
jpayne@69:             if ( (%$rref) ) {
jpayne@69:                 #-- skip reference sequence or set atts from hash
jpayne@69:                 if ( !exists ($rref->{$idR}) ) { next; }
jpayne@69:                 else { ($refoff, $reflen, $refdir) = @{$rref->{$idR}}; }
jpayne@69:             }
jpayne@69:             else {
jpayne@69:                 #-- no reference hash, so default atts
jpayne@69:                 ($refoff, $reflen, $refdir) = (0, $lenR, 1);
jpayne@69:             }
jpayne@69:             
jpayne@69:             if ( (%$qref) ) {
jpayne@69:                 #-- skip query sequence or set atts from hash
jpayne@69:                 if ( !exists ($qref->{$idQ}) ) { next; }
jpayne@69:                 else { ($qryoff, $qrylen, $qrydir) = @{$qref->{$idQ}}; }
jpayne@69:             }
jpayne@69:             else {
jpayne@69:                 #-- no query hash, so default atts
jpayne@69:                 ($qryoff, $qrylen, $qrydir) = (0, $lenQ, 1);
jpayne@69:             }
jpayne@69: 
jpayne@69:             #-- get the orientation right
jpayne@69:             if ( $refdir == -1 ) { $pR = $reflen - $pR + 1; }
jpayne@69:             if ( $qrydir == -1 ) { $pQ = $qrylen - $pQ + 1; }
jpayne@69: 
jpayne@69:             #-- plot it
jpayne@69:             $pR += $refoff;
jpayne@69:             $pQ += $qryoff;
jpayne@69: 
jpayne@69:             if ( $OPT_coverage ) {
jpayne@69:                 print HFILE "$pR 10 0\n", "$pR 10 0\n\n\n",
jpayne@69:             }
jpayne@69:             else {
jpayne@69:                 print HFILE "$pR $pQ 0\n", "$pR $pQ 0\n\n\n";
jpayne@69:             }            
jpayne@69:         }
jpayne@69: 
jpayne@69:         close (SNPS)
jpayne@69:             or print STDERR "WARNING: Trouble closing show-snps pipe, $!\n";
jpayne@69:     }
jpayne@69: 
jpayne@69: 
jpayne@69:     close (FFILE)
jpayne@69:         or print STDERR "WARNING: Trouble closing $OPT_Ffile, $!\n";
jpayne@69: 
jpayne@69:     close (RFILE)
jpayne@69:         or print STDERR "WARNING: Trouble closing $OPT_Rfile, $!\n";
jpayne@69: 
jpayne@69:     if ( defined $OPT_Hfile ) {
jpayne@69:         close (HFILE)
jpayne@69:             or print STDERR "WARNING: Trouble closing $OPT_Hfile, $!\n";
jpayne@69:     }
jpayne@69: 
jpayne@69: 
jpayne@69:     if ( !(%$rref) ) {
jpayne@69:         if ( $ismultiref ) {
jpayne@69:             print STDERR
jpayne@69:                 "WARNING: Multiple ref sequences overlaid, try -R or -r\n";
jpayne@69:         }
jpayne@69:         elsif ( defined $pidR ) {
jpayne@69:             $rref->{$pidR} = [ 0, $plenR, 1 ];
jpayne@69:         }
jpayne@69:     }
jpayne@69: 
jpayne@69:     if ( !(%$qref) ) {
jpayne@69:         if ( $ismultiqry && !$OPT_coverage ) {
jpayne@69:             print STDERR
jpayne@69:                 "WARNING: Multiple qry sequences overlaid, try -Q, -q or -c\n";
jpayne@69:         }
jpayne@69:         elsif ( defined $pidQ ) {
jpayne@69:             $qref->{$pidQ} = [ 0, $plenQ, 1 ];
jpayne@69:         }
jpayne@69:     }
jpayne@69: 
jpayne@69:     if ( !$isplotted ) {
jpayne@69:         die "ERROR: No alignment data to plot\n";
jpayne@69:     }
jpayne@69: }
jpayne@69: 
jpayne@69: 
jpayne@69: #----------------------------------------------------------------- WriteGP ----#
jpayne@69: sub WriteGP ($$)
jpayne@69: {
jpayne@69:     my $rref = shift;
jpayne@69:     my $qref = shift;
jpayne@69: 
jpayne@69:     print STDERR "Writing gnuplot script $OPT_Gfile\n";
jpayne@69: 
jpayne@69:     open (GFILE, ">$OPT_Gfile")
jpayne@69:         or die "ERROR: Could not open $OPT_Gfile, $!\n";
jpayne@69: 
jpayne@69:     my ($FWD, $REV, $HLT) = (1, 2, 3);
jpayne@69:     my $SIZE = $TERMSIZE{$OPT_terminal}{$OPT_size};
jpayne@69: 
jpayne@69:     #-- terminal specific stuff
jpayne@69:     my ($P_TERM, $P_SIZE, %P_PS, %P_LW);
jpayne@69:     foreach ( $OPT_terminal ) {
jpayne@69:         /^$X11/    and do {
jpayne@69:             $P_TERM = $OPT_gpstatus == 0 ?
jpayne@69:                 "$X11 font \"$FFACE,$FSIZE\"" : "$X11";
jpayne@69: 
jpayne@69:             %P_PS = ( $FWD => 1.0, $REV => 1.0, $HLT => 1.0 );
jpayne@69: 
jpayne@69:             %P_LW = $OPT_coverage || $OPT_color ?
jpayne@69:                 ( $FWD => 3.0, $REV => 3.0, $HLT => 3.0 ) :
jpayne@69:                 ( $FWD => 2.0, $REV => 2.0, $HLT => 2.0 );
jpayne@69: 
jpayne@69:             $P_SIZE = $OPT_coverage ?
jpayne@69:                 "set size 1,1" :
jpayne@69:                 "set size 1,1";
jpayne@69: 
jpayne@69:             last;
jpayne@69:         };
jpayne@69: 
jpayne@69:         /^$PS/     and do {
jpayne@69:             $P_TERM = defined $OPT_color && $OPT_color == 0 ?
jpayne@69:                 "$PS monochrome" : "$PS color";
jpayne@69:             $P_TERM .= $OPT_gpstatus == 0 ?
jpayne@69:                 " solid \"$FFACE\" $FSIZE" : " solid \"$FFACE\" $FSIZE";
jpayne@69: 
jpayne@69:             %P_PS = ( $FWD => 0.5, $REV => 0.5, $HLT => 0.5 );
jpayne@69: 
jpayne@69:             %P_LW = $OPT_coverage || $OPT_color ?
jpayne@69:                 ( $FWD => 4.0, $REV => 4.0, $HLT => 4.0 ) :
jpayne@69:                 ( $FWD => 2.0, $REV => 2.0, $HLT => 2.0 );
jpayne@69: 
jpayne@69:             $P_SIZE = $OPT_coverage ?
jpayne@69:                 "set size ".(1.0 * $SIZE).",".(0.5 * $SIZE) :
jpayne@69:                 "set size ".(1.0 * $SIZE).",".(1.0 * $SIZE);
jpayne@69: 
jpayne@69:             last;
jpayne@69:         };
jpayne@69: 
jpayne@69:         /^$PNG/    and do {
jpayne@69:             $P_TERM = $OPT_gpstatus == 0 ?
jpayne@69:                 "$PNG tiny size $SIZE,$SIZE" : "$PNG small";
jpayne@69:             if ( defined $OPT_color && $OPT_color == 0 ) {
jpayne@69:                 $P_TERM .= " xffffff x000000 x000000";
jpayne@69:                 $P_TERM .= " x000000 x000000 x000000";
jpayne@69:                 $P_TERM .= " x000000 x000000 x000000";
jpayne@69:             }
jpayne@69:             
jpayne@69:             %P_PS = ( $FWD => 1.0, $REV => 1.0, $HLT => 1.0 );
jpayne@69: 
jpayne@69:             %P_LW = $OPT_coverage || $OPT_color ?
jpayne@69:                 ( $FWD => 3.0, $REV => 3.0, $HLT => 3.0 ) :
jpayne@69:                 ( $FWD => 3.0, $REV => 3.0, $HLT => 3.0 );
jpayne@69: 
jpayne@69:             $P_SIZE = $OPT_coverage ?
jpayne@69:                 "set size 1,.375" :
jpayne@69:                 "set size 1,1";
jpayne@69: 
jpayne@69:             last;
jpayne@69:         };
jpayne@69: 
jpayne@69:         die "ERROR: Don't know how to initialize terminal, $OPT_terminal\n";
jpayne@69:     }
jpayne@69: 
jpayne@69:     #-- plot commands
jpayne@69:     my ($P_WITH, $P_FORMAT, $P_LS, $P_KEY, %P_PT, %P_LT);
jpayne@69: 
jpayne@69:     %P_PT = ( $FWD => 6, $REV => 6, $HLT => 6 );
jpayne@69:     %P_LT = defined $OPT_Hfile ?
jpayne@69:         ( $FWD => 2, $REV => 2, $HLT => 1 ) :
jpayne@69:         ( $FWD => 1, $REV => 3, $HLT => 2 );
jpayne@69: 
jpayne@69:     $P_WITH = $OPT_coverage || $OPT_color ? "w l" : "w lp";
jpayne@69: 
jpayne@69:     $P_FORMAT = "set format \"$TFORMAT\"";
jpayne@69:     if ( $OPT_gpstatus == 0 ) {
jpayne@69:         $P_LS = "set style line";
jpayne@69:         $P_KEY = "unset key";
jpayne@69:         $P_FORMAT .= "\nset mouse format \"$TFORMAT\"";
jpayne@69:         $P_FORMAT .= "\nset mouse mouseformat \"$MFORMAT\"";
jpayne@69:         $P_FORMAT .= "\nif(GPVAL_VERSION < 5) set mouse clipboardformat \"$MFORMAT\"";
jpayne@69:     }
jpayne@69:     else {
jpayne@69:         $P_LS = "set linestyle";
jpayne@69:         $P_KEY = "set nokey";
jpayne@69:     }
jpayne@69: 
jpayne@69: 
jpayne@69:     my @refk = keys (%$rref);
jpayne@69:     my @qryk = keys (%$qref);
jpayne@69:     my ($xrange, $yrange);
jpayne@69:     my ($xlabel, $ylabel);
jpayne@69:     my ($tic, $dir);
jpayne@69:     my $border = 0;
jpayne@69: 
jpayne@69:     #-- terminal header and output
jpayne@69:     print GFILE "set terminal $P_TERM\n";
jpayne@69: 
jpayne@69:     if ( defined $OPT_Pfile ) {
jpayne@69:         print GFILE "set output \"$OPT_Pfile\"\n";
jpayne@69:     }
jpayne@69: 
jpayne@69:     if ( defined $OPT_title ) {
jpayne@69:         print GFILE "set title \"$OPT_title\"\n";
jpayne@69:     }
jpayne@69: 
jpayne@69:     #-- set tics, determine labels, ranges (ref)
jpayne@69:     if ( scalar (@refk) == 1 ) {
jpayne@69:         $xlabel = $refk[0];
jpayne@69:         $xrange = $rref->{$xlabel}[1];
jpayne@69:     }
jpayne@69:     else {
jpayne@69:         $xrange = 0;
jpayne@69:         print GFILE "set xtics rotate \( \\\n";
jpayne@69:         foreach $xlabel ( sort { $rref->{$a}[0] <=> $rref->{$b}[0] } @refk ) {
jpayne@69:             $xrange += $rref->{$xlabel}[1];
jpayne@69:             $tic = $rref->{$xlabel}[0] + 1;
jpayne@69:             $dir = ($rref->{$xlabel}[2] == 1) ? "" : "*";
jpayne@69:             print GFILE " \"$dir$xlabel\" $tic, \\\n";
jpayne@69:         }
jpayne@69:         print GFILE " \"\" $xrange \\\n\)\n";
jpayne@69:         $xlabel = "REF";
jpayne@69:     }
jpayne@69:     if ( $xrange == 0 ) { $xrange = "*"; }
jpayne@69: 
jpayne@69:     #-- set tics, determine labels, ranges (qry)
jpayne@69:     if ( $OPT_coverage ) {
jpayne@69:         $ylabel = "%SIM";
jpayne@69:         $yrange = 110;
jpayne@69:     }
jpayne@69:     elsif ( scalar (@qryk) == 1 ) {
jpayne@69:         $ylabel = $qryk[0];
jpayne@69:         $yrange = $qref->{$ylabel}[1];
jpayne@69:     }
jpayne@69:     else {
jpayne@69:         $yrange = 0;
jpayne@69:         print GFILE "set ytics \( \\\n";
jpayne@69:         foreach $ylabel ( sort { $qref->{$a}[0] <=> $qref->{$b}[0] } @qryk ) {
jpayne@69:             $yrange += $qref->{$ylabel}[1];
jpayne@69:             $tic = $qref->{$ylabel}[0] + 1;
jpayne@69:             $dir = ($qref->{$ylabel}[2] == 1) ? "" : "*";
jpayne@69:             print GFILE " \"$dir$ylabel\" $tic, \\\n";
jpayne@69:         }
jpayne@69:         print GFILE " \"\" $yrange \\\n\)\n";
jpayne@69:         $ylabel = "QRY";
jpayne@69:     }
jpayne@69:     if ( $yrange == 0 ) { $yrange = "*"; }
jpayne@69: 
jpayne@69:     #-- determine borders
jpayne@69:     if ( $xrange ne "*" && scalar (@refk) == 1 ) { $border |= 10; }
jpayne@69:     if ( $yrange ne "*" && scalar (@qryk) == 1 ) { $border |= 5; }
jpayne@69:     if ( $OPT_coverage ) { $border |= 5; }
jpayne@69: 
jpayne@69:     #-- grid, labels, border
jpayne@69:     print GFILE
jpayne@69:         "$P_SIZE\n",
jpayne@69:         "set grid\n",
jpayne@69:         "$P_KEY\n",
jpayne@69:         "set border $border\n",
jpayne@69:         "set tics scale 0\n",
jpayne@69:         "set xlabel \"$xlabel\"\n",
jpayne@69:         "set ylabel \"$ylabel\"\n",
jpayne@69:         "$P_FORMAT\n";
jpayne@69: 
jpayne@69:     #-- ranges
jpayne@69:     if ( defined $OPT_xrange ) { print GFILE "set xrange $OPT_xrange\n"; }
jpayne@69:     else                       { print GFILE "set xrange [1:$xrange]\n"; }
jpayne@69: 
jpayne@69:     if ( defined $OPT_yrange ) { print GFILE "set yrange $OPT_yrange\n"; }
jpayne@69:     else                       { print GFILE "set yrange [1:$yrange]\n"; }
jpayne@69: 
jpayne@69:     #-- if %sim plot
jpayne@69:     if ( $OPT_color ) {
jpayne@69:         print GFILE
jpayne@69:             "set zrange [0:100]\n",
jpayne@69:             "set colorbox default\n",
jpayne@69:             "set cblabel \"%similarity\"\n",
jpayne@69:             "set cbrange [0:100]\n",
jpayne@69:             "set cbtics 20\n",
jpayne@69:             "set pm3d map\n",
jpayne@69:             "set palette model RGB defined ( \\\n",
jpayne@69:             "  0 \"#000000\", \\\n",
jpayne@69:             "  4 \"#DD00DD\", \\\n",
jpayne@69:             "  6 \"#0000DD\", \\\n",
jpayne@69:             "  7 \"#00DDDD\", \\\n",
jpayne@69:             "  8 \"#00DD00\", \\\n",
jpayne@69:             "  9 \"#DDDD00\", \\\n",
jpayne@69:             " 10 \"#DD0000\"  \\\n)\n";
jpayne@69:     }
jpayne@69: 
jpayne@69:     foreach my $s ( ($FWD, $REV, $HLT) ) {
jpayne@69:         my $ss = "$P_LS $s ";
jpayne@69:         $ss .= $OPT_color ? " palette" : " lt $P_LT{$s}";
jpayne@69:         $ss .= " lw $P_LW{$s}";
jpayne@69:         if ( ! $OPT_coverage || $s == $HLT ) {
jpayne@69:             $ss .= " pt $P_PT{$s} ps $P_PS{$s}";
jpayne@69:         }
jpayne@69:         print GFILE "$ss\n";
jpayne@69:     }
jpayne@69: 
jpayne@69:     #-- plot it
jpayne@69:     print GFILE
jpayne@69:         ($OPT_color ? "splot \\\n" : "plot \\\n");
jpayne@69:     print GFILE
jpayne@69:         " \"$OPT_Ffile\" title \"FWD\" $P_WITH ls $FWD, \\\n",
jpayne@69:         " \"$OPT_Rfile\" title \"REV\" $P_WITH ls $REV",
jpayne@69:         (! defined $OPT_Hfile ? "\n" :
jpayne@69:          ", \\\n \"$OPT_Hfile\" title \"HLT\" w lp ls $HLT");
jpayne@69:     
jpayne@69:     #-- interactive mode
jpayne@69:     if ( $OPT_terminal eq $X11 ) {
jpayne@69:         print GFILE "\n",
jpayne@69:         "print \"-- INTERACTIVE MODE --\"\n",
jpayne@69:         "print \"consult gnuplot docs for command list\"\n",
jpayne@69:         "print \"mouse 1: coords to clipboard\"\n",
jpayne@69:         "print \"mouse 2: mark on plot\"\n",
jpayne@69:         "print \"mouse 3: zoom box\"\n",
jpayne@69:         "print \"'h' for help in plot window\"\n",
jpayne@69:         "print \"enter to exit\"\n",
jpayne@69:         "pause -1\n";
jpayne@69:     }
jpayne@69: 
jpayne@69:     close (GFILE)
jpayne@69:         or print STDERR "WARNING: Trouble closing $OPT_Gfile, $!\n";
jpayne@69: }
jpayne@69: 
jpayne@69: 
jpayne@69: #------------------------------------------------------------------- RunGP ----#
jpayne@69: sub RunGP ( )
jpayne@69: {
jpayne@69:     if ( defined $OPT_Pfile ) {
jpayne@69:         print STDERR "Rendering plot $OPT_Pfile\n";
jpayne@69:     }
jpayne@69:     else {
jpayne@69:         print STDERR "Rendering plot to screen\n";
jpayne@69:     }
jpayne@69: 
jpayne@69:     my $cmd = "gnuplot";
jpayne@69: 
jpayne@69:     #-- x11 specifics
jpayne@69:     if ( $OPT_terminal eq $X11 ) {
jpayne@69:         my $size = $TERMSIZE{$OPT_terminal}{$OPT_size};
jpayne@69:         $cmd .= " -geometry ${size}x";
jpayne@69:         if ( $OPT_coverage ) { $size = sprintf ("%.0f", $size * .375); }
jpayne@69:         $cmd .= "${size}+0+0 -title mummerplot";
jpayne@69: 
jpayne@69:         if ( defined $OPT_color && $OPT_color == 0 ) {
jpayne@69:             $cmd .= " -mono";
jpayne@69:             $cmd .= " -xrm 'gnuplot*line1Dashes: 0'";
jpayne@69:             $cmd .= " -xrm 'gnuplot*line2Dashes: 0'";
jpayne@69:             $cmd .= " -xrm 'gnuplot*line3Dashes: 0'";
jpayne@69:         }
jpayne@69: 
jpayne@69:         if ( $OPT_rv ) {
jpayne@69:             $cmd .= " -rv";
jpayne@69:             $cmd .= " -xrm 'gnuplot*background: black'";
jpayne@69:             $cmd .= " -xrm 'gnuplot*textColor: white'";
jpayne@69:             $cmd .= " -xrm 'gnuplot*borderColor: white'";
jpayne@69:             $cmd .= " -xrm 'gnuplot*axisColor: white'";
jpayne@69:         }
jpayne@69:     }
jpayne@69: 
jpayne@69:     $cmd .= " $OPT_Gfile";
jpayne@69: 
jpayne@69:     system ($cmd)
jpayne@69:         and print STDERR "WARNING: Unable to run '$cmd', $!\n";
jpayne@69: }
jpayne@69: 
jpayne@69: 
jpayne@69: #---------------------------------------------------------------- ListenGP ----#
jpayne@69: sub ListenGP($$)
jpayne@69: {
jpayne@69:     my $rref = shift;
jpayne@69:     my $qref = shift;
jpayne@69: 
jpayne@69:     my ($refc, $qryc);
jpayne@69:     my ($refid, $qryid);
jpayne@69:     my ($rsock, $qsock);
jpayne@69:     my $oldclip = "";
jpayne@69: 
jpayne@69:     #-- get IDs sorted by offset
jpayne@69:     my @refo = sort { $rref->{$a}[0] <=> $rref->{$b}[0] } keys %$rref;
jpayne@69:     my @qryo = sort { $qref->{$a}[0] <=> $qref->{$b}[0] } keys %$qref;
jpayne@69: 
jpayne@69:     #-- attempt to connect sockets
jpayne@69:     if ( $OPT_rport ) {
jpayne@69:         $rsock = IO::Socket::INET->new("localhost:$OPT_rport")
jpayne@69:             or print STDERR "WARNING: Could not connect to rport $OPT_rport\n";
jpayne@69:     }
jpayne@69: 
jpayne@69:     if ( $OPT_qport ) {
jpayne@69:         $qsock = IO::Socket::INET->new("localhost:$OPT_qport")
jpayne@69:             or print STDERR "WARNING: Could not connect to qport $OPT_qport\n";
jpayne@69:     }
jpayne@69: 
jpayne@69:     #-- while parent still exists
jpayne@69:     while ( getppid != 1 ) {
jpayne@69: 
jpayne@69:         #-- query the clipboard
jpayne@69:         $_ = `xclip -o -silent -selection primary`;
jpayne@69:         if ( $? >> 8 ) {
jpayne@69:             die "WARNING: Unable to query clipboard with xclip\n";
jpayne@69:         }
jpayne@69: 
jpayne@69:         #-- if cliboard has changed and contains a coordinate
jpayne@69:         if ( $_ ne $oldclip && (($refc, $qryc) = /^\[(\d+), (\d+)\]/) ) {
jpayne@69: 
jpayne@69:             $oldclip = $_;
jpayne@69: 
jpayne@69:             #-- translate the reference position
jpayne@69:             $refid = "NULL";
jpayne@69:             for ( my $i = 0; $i < (scalar @refo); ++ $i ) {
jpayne@69:                 my $aref = $rref->{$refo[$i]};
jpayne@69:                 if ( $i == $#refo || $aref->[0] + $aref->[1] > $refc ) {
jpayne@69:                     $refid = $refo[$i];
jpayne@69:                     $refc -= $aref->[0];
jpayne@69:                     if ( $aref->[2] == -1 ) {
jpayne@69:                         $refc = $aref->[1] - $refc + 1;
jpayne@69:                     }
jpayne@69:                     last;
jpayne@69:                 }
jpayne@69:             }
jpayne@69: 
jpayne@69:             #-- translate the query position
jpayne@69:             $qryid = "NULL";
jpayne@69:             for ( my $i = 0; $i < (scalar @qryo); ++ $i ) {
jpayne@69:                 my $aref = $qref->{$qryo[$i]};
jpayne@69:                 if ( $i == $#qryo || $aref->[0] + $aref->[1] > $qryc ) {
jpayne@69:                     $qryid = $qryo[$i];
jpayne@69:                     $qryc -= $aref->[0];
jpayne@69:                     if ( $aref->[2] == -1 ) {
jpayne@69:                         $qryc = $aref->[1] - $qryc + 1;
jpayne@69:                     }
jpayne@69:                     last;
jpayne@69:                 }
jpayne@69:             }
jpayne@69: 
jpayne@69:             #-- print the info to stdout and socket
jpayne@69:             print "$refid\t$qryid\t$refc\t$qryc\n";
jpayne@69: 
jpayne@69:             if ( $rsock ) {
jpayne@69:                 print $rsock "contig I$refid $refc\n";
jpayne@69:                 print "sent \"contig I$refid $refc\" to $OPT_rport\n";
jpayne@69:             }
jpayne@69:             if ( $qsock ) {
jpayne@69:                 print $qsock "contig I$qryid $qryc\n";
jpayne@69:                 print "sent \"contig I$qryid $qryc\" to $OPT_qport\n";
jpayne@69:             }
jpayne@69:         }
jpayne@69: 
jpayne@69:         #-- sleep for half second
jpayne@69:         select undef, undef, undef, .5;
jpayne@69:     }
jpayne@69: 
jpayne@69:     exit (0);
jpayne@69: }
jpayne@69: 
jpayne@69: 
jpayne@69: #------------------------------------------------------------ ParseOptions ----#
jpayne@69: sub ParseOptions ( )
jpayne@69: {
jpayne@69:     my ($opt_small, $opt_medium, $opt_large);
jpayne@69:     my ($opt_ps, $opt_x11, $opt_png);
jpayne@69:     my $cnt;
jpayne@69: 
jpayne@69:     #-- Get options
jpayne@69:     my $err = $tigr -> TIGR_GetOptions
jpayne@69:         (
jpayne@69:          "b|breaklen:i" => \$OPT_breaklen,
jpayne@69:          "color!"       => \$OPT_color,
jpayne@69:          "c|coverage!"  => \$OPT_coverage,
jpayne@69:          "f|filter!"    => \$OPT_filter,
jpayne@69:          "l|layout!"    => \$OPT_layout,
jpayne@69:          "p|prefix=s"   => \$OPT_prefix,
jpayne@69:          "rv"           => \$OPT_rv,
jpayne@69:          "r|IdR=s"      => \$OPT_IdR,
jpayne@69:          "q|IdQ=s"      => \$OPT_IdQ,
jpayne@69:          "R|Rfile=s"    => \$OPT_IDRfile,
jpayne@69:          "Q|Qfile=s"    => \$OPT_IDQfile,
jpayne@69:          "rport=i"      => \$OPT_rport,
jpayne@69:          "qport=i"      => \$OPT_qport,
jpayne@69:          "s|size=s"     => \$OPT_size,
jpayne@69:          "S|SNP"        => \$OPT_SNP,
jpayne@69:          "t|terminal=s" => \$OPT_terminal,
jpayne@69:          "title=s"      => \$OPT_title,
jpayne@69:          "x|xrange=s"   => \$OPT_xrange,
jpayne@69:          "y|yrange=s"   => \$OPT_yrange,
jpayne@69:          "x11"          => \$opt_x11,
jpayne@69:          "postscript"   => \$opt_ps,
jpayne@69:          "png"          => \$opt_png,
jpayne@69:          "small"        => \$opt_small,
jpayne@69:          "medium"       => \$opt_medium,
jpayne@69:          "large"        => \$opt_large,
jpayne@69:          "fat"          => \$OPT_ONLY_USE_FATTEST,
jpayne@69:          );
jpayne@69: 
jpayne@69:     if ( !$err  ||  scalar (@ARGV) != 1 ) {
jpayne@69:         $tigr -> printUsageInfo( );
jpayne@69:         die "Try '$0 -h' for more information.\n";
jpayne@69:     }
jpayne@69: 
jpayne@69:     $cnt = 0;
jpayne@69:     if ( $opt_png ) { $OPT_terminal = $PNG; $cnt ++; }
jpayne@69:     if ( $opt_ps  ) { $OPT_terminal = $PS;  $cnt ++; }
jpayne@69:     if ( $opt_x11 ) { $OPT_terminal = $X11; $cnt ++; }
jpayne@69:     if ( $cnt > 1 ) {
jpayne@69:         print STDERR
jpayne@69:             "WARNING: Multiple terminals not allowed, using '$OPT_terminal'\n";
jpayne@69:     }
jpayne@69: 
jpayne@69:     $cnt = 0;
jpayne@69:     if ( $opt_large  ) { $OPT_size = $LARGE;  $cnt ++; }
jpayne@69:     if ( $opt_medium ) { $OPT_size = $MEDIUM; $cnt ++; }
jpayne@69:     if ( $opt_small  ) { $OPT_size = $SMALL;  $cnt ++; }
jpayne@69:     if ( $cnt > 1 ) {
jpayne@69:         print STDERR
jpayne@69:             "WARNING: Multiple sizes now allowed, using '$OPT_size'\n";
jpayne@69:     }
jpayne@69: 
jpayne@69:     #-- Check that status of gnuplot
jpayne@69:     $OPT_gpstatus = system ("gnuplot --version");
jpayne@69: 
jpayne@69:     if ( $OPT_gpstatus == -1 ) {
jpayne@69:         print STDERR
jpayne@69:             "WARNING: Could not find gnuplot, plot will not be rendered\n";
jpayne@69:     }
jpayne@69:     elsif ( $OPT_gpstatus ) {
jpayne@69:         print STDERR
jpayne@69:             "WARNING: Using outdated gnuplot, use v4.0 for best results\n";
jpayne@69: 
jpayne@69:         if ( $OPT_color ) {
jpayne@69:             print STDERR
jpayne@69:                 "WARNING: Turning of --color option for compatibility\n";
jpayne@69:             undef $OPT_color;
jpayne@69:         }
jpayne@69: 
jpayne@69:         if ( $OPT_terminal eq $PNG  &&  $OPT_size ne $SMALL ) { 
jpayne@69:             print STDERR
jpayne@69:                 "WARNING: Turning of --size option for compatibility\n";
jpayne@69:             $OPT_size = $SMALL;
jpayne@69:         }
jpayne@69:     }
jpayne@69: 
jpayne@69:     #-- Check options
jpayne@69:     if ( !exists $TERMSIZE{$OPT_terminal} ) {
jpayne@69:         die "ERROR: Invalid terminal type, $OPT_terminal\n";
jpayne@69:     }
jpayne@69: 
jpayne@69:     if ( !exists $TERMSIZE{$OPT_terminal}{$OPT_size} ) {
jpayne@69:         die "ERROR: Invalid terminal size, $OPT_size\n";
jpayne@69:     }
jpayne@69: 
jpayne@69:     if ( $OPT_xrange ) {
jpayne@69:         $OPT_xrange =~ tr/,/:/;
jpayne@69:         $OPT_xrange =~ /^\[\d+:\d+\]$/
jpayne@69:             or die "ERROR: Invalid xrange format, $OPT_xrange\n";
jpayne@69:     }
jpayne@69: 
jpayne@69:     if ( $OPT_yrange ) {
jpayne@69:         $OPT_yrange =~ tr/,/:/;
jpayne@69:         $OPT_yrange =~ /^\[\d+:\d+\]$/
jpayne@69:             or die "ERROR: Invalid yrange format, $OPT_yrange\n";
jpayne@69:     }
jpayne@69: 
jpayne@69:     #-- Set file names
jpayne@69:     $OPT_Mfile = $ARGV[0];
jpayne@69:     $tigr->isReadableFile ($OPT_Mfile)
jpayne@69:         or die "ERROR: Could not read $OPT_Mfile, $!\n";
jpayne@69: 
jpayne@69:     $OPT_Ffile = $OPT_prefix . $SUFFIX{$FWDPLOT};
jpayne@69:     $tigr->isWritableFile ($OPT_Ffile) or $tigr->isCreatableFile ($OPT_Ffile)
jpayne@69:         or die "ERROR: Could not write $OPT_Ffile, $!\n";
jpayne@69: 
jpayne@69:     $OPT_Rfile = $OPT_prefix . $SUFFIX{$REVPLOT};
jpayne@69:     $tigr->isWritableFile ($OPT_Rfile) or $tigr->isCreatableFile ($OPT_Rfile)
jpayne@69:         or die "ERROR: Could not write $OPT_Rfile, $!\n";
jpayne@69: 
jpayne@69:     if ( defined $OPT_breaklen || defined $OPT_SNP ) {
jpayne@69:         $OPT_Hfile = $OPT_prefix . $SUFFIX{$HLTPLOT};
jpayne@69:         $tigr->isWritableFile($OPT_Hfile) or $tigr->isCreatableFile($OPT_Hfile)
jpayne@69:             or die "ERROR: Could not write $OPT_Hfile, $!\n";
jpayne@69:     }
jpayne@69: 
jpayne@69:     if ($OPT_ONLY_USE_FATTEST)
jpayne@69:     {
jpayne@69:       $OPT_layout = 1;
jpayne@69:     }
jpayne@69: 
jpayne@69:     if ( $OPT_filter || $OPT_layout ) {
jpayne@69:         $OPT_Dfile = $OPT_prefix . $SUFFIX{$FILTER};
jpayne@69:         $tigr->isWritableFile($OPT_Dfile) or $tigr->isCreatableFile($OPT_Dfile)
jpayne@69:             or die "ERROR: Could not write $OPT_Dfile, $!\n";
jpayne@69:     }
jpayne@69: 
jpayne@69:     $OPT_Gfile = $OPT_prefix . $SUFFIX{$GNUPLOT};
jpayne@69:     $tigr->isWritableFile ($OPT_Gfile) or $tigr->isCreatableFile ($OPT_Gfile)
jpayne@69:         or die "ERROR: Could not write $OPT_Gfile, $!\n";
jpayne@69: 
jpayne@69:     if ( exists $SUFFIX{$OPT_terminal} ) {
jpayne@69:         $OPT_Pfile = $OPT_prefix . $SUFFIX{$OPT_terminal};
jpayne@69:         $tigr->isWritableFile($OPT_Pfile) or $tigr->isCreatableFile($OPT_Pfile)
jpayne@69:             or die "ERROR: Could not write $OPT_Pfile, $!\n";
jpayne@69:     }
jpayne@69: 
jpayne@69:     if ( defined $OPT_IDRfile ) {
jpayne@69:         $tigr->isReadableFile ($OPT_IDRfile)
jpayne@69:             or die "ERROR: Could not read $OPT_IDRfile, $!\n";
jpayne@69:     }
jpayne@69: 
jpayne@69:     if ( defined $OPT_IDQfile ) {
jpayne@69:         $tigr->isReadableFile ($OPT_IDQfile)
jpayne@69:             or die "ERROR: Could not read $OPT_IDQfile, $!\n";
jpayne@69:     }
jpayne@69: 
jpayne@69:     if ( (defined $OPT_rport || defined $OPT_qport) &&
jpayne@69:          ($OPT_terminal ne $X11 || $OPT_gpstatus ) ) {
jpayne@69:         print STDERR
jpayne@69:             "WARNING: Port options available only for v4.0 X11 plots\n";
jpayne@69:         undef $OPT_rport;
jpayne@69:         undef $OPT_qport;
jpayne@69:     }
jpayne@69: 
jpayne@69: 
jpayne@69:     if ( defined $OPT_color && defined $OPT_Hfile ) {
jpayne@69:         print STDERR
jpayne@69:             "WARNING: Turning off --color option so highlighting is visible\n";
jpayne@69:         undef $OPT_color;
jpayne@69:     }
jpayne@69: }