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