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