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