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