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