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