Mercurial > repos > rliterman > csp2

annotate CSP2/CSP2_env/env-d9b9114564458d9d-741b3de822f2aaca6c6caa4325c4afce/opt/mummer-3.23/scripts/mummerplot.pl @ 69:33d812a61356

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