code_perl/Maasha/Biopieces.pm

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