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