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