code_perl/Maasha/Biopieces.pm

   1 package Maasha::Biopieces;
   2
   3
   4 # Copyright (C) 2007-2008 Martin A. Hansen.
   5
   6 # This program is free software; you can redistribute it and/or
   7 # modify it under the terms of the GNU General Public License
   8 # as published by the Free Software Foundation; either version 2
   9 # of the License, or (at your option) any later version.
  10
  11 # This program is distributed in the hope that it will be useful,
  12 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14 # GNU General Public License for more details.
  15
  16 # You should have received a copy of the GNU General Public License
  17 # along with this program; if not, write to the Free Software
  18 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  19
  20 # http://www.gnu.org/copyleft/gpl.html
  21
  22
  23 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DESCRIPTION <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  24
  25
  26 # Routines for manipulation, parsing and emitting of human/machine readable biopieces records.
  27
  28
  29 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  30
  31
  32 use strict;
  33 use Data::Dumper;
  34 use Getopt::Long qw( :config bundling );
  35 use Time::HiRes qw( gettimeofday );
  36 use Storable qw( dclone );
  37 use Maasha::Config;
  38 use Maasha::Common;
  39 use Maasha::Fasta;
  40 use Maasha::Align;
  41 use Maasha::Matrix;
  42 use Maasha::Match;
  43 use Maasha::EMBL;
  44 use Maasha::Stockholm;
  45 use Maasha::Seq;
  46 use Maasha::Patscan;
  47 use Maasha::Plot;
  48 use Maasha::Calc;
  49 use Maasha::UCSC;
  50 use Maasha::NCBI;
  51 use Maasha::GFF;
  52 use Maasha::TwoBit;
  53 use Maasha::Solid;
  54 use Maasha::SQL;
  55
  56 use vars qw( @ISA @EXPORT_OK );
  57
  58 require Exporter;
  59
  60 @ISA = qw( Exporter );
  61
  62 @EXPORT_OK = qw(
  63     read_stream
  64     write_stream
  65     get_record
  66     put_record
  67 );
  68
  69 use constant {
  70     SEQ_NAME => 0,
  71     SEQ      => 1,
  72 };
  73
  74
  75 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> SIGNAL HANDLER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  76
  77
  78 $SIG{ '__DIE__' } = \&sig_handler;
  79 $SIG{ 'INT' }     = \&sig_handler;
  80 $SIG{ 'TERM' }    = \&sig_handler;
  81
  82
  83 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> GLOBALS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  84
  85
  86 my ( $script, $BP_TMP );
  87
  88 $script  = Maasha::Common::get_scriptname();
  89 $BP_TMP = Maasha::Common::get_tmpdir();
  90
  91
  92 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> LOG <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  93
  94
  95 my $log_fh = Maasha::Common::append_open( $ENV{ "LOG_DIR" } . "/biopieces.log" );
  96
  97 $log_fh->autoflush( 1 );
  98
  99 log( $log_fh, $script, \@ARGV );
 100
 101 close $log_fh;
 102
 103
 104 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> RUN SCRIPT <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 105
 106
 107 my $t0 = gettimeofday();
 108
 109 run_script( $script );
 110
 111 my $t1 = gettimeofday();
 112
 113 print STDERR "Program: $script" . ( " " x ( 25 - length( $script ) ) ) . sprintf( "Run time: %.4f\n", ( $t1 - $t0 ) );
 114
 115
 116 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> SUBROUTINES <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 117
 118
 119 sub log
 120 {
 121     # Martin A. Hansen, January 2008.
 122
 123     # Log messages to logfile.
 124
 125     my ( $fh,       # filehandle to logfile
 126          $script,   # script name
 127          $argv,     # reference to @ARGV
 128        ) = @_;
 129
 130     # Returns nothing.
 131
 132     my ( $time_stamp, $user );
 133
 134     $time_stamp = Maasha::Common::time_stamp();
 135
 136     $user = $ENV{ "USER" };
 137
 138     $script = "biopieces" if $script eq "-e";
 139
 140     print $fh "$time_stamp\t$user\t$script ", join( " ", @{ $argv } ), "\n";
 141 }
 142
 143
 144 sub run_script
 145 {
 146     # Martin A. Hansen, August 2007.
 147
 148     # Run a specific script.
 149
 150     my ( $script,   # script name
 151        ) = @_;
 152
 153     # Returns nothing.
 154
 155     my ( $options, $in, $out );
 156
 157     script_list_biopieces( $ENV{ 'INST_DIR'} . "/biopieces/usage/" ) if $script eq "list_biopieces";
 158
 159     script_print_usage( $ENV{ 'INST_DIR'} . "/biopieces/usage/$script" ) if -t STDIN and not @ARGV;
 160
 161     $options = get_options( $script );
 162
 163     $in  = read_stream( $options->{ "stream_in" } );
 164     $out = write_stream( $options->{ "stream_out" } );
 165
 166     if    ( $script eq "read_fasta" )               { script_read_fasta( $in, $out, $options ) }
 167     elsif ( $script eq "read_tab" )                 { script_read_tab( $in, $out, $options ) }
 168     elsif ( $script eq "read_psl" )                 { script_read_psl( $in, $out, $options ) }
 169     elsif ( $script eq "read_bed" )                 { script_read_bed( $in, $out, $options ) }
 170     elsif ( $script eq "read_blast_tab" )           { script_read_blast_tab( $in, $out, $options ) }
 171     elsif ( $script eq "read_embl" )                { script_read_embl( $in, $out, $options ) }
 172     elsif ( $script eq "read_stockholm" )           { script_read_stockholm( $in, $out, $options ) }
 173     elsif ( $script eq "read_phastcons" )           { script_read_phastcons( $in, $out, $options ) }
 174     elsif ( $script eq "read_soft" )                { script_read_soft( $in, $out, $options ) }
 175     elsif ( $script eq "read_gff" )                 { script_read_gff( $in, $out, $options ) }
 176     elsif ( $script eq "read_2bit" )                { script_read_2bit( $in, $out, $options ) }
 177     elsif ( $script eq "read_solexa" )              { script_read_solexa( $in, $out, $options ) }
 178     elsif ( $script eq "read_solid" )               { script_read_solid( $in, $out, $options ) }
 179     elsif ( $script eq "read_mysql" )               { script_read_mysql( $in, $out, $options ) }
 180     elsif ( $script eq "length_seq" )               { script_length_seq( $in, $out, $options ) }
 181     elsif ( $script eq "uppercase_seq" )            { script_uppercase_seq( $in, $out, $options ) }
 182     elsif ( $script eq "shuffle_seq" )              { script_shuffle_seq( $in, $out, $options ) }
 183     elsif ( $script eq "analyze_seq" )              { script_analyze_seq( $in, $out, $options ) }
 184     elsif ( $script eq "analyze_tags" )             { script_analyze_tags( $in, $out, $options ) }
 185     elsif ( $script eq "complexity_seq" )           { script_complexity_seq( $in, $out, $options ) }
 186     elsif ( $script eq "oligo_freq" )               { script_oligo_freq( $in, $out, $options ) }
 187     elsif ( $script eq "create_weight_matrix" )     { script_create_weight_matrix( $in, $out, $options ) }
 188     elsif ( $script eq "calc_bit_scores" )          { script_calc_bit_scores( $in, $out, $options ) }
 189     elsif ( $script eq "reverse_seq" )              { script_reverse_seq( $in, $out, $options ) }
 190     elsif ( $script eq "complement_seq" )           { script_complement_seq( $in, $out, $options ) }
 191     elsif ( $script eq "remove_indels" )            { script_remove_indels( $in, $out, $options ) }
 192     elsif ( $script eq "transliterate_seq" )        { script_transliterate_seq( $in, $out, $options ) }
 193     elsif ( $script eq "transliterate_vals" )       { script_transliterate_vals( $in, $out, $options ) }
 194     elsif ( $script eq "translate_seq" )            { script_translate_seq( $in, $out, $options ) }
 195     elsif ( $script eq "extract_seq" )              { script_extract_seq( $in, $out, $options ) }
 196     elsif ( $script eq "get_genome_seq" )           { script_get_genome_seq( $in, $out, $options ) }
 197     elsif ( $script eq "get_genome_align" )         { script_get_genome_align( $in, $out, $options ) }
 198     elsif ( $script eq "get_genome_phastcons" )     { script_get_genome_phastcons( $in, $out, $options ) }
 199     elsif ( $script eq "fold_seq" )                 { script_fold_seq( $in, $out, $options ) }
 200     elsif ( $script eq "split_seq" )                { script_split_seq( $in, $out, $options ) }
 201     elsif ( $script eq "split_bed" )                { script_split_bed( $in, $out, $options ) }
 202     elsif ( $script eq "align_seq" )                { script_align_seq( $in, $out, $options ) }
 203     elsif ( $script eq "tile_seq" )                 { script_tile_seq( $in, $out, $options ) }
 204     elsif ( $script eq "invert_align" )             { script_invert_align( $in, $out, $options ) }
 205     elsif ( $script eq "patscan_seq" )              { script_patscan_seq( $in, $out, $options ) }
 206     elsif ( $script eq "create_blast_db" )          { script_create_blast_db( $in, $out, $options ) }
 207     elsif ( $script eq "blast_seq" )                { script_blast_seq( $in, $out, $options ) }
 208     elsif ( $script eq "blat_seq" )                 { script_blat_seq( $in, $out, $options ) }
 209     elsif ( $script eq "match_seq" )                { script_match_seq( $in, $out, $options ) }
 210     elsif ( $script eq "create_vmatch_index" )      { script_create_vmatch_index( $in, $out, $options ) }
 211     elsif ( $script eq "vmatch_seq" )               { script_vmatch_seq( $in, $out, $options ) }
 212     elsif ( $script eq "write_fasta" )              { script_write_fasta( $in, $out, $options, $options ) }
 213     elsif ( $script eq "write_align" )              { script_write_align( $in, $out, $options ) }
 214     elsif ( $script eq "write_blast" )              { script_write_blast( $in, $out, $options ) }
 215     elsif ( $script eq "write_tab" )                { script_write_tab( $in, $out, $options ) }
 216     elsif ( $script eq "write_bed" )                { script_write_bed( $in, $out, $options ) }
 217     elsif ( $script eq "write_psl" )                { script_write_psl( $in, $out, $options ) }
 218     elsif ( $script eq "write_2bit" )               { script_write_2bit( $in, $out, $options, $options ) }
 219     elsif ( $script eq "write_solid" )              { script_write_solid( $in, $out, $options, $options ) }
 220     elsif ( $script eq "head_records" )             { script_head_records( $in, $out, $options ) }
 221     elsif ( $script eq "remove_keys" )              { script_remove_keys( $in, $out, $options ) }
 222     elsif ( $script eq "rename_keys" )              { script_rename_keys( $in, $out, $options ) }
 223     elsif ( $script eq "uniq_vals" )                { script_uniq_vals( $in, $out, $options ) }
 224     elsif ( $script eq "merge_vals" )               { script_merge_vals( $in, $out, $options ) }
 225     elsif ( $script eq "grab" )                     { script_grab( $in, $out, $options ) }
 226     elsif ( $script eq "compute" )                  { script_compute( $in, $out, $options ) }
 227     elsif ( $script eq "flip_tab" )                 { script_flip_tab( $in, $out, $options ) }
 228     elsif ( $script eq "add_ident" )                { script_add_ident( $in, $out, $options ) }
 229     elsif ( $script eq "count_records" )            { script_count_records( $in, $out, $options ) }
 230     elsif ( $script eq "random_records" )           { script_random_records( $in, $out, $options ) }
 231     elsif ( $script eq "sort_records" )             { script_sort_records( $in, $out, $options ) }
 232     elsif ( $script eq "count_vals" )               { script_count_vals( $in, $out, $options ) }
 233     elsif ( $script eq "plot_histogram" )           { script_plot_histogram( $in, $out, $options ) }
 234     elsif ( $script eq "plot_lendist" )             { script_plot_lendist( $in, $out, $options ) }
 235     elsif ( $script eq "plot_chrdist" )             { script_plot_chrdist( $in, $out, $options ) }
 236     elsif ( $script eq "plot_karyogram" )           { script_plot_karyogram( $in, $out, $options ) }
 237     elsif ( $script eq "plot_matches" )             { script_plot_matches( $in, $out, $options ) }
 238     elsif ( $script eq "plot_seqlogo" )             { script_plot_seqlogo( $in, $out, $options ) }
 239     elsif ( $script eq "plot_phastcons_profiles" )  { script_plot_phastcons_profiles( $in, $out, $options ) }
 240     elsif ( $script eq "analyze_bed" )              { script_analyze_bed( $in, $out, $options ) }
 241     elsif ( $script eq "analyze_vals" )             { script_analyze_vals( $in, $out, $options ) }
 242     elsif ( $script eq "length_vals" )              { script_length_vals( $in, $out, $options ) }
 243     elsif ( $script eq "sum_vals" )                 { script_sum_vals( $in, $out, $options ) }
 244     elsif ( $script eq "mean_vals" )                { script_mean_vals( $in, $out, $options ) }
 245     elsif ( $script eq "median_vals" )              { script_median_vals( $in, $out, $options ) }
 246     elsif ( $script eq "max_vals" )                 { script_max_vals( $in, $out, $options ) }
 247     elsif ( $script eq "min_vals" )                 { script_min_vals( $in, $out, $options ) }
 248     elsif ( $script eq "upload_to_ucsc" )           { script_upload_to_ucsc( $in, $out, $options ) }
 249
 250     close $in if defined $in;
 251     close $out;
 252
 253     # unset status   - missing
 254     # write log file - missing
 255 }
 256
 257
 258 sub get_options
 259 {
 260     # Martin A. Hansen, February 2008.
 261
 262     # Gets options from commandline and checks these vigerously.
 263
 264     my ( $script,   # name of script
 265        ) = @_;
 266
 267     # Returns hash
 268
 269     my ( %options, @options, $opt, @genomes );
 270
 271     if ( $script eq "read_fasta" )
 272     {
 273         @options = qw(
 274             data_in|i=s
 275             num|n=s
 276         );
 277     }
 278     elsif ( $script eq "read_tab" )
 279     {
 280         @options = qw(
 281             data_in|i=s
 282             delimit|d=s
 283             cols|c=s
 284             keys|k=s
 285             skip|s=s
 286             num|n=s
 287         );
 288     }
 289     elsif ( $script eq "read_psl" )
 290     {
 291         @options = qw(
 292             data_in|i=s
 293             num|n=s
 294         );
 295     }
 296     elsif ( $script eq "read_bed" )
 297     {
 298         @options = qw(
 299             data_in|i=s
 300             num|n=s
 301         );
 302     }
 303     elsif ( $script eq "read_blast_tab" )
 304     {
 305         @options = qw(
 306             data_in|i=s
 307             num|n=s
 308         );
 309     }
 310     elsif ( $script eq "read_embl" )
 311     {
 312         @options = qw(
 313             data_in|i=s
 314             num|n=s
 315             keys|k=s
 316             feats|f=s
 317             quals|q=s
 318         );
 319     }
 320     elsif ( $script eq "read_stockholm" )
 321     {
 322         @options = qw(
 323             data_in|i=s
 324             num|n=s
 325         );
 326     }
 327     elsif ( $script eq "read_phastcons" )
 328     {
 329         @options = qw(
 330             data_in|i=s
 331             num|n=s
 332             min|m=s
 333             dist|d=s
 334             threshold|t=f
 335             gap|g=s
 336         );
 337     }
 338     elsif ( $script eq "read_soft" )
 339     {
 340         @options = qw(
 341             data_in|i=s
 342             num|n=s
 343         );
 344     }
 345     elsif ( $script eq "read_gff" )
 346     {
 347         @options = qw(
 348             data_in|i=s
 349             num|n=s
 350         );
 351     }
 352     elsif ( $script eq "read_2bit" )
 353     {
 354         @options = qw(
 355             data_in|i=s
 356             num|n=s
 357             no_mask|N
 358         );
 359     }
 360     elsif ( $script eq "read_solexa" )
 361     {
 362         @options = qw(
 363             data_in|i=s
 364             num|n=s
 365             quality|q=s
 366         );
 367     }
 368     elsif ( $script eq "read_solid" )
 369     {
 370         @options = qw(
 371             data_in|i=s
 372             num|n=s
 373             quality|q=s
 374         );
 375     }
 376     elsif ( $script eq "read_mysql" )
 377     {
 378         @options = qw(
 379             database|d=s
 380             query|q=s
 381             user|u=s
 382             password|p=s
 383         );
 384     }
 385     elsif ( $script eq "length_seq" )
 386     {
 387         @options = qw(
 388             no_stream|x
 389             data_out|o=s
 390         );
 391     }
 392     elsif ( $script eq "oligo_freq" )
 393     {
 394         @options = qw(
 395             word_size|w=s
 396             all|a
 397         );
 398     }
 399     elsif ( $script eq "create_weight_matrix" )
 400     {
 401         @options = qw(
 402             percent|p
 403         );
 404     }
 405     elsif ( $script eq "transliterate_seq" )
 406     {
 407         @options = qw(
 408             search|s=s
 409             replace|r=s
 410             delete|d=s
 411         );
 412     }
 413     elsif ( $script eq "transliterate_vals" )
 414     {
 415         @options = qw(
 416             keys|k=s
 417             search|s=s
 418             replace|r=s
 419             delete|d=s
 420         );
 421     }
 422     elsif ( $script eq "translate_seq" )
 423     {
 424         @options = qw(
 425             frames|f=s
 426         );
 427     }
 428     elsif ( $script eq "extract_seq" )
 429     {
 430         @options = qw(
 431             beg|b=s
 432             end|e=s
 433             len|l=s
 434         );
 435     }
 436     elsif ( $script eq "get_genome_seq" )
 437     {
 438         @options = qw(
 439             genome|g=s
 440             chr|c=s
 441             beg|b=s
 442             end|e=s
 443             len|l=s
 444             flank|f=s
 445             mask|m
 446         );
 447     }
 448     elsif ( $script eq "get_genome_align" )
 449     {
 450         @options = qw(
 451             genome|g=s
 452             chr|c=s
 453             beg|b=s
 454             end|e=s
 455             len|l=s
 456             strand|s=s
 457         );
 458     }
 459     elsif ( $script eq "get_genome_phastcons" )
 460     {
 461         @options = qw(
 462             genome|g=s
 463             chr|c=s
 464             beg|b=s
 465             end|e=s
 466             len|l=s
 467             flank|f=s
 468         );
 469     }
 470     elsif ( $script eq "split_seq" )
 471     {
 472         @options = qw(
 473             word_size|w=s
 474             uniq|u
 475         );
 476     }
 477     elsif ( $script eq "split_bed" )
 478     {
 479         @options = qw(
 480             window_size|w=s
 481             step_size|s=s
 482         );
 483     }
 484     elsif ( $script eq "tile_seq" )
 485     {
 486         @options = qw(
 487             identity|i=s
 488             supress_indels|s
 489         );
 490     }
 491     elsif ( $script eq "invert_align" )
 492     {
 493         @options = qw(
 494             soft|s
 495         );
 496     }
 497     elsif ( $script eq "patscan_seq" )
 498     {
 499         @options = qw(
 500             patterns|p=s
 501             patterns_in|P=s
 502             comp|c
 503             max_hits|h=s
 504             max_misses|m=s
 505             genome|g=s
 506         );
 507     }
 508     elsif ( $script eq "create_blast_db" )
 509     {
 510         @options = qw(
 511             no_stream|x
 512             database|d=s
 513         );
 514     }
 515     elsif ( $script eq "blast_seq" )
 516     {
 517         @options = qw(
 518             database|d=s
 519             genome|g=s
 520             program|p=s
 521             e_val|e=f
 522             filter|f
 523             cpus|c=s
 524             no_filter|F
 525         );
 526     }
 527     elsif ( $script eq "blat_seq" )
 528     {
 529         @options = qw(
 530             genome|g=s
 531             tile_size|t=s
 532             step_size|s=s
 533             min_identity|m=s
 534             min_score|M=s
 535             one_off|o=s
 536             ooc|c
 537         );
 538     }
 539     elsif ( $script eq "match_seq" )
 540     {
 541         @options = qw(
 542             word_size|w=s
 543             direction|d=s
 544         );
 545     }
 546     elsif ( $script eq "create_vmatch_index" )
 547     {
 548         @options = qw(
 549             index_name|i=s
 550             prefix_length|p=s
 551             no_stream|x
 552         );
 553     }
 554     elsif ( $script eq "vmatch_seq" )
 555     {
 556         @options = qw(
 557             genome|g=s
 558             index_name|i=s
 559             count|c
 560             max_hits|m=s
 561             hamming_dist|h=s
 562             edit_dist|e=s
 563         );
 564     }
 565     elsif ( $script eq "write_fasta" )
 566     {
 567         @options = qw(
 568             wrap|w=s
 569             no_stream|x
 570             data_out|o=s
 571             compress|Z
 572         );
 573     }
 574     elsif ( $script eq "write_align" )
 575     {
 576         @options = qw(
 577             wrap|w=s
 578             no_stream|x
 579             no_ruler|R
 580             no_consensus|C
 581             data_out|o=s
 582         );
 583     }
 584     elsif ( $script eq "write_blast" )
 585     {
 586         @options = qw(
 587             no_stream|x
 588             data_out|o=s
 589             comment|c
 590             compress|Z
 591         );
 592     }
 593     elsif ( $script eq "write_tab" )
 594     {
 595         @options = qw(
 596             no_stream|x
 597             data_out|o=s
 598             delimit|d=s
 599             keys|k=s
 600             no_keys|K=s
 601             comment|c
 602             compress|Z
 603         );
 604     }
 605     elsif ( $script eq "write_bed" )
 606     {
 607         @options = qw(
 608             no_stream|x
 609             data_out|o=s
 610             compress|Z
 611         );
 612     }
 613     elsif ( $script eq "write_psl" )
 614     {
 615         @options = qw(
 616             no_stream|x
 617             data_out|o=s
 618             compress|Z
 619         );
 620     }
 621     elsif ( $script eq "write_2bit" )
 622     {
 623         @options = qw(
 624             no_stream|x
 625             data_out|o=s
 626             no_mask|N
 627         );
 628     }
 629     elsif ( $script eq "write_solid" )
 630     {
 631         @options = qw(
 632             wrap|w=s
 633             no_stream|x
 634             data_out|o=s
 635             compress|Z
 636         );
 637     }
 638     elsif ( $script eq "plot_seqlogo" )
 639     {
 640         @options = qw(
 641             no_stream|x
 642             data_out|o=s
 643         );
 644     }
 645     elsif ( $script eq "plot_phastcons_profiles" )
 646     {
 647         @options = qw(
 648             no_stream|x
 649             data_out|o=s
 650             genome|g=s
 651             mean|m
 652             median|M
 653             flank|f=s
 654             terminal|t=s
 655             title|T=s
 656             xlabel|X=s
 657             ylabel|Y=s
 658         );
 659     }
 660     elsif ( $script eq "analyze_vals" )
 661     {
 662         @options = qw(
 663             no_stream|x
 664             keys|k=s
 665         );
 666     }
 667     elsif ( $script eq "head_records" )
 668     {
 669         @options = qw(
 670             num|n=s
 671         );
 672     }
 673     elsif ( $script eq "remove_keys" )
 674     {
 675         @options = qw(
 676             keys|k=s
 677             save_keys|K=s
 678         );
 679     }
 680     elsif ( $script eq "rename_keys" )
 681     {
 682         @options = qw(
 683             keys|k=s
 684         );
 685     }
 686     elsif ( $script eq "uniq_vals" )
 687     {
 688         @options = qw(
 689             key|k=s
 690             invert|i
 691         );
 692     }
 693     elsif ( $script eq "merge_vals" )
 694     {
 695         @options = qw(
 696             keys|k=s
 697             delimit|d=s
 698         );
 699     }
 700     elsif ( $script eq "grab" )
 701     {
 702         @options = qw(
 703             patterns|p=s
 704             patterns_in|P=s
 705             regex|r=s
 706             eval|e=s
 707             exact_in|E=s
 708             invert|i
 709             case_insensitive|c
 710             keys|k=s
 711             keys_only|K
 712             vals_only|V
 713         );
 714     }
 715     elsif ( $script eq "compute" )
 716     {
 717         @options = qw(
 718             eval|e=s
 719         );
 720     }
 721     elsif ( $script eq "add_ident" )
 722     {
 723         @options = qw(
 724             prefix|p=s
 725             key|k=s
 726         );
 727     }
 728     elsif ( $script eq "count_records" )
 729     {
 730         @options = qw(
 731             no_stream|x
 732             data_out|o=s
 733         );
 734     }
 735     elsif ( $script eq "random_records" )
 736     {
 737         @options = qw(
 738             num|n=s
 739         );
 740     }
 741     elsif ( $script eq "sort_records" )
 742     {
 743         @options = qw(
 744             reverse|r
 745             keys|k=s
 746         );
 747     }
 748     elsif ( $script eq "count_vals" )
 749     {
 750         @options = qw(
 751             keys|k=s
 752         );
 753     }
 754     elsif ( $script eq "plot_histogram" )
 755     {
 756         @options = qw(
 757             no_stream|x
 758             data_out|o=s
 759             terminal|t=s
 760             title|T=s
 761             xlabel|X=s
 762             ylabel|Y=s
 763             key|k=s
 764             sort|s=s
 765         );
 766     }
 767     elsif ( $script eq "plot_lendist" )
 768     {
 769         @options = qw(
 770             no_stream|x
 771             data_out|o=s
 772             terminal|t=s
 773             title|T=s
 774             xlabel|X=s
 775             ylabel|Y=s
 776             key|k=s
 777         );
 778     }
 779     elsif ( $script eq "plot_chrdist" )
 780     {
 781         @options = qw(
 782             no_stream|x
 783             data_out|o=s
 784             terminal|t=s
 785             title|T=s
 786             xlabel|X=s
 787             ylabel|Y=s
 788         );
 789     }
 790     elsif ( $script eq "plot_karyogram" )
 791     {
 792         @options = qw(
 793             no_stream|x
 794             data_out|o=s
 795             genome|g=s
 796             feat_color|f=s
 797         );
 798     }
 799     elsif ( $script eq "plot_matches" )
 800     {
 801         @options = qw(
 802             no_stream|x
 803             data_out|o=s
 804             terminal|t=s
 805             title|T=s
 806             xlabel|X=s
 807             ylabel|Y=s
 808             direction|d=s
 809         );
 810     }
 811     elsif ( $script eq "length_vals" )
 812     {
 813         @options = qw(
 814             keys|k=s
 815         );
 816     }
 817     elsif ( $script eq "sum_vals" )
 818     {
 819         @options = qw(
 820             no_stream|x
 821             data_out|o=s
 822             keys|k=s
 823         );
 824     }
 825     elsif ( $script eq "mean_vals" )
 826     {
 827         @options = qw(
 828             no_stream|x
 829             data_out|o=s
 830             keys|k=s
 831         );
 832     }
 833     elsif ( $script eq "median_vals" )
 834     {
 835         @options = qw(
 836             no_stream|x
 837             data_out|o=s
 838             keys|k=s
 839         );
 840     }
 841     elsif ( $script eq "max_vals" )
 842     {
 843         @options = qw(
 844             no_stream|x
 845             data_out|o=s
 846             keys|k=s
 847         );
 848     }
 849     elsif ( $script eq "min_vals" )
 850     {
 851         @options = qw(
 852             no_stream|x
 853             data_out|o=s
 854             keys|k=s
 855         );
 856     }
 857     elsif ( $script eq "upload_to_ucsc" )
 858     {
 859         @options = qw(
 860             no_stream|x
 861             database|d=s
 862             table|t=s
 863             short_label|s=s
 864             long_label|l=s
 865             group|g=s
 866             priority|p=f
 867             use_score|u
 868             visibility|v=s
 869             wiggle|w
 870             color|c=s
 871             chunk_size|C=s
 872         );
 873     }
 874
 875     push @options, qw(
 876         stream_in|I=s
 877         stream_out|O=s
 878         verbose
 879     );
 880
 881 #    print STDERR Dumper( \@options );
 882
 883     GetOptions(
 884         \%options,
 885         @options,
 886     );
 887
 888     $options{ "cols" }      = [ split ",", $options{ "cols" } ]      if defined $options{ "cols" };
 889     $options{ "keys" }      = [ split ",", $options{ "keys" } ]      if defined $options{ "keys" };
 890     $options{ "no_keys" }   = [ split ",", $options{ "no_keys" } ]   if defined $options{ "no_keys" };
 891     $options{ "save_keys" } = [ split ",", $options{ "save_keys" } ] if defined $options{ "save_keys" };
 892     $options{ "quals" }     = [ split ",", $options{ "quals" } ]     if defined $options{ "quals" };
 893     $options{ "feats" }     = [ split ",", $options{ "feats" } ]     if defined $options{ "feats" };
 894     $options{ "frames" }    = [ split ",", $options{ "frames" } ]    if defined $options{ "frames" };
 895
 896     # ---- check arguments ----
 897
 898     if ( $options{ 'data_in' } )
 899     {
 900         $options{ "files" } = getopt_files( $options{ 'data_in' } );
 901
 902         Maasha::Common::error( qq(Argument to --data_in must be a valid file or fileglob expression) ) if scalar @{ $options{ "files" } } == 0;
 903     }
 904
 905     map { Maasha::Common::error( qq(Argument to --cols must be a whole numbers - not "$_") ) if $_ !~ /^\d+$/ } @{ $options{ "cols" } } if $options{ "cols" };
 906
 907 #    print STDERR Dumper( \%options );
 908
 909     foreach $opt ( keys %options )
 910     {
 911         if ( $opt =~ /stream_in|pattern_in|exact_in/ and not -f $options{ $opt } )
 912         {
 913             Maasha::Common::error( qq(Argument to --$opt must be a valid file or fileglob expression - not "$options{ $opt }") );
 914         }
 915         elsif ( $opt =~ /beg|end|word_size|wrap|chunk_size|tile_size|len|prefix_length|num|skip|cpus|window_size|step_size/ and $options{ $opt } !~ /^\d+$/ )
 916         {
 917             Maasha::Common::error( qq(Argument to --$opt must be a whole number - not "$options{ $opt }") );
 918         }
 919         elsif ( $opt =~ /max_hits|max_hits|max_misses|dist|edit_dist|flank|gap|hamming_dist|priority/ and $options{ $opt } !~ /^-?\d+$/ )
 920         {
 921             Maasha::Common::error( qq(Argument to --$opt must be an integer - not "$options{ $opt }") );
 922         }
 923         elsif ( $opt =~ /identity|threshold/ and $options{ $opt } !~ /^-?(?:\d+(?:\.\d*)?|\.\d+)$/ )
 924         {
 925             Maasha::Common::error( qq(Argument to --$opt must be a decimal number - not "$options{ $opt }") );
 926         }
 927         elsif ( $opt =~ /e_val/ and $options{ $opt } !~ /^([+-]?)(?=\d|\.\d)\d*(\.\d*)?([Ee]([+-]?\d+))?$/ )
 928         {
 929             Maasha::Common::error( qq(Argument to --$opt must be a float - not "$options{ $opt }") );
 930         }
 931         elsif ( $opt =~ /strand/ and $options{ $opt } !~ /^(\+|-)$/ )
 932         {
 933             Maasha::Common::error( qq(Argument to --$opt must be "+" or "-" - not "$options{ $opt }") );
 934         }
 935         elsif ( $opt eq "genome" )
 936         {
 937             @genomes = Maasha::Config::genomes();
 938
 939             if ( not grep $options{ $opt }, @genomes ) {
 940                 Maasha::Common::error( qq(Genome $options{ $opt } not found in "$ENV{ 'INST_DIR' }/conf/genomes.conf") );
 941             }
 942         }
 943         elsif ( $opt eq "terminal" and not $options{ $opt } =~ /^(svg|post|dumb)/ )
 944         {
 945             Maasha::Common::error( qq(Bad --$opt argument "$options{ $opt }") );
 946         }
 947         elsif ( $opt eq "table" and $options{ $opt } =~ /-\./ )
 948         {
 949             Maasha::Common::error( qq(Character '$options{ $opt }' is not allowed in table names) );
 950         }
 951     }
 952
 953     Maasha::Common::error( qq(no --database specified) )                if $script eq "create_blast_db"     and not $options{ "database" };
 954     Maasha::Common::error( qq(no --index_name specified) )              if $script eq "create_vmatch_index" and not $options{ "index_name" };
 955     Maasha::Common::error( qq(no --database or --genome specified) )    if $script eq "blast_seq" and not $options{ "genome" } and not $options{ "database" };
 956     Maasha::Common::error( qq(both --database and --genome specified) ) if $script eq "blast_seq" and $options{ "genome" } and $options{ "database" };
 957     Maasha::Common::error( qq(no --index_name or --genome specified) )  if $script eq "vmatch_seq" and not $options{ "genome" } and not $options{ "index_name" };
 958     Maasha::Common::error( qq(both --index and --genome specified) )    if $script eq "vmatch_seq" and $options{ "genome" } and $options{ "index" };
 959     Maasha::Common::error( qq(no --genome specified) )                  if $script =~ /get_genome_seq|get_genome_align|get_genome_phastcons|blat_seq|plot_phastcons_profiles|plot_karyogram/ and not $options{ "genome" };
 960     Maasha::Common::error( qq(no --key specified) )                     if $script =~ /plot_lendist|plot_histogram/ and not $options{ "key" };
 961     Maasha::Common::error( qq(no --keys speficied) )                    if $script =~ /sort_records|count_vals|sum_vals|mean_vals|median_vals|length_vals/ and not $options{ "keys" };
 962
 963     if ( $script eq "upload_to_ucsc" )
 964     {
 965         Maasha::Common::error( qq(no --database specified) ) if not $options{ "database" };
 966         Maasha::Common::error( qq(no --table specified) )    if not $options{ "table" };
 967     }
 968
 969     return wantarray ? %options : \%options;
 970 }
 971
 972
 973 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> SCRIPTS  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 974
 975
 976 sub script_print_usage
 977 {
 978     # Martin A. Hansen, January 2008.
 979
 980     # Retrieves usage information from file and
 981     # prints this nicely formatted.
 982
 983     my ( $path,   #   full path to usage file
 984        ) = @_;
 985
 986     # Returns nothing.
 987
 988     my ( $script, $fh, $line, @lines, @list, %hash, $key );
 989
 990     $script = ( split "/", $path )[ -1 ];
 991
 992     $fh = Maasha::Common::read_open( $path );
 993
 994     push @list, "Program name";
 995
 996     $hash{ "Program name" } = [ $script ];
 997
 998     while ( $line = <$fh> )
 999     {
1000         chomp $line;
1001
1002         $line =~ s/\$script/$script/g;
1003
1004         if ( $line =~ /^([^:]+):\s+(.+)$/ )
1005         {
1006             push @list, $1 if not exists $hash{ $1 };
1007             push @{ $hash{ $1 } }, $2;
1008         }
1009     }
1010
1011     close $fh;
1012
1013     print "\n";
1014
1015     foreach $key ( @list )
1016     {
1017         if ( scalar @{ $hash{ $key } } == 1 )
1018         {
1019             @lines = Maasha::Common::wrap_line( $hash{ $key }->[ 0 ], 80 );
1020
1021             printf( "%-15s%s\n", "$key:", shift @lines );
1022
1023             map { printf( "%-15s%s\n", "", $_ ) } @lines;
1024
1025             print "\n";
1026         }
1027         else
1028         {
1029             print "$key:\n";
1030
1031             map { print "   $_\n" } @{ $hash{ $key } };
1032
1033             print "\n";
1034         }
1035     }
1036
1037     exit;
1038 }
1039
1040
1041 sub script_list_biopieces
1042 {
1043     # Martin A. Hansen, January 2008.
1044
1045     # Prints the description from the usage for each of the biopieces.
1046
1047     my ( $path,    # full path to usage directory
1048        ) = @_;
1049
1050     # Returns nothing.
1051
1052     my ( @files, $file, $fh, $line, @lines, $program );
1053
1054     @files = Maasha::Common::ls_files( $path );
1055
1056     foreach $file ( sort @files )
1057     {
1058         $program = ( split "/", $file )[ -1 ];
1059
1060         $fh = Maasha::Common::read_open( $file );
1061
1062         while ( $line = <$fh> )
1063         {
1064             chomp $line;
1065
1066             if ( $line =~ /^Description:\s+(.+)/ )
1067             {
1068                 @lines = Maasha::Common::wrap_line( $1, 60 );
1069
1070                 printf( "%-30s%s\n", $program, shift @lines );
1071
1072                 map { printf( "%-30s%s\n", "", $_ ) } @lines;
1073             }
1074         }
1075
1076         close $fh;
1077     }
1078
1079     exit;
1080 }
1081
1082
1083 sub script_read_fasta
1084 {
1085     # Martin A. Hansen, August 2007.
1086
1087     # Read sequences from FASTA file.
1088
1089     my ( $in,        # handle to in stream
1090          $out,       # handle to out stream
1091          $options,   # options hash
1092        ) = @_;
1093
1094     # Returns nothing.
1095
1096     my ( $record, $file, $data_in, $entry, $num );
1097
1098     while ( $record = get_record( $in ) ) {
1099         put_record( $record, $out );
1100     }
1101
1102     $num = 1;
1103
1104     foreach $file ( @{ $options->{ "files" } } )
1105     {
1106         $data_in = Maasha::Common::read_open( $file );
1107
1108         while ( $entry = Maasha::Fasta::get_entry( $data_in ) )
1109         {
1110             if ( defined $entry->[ SEQ_NAME ] and $entry->[ SEQ ] )
1111             {
1112                 $record = {
1113                     SEQ_NAME => $entry->[ SEQ_NAME ],
1114                     SEQ      => $entry->[ SEQ ],
1115                     SEQ_LEN  => length $entry->[ SEQ ],
1116                 };
1117
1118                 put_record( $record, $out );
1119             }
1120
1121             goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1122
1123             $num++;
1124         }
1125
1126         close $data_in;
1127     }
1128
1129     NUM:
1130
1131     close $data_in if $data_in;
1132 }
1133
1134
1135 sub script_read_tab
1136 {
1137     # Martin A. Hansen, August 2007.
1138
1139     # Read table or table columns from stream or file.
1140
1141     my ( $in,        # handle to in stream
1142          $out,       # handle to out stream
1143          $options,   # options hash
1144        ) = @_;
1145
1146     # Returns nothing.
1147
1148     my ( $file, $line, @fields, @fields2, $i, $record, $data_in, $skip, $num );
1149
1150     $options->{ 'delimit' } ||= '\s+';
1151
1152     while ( $record = get_record( $in ) ) {
1153         put_record( $record, $out );
1154     }
1155
1156     $skip = $options->{ 'skip' } ||= 0;
1157     $num = 1;
1158
1159     foreach $file ( @{ $options->{ "files" } } )
1160     {
1161         $data_in = Maasha::Common::read_open( $file );
1162
1163         while ( $line = <$data_in> )
1164         {
1165             if ( $skip )
1166             {
1167                 $skip--;
1168                 next;
1169             }
1170
1171             next if $line =~ /^#|^$/;
1172
1173             chomp $line;
1174
1175             undef $record;
1176             undef @fields2;
1177
1178             @fields = split /$options->{'delimit'}/, $line;
1179
1180             if ( $options->{ "cols" } ) {
1181                 map { push @fields2, $fields[ $_ ] } @{ $options->{ "cols" } };
1182             } else {
1183                 @fields2 = @fields;
1184             }
1185
1186             for ( $i = 0; $i < @fields2; $i++ )
1187             {
1188                 if ( $options->{ "keys" }->[ $i ] ) {
1189                     $record->{ $options->{ "keys" }->[ $i ] } = $fields2[ $i ];
1190                 } else {
1191                     $record->{ "V" . $i } = $fields2[ $i ];
1192                 }
1193             }
1194
1195             put_record( $record, $out );
1196
1197             goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1198
1199             $num++;
1200         }
1201
1202         close $data_in;
1203     }
1204
1205     NUM:
1206
1207     close $data_in if $data_in;
1208 }
1209
1210
1211 sub script_read_psl
1212 {
1213     # Martin A. Hansen, August 2007.
1214
1215     # Read psl table from stream or file.
1216
1217     my ( $in,        # handle to in stream
1218          $out,       # handle to out stream
1219          $options,   # options hash
1220        ) = @_;
1221
1222     # Returns nothing.
1223
1224     my ( $record, @files, $file, $entries, $entry, $num );
1225
1226     while ( $record = get_record( $in ) ) {
1227         put_record( $record, $out );
1228     }
1229
1230     $num = 1;
1231
1232     foreach $file ( @{ $options->{ "files" } } )
1233     {
1234         $entries = Maasha::UCSC::psl_get_entries( $file );
1235
1236         foreach $entry ( @{ $entries } )
1237         {
1238             put_record( $entry, $out );
1239
1240             goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1241
1242             $num++;
1243         }
1244     }
1245
1246     NUM:
1247 }
1248
1249
1250 sub script_read_bed
1251 {
1252     # Martin A. Hansen, August 2007.
1253
1254     # Read bed table from stream or file.
1255
1256     my ( $in,        # handle to in stream
1257          $out,       # handle to out stream
1258          $options,   # options hash
1259        ) = @_;
1260
1261     # Returns nothing.
1262
1263     my ( $file, $record, $entry, $data_in, $num );
1264
1265     while ( $record = get_record( $in ) ) {
1266         put_record( $record, $out );
1267     }
1268
1269     $num = 1;
1270
1271     foreach $file ( @{ $options->{ "files" } } )
1272     {
1273         $data_in = Maasha::Common::read_open( $file );
1274
1275         while ( $entry = Maasha::UCSC::bed_get_entry( $data_in ) )
1276         {
1277             put_record( $entry, $out );
1278
1279             goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1280
1281             $num++;
1282         }
1283
1284         close $data_in;
1285     }
1286
1287     NUM:
1288
1289     close $data_in if $data_in;
1290 }
1291
1292
1293 sub script_read_blast_tab
1294 {
1295     # Martin A. Hansen, September 2007.
1296
1297     # Read tabular BLAST output from NCBI blast run with -m8 or -m9.
1298
1299     my ( $in,        # handle to in stream
1300          $out,       # handle to out stream
1301          $options,   # options hash
1302        ) = @_;
1303
1304     # Returns nothing.
1305
1306     my ( $file, $line, @fields, $strand, $record, $data_in, $num );
1307
1308     while ( $record = get_record( $in ) ) {
1309         put_record( $record, $out );
1310     }
1311
1312     $num = 1;
1313
1314     foreach $file ( @{ $options->{ "files" } } )
1315     {
1316         $data_in = Maasha::Common::read_open( $file );
1317
1318         while ( $line = <$data_in> )
1319         {
1320             chomp $line;
1321
1322             next if $line =~ /^#/;
1323
1324             @fields = split /\t/, $line;
1325
1326             $record->{ "REC_TYPE" }   = "BLAST";
1327             $record->{ "Q_ID" }       = $fields[ 0 ];
1328             $record->{ "S_ID" }       = $fields[ 1 ];
1329             $record->{ "IDENT" }      = $fields[ 2 ];
1330             $record->{ "ALIGN_LEN" }  = $fields[ 3 ];
1331             $record->{ "MISMATCHES" } = $fields[ 4 ];
1332             $record->{ "GAPS" }       = $fields[ 5 ];
1333             $record->{ "Q_BEG" }      = $fields[ 6 ] - 1; # BLAST is 1-based
1334             $record->{ "Q_END" }      = $fields[ 7 ] - 1; # BLAST is 1-based
1335             $record->{ "S_BEG" }      = $fields[ 8 ] - 1; # BLAST is 1-based
1336             $record->{ "S_END" }      = $fields[ 9 ] - 1; # BLAST is 1-based
1337             $record->{ "E_VAL" }      = $fields[ 10 ];
1338             $record->{ "BIT_SCORE" }  = $fields[ 11 ];
1339
1340             if ( $record->{ "S_BEG" } > $record->{ "S_END" } )
1341             {
1342                 $record->{ "STRAND" } = '-';
1343
1344                 ( $record->{ "S_BEG" }, $record->{ "S_END" } ) = ( $record->{ "S_END" }, $record->{ "S_BEG" } );
1345             }
1346             else
1347             {
1348                 $record->{ "STRAND" } = '+';
1349             }
1350
1351             put_record( $record, $out );
1352
1353             goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1354
1355             $num++;
1356         }
1357
1358         close $data_in;
1359     }
1360
1361     NUM:
1362
1363     close $data_in if $data_in;
1364 }
1365
1366
1367 sub script_read_embl
1368 {
1369     # Martin A. Hansen, August 2007.
1370
1371     # Read EMBL format.
1372
1373     my ( $in,        # handle to in stream
1374          $out,       # handle to out stream
1375          $options,   # options hash
1376        ) = @_;
1377
1378     # Returns nothing.
1379
1380     my ( %options2, $file, $data_in, $num, $entry, $record );
1381
1382     map { $options2{ "keys" }{ $_ } = 1 }  @{ $options->{ "keys" } };
1383     map { $options2{ "feats" }{ $_ } = 1 } @{ $options->{ "feats" } };
1384     map { $options2{ "quals" }{ $_ } = 1 } @{ $options->{ "quals" } };
1385
1386     while ( $record = get_record( $in ) ) {
1387         put_record( $record, $out );
1388     }
1389
1390     $num = 1;
1391
1392     foreach $file ( @{ $options->{ "files" } } )
1393     {
1394         $data_in = Maasha::Common::read_open( $file );
1395
1396         while ( $entry = Maasha::EMBL::get_embl_entry( $data_in ) )
1397         {
1398             $record = Maasha::EMBL::parse_embl_entry( $entry, \%options2 );
1399
1400             my ( $feat, $feat2, $qual, $qual_val, $record_copy );
1401
1402             $record_copy = dclone $record;
1403
1404             delete $record_copy->{ "FT" };
1405
1406             put_record( $record_copy, $out );
1407
1408             delete $record_copy->{ "SEQ" };
1409
1410             foreach $feat ( keys %{ $record->{ "FT" } } )
1411             {
1412                 $record_copy->{ "FEAT_TYPE" } = $feat;
1413
1414                 foreach $feat2 ( @{ $record->{ "FT" }->{ $feat } } )
1415                 {
1416                     foreach $qual ( keys %{ $feat2 } )
1417                     {
1418                         $qual_val = join "; ", @{ $feat2->{ $qual } };
1419
1420                         $qual =~ s/^_//;
1421                         $qual = uc $qual;
1422
1423                         $record_copy->{ $qual } = $qual_val;
1424                     }
1425
1426                     put_record( $record_copy, $out );
1427                 }
1428             }
1429
1430             goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1431
1432             $num++;
1433         }
1434
1435         close $data_in;
1436     }
1437
1438     NUM:
1439
1440     close $data_in if $data_in;
1441 }
1442
1443
1444 sub script_read_stockholm
1445 {
1446     # Martin A. Hansen, August 2007.
1447
1448     # Read Stockholm format.
1449
1450     my ( $in,        # handle to in stream
1451          $out,       # handle to out stream
1452          $options,   # options hash
1453        ) = @_;
1454
1455     # Returns nothing.
1456
1457     my ( $data_in, $file, $num, $entry, $record, $record_anno, $record_align, $key, $seq );
1458
1459     while ( $record = get_record( $in ) ) {
1460         put_record( $record, $out );
1461     }
1462
1463     $num = 1;
1464
1465     foreach $file ( @{ $options->{ "files" } } )
1466     {
1467         $data_in = Maasha::Common::read_open( $file );
1468
1469         while ( $entry = Maasha::Stockholm::get_stockholm_entry( $data_in ) )
1470         {
1471             $record = Maasha::Stockholm::parse_stockholm_entry( $entry );
1472
1473             undef $record_anno;
1474
1475             foreach $key ( keys %{ $record->{ "GF" } } ) {
1476                 $record_anno->{ $key } = $record->{ "GF" }->{ $key };
1477             }
1478
1479             $record_anno->{ "ALIGN" } = $num;
1480
1481             put_record( $record_anno, $out );
1482
1483             foreach $seq ( @{ $record->{ "ALIGN" } } )
1484             {
1485                 undef $record_align;
1486
1487                 $record_align = {
1488                     SEQ_NAME  => $seq->[ 0 ],
1489                     SEQ       => $seq->[ 1 ],
1490                 };
1491
1492                 put_record( $record_align, $out );
1493             }
1494
1495             goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1496
1497             $num++;
1498         }
1499
1500         close $data_in;
1501     }
1502
1503     NUM:
1504
1505     close $data_in if $data_in;
1506 }
1507
1508
1509 sub script_read_phastcons
1510 {
1511     # Martin A. Hansen, December 2007.
1512
1513     # Read PhastCons format.
1514
1515     my ( $in,        # handle to in stream
1516          $out,       # handle to out stream
1517          $options,   # options hash
1518        ) = @_;
1519
1520     # Returns nothing.
1521
1522     my ( $data_in, $file, $num, $entry, @records, $record );
1523
1524     $options->{ "min" }       ||= 10;
1525     $options->{ "dist" }      ||= 25;
1526     $options->{ "threshold" } ||= 0.8;
1527     $options->{ "gap" }       ||= 5;
1528
1529     while ( $record = get_record( $in ) ) {
1530         put_record( $record, $out );
1531     }
1532
1533     $num = 1;
1534
1535     foreach $file ( @{ $options->{ "files" } } )
1536     {
1537         $data_in = Maasha::Common::read_open( $file );
1538
1539         while ( $entry = Maasha::UCSC::phastcons_get_entry( $data_in ) )
1540         {
1541             @records = Maasha::UCSC::phastcons_parse_entry( $entry, $options );
1542
1543             foreach $record ( @records )
1544             {
1545                 $record->{ "REC_TYPE" } = "BED";
1546                 $record->{ "BED_LEN" }  = $record->{ "CHR_END" } - $record->{ "CHR_BEG" } + 1;
1547
1548                 put_record( $record, $out );
1549
1550                 goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1551
1552                 $num++;
1553             }
1554         }
1555
1556         close $data_in;
1557     }
1558
1559     NUM:
1560
1561     close $data_in if $data_in;
1562 }
1563
1564
1565 sub script_read_soft
1566 {
1567     # Martin A. Hansen, December 2007.
1568
1569     # Read soft format.
1570     # http://www.ncbi.nlm.nih.gov/geo/info/soft2.html
1571
1572     my ( $in,        # handle to in stream
1573          $out,       # handle to out stream
1574          $options,   # options hash
1575        ) = @_;
1576
1577     # Returns nothing.
1578
1579     my ( $data_in, $file, $num, $records, $record, $soft_index, $fh, @platforms, $plat_table, @samples, $sample, $old_end );
1580
1581     while ( $record = get_record( $in ) ) {
1582         put_record( $record, $out );
1583     }
1584
1585     $num = 1;
1586
1587     foreach $file ( @{ $options->{ "files" } } )
1588     {
1589         $soft_index = Maasha::NCBI::soft_index_file( $file );
1590
1591         $fh         = Maasha::Common::read_open( $file );
1592
1593         @platforms  = grep { $_->[ 0 ] =~ /PLATFORM/ } @{ $soft_index };
1594
1595         $plat_table = Maasha::NCBI::soft_get_platform( $fh, $platforms[ 0 ]->[ 1 ], $platforms[ -1 ]->[ 2 ] );
1596
1597         @samples    = grep { $_->[ 0 ] =~ /SAMPLE/ } @{ $soft_index };
1598
1599         $old_end    = $platforms[ -1 ]->[ 2 ];
1600
1601         foreach $sample ( @samples )
1602         {
1603             $records = Maasha::NCBI::soft_get_sample( $fh, $plat_table, $sample->[ 1 ] - $old_end - 1, $sample->[ 2 ] - $old_end - 1 );
1604
1605             foreach $record ( @{ $records } )
1606             {
1607                 put_record( $record, $out );
1608
1609                 goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1610
1611                 $num++;
1612             }
1613
1614             $old_end = $sample->[ 2 ];
1615         }
1616
1617         close $fh;
1618     }
1619
1620     NUM:
1621
1622     close $data_in if $data_in;
1623     close $fh if $fh;
1624 }
1625
1626
1627 sub script_read_gff
1628 {
1629     # Martin A. Hansen, February 2008.
1630
1631     # Read soft format.
1632     # http://www.ncbi.nlm.nih.gov/geo/info/soft2.html
1633
1634     my ( $in,        # handle to in stream
1635          $out,       # handle to out stream
1636          $options,   # options hash
1637        ) = @_;
1638
1639     # Returns nothing.
1640
1641     my ( $data_in, $file, $fh, $num, $record, $entry );
1642
1643     while ( $record = get_record( $in ) ) {
1644         put_record( $record, $out );
1645     }
1646
1647     $num = 1;
1648
1649     foreach $file ( @{ $options->{ "files" } } )
1650     {
1651         $fh = Maasha::Common::read_open( $file );
1652
1653         while ( $entry = Maasha::GFF::get_entry( $fh ) )
1654         {
1655             put_record( $entry, $out );
1656
1657             goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1658
1659             $num++;
1660         }
1661
1662         close $fh;
1663     }
1664
1665     NUM:
1666
1667     close $data_in if $data_in;
1668 }
1669
1670
1671 sub script_read_2bit
1672 {
1673     # Martin A. Hansen, March 2008.
1674
1675     # Read sequences from 2bit file.
1676
1677     my ( $in,        # handle to in stream
1678          $out,       # handle to out stream
1679          $options,   # options hash
1680        ) = @_;
1681
1682     # Returns nothing.
1683
1684     my ( $record, $file, $data_in, $mask, $toc, $line, $num );
1685
1686     $mask = 1 if not $options->{ "no_mask" };
1687
1688     while ( $record = get_record( $in ) ) {
1689         put_record( $record, $out );
1690     }
1691
1692     $num = 1;
1693
1694     foreach $file ( @{ $options->{ "files" } } )
1695     {
1696         $data_in = Maasha::Common::read_open( $file );
1697
1698         $toc = Maasha::TwoBit::twobit_get_TOC( $data_in );
1699
1700         foreach $line ( @{ $toc } )
1701         {
1702             $record->{ "SEQ_NAME" } = $line->[ 0 ];
1703             $record->{ "SEQ" }      = Maasha::TwoBit::twobit_get_seq( $data_in, $line->[ 1 ], undef, undef, $mask );
1704             $record->{ "SEQ_LEN" }  = length $record->{ "SEQ" };
1705
1706             put_record( $record, $out );
1707
1708             goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1709
1710             $num++;
1711         }
1712
1713         close $data_in;
1714     }
1715
1716     NUM:
1717
1718     close $data_in if $data_in;
1719 }
1720
1721
1722 sub script_read_solexa
1723 {
1724     # Martin A. Hansen, March 2008.
1725
1726     # Read Solexa sequence reads from file.
1727
1728     my ( $in,        # handle to in stream
1729          $out,       # handle to out stream
1730          $options,   # options hash
1731        ) = @_;
1732
1733     # Returns nothing.
1734
1735     my ( $record, $file, $base_name, $data_in, $line, $num, @fields, @seqs, @scores, $i, $seq, $seq_count );
1736
1737     $options->{ "quality" } ||= 20;
1738
1739     while ( $record = get_record( $in ) ) {
1740         put_record( $record, $out );
1741     }
1742
1743     $num = 1;
1744
1745     foreach $file ( @{ $options->{ "files" } } )
1746     {
1747         $data_in   = Maasha::Common::read_open( $file );
1748         $base_name = Maasha::Common::get_basename( $file );
1749         $base_name =~ s/\..*//;
1750
1751         $seq_count = 0;
1752
1753         while ( $line = <$data_in> )
1754         {
1755             @fields = split /:/, $line;
1756             @seqs   = split //, $fields[ 5 ];
1757             @scores = split / /, $fields[ -1 ];
1758
1759             for ( $i = 0; $i < @scores; $i++ ) {
1760                 $seqs[ $i ] = lc $seqs[ $i ] if $scores[ $i ] < $options->{ "quality" };
1761             }
1762
1763             $seq = join "", @seqs;
1764
1765             $record->{ "SEQ_NAME" }     = sprintf( "%s_ID%08d", $base_name, $seq_count );
1766             $record->{ "SEQ" }          = $seq;
1767             $record->{ "SEQ_LEN" }      = length $seq;
1768             $record->{ "SCORE_MEAN" }   = sprintf ( "%.2f", Maasha::Calc::mean( \@scores ) );
1769
1770             put_record( $record, $out );
1771
1772             goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1773
1774             $seq_count++;
1775             $num++;
1776         }
1777
1778         close $data_in;
1779     }
1780
1781     NUM:
1782
1783     close $data_in if $data_in;
1784 }
1785
1786
1787 sub script_read_solid
1788 {
1789     # Martin A. Hansen, April 2008.
1790
1791     # Read Solid sequence from file.
1792
1793     my ( $in,        # handle to in stream
1794          $out,       # handle to out stream
1795          $options,   # options hash
1796        ) = @_;
1797
1798     # Returns nothing.
1799
1800     my ( $record, $file, $data_in, $line, $num, $seq_name, $seq_cs, $seq_qual, @scores, @seqs, $i );
1801
1802     $options->{ "quality" } ||= 15;
1803
1804     while ( $record = get_record( $in ) ) {
1805         put_record( $record, $out );
1806     }
1807
1808     $num = 1;
1809
1810     foreach $file ( @{ $options->{ "files" } } )
1811     {
1812         $data_in = Maasha::Common::read_open( $file );
1813
1814         while ( $line = <$data_in> )
1815         {
1816             chomp $line;
1817
1818             ( $seq_name, $seq_cs, $seq_qual ) = split /\t/, $line;
1819
1820             @scores = split /,/, $seq_qual;
1821             @seqs   = split //, Maasha::Solid::color_space2seq( $seq_cs );
1822
1823             for ( $i = 0; $i < @seqs; $i++ ) {
1824                 $seqs[ $i ] = lc $seqs[ $i ] if $scores[ $i ] < $options->{ "quality" };
1825             }
1826
1827             $record = {
1828                 SEQ_NAME   => $seq_name,
1829                 SEQ_CS     => $seq_cs,
1830                 SEQ_QUAL   => $seq_qual,
1831                 SEQ_LEN    => length $seq_cs,
1832                 SEQ        => join( "", @seqs ),
1833                 SCORE_MEAN => sprintf( "%.2f", Maasha::Calc::mean( \@scores ) ),
1834             };
1835
1836             put_record( $record, $out );
1837
1838             goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1839
1840             $num++;
1841         }
1842
1843         close $data_in;
1844     }
1845
1846     NUM:
1847
1848     close $data_in if $data_in;
1849 }
1850
1851
1852 sub script_read_mysql
1853 {
1854     # Martin A. Hansen, May 2008.
1855
1856     # Read a MySQL query into stream.
1857
1858     my ( $in,        # handle to in stream
1859          $out,       # handle to out stream
1860          $options,   # options hash
1861        ) = @_;
1862
1863     # Returns nothing.
1864
1865     my ( $record, $dbh, $results );
1866
1867     $options->{ "user" }     ||= Maasha::UCSC::ucsc_get_user();
1868     $options->{ "password" } ||= Maasha::UCSC::ucsc_get_password();
1869
1870     while ( $record = get_record( $in ) ) {
1871         put_record( $record, $out );
1872     }
1873
1874     $dbh = Maasha::SQL::connect( $options->{ "database" }, $options->{ "user" }, $options->{ "password" } );
1875
1876     $results = Maasha::SQL::query_hashref_list( $dbh, $options->{ "query" } );
1877
1878     Maasha::SQL::disconnect( $dbh );
1879
1880     map { put_record( $_ ) } @{ $results };
1881 }
1882
1883
1884 sub script_length_seq
1885 {
1886     # Martin A. Hansen, August 2007.
1887
1888     # Determine the length of sequences in stream.
1889
1890     my ( $in,        # handle to in stream
1891          $out,       # handle to out stream
1892          $options,   # options hash
1893        ) = @_;
1894
1895     # Returns nothing.
1896
1897     my ( $record, $total );
1898
1899     while ( $record = get_record( $in ) )
1900     {
1901         if ( $record->{ "SEQ" } )
1902         {
1903             $record->{ "SEQ_LEN" } = length $record->{ "SEQ" };
1904             $total += $record->{ "SEQ_LEN" };
1905         }
1906
1907         put_record( $record, $out ) if not $options->{ "no_stream" };
1908     }
1909
1910     put_record( { TOTAL_SEQ_LEN => $total }, $out );
1911 }
1912
1913
1914 sub script_uppercase_seq
1915 {
1916     # Martin A. Hansen, August 2007.
1917
1918     # Uppercases sequences in stream.
1919
1920     my ( $in,    # handle to in stream
1921          $out,   # handle to out stream
1922        ) = @_;
1923
1924     # Returns nothing.
1925
1926     my ( $record );
1927
1928     while ( $record = get_record( $in ) )
1929     {
1930         $record->{ "SEQ" } = uc $record->{ "SEQ" } if $record->{ "SEQ" };
1931
1932         put_record( $record, $out );
1933     }
1934 }
1935
1936
1937 sub script_shuffle_seq
1938 {
1939     # Martin A. Hansen, December 2007.
1940
1941     # Shuffle sequences in stream.
1942
1943     my ( $in,    # handle to in stream
1944          $out,   # handle to out stream
1945        ) = @_;
1946
1947     # Returns nothing.
1948
1949     my ( $record );
1950
1951     while ( $record = get_record( $in ) )
1952     {
1953         $record->{ "SEQ" } = Maasha::Seq::seq_shuffle( $record->{ "SEQ" } ) if $record->{ "SEQ" };
1954
1955         put_record( $record, $out );
1956     }
1957 }
1958
1959
1960 sub script_analyze_seq
1961 {
1962     # Martin A. Hansen, August 2007.
1963
1964     # Analyze sequence composition of sequences in stream.
1965
1966     my ( $in,     # handle to in stream
1967          $out,    # handle to out stream
1968        ) = @_;
1969
1970     # Returns nothing.
1971
1972     my ( $record, $analysis );
1973
1974     while ( $record = get_record( $in ) )
1975     {
1976         if ( $record->{ "SEQ" } )
1977         {
1978             $analysis = Maasha::Seq::seq_analyze( $record->{ "SEQ" } );
1979
1980             map { $record->{ $_ } = $analysis->{ $_ } } keys %{ $analysis };
1981         }
1982
1983         put_record( $record, $out );
1984     }
1985 }
1986
1987
1988 sub script_analyze_tags
1989 {
1990     # Martin A. Hansen, August 2008.
1991
1992     # Analyze sequence tags in stream.
1993
1994     my ( $in,     # handle to in stream
1995          $out,    # handle to out stream
1996        ) = @_;
1997
1998     # Returns nothing.
1999
2000     my ( $record, $analysis, %len_hash, %clone_hash, $clones, $key, $tag_record );
2001
2002     while ( $record = get_record( $in ) )
2003     {
2004         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } )
2005         {
2006             if ( $record->{ "SEQ_NAME" } =~ /_(\d+)$/ )
2007             {
2008                 $clones = $1;
2009
2010                 $len_hash{ length( $record->{ "SEQ" } ) }++;
2011                 $clone_hash{ length( $record->{ "SEQ" } ) } += $clones;
2012             }
2013         }
2014         elsif ( $record->{ "Q_ID" } and $record->{ "BED_LEN" } )
2015         {
2016             if ( $record->{ "Q_ID" } =~ /_(\d+)$/ )
2017             {
2018                 $clones = $1;
2019
2020                 $len_hash{ $record->{ "BED_LEN" } }++;
2021                 $clone_hash{ $record->{ "BED_LEN" } } += $clones;
2022             }
2023         }
2024     }
2025
2026     foreach $key ( sort { $a <=> $b } keys %len_hash )
2027     {
2028         $tag_record->{ "TAG_LEN" }    = $key;
2029         $tag_record->{ "TAG_COUNT" }  = $len_hash{ $key };
2030         $tag_record->{ "TAG_CLONES" } = $clone_hash{ $key };
2031
2032         put_record( $tag_record, $out );
2033     }
2034 }
2035
2036
2037 sub script_complexity_seq
2038 {
2039     # Martin A. Hansen, May 2008.
2040
2041     # Generates an index calculated as the most common di-residue over
2042     # the sequence length for all sequences in stream.
2043
2044     my ( $in,     # handle to in stream
2045          $out,    # handle to out stream
2046        ) = @_;
2047
2048     # Returns nothing.
2049
2050     my ( $record, $index );
2051
2052     while ( $record = get_record( $in ) )
2053     {
2054         $record->{ "SEQ_COMPLEXITY" } = sprintf( "%.2f", Maasha::Seq::seq_complexity( $record->{ "SEQ" } ) ) if $record->{ "SEQ" };
2055
2056         put_record( $record, $out );
2057     }
2058 }
2059
2060
2061 sub script_oligo_freq
2062 {
2063     # Martin A. Hansen, August 2007.
2064
2065     # Determine the length of sequences in stream.
2066
2067     my ( $in,        # handle to in stream
2068          $out,       # handle to out stream
2069          $options,   # options hash
2070        ) = @_;
2071
2072     # Returns nothing.
2073
2074     my ( $record, %oligos, @freq_table );
2075
2076     $options->{ "word_size" } ||= 7;
2077
2078     while ( $record = get_record( $in ) )
2079     {
2080         if ( $record->{ "SEQ" } )
2081         {
2082             map { $oligos{ $_ }++ } Maasha::Seq::seq2oligos( \$record->{ "SEQ" }, $options->{ "word_size" } );
2083
2084             if ( not $options->{ "all" } )
2085             {
2086                 @freq_table = Maasha::Seq::oligo_freq( \%oligos );
2087
2088                 map { put_record( $_, $out ) } @freq_table;
2089
2090                 undef %oligos;
2091             }
2092         }
2093
2094         put_record( $record, $out );
2095     }
2096
2097     if ( $options->{ "all" } )
2098     {
2099         @freq_table = Maasha::Seq::oligo_freq( \%oligos );
2100
2101         map { put_record( $_, $out ) } @freq_table;
2102     }
2103 }
2104
2105
2106 sub script_create_weight_matrix
2107 {
2108     # Martin A. Hansen, August 2007.
2109
2110     # Creates a weight matrix from an alignmnet.
2111
2112     my ( $in,        # handle to in stream
2113          $out,       # handle to out stream
2114          $options,   # options hash
2115        ) = @_;
2116
2117     # Returns nothing.
2118
2119     my ( $record, $count, $i, $res, %freq_hash, %res_hash, $freq );
2120
2121     $count = 0;
2122
2123     while ( $record = get_record( $in ) )
2124     {
2125         if ( $record->{ "SEQ" } )
2126         {
2127             for ( $i = 0; $i < length $record->{ "SEQ" }; $i++ )
2128             {
2129                 $res = substr $record->{ "SEQ" }, $i, 1;
2130
2131                 $freq_hash{ $i }{ $res }++;
2132                 $res_hash{ $res } = 1;
2133             }
2134
2135             $count++;
2136         }
2137         else
2138         {
2139             put_record( $record, $out );
2140         }
2141     }
2142
2143     foreach $res ( sort keys %res_hash )
2144     {
2145         undef $record;
2146
2147         $record->{ "V0" } = $res;
2148
2149         for ( $i = 0; $i < keys %freq_hash; $i++ )
2150         {
2151             $freq = $freq_hash{ $i }{ $res } || 0;
2152
2153             if ( $options->{ "percent" } ) {
2154                 $freq = sprintf( "%.0f", 100 * $freq / $count ) if $freq > 0;
2155             }
2156
2157             $record->{ "V" . ( $i + 1 ) } = $freq;
2158         }
2159
2160         put_record( $record, $out );
2161     }
2162 }
2163
2164
2165 sub script_calc_bit_scores
2166 {
2167     # Martin A. Hansen, March 2007.
2168
2169     # Calculates the bit scores for each position from an alignmnet in the stream.
2170
2171     my ( $in,        # handle to in stream
2172          $out,       # handle to out stream
2173        ) = @_;
2174
2175     # Returns nothing.
2176
2177     my ( $record, $type, $count, $i, $res, %freq_hash, $bit_max, $bit_height, $bit_diff );
2178
2179     $count = 0;
2180
2181     while ( $record = get_record( $in ) )
2182     {
2183         if ( $record->{ "SEQ" } )
2184         {
2185             $type = Maasha::Seq::seq_guess_type( $record->{ "SEQ" } ) if not $type;
2186
2187             for ( $i = 0; $i < length $record->{ "SEQ" }; $i++ )
2188             {
2189                 $res = substr $record->{ "SEQ" }, $i, 1;
2190
2191                 next if $res =~ /-|_|~|\./;
2192
2193                 $freq_hash{ $i }{ $res }++;
2194             }
2195
2196             $count++;
2197         }
2198         else
2199         {
2200             put_record( $record, $out );
2201         }
2202     }
2203
2204     undef $record;
2205
2206     if ( $type eq "protein" ) {
2207         $bit_max = 4;
2208     } else {
2209         $bit_max = 2;
2210     }
2211
2212     for ( $i = 0; $i < keys %freq_hash; $i++ )
2213     {
2214         $bit_height = Maasha::Seq::seqlogo_calc_bit_height( $freq_hash{ $i }, $count );
2215
2216         $bit_diff = $bit_max - $bit_height;
2217
2218         $record->{ "V" . ( $i ) } = sprintf( "%.2f", $bit_diff );
2219     }
2220
2221     put_record( $record, $out );
2222 }
2223
2224
2225 sub script_reverse_seq
2226 {
2227     # Martin A. Hansen, August 2007.
2228
2229     # Reverse sequence in record.
2230
2231     my ( $in,    # handle to in stream
2232          $out,   # handle to out stream
2233        ) = @_;
2234
2235     # Returns nothing.
2236
2237     my ( $record );
2238
2239     while ( $record = get_record( $in ) )
2240     {
2241         if ( $record->{ "SEQ" } ) {
2242             $record->{ "SEQ" } = reverse $record->{ "SEQ" };
2243         }
2244
2245         put_record( $record, $out );
2246     }
2247 }
2248
2249
2250 sub script_complement_seq
2251 {
2252     # Martin A. Hansen, August 2007.
2253
2254     # Complement sequence in record.
2255
2256     my ( $in,     # handle to in stream
2257          $out,    # handle to out stream
2258        ) = @_;
2259
2260     # Returns nothing.
2261
2262     my ( $record, $type );
2263
2264     while ( $record = get_record( $in ) )
2265     {
2266         if ( $record->{ "SEQ" } )
2267         {
2268             if ( not $type ) {
2269                 $type = Maasha::Seq::seq_guess_type( $record->{ "SEQ" } );
2270             }
2271
2272             if ( $type eq "rna" ) {
2273                 Maasha::Seq::rna_comp( \$record->{ "SEQ" } );
2274             } elsif ( $type eq "dna" ) {
2275                 Maasha::Seq::dna_comp( \$record->{ "SEQ" } );
2276             }
2277         }
2278
2279         put_record( $record, $out );
2280     }
2281 }
2282
2283
2284 sub script_remove_indels
2285 {
2286     # Martin A. Hansen, August 2007.
2287
2288     # Remove indels from sequences in stream.
2289
2290     my ( $in,     # handle to in stream
2291          $out,    # handle to out stream
2292        ) = @_;
2293
2294     # Returns nothing.
2295
2296     my ( $record );
2297
2298     while ( $record = get_record( $in ) )
2299     {
2300         $record->{ 'SEQ' } =~ tr/-~.//d if $record->{ "SEQ" };
2301
2302         put_record( $record, $out );
2303     }
2304 }
2305
2306
2307 sub script_transliterate_seq
2308 {
2309     # Martin A. Hansen, August 2007.
2310
2311     # Transliterate chars from sequence in record.
2312
2313     my ( $in,        # handle to in stream
2314          $out,       # handle to out stream
2315          $options,   # options hash
2316        ) = @_;
2317
2318     # Returns nothing.
2319
2320     my ( $record, $search, $replace, $delete );
2321
2322     $search  = $options->{ "search" }  || "";
2323     $replace = $options->{ "replace" } || "";
2324     $delete  = $options->{ "delete" }  || "";
2325
2326     while ( $record = get_record( $in ) )
2327     {
2328         if ( $record->{ "SEQ" } )
2329         {
2330             if ( $search and $replace ) {
2331                 eval "\$record->{ 'SEQ' } =~ tr/$search/$replace/";
2332             } elsif ( $delete ) {
2333                 eval "\$record->{ 'SEQ' } =~ tr/$delete//d";
2334             }
2335         }
2336
2337         put_record( $record, $out );
2338     }
2339 }
2340
2341
2342 sub script_transliterate_vals
2343 {
2344     # Martin A. Hansen, April 2008.
2345
2346     # Transliterate chars from values in record.
2347
2348     my ( $in,        # handle to in stream
2349          $out,       # handle to out stream
2350          $options,   # options hash
2351        ) = @_;
2352
2353     # Returns nothing.
2354
2355     my ( $record, $search, $replace, $delete, $key );
2356
2357     $search  = $options->{ "search" }  || "";
2358     $replace = $options->{ "replace" } || "";
2359     $delete  = $options->{ "delete" }  || "";
2360
2361     while ( $record = get_record( $in ) )
2362     {
2363         foreach $key ( @{ $options->{ "keys" } } )
2364         {
2365             if ( exists $record->{ $key } )
2366             {
2367                 if ( $search and $replace ) {
2368                     eval "\$record->{ $key } =~ tr/$search/$replace/";
2369                 } elsif ( $delete ) {
2370                     eval "\$record->{ $key } =~ tr/$delete//d";
2371                 }
2372             }
2373         }
2374
2375         put_record( $record, $out );
2376     }
2377 }
2378
2379
2380 sub script_translate_seq
2381 {
2382     # Martin A. Hansen, February 2008.
2383
2384     # Translate DNA sequence into protein sequence.
2385
2386     my ( $in,        # handle to in stream
2387          $out,       # handle to out stream
2388          $options,   # options hash
2389        ) = @_;
2390
2391     # Returns nothing.
2392
2393     my ( $record, $frame, %new_record );
2394
2395     $options->{ "frames" } ||= [ 1, 2, 3, -1, -2, -3 ];
2396
2397     while ( $record = get_record( $in ) )
2398     {
2399         if ( $record->{ "SEQ" } )
2400         {
2401             if ( Maasha::Seq::seq_guess_type( $record->{ "SEQ" } ) eq "dna" )
2402             {
2403                 foreach $frame ( @{ $options->{ "frames" } } )
2404                 {
2405                     %new_record = %{ $record };
2406
2407                     $new_record{ "SEQ" }     = Maasha::Seq::translate( $record->{ "SEQ" }, $frame );
2408                     $new_record{ "SEQ_LEN" } = length $record->{ "SEQ" };
2409                     $new_record{ "FRAME" }   = $frame;
2410
2411                     put_record( \%new_record, $out );
2412                 }
2413             }
2414         }
2415         else
2416         {
2417             put_record( $record, $out );
2418         }
2419     }
2420 }
2421
2422
2423 sub script_extract_seq
2424 {
2425     # Martin A. Hansen, August 2007.
2426
2427     # Extract subsequences from sequences in record.
2428
2429     my ( $in,        # handle to in stream
2430          $out,       # handle to out stream
2431          $options,   # options hash
2432        ) = @_;
2433
2434     # Returns nothing.
2435
2436     my ( $beg, $end, $len, $record );
2437
2438     if ( not defined $options->{ "beg" } or $options->{ "beg" } < 0 ) {
2439         $beg = 0;
2440     } else {
2441         $beg = $options->{ "beg" } - 1;   # correcting for start offset
2442     }
2443
2444     if ( defined $options->{ "end" } and $options->{ "end" } - 1 < $beg ) {
2445         $end = $beg - 1;
2446     } elsif ( defined $options->{ "end" } ) {
2447         $end = $options->{ "end" } - 1;   # correcting for start offset
2448     }
2449
2450     $len = $options->{ "len" };
2451
2452 #    print "beg->$beg,  end->$end,  len->$len\n";
2453
2454     while ( $record = get_record( $in ) )
2455     {
2456         if ( $record->{ "SEQ" } )
2457         {
2458             if ( defined $beg and defined $end )
2459             {
2460                 if ( $end - $beg + 1 > length $record->{ "SEQ" } ) {
2461                     $record->{ "SEQ" } = substr $record->{ "SEQ" }, $beg;
2462                 } else {
2463                     $record->{ "SEQ" } = substr $record->{ "SEQ" }, $beg, $end - $beg + 1;
2464                 }
2465             }
2466             elsif ( defined $beg and defined $len )
2467             {
2468                 if ( $len > length $record->{ "SEQ" } ) {
2469                     $record->{ "SEQ" } = substr $record->{ "SEQ" }, $beg;
2470                 } else {
2471                     $record->{ "SEQ" } = substr $record->{ "SEQ" }, $beg, $len;
2472                 }
2473             }
2474             elsif ( defined $beg )
2475             {
2476                 $record->{ "SEQ" } = substr $record->{ "SEQ" }, $beg;
2477             }
2478         }
2479
2480         put_record( $record, $out );
2481     }
2482 }
2483
2484
2485 sub script_get_genome_seq
2486 {
2487     # Martin A. Hansen, December 2007.
2488
2489     # Gets a subsequence from a genome.
2490
2491     my ( $in,        # handle to in stream
2492          $out,       # handle to out stream
2493          $options,   # options hash
2494        ) = @_;
2495
2496     # Returns nothing.
2497
2498     my ( $record, $genome_file, $index_file, $index, $fh, $index_head, $index_beg, $index_len, $beg, $len, %lookup_hash, @begs, @lens, $i );
2499
2500     $options->{ "flank" } ||= 0;
2501
2502     if ( $options->{ "genome" } )
2503     {
2504         $genome_file = Maasha::Config::genome_fasta( $options->{ 'genome' } );
2505         $index_file  = Maasha::Config::genome_fasta_index( $options->{ 'genome' } );
2506
2507         $fh          = Maasha::Common::read_open( $genome_file );
2508         $index       = Maasha::Fasta::index_retrieve( $index_file );
2509
2510         shift @{ $index }; # Get rid of the file size info
2511
2512         map { $lookup_hash{ $_->[ 0 ] } = [ $_->[ 1 ], $_->[ 2 ] ] } @{ $index };
2513
2514         if ( exists $lookup_hash{ $options->{ "chr" } } and defined $options->{ "beg" } and ( defined $options->{ "end" } or defined $options->{ "len" } ) )
2515         {
2516             ( $index_beg, $index_len ) = @{ $lookup_hash{ $options->{ "chr" } } };
2517
2518             $beg = $index_beg + $options->{ "beg" } - 1;
2519
2520             if ( $options->{ "len" } ) {
2521                 $len = $options->{ "len" };
2522             } elsif ( $options->{ "end" } ) {
2523                 $len = ( $options->{ "end" } - $options->{ "beg" } + 1 );
2524             }
2525
2526             $beg -= $options->{ "flank" };
2527             $len += 2 * $options->{ "flank" };
2528
2529             if ( $beg <= $index_beg )
2530             {
2531                 $len -= $index_beg - $beg;
2532                 $beg = $index_beg;
2533             }
2534
2535             $len = $index_beg + $index_len - $beg if $beg + $len > $index_beg + $index_len;
2536
2537             next if $beg > $index_beg + $index_len;
2538
2539             $record->{ "CHR" }     = $options->{ "chr" };
2540             $record->{ "CHR_BEG" } = $beg - $index_beg;
2541             $record->{ "CHR_END" } = $record->{ "CHR_BEG" } + $len - 1;
2542
2543             $record->{ "SEQ" }     = Maasha::Common::file_read( $fh, $beg, $len );
2544             $record->{ "SEQ_LEN" } = $len;
2545
2546             put_record( $record, $out );
2547         }
2548     }
2549
2550     while ( $record = get_record( $in ) )
2551     {
2552         if ( $options->{ "genome" } and not $record->{ "SEQ" } )
2553         {
2554             if ( $record->{ "REC_TYPE" } eq "BED" and exists $lookup_hash{ $record->{ "CHR" } } )
2555             {
2556                 ( $index_beg, $index_len ) = @{ $lookup_hash{ $record->{ "CHR" } } };
2557
2558                 $beg = $record->{ "CHR_BEG" } + $index_beg;
2559                 $len = $record->{ "CHR_END" } - $record->{ "CHR_BEG" } + 1;
2560             }
2561             elsif ( $record->{ "REC_TYPE" } eq "PSL" and exists $lookup_hash{ $record->{ "S_ID" } } )
2562             {
2563                 ( $index_beg, $index_len ) = @{ $lookup_hash{ $record->{ "S_ID" } } };
2564
2565                 $beg = $record->{ "S_BEG" } + $index_beg;
2566                 $len = $record->{ "S_END" } - $record->{ "S_BEG" } + 1;
2567             }
2568             elsif ( $record->{ "REC_TYPE" } eq "BLAST" and exists $lookup_hash{ $record->{ "S_ID" } } )
2569             {
2570                 ( $index_beg, $index_len ) = @{ $lookup_hash{ $record->{ "S_ID" } } };
2571
2572                 $beg = $record->{ "S_BEG" } + $index_beg;
2573                 $len = $record->{ "S_END" } - $record->{ "S_BEG" } + 1;
2574             }
2575
2576             $beg -= $options->{ "flank" };
2577             $len += 2 * $options->{ "flank" };
2578
2579             if ( $beg <= $index_beg )
2580             {
2581                 $len -= $index_beg - $beg;
2582                 $beg = $index_beg;
2583             }
2584
2585             $len = $index_beg + $index_len - $beg if $beg + $len > $index_beg + $index_len;
2586
2587             next if $beg > $index_beg + $index_len;
2588
2589             $record->{ "CHR_BEG" } = $beg - $index_beg;
2590             $record->{ "CHR_END" } = $record->{ "CHR_BEG" } + $len - 1;
2591
2592             $record->{ "SEQ" } = Maasha::Common::file_read( $fh, $beg, $len );
2593
2594             if ( $record->{ "STRAND" } and $record->{ "STRAND" } eq "-" )
2595             {
2596                 Maasha::Seq::dna_comp( \$record->{ "SEQ" } );
2597                 $record->{ "SEQ" } = reverse $record->{ "SEQ" };
2598             }
2599
2600             if ( $options->{ "mask" } )
2601             {
2602                 if ( $record->{ "BLOCKCOUNT" } > 1 ) # uppercase hit block segments and lowercase the rest.
2603                 {
2604                     $record->{ "SEQ" } = lc $record->{ "SEQ" };
2605
2606                     @begs = split ",", $record->{ "Q_BEGS" };
2607                     @lens = split ",", $record->{ "BLOCKSIZES" };
2608
2609                     for ( $i = 0; $i < @begs; $i++ ) {
2610                         substr $record->{ "SEQ" }, $begs[ $i ], $lens[ $i ], uc substr $record->{ "SEQ" }, $begs[ $i ], $lens[ $i ];
2611                     }
2612                 }
2613             }
2614         }
2615
2616         put_record( $record, $out );
2617     }
2618
2619     close $fh if $fh;
2620 }
2621
2622
2623 sub script_get_genome_align
2624 {
2625     # Martin A. Hansen, April 2008.
2626
2627     # Gets a subalignment from a multiple genome alignment.
2628
2629     my ( $in,        # handle to in stream
2630          $out,       # handle to out stream
2631          $options,   # options hash
2632        ) = @_;
2633
2634     # Returns nothing.
2635
2636     my ( $record, $maf_track, $align, $align_num, $beg, $end, $len, $entry );
2637
2638     $options->{ "strand" } ||= "+";
2639
2640     $align_num = 1;
2641
2642     $maf_track = Maasha::Config::maf_track( $options->{ "genome" } );
2643
2644     if ( $options->{ "chr" } and $options->{ "beg" } and ( $options->{ "end" } or $options->{ "len" } ) )
2645     {
2646         $beg = $options->{ "beg" } - 1;
2647
2648         if ( $options->{ "end" } ) {
2649             $end = $options->{ "end" };
2650         } elsif ( $options->{ "len" } ) {
2651             $end = $beg + $options->{ "len" };
2652         }
2653
2654         $align = Maasha::UCSC::maf_extract( $BP_TMP, $options->{ "genome" }, $maf_track, $options->{ "chr" }, $beg, $end, $options->{ "strand" } );
2655
2656         foreach $entry ( @{ $align } )
2657         {
2658             $entry->{ "CHR" }     = $record->{ "CHR" };
2659             $entry->{ "CHR_BEG" } = $record->{ "CHR_BEG" };
2660             $entry->{ "CHR_END" } = $record->{ "CHR_END" };
2661             $entry->{ "STRAND" }  = $record->{ "STRAND" } || '+';
2662             $entry->{ "Q_ID" }    = $record->{ "Q_ID" };
2663             $entry->{ "SCORE" }   = $record->{ "SCORE" };
2664
2665             put_record( $entry, $out );
2666         }
2667     }
2668
2669     while ( $record = get_record( $in ) )
2670     {
2671         if ( $record->{ "REC_TYPE" } eq "BED" )
2672         {
2673             $align = Maasha::UCSC::maf_extract( $BP_TMP, $options->{ "genome" }, $maf_track, $record->{ "CHR" }, $record->{ "CHR_BEG" }, $record->{ "CHR_END" }, $record->{ "STRAND" } );
2674         }
2675         elsif ( $record->{ "REC_TYPE" } eq "PSL" )
2676         {
2677             $align = Maasha::UCSC::maf_extract( $BP_TMP, $options->{ "genome" }, $maf_track, $record->{ "S_ID" }, $record->{ "S_BEG" }, $record->{ "S_END" }, $record->{ "STRAND" } );
2678         }
2679         elsif ( $record->{ "REC_TYPE" } eq "BLAST" )
2680         {
2681             $align = Maasha::UCSC::maf_extract( $BP_TMP, $options->{ "genome" }, $maf_track, $record->{ "S_ID" }, $record->{ "S_BEG" }, $record->{ "S_END" }, $record->{ "STRAND" } );
2682         }
2683
2684         foreach $entry ( @{ $align } )
2685         {
2686             $entry->{ "CHR" }     = $record->{ "CHR" };
2687             $entry->{ "CHR_BEG" } = $record->{ "CHR_BEG" };
2688             $entry->{ "CHR_END" } = $record->{ "CHR_END" };
2689             $entry->{ "STRAND" }  = $record->{ "STRAND" };
2690             $entry->{ "Q_ID" }    = $record->{ "Q_ID" };
2691             $entry->{ "SCORE" }   = $record->{ "SCORE" };
2692
2693             put_record( $entry, $out );
2694         }
2695
2696         $align_num++;
2697     }
2698 }
2699
2700
2701 sub script_get_genome_phastcons
2702 {
2703     # Martin A. Hansen, February 2008.
2704
2705     # Get phastcons scores from genome intervals.
2706
2707     my ( $in,        # handle to in stream
2708          $out,       # handle to out stream
2709          $options,   # options hash
2710        ) = @_;
2711
2712     # Returns nothing.
2713
2714     my ( $phastcons_file, $phastcons_index, $index, $fh_phastcons, $scores, $record );
2715
2716     $options->{ "flank" } ||= 0;
2717
2718     $phastcons_file  = Maasha::Config::genome_phastcons( $options->{ "genome" } );
2719     $phastcons_index = Maasha::Config::genome_phastcons_index( $options->{ "genome" } );
2720
2721     $index           = Maasha::UCSC::phastcons_index_retrieve( $phastcons_index );
2722     $fh_phastcons    = Maasha::Common::read_open( $phastcons_file );
2723
2724     if ( defined $options->{ "chr" } and defined $options->{ "beg" } and ( defined $options->{ "end" } or defined $options->{ "len" } ) )
2725     {
2726         $options->{ "beg" } -= 1;   # request is 1-based
2727         $options->{ "end" } -= 1;   # request is 1-based
2728
2729         if ( $options->{ "len" } ) {
2730             $options->{ "end" } = $options->{ "beg" } + $options->{ "len" } - 1;
2731         }
2732
2733         $scores = Maasha::UCSC::phastcons_index_lookup( $index, $fh_phastcons, $options->{ "chr" }, $options->{ "beg" }, $options->{ "end" }, $options->{ "flank" } );
2734
2735         $record->{ "CHR" }       = $options->{ "chr" };
2736         $record->{ "CHR_BEG" }   = $options->{ "beg" } - $options->{ "flank" };
2737         $record->{ "CHR_END" }   = $options->{ "end" } + $options->{ "flank" };
2738
2739         $record->{ "PHASTCONS" }   = join ",", @{ $scores };
2740         $record->{ "PHAST_COUNT" } = scalar @{ $scores };  # DEBUG
2741
2742         put_record( $record, $out );
2743     }
2744
2745     while ( $record = get_record( $in ) )
2746     {
2747         if ( $record->{ "REC_TYPE" } eq "BED" )
2748         {
2749             $scores = Maasha::UCSC::phastcons_index_lookup( $index, $fh_phastcons, $record->{ "CHR" }, $record->{ "CHR_BEG" }, $record->{ "CHR_END" }, $options->{ "flank" } );
2750         }
2751         elsif ( $record->{ "REC_TYPE" } eq "PSL" )
2752         {
2753             $scores = Maasha::UCSC::phastcons_index_lookup( $index, $fh_phastcons, $record->{ "S_ID" }, $record->{ "S_BEG" }, $record->{ "S_END" }, $options->{ "flank" } );
2754         }
2755         elsif ( $record->{ "REC_TYPE" } eq "BLAST" )
2756         {
2757             $scores = Maasha::UCSC::phastcons_index_lookup( $index, $fh_phastcons, $record->{ "S_ID" }, $record->{ "S_BEG" }, $record->{ "S_END" }, $options->{ "flank" } );
2758         }
2759
2760         $record->{ "PHASTCONS" } = join ",", @{ $scores } if @{ $scores };
2761 #        $record->{ "PHAST_COUNT" } = @{ $scores } if @{ $scores };  # DEBUG
2762
2763         put_record( $record, $out );
2764     }
2765
2766     close $fh_phastcons if $fh_phastcons;
2767 }
2768
2769
2770 sub script_fold_seq
2771 {
2772     # Martin A. Hansen, December 2007.
2773
2774     # Folds sequences in stream into secondary structures.
2775
2776     my ( $in,     # handle to in stream
2777          $out,    # handle to out stream
2778        ) = @_;
2779
2780     # Returns nothing.
2781
2782     my ( $record, $type, $struct, $index );
2783
2784     while ( $record = get_record( $in ) )
2785     {
2786         if ( $record->{ "SEQ" } )
2787         {
2788             if ( not $type ) {
2789                 $type = Maasha::Seq::seq_guess_type( $record->{ "SEQ" } );
2790             }
2791
2792             if ( $type ne "protein" )
2793             {
2794                 ( $struct, $index ) = Maasha::Seq::fold_struct_rnafold( $record->{ "SEQ" } );
2795                 $record->{ "SEC_STRUCT" }  = $struct;
2796                 $record->{ "FREE_ENERGY" } = $index;
2797                 $record->{ "SCORE" }       = abs int $index;
2798                 $record->{ "SIZE" }        = length $struct;
2799                 $record->{ "CONF" }        = "1," x $record->{ "SIZE" };
2800             }
2801         }
2802
2803         put_record( $record, $out );
2804     }
2805 }
2806
2807
2808 sub script_split_seq
2809 {
2810     # Martin A. Hansen, August 2007.
2811
2812     # Split a sequence in stream into words.
2813
2814     my ( $in,        # handle to in stream
2815          $out,       # handle to out stream
2816          $options,   # options hash
2817        ) = @_;
2818
2819     # Returns nothing.
2820
2821     my ( $record, $new_record, $i, $subseq, %lookup );
2822
2823     $options->{ "word_size" } ||= 7;
2824
2825     while ( $record = get_record( $in ) )
2826     {
2827         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } )
2828         {
2829             for ( $i = 0; $i < length( $record->{ "SEQ" } ) - $options->{ "word_size" } + 1; $i++ )
2830             {
2831                 $subseq = substr $record->{ "SEQ" }, $i, $options->{ "word_size" };
2832
2833                 if ( $options->{ "uniq" } and not $lookup{ $subseq } )
2834                 {
2835                     $new_record->{ "REC_TYPE" } = "SPLIT";
2836                     $new_record->{ "SEQ_NAME" } = $record->{ "SEQ_NAME" } . "[" . ( $i + 1 ) . "-" . ( $i + $options->{ "word_size" } ) . "]";
2837                     $new_record->{ "SEQ" }      = $subseq;
2838
2839                     put_record( $new_record, $out );
2840
2841                     $lookup{ $subseq } = 1;
2842                 }
2843                 else
2844                 {
2845                     $new_record->{ "REC_TYPE" } = "SPLIT";
2846                     $new_record->{ "SEQ_NAME" } = $record->{ "SEQ_NAME" } . "[" . ( $i + 1 ) . "-" . ( $i + $options->{ "word_size" } ) . "]";
2847                     $new_record->{ "SEQ" }      = $subseq;
2848
2849                     put_record( $new_record, $out );
2850                 }
2851             }
2852         }
2853         else
2854         {
2855             put_record( $record, $out );
2856         }
2857     }
2858 }
2859
2860
2861 sub script_split_bed
2862 {
2863     # Martin A. Hansen, June 2008.
2864
2865     # Split a BED record into overlapping windows.
2866
2867     my ( $in,        # handle to in stream
2868          $out,       # handle to out stream
2869          $options,   # options hash
2870        ) = @_;
2871
2872     # Returns nothing.
2873
2874     my ( $record, $new_record, $i );
2875
2876     $options->{ "window_size" } ||= 20;
2877     $options->{ "step_size" }   ||= 1;
2878
2879     while ( $record = get_record( $in ) )
2880     {
2881         if ( $record->{ "CHR" } and $record->{ "CHR_BEG" } and $record->{ "CHR_END" } )
2882         {
2883             $record->{ "BED_LEN" } = $record->{ "CHR_END" } - $record->{ "CHR_BEG" } + 1;
2884
2885             for ( $i = 0; $i < $record->{ "BED_LEN" } - $options->{ "window_size" }; $i += $options->{ "step_size" } )
2886             {
2887                 $new_record->{ "REC_TYPE" } = "BED";
2888                 $new_record->{ "CHR" }      = $record->{ "CHR" };
2889                 $new_record->{ "CHR_BEG" }  = $record->{ "CHR_BEG" } + $i;
2890                 $new_record->{ "CHR_END" }  = $record->{ "CHR_BEG" } + $i + $options->{ "window_size" };
2891                 $new_record->{ "BED_LEN" }  = $options->{ "window_size" };
2892                 $new_record->{ "Q_ID" }     = $record->{ "Q_ID" } . "_$i";
2893                 $new_record->{ "SCORE" }    = $record->{ "SCORE" };
2894                 $new_record->{ "STRAND" }   = $record->{ "STRAND" };
2895
2896                 put_record( $new_record, $out );
2897             }
2898         }
2899         else
2900         {
2901             put_record( $record, $out );
2902         }
2903     }
2904 }
2905
2906
2907 sub script_align_seq
2908 {
2909     # Martin A. Hansen, August 2007.
2910
2911     # Align sequences in stream.
2912
2913     my ( $in,    # handle to in stream
2914          $out,   # handle to out stream
2915        ) = @_;
2916
2917     # Returns nothing.
2918
2919     my ( $record, @entries, $entry );
2920
2921     while ( $record = get_record( $in ) )
2922     {
2923         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } ) {
2924             push @entries, [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ];
2925         } elsif ( $record->{ "Q_ID" } and $record->{ "SEQ" } ) {
2926             push @entries, [ $record->{ "Q_ID" }, $record->{ "SEQ" } ];
2927         } else {
2928             put_record( $record, $out );
2929         }
2930     }
2931
2932     @entries = Maasha::Align::align( \@entries );
2933
2934     foreach $entry ( @entries )
2935     {
2936         if ( $entry->[ SEQ_NAME ] and $entry->[ SEQ ] )
2937         {
2938             $record = {
2939                 SEQ_NAME => $entry->[ SEQ_NAME ],
2940                 SEQ      => $entry->[ SEQ ],
2941             };
2942
2943             put_record( $record, $out );
2944         }
2945     }
2946 }
2947
2948
2949 sub script_tile_seq
2950 {
2951     # Martin A. Hansen, February 2008.
2952
2953     # Using the first sequence in stream as reference, tile
2954     # all subsequent sequences based on pairwise alignments.
2955
2956     my ( $in,        # handle to in stream
2957          $out,       # handle to out stream
2958          $options,   # options hash
2959        ) = @_;
2960
2961     # Returns nothing.
2962
2963     my ( $record, $first, $ref_entry, @entries );
2964
2965     $first = 1;
2966
2967     while ( $record = get_record( $in ) )
2968     {
2969         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } )
2970         {
2971             if ( $first )
2972             {
2973                 $ref_entry = [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ];
2974
2975                 $first = 0;
2976             }
2977             else
2978             {
2979                 push @entries, [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ];
2980             }
2981         }
2982         else
2983         {
2984             put_record( $record, $out );
2985         }
2986     }
2987
2988     @entries = Maasha::Align::align_tile( $ref_entry, \@entries, $options );
2989
2990     map { put_record( { SEQ_NAME => $_->[ SEQ_NAME ], SEQ => $_->[ SEQ ] }, $out ) } @entries;
2991 }
2992
2993
2994 sub script_invert_align
2995 {
2996     # Martin A. Hansen, February 2008.
2997
2998     # Inverts an alignment showing only non-mathing residues
2999     # using the first sequence as reference.
3000
3001     my ( $in,        # handle to in stream
3002          $out,       # handle to out stream
3003          $options,   # options hash
3004        ) = @_;
3005
3006     # Returns nothing.
3007
3008     my ( $record, @entries );
3009
3010     while ( $record = get_record( $in ) )
3011     {
3012         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } )
3013         {
3014             push @entries, [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ];
3015         }
3016         else
3017         {
3018             put_record( $record, $out );
3019         }
3020     }
3021
3022     Maasha::Align::align_invert( \@entries, $options->{ "soft" } );
3023
3024     map { put_record( { SEQ_NAME => $_->[ SEQ_NAME ], SEQ => $_->[ SEQ ] }, $out ) } @entries;
3025 }
3026
3027
3028 sub script_patscan_seq
3029 {
3030     # Martin A. Hansen, August 2007.
3031
3032     # Locates patterns in sequences using scan_for_matches.
3033
3034     my ( $in,        # handle to in stream
3035          $out,       # handle to out stream
3036          $options,   # options hash
3037        ) = @_;
3038
3039     # Returns nothing.
3040
3041     my ( $genome_file, @args, $arg, $type, $seq_file, $pat_file, $out_file, $fh_in, $fh_out, $record, $patterns, $pattern, $entry, $result, %head_hash, $i );
3042
3043     if ( $options->{ "patterns" } ) {
3044         $patterns = Maasha::Patscan::parse_patterns( $options->{ "patterns" } );
3045     } elsif ( -f $options->{ "patterns_in" } ) {
3046         $patterns = Maasha::Patscan::read_patterns( $options->{ "patterns_in" } );
3047     }
3048
3049     $genome_file = Maasha::Config::genome_fasta( $options->{ 'genome' } ) if $options->{ 'genome' };
3050
3051     push @args, "-c"                            if $options->{ "comp" };
3052     push @args, "-m $options->{ 'max_hits' }"   if $options->{ 'max_hits' };
3053     push @args, "-n $options->{ 'max_misses' }" if $options->{ 'max_hits' };
3054
3055     $seq_file = "$BP_TMP/patscan.seq";
3056     $pat_file = "$BP_TMP/patscan.pat";
3057     $out_file = "$BP_TMP/patscan.out";
3058
3059     $fh_out = Maasha::Common::write_open( $seq_file );
3060
3061     $i = 0;
3062
3063     while ( $record = get_record( $in ) )
3064     {
3065         if ( $record->{ "SEQ" } and $record->{ "SEQ_NAME" } )
3066         {
3067             $type = Maasha::Seq::seq_guess_type( $record->{ "SEQ" } ) if not $type;
3068
3069             Maasha::Fasta::put_entry( [ $i, $record->{ "SEQ" } ], $fh_out );
3070
3071             $head_hash{ $i } = $record->{ "SEQ_NAME" };
3072
3073             $i++;
3074         }
3075
3076 #        put_record( $record, $out );
3077     }
3078
3079     close $fh_out;
3080
3081     $arg  = join " ", @args;
3082     $arg .= " -p" if $type eq "protein";
3083
3084     foreach $pattern ( @{ $patterns } )
3085     {
3086         $fh_out = Maasha::Common::write_open( $pat_file );
3087
3088         print $fh_out "$pattern\n";
3089
3090         close $fh_out;
3091
3092         if ( $options->{ 'genome' } ) {
3093             `scan_for_matches $arg $pat_file < $genome_file > $out_file`;
3094             # Maasha::Common::run( "scan_for_matches", "$arg $pat_file < $genome_file > $out_file" );
3095         } else {
3096             `scan_for_matches $arg $pat_file < $seq_file > $out_file`;
3097             # Maasha::Common::run( "scan_for_matches", "$arg $pat_file < $seq_file > $out_file" );
3098         }
3099
3100         $fh_in = Maasha::Common::read_open( $out_file );
3101
3102         while ( $entry = Maasha::Fasta::get_entry( $fh_in ) )
3103         {
3104             $result = Maasha::Patscan::parse_scan_result( $entry, $pattern );
3105
3106             if ( $options->{ 'genome' } )
3107             {
3108                 $result->{ "CHR" }     = $result->{ "S_ID" };
3109                 $result->{ "CHR_BEG" } = $result->{ "S_BEG" };
3110                 $result->{ "CHR_END" } = $result->{ "S_END" };
3111
3112                 delete $result->{ "S_ID" };
3113                 delete $result->{ "S_BEG" };
3114                 delete $result->{ "S_END" };
3115             }
3116             else
3117             {
3118                 $result->{ "S_ID" } = $head_hash{ $result->{ "S_ID" } };
3119             }
3120
3121             put_record( $result, $out );
3122         }
3123
3124         close $fh_in;
3125     }
3126
3127     unlink $pat_file;
3128     unlink $seq_file;
3129     unlink $out_file;
3130 }
3131
3132
3133 sub script_create_blast_db
3134 {
3135     # Martin A. Hansen, September 2007.
3136
3137     # Creates a NCBI BLAST database with formatdb
3138
3139     my ( $in,        # handle to in stream
3140          $out,       # handle to out stream
3141          $options,   # options hash
3142        ) = @_;
3143
3144     # Returns nothing.
3145
3146     my ( $fh, $seq_type, $path, $record );
3147
3148     $path = $options->{ "database" };
3149
3150     $fh = Maasha::Common::write_open( $path );
3151
3152     while ( $record = get_record( $in ) )
3153     {
3154         put_record( $record, $out ) if not $options->{ "no_stream" };
3155
3156         if ( $record->{ "SEQ" } and $record->{ "SEQ_NAME" } )
3157         {
3158             $seq_type = Maasha::Seq::seq_guess_type( $record->{ "SEQ" } ) if not $seq_type;
3159
3160             Maasha::Fasta::put_entry( [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ], $fh );
3161         }
3162     }
3163
3164     close $fh;
3165
3166     if ( $seq_type eq "protein" ) {
3167         Maasha::Common::run( "formatdb", "-p T -i $path -t $options->{ 'database' }" );
3168     } else {
3169         Maasha::Common::run( "formatdb", "-p F -i $path -t $options->{ 'database' }" );
3170     }
3171
3172     unlink $path;
3173 }
3174
3175
3176 sub script_blast_seq
3177 {
3178     # Martin A. Hansen, September 2007.
3179
3180     # BLASTs sequences in stream against a given database.
3181
3182     my ( $in,        # handle to in stream
3183          $out,       # handle to out stream
3184          $options,   # options hash
3185        ) = @_;
3186
3187     # Returns nothing.
3188
3189     my ( $genome, $q_type, $s_type, $tmp_in, $tmp_out, $fh_in, $fh_out, $record, $line, @fields );
3190
3191     $options->{ "e_val" }  = 10 if not defined $options->{ "e_val" };
3192     $options->{ "filter" } = "F";
3193     $options->{ "filter" } = "T" if $options->{ "filter" };
3194     $options->{ "cpus" } ||= 1;
3195
3196     $options->{ "database" } = Maasha::Config::genome_blast( $options->{ 'genome' } ) if $options->{ 'genome' };
3197
3198     $tmp_in  = "$BP_TMP/blast_query.seq";
3199     $tmp_out = "$BP_TMP/blast.result";
3200
3201     $fh_out = Maasha::Common::write_open( $tmp_in );
3202
3203     while ( $record = get_record( $in ) )
3204     {
3205         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } )
3206         {
3207             $q_type = Maasha::Seq::seq_guess_type( $record->{ "SEQ" } ) if not $q_type;
3208
3209             Maasha::Fasta::put_entry( [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ], $fh_out );
3210         }
3211
3212         put_record( $record, $out );
3213     }
3214
3215     close $fh_out;
3216
3217     if ( -f $options->{ 'database' } . ".phr" ) {
3218         $s_type = "protein";
3219     } else {
3220         $s_type = "nucleotide";
3221     }
3222
3223     if ( not $options->{ 'program' } )
3224     {
3225         if ( $q_type ne "protein" and $s_type ne "protein" ) {
3226             $options->{ 'program' } = "blastn";
3227         } elsif ( $q_type eq "protein" and $s_type eq "protein" ) {
3228             $options->{ 'program' } = "blastp";
3229         } elsif ( $q_type ne "protein" and $s_type eq "protein" ) {
3230             $options->{ 'program' } = "blastx";
3231         } elsif ( $q_type eq "protein" and $s_type ne "protein" ) {
3232             $options->{ 'program' } = "tblastn";
3233         }
3234     }
3235
3236     Maasha::Common::run( "blastall", "-p $options->{ 'program' } -e $options->{ 'e_val' } -a $options->{ 'cpus' } -m 8 -i $tmp_in -d $options->{ 'database' } -F $options->{ 'filter' } -o $tmp_out > /dev/null 2>&1", 1 );
3237
3238     unlink $tmp_in;
3239
3240     $fh_out = Maasha::Common::read_open( $tmp_out );
3241
3242     undef $record;
3243
3244     while ( $line = <$fh_out> )
3245     {
3246         chomp $line;
3247
3248         next if $line =~ /^#/;
3249
3250         @fields = split /\s+/, $line;
3251
3252         $record->{ "REC_TYPE" }   = "BLAST";
3253         $record->{ "Q_ID" }       = $fields[ 0 ];
3254         $record->{ "S_ID" }       = $fields[ 1 ];
3255         $record->{ "IDENT" }      = $fields[ 2 ];
3256         $record->{ "ALIGN_LEN" }  = $fields[ 3 ];
3257         $record->{ "MISMATCHES" } = $fields[ 4 ];
3258         $record->{ "GAPS" }       = $fields[ 5 ];
3259         $record->{ "Q_BEG" }      = $fields[ 6 ] - 1; # BLAST is 1-based
3260         $record->{ "Q_END" }      = $fields[ 7 ] - 1; # BLAST is 1-based
3261         $record->{ "S_BEG" }      = $fields[ 8 ] - 1; # BLAST is 1-based
3262         $record->{ "S_END" }      = $fields[ 9 ] - 1; # BLAST is 1-based
3263         $record->{ "E_VAL" }      = $fields[ 10 ];
3264         $record->{ "BIT_SCORE" }  = $fields[ 11 ];
3265
3266         if ( $record->{ "S_BEG" } > $record->{ "S_END" } )
3267         {
3268             $record->{ "STRAND" } = '-';
3269
3270             ( $record->{ "S_BEG" }, $record->{ "S_END" } ) = ( $record->{ "S_END" }, $record->{ "S_BEG" } );
3271         }
3272         else
3273         {
3274             $record->{ "STRAND" } = '+';
3275         }
3276
3277         put_record( $record, $out );
3278     }
3279
3280     close $fh_out;
3281
3282     unlink $tmp_out;
3283 }
3284
3285
3286 sub script_blat_seq
3287 {
3288     # Martin A. Hansen, August 2007.
3289
3290     # BLATs sequences in stream against a given genome.
3291
3292     my ( $in,        # handle to in stream
3293          $out,       # handle to out stream
3294          $options,   # options hash
3295        ) = @_;
3296
3297     # Returns nothing.
3298
3299     my ( $blat_args, $genome_file, $query_file, $fh_in, $fh_out, $type, $record, $result_file, $entries );
3300
3301     $genome_file = Maasha::Config::genome_fasta( $options->{ "genome" } );
3302
3303     $options->{ 'tile_size' }    ||= 11;
3304     $options->{ 'one_off' }      ||= 0;
3305     $options->{ 'min_identity' } ||= 90;
3306     $options->{ 'min_score' }    ||= 0;
3307     $options->{ 'step_size' }    ||= $options->{ 'tile_size' };
3308
3309     $blat_args .= " -tileSize=$options->{ 'tile_size' }";
3310     $blat_args .= " -oneOff=$options->{ 'one_off' }";
3311     $blat_args .= " -minIdentity=$options->{ 'min_identity' }";
3312     $blat_args .= " -minScore=$options->{ 'min_score' }";
3313     $blat_args .= " -stepSize=$options->{ 'step_size' }";
3314     $blat_args .= " -ooc=" . Maasha::Config::genome_blat_ooc( $options->{ "genome" }, 11 ) if $options->{ 'ooc' };
3315
3316     $query_file = "$BP_TMP/blat.seq";
3317
3318     $fh_out = Maasha::Common::write_open( $query_file );
3319
3320     while ( $record = get_record( $in ) )
3321     {
3322         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } )
3323         {
3324             Maasha::Fasta::put_entry( [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ], $fh_out, 80 );
3325             $type = Maasha::Seq::seq_guess_type( $record->{ "SEQ" } ) if not $type;
3326         }
3327
3328         put_record( $record, $out );
3329     }
3330
3331     close $fh_out;
3332
3333     $blat_args .= " -t=dnax" if $type eq "protein";
3334     $blat_args .= " -q=$type";
3335
3336     $result_file = "$BP_TMP/blat.psl";
3337
3338     Maasha::Common::run( "blat", "$genome_file $query_file $blat_args $result_file > /dev/null 2>&1" );
3339
3340     unlink $query_file;
3341
3342     $entries = Maasha::UCSC::psl_get_entries( $result_file );
3343
3344     map { put_record( $_, $out ) } @{ $entries };
3345
3346     unlink $result_file;
3347 }
3348
3349
3350 sub script_match_seq
3351 {
3352     # Martin A. Hansen, August 2007.
3353
3354     # BLATs sequences in stream against a given genome.
3355
3356     my ( $in,        # handle to in stream
3357          $out,       # handle to out stream
3358          $options,   # options hash
3359        ) = @_;
3360
3361     # Returns nothing.
3362
3363     my ( $record, @entries, $results );
3364
3365     $options->{ "word_size" } ||= 20;
3366     $options->{ "direction" } ||= "both";
3367
3368     while ( $record = get_record( $in ) )
3369     {
3370         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } ) {
3371             push @entries, [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ];
3372         }
3373
3374         put_record( $record, $out );
3375     }
3376
3377     if ( @entries == 1 )
3378     {
3379         $results = Maasha::Match::match_mummer( [ $entries[ 0 ] ], [ $entries[ 0 ] ], $options, $BP_TMP );
3380
3381         map { put_record( $_, $out ) } @{ $results };
3382     }
3383     elsif ( @entries == 2 )
3384     {
3385         $results = Maasha::Match::match_mummer( [ $entries[ 0 ] ], [ $entries[ 1 ] ], $options, $BP_TMP );
3386
3387         map { put_record( $_, $out ) } @{ $results };
3388     }
3389 }
3390
3391
3392 sub script_create_vmatch_index
3393 {
3394     # Martin A. Hansen, January 2008.
3395
3396     # Create a vmatch index from sequences in the stream.
3397
3398     my ( $in,        # handle to in stream
3399          $out,       # handle to out stream
3400          $options,   # options hash
3401        ) = @_;
3402
3403     # Returns nothing.
3404
3405     my ( $record, $file_tmp, $fh_tmp, $type );
3406
3407     if ( $options->{ "index_name" } )
3408     {
3409         $file_tmp = $options->{ 'index_name' };
3410         $fh_tmp   = Maasha::Common::write_open( $file_tmp );
3411     }
3412
3413     while ( $record = get_record( $in ) )
3414     {
3415         if ( $options->{ "index_name" } and $record->{ "SEQ_NAME" } and $record->{ "SEQ" } )
3416         {
3417             Maasha::Fasta::put_entry( [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ], $fh_tmp );
3418
3419             $type = Maasha::Seq::seq_guess_type( $record->{ "SEQ" } ) if not $type;
3420         }
3421
3422         put_record( $record, $out ) if not $options->{ "no_stream" };
3423     }
3424
3425     if ( $options->{ "index_name" } )
3426     {
3427         close $fh_tmp;
3428
3429         if ( $type eq "protein" ) {
3430             Maasha::Common::run( "mkvtree", "-db $file_tmp -protein -pl $options->{ 'prefix_length' } -allout -indexname $file_tmp > /dev/null 2>&1" );
3431         } else {
3432             Maasha::Common::run( "mkvtree", "-db $file_tmp -dna -pl $options->{ 'prefix_length' } -allout -indexname $file_tmp > /dev/null 2>&1" );
3433         }
3434
3435         unlink $file_tmp;
3436     }
3437 }
3438
3439
3440 sub script_vmatch_seq
3441 {
3442     # Martin A. Hansen, August 2007.
3443
3444     # Vmatches sequences in stream against a given genome.
3445
3446     my ( $in,        # handle to in stream
3447          $out,       # handle to out stream
3448          $options,   # options hash
3449        ) = @_;
3450
3451     # Returns nothing.
3452
3453     my ( @index_files, @records, $result_file, $fh_in, $record );
3454
3455     $options->{ 'count' } = 1 if $options->{ 'max_hits' };
3456
3457     if ( $options->{ "index_name" } ) {
3458         @index_files = $options->{ "index_name" };
3459     } else {
3460         @index_files = Maasha::Config::genome_vmatch( $options->{ "genome" } );
3461     }
3462
3463     while ( $record = get_record( $in ) )
3464     {
3465         push @records, $record;
3466
3467         put_record( $record, $out );
3468     }
3469
3470     $result_file = Maasha::Match::match_vmatch( $BP_TMP, \@records, \@index_files, $options );
3471
3472     undef @records;
3473
3474     $fh_in = Maasha::Common::read_open( $result_file );
3475
3476     while ( $record = Maasha::Match::vmatch_get_entry( $fh_in ) ) {
3477         put_record( $record, $out );
3478     }
3479
3480     close $fh_in;
3481
3482     unlink $result_file;
3483 }
3484
3485
3486 sub script_write_fasta
3487 {
3488     # Martin A. Hansen, August 2007.
3489
3490     # Write FASTA entries from sequences in stream.
3491
3492     my ( $in,        # handle to in stream
3493          $out,       # handle to out stream
3494          $options,   # options hash
3495        ) = @_;
3496
3497     # Returns nothing.
3498
3499     my ( $record, $fh );
3500
3501     $fh = write_stream( $options->{ "data_out" }, $options->{ "compress" } );
3502
3503     while ( $record = get_record( $in ) )
3504     {
3505         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } ) {
3506             Maasha::Fasta::put_entry( [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ], $fh, $options->{ "wrap" } );
3507         }
3508
3509         put_record( $record, $out ) if not $options->{ "no_stream" };
3510     }
3511
3512     close $fh;
3513 }
3514
3515
3516 sub script_write_align
3517 {
3518     # Martin A. Hansen, August 2007.
3519
3520     # Write pretty alignments aligned sequences in stream.
3521
3522     my ( $in,        # handle to in stream
3523          $out,       # handle to out stream
3524          $options,   # options hash
3525        ) = @_;
3526
3527     # Returns nothing.
3528
3529     my ( $fh, $record, @entries );
3530
3531     $fh = write_stream( $options->{ "data_out" } ) ;
3532
3533     while ( $record = get_record( $in ) )
3534     {
3535         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } ) {
3536             push @entries, [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ];
3537         }
3538
3539         put_record( $record, $out ) if not $options->{ "no_stream" };
3540     }
3541
3542     if ( scalar( @entries ) == 2 ) {
3543         Maasha::Align::align_print_pairwise( $entries[ 0 ], $entries[ 1 ], $fh, $options->{ "wrap" } );
3544     } elsif ( scalar ( @entries ) > 2 ) {
3545         Maasha::Align::align_print_multi( \@entries, $fh, $options->{ "wrap" }, $options->{ "no_ruler" }, $options->{ "no_consensus" } );
3546     }
3547
3548     close $fh if $fh;
3549 }
3550
3551
3552 sub script_write_blast
3553 {
3554     # Martin A. Hansen, November 2007.
3555
3556     # Write data in blast table format (-m8 and 9).
3557
3558     my ( $in,        # handle to in stream
3559          $out,       # handle to out stream
3560          $options,   # options hash
3561        ) = @_;
3562
3563     # Returns nothing.
3564
3565     my ( $fh, $record, $first );
3566
3567     $fh = write_stream( $options->{ "data_out" }, $options->{ "compress" } ) ;
3568
3569     $first = 1;
3570
3571     while ( $record = get_record( $in ) )
3572     {
3573         if ( $record->{ "REC_TYPE" } eq "BLAST" )
3574         {
3575             if ( $options->{ "comment" } and $first )
3576             {
3577                 print "# Fields: Query id, Subject id, % identity, alignment length, mismatches, gap openings, q. start, q. end, s. start, s. end, e-value, bit score\n";
3578
3579                 $first = 0;
3580             }
3581
3582             if ( $record->{ "STRAND" } eq "-" ) {
3583                 ( $record->{ "S_BEG" }, $record->{ "S_END" } ) = ( $record->{ "S_END" }, $record->{ "S_BEG" } );
3584             }
3585
3586             print $fh join( "\t",
3587                 $record->{ "Q_ID" },
3588                 $record->{ "S_ID" },
3589                 $record->{ "IDENT" },
3590                 $record->{ "ALIGN_LEN" },
3591                 $record->{ "MISMATCHES" },
3592                 $record->{ "GAPS" },
3593                 $record->{ "Q_BEG" } + 1,
3594                 $record->{ "Q_END" } + 1,
3595                 $record->{ "S_BEG" } + 1,
3596                 $record->{ "S_END" } + 1,
3597                 $record->{ "E_VAL" },
3598                 $record->{ "BIT_SCORE" }
3599             ), "\n";
3600         }
3601
3602         put_record( $record, $out ) if not $options->{ "no_stream" };
3603     }
3604
3605     close $fh;
3606 }
3607
3608
3609 sub script_write_tab
3610 {
3611     # Martin A. Hansen, August 2007.
3612
3613     # Write data as table.
3614
3615     my ( $in,        # handle to in stream
3616          $out,       # handle to out stream
3617          $options,   # options hash
3618        ) = @_;
3619
3620     # Returns nothing.
3621
3622     my ( $fh, $record, $key, @keys, @vals, $ok, %no_keys, $A, $B );
3623
3624     $options->{ "delimit" } ||= "\t";
3625
3626     map { $no_keys{ $_ } = 1 } @{ $options->{ "no_keys" } };
3627
3628     $fh = write_stream( $options->{ "data_out" }, $options->{ "compress" } );
3629
3630     while ( $record = get_record( $in ) )
3631     {
3632         undef @vals;
3633         $ok = 1;
3634
3635         if ( $options->{ "keys" } )
3636         {
3637             map { $ok = 0 if not exists $record->{ $_ } } @{ $options->{ "keys" } };
3638
3639             if ( $ok )
3640             {
3641                 foreach $key ( @{ $options->{ "keys" }  } )
3642                 {
3643                     if ( exists $record->{ $key } )
3644                     {
3645                         push @keys, $key if $options->{ "comment" };
3646                         push @vals, $record->{ $key };
3647                     }
3648                 }
3649              }
3650         }
3651         else
3652         {
3653             foreach $key ( sort { $A = $a; $B = $b; $A =~ s/^V(\d+)$/$1/; $B =~ s/^V(\d+)$/$1/; $A <=> $B } keys %{ $record } )
3654             {
3655                 next if exists $no_keys{ $key };
3656
3657                 push @keys, $key if $options->{ "comment" };
3658                 push @vals, $record->{ $key };
3659             }
3660         }
3661
3662         if ( @keys and $options->{ "comment" } )
3663         {
3664             print $fh "#", join( $options->{ "delimit" }, @keys ), "\n";
3665
3666             delete $options->{ "comment" };
3667         }
3668
3669         print $fh join( $options->{ "delimit" }, @vals ), "\n" if @vals;
3670
3671         put_record( $record, $out ) if not $options->{ "no_stream" };
3672     }
3673
3674     close $fh;
3675 }
3676
3677
3678 sub script_write_bed
3679 {
3680     # Martin A. Hansen, August 2007.
3681
3682     # Write BED format for the UCSC genome browser using records in stream.
3683
3684     # Crude - needs lots of work!
3685
3686     my ( $in,        # handle to in stream
3687          $out,       # handle to out stream
3688          $options,   # options hash
3689        ) = @_;
3690
3691     # Returns nothing.
3692
3693     my ( $fh, $record, $new_record );
3694
3695     $fh = write_stream( $options->{ "data_out" }, $options->{ "compress" } );
3696
3697     while ( $record = get_record( $in ) )
3698     {
3699         if ( $record->{ "REC_TYPE" } eq "BED" )                                             # ---- Hits from BED ----
3700         {
3701             Maasha::UCSC::bed_put_entry( $record, $fh, $record->{ "BED_COLS" } );
3702         }
3703         elsif ( $record->{ "REC_TYPE" } eq "PSL" and $record->{ "S_ID" } =~ /^chr/i )       # ---- Hits from BLAT (PSL) ----
3704         {
3705             $new_record->{ "CHR" }     = $record->{ "S_ID" };
3706             $new_record->{ "CHR_BEG" } = $record->{ "S_BEG" };
3707             $new_record->{ "CHR_END" } = $record->{ "S_END" };
3708             $new_record->{ "Q_ID" }    = $record->{ "Q_ID" };
3709             $new_record->{ "SCORE" }   = $record->{ "SCORE" } || 999;
3710             $new_record->{ "STRAND" }  = $record->{ "STRAND" };
3711
3712             Maasha::UCSC::bed_put_entry( $new_record, $fh, 6 );
3713         }
3714         elsif ( $record->{ "REC_TYPE" } eq "PATSCAN" and $record->{ "CHR" } )               # ---- Hits from patscan_seq ----
3715         {
3716             Maasha::UCSC::bed_put_entry( $record, $fh, 6 );
3717         }
3718         elsif ( $record->{ "REC_TYPE" } eq "BLAST" and $record->{ "S_ID" } =~ /^chr/i )     # ---- Hits from BLAST ----
3719         {
3720             $new_record->{ "CHR" }     = $record->{ "S_ID" };
3721             $new_record->{ "CHR_BEG" } = $record->{ "S_BEG" };
3722             $new_record->{ "CHR_END" } = $record->{ "S_END" };
3723             $new_record->{ "Q_ID" }    = $record->{ "Q_ID" };
3724             $new_record->{ "SCORE" }   = $record->{ "SCORE" } || 999; # or use E_VAL somehow
3725             $new_record->{ "STRAND" }  = $record->{ "STRAND" };
3726
3727             Maasha::UCSC::bed_put_entry( $new_record, $fh, 6 );
3728         }
3729         elsif ( $record->{ "REC_TYPE" } eq "VMATCH" and $record->{ "S_ID" } =~ /^chr/i )    # ---- Hits from Vmatch ----
3730         {
3731             $new_record->{ "CHR" }     = $record->{ "S_ID" };
3732             $new_record->{ "CHR_BEG" } = $record->{ "S_BEG" };
3733             $new_record->{ "CHR_END" } = $record->{ "S_END" };
3734             $new_record->{ "Q_ID" }    = $record->{ "Q_ID" };
3735             $new_record->{ "SCORE" }   = $record->{ "SCORE" } || 999; # or use E_VAL somehow
3736             $new_record->{ "STRAND" }  = $record->{ "STRAND" };
3737
3738             Maasha::UCSC::bed_put_entry( $new_record, $fh, 6 );
3739         }
3740         elsif ( $record->{ "CHR" } and defined $record->{ "CHR_BEG" } and $record->{ "CHR_END" } )  # ---- Generic data from tables ----
3741         {
3742             Maasha::UCSC::bed_put_entry( $record, $fh );
3743         }
3744
3745         put_record( $record, $out ) if not $options->{ "no_stream" };
3746     }
3747
3748     close $fh;
3749 }
3750
3751
3752 sub script_write_psl
3753 {
3754     # Martin A. Hansen, August 2007.
3755
3756     # Write PSL output from stream.
3757
3758     my ( $in,        # handle to in stream
3759          $out,       # handle to out stream
3760          $options,   # options hash
3761        ) = @_;
3762
3763     # Returns nothing.
3764
3765     my ( $fh, $record, @output, $first );
3766
3767     $first = 1;
3768
3769     $fh = write_stream( $options->{ "data_out" }, $options->{ "compress" } );
3770
3771     while ( $record = get_record( $in ) )
3772     {
3773         put_record( $record, $out ) if not $options->{ "no_stream" };
3774
3775         if ( $record->{ "REC_TYPE" } and $record->{ "REC_TYPE" } eq "PSL" )
3776         {
3777             Maasha::UCSC::psl_put_header( $fh ) if $first;
3778             Maasha::UCSC::psl_put_entry( $record, $fh );
3779             $first = 0;
3780         }
3781     }
3782
3783     close $fh;
3784 }
3785
3786
3787 sub script_write_2bit
3788 {
3789     # Martin A. Hansen, March 2008.
3790
3791     # Write sequence entries from stream in 2bit format.
3792
3793     my ( $in,        # handle to in stream
3794          $out,       # handle to out stream
3795          $options,   # options hash
3796        ) = @_;
3797
3798     # Returns nothing.
3799
3800     my ( $record, $mask, $tmp_file, $fh_tmp, $fh_in, $fh_out );
3801
3802     $mask = 1 if not $options->{ "no_mask" };
3803
3804     $tmp_file = "$BP_TMP/write_2bit.fna";
3805     $fh_tmp   = Maasha::Common::write_open( $tmp_file );
3806
3807     $fh_out = write_stream( $options->{ "data_out" } );
3808
3809     while ( $record = get_record( $in ) )
3810     {
3811         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } ) {
3812             Maasha::Fasta::put_entry( [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ], $fh_tmp );
3813         }
3814
3815         put_record( $record, $out ) if not $options->{ "no_stream" };
3816     }
3817
3818     close $fh_tmp;
3819
3820     $fh_in = Maasha::Common::read_open( $tmp_file );
3821
3822     Maasha::TwoBit::fasta2twobit( $fh_in, $fh_out, $mask );
3823
3824     close $fh_in;
3825     close $fh_out;
3826
3827     unlink $tmp_file;
3828 }
3829
3830
3831 sub script_write_solid
3832 {
3833     # Martin A. Hansen, April 2008.
3834
3835     # Write di-base encoded Solid sequence from entries in stream.
3836
3837     my ( $in,        # handle to in stream
3838          $out,       # handle to out stream
3839          $options,   # options hash
3840        ) = @_;
3841
3842     # Returns nothing.
3843
3844     my ( $record, $fh, $seq_cs );
3845
3846     $fh = write_stream( $options->{ "data_out" }, $options->{ "compress" } );
3847
3848     while ( $record = get_record( $in ) )
3849     {
3850         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } )
3851         {
3852             $seq_cs = Maasha::Solid::seq2color_space( $record->{ "SEQ" } );
3853
3854             Maasha::Fasta::put_entry( [ $record->{ "SEQ_NAME" }, $seq_cs ], $fh, $options->{ "wrap" } );
3855         }
3856
3857         put_record( $record, $out ) if not $options->{ "no_stream" };
3858     }
3859
3860     close $fh;
3861 }
3862
3863
3864 sub script_plot_seqlogo
3865 {
3866     # Martin A. Hansen, August 2007.
3867
3868     # Calculates and writes a sequence logo for alignments.
3869
3870     my ( $in,        # handle to in stream
3871          $out,       # handle to out stream
3872          $options,   # options hash
3873        ) = @_;
3874
3875     # Returns nothing.
3876
3877     my ( $record, @entries, $logo, $fh );
3878
3879     while ( $record = get_record( $in ) )
3880     {
3881         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } ) {
3882             push @entries, [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ];
3883         }
3884
3885         put_record( $record, $out ) if not $options->{ "no_stream" };
3886     }
3887
3888     $logo = Maasha::Plot::seq_logo( \@entries );
3889
3890     $fh = write_stream( $options->{ "data_out" } );
3891
3892     print $fh $logo;
3893
3894     close $fh;
3895 }
3896
3897
3898 sub script_plot_phastcons_profiles
3899 {
3900     # Martin A. Hansen, January 2008.
3901
3902     # Plots PhastCons profiles.
3903
3904     my ( $in,        # handle to in stream
3905          $out,       # handle to out stream
3906          $options,   # options hash
3907        ) = @_;
3908
3909     # Returns nothing.
3910
3911     my ( $phastcons_file, $phastcons_index, $index, $fh_phastcons, $record, $scores, $AoA, $plot, $fh );
3912
3913     $options->{ "title" } ||= "PhastCons Profiles";
3914
3915     $phastcons_file  = Maasha::Config::genome_phastcons( $options->{ "genome" } );
3916     $phastcons_index = Maasha::Config::genome_phastcons_index( $options->{ "genome" } );
3917
3918     $index           = Maasha::UCSC::phastcons_index_retrieve( $phastcons_index );
3919     $fh_phastcons    = Maasha::Common::read_open( $phastcons_file );
3920
3921     while ( $record = get_record( $in ) )
3922     {
3923         if ( $record->{ "CHR" } and $record->{ "CHR_BEG" } and $record->{ "CHR_END" } )
3924         {
3925             $scores = Maasha::UCSC::phastcons_index_lookup( $index, $fh_phastcons, $record->{ "CHR" }, $record->{ "CHR_BEG" }, $record->{ "CHR_END" }, $options->{ "flank" } );
3926
3927             push @{ $AoA }, [ @{ $scores } ];
3928         }
3929
3930         put_record( $record, $out ) if not $options->{ "no_stream" };
3931     }
3932
3933     Maasha::UCSC::phastcons_normalize( $AoA );
3934
3935     $AoA = [ [ Maasha::UCSC::phastcons_mean( $AoA ) ] ] if $options->{ "mean" };
3936     $AoA = [ [ Maasha::UCSC::phastcons_median( $AoA ) ] ] if $options->{ "median" };
3937
3938     $AoA = Maasha::Matrix::matrix_flip( $AoA );
3939
3940     $plot = Maasha::Plot::lineplot_simple( $AoA, $options, $BP_TMP );
3941
3942     $fh = write_stream( $options->{ "data_out" } );
3943
3944     print $fh "$_\n" foreach @{ $plot };
3945
3946     close $fh;
3947 }
3948
3949
3950 sub script_analyze_bed
3951 {
3952     # Martin A. Hansen, March 2008.
3953
3954     # Analyze BED entries in stream.
3955
3956     my ( $in,        # handle to in stream
3957          $out,       # handle to out stream
3958          $options,   # options hash
3959        ) = @_;
3960
3961     # Returns nothing.
3962
3963     my ( $record );
3964
3965     while ( $record = get_record( $in ) )
3966     {
3967         $record = Maasha::UCSC::bed_analyze( $record ) if $record->{ "REC_TYPE" } eq "BED";
3968
3969         put_record( $record, $out );
3970     }
3971 }
3972
3973
3974 sub script_analyze_vals
3975 {
3976     # Martin A. Hansen, August 2007.
3977
3978     # Analyze values for given keys in stream.
3979
3980     my ( $in,        # handle to in stream
3981          $out,       # handle to out stream
3982          $options,   # options hash
3983        ) = @_;
3984
3985     # Returns nothing.
3986
3987     my ( $record, $key, @keys, %key_hash, $analysis, $len );
3988
3989     map { $key_hash{ $_ } = 1 } @{ $options->{ "keys" } };
3990
3991     while ( $record = get_record( $in ) )
3992     {
3993         foreach $key ( keys %{ $record } )
3994         {
3995             next if $options->{ "keys" } and not exists $key_hash{ $key };
3996
3997             $analysis->{ $key }->{ "COUNT" }++;
3998
3999             if ( Maasha::Calc::is_a_number( $record->{ $key } ) )
4000             {
4001                 $analysis->{ $key }->{ "TYPE" } = "num";
4002                 $analysis->{ $key }->{ "SUM" } += $record->{ $key };
4003                 $analysis->{ $key }->{ "MAX" } = $record->{ $key } if $record->{ $key } > $analysis->{ $key }->{ "MAX" } or not $analysis->{ $key }->{ "MAX" };
4004                 $analysis->{ $key }->{ "MIN" } = $record->{ $key } if $record->{ $key } < $analysis->{ $key }->{ "MIN" } or not $analysis->{ $key }->{ "MIN" };
4005             }
4006             else
4007             {
4008                 $len = length $record->{ $key };
4009
4010                 $analysis->{ $key }->{ "TYPE" } = "alph";
4011                 $analysis->{ $key }->{ "SUM" } += $len;
4012                 $analysis->{ $key }->{ "MAX" } = $len if $len > $analysis->{ $key }->{ "MAX" } or not $analysis->{ $key }->{ "MAX" };
4013                 $analysis->{ $key }->{ "MIN" } = $len if $len < $analysis->{ $key }->{ "MIM" } or not $analysis->{ $key }->{ "MIN" };
4014             }
4015         }
4016
4017         put_record( $record, $out ) if not $options->{ "no_stream" };
4018     }
4019
4020     foreach $key ( keys %{ $analysis } )
4021     {
4022         $analysis->{ $key }->{ "MEAN" } = sprintf "%.2f", $analysis->{ $key }->{ "SUM" } / $analysis->{ $key }->{ "COUNT" };
4023         $analysis->{ $key }->{ "SUM" }  = sprintf "%.2f", $analysis->{ $key }->{ "SUN" };
4024     }
4025
4026     my ( $keys, $types, $counts, $mins, $maxs, $sums, $means );
4027
4028     $keys   = "KEY  ";
4029     $types  = "TYPE ";
4030     $counts = "COUNT";
4031     $mins   = "MIN  ";
4032     $maxs   = "MAX  ";
4033     $sums   = "SUM  ";
4034     $means  = "MEAN ";
4035
4036     if ( $options->{ "keys" } ) {
4037         @keys = @{ $options->{ "keys" } };
4038     } else {
4039         @keys = keys %{ $analysis };
4040     }
4041
4042     foreach $key ( @keys )
4043     {
4044         $keys   .= sprintf "% 15s", $key;
4045         $types  .= sprintf "% 15s", $analysis->{ $key }->{ "TYPE" };
4046         $counts .= sprintf "% 15s", $analysis->{ $key }->{ "COUNT" };
4047         $mins   .= sprintf "% 15s", $analysis->{ $key }->{ "MIN" };
4048         $maxs   .= sprintf "% 15s", $analysis->{ $key }->{ "MAX" };
4049         $sums   .= sprintf "% 15s", $analysis->{ $key }->{ "SUM" };
4050         $means  .= sprintf "% 15s", $analysis->{ $key }->{ "MEAN" };
4051     }
4052
4053     print $out "$keys\n";
4054     print $out "$types\n";
4055     print $out "$counts\n";
4056     print $out "$mins\n";
4057     print $out "$maxs\n";
4058     print $out "$sums\n";
4059     print $out "$means\n";
4060 }
4061
4062
4063 sub script_head_records
4064 {
4065     # Martin A. Hansen, August 2007.
4066
4067     # Display the first sequences in stream.
4068
4069     my ( $in,        # handle to in stream
4070          $out,       # handle to out stream
4071          $options,   # options hash
4072        ) = @_;
4073
4074     # Returns nothing.
4075
4076     my ( $record, $count );
4077
4078     $options->{ "num" } ||= 10;
4079
4080     $count = 0;
4081
4082     while ( $record = get_record( $in ) )
4083     {
4084         $count++;
4085
4086         put_record( $record, $out );
4087
4088         last if $count == $options->{ "num" };
4089     }
4090 }
4091
4092
4093 sub script_remove_keys
4094 {
4095     # Martin A. Hansen, August 2007.
4096
4097     # Remove keys from stream.
4098
4099     my ( $in,        # handle to in stream
4100          $out,       # handle to out stream
4101          $options,   # options hash
4102        ) = @_;
4103
4104     # Returns nothing.
4105
4106     my ( $record, $new_record );
4107
4108     while ( $record = get_record( $in ) )
4109     {
4110         if ( $options->{ "keys" } )
4111         {
4112             map { delete $record->{ $_ } } @{ $options->{ "keys" } };
4113         }
4114         elsif ( $options->{ "save_keys" } )
4115         {
4116             map { $new_record->{ $_ } = $record->{ $_ } if exists $record->{ $_ } } @{ $options->{ "save_keys" } };
4117
4118             $record = $new_record;
4119         }
4120
4121         put_record( $record, $out ) if keys %{ $record };
4122     }
4123 }
4124
4125
4126 sub script_rename_keys
4127 {
4128     # Martin A. Hansen, August 2007.
4129
4130     # Rename keys in stream.
4131
4132     my ( $in,        # handle to in stream
4133          $out,       # handle to out stream
4134          $options,   # options hash
4135        ) = @_;
4136
4137     # Returns nothing.
4138
4139     my ( $record );
4140
4141     while ( $record = get_record( $in ) )
4142     {
4143         if ( exists $record->{ $options->{ "keys" }->[ 0 ] } )
4144         {
4145             $record->{ $options->{ "keys" }->[ 1 ] } = $record->{ $options->{ "keys" }->[ 0 ] };
4146
4147             delete $record->{ $options->{ "keys" }->[ 0 ] };
4148         }
4149
4150         put_record( $record, $out );
4151     }
4152 }
4153
4154
4155 sub script_uniq_vals
4156 {
4157     # Martin A. Hansen, August 2007.
4158
4159     # Find unique values in stream.
4160
4161     my ( $in,        # handle to in stream
4162          $out,       # handle to out stream
4163          $options,   # options hash
4164        ) = @_;
4165
4166     # Returns nothing.
4167
4168     my ( %hash, $record );
4169
4170     while ( $record = get_record( $in ) )
4171     {
4172         if ( $record->{ $options->{ "key" } } )
4173         {
4174             if ( not $hash{ $record->{ $options->{ "key" } } } and not $options->{ "invert" } )
4175             {
4176                 put_record( $record, $out );
4177
4178                 $hash{ $record->{ $options->{ "key" } } } = 1;
4179             }
4180             elsif ( $hash{ $record->{ $options->{ "key" } } } and $options->{ "invert" } )
4181             {
4182                 put_record( $record, $out );
4183             }
4184             else
4185             {
4186                 $hash{ $record->{ $options->{ "key" } } } = 1;
4187             }
4188         }
4189         else
4190         {
4191             put_record( $record, $out );
4192         }
4193     }
4194 }
4195
4196
4197 sub script_merge_vals
4198 {
4199     # Martin A. Hansen, August 2007.
4200
4201     # Rename keys in stream.
4202
4203     my ( $in,        # handle to in stream
4204          $out,       # handle to out stream
4205          $options,   # options hash
4206        ) = @_;
4207
4208     # Returns nothing.
4209
4210     my ( $record, @join, $i );
4211
4212     $options->{ "delimit" } ||= '_';
4213
4214     while ( $record = get_record( $in ) )
4215     {
4216         if ( exists $record->{ $options->{ "keys" }->[ 0 ] } )
4217         {
4218             @join = $record->{ $options->{ "keys" }->[ 0 ] };
4219
4220             for ( $i = 1; $i < @{ $options->{ "keys" } }; $i++ ) {
4221                 push @join, $record->{ $options->{ "keys" }->[ $i ] } if exists $record->{ $options->{ "keys" }->[ $i ] };
4222             }
4223
4224             $record->{ $options->{ "keys" }->[ 0 ] } = join $options->{ "delimit" }, @join;
4225         }
4226
4227         put_record( $record, $out );
4228     }
4229 }
4230
4231
4232 sub script_grab
4233 {
4234     # Martin A. Hansen, August 2007.
4235
4236     # Grab for records in stream.
4237
4238     my ( $in,        # handle to in stream
4239          $out,       # handle to out stream
4240          $options,   # options hash
4241        ) = @_;
4242
4243     # Returns nothing.
4244
4245     my ( $patterns, $pattern, $record, $key, $pos, $op, $val, %lookup_hash );
4246
4247     if ( $options->{ "patterns" } )
4248     {
4249         $patterns = [ split ",", $options->{ "patterns" } ];
4250     }
4251     elsif ( -f $options->{ "patterns_in" } )
4252     {
4253         $patterns = Maasha::Patscan::read_patterns( $options->{ "patterns_in" } );
4254     }
4255     elsif ( -f $options->{ "exact_in" } )
4256     {
4257         $patterns = Maasha::Patscan::read_patterns( $options->{ "exact_in" } );
4258
4259         map { $lookup_hash{ $_ } = 1 } @{ $patterns };
4260
4261         undef $patterns;
4262     }
4263
4264     if ( $options->{ "eval" } )
4265     {
4266         if ( $options->{ "eval" } =~ /^([^><=! ]+)\s*(>=|<=|>|<|=|!=|eq|ne)\s*(.+)$/ )
4267         {
4268             $key = $1;
4269             $op  = $2;
4270             $val = $3;
4271         }
4272     }
4273
4274     while ( $record = get_record( $in ) )
4275     {
4276         $pos = -1;
4277
4278         if ( %lookup_hash )
4279         {
4280             if ( $options->{ "keys" } )
4281             {
4282                 foreach $key ( @{ $options->{ "keys" } } )
4283                 {
4284                     if ( exists $lookup_hash{ $record->{ $key } } )
4285                     {
4286                         $pos = 1;
4287                         goto FOUND;
4288                     }
4289                 }
4290             }
4291             else
4292             {
4293                 foreach $key ( keys %{ $record } )
4294                 {
4295                     if ( not $options->{ "vals_only" } )
4296                     {
4297                         if ( exists $lookup_hash{ $key } )
4298                         {
4299                             $pos = 1;
4300                             goto FOUND;
4301                         }
4302                     }
4303
4304                     if ( not $options->{ "keys_only" } )
4305                     {
4306                         if ( exists $lookup_hash{ $record->{ $key } } )
4307                         {
4308                             $pos = 1;
4309                             goto FOUND;
4310                         }
4311                     }
4312                 }
4313             }
4314         }
4315         elsif ( $patterns )
4316         {
4317             foreach $pattern ( @{ $patterns } )
4318             {
4319                 if ( $options->{ "keys" } )
4320                 {
4321                     foreach $key ( @{ $options->{ "keys" } } )
4322                     {
4323                         $pos = index $record->{ $key }, $pattern;
4324
4325                         goto FOUND if $pos >= 0;
4326                     }
4327                 }
4328                 else
4329                 {
4330                     foreach $key ( keys %{ $record } )
4331                     {
4332                         if ( not $options->{ "vals_only" } )
4333                         {
4334                             $pos = index $key, $pattern;
4335
4336                             goto FOUND if $pos >= 0;
4337                         }
4338
4339                         if ( not $options->{ "keys_only" } )
4340                         {
4341                             $pos = index $record->{ $key }, $pattern;
4342
4343                             goto FOUND if $pos >= 0;
4344                         }
4345                     }
4346                 }
4347             }
4348         }
4349         elsif ( $options->{ "regex" } )
4350         {
4351             if ( $options->{ "keys" } )
4352             {
4353                 foreach $key ( @{ $options->{ "keys" } } )
4354                 {
4355                     if ( $options->{ "case_insensitive" } ) {
4356                         $pos = 1 if $record->{ $key } =~ /$options->{'regex'}/i;
4357                     } else {
4358                         $pos = 1 if $record->{ $key } =~ /$options->{'regex'}/;
4359                     }
4360
4361                     goto FOUND if $pos >= 0;
4362                 }
4363             }
4364             else
4365             {
4366                 foreach $key ( keys %{ $record } )
4367                 {
4368                     if ( not $options->{ "vals_only" } )
4369                     {
4370                         if ( $options->{ "case_insensitive" } ) {
4371                             $pos = 1 if $key =~ /$options->{'regex'}/i;
4372                         } else {
4373                             $pos = 1 if $key =~ /$options->{'regex'}/;
4374                         }
4375
4376                         goto FOUND if $pos >= 0;
4377                     }
4378
4379                     if ( not $options->{ "keys_only" } )
4380                     {
4381                         if ( $options->{ "case_insensitive" } ) {
4382                             $pos = 1 if $record->{ $key } =~ /$options->{'regex'}/i;
4383                         } else {
4384                             $pos = 1 if $record->{ $key } =~ /$options->{'regex'}/;
4385                         }
4386
4387                         goto FOUND if $pos >= 0;
4388                     }
4389                 }
4390             }
4391         }
4392         elsif ( $options->{ "eval" } )
4393         {
4394             if ( defined $record->{ $key } )
4395             {
4396                 if ( $op eq "<" and $record->{ $key } < $val ) {
4397                     $pos = 1 and goto FOUND;
4398                 } elsif ( $op eq ">" and $record->{ $key } > $val ) {
4399                     $pos = 1 and goto FOUND;
4400                 } elsif ( $op eq ">=" and $record->{ $key } >= $val ) {
4401                     $pos = 1 and goto FOUND;
4402                 } elsif ( $op eq "<=" and $record->{ $key } <= $val ) {
4403                     $pos = 1 and goto FOUND;
4404                 } elsif ( $op eq "=" and $record->{ $key } == $val ) {
4405                     $pos = 1 and goto FOUND;
4406                 } elsif ( $op eq "!=" and $record->{ $key } != $val ) {
4407                     $pos = 1 and goto FOUND;
4408                 } elsif ( $op eq "eq" and $record->{ $key } eq $val ) {
4409                     $pos = 1 and goto FOUND;
4410                 } elsif ( $op eq "ne" and $record->{ $key } ne $val ) {
4411                     $pos = 1 and goto FOUND;
4412                 }
4413             }
4414         }
4415
4416         FOUND:
4417
4418         if ( $pos >= 0 and not $options->{ "invert" } ) {
4419             put_record( $record, $out );
4420         } elsif ( $pos < 0 and $options->{ "invert" } ) {
4421             put_record( $record, $out );
4422         }
4423     }
4424 }
4425
4426
4427 sub script_compute
4428 {
4429     # Martin A. Hansen, August 2007.
4430
4431     # Evaluate extression for records in stream.
4432
4433     my ( $in,        # handle to in stream
4434          $out,       # handle to out stream
4435          $options,   # options hash
4436        ) = @_;
4437
4438     # Returns nothing.
4439
4440     my ( $record, $eval_key, $eval_val, $check, @keys );
4441
4442     while ( $record = get_record( $in ) )
4443     {
4444         if ( $options->{ "eval" } )
4445         {
4446             if ( $options->{ "eval" } =~ /^(.+)\s*=\s*(.+)$/ )
4447             {
4448                 $eval_key = $1;
4449                 $eval_val = $2;
4450             }
4451
4452             if ( not $check )
4453             {
4454                 @keys = split /\W+/, $eval_val;
4455                 @keys = grep { ! /^\d+$/ } @keys;
4456
4457                 $check = 1;
4458             }
4459
4460             map { $eval_val =~ s/$_/$record->{ $_ }/g } @keys;
4461
4462             $record->{ $eval_key } = eval "$eval_val" or Maasha::Common::error( "eval failed -> $@" );
4463         }
4464
4465         put_record( $record, $out );
4466     }
4467 }
4468
4469
4470 sub script_flip_tab
4471 {
4472     # Martin A. Hansen, June 2008.
4473
4474     # Flip a table.
4475
4476     my ( $in,        # handle to in stream
4477          $out,       # handle to out stream
4478          $options,   # options hash
4479        ) = @_;
4480
4481     # Returns nothing.
4482
4483     my ( $record, $key, $A, $B, @rows, @matrix, $row, $i );
4484
4485     while ( $record = get_record( $in ) )
4486     {
4487         undef @rows;
4488
4489         foreach $key ( sort { $A = $a; $B = $b; $A =~ s/^V(\d+)$/$1/; $B =~ s/^V(\d+)$/$1/; $A <=> $B } keys %{ $record } )
4490         {
4491             push @rows, $record->{ $key };
4492
4493         }
4494
4495         push @matrix, [ @rows ];
4496     }
4497
4498     undef $record;
4499
4500     @matrix = Maasha::Matrix::matrix_flip( \@matrix );
4501
4502     foreach $row ( @matrix )
4503     {
4504         for ( $i = 0; $i < @{ $row }; $i++ ) {
4505             $record->{ "V$i" } = $row->[ $i ];
4506         }
4507
4508         put_record( $record, $out );
4509     }
4510 }
4511
4512
4513 sub script_add_ident
4514 {
4515     # Martin A. Hansen, May 2008.
4516
4517     # Add a unique identifier to each record in stream.
4518
4519     my ( $in,        # handle to in stream
4520          $out,       # handle to out stream
4521          $options,   # options hash
4522        ) = @_;
4523
4524     # Returns nothing.
4525
4526     my ( $record, $key, $prefix, $i );
4527
4528     $key    = $options->{ "key" }    || "ID";
4529     $prefix = $options->{ "prefix" } || "ID";
4530
4531     $i = 0;
4532
4533     while ( $record = get_record( $in ) )
4534     {
4535         $record->{ $key } = sprintf( "$prefix%08d", $i );
4536
4537         put_record( $record, $out );
4538
4539         $i++;
4540     }
4541 }
4542
4543
4544 sub script_count_records
4545 {
4546     # Martin A. Hansen, August 2007.
4547
4548     # Count records in stream.
4549
4550     my ( $in,        # handle to in stream
4551          $out,       # handle to out stream
4552          $options,   # options hash
4553        ) = @_;
4554
4555     # Returns nothing.
4556
4557     my ( $record, $count, $result, $fh, $line );
4558
4559     $count = 0;
4560
4561     if ( $options->{ "no_stream" } )
4562     {
4563         while ( $line = <$in> )
4564         {
4565             chomp $line;
4566
4567             $count++ if $line eq "---";
4568         }
4569     }
4570     else
4571     {
4572         while ( $record = get_record( $in ) )
4573         {
4574             put_record( $record, $out );
4575
4576             $count++;
4577         }
4578     }
4579
4580     $result = { "count_records" => $count };
4581
4582     $fh = write_stream( $options->{ "data_out" } );
4583
4584     put_record( $result, $fh );
4585
4586     close $fh;
4587 }
4588
4589
4590 sub script_random_records
4591 {
4592     # Martin A. Hansen, August 2007.
4593
4594     # Pick a number or random records from stream.
4595
4596     my ( $in,        # handle to in stream
4597          $out,       # handle to out stream
4598          $options,   # options hash
4599        ) = @_;
4600
4601     # Returns nothing.
4602
4603     my ( $record, $tmp_file, $fh_out, $fh_in, $count, $i, %rand_hash, $rand, $max );
4604
4605     $options->{ "num" } ||= 10;
4606
4607     $tmp_file = "$BP_TMP/random_records.tmp";
4608
4609     $fh_out = Maasha::Common::write_open( $tmp_file );
4610
4611     $count = 0;
4612
4613     while ( $record = get_record( $in ) )
4614     {
4615         put_record( $record, $fh_out );
4616
4617         $count++;
4618     }
4619
4620     close $fh_out;
4621
4622     $max = 0;
4623     $i   = 0;
4624
4625     Maasha::Common::error( qq(Requested random records > records in stream) ) if $options->{ "num" } > $count;
4626
4627     while ( $i < $options->{ "num" } )
4628     {
4629         $rand = int( rand( $count ) );
4630
4631         if ( not exists $rand_hash{ $rand } )
4632         {
4633             $rand_hash{ $rand } = 1;
4634
4635             $max = $rand if $rand > $max;
4636
4637             $i++;
4638         }
4639     }
4640
4641     $fh_in = Maasha::Common::read_open( $tmp_file );
4642
4643     $count = 0;
4644
4645     while ( $record = get_record( $fh_in ) )
4646     {
4647         put_record( $record, $out ) if exists $rand_hash{ $count };
4648
4649         last if $count == $max;
4650
4651         $count++;
4652     }
4653
4654     close $fh_in;
4655
4656     unlink $tmp_file;
4657 }
4658
4659
4660 sub script_sort_records
4661 {
4662     # Martin A. Hansen, August 2007.
4663
4664     # Sort to sort records according to keys.
4665
4666     my ( $in,        # handle to in stream
4667          $out,       # handle to out stream
4668          $options,   # options hash
4669        ) = @_;
4670
4671     # Returns nothing.
4672
4673     my ( @keys, $key, @sort_cmd, $sort_str, $sort_sub, @records, $record, $i );
4674
4675     foreach $key ( @{ $options->{ "keys" } } )
4676     {
4677         if ( $key =~ s/n$// ) {
4678             push @sort_cmd, qq(\$a->{ "$key" } <=> \$b->{ "$key" });
4679         } else {
4680             push @sort_cmd, qq(\$a->{ "$key" } cmp \$b->{ "$key" });
4681         }
4682     }
4683
4684     $sort_str = join " or ", @sort_cmd;
4685     $sort_sub = eval "sub { $sort_str }";   # NB security issue!
4686
4687     while ( $record = get_record( $in ) ) {
4688         push @records, $record;
4689     }
4690
4691     @records = sort $sort_sub @records;
4692
4693     if ( $options->{ "reverse" } )
4694     {
4695         for ( $i = scalar @records - 1; $i >= 0; $i-- ) {
4696             put_record( $records[ $i ], $out );
4697         }
4698     }
4699     else
4700     {
4701         for ( $i = 0; $i < scalar @records; $i++ ) {
4702             put_record( $records[ $i ], $out );
4703         }
4704     }
4705 }
4706
4707
4708 sub script_count_vals
4709 {
4710     # Martin A. Hansen, August 2007.
4711
4712     # Count records in stream.
4713
4714     my ( $in,        # handle to in stream
4715          $out,       # handle to out stream
4716          $options,   # options hash
4717        ) = @_;
4718
4719     # Returns nothing.
4720
4721     my ( $num, $record, %count_hash, @records, $tmp_file, $fh_out, $fh_in, $cache );
4722
4723     $tmp_file = "$BP_TMP/count_cache.tmp";
4724
4725     $fh_out   = Maasha::Common::write_open( $tmp_file );
4726
4727     $num      = 0;
4728
4729     while ( $record = get_record( $in ) )
4730     {
4731         map { $count_hash{ $_ }{ $record->{ $_ } }++ if exists $record->{ $_ } } @{ $options->{ "keys" } };
4732
4733         push @records, $record;
4734
4735         if ( scalar @records > 5_000_000 )   # too many records to hold in memory - use disk cache
4736         {
4737             map { put_record( $_, $fh_out ) } @records;
4738
4739             undef @records;
4740
4741             $cache = 1;
4742         }
4743
4744         print STDERR "verbose: records read $num\n" if ( $options->{ 'verbose' } and ( $num % 1_000_000 ) == 0 );
4745
4746         $num++;
4747     }
4748
4749     close $fh_out;
4750
4751     if ( $cache )
4752     {
4753         $num      = 0;
4754
4755         $fh_in = Maasha::Common::read_open( $tmp_file );
4756
4757         while ( $record = get_record( $fh_in ) )
4758         {
4759             map { $record->{ $_ . "_COUNT" } = $count_hash{ $_ }{ $record->{ $_ } } if exists $record->{ $_ } } @{ $options->{ "keys" } };
4760
4761             put_record( $record, $out );
4762
4763             print STDERR "verbose: cache read $num\n" if ( $options->{ 'verbose' } and ( $num % 1_000_000 ) == 0 );
4764
4765             $num++;
4766         }
4767
4768         close $fh_in;
4769     }
4770
4771     foreach $record ( @records )
4772     {
4773         map { $record->{ $_ . "_COUNT" } = $count_hash{ $_ }{ $record->{ $_ } } if exists $record->{ $_ } } @{ $options->{ "keys" } };
4774
4775         put_record( $record, $out );
4776     }
4777
4778     unlink $tmp_file;
4779 }
4780
4781
4782 sub script_plot_histogram
4783 {
4784     # Martin A. Hansen, September 2007.
4785
4786     # Plot a simple histogram for a given key using GNU plot.
4787
4788     my ( $in,        # handle to in stream
4789          $out,       # handle to out stream
4790          $options,   # options hash
4791        ) = @_;
4792
4793     # Returns nothing.
4794
4795     my ( $record, %data_hash, $max, @data_list, $i, $result, $fh );
4796
4797     $options->{ "title" } ||= "Histogram";
4798     $options->{ "sort" }  ||= "num";
4799
4800     while ( $record = get_record( $in ) )
4801     {
4802         $data_hash{ $record->{ $options->{ "key" } } }++ if $record->{ $options->{ "key" } };
4803
4804         put_record( $record, $out ) if not $options->{ "no_stream" };
4805     }
4806
4807     if ( $options->{ "sort" } eq "num" ) {
4808         map { push @data_list, [ $_, $data_hash{ $_ } ] } sort { $a <=> $b } keys %data_hash;
4809     } else {
4810         map { push @data_list, [ $_, $data_hash{ $_ } ] } sort keys %data_hash;
4811     }
4812
4813     $result = Maasha::Plot::histogram_simple( \@data_list, $options );
4814
4815     $fh = write_stream( $options->{ "data_out" } );
4816
4817     print $fh "$_\n" foreach @{ $result };
4818
4819     close $fh;
4820 }
4821
4822
4823 sub script_plot_lendist
4824 {
4825     # Martin A. Hansen, August 2007.
4826
4827     # Plot length distribution using GNU plot.
4828
4829     my ( $in,        # handle to in stream
4830          $out,       # handle to out stream
4831          $options,   # options hash
4832        ) = @_;
4833
4834     # Returns nothing.
4835
4836     my ( $record, %data_hash, $max, @data_list, $i, $result, $fh );
4837
4838     $options->{ "title" } ||= "Length Distribution";
4839
4840     while ( $record = get_record( $in ) )
4841     {
4842         $data_hash{ $record->{ $options->{ "key" } } }++ if $record->{ $options->{ "key" } };
4843
4844         put_record( $record, $out ) if not $options->{ "no_stream" };
4845     }
4846
4847     $max = Maasha::Calc::list_max( [ keys %data_hash ] );
4848
4849     for ( $i = 0; $i < $max; $i++ ) {
4850         push @data_list, [ $i, $data_hash{ $i } || 0 ];
4851     }
4852
4853     $result = Maasha::Plot::histogram_lendist( \@data_list, $options );
4854
4855     $fh = write_stream( $options->{ "data_out" } );
4856
4857     print $fh "$_\n" foreach @{ $result };
4858
4859     close $fh;
4860 }
4861
4862
4863 sub script_plot_chrdist
4864 {
4865     # Martin A. Hansen, August 2007.
4866
4867     # Plot chromosome distribution using GNU plot.
4868
4869     my ( $in,        # handle to in stream
4870          $out,       # handle to out stream
4871          $options,   # options hash
4872        ) = @_;
4873
4874     # Returns nothing.
4875
4876     my ( $record, %data_hash, @data_list, $elem, $sort_key, $count, $result, $fh );
4877
4878     $options->{ "title" } ||= "Chromosome Distribution";
4879
4880     while ( $record = get_record( $in ) )
4881     {
4882         if ( $record->{ "CHR" } ) {                                                             # generic
4883             $data_hash{ $record->{ "CHR" } }++;
4884         } elsif ( $record->{ "REC_TYPE" } eq "PATSCAN" and $record->{ "S_ID" } =~ /^chr/i ) {   # patscan
4885             $data_hash{ $record->{ "S_ID" } }++;
4886         } elsif ( $record->{ "REC_TYPE" } eq "PSL" and $record->{ "S_ID" } =~ /^chr/i ) {       # BLAT / PSL
4887             $data_hash{ $record->{ "S_ID" } }++;
4888         } elsif ( $record->{ "REC_TYPE" } eq "BLAST" and $record->{ "S_ID" } =~ /^chr/i ) {     # BLAST
4889             $data_hash{ $record->{ "S_ID" } }++;
4890         }
4891
4892         put_record( $record, $out ) if not $options->{ "no_stream" };
4893     }
4894
4895     foreach $elem ( keys %data_hash )
4896     {
4897         $sort_key = $elem;
4898
4899         $sort_key =~ s/chr//i;
4900
4901         $sort_key =~ s/^X(.*)/99$1/;
4902         $sort_key =~ s/^Y(.*)/99$1/;
4903         $sort_key =~ s/^Z(.*)/999$1/;
4904         $sort_key =~ s/^M(.*)/9999$1/;
4905         $sort_key =~ s/^U(.*)/99999$1/;
4906
4907         $count = $sort_key =~ tr/_//;
4908
4909         $sort_key =~ s/_.*/"999999" x $count/ex;
4910
4911         push @data_list, [ $elem, $data_hash{ $elem }, $sort_key ];
4912     }
4913
4914     @data_list = sort { $a->[ 2 ] <=> $b->[ 2 ] } @data_list;
4915
4916     $result = Maasha::Plot::histogram_chrdist( \@data_list, $options );
4917
4918     $fh = write_stream( $options->{ "data_out" } );
4919
4920     print $fh "$_\n" foreach @{ $result };
4921
4922     close $fh;
4923 }
4924
4925
4926 sub script_plot_karyogram
4927 {
4928     # Martin A. Hansen, August 2007.
4929
4930     # Plot hits on karyogram.
4931
4932     my ( $in,        # handle to in stream
4933          $out,       # handle to out stream
4934          $options,   # options hash
4935        ) = @_;
4936
4937     # Returns nothing.
4938
4939     my ( %options, $record, @data, $fh, $result, %data_hash );
4940
4941     $options->{ "genome" }     ||= "human";
4942     $options->{ "feat_color" } ||= "black";
4943
4944     while ( $record = get_record( $in ) )
4945     {
4946         if ( $record->{ "CHR" } and $record->{ "CHR_BEG" } and $record->{ "CHR_END" } )
4947         {
4948             push @{ $data_hash{ $record->{ "CHR" } } }, [ $record->{ "CHR_BEG" }, $record->{ "CHR_END" }, $options->{ "feat_color" } ];
4949         }
4950
4951         put_record( $record, $out ) if not $options->{ "no_stream" };
4952     }
4953
4954     $result = Maasha::Plot::karyogram( \%data_hash, \%options );
4955
4956     $fh = write_stream( $options->{ "data_out" } );
4957
4958     print $fh $result;
4959
4960     close $fh;
4961 }
4962
4963
4964 sub script_plot_matches
4965 {
4966     # Martin A. Hansen, August 2007.
4967
4968     # Plot matches in 2D generating a dotplot.
4969
4970     my ( $in,        # handle to in stream
4971          $out,       # handle to out stream
4972          $options,   # options hash
4973        ) = @_;
4974
4975     # Returns nothing.
4976
4977     my ( $record, @data, $fh, $result, %data_hash );
4978
4979     $options->{ "direction" } ||= "both";
4980
4981     while ( $record = get_record( $in ) )
4982     {
4983         if ( defined $record->{ "Q_BEG" } and defined $record->{ "S_BEG" } and $record->{ "Q_END" } and $record->{ "S_END" } ) {
4984             push @data, $record;
4985         }
4986
4987         put_record( $record, $out ) if not $options->{ "no_stream" };
4988     }
4989
4990     $options->{ "title" }  ||= "plot_matches";
4991     $options->{ "xlabel" } ||= $data[ 0 ]->{ "Q_ID" };
4992     $options->{ "ylabel" } ||= $data[ 0 ]->{ "S_ID" };
4993
4994     $result = Maasha::Plot::dotplot_matches( \@data, $options, $BP_TMP );
4995
4996     $fh = write_stream( $options->{ "data_out" } );
4997
4998     print $fh "$_\n" foreach @{ $result };
4999
5000     close $fh;
5001 }
5002
5003
5004 sub script_length_vals
5005 {
5006     # Martin A. Hansen, August 2007.
5007
5008     # Determine the length of the value for given keys.
5009
5010     my ( $in,        # handle to in stream
5011          $out,       # handle to out stream
5012          $options,   # options hash
5013        ) = @_;
5014
5015     # Returns nothing.
5016
5017     my ( $record, $key );
5018
5019     while ( $record = get_record( $in ) )
5020     {
5021         foreach $key ( @{ $options->{ "keys" } } )
5022         {
5023             if ( $record->{ $key } ) {
5024                 $record->{ $key . "_LEN" } = length $record->{ $key };
5025             }
5026         }
5027
5028         put_record( $record, $out );
5029     }
5030 }
5031
5032
5033 sub script_sum_vals
5034 {
5035     # Martin A. Hansen, August 2007.
5036
5037     # Calculates the sums for values of given keys.
5038
5039     my ( $in,        # handle to in stream
5040          $out,       # handle to out stream
5041          $options,   # options hash
5042        ) = @_;
5043
5044     # Returns nothing.
5045
5046     my ( $record, $key, %sum_hash, $fh );
5047
5048     while ( $record = get_record( $in ) )
5049     {
5050         foreach $key ( @{ $options->{ "keys" } } )
5051         {
5052             if ( $record->{ $key } ) {
5053                 $sum_hash{ $key } += $record->{ $key };
5054             }
5055         }
5056
5057         put_record( $record, $out ) if not $options->{ "no_stream" };
5058     }
5059
5060     $fh = write_stream( $options->{ "data_out" } );
5061
5062     foreach $key ( @{ $options->{ "keys" } } ) {
5063         put_record( { $key . "_SUM" => $sum_hash{ $key } || 0 } , $fh );
5064     }
5065
5066     close $fh;
5067 }
5068
5069
5070 sub script_mean_vals
5071 {
5072     # Martin A. Hansen, August 2007.
5073
5074     # Calculate the mean of values of given keys.
5075
5076     my ( $in,        # handle to in stream
5077          $out,       # handle to out stream
5078          $options,   # options hash
5079        ) = @_;
5080
5081     # Returns nothing.
5082
5083     my ( $record, $key, %sum_hash, %count_hash, $mean, $fh );
5084
5085     while ( $record = get_record( $in ) )
5086     {
5087         foreach $key ( @{ $options->{ "keys" } } )
5088         {
5089             if ( $record->{ $key } )
5090             {
5091                 $sum_hash{ $key } += $record->{ $key };
5092                 $count_hash{ $key }++;
5093             }
5094         }
5095
5096         put_record( $record, $out ) if not $options->{ "no_stream" };
5097     }
5098
5099     $fh = write_stream( $options->{ "data_out" } );
5100
5101     foreach $key ( @{ $options->{ "keys" } } )
5102     {
5103         if ( $count_hash{ $key } ) {
5104             $mean = sprintf( "%.2f", ( $sum_hash{ $key } / $count_hash{ $key } ) );
5105         } else {
5106             $mean = "N/A";
5107         }
5108
5109         put_record( { $key . "_MEAN" => $mean } , $fh );
5110     }
5111
5112     close $fh;
5113 }
5114
5115
5116 sub script_median_vals
5117 {
5118     # Martin A. Hansen, March 2008.
5119
5120     # Calculate the median values of given keys.
5121
5122     my ( $in,        # handle to in stream
5123          $out,       # handle to out stream
5124          $options,   # options hash
5125        ) = @_;
5126
5127     # Returns nothing.
5128
5129     my ( $record, $key, %median_hash, $median, $fh );
5130
5131     while ( $record = get_record( $in ) )
5132     {
5133         foreach $key ( @{ $options->{ "keys" } } ) {
5134             push @{ $median_hash{ $key } }, $record->{ $key } if defined $record->{ $key };
5135         }
5136
5137         put_record( $record, $out ) if not $options->{ "no_stream" };
5138     }
5139
5140     $fh = write_stream( $options->{ "data_out" } );
5141
5142     foreach $key ( @{ $options->{ "keys" } } )
5143     {
5144         if ( $median_hash{ $key } ) {
5145             $median = Maasha::Calc::median( $median_hash{ $key } );
5146         } else {
5147             $median = "N/A";
5148         }
5149
5150         put_record( { $key . "_MEDIAN" => $median } , $fh );
5151     }
5152
5153     close $fh;
5154 }
5155
5156
5157 sub script_max_vals
5158 {
5159     # Martin A. Hansen, February 2008.
5160
5161     # Determine the maximum values of given keys.
5162
5163     my ( $in,        # handle to in stream
5164          $out,       # handle to out stream
5165          $options,   # options hash
5166        ) = @_;
5167
5168     # Returns nothing.
5169
5170     my ( $record, $key, $fh, %max_hash, $max_record );
5171
5172     while ( $record = get_record( $in ) )
5173     {
5174         foreach $key ( @{ $options->{ "keys" } } )
5175         {
5176             if ( $record->{ $key } )
5177             {
5178                 $max_hash{ $key } = $record->{ $key } if $record->{ $key } > $max_hash{ $key };
5179             }
5180         }
5181
5182         put_record( $record, $out ) if not $options->{ "no_stream" };
5183     }
5184
5185     $fh = write_stream( $options->{ "data_out" } );
5186
5187     foreach $key ( @{ $options->{ "keys" } } )
5188     {
5189         $max_record->{ $key . "_MAX" } = $max_hash{ $key };
5190     }
5191
5192     put_record( $max_record, $fh );
5193
5194     close $fh;
5195 }
5196
5197
5198 sub script_min_vals
5199 {
5200     # Martin A. Hansen, February 2008.
5201
5202     # Determine the minimum values of given keys.
5203
5204     my ( $in,        # handle to in stream
5205          $out,       # handle to out stream
5206          $options,   # options hash
5207        ) = @_;
5208
5209     # Returns nothing.
5210
5211     my ( $record, $key, $fh, %min_hash, $min_record );
5212
5213     while ( $record = get_record( $in ) )
5214     {
5215         foreach $key ( @{ $options->{ "keys" } } )
5216         {
5217             if ( defined $record->{ $key } )
5218             {
5219                 if ( exists $min_hash{ $key } ) {
5220                     $min_hash{ $key } = $record->{ $key } if $record->{ $key } < $min_hash{ $key };
5221                 } else {
5222                     $min_hash{ $key } = $record->{ $key };
5223                 }
5224             }
5225         }
5226
5227         put_record( $record, $out ) if not $options->{ "no_stream" };
5228     }
5229
5230     $fh = write_stream( $options->{ "data_out" } );
5231
5232     foreach $key ( @{ $options->{ "keys" } } )
5233     {
5234         $min_record->{ $key . "_MIN" } = $min_hash{ $key };
5235     }
5236
5237     put_record( $min_record, $fh );
5238
5239     close $fh;
5240 }
5241
5242
5243 sub script_upload_to_ucsc
5244 {
5245     # Martin A. Hansen, August 2007.
5246
5247     # Calculate the mean of values of given keys.
5248
5249     my ( $in,        # handle to in stream
5250          $out,       # handle to out stream
5251          $options,   # options hash
5252        ) = @_;
5253
5254     # Returns nothing.
5255
5256     my ( $record, $file, $wib_file, $wig_file, $wib_dir, $fh_in, $fh_out, $i, $first, $format, $args, $type, $columns, $append, %fh_hash,
5257          $chr, $beg, $end, $block, $line, $max, $beg_block, $entry, $q_id, $clones );
5258
5259     $options->{ "short_label" } ||= $options->{ 'table' };
5260     $options->{ "long_label" }  ||= $options->{ 'table' };
5261     $options->{ "group" }       ||= $ENV{ "LOGNAME" };
5262     $options->{ "priority" }    ||= 1;
5263     $options->{ "visibility" }  ||= "pack";
5264     $options->{ "color" }       ||= join( ",", int( rand( 255 ) ), int( rand( 255 ) ), int( rand( 255 ) ) );
5265     $options->{ "chunk_size" }  ||= 10_000_000_000;    # Due to 32-bit UCSC compilation really large tables cannot be loaded in one go.
5266
5267     $file = "$BP_TMP/ucsc_upload.tmp";
5268
5269     $append = 0;
5270
5271     $first = 1;
5272
5273     $i = 0;
5274
5275     if ( $options->{ 'wiggle' } )
5276     {
5277         $options->{ "visibility" } = "full";
5278
5279         while ( $record = get_record( $in ) )
5280         {
5281             put_record( $record, $out ) if not $options->{ "no_stream" };
5282
5283             $record->{ "CHR" }     = $record->{ "S_ID" }  if not defined $record->{ "CHR" };
5284             $record->{ "CHR_BEG" } = $record->{ "S_BEG" } if not defined $record->{ "CHR_BEG" };
5285             $record->{ "CHR_END" } = $record->{ "S_END" } if not defined $record->{ "CHR_END" };
5286
5287             $fh_hash{ $record->{ "CHR" } } = Maasha::Common::write_open( "$BP_TMP/$record->{ 'CHR' }" ) if not exists $fh_hash{ $record->{ "CHR" } };
5288
5289             $fh_out = $fh_hash{ $record->{ "CHR" } };
5290
5291             Maasha::UCSC::bed_put_entry( $record, $fh_out, 5 );
5292         }
5293
5294         map { close $_ } keys %fh_hash;
5295
5296         $fh_out = Maasha::Common::write_open( $file );
5297
5298         foreach $chr ( sort keys %fh_hash )
5299         {
5300             Maasha::Common::run( "bedSort", "$BP_TMP/$chr $BP_TMP/$chr" );
5301
5302             $fh_in = Maasha::Common::read_open( "$BP_TMP/$chr" );
5303
5304             undef $block;
5305
5306             while ( $entry = Maasha::UCSC::bed_get_entry( $fh_in, 5 ) )
5307             {
5308                 $chr  = $entry->{ 'CHR' };
5309                 $beg  = $entry->{ 'CHR_BEG' };
5310                 $end  = $entry->{ 'CHR_END' };
5311                 $q_id = $entry->{ 'Q_ID' };
5312
5313                 if ( $q_id =~ /_(\d+)$/ ) {
5314                     $clones = $1;
5315                 } else {
5316                     $clones = 1;
5317                 }
5318
5319                 if ( $block )
5320                 {
5321                     if ( $beg > $max )
5322                     {
5323                         Maasha::UCSC::fixedstep_put_entry( $chr, $beg_block, $block, $fh_out );
5324                         undef $block;
5325                     }
5326                     else
5327                     {
5328                         for ( $i = $beg - $beg_block; $i < ( $beg - $beg_block ) + ( $end - $beg ); $i++ ) {
5329                             $block->[ $i ] += $clones;
5330                         }
5331
5332                         $max = Maasha::Calc::max( $max, $end );
5333                     }
5334                 }
5335
5336                 if ( not $block )
5337                 {
5338                     $beg_block = $beg;
5339                     $max       = $end;
5340
5341                     for ( $i = 0; $i < ( $end - $beg ); $i++ ) {
5342                         $block->[ $i ] += $clones;
5343                     }
5344                 }
5345             }
5346
5347             close $fh_in;
5348
5349             Maasha::UCSC::fixedstep_put_entry( $chr, $beg_block, $block, $fh_out );
5350
5351             unlink "$BP_TMP/$chr";
5352         }
5353
5354         close $fh_out;
5355
5356         $wig_file = "$options->{ 'table' }.wig";
5357         $wib_file = "$options->{ 'table' }.wib";
5358
5359         $wib_dir  = "$ENV{ 'BP_DATA' }/genomes/$options->{ 'database' }/wib";
5360
5361         Maasha::Common::dir_create_if_not_exists( $wib_dir );
5362
5363         # Maasha::Common::run( "wigEncode", "$file $wig_file $wib_file > /dev/null 2>&1" );
5364
5365         `cd $BP_TMP && wigEncode $file $wig_file $wib_file > /dev/null 2>&1`;
5366         Maasha::Common::run( "mv", "$BP_TMP/$wib_file $wib_dir" );
5367
5368         unlink $file;
5369
5370         $file = $wig_file;
5371
5372         $format = "WIGGLE";
5373     }
5374     else
5375     {
5376         $fh_out = Maasha::Common::write_open( $file );
5377
5378         while ( $record = get_record( $in ) )
5379         {
5380             put_record( $record, $out ) if not $options->{ "no_stream" };
5381
5382             if ( $record->{ "REC_TYPE" } eq "PSL" )
5383             {
5384                 Maasha::UCSC::psl_put_header( $fh_out ) if $first;
5385                 Maasha::UCSC::psl_put_entry( $record, $fh_out );
5386
5387                 $first = 0;
5388
5389                 $format = "PSL" if not $format;
5390             }
5391             elsif ( $record->{ "REC_TYPE" } eq "BED" and $record->{ "SEC_STRUCT" } )
5392             {
5393                 # chrom chromStart  chromEnd    name    score   strand  size    secStr  conf
5394
5395                 print $fh_out join ( "\t",
5396                     $record->{ "CHR" },
5397                     $record->{ "CHR_BEG" },
5398                     $record->{ "CHR_END" } + 1,
5399                     $record->{ "Q_ID" },
5400                     $record->{ "SCORE" },
5401                     $record->{ "STRAND" },
5402                     $record->{ "SIZE" },
5403                     $record->{ "SEC_STRUCT" },
5404                     $record->{ "CONF" },
5405                 ), "\n";
5406
5407                 $format  = "BED_SS" if not $format;
5408             }
5409             elsif ( $record->{ "REC_TYPE" } eq "BED" )
5410             {
5411                 Maasha::UCSC::bed_put_entry( $record, $fh_out, $record->{ "BED_COLS" } );
5412
5413                 $format  = "BED"                   if not $format;
5414                 $columns = $record->{ "BED_COLS" } if not $columns;
5415             }
5416             elsif ( $record->{ "REC_TYPE" } eq "PATSCAN" and $record->{ "CHR" } )
5417             {
5418                 Maasha::UCSC::bed_put_entry( $record, $fh_out, 6 );
5419
5420                 $format  = "BED" if not $format;
5421                 $columns = 6     if not $columns;
5422             }
5423             elsif ( $record->{ "REC_TYPE" } eq "BLAST" and $record->{ "S_ID" } =~ /^chr/ )
5424             {
5425                 $record->{ "CHR" }     = $record->{ "S_ID" };
5426                 $record->{ "CHR_BEG" } = $record->{ "S_BEG" };
5427                 $record->{ "CHR_END" } = $record->{ "S_END" };
5428                 $record->{ "SCORE" }   = $record->{ "BIT_SCORE" } * 1000;
5429
5430                 $format  = "BED" if not $format;
5431                 $columns = 6     if not $columns;
5432
5433                 Maasha::UCSC::bed_put_entry( $record, $fh_out );
5434             }
5435             elsif ( $record->{ "REC_TYPE" } eq "VMATCH" and $record->{ "S_ID" } =~ /^chr/i )
5436             {
5437                 $record->{ "CHR" }     = $record->{ "S_ID" };
5438                 $record->{ "CHR_BEG" } = $record->{ "S_BEG" };
5439                 $record->{ "CHR_END" } = $record->{ "S_END" };
5440                 $record->{ "SCORE" }   = $record->{ "SCORE" } || 999;
5441                 $record->{ "SCORE" }   = int( $record->{ "SCORE" } );
5442
5443                 $format  = "BED" if not $format;
5444                 $columns = 6     if not $columns;
5445
5446                 Maasha::UCSC::bed_put_entry( $record, $fh_out, 6 );
5447             }
5448
5449             if ( $i == $options->{ "chunk_size" } )
5450             {
5451                 close $fh_out;
5452
5453                 if ( $format eq "BED" ) {
5454                     Maasha::UCSC::bed_upload_to_ucsc( $BP_TMP, $file, $options, $append );
5455                 } elsif ( $format eq "PSL" ) {
5456                     Maasha::UCSC::psl_upload_to_ucsc( $file, $options, $append );
5457                 }
5458
5459                 unlink $file;
5460
5461                 $first = 1;
5462
5463                 $append = 1;
5464
5465                 $fh_out = Maasha::Common::write_open( $file );
5466             }
5467
5468             $i++;
5469         }
5470     }
5471
5472     close $fh_out;
5473
5474     if ( exists $options->{ "database" } and $options->{ "table" } )
5475     {
5476         if ( $format eq "BED" )
5477         {
5478             $type = "bed $columns";
5479
5480             Maasha::UCSC::bed_upload_to_ucsc( $BP_TMP, $file, $options, $append );
5481         }
5482         elsif ( $format eq "BED_SS" )
5483         {
5484             $options->{ "sec_struct" } = 1;
5485
5486             $type = "sec_struct";
5487
5488             Maasha::UCSC::bed_upload_to_ucsc( $BP_TMP, $file, $options, $append );
5489         }
5490         elsif ( $format eq "PSL" )
5491         {
5492             $type = "psl";
5493
5494             Maasha::UCSC::psl_upload_to_ucsc( $file, $options, $append );
5495         }
5496         elsif ( $format eq "WIGGLE" )
5497         {
5498             $type = "wig 0";
5499
5500             Maasha::UCSC::wiggle_upload_to_ucsc( $BP_TMP, $wib_dir, $file, $options );
5501         }
5502
5503         unlink $file;
5504
5505         Maasha::UCSC::update_my_tracks( $options, $type );
5506     }
5507 }
5508
5509
5510 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
5511
5512
5513 sub read_stream
5514 {
5515     # Martin A. Hansen, July 2007.
5516
5517     # Opens a stream to STDIN or a file,
5518
5519     my ( $path,   # path - OPTIONAL
5520        ) = @_;
5521
5522     # Returns filehandle.
5523
5524     my ( $fh );
5525
5526     if ( not -t STDIN ) {
5527         $fh = Maasha::Common::read_stdin();
5528     } elsif ( not $path ) {
5529         Maasha::Common::error( qq(no data stream) );
5530     } else {
5531         $fh = Maasha::Common::read_open( $path );
5532     }
5533
5534 #    $fh->autoflush(1) if $fh;  # Disable file buffer for debugging.
5535
5536     return $fh;
5537 }
5538
5539
5540 sub write_stream
5541 {
5542     # Martin A. Hansen, August 2007.
5543
5544     # Opens a stream to STDOUT or a file.
5545
5546     my ( $path,   # path          - OPTIONAL
5547          $gzip,   # compress data - OPTIONAL
5548        ) = @_;
5549
5550     # Returns filehandle.
5551
5552     my ( $fh );
5553
5554     if ( $path ) {
5555         $fh = Maasha::Common::write_open( $path, $gzip );
5556     } else {
5557         $fh = Maasha::Common::write_stdout();
5558     }
5559
5560     return $fh;
5561 }
5562
5563
5564 sub get_record
5565 {
5566     # Martin A. Hansen, July 2007.
5567
5568     # Reads one record at a time and converts that record
5569     # to a Perl data structure (a hash) which is returned.
5570
5571     my ( $fh,
5572        ) = @_;
5573
5574     # Returns data structure.
5575
5576     my ( $block, @lines, $line, $key, $value, %record );
5577
5578     local $/ = "\n---\n";
5579
5580     $block = <$fh>;
5581
5582     chomp $block;
5583
5584     return if not defined $block;
5585
5586     @lines = split "\n", $block;
5587
5588     foreach $line ( @lines )
5589     {
5590         ( $key, $value ) = split ": ", $line;
5591
5592         $record{ $key } = $value;
5593     }
5594
5595     return wantarray ? %record : \%record;
5596 }
5597
5598
5599 sub put_record
5600 {
5601     # Martin A. Hansen, July 2007.
5602
5603     # Given a Perl datastructure (a hash ref) emits this to STDOUT or a filehandle.
5604
5605     my ( $data,   # data structure
5606          $fh,     # file handle - OPTIONAL
5607        ) = @_;
5608
5609     # Returns nothing.
5610
5611     if ( scalar keys %{ $data } )
5612     {
5613         if ( $fh )
5614         {
5615             map { print $fh "$_: $data->{ $_ }\n" } keys %{ $data };
5616             print $fh "---\n";
5617         }
5618         else
5619         {
5620             map { print "$_: $data->{ $_ }\n" } keys %{ $data };
5621             print "---\n";
5622         }
5623     }
5624
5625     undef $data;
5626 }
5627
5628
5629 sub getopt_files
5630 {
5631     # Martin A. Hansen, November 2007.
5632
5633     # Extracts files from an explicit GetOpt::Long argument
5634     # allowing for the use of glob. E.g.
5635     # --data_in=test.fna
5636     # --data_in=test.fna,test2.fna
5637     # --data_in=*.fna
5638     # --data_in=test.fna,/dir/*.fna
5639
5640     my ( $option,   # option from GetOpt::Long
5641        ) = @_;
5642
5643     # Returns a list.
5644
5645     my ( $elem, @files );
5646
5647     foreach $elem ( split ",", $option )
5648     {
5649         if ( -f $elem ) {
5650             push @files, $elem;
5651         } elsif ( $elem =~ /\*/ ) {
5652             push @files, glob( $elem );
5653         }
5654     }
5655
5656     return wantarray ? @files : \@files;
5657 }
5658
5659
5660 sub sig_handler
5661 {
5662     # Martin A. Hansen, April 2008.
5663
5664     # Removes temporary directory and exits gracefully.
5665     # This subroutine is meant to be run always as the last
5666     # thing even if a script is dies or is interrupted
5667     # or killed.
5668
5669     my ( $sig,   # signal from the %SIG
5670        ) = @_;
5671
5672     # print STDERR "signal->$sig<-\n";
5673
5674     chomp $sig;
5675
5676     sleep 1;
5677
5678     if ( -d $BP_TMP )
5679     {
5680         if ( $sig =~ /MAASHA_ERROR/ ) {
5681             print STDERR "\nProgram '$script' had an error"                     . "  -  Please wait for temporary data to be removed\n";
5682         } elsif ( $sig eq "INT" ) {
5683             print STDERR "\nProgram '$script' interrupted (ctrl-c was pressed)" . "  -  Please wait for temporary data to be removed\n";
5684         } elsif ( $sig eq "TERM" ) {
5685             print STDERR "\nProgram '$script' terminated (someone used kill?)"  . "  -  Please wait for temporary data to be removed\n";
5686         } else {
5687             print STDERR "\nProgram '$script' died->$sig"                       . "  -  Please wait for temporary data to be removed\n";
5688         }
5689
5690         # This is a really bad solution, potentially, anyone can include this module and set
5691         # the BP_TMP to point at any dir and thus take out the machine !!!
5692
5693         Maasha::Common::dir_remove( $BP_TMP );
5694     }
5695
5696     exit( 0 );
5697 }
5698
5699
5700 END
5701 {
5702     # This is a really bad solution, potentially, anyone can include this module and set
5703     # the BP_TMP to point at any dir and thus take out the machine !!!
5704
5705     Maasha::Common::dir_remove( $BP_TMP );
5706 }
5707
5708
5709 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
5710
5711 1;
5712
5713 __END__
5714
5715
5716 sub script_read_soft
5717 {
5718     # Martin A. Hansen, December 2007.
5719
5720     # Read soft format.
5721     # http://www.ncbi.nlm.nih.gov/geo/info/soft2.html
5722
5723     my ( $in,        # handle to in stream
5724          $out,       # handle to out stream
5725          $options,   # options hash
5726        ) = @_;
5727
5728     # Returns nothing.
5729
5730     my ( $data_in, $file, $num, $records, $record );
5731
5732     while ( $record = get_record( $in ) ) {
5733         put_record( $record, $out );
5734     }
5735
5736     $num = 1;
5737
5738     foreach $file ( @{ $options->{ "files" } } )
5739     {
5740         $records = Maasha::NCBI::soft_parse( $file );
5741
5742         foreach $record ( @{ $records } )
5743         {
5744             put_record( $record, $out );
5745
5746             goto NUM if $options->{ "num" } and $num == $options->{ "num" };
5747
5748             $num++;
5749         }
5750     }
5751
5752     NUM:
5753
5754     close $data_in if $data_in;
5755 }