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     Maasha::Common::run( "blastall", "-p $options->{ 'program' } -e $options->{ 'e_val' } -a $options->{ 'cpus' } -m 8 -i $tmp_in -d $options->{ 'database' } -F $options->{ 'filter' } -o $tmp_out > /dev/null 2>&1", 1 );
3477
3478     unlink $tmp_in;
3479
3480     $fh_out = Maasha::Common::read_open( $tmp_out );
3481
3482     undef $record;
3483
3484     while ( $line = <$fh_out> )
3485     {
3486         chomp $line;
3487
3488         next if $line =~ /^#/;
3489
3490         @fields = split /\s+/, $line;
3491
3492         $record->{ "REC_TYPE" }   = "BLAST";
3493         $record->{ "Q_ID" }       = $fields[ 0 ];
3494         $record->{ "S_ID" }       = $fields[ 1 ];
3495         $record->{ "IDENT" }      = $fields[ 2 ];
3496         $record->{ "ALIGN_LEN" }  = $fields[ 3 ];
3497         $record->{ "MISMATCHES" } = $fields[ 4 ];
3498         $record->{ "GAPS" }       = $fields[ 5 ];
3499         $record->{ "Q_BEG" }      = $fields[ 6 ] - 1; # BLAST is 1-based
3500         $record->{ "Q_END" }      = $fields[ 7 ] - 1; # BLAST is 1-based
3501         $record->{ "S_BEG" }      = $fields[ 8 ] - 1; # BLAST is 1-based
3502         $record->{ "S_END" }      = $fields[ 9 ] - 1; # BLAST is 1-based
3503         $record->{ "E_VAL" }      = $fields[ 10 ];
3504         $record->{ "BIT_SCORE" }  = $fields[ 11 ];
3505
3506         if ( $record->{ "S_BEG" } > $record->{ "S_END" } )
3507         {
3508             $record->{ "STRAND" } = '-';
3509
3510             ( $record->{ "S_BEG" }, $record->{ "S_END" } ) = ( $record->{ "S_END" }, $record->{ "S_BEG" } );
3511         }
3512         else
3513         {
3514             $record->{ "STRAND" } = '+';
3515         }
3516
3517         put_record( $record, $out );
3518     }
3519
3520     close $fh_out;
3521
3522     unlink $tmp_out;
3523 }
3524
3525
3526 sub script_blat_seq
3527 {
3528     # Martin A. Hansen, August 2007.
3529
3530     # BLATs sequences in stream against a given genome.
3531
3532     my ( $in,        # handle to in stream
3533          $out,       # handle to out stream
3534          $options,   # options hash
3535        ) = @_;
3536
3537     # Returns nothing.
3538
3539     my ( $blat_args, $genome_file, $query_file, $fh_in, $fh_out, $type, $record, $result_file, $entries, $entry );
3540
3541     $genome_file = "$ENV{ 'BP_DATA' }/genomes/$options->{ 'genome' }/fasta/$options->{ 'genome' }.fna";
3542
3543     $options->{ 'tile_size' }    ||= 11;
3544     $options->{ 'one_off' }      ||= 0;
3545     $options->{ 'min_identity' } ||= 90;
3546     $options->{ 'min_score' }    ||= 0;
3547     $options->{ 'step_size' }    ||= $options->{ 'tile_size' };
3548
3549     $blat_args .= " -tileSize=$options->{ 'tile_size' }";
3550     $blat_args .= " -oneOff=$options->{ 'one_off' }";
3551     $blat_args .= " -minIdentity=$options->{ 'min_identity' }";
3552     $blat_args .= " -minScore=$options->{ 'min_score' }";
3553     $blat_args .= " -stepSize=$options->{ 'step_size' }";
3554 #    $blat_args .= " -ooc=" . Maasha::Config::genome_blat_ooc( $options->{ "genome" }, 11 ) if $options->{ 'ooc' };
3555
3556     $query_file = "$BP_TMP/blat.seq";
3557
3558     $fh_out = Maasha::Common::write_open( $query_file );
3559
3560     while ( $record = get_record( $in ) )
3561     {
3562         if ( $entry = record2fasta( $record ) )
3563         {
3564             $type = Maasha::Seq::seq_guess_type( $entry->[ SEQ ] ) if not $type;
3565             Maasha::Fasta::put_entry( $entry, $fh_out, 80 );
3566         }
3567
3568         put_record( $record, $out );
3569     }
3570
3571     close $fh_out;
3572
3573     $blat_args .= " -t=dnax" if $type eq "protein";
3574     $blat_args .= " -q=$type";
3575
3576     $result_file = "$BP_TMP/blat.psl";
3577
3578     Maasha::Common::run( "blat", "$genome_file $query_file $blat_args $result_file > /dev/null 2>&1" );
3579
3580     unlink $query_file;
3581
3582     $entries = Maasha::UCSC::psl_get_entries( $result_file );
3583
3584     map { put_record( $_, $out ) } @{ $entries };
3585
3586     unlink $result_file;
3587 }
3588
3589
3590 sub script_soap_seq
3591 {
3592     # Martin A. Hansen, July 2008.
3593
3594     # soap sequences in stream against a given file or genome.
3595
3596     my ( $in,        # handle to in stream
3597          $out,       # handle to out stream
3598          $options,   # options hash
3599        ) = @_;
3600
3601     # Returns nothing.
3602
3603     my ( $genome, $tmp_in, $tmp_out, $fh_in, $fh_out, $record, $line, @fields, $entry, $count, $args );
3604
3605     $options->{ "seed_size" }  ||= 10;
3606     $options->{ "mismatches" } ||= 2;
3607     $options->{ "gap_size" }   ||= 0;
3608     $options->{ "cpus" }       ||= 1;
3609
3610     if ( $options->{ "genome" } ) {
3611         $options->{ "in_file" } = "$ENV{ 'BP_DATA' }/genomes/$options->{ 'genome' }/fasta/$options->{ 'genome' }.fna";
3612     }
3613
3614     $tmp_in  = "$BP_TMP/soap_query.seq";
3615     $tmp_out = "$BP_TMP/soap.result";
3616
3617     $fh_out = Maasha::Common::write_open( $tmp_in );
3618
3619     $count = 0;
3620
3621     while ( $record = get_record( $in ) )
3622     {
3623         if ( $entry = record2fasta( $record ) )
3624         {
3625             Maasha::Fasta::put_entry( $entry, $fh_out );
3626
3627             $count++;
3628         }
3629
3630         put_record( $record, $out );
3631     }
3632
3633     close $fh_out;
3634
3635     if ( $count > 0 )
3636     {
3637         $args = join( " ",
3638             "-s $options->{ 'seed_size' }",
3639             "-r 2",
3640             "-a $tmp_in",
3641             "-v $options->{ 'mismatches' }",
3642             "-g $options->{ 'gap_size' }",
3643             "-p $options->{ 'cpus' }",
3644             "-d $options->{ 'in_file' }",
3645             "-o $tmp_out",
3646         );
3647
3648         $args .= " > /dev/null 2>&1" if not $options->{ 'verbose' };
3649
3650         Maasha::Common::run( "soap", $args, 1 );
3651
3652         unlink $tmp_in;
3653
3654         $fh_out = Maasha::Common::read_open( $tmp_out );
3655
3656         undef $record;
3657
3658         while ( $line = <$fh_out> )
3659         {
3660             chomp $line;
3661
3662             @fields = split /\t/, $line;
3663
3664             $record->{ "REC_TYPE" }   = "SOAP";
3665             $record->{ "Q_ID" }       = $fields[ 0 ];
3666             $record->{ "SCORE" }      = $fields[ 3 ];
3667             $record->{ "STRAND" }     = $fields[ 6 ];
3668             $record->{ "S_ID" }       = $fields[ 7 ];
3669             $record->{ "S_BEG" }      = $fields[ 8 ] - 1; # soap is 1-based
3670             $record->{ "S_END" }      = $fields[ 8 ] + $fields[ 5 ] - 2;
3671
3672             put_record( $record, $out );
3673         }
3674
3675         close $fh_out;
3676     }
3677
3678     unlink $tmp_out;
3679 }
3680
3681
3682 sub script_match_seq
3683 {
3684     # Martin A. Hansen, August 2007.
3685
3686     # BLATs sequences in stream against a given genome.
3687
3688     my ( $in,        # handle to in stream
3689          $out,       # handle to out stream
3690          $options,   # options hash
3691        ) = @_;
3692
3693     # Returns nothing.
3694
3695     my ( $record, @entries, $results );
3696
3697     $options->{ "word_size" } ||= 20;
3698     $options->{ "direction" } ||= "both";
3699
3700     while ( $record = get_record( $in ) )
3701     {
3702         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } ) {
3703             push @entries, [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ];
3704         }
3705
3706         put_record( $record, $out );
3707     }
3708
3709     if ( @entries == 1 )
3710     {
3711         $results = Maasha::Match::match_mummer( [ $entries[ 0 ] ], [ $entries[ 0 ] ], $options, $BP_TMP );
3712
3713         map { put_record( $_, $out ) } @{ $results };
3714     }
3715     elsif ( @entries == 2 )
3716     {
3717         $results = Maasha::Match::match_mummer( [ $entries[ 0 ] ], [ $entries[ 1 ] ], $options, $BP_TMP );
3718
3719         map { put_record( $_, $out ) } @{ $results };
3720     }
3721 }
3722
3723
3724 sub script_create_vmatch_index
3725 {
3726     # Martin A. Hansen, January 2008.
3727
3728     # Create a vmatch index from sequences in the stream.
3729
3730     my ( $in,        # handle to in stream
3731          $out,       # handle to out stream
3732          $options,   # options hash
3733        ) = @_;
3734
3735     # Returns nothing.
3736
3737     my ( $record, $file_tmp, $fh_tmp, $type, $entry );
3738
3739     if ( $options->{ "index_name" } )
3740     {
3741         $file_tmp = $options->{ 'index_name' };
3742         $fh_tmp   = Maasha::Common::write_open( $file_tmp );
3743     }
3744
3745     while ( $record = get_record( $in ) )
3746     {
3747         if ( $options->{ "index_name" } and $entry = record2fasta( $record ) )
3748         {
3749             Maasha::Fasta::put_entry( $entry, $fh_tmp );
3750
3751             $type = Maasha::Seq::seq_guess_type( $entry->[ SEQ ] ) if not defined $type;
3752         }
3753
3754         put_record( $record, $out ) if not $options->{ "no_stream" };
3755     }
3756
3757     if ( $options->{ "index_name" } )
3758     {
3759         close $fh_tmp;
3760
3761         if ( $type eq "protein" ) {
3762             Maasha::Common::run( "mkvtree", "-db $file_tmp -protein -pl $options->{ 'prefix_length' } -allout -indexname $file_tmp > /dev/null 2>&1" );
3763         } else {
3764             Maasha::Common::run( "mkvtree", "-db $file_tmp -dna -pl $options->{ 'prefix_length' } -allout -indexname $file_tmp > /dev/null 2>&1" );
3765         }
3766
3767         unlink $file_tmp;
3768     }
3769 }
3770
3771
3772 sub script_vmatch_seq
3773 {
3774     # Martin A. Hansen, August 2007.
3775
3776     # Vmatches sequences in stream against a given genome.
3777
3778     my ( $in,        # handle to in stream
3779          $out,       # handle to out stream
3780          $options,   # options hash
3781        ) = @_;
3782
3783     # Returns nothing.
3784
3785     my ( @index_files, @records, $result_file, $fh_in, $record, %hash );
3786
3787     $options->{ 'count' } = 1 if $options->{ 'max_hits' };
3788
3789     if ( $options->{ "index_name" } )
3790     {
3791         @index_files = $options->{ "index_name" };
3792     }
3793     else
3794     {
3795         @index_files = Maasha::Common::ls_files( "$ENV{ 'BP_DATA' }/genomes/$options->{ 'genome' }/vmatch" );
3796
3797         map { $_ =~ /^(.+)\.[a-z1]{3,4}$/; $hash{ $1 } = 1 } @index_files;
3798
3799         @index_files = sort keys %hash;
3800     }
3801
3802     while ( $record = get_record( $in ) )
3803     {
3804         push @records, $record;
3805
3806         put_record( $record, $out );
3807     }
3808
3809     $result_file = Maasha::Match::match_vmatch( $BP_TMP, \@records, \@index_files, $options );
3810
3811     undef @records;
3812
3813     $fh_in = Maasha::Common::read_open( $result_file );
3814
3815     while ( $record = Maasha::Match::vmatch_get_entry( $fh_in ) ) {
3816         put_record( $record, $out );
3817     }
3818
3819     close $fh_in;
3820
3821     unlink $result_file;
3822 }
3823
3824
3825 sub script_write_fasta
3826 {
3827     # Martin A. Hansen, August 2007.
3828
3829     # Write FASTA entries from sequences in stream.
3830
3831     my ( $in,        # handle to in stream
3832          $out,       # handle to out stream
3833          $options,   # options hash
3834        ) = @_;
3835
3836     # Returns nothing.
3837
3838     my ( $record, $fh, $entry );
3839
3840     $fh = write_stream( $options->{ "data_out" }, $options->{ "compress" } );
3841
3842     while ( $record = get_record( $in ) )
3843     {
3844         if ( $entry = record2fasta( $record ) ) {
3845             Maasha::Fasta::put_entry( $entry, $fh, $options->{ "wrap" } );
3846         }
3847
3848         put_record( $record, $out ) if not $options->{ "no_stream" };
3849     }
3850
3851     close $fh;
3852 }
3853
3854
3855 sub script_write_align
3856 {
3857     # Martin A. Hansen, August 2007.
3858
3859     # Write pretty alignments aligned sequences in stream.
3860
3861     my ( $in,        # handle to in stream
3862          $out,       # handle to out stream
3863          $options,   # options hash
3864        ) = @_;
3865
3866     # Returns nothing.
3867
3868     my ( $fh, $record, @entries );
3869
3870     $fh = write_stream( $options->{ "data_out" } ) ;
3871
3872     while ( $record = get_record( $in ) )
3873     {
3874         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } ) {
3875             push @entries, [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ];
3876         }
3877
3878         put_record( $record, $out ) if not $options->{ "no_stream" };
3879     }
3880
3881     if ( scalar( @entries ) == 2 ) {
3882         Maasha::Align::align_print_pairwise( $entries[ 0 ], $entries[ 1 ], $fh, $options->{ "wrap" } );
3883     } elsif ( scalar ( @entries ) > 2 ) {
3884         Maasha::Align::align_print_multi( \@entries, $fh, $options->{ "wrap" }, $options->{ "no_ruler" }, $options->{ "no_consensus" } );
3885     }
3886
3887     close $fh if $fh;
3888 }
3889
3890
3891 sub script_write_blast
3892 {
3893     # Martin A. Hansen, November 2007.
3894
3895     # Write data in blast table format (-m8 and 9).
3896
3897     my ( $in,        # handle to in stream
3898          $out,       # handle to out stream
3899          $options,   # options hash
3900        ) = @_;
3901
3902     # Returns nothing.
3903
3904     my ( $fh, $record, $first );
3905
3906     $fh = write_stream( $options->{ "data_out" }, $options->{ "compress" } ) ;
3907
3908     $first = 1;
3909
3910     while ( $record = get_record( $in ) )
3911     {
3912         if ( $record->{ "REC_TYPE" } eq "BLAST" )
3913         {
3914             if ( $options->{ "comment" } and $first )
3915             {
3916                 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";
3917
3918                 $first = 0;
3919             }
3920
3921             if ( $record->{ "STRAND" } eq "-" ) {
3922                 ( $record->{ "S_BEG" }, $record->{ "S_END" } ) = ( $record->{ "S_END" }, $record->{ "S_BEG" } );
3923             }
3924
3925             print $fh join( "\t",
3926                 $record->{ "Q_ID" },
3927                 $record->{ "S_ID" },
3928                 $record->{ "IDENT" },
3929                 $record->{ "ALIGN_LEN" },
3930                 $record->{ "MISMATCHES" },
3931                 $record->{ "GAPS" },
3932                 $record->{ "Q_BEG" } + 1,
3933                 $record->{ "Q_END" } + 1,
3934                 $record->{ "S_BEG" } + 1,
3935                 $record->{ "S_END" } + 1,
3936                 $record->{ "E_VAL" },
3937                 $record->{ "BIT_SCORE" }
3938             ), "\n";
3939         }
3940
3941         put_record( $record, $out ) if not $options->{ "no_stream" };
3942     }
3943
3944     close $fh;
3945 }
3946
3947
3948 sub script_write_tab
3949 {
3950     # Martin A. Hansen, August 2007.
3951
3952     # Write data as table.
3953
3954     my ( $in,        # handle to in stream
3955          $out,       # handle to out stream
3956          $options,   # options hash
3957        ) = @_;
3958
3959     # Returns nothing.
3960
3961     my ( $fh, $record, $key, @keys, @vals, $ok, %no_keys, $A, $B );
3962
3963     $options->{ "delimit" } ||= "\t";
3964
3965     map { $no_keys{ $_ } = 1 } @{ $options->{ "no_keys" } };
3966
3967     $fh = write_stream( $options->{ "data_out" }, $options->{ "compress" } );
3968
3969     while ( $record = get_record( $in ) )
3970     {
3971         undef @vals;
3972         $ok = 1;
3973
3974         if ( $options->{ "keys" } )
3975         {
3976             map { $ok = 0 if not exists $record->{ $_ } } @{ $options->{ "keys" } };
3977
3978             if ( $ok )
3979             {
3980                 foreach $key ( @{ $options->{ "keys" }  } )
3981                 {
3982                     if ( exists $record->{ $key } )
3983                     {
3984                         push @keys, $key if $options->{ "comment" };
3985                         push @vals, $record->{ $key };
3986                     }
3987                 }
3988              }
3989         }
3990         else
3991         {
3992             foreach $key ( sort { $A = $a; $B = $b; $A =~ s/^V(\d+)$/$1/; $B =~ s/^V(\d+)$/$1/; $A <=> $B } keys %{ $record } )
3993             {
3994                 next if exists $no_keys{ $key };
3995
3996                 push @keys, $key if $options->{ "comment" };
3997                 push @vals, $record->{ $key };
3998             }
3999         }
4000
4001         if ( @keys and $options->{ "comment" } )
4002         {
4003             print $fh "#", join( $options->{ "delimit" }, @keys ), "\n";
4004
4005             delete $options->{ "comment" };
4006         }
4007
4008         print $fh join( $options->{ "delimit" }, @vals ), "\n" if @vals;
4009
4010         put_record( $record, $out ) if not $options->{ "no_stream" };
4011     }
4012
4013     close $fh;
4014 }
4015
4016
4017 sub script_write_bed
4018 {
4019     # Martin A. Hansen, August 2007.
4020
4021     # Write BED format for the UCSC genome browser using records in stream.
4022
4023     my ( $in,        # handle to in stream
4024          $out,       # handle to out stream
4025          $options,   # options hash
4026        ) = @_;
4027
4028     # Returns nothing.
4029
4030     my ( $fh, $record, $new_record );
4031
4032     $fh = write_stream( $options->{ "data_out" }, $options->{ "compress" } );
4033
4034     while ( $record = get_record( $in ) )
4035     {
4036         if ( $record->{ "REC_TYPE" } eq "BED" )                                             # ---- Hits from BED ----
4037         {
4038             Maasha::UCSC::bed_put_entry( $record, $fh, $record->{ "BED_COLS" } );
4039         }
4040         elsif ( $record->{ "REC_TYPE" } eq "PSL" and $record->{ "S_ID" } =~ /^chr/i )       # ---- Hits from BLAT (PSL) ----
4041         {
4042             $new_record->{ "CHR" }     = $record->{ "S_ID" };
4043             $new_record->{ "CHR_BEG" } = $record->{ "S_BEG" };
4044             $new_record->{ "CHR_END" } = $record->{ "S_END" };
4045             $new_record->{ "Q_ID" }    = $record->{ "Q_ID" };
4046             $new_record->{ "SCORE" }   = $record->{ "SCORE" } || 999;
4047             $new_record->{ "STRAND" }  = $record->{ "STRAND" };
4048
4049             Maasha::UCSC::bed_put_entry( $new_record, $fh, 6 );
4050         }
4051         elsif ( $record->{ "REC_TYPE" } eq "PATSCAN" and $record->{ "CHR" } )               # ---- Hits from patscan_seq ----
4052         {
4053             Maasha::UCSC::bed_put_entry( $record, $fh, 6 );
4054         }
4055         elsif ( $record->{ "REC_TYPE" } eq "BLAST" and $record->{ "S_ID" } =~ /^chr/i )     # ---- Hits from BLAST ----
4056         {
4057             $new_record->{ "CHR" }     = $record->{ "S_ID" };
4058             $new_record->{ "CHR_BEG" } = $record->{ "S_BEG" };
4059             $new_record->{ "CHR_END" } = $record->{ "S_END" };
4060             $new_record->{ "Q_ID" }    = $record->{ "Q_ID" };
4061             $new_record->{ "SCORE" }   = $record->{ "SCORE" } || 999; # or use E_VAL somehow
4062             $new_record->{ "STRAND" }  = $record->{ "STRAND" };
4063
4064             Maasha::UCSC::bed_put_entry( $new_record, $fh, 6 );
4065         }
4066         elsif ( $record->{ "REC_TYPE" } eq "VMATCH" and $record->{ "S_ID" } =~ /^chr/i )    # ---- Hits from Vmatch ----
4067         {
4068             $new_record->{ "CHR" }     = $record->{ "S_ID" };
4069             $new_record->{ "CHR_BEG" } = $record->{ "S_BEG" };
4070             $new_record->{ "CHR_END" } = $record->{ "S_END" };
4071             $new_record->{ "Q_ID" }    = $record->{ "Q_ID" };
4072             $new_record->{ "SCORE" }   = $record->{ "SCORE" } || 999; # or use E_VAL somehow
4073             $new_record->{ "STRAND" }  = $record->{ "STRAND" };
4074
4075             Maasha::UCSC::bed_put_entry( $new_record, $fh, 6 );
4076         }
4077         elsif ( $record->{ "REC_TYPE" } eq "SOAP" and $record->{ "S_ID" } =~ /^chr/i )    # ---- Hits from Vmatch ----
4078         {
4079             $new_record->{ "CHR" }     = $record->{ "S_ID" };
4080             $new_record->{ "CHR_BEG" } = $record->{ "S_BEG" };
4081             $new_record->{ "CHR_END" } = $record->{ "S_END" };
4082             $new_record->{ "Q_ID" }    = $record->{ "Q_ID" };
4083             $new_record->{ "SCORE" }   = $record->{ "SCORE" } || 999;
4084             $new_record->{ "STRAND" }  = $record->{ "STRAND" };
4085
4086             Maasha::UCSC::bed_put_entry( $new_record, $fh, 6 );
4087         }
4088         elsif ( $record->{ "CHR" } and defined $record->{ "CHR_BEG" } and $record->{ "CHR_END" } )  # ---- Generic data from tables ----
4089         {
4090             Maasha::UCSC::bed_put_entry( $record, $fh );
4091         }
4092
4093         put_record( $record, $out ) if not $options->{ "no_stream" };
4094     }
4095
4096     close $fh;
4097 }
4098
4099
4100 sub script_write_psl
4101 {
4102     # Martin A. Hansen, August 2007.
4103
4104     # Write PSL output from stream.
4105
4106     my ( $in,        # handle to in stream
4107          $out,       # handle to out stream
4108          $options,   # options hash
4109        ) = @_;
4110
4111     # Returns nothing.
4112
4113     my ( $fh, $record, @output, $first );
4114
4115     $first = 1;
4116
4117     $fh = write_stream( $options->{ "data_out" }, $options->{ "compress" } );
4118
4119     while ( $record = get_record( $in ) )
4120     {
4121         put_record( $record, $out ) if not $options->{ "no_stream" };
4122
4123         if ( $record->{ "REC_TYPE" } and $record->{ "REC_TYPE" } eq "PSL" )
4124         {
4125             Maasha::UCSC::psl_put_header( $fh ) if $first;
4126             Maasha::UCSC::psl_put_entry( $record, $fh );
4127             $first = 0;
4128         }
4129     }
4130
4131     close $fh;
4132 }
4133
4134
4135 sub script_write_fixedstep
4136 {
4137     # Martin A. Hansen, Juli 2008.
4138
4139     # Write fixedStep entries from recrods in the stream.
4140
4141     my ( $in,        # handle to in stream
4142          $out,       # handle to out stream
4143          $options,   # options hash
4144        ) = @_;
4145
4146     # Returns nothing.
4147
4148     my ( $fh, $record, $vals );
4149
4150     $fh = write_stream( $options->{ "data_out" }, $options->{ "compress" } );
4151
4152     while ( $record = get_record( $in ) )
4153     {
4154         if ( $record->{ "CHR" } and $record->{ "CHR_BEG" } and $record->{ "STEP" } and $record->{ "VALS" } )
4155         {
4156             print $fh "fixedStep chrom=$record->{ 'CHR' } start=$record->{ 'CHR_BEG' } step=$record->{ 'STEP' }\n";
4157
4158             $vals = $record->{ 'VALS' };
4159
4160             $vals =~ tr/,/\n/;
4161
4162             print $fh "$vals\n";
4163         }
4164
4165         put_record( $record, $out ) if not $options->{ "no_stream" };
4166     }
4167
4168     close $fh;
4169 }
4170
4171
4172 sub script_write_2bit
4173 {
4174     # Martin A. Hansen, March 2008.
4175
4176     # Write sequence entries from stream in 2bit format.
4177
4178     my ( $in,        # handle to in stream
4179          $out,       # handle to out stream
4180          $options,   # options hash
4181        ) = @_;
4182
4183     # Returns nothing.
4184
4185     my ( $record, $mask, $tmp_file, $fh_tmp, $fh_in, $fh_out, $entry );
4186
4187     $mask = 1 if not $options->{ "no_mask" };
4188
4189     $tmp_file = "$BP_TMP/write_2bit.fna";
4190     $fh_tmp   = Maasha::Common::write_open( $tmp_file );
4191
4192     $fh_out = write_stream( $options->{ "data_out" } );
4193
4194     while ( $record = get_record( $in ) )
4195     {
4196         if ( $entry = record2fasta( $record ) ) {
4197             Maasha::Fasta::put_entry( $entry, $fh_tmp );
4198         }
4199
4200         put_record( $record, $out ) if not $options->{ "no_stream" };
4201     }
4202
4203     close $fh_tmp;
4204
4205     $fh_in = Maasha::Common::read_open( $tmp_file );
4206
4207     Maasha::TwoBit::fasta2twobit( $fh_in, $fh_out, $mask );
4208
4209     close $fh_in;
4210     close $fh_out;
4211
4212     unlink $tmp_file;
4213 }
4214
4215
4216 sub script_write_solid
4217 {
4218     # Martin A. Hansen, April 2008.
4219
4220     # Write di-base encoded Solid sequence from entries in stream.
4221
4222     my ( $in,        # handle to in stream
4223          $out,       # handle to out stream
4224          $options,   # options hash
4225        ) = @_;
4226
4227     # Returns nothing.
4228
4229     my ( $record, $fh, $entry );
4230
4231     $fh = write_stream( $options->{ "data_out" }, $options->{ "compress" } );
4232
4233     while ( $record = get_record( $in ) )
4234     {
4235         if ( $entry = record2fasta( $record ) )
4236         {
4237             $entry->[ SEQ ] = Maasha::Solid::seq2color_space( uc $entry->[ SEQ ] );
4238
4239             Maasha::Fasta::put_entry( $entry, $fh, $options->{ "wrap" } );
4240         }
4241
4242         put_record( $record, $out ) if not $options->{ "no_stream" };
4243     }
4244
4245     close $fh;
4246 }
4247
4248
4249 sub script_plot_seqlogo
4250 {
4251     # Martin A. Hansen, August 2007.
4252
4253     # Calculates and writes a sequence logo for alignments.
4254
4255     my ( $in,        # handle to in stream
4256          $out,       # handle to out stream
4257          $options,   # options hash
4258        ) = @_;
4259
4260     # Returns nothing.
4261
4262     my ( $record, @entries, $logo, $fh );
4263
4264     while ( $record = get_record( $in ) )
4265     {
4266         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } ) {
4267             push @entries, [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ];
4268         }
4269
4270         put_record( $record, $out ) if not $options->{ "no_stream" };
4271     }
4272
4273     $logo = Maasha::Plot::seq_logo( \@entries );
4274
4275     $fh = write_stream( $options->{ "data_out" } );
4276
4277     print $fh $logo;
4278
4279     close $fh;
4280 }
4281
4282
4283 sub script_plot_phastcons_profiles
4284 {
4285     # Martin A. Hansen, January 2008.
4286
4287     # Plots PhastCons profiles.
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 ( $phastcons_file, $phastcons_index, $index, $fh_phastcons, $record, $scores, $AoA, $plot, $fh );
4297
4298     $options->{ "title" } ||= "PhastCons Profiles";
4299
4300     $phastcons_file  = Maasha::Config::genome_phastcons( $options->{ "genome" } );
4301     $phastcons_index = Maasha::Config::genome_phastcons_index( $options->{ "genome" } );
4302
4303     $index           = Maasha::UCSC::phastcons_index_retrieve( $phastcons_index );
4304     $fh_phastcons    = Maasha::Common::read_open( $phastcons_file );
4305
4306     while ( $record = get_record( $in ) )
4307     {
4308         if ( $record->{ "CHR" } and $record->{ "CHR_BEG" } and $record->{ "CHR_END" } )
4309         {
4310             $scores = Maasha::UCSC::phastcons_index_lookup( $index, $fh_phastcons, $record->{ "CHR" }, $record->{ "CHR_BEG" }, $record->{ "CHR_END" }, $options->{ "flank" } );
4311
4312             push @{ $AoA }, [ @{ $scores } ];
4313         }
4314
4315         put_record( $record, $out ) if not $options->{ "no_stream" };
4316     }
4317
4318     Maasha::UCSC::phastcons_normalize( $AoA );
4319
4320     $AoA = [ [ Maasha::UCSC::phastcons_mean( $AoA ) ] ] if $options->{ "mean" };
4321     $AoA = [ [ Maasha::UCSC::phastcons_median( $AoA ) ] ] if $options->{ "median" };
4322
4323     $AoA = Maasha::Matrix::matrix_flip( $AoA );
4324
4325     $plot = Maasha::Plot::lineplot_simple( $AoA, $options, $BP_TMP );
4326
4327     $fh = write_stream( $options->{ "data_out" } );
4328
4329     print $fh "$_\n" foreach @{ $plot };
4330
4331     close $fh;
4332 }
4333
4334
4335 sub script_analyze_bed
4336 {
4337     # Martin A. Hansen, March 2008.
4338
4339     # Analyze BED entries in stream.
4340
4341     my ( $in,        # handle to in stream
4342          $out,       # handle to out stream
4343          $options,   # options hash
4344        ) = @_;
4345
4346     # Returns nothing.
4347
4348     my ( $record );
4349
4350     while ( $record = get_record( $in ) )
4351     {
4352         $record = Maasha::UCSC::bed_analyze( $record ) if $record->{ "REC_TYPE" } eq "BED";
4353
4354         put_record( $record, $out );
4355     }
4356 }
4357
4358
4359 sub script_analyze_vals
4360 {
4361     # Martin A. Hansen, August 2007.
4362
4363     # Analyze values for given keys in stream.
4364
4365     my ( $in,        # handle to in stream
4366          $out,       # handle to out stream
4367          $options,   # options hash
4368        ) = @_;
4369
4370     # Returns nothing.
4371
4372     my ( $record, $key, @keys, %key_hash, $analysis, $len );
4373
4374     map { $key_hash{ $_ } = 1 } @{ $options->{ "keys" } };
4375
4376     while ( $record = get_record( $in ) )
4377     {
4378         foreach $key ( keys %{ $record } )
4379         {
4380             next if $options->{ "keys" } and not exists $key_hash{ $key };
4381
4382             $analysis->{ $key }->{ "COUNT" }++;
4383
4384             if ( Maasha::Calc::is_a_number( $record->{ $key } ) )
4385             {
4386                 $analysis->{ $key }->{ "TYPE" } = "num";
4387                 $analysis->{ $key }->{ "SUM" } += $record->{ $key };
4388                 $analysis->{ $key }->{ "MAX" } = $record->{ $key } if $record->{ $key } > $analysis->{ $key }->{ "MAX" } or not $analysis->{ $key }->{ "MAX" };
4389                 $analysis->{ $key }->{ "MIN" } = $record->{ $key } if $record->{ $key } < $analysis->{ $key }->{ "MIN" } or not $analysis->{ $key }->{ "MIN" };
4390             }
4391             else
4392             {
4393                 $len = length $record->{ $key };
4394
4395                 $analysis->{ $key }->{ "TYPE" } = "alph";
4396                 $analysis->{ $key }->{ "SUM" } += $len;
4397                 $analysis->{ $key }->{ "MAX" } = $len if $len > $analysis->{ $key }->{ "MAX" } or not $analysis->{ $key }->{ "MAX" };
4398                 $analysis->{ $key }->{ "MIN" } = $len if $len < $analysis->{ $key }->{ "MIM" } or not $analysis->{ $key }->{ "MIN" };
4399             }
4400         }
4401
4402         put_record( $record, $out ) if not $options->{ "no_stream" };
4403     }
4404
4405     foreach $key ( keys %{ $analysis } )
4406     {
4407         $analysis->{ $key }->{ "MEAN" } = sprintf "%.2f", $analysis->{ $key }->{ "SUM" } / $analysis->{ $key }->{ "COUNT" };
4408         $analysis->{ $key }->{ "SUM" }  = sprintf "%.2f", $analysis->{ $key }->{ "SUM" };
4409     }
4410
4411     my ( $keys, $types, $counts, $mins, $maxs, $sums, $means );
4412
4413     $keys   = "KEY  ";
4414     $types  = "TYPE ";
4415     $counts = "COUNT";
4416     $mins   = "MIN  ";
4417     $maxs   = "MAX  ";
4418     $sums   = "SUM  ";
4419     $means  = "MEAN ";
4420
4421     if ( $options->{ "keys" } ) {
4422         @keys = @{ $options->{ "keys" } };
4423     } else {
4424         @keys = keys %{ $analysis };
4425     }
4426
4427     foreach $key ( @keys )
4428     {
4429         $keys   .= sprintf "% 15s", $key;
4430         $types  .= sprintf "% 15s", $analysis->{ $key }->{ "TYPE" };
4431         $counts .= sprintf "% 15s", $analysis->{ $key }->{ "COUNT" };
4432         $mins   .= sprintf "% 15s", $analysis->{ $key }->{ "MIN" };
4433         $maxs   .= sprintf "% 15s", $analysis->{ $key }->{ "MAX" };
4434         $sums   .= sprintf "% 15s", $analysis->{ $key }->{ "SUM" };
4435         $means  .= sprintf "% 15s", $analysis->{ $key }->{ "MEAN" };
4436     }
4437
4438     print $out "$keys\n";
4439     print $out "$types\n";
4440     print $out "$counts\n";
4441     print $out "$mins\n";
4442     print $out "$maxs\n";
4443     print $out "$sums\n";
4444     print $out "$means\n";
4445 }
4446
4447
4448 sub script_head_records
4449 {
4450     # Martin A. Hansen, August 2007.
4451
4452     # Display the first sequences in stream.
4453
4454     my ( $in,        # handle to in stream
4455          $out,       # handle to out stream
4456          $options,   # options hash
4457        ) = @_;
4458
4459     # Returns nothing.
4460
4461     my ( $record, $count );
4462
4463     $options->{ "num" } ||= 10;
4464
4465     $count = 0;
4466
4467     while ( $record = get_record( $in ) )
4468     {
4469         $count++;
4470
4471         put_record( $record, $out );
4472
4473         last if $count == $options->{ "num" };
4474     }
4475 }
4476
4477
4478 sub script_remove_keys
4479 {
4480     # Martin A. Hansen, August 2007.
4481
4482     # Remove keys from stream.
4483
4484     my ( $in,        # handle to in stream
4485          $out,       # handle to out stream
4486          $options,   # options hash
4487        ) = @_;
4488
4489     # Returns nothing.
4490
4491     my ( $record, $new_record );
4492
4493     while ( $record = get_record( $in ) )
4494     {
4495         if ( $options->{ "keys" } )
4496         {
4497             map { delete $record->{ $_ } } @{ $options->{ "keys" } };
4498         }
4499         elsif ( $options->{ "save_keys" } )
4500         {
4501             map { $new_record->{ $_ } = $record->{ $_ } if exists $record->{ $_ } } @{ $options->{ "save_keys" } };
4502
4503             $record = $new_record;
4504         }
4505
4506         put_record( $record, $out ) if keys %{ $record };
4507     }
4508 }
4509
4510
4511 sub script_rename_keys
4512 {
4513     # Martin A. Hansen, August 2007.
4514
4515     # Rename keys in stream.
4516
4517     my ( $in,        # handle to in stream
4518          $out,       # handle to out stream
4519          $options,   # options hash
4520        ) = @_;
4521
4522     # Returns nothing.
4523
4524     my ( $record );
4525
4526     while ( $record = get_record( $in ) )
4527     {
4528         if ( exists $record->{ $options->{ "keys" }->[ 0 ] } )
4529         {
4530             $record->{ $options->{ "keys" }->[ 1 ] } = $record->{ $options->{ "keys" }->[ 0 ] };
4531
4532             delete $record->{ $options->{ "keys" }->[ 0 ] };
4533         }
4534
4535         put_record( $record, $out );
4536     }
4537 }
4538
4539
4540 sub script_uniq_vals
4541 {
4542     # Martin A. Hansen, August 2007.
4543
4544     # Find unique values in stream.
4545
4546     my ( $in,        # handle to in stream
4547          $out,       # handle to out stream
4548          $options,   # options hash
4549        ) = @_;
4550
4551     # Returns nothing.
4552
4553     my ( %hash, $record );
4554
4555     while ( $record = get_record( $in ) )
4556     {
4557         if ( $record->{ $options->{ "key" } } )
4558         {
4559             if ( not $hash{ $record->{ $options->{ "key" } } } and not $options->{ "invert" } )
4560             {
4561                 put_record( $record, $out );
4562
4563                 $hash{ $record->{ $options->{ "key" } } } = 1;
4564             }
4565             elsif ( $hash{ $record->{ $options->{ "key" } } } and $options->{ "invert" } )
4566             {
4567                 put_record( $record, $out );
4568             }
4569             else
4570             {
4571                 $hash{ $record->{ $options->{ "key" } } } = 1;
4572             }
4573         }
4574         else
4575         {
4576             put_record( $record, $out );
4577         }
4578     }
4579 }
4580
4581
4582 sub script_merge_vals
4583 {
4584     # Martin A. Hansen, August 2007.
4585
4586     # Rename keys in stream.
4587
4588     my ( $in,        # handle to in stream
4589          $out,       # handle to out stream
4590          $options,   # options hash
4591        ) = @_;
4592
4593     # Returns nothing.
4594
4595     my ( $record, @join, $i );
4596
4597     $options->{ "delimit" } ||= '_';
4598
4599     while ( $record = get_record( $in ) )
4600     {
4601         if ( exists $record->{ $options->{ "keys" }->[ 0 ] } )
4602         {
4603             @join = $record->{ $options->{ "keys" }->[ 0 ] };
4604
4605             for ( $i = 1; $i < @{ $options->{ "keys" } }; $i++ ) {
4606                 push @join, $record->{ $options->{ "keys" }->[ $i ] } if exists $record->{ $options->{ "keys" }->[ $i ] };
4607             }
4608
4609             $record->{ $options->{ "keys" }->[ 0 ] } = join $options->{ "delimit" }, @join;
4610         }
4611
4612         put_record( $record, $out );
4613     }
4614 }
4615
4616
4617 sub script_merge_records
4618 {
4619     # Martin A. Hansen, July 2008.
4620
4621     # Merges records in the stream based on identical values of two given keys.
4622
4623     my ( $in,        # handle to in stream
4624          $out,       # handle to out stream
4625          $options,   # options hash
4626        ) = @_;
4627
4628     # Returns nothing.
4629
4630     my ( $merge, $record, $file1, $file2, $fh1, $fh2, $key1, $key2, @keys1, @keys2, @vals1, @vals2,
4631          $num1, $num2, $num, $cmp, $i );
4632
4633     $merge = $options->{ "merge" } || "AandB";
4634
4635     $file1 = "$BP_TMP/merge_records1.tmp";
4636     $file2 = "$BP_TMP/merge_records2.tmp";
4637
4638     $fh1   = Maasha::Common::write_open( $file1 );
4639     $fh2   = Maasha::Common::write_open( $file2 );
4640
4641     $key1  = $options->{ "keys" }->[ 0 ];
4642     $key2  = $options->{ "keys" }->[ 1 ];
4643
4644     $num   = $key2 =~ s/n$//;
4645     $num1  = 0;
4646     $num2  = 0;
4647
4648     while ( $record = get_record( $in ) )
4649     {
4650         if ( exists $record->{ $key1 } )
4651         {
4652             @keys1 = $key1;
4653             @vals1 = $record->{ $key1 };
4654
4655             delete $record->{ $key1 };
4656
4657             map { push @keys1, $_; push @vals1, $record->{ $_ } } keys %{ $record };
4658
4659             print $fh1 join( "\t", @vals1 ), "\n";
4660
4661             $num1++;
4662         }
4663         elsif ( exists $record->{ $key2 } )
4664         {
4665             @keys2 = $key2;
4666             @vals2 = $record->{ $key2 };
4667
4668             delete $record->{ $key2 };
4669
4670             map { push @keys2, $_; push @vals2, $record->{ $_ } } keys %{ $record };
4671
4672             print $fh2 join( "\t", @vals2 ), "\n";
4673
4674             $num2++;
4675         }
4676     }
4677
4678     close $fh1;
4679     close $fh2;
4680
4681     if ( $num )
4682     {
4683         Maasha::Common::run( "sort", "-k 1,1n $file1 > $file1.sort" ) and rename "$file1.sort", $file1;
4684         Maasha::Common::run( "sort", "-k 1,1n $file2 > $file2.sort" ) and rename "$file2.sort", $file2;
4685     }
4686     else
4687     {
4688         Maasha::Common::run( "sort", "-k 1,1 $file1 > $file1.sort" ) and rename "$file1.sort", $file1;
4689         Maasha::Common::run( "sort", "-k 1,1 $file2 > $file2.sort" ) and rename "$file2.sort", $file2;
4690     }
4691
4692     $fh1 = Maasha::Common::read_open( $file1 );
4693     $fh2 = Maasha::Common::read_open( $file2 );
4694
4695     @vals1 = Maasha::Common::get_fields( $fh1 );
4696     @vals2 = Maasha::Common::get_fields( $fh2 );
4697
4698     while ( $num1 > 0 and $num2 > 0 )
4699     {
4700         undef $record;
4701
4702         if ( $num ) {
4703             $cmp = $vals1[ 0 ] <=> $vals2[ 0 ];
4704         } else {
4705             $cmp = $vals1[ 0 ] cmp $vals2[ 0 ];
4706         }
4707
4708         if ( $cmp < 0 )
4709         {
4710             if ( $merge =~ /^(AorB|AnotB)$/ )
4711             {
4712                 for ( $i = 0; $i < @keys1; $i++ ) {
4713                     $record->{ $keys1[ $i ] } = $vals1[ $i ];
4714                 }
4715
4716                 put_record( $record, $out );
4717             }
4718
4719             @vals1 = Maasha::Common::get_fields( $fh1 );
4720             $num1--;
4721         }
4722         elsif ( $cmp > 0 )
4723         {
4724             if ( $merge =~ /^(BorA|BnotA)$/ )
4725             {
4726                 for ( $i = 0; $i < @keys2; $i++ ) {
4727                     $record->{ $keys2[ $i ] } = $vals2[ $i ];
4728                 }
4729
4730                 put_record( $record, $out );
4731             }
4732
4733             @vals2 = Maasha::Common::get_fields( $fh2 );
4734             $num2--;
4735         }
4736         else
4737         {
4738             if ( $merge =~ /^(AandB|AorB|BorA)$/ )
4739             {
4740                 for ( $i = 0; $i < @keys1; $i++ ) {
4741                     $record->{ $keys1[ $i ] } = $vals1[ $i ];
4742                 }
4743
4744                 for ( $i = 1; $i < @keys2; $i++ ) {
4745                     $record->{ $keys2[ $i ] } = $vals2[ $i ];
4746                 }
4747
4748                 put_record( $record, $out );
4749             }
4750
4751             @vals1 = Maasha::Common::get_fields( $fh1 );
4752             @vals2 = Maasha::Common::get_fields( $fh2 );
4753             $num1--;
4754             $num2--;
4755         }
4756     }
4757
4758     close $fh1;
4759     close $fh2;
4760
4761     unlink $file1;
4762     unlink $file2;
4763
4764     if ( $num1 > 0 and $merge =~ /^(AorB|AnotB)$/ )
4765     {
4766         undef $record;
4767
4768         for ( $i = 0; $i < @keys1; $i++ ) {
4769             $record->{ $keys1[ $i ] } = $vals1[ $i ];
4770         }
4771
4772         put_record( $record, $out );
4773     }
4774
4775     if ( $num2 > 0 and $merge =~ /^(BorA|BnotA)$/ )
4776     {
4777         undef $record;
4778
4779         for ( $i = 0; $i < @keys2; $i++ ) {
4780             $record->{ $keys2[ $i ] } = $vals2[ $i ];
4781         }
4782
4783         put_record( $record, $out );
4784     }
4785 }
4786
4787
4788 sub script_grab
4789 {
4790     # Martin A. Hansen, August 2007.
4791
4792     # Grab for records in stream.
4793
4794     my ( $in,        # handle to in stream
4795          $out,       # handle to out stream
4796          $options,   # options hash
4797        ) = @_;
4798
4799     # Returns nothing.
4800
4801     my ( $patterns, $pattern, $record, $key, $pos, $op, $val, %lookup_hash );
4802
4803     if ( $options->{ "patterns" } )
4804     {
4805         $patterns = [ split ",", $options->{ "patterns" } ];
4806     }
4807     elsif ( -f $options->{ "patterns_in" } )
4808     {
4809         $patterns = Maasha::Patscan::read_patterns( $options->{ "patterns_in" } );
4810     }
4811     elsif ( -f $options->{ "exact_in" } )
4812     {
4813         $patterns = Maasha::Patscan::read_patterns( $options->{ "exact_in" } );
4814
4815         map { $lookup_hash{ $_ } = 1 } @{ $patterns };
4816
4817         undef $patterns;
4818     }
4819
4820     if ( $options->{ "eval" } )
4821     {
4822         if ( $options->{ "eval" } =~ /^([^><=! ]+)\s*(>=|<=|>|<|=|!=|eq|ne)\s*(.+)$/ )
4823         {
4824             $key = $1;
4825             $op  = $2;
4826             $val = $3;
4827         }
4828     }
4829
4830     while ( $record = get_record( $in ) )
4831     {
4832         $pos = -1;
4833
4834         if ( %lookup_hash )
4835         {
4836             if ( $options->{ "keys" } )
4837             {
4838                 foreach $key ( @{ $options->{ "keys" } } )
4839                 {
4840                     if ( exists $lookup_hash{ $record->{ $key } } )
4841                     {
4842                         $pos = 1;
4843                         goto FOUND;
4844                     }
4845                 }
4846             }
4847             else
4848             {
4849                 foreach $key ( keys %{ $record } )
4850                 {
4851                     if ( not $options->{ "vals_only" } )
4852                     {
4853                         if ( exists $lookup_hash{ $key } )
4854                         {
4855                             $pos = 1;
4856                             goto FOUND;
4857                         }
4858                     }
4859
4860                     if ( not $options->{ "keys_only" } )
4861                     {
4862                         if ( exists $lookup_hash{ $record->{ $key } } )
4863                         {
4864                             $pos = 1;
4865                             goto FOUND;
4866                         }
4867                     }
4868                 }
4869             }
4870         }
4871         elsif ( $patterns )
4872         {
4873             foreach $pattern ( @{ $patterns } )
4874             {
4875                 if ( $options->{ "keys" } )
4876                 {
4877                     foreach $key ( @{ $options->{ "keys" } } )
4878                     {
4879                         $pos = index $record->{ $key }, $pattern;
4880
4881                         goto FOUND if $pos >= 0;
4882                     }
4883                 }
4884                 else
4885                 {
4886                     foreach $key ( keys %{ $record } )
4887                     {
4888                         if ( not $options->{ "vals_only" } )
4889                         {
4890                             $pos = index $key, $pattern;
4891
4892                             goto FOUND if $pos >= 0;
4893                         }
4894
4895                         if ( not $options->{ "keys_only" } )
4896                         {
4897                             $pos = index $record->{ $key }, $pattern;
4898
4899                             goto FOUND if $pos >= 0;
4900                         }
4901                     }
4902                 }
4903             }
4904         }
4905         elsif ( $options->{ "regex" } )
4906         {
4907             if ( $options->{ "keys" } )
4908             {
4909                 foreach $key ( @{ $options->{ "keys" } } )
4910                 {
4911                     if ( $options->{ "case_insensitive" } ) {
4912                         $pos = 1 if $record->{ $key } =~ /$options->{'regex'}/i;
4913                     } else {
4914                         $pos = 1 if $record->{ $key } =~ /$options->{'regex'}/;
4915                     }
4916
4917                     goto FOUND if $pos >= 0;
4918                 }
4919             }
4920             else
4921             {
4922                 foreach $key ( keys %{ $record } )
4923                 {
4924                     if ( not $options->{ "vals_only" } )
4925                     {
4926                         if ( $options->{ "case_insensitive" } ) {
4927                             $pos = 1 if $key =~ /$options->{'regex'}/i;
4928                         } else {
4929                             $pos = 1 if $key =~ /$options->{'regex'}/;
4930                         }
4931
4932                         goto FOUND if $pos >= 0;
4933                     }
4934
4935                     if ( not $options->{ "keys_only" } )
4936                     {
4937                         if ( $options->{ "case_insensitive" } ) {
4938                             $pos = 1 if $record->{ $key } =~ /$options->{'regex'}/i;
4939                         } else {
4940                             $pos = 1 if $record->{ $key } =~ /$options->{'regex'}/;
4941                         }
4942
4943                         goto FOUND if $pos >= 0;
4944                     }
4945                 }
4946             }
4947         }
4948         elsif ( $options->{ "eval" } )
4949         {
4950             if ( defined $record->{ $key } )
4951             {
4952                 if ( $op eq "<" and $record->{ $key } < $val ) {
4953                     $pos = 1 and goto FOUND;
4954                 } elsif ( $op eq ">" and $record->{ $key } > $val ) {
4955                     $pos = 1 and goto FOUND;
4956                 } elsif ( $op eq ">=" and $record->{ $key } >= $val ) {
4957                     $pos = 1 and goto FOUND;
4958                 } elsif ( $op eq "<=" and $record->{ $key } <= $val ) {
4959                     $pos = 1 and goto FOUND;
4960                 } elsif ( $op eq "=" and $record->{ $key } == $val ) {
4961                     $pos = 1 and goto FOUND;
4962                 } elsif ( $op eq "!=" and $record->{ $key } != $val ) {
4963                     $pos = 1 and goto FOUND;
4964                 } elsif ( $op eq "eq" and $record->{ $key } eq $val ) {
4965                     $pos = 1 and goto FOUND;
4966                 } elsif ( $op eq "ne" and $record->{ $key } ne $val ) {
4967                     $pos = 1 and goto FOUND;
4968                 }
4969             }
4970         }
4971
4972         FOUND:
4973
4974         if ( $pos >= 0 and not $options->{ "invert" } ) {
4975             put_record( $record, $out );
4976         } elsif ( $pos < 0 and $options->{ "invert" } ) {
4977             put_record( $record, $out );
4978         }
4979     }
4980 }
4981
4982
4983 sub script_compute
4984 {
4985     # Martin A. Hansen, August 2007.
4986
4987     # Evaluate extression for records in stream.
4988
4989     my ( $in,        # handle to in stream
4990          $out,       # handle to out stream
4991          $options,   # options hash
4992        ) = @_;
4993
4994     # Returns nothing.
4995
4996     my ( $record, $eval_key, $eval_val, $check, @keys );
4997
4998     while ( $record = get_record( $in ) )
4999     {
5000         if ( $options->{ "eval" } )
5001         {
5002             if ( $options->{ "eval" } =~ /^(.+)\s*=\s*(.+)$/ )
5003             {
5004                 $eval_key = $1;
5005                 $eval_val = $2;
5006             }
5007
5008             if ( not $check )
5009             {
5010                 @keys = split /\W+/, $eval_val;
5011                 @keys = grep { ! /^\d+$/ } @keys;
5012
5013                 $check = 1;
5014             }
5015
5016             map { $eval_val =~ s/$_/$record->{ $_ }/g } @keys;
5017
5018             $record->{ $eval_key } = eval "$eval_val" or Maasha::Common::error( "eval failed -> $@" );
5019         }
5020
5021         put_record( $record, $out );
5022     }
5023 }
5024
5025
5026 sub script_flip_tab
5027 {
5028     # Martin A. Hansen, June 2008.
5029
5030     # Flip a table.
5031
5032     my ( $in,        # handle to in stream
5033          $out,       # handle to out stream
5034          $options,   # options hash
5035        ) = @_;
5036
5037     # Returns nothing.
5038
5039     my ( $record, $key, $A, $B, @rows, @matrix, $row, $i );
5040
5041     while ( $record = get_record( $in ) )
5042     {
5043         undef @rows;
5044
5045         foreach $key ( sort { $A = $a; $B = $b; $A =~ s/^V(\d+)$/$1/; $B =~ s/^V(\d+)$/$1/; $A <=> $B } keys %{ $record } )
5046         {
5047             push @rows, $record->{ $key };
5048
5049         }
5050
5051         push @matrix, [ @rows ];
5052     }
5053
5054     undef $record;
5055
5056     @matrix = Maasha::Matrix::matrix_flip( \@matrix );
5057
5058     foreach $row ( @matrix )
5059     {
5060         for ( $i = 0; $i < @{ $row }; $i++ ) {
5061             $record->{ "V$i" } = $row->[ $i ];
5062         }
5063
5064         put_record( $record, $out );
5065     }
5066 }
5067
5068
5069 sub script_add_ident
5070 {
5071     # Martin A. Hansen, May 2008.
5072
5073     # Add a unique identifier to each record in stream.
5074
5075     my ( $in,        # handle to in stream
5076          $out,       # handle to out stream
5077          $options,   # options hash
5078        ) = @_;
5079
5080     # Returns nothing.
5081
5082     my ( $record, $key, $prefix, $i );
5083
5084     $key    = $options->{ "key" }    || "ID";
5085     $prefix = $options->{ "prefix" } || "ID";
5086
5087     $i = 0;
5088
5089     while ( $record = get_record( $in ) )
5090     {
5091         $record->{ $key } = sprintf( "$prefix%08d", $i );
5092
5093         put_record( $record, $out );
5094
5095         $i++;
5096     }
5097 }
5098
5099
5100 sub script_count_records
5101 {
5102     # Martin A. Hansen, August 2007.
5103
5104     # Count records in stream.
5105
5106     my ( $in,        # handle to in stream
5107          $out,       # handle to out stream
5108          $options,   # options hash
5109        ) = @_;
5110
5111     # Returns nothing.
5112
5113     my ( $record, $count, $result, $fh, $line );
5114
5115     $count = 0;
5116
5117     if ( $options->{ "no_stream" } )
5118     {
5119         while ( $line = <$in> )
5120         {
5121             chomp $line;
5122
5123             $count++ if $line eq "---";
5124         }
5125     }
5126     else
5127     {
5128         while ( $record = get_record( $in ) )
5129         {
5130             put_record( $record, $out );
5131
5132             $count++;
5133         }
5134     }
5135
5136     $result = { "RECORDS_COUNT" => $count };
5137
5138     $fh = write_stream( $options->{ "data_out" } );
5139
5140     put_record( $result, $fh );
5141
5142     close $fh;
5143 }
5144
5145
5146 sub script_random_records
5147 {
5148     # Martin A. Hansen, August 2007.
5149
5150     # Pick a number or random records from stream.
5151
5152     my ( $in,        # handle to in stream
5153          $out,       # handle to out stream
5154          $options,   # options hash
5155        ) = @_;
5156
5157     # Returns nothing.
5158
5159     my ( $record, $tmp_file, $fh_out, $fh_in, $count, $i, %rand_hash, $rand, $max );
5160
5161     $options->{ "num" } ||= 10;
5162
5163     $tmp_file = "$BP_TMP/random_records.tmp";
5164
5165     $fh_out = Maasha::Common::write_open( $tmp_file );
5166
5167     $count = 0;
5168
5169     while ( $record = get_record( $in ) )
5170     {
5171         put_record( $record, $fh_out );
5172
5173         $count++;
5174     }
5175
5176     close $fh_out;
5177
5178     $max = 0;
5179     $i   = 0;
5180
5181     Maasha::Common::error( qq(Requested random records > records in stream) ) if $options->{ "num" } > $count;
5182
5183     while ( $i < $options->{ "num" } )
5184     {
5185         $rand = int( rand( $count ) );
5186
5187         if ( not exists $rand_hash{ $rand } )
5188         {
5189             $rand_hash{ $rand } = 1;
5190
5191             $max = $rand if $rand > $max;
5192
5193             $i++;
5194         }
5195     }
5196
5197     $fh_in = Maasha::Common::read_open( $tmp_file );
5198
5199     $count = 0;
5200
5201     while ( $record = get_record( $fh_in ) )
5202     {
5203         put_record( $record, $out ) if exists $rand_hash{ $count };
5204
5205         last if $count == $max;
5206
5207         $count++;
5208     }
5209
5210     close $fh_in;
5211
5212     unlink $tmp_file;
5213 }
5214
5215
5216 sub script_sort_records
5217 {
5218     # Martin A. Hansen, August 2007.
5219
5220     # Sort to sort records according to keys.
5221
5222     my ( $in,        # handle to in stream
5223          $out,       # handle to out stream
5224          $options,   # options hash
5225        ) = @_;
5226
5227     # Returns nothing.
5228
5229     my ( @keys, $key, @sort_cmd, $sort_str, $sort_sub, @records, $record, $i );
5230
5231     foreach $key ( @{ $options->{ "keys" } } )
5232     {
5233         if ( $key =~ s/n$// ) {
5234             push @sort_cmd, qq(\$a->{ "$key" } <=> \$b->{ "$key" });
5235         } else {
5236             push @sort_cmd, qq(\$a->{ "$key" } cmp \$b->{ "$key" });
5237         }
5238     }
5239
5240     $sort_str = join " or ", @sort_cmd;
5241     $sort_sub = eval "sub { $sort_str }";   # NB security issue!
5242
5243     while ( $record = get_record( $in ) ) {
5244         push @records, $record;
5245     }
5246
5247     @records = sort $sort_sub @records;
5248
5249     if ( $options->{ "reverse" } )
5250     {
5251         for ( $i = scalar @records - 1; $i >= 0; $i-- ) {
5252             put_record( $records[ $i ], $out );
5253         }
5254     }
5255     else
5256     {
5257         for ( $i = 0; $i < scalar @records; $i++ ) {
5258             put_record( $records[ $i ], $out );
5259         }
5260     }
5261 }
5262
5263
5264 sub script_count_vals
5265 {
5266     # Martin A. Hansen, August 2007.
5267
5268     # Count records in stream.
5269
5270     my ( $in,        # handle to in stream
5271          $out,       # handle to out stream
5272          $options,   # options hash
5273        ) = @_;
5274
5275     # Returns nothing.
5276
5277     my ( $num, $record, %count_hash, @records, $tmp_file, $fh_out, $fh_in, $cache );
5278
5279     $tmp_file = "$BP_TMP/count_cache.tmp";
5280
5281     $fh_out   = Maasha::Common::write_open( $tmp_file );
5282
5283     $cache    = 0;
5284     $num      = 0;
5285
5286     while ( $record = get_record( $in ) )
5287     {
5288         map { $count_hash{ $_ }{ $record->{ $_ } }++ if exists $record->{ $_ } } @{ $options->{ "keys" } };
5289
5290         push @records, $record;
5291
5292         if ( scalar @records > 5_000_000 )   # too many records to hold in memory - use disk cache
5293         {
5294             map { put_record( $_, $fh_out ) } @records;
5295
5296             undef @records;
5297
5298             $cache = 1;
5299         }
5300
5301         print STDERR "verbose: records read $num\n" if ( $options->{ 'verbose' } and ( $num % 1_000_000 ) == 0 );
5302
5303         $num++;
5304     }
5305
5306     close $fh_out;
5307
5308     if ( $cache )
5309     {
5310         $num   = 0;
5311
5312         $fh_in = Maasha::Common::read_open( $tmp_file );
5313
5314         while ( $record = get_record( $fh_in ) )
5315         {
5316             map { $record->{ $_ . "_COUNT" } = $count_hash{ $_ }{ $record->{ $_ } } if exists $record->{ $_ } } @{ $options->{ "keys" } };
5317
5318             put_record( $record, $out );
5319
5320             print STDERR "verbose: cache read $num\n" if ( $options->{ 'verbose' } and ( $num % 1_000_000 ) == 0 );
5321
5322             $num++;
5323         }
5324
5325         close $fh_in;
5326     }
5327
5328     foreach $record ( @records )
5329     {
5330         map { $record->{ $_ . "_COUNT" } = $count_hash{ $_ }{ $record->{ $_ } } if exists $record->{ $_ } } @{ $options->{ "keys" } };
5331
5332         put_record( $record, $out );
5333     }
5334
5335     unlink $tmp_file;
5336 }
5337
5338
5339 sub script_plot_histogram
5340 {
5341     # Martin A. Hansen, September 2007.
5342
5343     # Plot a simple histogram for a given key using GNU plot.
5344
5345     my ( $in,        # handle to in stream
5346          $out,       # handle to out stream
5347          $options,   # options hash
5348        ) = @_;
5349
5350     # Returns nothing.
5351
5352     my ( $record, %data_hash, $max, @data_list, $i, $result, $fh );
5353
5354     $options->{ "title" } ||= "Histogram";
5355     $options->{ "sort" }  ||= "num";
5356
5357     while ( $record = get_record( $in ) )
5358     {
5359         $data_hash{ $record->{ $options->{ "key" } } }++ if $record->{ $options->{ "key" } };
5360
5361         put_record( $record, $out ) if not $options->{ "no_stream" };
5362     }
5363
5364     if ( $options->{ "sort" } eq "num" ) {
5365         map { push @data_list, [ $_, $data_hash{ $_ } ] } sort { $a <=> $b } keys %data_hash;
5366     } else {
5367         map { push @data_list, [ $_, $data_hash{ $_ } ] } sort keys %data_hash;
5368     }
5369
5370     $result = Maasha::Plot::histogram_simple( \@data_list, $options );
5371
5372     $fh = write_stream( $options->{ "data_out" } );
5373
5374     print $fh "$_\n" foreach @{ $result };
5375
5376     close $fh;
5377 }
5378
5379
5380 sub script_plot_lendist
5381 {
5382     # Martin A. Hansen, August 2007.
5383
5384     # Plot length distribution using GNU plot.
5385
5386     my ( $in,        # handle to in stream
5387          $out,       # handle to out stream
5388          $options,   # options hash
5389        ) = @_;
5390
5391     # Returns nothing.
5392
5393     my ( $record, %data_hash, $max, @data_list, $i, $result, $fh );
5394
5395     $options->{ "title" } ||= "Length Distribution";
5396
5397     while ( $record = get_record( $in ) )
5398     {
5399         $data_hash{ $record->{ $options->{ "key" } } }++ if $record->{ $options->{ "key" } };
5400
5401         put_record( $record, $out ) if not $options->{ "no_stream" };
5402     }
5403
5404     $max = Maasha::Calc::list_max( [ keys %data_hash ] );
5405
5406     for ( $i = 0; $i < $max; $i++ ) {
5407         push @data_list, [ $i, $data_hash{ $i } || 0 ];
5408     }
5409
5410     $result = Maasha::Plot::histogram_lendist( \@data_list, $options );
5411
5412     $fh = write_stream( $options->{ "data_out" } );
5413
5414     print $fh "$_\n" foreach @{ $result };
5415
5416     close $fh;
5417 }
5418
5419
5420 sub script_plot_chrdist
5421 {
5422     # Martin A. Hansen, August 2007.
5423
5424     # Plot chromosome distribution using GNU plot.
5425
5426     my ( $in,        # handle to in stream
5427          $out,       # handle to out stream
5428          $options,   # options hash
5429        ) = @_;
5430
5431     # Returns nothing.
5432
5433     my ( $record, %data_hash, @data_list, $elem, $sort_key, $count, $result, $fh );
5434
5435     $options->{ "title" } ||= "Chromosome Distribution";
5436
5437     while ( $record = get_record( $in ) )
5438     {
5439         if ( $record->{ "CHR" } ) {                                                             # generic
5440             $data_hash{ $record->{ "CHR" } }++;
5441         } elsif ( $record->{ "REC_TYPE" } eq "PATSCAN" and $record->{ "S_ID" } =~ /^chr/i ) {   # patscan
5442             $data_hash{ $record->{ "S_ID" } }++;
5443         } elsif ( $record->{ "REC_TYPE" } eq "PSL" and $record->{ "S_ID" } =~ /^chr/i ) {       # BLAT / PSL
5444             $data_hash{ $record->{ "S_ID" } }++;
5445         } elsif ( $record->{ "REC_TYPE" } eq "BLAST" and $record->{ "S_ID" } =~ /^chr/i ) {     # BLAST
5446             $data_hash{ $record->{ "S_ID" } }++;
5447         }
5448
5449         put_record( $record, $out ) if not $options->{ "no_stream" };
5450     }
5451
5452     foreach $elem ( keys %data_hash )
5453     {
5454         $sort_key = $elem;
5455
5456         $sort_key =~ s/chr//i;
5457
5458         $sort_key =~ s/^X(.*)/99$1/;
5459         $sort_key =~ s/^Y(.*)/99$1/;
5460         $sort_key =~ s/^Z(.*)/999$1/;
5461         $sort_key =~ s/^M(.*)/9999$1/;
5462         $sort_key =~ s/^U(.*)/99999$1/;
5463
5464         $count = $sort_key =~ tr/_//;
5465
5466         $sort_key =~ s/_.*/"999999" x $count/ex;
5467
5468         push @data_list, [ $elem, $data_hash{ $elem }, $sort_key ];
5469     }
5470
5471     @data_list = sort { $a->[ 2 ] <=> $b->[ 2 ] } @data_list;
5472
5473     $result = Maasha::Plot::histogram_chrdist( \@data_list, $options );
5474
5475     $fh = write_stream( $options->{ "data_out" } );
5476
5477     print $fh "$_\n" foreach @{ $result };
5478
5479     close $fh;
5480 }
5481
5482
5483 sub script_plot_karyogram
5484 {
5485     # Martin A. Hansen, August 2007.
5486
5487     # Plot hits on karyogram.
5488
5489     my ( $in,        # handle to in stream
5490          $out,       # handle to out stream
5491          $options,   # options hash
5492        ) = @_;
5493
5494     # Returns nothing.
5495
5496     my ( %options, $record, @data, $fh, $result, %data_hash );
5497
5498     $options->{ "genome" }     ||= "human";
5499     $options->{ "feat_color" } ||= "black";
5500
5501     while ( $record = get_record( $in ) )
5502     {
5503         if ( $record->{ "CHR" } and $record->{ "CHR_BEG" } and $record->{ "CHR_END" } )
5504         {
5505             push @{ $data_hash{ $record->{ "CHR" } } }, [ $record->{ "CHR_BEG" }, $record->{ "CHR_END" }, $options->{ "feat_color" } ];
5506         }
5507
5508         put_record( $record, $out ) if not $options->{ "no_stream" };
5509     }
5510
5511     $result = Maasha::Plot::karyogram( \%data_hash, \%options );
5512
5513     $fh = write_stream( $options->{ "data_out" } );
5514
5515     print $fh $result;
5516
5517     close $fh;
5518 }
5519
5520
5521 sub script_plot_matches
5522 {
5523     # Martin A. Hansen, August 2007.
5524
5525     # Plot matches in 2D generating a dotplot.
5526
5527     my ( $in,        # handle to in stream
5528          $out,       # handle to out stream
5529          $options,   # options hash
5530        ) = @_;
5531
5532     # Returns nothing.
5533
5534     my ( $record, @data, $fh, $result, %data_hash );
5535
5536     $options->{ "direction" } ||= "both";
5537
5538     while ( $record = get_record( $in ) )
5539     {
5540         if ( defined $record->{ "Q_BEG" } and defined $record->{ "S_BEG" } and $record->{ "Q_END" } and $record->{ "S_END" } ) {
5541             push @data, $record;
5542         }
5543
5544         put_record( $record, $out ) if not $options->{ "no_stream" };
5545     }
5546
5547     $options->{ "title" }  ||= "plot_matches";
5548     $options->{ "xlabel" } ||= $data[ 0 ]->{ "Q_ID" };
5549     $options->{ "ylabel" } ||= $data[ 0 ]->{ "S_ID" };
5550
5551     $result = Maasha::Plot::dotplot_matches( \@data, $options, $BP_TMP );
5552
5553     $fh = write_stream( $options->{ "data_out" } );
5554
5555     print $fh "$_\n" foreach @{ $result };
5556
5557     close $fh;
5558 }
5559
5560
5561 sub script_length_vals
5562 {
5563     # Martin A. Hansen, August 2007.
5564
5565     # Determine the length of the value for given keys.
5566
5567     my ( $in,        # handle to in stream
5568          $out,       # handle to out stream
5569          $options,   # options hash
5570        ) = @_;
5571
5572     # Returns nothing.
5573
5574     my ( $record, $key );
5575
5576     while ( $record = get_record( $in ) )
5577     {
5578         foreach $key ( @{ $options->{ "keys" } } )
5579         {
5580             if ( $record->{ $key } ) {
5581                 $record->{ $key . "_LEN" } = length $record->{ $key };
5582             }
5583         }
5584
5585         put_record( $record, $out );
5586     }
5587 }
5588
5589
5590 sub script_sum_vals
5591 {
5592     # Martin A. Hansen, August 2007.
5593
5594     # Calculates the sums for values of given keys.
5595
5596     my ( $in,        # handle to in stream
5597          $out,       # handle to out stream
5598          $options,   # options hash
5599        ) = @_;
5600
5601     # Returns nothing.
5602
5603     my ( $record, $key, %sum_hash, $fh );
5604
5605     while ( $record = get_record( $in ) )
5606     {
5607         foreach $key ( @{ $options->{ "keys" } } )
5608         {
5609             if ( $record->{ $key } ) {
5610                 $sum_hash{ $key } += $record->{ $key };
5611             }
5612         }
5613
5614         put_record( $record, $out ) if not $options->{ "no_stream" };
5615     }
5616
5617     $fh = write_stream( $options->{ "data_out" } );
5618
5619     foreach $key ( @{ $options->{ "keys" } } ) {
5620         put_record( { $key . "_SUM" => $sum_hash{ $key } || 0 } , $fh );
5621     }
5622
5623     close $fh;
5624 }
5625
5626
5627 sub script_mean_vals
5628 {
5629     # Martin A. Hansen, August 2007.
5630
5631     # Calculate the mean of values of given keys.
5632
5633     my ( $in,        # handle to in stream
5634          $out,       # handle to out stream
5635          $options,   # options hash
5636        ) = @_;
5637
5638     # Returns nothing.
5639
5640     my ( $record, $key, %sum_hash, %count_hash, $mean, $fh );
5641
5642     while ( $record = get_record( $in ) )
5643     {
5644         foreach $key ( @{ $options->{ "keys" } } )
5645         {
5646             if ( $record->{ $key } )
5647             {
5648                 $sum_hash{ $key } += $record->{ $key };
5649                 $count_hash{ $key }++;
5650             }
5651         }
5652
5653         put_record( $record, $out ) if not $options->{ "no_stream" };
5654     }
5655
5656     $fh = write_stream( $options->{ "data_out" } );
5657
5658     foreach $key ( @{ $options->{ "keys" } } )
5659     {
5660         if ( $count_hash{ $key } ) {
5661             $mean = sprintf( "%.2f", ( $sum_hash{ $key } / $count_hash{ $key } ) );
5662         } else {
5663             $mean = "N/A";
5664         }
5665
5666         put_record( { $key . "_MEAN" => $mean } , $fh );
5667     }
5668
5669     close $fh;
5670 }
5671
5672
5673 sub script_median_vals
5674 {
5675     # Martin A. Hansen, March 2008.
5676
5677     # Calculate the median values of given keys.
5678
5679     my ( $in,        # handle to in stream
5680          $out,       # handle to out stream
5681          $options,   # options hash
5682        ) = @_;
5683
5684     # Returns nothing.
5685
5686     my ( $record, $key, %median_hash, $median, $fh );
5687
5688     while ( $record = get_record( $in ) )
5689     {
5690         foreach $key ( @{ $options->{ "keys" } } ) {
5691             push @{ $median_hash{ $key } }, $record->{ $key } if defined $record->{ $key };
5692         }
5693
5694         put_record( $record, $out ) if not $options->{ "no_stream" };
5695     }
5696
5697     $fh = write_stream( $options->{ "data_out" } );
5698
5699     foreach $key ( @{ $options->{ "keys" } } )
5700     {
5701         if ( $median_hash{ $key } ) {
5702             $median = Maasha::Calc::median( $median_hash{ $key } );
5703         } else {
5704             $median = "N/A";
5705         }
5706
5707         put_record( { $key . "_MEDIAN" => $median } , $fh );
5708     }
5709
5710     close $fh;
5711 }
5712
5713
5714 sub script_max_vals
5715 {
5716     # Martin A. Hansen, February 2008.
5717
5718     # Determine the maximum values of given keys.
5719
5720     my ( $in,        # handle to in stream
5721          $out,       # handle to out stream
5722          $options,   # options hash
5723        ) = @_;
5724
5725     # Returns nothing.
5726
5727     my ( $record, $key, $fh, %max_hash, $max_record );
5728
5729     while ( $record = get_record( $in ) )
5730     {
5731         foreach $key ( @{ $options->{ "keys" } } )
5732         {
5733             if ( $record->{ $key } )
5734             {
5735                 $max_hash{ $key } = $record->{ $key } if $record->{ $key } > $max_hash{ $key };
5736             }
5737         }
5738
5739         put_record( $record, $out ) if not $options->{ "no_stream" };
5740     }
5741
5742     $fh = write_stream( $options->{ "data_out" } );
5743
5744     foreach $key ( @{ $options->{ "keys" } } )
5745     {
5746         $max_record->{ $key . "_MAX" } = $max_hash{ $key };
5747     }
5748
5749     put_record( $max_record, $fh );
5750
5751     close $fh;
5752 }
5753
5754
5755 sub script_min_vals
5756 {
5757     # Martin A. Hansen, February 2008.
5758
5759     # Determine the minimum values of given keys.
5760
5761     my ( $in,        # handle to in stream
5762          $out,       # handle to out stream
5763          $options,   # options hash
5764        ) = @_;
5765
5766     # Returns nothing.
5767
5768     my ( $record, $key, $fh, %min_hash, $min_record );
5769
5770     while ( $record = get_record( $in ) )
5771     {
5772         foreach $key ( @{ $options->{ "keys" } } )
5773         {
5774             if ( defined $record->{ $key } )
5775             {
5776                 if ( exists $min_hash{ $key } ) {
5777                     $min_hash{ $key } = $record->{ $key } if $record->{ $key } < $min_hash{ $key };
5778                 } else {
5779                     $min_hash{ $key } = $record->{ $key };
5780                 }
5781             }
5782         }
5783
5784         put_record( $record, $out ) if not $options->{ "no_stream" };
5785     }
5786
5787     $fh = write_stream( $options->{ "data_out" } );
5788
5789     foreach $key ( @{ $options->{ "keys" } } )
5790     {
5791         $min_record->{ $key . "_MIN" } = $min_hash{ $key };
5792     }
5793
5794     put_record( $min_record, $fh );
5795
5796     close $fh;
5797 }
5798
5799
5800 sub script_upload_to_ucsc
5801 {
5802     # Martin A. Hansen, August 2007.
5803
5804     # Calculate the mean of values of given keys.
5805
5806     my ( $in,        # handle to in stream
5807          $out,       # handle to out stream
5808          $options,   # options hash
5809        ) = @_;
5810
5811     # Returns nothing.
5812
5813     my ( $record, $file, $wib_file, $wig_file, $wib_dir, $fh_in, $fh_out, $i, $first, $format, $args, $type, $columns, $append, %fh_hash,
5814          $chr, $beg, $end, $block, $line, $max, $beg_block, $entry, $q_id, $clones );
5815
5816     $options->{ "short_label" } ||= $options->{ 'table' };
5817     $options->{ "long_label" }  ||= $options->{ 'table' };
5818     $options->{ "group" }       ||= $ENV{ "LOGNAME" };
5819     $options->{ "priority" }    ||= 1;
5820     $options->{ "visibility" }  ||= "pack";
5821     $options->{ "color" }       ||= join( ",", int( rand( 255 ) ), int( rand( 255 ) ), int( rand( 255 ) ) );
5822     $options->{ "chunk_size" }  ||= 10_000_000_000;    # Due to 32-bit UCSC compilation really large tables cannot be loaded in one go.
5823
5824     $file = "$BP_TMP/ucsc_upload.tmp";
5825
5826     $append = 0;
5827
5828     $first = 1;
5829
5830     $i = 0;
5831
5832     if ( $options->{ 'wiggle' } )
5833     {
5834         $options->{ "visibility" } = "full";
5835
5836         while ( $record = get_record( $in ) )
5837         {
5838             put_record( $record, $out ) if not $options->{ "no_stream" };
5839
5840             $record->{ "CHR" }     = $record->{ "S_ID" }  if not defined $record->{ "CHR" };
5841             $record->{ "CHR_BEG" } = $record->{ "S_BEG" } if not defined $record->{ "CHR_BEG" };
5842             $record->{ "CHR_END" } = $record->{ "S_END" } if not defined $record->{ "CHR_END" };
5843
5844             $fh_hash{ $record->{ "CHR" } } = Maasha::Common::write_open( "$BP_TMP/$record->{ 'CHR' }" ) if not exists $fh_hash{ $record->{ "CHR" } };
5845
5846             $fh_out = $fh_hash{ $record->{ "CHR" } };
5847
5848             Maasha::UCSC::bed_put_entry( $record, $fh_out, 5 );
5849         }
5850
5851         map { close $_ } keys %fh_hash;
5852
5853         $fh_out = Maasha::Common::write_open( $file );
5854
5855         foreach $chr ( sort keys %fh_hash )
5856         {
5857             Maasha::Common::run( "bedSort", "$BP_TMP/$chr $BP_TMP/$chr" );
5858
5859             $fh_in = Maasha::Common::read_open( "$BP_TMP/$chr" );
5860
5861             undef $block;
5862
5863             while ( $entry = Maasha::UCSC::bed_get_entry( $fh_in, 5 ) )
5864             {
5865                 $chr  = $entry->{ 'CHR' };
5866                 $beg  = $entry->{ 'CHR_BEG' };
5867                 $end  = $entry->{ 'CHR_END' };
5868                 $q_id = $entry->{ 'Q_ID' };
5869
5870                 if ( $options->{ "score" } ) {
5871                     $clones = $entry->{ 'SCORE' };
5872                 } if ( $q_id =~ /_(\d+)$/ ) {
5873                     $clones = $1;
5874                 } else {
5875                     $clones = 1;
5876                 }
5877
5878                 if ( $block )
5879                 {
5880                     if ( $beg > $max )
5881                     {
5882                         Maasha::UCSC::fixedstep_put_entry( $chr, $beg_block, $block, $fh_out );
5883                         undef $block;
5884                     }
5885                     else
5886                     {
5887                         for ( $i = $beg - $beg_block; $i < ( $beg - $beg_block ) + ( $end - $beg ); $i++ ) {
5888                             $block->[ $i ] += $clones;
5889                         }
5890
5891                         $max = Maasha::Calc::max( $max, $end );
5892                     }
5893                 }
5894
5895                 if ( not $block )
5896                 {
5897                     $beg_block = $beg;
5898                     $max       = $end;
5899
5900                     for ( $i = 0; $i < ( $end - $beg ); $i++ ) {
5901                         $block->[ $i ] += $clones;
5902                     }
5903                 }
5904             }
5905
5906             close $fh_in;
5907
5908             Maasha::UCSC::fixedstep_put_entry( $chr, $beg_block, $block, $fh_out, $options->{ "log10" } );
5909
5910             unlink "$BP_TMP/$chr";
5911         }
5912
5913         close $fh_out;
5914
5915         $wig_file = "$options->{ 'table' }.wig";
5916         $wib_file = "$options->{ 'table' }.wib";
5917
5918         $wib_dir  = "$ENV{ 'HOME' }/ucsc/wib";
5919
5920         Maasha::Common::dir_create_if_not_exists( $wib_dir );
5921
5922         # Maasha::Common::run( "wigEncode", "$file $wig_file $wib_file > /dev/null 2>&1" );
5923
5924         `cd $BP_TMP && wigEncode $file $wig_file $wib_file > /dev/null 2>&1`;
5925         Maasha::Common::run( "mv", "$BP_TMP/$wib_file $wib_dir" );
5926
5927         unlink $file;
5928
5929         $file = $wig_file;
5930
5931         $format = "WIGGLE";
5932     }
5933     else
5934     {
5935         $fh_out = Maasha::Common::write_open( $file );
5936
5937         while ( $record = get_record( $in ) )
5938         {
5939             put_record( $record, $out ) if not $options->{ "no_stream" };
5940
5941             if ( $record->{ "REC_TYPE" } eq "PSL" )
5942             {
5943                 Maasha::UCSC::psl_put_header( $fh_out ) if $first;
5944                 Maasha::UCSC::psl_put_entry( $record, $fh_out );
5945
5946                 $first = 0;
5947
5948                 $format = "PSL" if not $format;
5949             }
5950             elsif ( $record->{ "REC_TYPE" } eq "BED" and $record->{ "SEC_STRUCT" } )
5951             {
5952                 # chrom chromStart  chromEnd    name    score   strand  size    secStr  conf
5953
5954                 print $fh_out join ( "\t",
5955                     $record->{ "CHR" },
5956                     $record->{ "CHR_BEG" },
5957                     $record->{ "CHR_END" } + 1,
5958                     $record->{ "Q_ID" },
5959                     $record->{ "SCORE" },
5960                     $record->{ "STRAND" },
5961                     $record->{ "SIZE" },
5962                     $record->{ "SEC_STRUCT" },
5963                     $record->{ "CONF" },
5964                 ), "\n";
5965
5966                 $format  = "BED_SS" if not $format;
5967             }
5968             elsif ( $record->{ "REC_TYPE" } eq "BED" )
5969             {
5970                 Maasha::UCSC::bed_put_entry( $record, $fh_out, $record->{ "BED_COLS" } );
5971
5972                 $format  = "BED"                   if not $format;
5973                 $columns = $record->{ "BED_COLS" } if not $columns;
5974             }
5975             elsif ( $record->{ "REC_TYPE" } eq "PATSCAN" and $record->{ "CHR" } )
5976             {
5977                 Maasha::UCSC::bed_put_entry( $record, $fh_out, 6 );
5978
5979                 $format  = "BED" if not $format;
5980                 $columns = 6     if not $columns;
5981             }
5982             elsif ( $record->{ "REC_TYPE" } eq "BLAST" and $record->{ "S_ID" } =~ /^chr/ )
5983             {
5984                 $record->{ "CHR" }     = $record->{ "S_ID" };
5985                 $record->{ "CHR_BEG" } = $record->{ "S_BEG" };
5986                 $record->{ "CHR_END" } = $record->{ "S_END" };
5987                 $record->{ "SCORE" }   = $record->{ "BIT_SCORE" } * 1000;
5988
5989                 $format  = "BED" if not $format;
5990                 $columns = 6     if not $columns;
5991
5992                 Maasha::UCSC::bed_put_entry( $record, $fh_out );
5993             }
5994             elsif ( $record->{ "REC_TYPE" } eq "VMATCH" and $record->{ "S_ID" } =~ /^chr/i )
5995             {
5996                 $record->{ "CHR" }     = $record->{ "S_ID" };
5997                 $record->{ "CHR_BEG" } = $record->{ "S_BEG" };
5998                 $record->{ "CHR_END" } = $record->{ "S_END" };
5999                 $record->{ "SCORE" }   = $record->{ "SCORE" } || 999;
6000                 $record->{ "SCORE" }   = int( $record->{ "SCORE" } );
6001
6002                 $format  = "BED" if not $format;
6003                 $columns = 6     if not $columns;
6004
6005                 Maasha::UCSC::bed_put_entry( $record, $fh_out, 6 );
6006             }
6007
6008             if ( $i == $options->{ "chunk_size" } )
6009             {
6010                 close $fh_out;
6011
6012                 if ( $format eq "BED" ) {
6013                     Maasha::UCSC::bed_upload_to_ucsc( $BP_TMP, $file, $options, $append );
6014                 } elsif ( $format eq "PSL" ) {
6015                     Maasha::UCSC::psl_upload_to_ucsc( $file, $options, $append );
6016                 }
6017
6018                 unlink $file;
6019
6020                 $first = 1;
6021
6022                 $append = 1;
6023
6024                 $fh_out = Maasha::Common::write_open( $file );
6025             }
6026
6027             $i++;
6028         }
6029     }
6030
6031     close $fh_out;
6032
6033     if ( exists $options->{ "database" } and $options->{ "table" } )
6034     {
6035         if ( $format eq "BED" )
6036         {
6037             $type = "bed $columns";
6038
6039             Maasha::UCSC::bed_upload_to_ucsc( $BP_TMP, $file, $options, $append );
6040         }
6041         elsif ( $format eq "BED_SS" )
6042         {
6043             $options->{ "sec_struct" } = 1;
6044
6045             $type = "sec_struct";
6046
6047             Maasha::UCSC::bed_upload_to_ucsc( $BP_TMP, $file, $options, $append );
6048         }
6049         elsif ( $format eq "PSL" )
6050         {
6051             $type = "psl";
6052
6053             Maasha::UCSC::psl_upload_to_ucsc( $file, $options, $append );
6054         }
6055         elsif ( $format eq "WIGGLE" )
6056         {
6057             $type = "wig 0";
6058
6059             Maasha::UCSC::wiggle_upload_to_ucsc( $BP_TMP, $wib_dir, $file, $options );
6060         }
6061
6062         unlink $file;
6063
6064         Maasha::UCSC::update_my_tracks( $options, $type );
6065     }
6066 }
6067
6068
6069 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
6070
6071
6072 sub record2fasta
6073 {
6074     # Martin A. Hansen, July 2008.
6075
6076     # Given a biopiece record converts it to a FASTA record.
6077     # If no generic SEQ or SEQ_NAME is found, the Q_* and S_* are
6078     # tried in that order.
6079
6080     my ( $record,    # record
6081        ) = @_;
6082
6083     # Returns a tuple.
6084
6085     my ( $seq_name, $seq );
6086
6087     $seq_name = $record->{ "SEQ_NAME" } || $record->{ "Q_ID" }  || $record->{ "S_ID" };
6088     $seq      = $record->{ "SEQ" }      || $record->{ "Q_SEQ" } || $record->{ "S_SEQ" };
6089
6090     if ( defined $seq_name and defined $seq ) {
6091         return wantarray ? ( $seq_name, $seq ) : [ $seq_name, $seq ];
6092     } else {
6093         return;
6094     }
6095 }
6096
6097
6098 sub read_stream
6099 {
6100     # Martin A. Hansen, July 2007.
6101
6102     # Opens a stream to STDIN or a file,
6103
6104     my ( $path,   # path - OPTIONAL
6105        ) = @_;
6106
6107     # Returns filehandle.
6108
6109     my ( $fh );
6110
6111     if ( not -t STDIN ) {
6112         $fh = Maasha::Common::read_stdin();
6113     } elsif ( not $path ) {
6114 #        Maasha::Common::error( qq(no data stream) );
6115     } else {
6116         $fh = Maasha::Common::read_open( $path );
6117     }
6118
6119 #    $fh->autoflush(1) if $fh;  # Disable file buffer for debugging.
6120
6121     return $fh;
6122 }
6123
6124
6125 sub write_stream
6126 {
6127     # Martin A. Hansen, August 2007.
6128
6129     # Opens a stream to STDOUT or a file.
6130
6131     my ( $path,   # path          - OPTIONAL
6132          $gzip,   # compress data - OPTIONAL
6133        ) = @_;
6134
6135     # Returns filehandle.
6136
6137     my ( $fh );
6138
6139     if ( $path ) {
6140         $fh = Maasha::Common::write_open( $path, $gzip );
6141     } else {
6142         $fh = Maasha::Common::write_stdout();
6143     }
6144
6145     return $fh;
6146 }
6147
6148
6149 sub get_record
6150 {
6151     # Martin A. Hansen, July 2007.
6152
6153     # Reads one record at a time and converts that record
6154     # to a Perl data structure (a hash) which is returned.
6155
6156     my ( $fh,   # handle to stream
6157        ) = @_;
6158
6159     # Returns a hash.
6160
6161     my ( $block, @lines, $line, $key, $value, %record );
6162
6163     local $/ = "\n---\n";
6164
6165     $block = <$fh>;
6166
6167     chomp $block;
6168
6169     return if not defined $block;
6170
6171     @lines = split "\n", $block;
6172
6173     foreach $line ( @lines )
6174     {
6175         ( $key, $value ) = split ": ", $line, 2;
6176
6177         $record{ $key } = $value;
6178     }
6179
6180     return wantarray ? %record : \%record;
6181 }
6182
6183
6184 sub put_record
6185 {
6186     # Martin A. Hansen, July 2007.
6187
6188     # Given a Perl datastructure (a hash ref) emits this to STDOUT or a filehandle.
6189
6190     my ( $data,   # data structure
6191          $fh,     # file handle - OPTIONAL
6192        ) = @_;
6193
6194     # Returns nothing.
6195
6196     if ( scalar keys %{ $data } )
6197     {
6198         if ( $fh )
6199         {
6200             map { print $fh "$_: $data->{ $_ }\n" } keys %{ $data };
6201             print $fh "---\n";
6202         }
6203         else
6204         {
6205             map { print "$_: $data->{ $_ }\n" } keys %{ $data };
6206             print "---\n";
6207         }
6208     }
6209
6210     undef $data;
6211 }
6212
6213
6214 sub getopt_files
6215 {
6216     # Martin A. Hansen, November 2007.
6217
6218     # Extracts files from an explicit GetOpt::Long argument
6219     # allowing for the use of glob. E.g.
6220     # --data_in=test.fna
6221     # --data_in=test.fna,test2.fna
6222     # --data_in=*.fna
6223     # --data_in=test.fna,/dir/*.fna
6224
6225     my ( $option,   # option from GetOpt::Long
6226        ) = @_;
6227
6228     # Returns a list.
6229
6230     my ( $elem, @files );
6231
6232     foreach $elem ( split ",", $option )
6233     {
6234         if ( -f $elem ) {
6235             push @files, $elem;
6236         } elsif ( $elem =~ /\*/ ) {
6237             push @files, glob( $elem );
6238         }
6239     }
6240
6241     return wantarray ? @files : \@files;
6242 }
6243
6244
6245 sub sig_handler
6246 {
6247     # Martin A. Hansen, April 2008.
6248
6249     # Removes temporary directory and exits gracefully.
6250     # This subroutine is meant to be run always as the last
6251     # thing even if a script is dies or is interrupted
6252     # or killed.
6253
6254     my ( $sig,   # signal from the %SIG
6255        ) = @_;
6256
6257     # print STDERR "signal->$sig<-\n";
6258
6259     chomp $sig;
6260
6261     sleep 1;
6262
6263     if ( -d $BP_TMP )
6264     {
6265         if ( $sig =~ /MAASHA_ERROR/ ) {
6266             print STDERR "\nProgram '$script' had an error"                     . "  -  Please wait for temporary data to be removed\n";
6267         } elsif ( $sig eq "INT" ) {
6268             print STDERR "\nProgram '$script' interrupted (ctrl-c was pressed)" . "  -  Please wait for temporary data to be removed\n";
6269         } elsif ( $sig eq "TERM" ) {
6270             print STDERR "\nProgram '$script' terminated (someone used kill?)"  . "  -  Please wait for temporary data to be removed\n";
6271         } else {
6272             print STDERR "\nProgram '$script' died->$sig"                       . "  -  Please wait for temporary data to be removed\n";
6273         }
6274
6275         clean_tmp();
6276     }
6277
6278     exit( 0 );
6279 }
6280
6281
6282 sub clean_tmp
6283 {
6284     # Martin A. Hansen, July 2008.
6285
6286     # Cleans out any unused temporary files and direcotries in BP_TMP.
6287
6288     # Returns nothing.
6289
6290     my ( $tmpdir, @dirs, $curr_pid, $dir, $user, $sid, $pid );
6291
6292     $tmpdir = $ENV{ 'BP_TMP' } || Maasha::Common::error( 'No BP_TMP variable in environment.' );
6293
6294     $curr_pid = Maasha::Common::get_processid();
6295
6296     @dirs = Maasha::Common::ls_dirs( $tmpdir );
6297
6298     foreach $dir ( @dirs )
6299     {
6300         if ( $dir =~ /^$tmpdir\/(.+)_(\d+)_(\d+)_bp_tmp$/ )
6301         {
6302             $user = $1;
6303             $sid  = $2;
6304             $pid  = $3;
6305
6306             if ( $user eq Maasha::Common::get_user() )
6307             {
6308                 if ( not Maasha::Common::process_running( $pid ) )
6309                 {
6310                     # print STDERR "Removing stale dir: $dir\n";
6311                     Maasha::Common::dir_remove( $dir );
6312                 }
6313                 elsif ( $pid == $curr_pid )
6314                 {
6315                     # print STDERR "Removing current dir: $dir\n";
6316                     Maasha::Common::dir_remove( $dir );
6317                 }
6318             }
6319         }
6320     }
6321 }
6322
6323
6324 END
6325 {
6326     clean_tmp();
6327 }
6328
6329
6330 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
6331
6332 1;
6333
6334 __END__