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