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, $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
2679     while ( $record = get_record( $in ) )
2680     {
2681         if ( $bed_entry = Maasha::UCSC::BED::biopiece2bed( $record ) ) {
2682             Maasha::UCSC::BED::bed_entry_put( $bed_entry, $fh_out, undef, $options->{ 'check' } );
2683         }
2684     }
2685
2686     close $fh_out;
2687
2688     $file_hash = Maasha::UCSC::BED::bed_file_split_on_chr( $bed_file, $BP_TMP );
2689
2690     unlink $bed_file;
2691
2692     foreach $chr ( sort keys %{ $file_hash } )
2693     {
2694         $bed_file       = $file_hash->{ $chr };
2695
2696         $fixedstep_file = Maasha::UCSC::Wiggle::fixedstep_calc( $bed_file, $chr, $options->{ 'use_score' }, $options->{ 'use_log10' } );
2697
2698         #$fixedstep_file = "$bed_file.fixedstep";
2699
2700         # Maasha::Common::run( "bed2fixedstep", "< $bed_file > $fixedstep_file" );
2701
2702
2703         $fh_in = Maasha::Filesys::file_read_open( $fixedstep_file );
2704
2705         while ( $fixedstep_entry = Maasha::UCSC::Wiggle::fixedstep_entry_get( $fh_in ) )
2706         {
2707             if ( $record = Maasha::UCSC::Wiggle::fixedstep2biopiece( $fixedstep_entry ) ) {
2708                 put_record( $record, $out );
2709             }
2710         }
2711
2712         close $fh_in;
2713
2714         unlink $bed_file;
2715         unlink $fixedstep_file;
2716     }
2717 }
2718
2719
2720 sub script_reverse_seq
2721 {
2722     # Martin A. Hansen, August 2007.
2723
2724     # Reverse sequence in record.
2725
2726     my ( $in,    # handle to in stream
2727          $out,   # handle to out stream
2728        ) = @_;
2729
2730     # Returns nothing.
2731
2732     my ( $record );
2733
2734     while ( $record = get_record( $in ) )
2735     {
2736         if ( $record->{ "SEQ" } ) {
2737             $record->{ "SEQ" } = reverse $record->{ "SEQ" };
2738         }
2739
2740         put_record( $record, $out );
2741     }
2742 }
2743
2744
2745 sub script_complement_seq
2746 {
2747     # Martin A. Hansen, August 2007.
2748
2749     # Complement sequence in record.
2750
2751     my ( $in,     # handle to in stream
2752          $out,    # handle to out stream
2753        ) = @_;
2754
2755     # Returns nothing.
2756
2757     my ( $record, $type );
2758
2759     while ( $record = get_record( $in ) )
2760     {
2761         if ( $record->{ "SEQ" } )
2762         {
2763             if ( not $type ) {
2764                 $type = Maasha::Seq::seq_guess_type( $record->{ "SEQ" } );
2765             }
2766
2767             if ( $type eq "rna" ) {
2768                 Maasha::Seq::rna_comp( \$record->{ "SEQ" } );
2769             } elsif ( $type eq "dna" ) {
2770                 Maasha::Seq::dna_comp( \$record->{ "SEQ" } );
2771             }
2772         }
2773
2774         put_record( $record, $out );
2775     }
2776 }
2777
2778
2779 sub script_remove_indels
2780 {
2781     # Martin A. Hansen, August 2007.
2782
2783     # Remove indels from sequences in stream.
2784
2785     my ( $in,     # handle to in stream
2786          $out,    # handle to out stream
2787        ) = @_;
2788
2789     # Returns nothing.
2790
2791     my ( $record );
2792
2793     while ( $record = get_record( $in ) )
2794     {
2795         $record->{ 'SEQ' } =~ tr/-~.//d if $record->{ "SEQ" };
2796
2797         put_record( $record, $out );
2798     }
2799 }
2800
2801
2802 sub script_transliterate_seq
2803 {
2804     # Martin A. Hansen, August 2007.
2805
2806     # Transliterate chars from sequence in record.
2807
2808     my ( $in,        # handle to in stream
2809          $out,       # handle to out stream
2810          $options,   # options hash
2811        ) = @_;
2812
2813     # Returns nothing.
2814
2815     my ( $record, $search, $replace, $delete );
2816
2817     $search  = $options->{ "search" }  || "";
2818     $replace = $options->{ "replace" } || "";
2819     $delete  = $options->{ "delete" }  || "";
2820
2821     while ( $record = get_record( $in ) )
2822     {
2823         if ( $record->{ "SEQ" } )
2824         {
2825             if ( $search and $replace ) {
2826                 eval "\$record->{ 'SEQ' } =~ tr/$search/$replace/";
2827             } elsif ( $delete ) {
2828                 eval "\$record->{ 'SEQ' } =~ tr/$delete//d";
2829             }
2830         }
2831
2832         put_record( $record, $out );
2833     }
2834 }
2835
2836
2837 sub script_transliterate_vals
2838 {
2839     # Martin A. Hansen, April 2008.
2840
2841     # Transliterate chars from values in record.
2842
2843     my ( $in,        # handle to in stream
2844          $out,       # handle to out stream
2845          $options,   # options hash
2846        ) = @_;
2847
2848     # Returns nothing.
2849
2850     my ( $record, $search, $replace, $delete, $key );
2851
2852     $search  = $options->{ "search" }  || "";
2853     $replace = $options->{ "replace" } || "";
2854     $delete  = $options->{ "delete" }  || "";
2855
2856     while ( $record = get_record( $in ) )
2857     {
2858         foreach $key ( @{ $options->{ "keys" } } )
2859         {
2860             if ( exists $record->{ $key } )
2861             {
2862                 if ( $search and $replace ) {
2863                     eval "\$record->{ $key } =~ tr/$search/$replace/";
2864                 } elsif ( $delete ) {
2865                     eval "\$record->{ $key } =~ tr/$delete//d";
2866                 }
2867             }
2868         }
2869
2870         put_record( $record, $out );
2871     }
2872 }
2873
2874
2875 sub script_translate_seq
2876 {
2877     # Martin A. Hansen, February 2008.
2878
2879     # Translate DNA sequence into protein sequence.
2880
2881     my ( $in,        # handle to in stream
2882          $out,       # handle to out stream
2883          $options,   # options hash
2884        ) = @_;
2885
2886     # Returns nothing.
2887
2888     my ( $record, $frame, %new_record );
2889
2890     $options->{ "frames" } ||= [ 1, 2, 3, -1, -2, -3 ];
2891
2892     while ( $record = get_record( $in ) )
2893     {
2894         if ( $record->{ "SEQ" } )
2895         {
2896             if ( Maasha::Seq::seq_guess_type( $record->{ "SEQ" } ) eq "dna" )
2897             {
2898                 foreach $frame ( @{ $options->{ "frames" } } )
2899                 {
2900                     %new_record = %{ $record };
2901
2902                     $new_record{ "SEQ" }     = Maasha::Seq::translate( $record->{ "SEQ" }, $frame );
2903                     $new_record{ "SEQ_LEN" } = length $new_record{ "SEQ" };
2904                     $new_record{ "FRAME" }   = $frame;
2905
2906                     put_record( \%new_record, $out );
2907                 }
2908             }
2909         }
2910         else
2911         {
2912             put_record( $record, $out );
2913         }
2914     }
2915 }
2916
2917
2918 sub script_extract_seq
2919 {
2920     # Martin A. Hansen, August 2007.
2921
2922     # Extract subsequences from sequences in record.
2923
2924     my ( $in,        # handle to in stream
2925          $out,       # handle to out stream
2926          $options,   # options hash
2927        ) = @_;
2928
2929     # Returns nothing.
2930
2931     my ( $beg, $end, $len, $record );
2932
2933     if ( not defined $options->{ "beg" } or $options->{ "beg" } < 0 ) {
2934         $beg = 0;
2935     } else {
2936         $beg = $options->{ "beg" } - 1;   # correcting for start offset
2937     }
2938
2939     if ( defined $options->{ "end" } and $options->{ "end" } - 1 < $beg ) {
2940         $end = $beg - 1;
2941     } elsif ( defined $options->{ "end" } ) {
2942         $end = $options->{ "end" } - 1;   # correcting for start offset
2943     }
2944
2945     $len = $options->{ "len" };
2946
2947 #    print "beg->$beg,  end->$end,  len->$len\n";
2948
2949     while ( $record = get_record( $in ) )
2950     {
2951         if ( $record->{ "SEQ" } )
2952         {
2953             if ( defined $beg and defined $end )
2954             {
2955                 if ( $end - $beg + 1 > length $record->{ "SEQ" } ) {
2956                     $record->{ "SEQ" } = substr $record->{ "SEQ" }, $beg;
2957                 } else {
2958                     $record->{ "SEQ" } = substr $record->{ "SEQ" }, $beg, $end - $beg + 1;
2959                 }
2960             }
2961             elsif ( defined $beg and defined $len )
2962             {
2963                 if ( $len > length $record->{ "SEQ" } ) {
2964                     $record->{ "SEQ" } = substr $record->{ "SEQ" }, $beg;
2965                 } else {
2966                     $record->{ "SEQ" } = substr $record->{ "SEQ" }, $beg, $len;
2967                 }
2968             }
2969             elsif ( defined $beg )
2970             {
2971                 $record->{ "SEQ" } = substr $record->{ "SEQ" }, $beg;
2972             }
2973         }
2974
2975         $record->{ "SEQ_LEN" } = length $record->{ "SEQ" };
2976
2977         put_record( $record, $out );
2978     }
2979 }
2980
2981
2982 sub script_get_genome_seq
2983 {
2984     # Martin A. Hansen, December 2007.
2985
2986     # Gets a subsequence from a genome.
2987
2988     my ( $in,        # handle to in stream
2989          $out,       # handle to out stream
2990          $options,   # options hash
2991        ) = @_;
2992
2993     # Returns nothing.
2994
2995     my ( $record, $genome, $genome_file, $index_file, $index, $fh, $index_head, $index_beg, $index_len, $beg, $len, %lookup_hash, @begs, @lens, $i, $seq );
2996
2997     $options->{ "flank" } ||= 0;
2998
2999     if ( $options->{ "genome" } )
3000     {
3001         $genome      = $options->{ "genome" };
3002
3003         $genome_file = "$ENV{ 'BP_DATA' }/genomes/$genome/fasta/$genome.fna";
3004         $index_file  = "$ENV{ 'BP_DATA' }/genomes/$genome/fasta/$genome.index";
3005
3006         $fh          = Maasha::Common::read_open( $genome_file );
3007         $index       = Maasha::Fasta::index_retrieve( $index_file );
3008
3009         shift @{ $index }; # Get rid of the file size info
3010
3011         map { $lookup_hash{ $_->[ 0 ] } = [ $_->[ 1 ], $_->[ 2 ] ] } @{ $index };
3012
3013         if ( exists $lookup_hash{ $options->{ "chr" } } and defined $options->{ "beg" } and ( defined $options->{ "end" } or defined $options->{ "len" } ) )
3014         {
3015             ( $index_beg, $index_len ) = @{ $lookup_hash{ $options->{ "chr" } } };
3016
3017             $beg = $index_beg + $options->{ "beg" } - 1;
3018
3019             if ( $options->{ "len" } ) {
3020                 $len = $options->{ "len" };
3021             } elsif ( $options->{ "end" } ) {
3022                 $len = ( $options->{ "end" } - $options->{ "beg" } + 1 );
3023             }
3024
3025             $beg -= $options->{ "flank" };
3026             $len += 2 * $options->{ "flank" };
3027
3028             if ( $beg <= $index_beg )
3029             {
3030                 $len -= $index_beg - $beg;
3031                 $beg = $index_beg;
3032             }
3033
3034             $len = $index_beg + $index_len - $beg if $beg + $len > $index_beg + $index_len;
3035
3036             next if $beg > $index_beg + $index_len;
3037
3038             $record->{ "CHR" }     = $options->{ "chr" };
3039             $record->{ "CHR_BEG" } = $beg - $index_beg;
3040             $record->{ "CHR_END" } = $record->{ "CHR_BEG" } + $len - 1;
3041
3042             $record->{ "SEQ" }     = Maasha::Common::file_read( $fh, $beg, $len );
3043             $record->{ "SEQ_LEN" } = $len;
3044
3045             put_record( $record, $out );
3046         }
3047     }
3048
3049     while ( $record = get_record( $in ) )
3050     {
3051         if ( $options->{ "genome" } and not $record->{ "SEQ" } )
3052         {
3053             if ( $record->{ "REC_TYPE" } eq "BED" and exists $lookup_hash{ $record->{ "CHR" } } )
3054             {
3055                 ( $index_beg, $index_len ) = @{ $lookup_hash{ $record->{ "CHR" } } };
3056
3057                 $beg = $record->{ "CHR_BEG" } + $index_beg;
3058                 $len = $record->{ "CHR_END" } - $record->{ "CHR_BEG" } + 1;
3059             }
3060             elsif ( $record->{ "REC_TYPE" } eq "PSL" and exists $lookup_hash{ $record->{ "S_ID" } } )
3061             {
3062                 ( $index_beg, $index_len ) = @{ $lookup_hash{ $record->{ "S_ID" } } };
3063
3064                 $beg = $record->{ "S_BEG" } + $index_beg;
3065                 $len = $record->{ "S_END" } - $record->{ "S_BEG" } + 1;
3066             }
3067             elsif ( $record->{ "REC_TYPE" } eq "BLAST" and exists $lookup_hash{ $record->{ "S_ID" } } )
3068             {
3069                 ( $index_beg, $index_len ) = @{ $lookup_hash{ $record->{ "S_ID" } } };
3070
3071                 $beg = $record->{ "S_BEG" } + $index_beg;
3072                 $len = $record->{ "S_END" } - $record->{ "S_BEG" } + 1;
3073             }
3074
3075             $beg -= $options->{ "flank" };
3076             $len += 2 * $options->{ "flank" };
3077
3078             if ( $beg <= $index_beg )
3079             {
3080                 $len -= $index_beg - $beg;
3081                 $beg = $index_beg;
3082             }
3083
3084             $len = $index_beg + $index_len - $beg if $beg + $len > $index_beg + $index_len;
3085
3086             next if $beg > $index_beg + $index_len;
3087
3088             $record->{ "CHR_BEG" } = $beg - $index_beg;
3089             $record->{ "CHR_END" } = $record->{ "CHR_BEG" } + $len - 1;
3090
3091             $record->{ "SEQ" } = Maasha::Common::file_read( $fh, $beg, $len );
3092
3093             if ( $record->{ "STRAND" } and $record->{ "STRAND" } eq "-" )
3094             {
3095                 Maasha::Seq::dna_comp( \$record->{ "SEQ" } );
3096                 $record->{ "SEQ" } = reverse $record->{ "SEQ" };
3097             }
3098
3099             if ( $options->{ "mask" } )
3100             {
3101                 if ( $record->{ "BLOCK_COUNT" } > 1 ) # uppercase hit block segments and lowercase the rest.
3102                 {
3103                     $record->{ "SEQ" } = lc $record->{ "SEQ" };
3104
3105                     @begs = split ",", $record->{ "Q_BEGS" };
3106                     @lens = split ",", $record->{ "BLOCK_LENS" };
3107
3108                     for ( $i = 0; $i < @begs; $i++ ) {
3109                         substr $record->{ "SEQ" }, $begs[ $i ], $lens[ $i ], uc substr $record->{ "SEQ" }, $begs[ $i ], $lens[ $i ];
3110                     }
3111                 }
3112             }
3113             elsif ( $options->{ "splice" } )
3114             {
3115                 if ( $record->{ "BLOCK_COUNT" } > 1 ) # splice block sequences
3116                 {
3117                     $seq  = "";
3118                     @begs = split ",", $record->{ "Q_BEGS" };
3119                     @lens = split ",", $record->{ "BLOCK_LENS" };
3120
3121                     for ( $i = 0; $i < @begs; $i++ ) {
3122                         $seq .= substr $record->{ "SEQ" }, $begs[ $i ], $lens[ $i ];
3123                     }
3124
3125                     $record->{ "SEQ" } = $seq;
3126                 }
3127             }
3128
3129             $record->{ "SEQ_LEN" } = length $record->{ "SEQ" };
3130         }
3131
3132         put_record( $record, $out );
3133     }
3134
3135     close $fh if $fh;
3136 }
3137
3138
3139 sub script_get_genome_align
3140 {
3141     # Martin A. Hansen, April 2008.
3142
3143     # Gets a subalignment from a multiple genome alignment.
3144
3145     my ( $in,        # handle to in stream
3146          $out,       # handle to out stream
3147          $options,   # options hash
3148        ) = @_;
3149
3150     # Returns nothing.
3151
3152     my ( $record, $maf_track, $align, $align_num, $beg, $end, $len, $entry );
3153
3154     $options->{ "strand" } ||= "+";
3155
3156     $align_num = 1;
3157
3158     $maf_track = Maasha::Config::maf_track( $options->{ "genome" } );
3159
3160     if ( $options->{ "chr" } and $options->{ "beg" } and ( $options->{ "end" } or $options->{ "len" } ) )
3161     {
3162         $beg = $options->{ "beg" } - 1;
3163
3164         if ( $options->{ "end" } ) {
3165             $end = $options->{ "end" };
3166         } elsif ( $options->{ "len" } ) {
3167             $end = $beg + $options->{ "len" };
3168         }
3169
3170         $align = Maasha::UCSC::maf_extract( $BP_TMP, $options->{ "genome" }, $maf_track, $options->{ "chr" }, $beg, $end, $options->{ "strand" } );
3171
3172         foreach $entry ( @{ $align } )
3173         {
3174             $entry->{ "CHR" }     = $record->{ "CHR" };
3175             $entry->{ "CHR_BEG" } = $record->{ "CHR_BEG" };
3176             $entry->{ "CHR_END" } = $record->{ "CHR_END" };
3177             $entry->{ "STRAND" }  = $record->{ "STRAND" } || '+';
3178             $entry->{ "Q_ID" }    = $record->{ "Q_ID" };
3179             $entry->{ "SCORE" }   = $record->{ "SCORE" };
3180
3181             put_record( $entry, $out );
3182         }
3183     }
3184
3185     while ( $record = get_record( $in ) )
3186     {
3187         if ( $record->{ "REC_TYPE" } eq "BED" )
3188         {
3189             $align = Maasha::UCSC::maf_extract( $BP_TMP, $options->{ "genome" }, $maf_track, $record->{ "CHR" }, $record->{ "CHR_BEG" }, $record->{ "CHR_END" }, $record->{ "STRAND" } );
3190         }
3191         elsif ( $record->{ "REC_TYPE" } eq "VMATCH" )
3192         {
3193             $align = Maasha::UCSC::maf_extract( $BP_TMP, $options->{ "genome" }, $maf_track, $record->{ "S_ID" }, $record->{ "S_BEG" }, $record->{ "S_END" } + 1, $record->{ "STRAND" } );
3194         }
3195         elsif ( $record->{ "REC_TYPE" } eq "PSL" )
3196         {
3197             $align = Maasha::UCSC::maf_extract( $BP_TMP, $options->{ "genome" }, $maf_track, $record->{ "S_ID" }, $record->{ "S_BEG" }, $record->{ "S_END" }, $record->{ "STRAND" } );
3198         }
3199         elsif ( $record->{ "REC_TYPE" } eq "BLAST" )
3200         {
3201             $align = Maasha::UCSC::maf_extract( $BP_TMP, $options->{ "genome" }, $maf_track, $record->{ "S_ID" }, $record->{ "S_BEG" }, $record->{ "S_END" }, $record->{ "STRAND" } );
3202         }
3203
3204         foreach $entry ( @{ $align } )
3205         {
3206             $entry->{ "CHR" }     = $record->{ "CHR" };
3207             $entry->{ "CHR_BEG" } = $record->{ "CHR_BEG" };
3208             $entry->{ "CHR_END" } = $record->{ "CHR_END" };
3209             $entry->{ "STRAND" }  = $record->{ "STRAND" };
3210             $entry->{ "Q_ID" }    = $record->{ "Q_ID" };
3211             $entry->{ "SCORE" }   = $record->{ "SCORE" };
3212
3213             put_record( $entry, $out );
3214         }
3215
3216         $align_num++;
3217     }
3218 }
3219
3220
3221 sub script_get_genome_phastcons
3222 {
3223     # Martin A. Hansen, February 2008.
3224
3225     # Get phastcons scores from genome intervals.
3226
3227     my ( $in,        # handle to in stream
3228          $out,       # handle to out stream
3229          $options,   # options hash
3230        ) = @_;
3231
3232     # Returns nothing.
3233
3234     my ( $phastcons_file, $phastcons_index, $index, $fh_phastcons, $scores, $record );
3235
3236     $options->{ "flank" } ||= 0;
3237
3238     $phastcons_file  = Maasha::Config::genome_phastcons( $options->{ "genome" } );
3239     $phastcons_index = Maasha::Config::genome_phastcons_index( $options->{ "genome" } );
3240
3241     $index           = Maasha::UCSC::fixedstep_index_retrieve( $phastcons_index );
3242     $fh_phastcons    = Maasha::Common::read_open( $phastcons_file );
3243
3244     if ( defined $options->{ "chr" } and defined $options->{ "beg" } and ( defined $options->{ "end" } or defined $options->{ "len" } ) )
3245     {
3246         $options->{ "beg" } -= 1;   # request is 1-based
3247         $options->{ "end" } -= 1;   # request is 1-based
3248
3249         if ( $options->{ "len" } ) {
3250             $options->{ "end" } = $options->{ "beg" } + $options->{ "len" } - 1;
3251         }
3252
3253         $scores = Maasha::UCSC::fixedstep_index_lookup( $index, $fh_phastcons, $options->{ "chr" }, $options->{ "beg" }, $options->{ "end" }, $options->{ "flank" } );
3254
3255         $record->{ "CHR" }       = $options->{ "chr" };
3256         $record->{ "CHR_BEG" }   = $options->{ "beg" } - $options->{ "flank" };
3257         $record->{ "CHR_END" }   = $options->{ "end" } + $options->{ "flank" };
3258
3259         $record->{ "PHASTCONS" }   = join ",", @{ $scores };
3260         $record->{ "PHAST_COUNT" } = scalar @{ $scores };  # DEBUG
3261
3262         put_record( $record, $out );
3263     }
3264
3265     while ( $record = get_record( $in ) )
3266     {
3267         if ( $record->{ "REC_TYPE" } eq "BED" )
3268         {
3269             $scores = Maasha::UCSC::fixedstep_index_lookup( $index, $fh_phastcons, $record->{ "CHR" }, $record->{ "CHR_BEG" }, $record->{ "CHR_END" }, $options->{ "flank" } );
3270         }
3271         elsif ( $record->{ "REC_TYPE" } eq "PSL" )
3272         {
3273             $scores = Maasha::UCSC::fixedstep_index_lookup( $index, $fh_phastcons, $record->{ "S_ID" }, $record->{ "S_BEG" }, $record->{ "S_END" }, $options->{ "flank" } );
3274         }
3275         elsif ( $record->{ "REC_TYPE" } eq "BLAST" )
3276         {
3277             $scores = Maasha::UCSC::fixedstep_index_lookup( $index, $fh_phastcons, $record->{ "S_ID" }, $record->{ "S_BEG" }, $record->{ "S_END" }, $options->{ "flank" } );
3278         }
3279
3280         $record->{ "PHASTCONS" } = join ",", @{ $scores } if @{ $scores };
3281 #        $record->{ "PHAST_COUNT" } = @{ $scores } if @{ $scores };  # DEBUG
3282
3283         put_record( $record, $out );
3284     }
3285
3286     close $fh_phastcons if $fh_phastcons;
3287 }
3288
3289
3290 sub script_fold_seq
3291 {
3292     # Martin A. Hansen, December 2007.
3293
3294     # Folds sequences in stream into secondary structures.
3295
3296     my ( $in,     # handle to in stream
3297          $out,    # handle to out stream
3298        ) = @_;
3299
3300     # Returns nothing.
3301
3302     my ( $record, $type, $struct, $index );
3303
3304     while ( $record = get_record( $in ) )
3305     {
3306         if ( $record->{ "SEQ" } )
3307         {
3308             if ( not $type ) {
3309                 $type = Maasha::Seq::seq_guess_type( $record->{ "SEQ" } );
3310             }
3311
3312             if ( $type ne "protein" )
3313             {
3314                 ( $struct, $index ) = Maasha::Seq::fold_struct_rnafold( $record->{ "SEQ" } );
3315                 $record->{ "SEC_STRUCT" }  = $struct;
3316                 $record->{ "FREE_ENERGY" } = $index;
3317                 $record->{ "SCORE" }       = abs int $index;
3318                 $record->{ "SIZE" }        = length $struct;
3319                 $record->{ "CONF" }        = "1," x $record->{ "SIZE" };
3320             }
3321         }
3322
3323         put_record( $record, $out );
3324     }
3325 }
3326
3327
3328 sub script_split_seq
3329 {
3330     # Martin A. Hansen, August 2007.
3331
3332     # Split a sequence in stream into words.
3333
3334     my ( $in,        # handle to in stream
3335          $out,       # handle to out stream
3336          $options,   # options hash
3337        ) = @_;
3338
3339     # Returns nothing.
3340
3341     my ( $record, $new_record, $i, $subseq, %lookup );
3342
3343     $options->{ "word_size" } ||= 7;
3344
3345     while ( $record = get_record( $in ) )
3346     {
3347         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } )
3348         {
3349             for ( $i = 0; $i < length( $record->{ "SEQ" } ) - $options->{ "word_size" } + 1; $i++ )
3350             {
3351                 $subseq = substr $record->{ "SEQ" }, $i, $options->{ "word_size" };
3352
3353                 if ( $options->{ "uniq" } and not $lookup{ $subseq } )
3354                 {
3355                     $new_record->{ "SEQ_NAME" } = $record->{ "SEQ_NAME" } . "[" . ( $i + 1 ) . "-" . ( $i + $options->{ "word_size" } ) . "]";
3356                     $new_record->{ "SEQ" }      = $subseq;
3357
3358                     put_record( $new_record, $out );
3359
3360                     $lookup{ $subseq } = 1;
3361                 }
3362                 else
3363                 {
3364                     $new_record->{ "SEQ_NAME" } = $record->{ "SEQ_NAME" } . "[" . ( $i + 1 ) . "-" . ( $i + $options->{ "word_size" } ) . "]";
3365                     $new_record->{ "SEQ" }      = $subseq;
3366
3367                     put_record( $new_record, $out );
3368                 }
3369             }
3370         }
3371         else
3372         {
3373             put_record( $record, $out );
3374         }
3375     }
3376 }
3377
3378
3379 sub script_split_bed
3380 {
3381     # Martin A. Hansen, June 2008.
3382
3383     # Split a BED record into overlapping windows.
3384
3385     my ( $in,        # handle to in stream
3386          $out,       # handle to out stream
3387          $options,   # options hash
3388        ) = @_;
3389
3390     # Returns nothing.
3391
3392     my ( $record, $new_record, $i );
3393
3394     $options->{ "window_size" } ||= 20;
3395     $options->{ "step_size" }   ||= 1;
3396
3397     while ( $record = get_record( $in ) )
3398     {
3399         if ( $record->{ "CHR" } and $record->{ "CHR_BEG" } and $record->{ "CHR_END" } )
3400         {
3401             $record->{ "BED_LEN" } = $record->{ "CHR_END" } - $record->{ "CHR_BEG" } + 1;
3402
3403             for ( $i = 0; $i < $record->{ "BED_LEN" } - $options->{ "window_size" }; $i += $options->{ "step_size" } )
3404             {
3405                 $new_record->{ "REC_TYPE" } = "BED";
3406                 $new_record->{ "CHR" }      = $record->{ "CHR" };
3407                 $new_record->{ "CHR_BEG" }  = $record->{ "CHR_BEG" } + $i;
3408                 $new_record->{ "CHR_END" }  = $record->{ "CHR_BEG" } + $i + $options->{ "window_size" };
3409                 $new_record->{ "BED_LEN" }  = $options->{ "window_size" };
3410                 $new_record->{ "Q_ID" }     = $record->{ "Q_ID" } . "_$i";
3411                 $new_record->{ "SCORE" }    = $record->{ "SCORE" };
3412                 $new_record->{ "STRAND" }   = $record->{ "STRAND" };
3413
3414                 put_record( $new_record, $out );
3415             }
3416         }
3417         else
3418         {
3419             put_record( $record, $out );
3420         }
3421     }
3422 }
3423
3424
3425 sub script_align_seq
3426 {
3427     # Martin A. Hansen, August 2007.
3428
3429     # Align sequences in stream.
3430
3431     my ( $in,    # handle to in stream
3432          $out,   # handle to out stream
3433        ) = @_;
3434
3435     # Returns nothing.
3436
3437     my ( $record, @entries, $entry );
3438
3439     while ( $record = get_record( $in ) )
3440     {
3441         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } ) {
3442             push @entries, [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ];
3443         } elsif ( $record->{ "Q_ID" } and $record->{ "SEQ" } ) {
3444             push @entries, [ $record->{ "Q_ID" }, $record->{ "SEQ" } ];
3445         } else {
3446             put_record( $record, $out );
3447         }
3448     }
3449
3450     @entries = Maasha::Align::align( \@entries );
3451
3452     foreach $entry ( @entries )
3453     {
3454         if ( $entry->[ SEQ_NAME ] and $entry->[ SEQ ] )
3455         {
3456             $record = {
3457                 SEQ_NAME => $entry->[ SEQ_NAME ],
3458                 SEQ      => $entry->[ SEQ ],
3459             };
3460
3461             put_record( $record, $out );
3462         }
3463     }
3464 }
3465
3466
3467 sub script_tile_seq
3468 {
3469     # Martin A. Hansen, February 2008.
3470
3471     # Using the first sequence in stream as reference, tile
3472     # all subsequent sequences based on pairwise alignments.
3473
3474     my ( $in,        # handle to in stream
3475          $out,       # handle to out stream
3476          $options,   # options hash
3477        ) = @_;
3478
3479     # Returns nothing.
3480
3481     my ( $record, $first, $ref_entry, @entries );
3482
3483     $first = 1;
3484
3485     while ( $record = get_record( $in ) )
3486     {
3487         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } )
3488         {
3489             if ( $first )
3490             {
3491                 $ref_entry = [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ];
3492
3493                 $first = 0;
3494             }
3495             else
3496             {
3497                 push @entries, [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ];
3498             }
3499         }
3500         else
3501         {
3502             put_record( $record, $out );
3503         }
3504     }
3505
3506     @entries = Maasha::Align::align_tile( $ref_entry, \@entries, $options );
3507
3508     map { put_record( { SEQ_NAME => $_->[ SEQ_NAME ], SEQ => $_->[ SEQ ] }, $out ) } @entries;
3509 }
3510
3511
3512 sub script_invert_align
3513 {
3514     # Martin A. Hansen, February 2008.
3515
3516     # Inverts an alignment showing only non-mathing residues
3517     # using the first sequence as reference.
3518
3519     my ( $in,        # handle to in stream
3520          $out,       # handle to out stream
3521          $options,   # options hash
3522        ) = @_;
3523
3524     # Returns nothing.
3525
3526     my ( $record, @entries );
3527
3528     while ( $record = get_record( $in ) )
3529     {
3530         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } )
3531         {
3532             push @entries, [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ];
3533         }
3534         else
3535         {
3536             put_record( $record, $out );
3537         }
3538     }
3539
3540     Maasha::Align::align_invert( \@entries, $options->{ "soft" } );
3541
3542     map { put_record( { SEQ_NAME => $_->[ SEQ_NAME ], SEQ => $_->[ SEQ ] }, $out ) } @entries;
3543 }
3544
3545
3546 sub script_patscan_seq
3547 {
3548     # Martin A. Hansen, August 2007.
3549
3550     # Locates patterns in sequences using scan_for_matches.
3551
3552     my ( $in,        # handle to in stream
3553          $out,       # handle to out stream
3554          $options,   # options hash
3555        ) = @_;
3556
3557     # Returns nothing.
3558
3559     my ( $genome_file, @args, $arg, $type, $seq_file, $pat_file, $out_file, $fh_in, $fh_out, $record, $patterns, $pattern, $entry, $result, %head_hash, $i );
3560
3561     if ( $options->{ "patterns" } ) {
3562         $patterns = Maasha::Patscan::parse_patterns( $options->{ "patterns" } );
3563     } elsif ( -f $options->{ "patterns_in" } ) {
3564         $patterns = Maasha::Patscan::read_patterns( $options->{ "patterns_in" } );
3565     }
3566
3567     $genome_file = "$ENV{ 'BP_DATA' }/genomes/$options->{ 'genome' }/fasta/$options->{ 'genome' }.fna" if $options->{ 'genome' };
3568
3569     push @args, "-c"                            if $options->{ "comp" };
3570     push @args, "-m $options->{ 'max_hits' }"   if $options->{ 'max_hits' };
3571     push @args, "-n $options->{ 'max_misses' }" if $options->{ 'max_hits' };
3572
3573     $seq_file = "$BP_TMP/patscan.seq";
3574     $pat_file = "$BP_TMP/patscan.pat";
3575     $out_file = "$BP_TMP/patscan.out";
3576
3577     $fh_out = Maasha::Common::write_open( $seq_file );
3578
3579     $i = 0;
3580
3581     while ( $record = get_record( $in ) )
3582     {
3583         if ( $record->{ "SEQ" } and $record->{ "SEQ_NAME" } )
3584         {
3585             $type = Maasha::Seq::seq_guess_type( $record->{ "SEQ" } ) if not $type;
3586
3587             Maasha::Fasta::put_entry( [ $i, $record->{ "SEQ" } ], $fh_out );
3588
3589             $head_hash{ $i } = $record->{ "SEQ_NAME" };
3590
3591             $i++;
3592         }
3593     }
3594
3595     close $fh_out;
3596
3597     $arg  = join " ", @args;
3598     $arg .= " -p" if $type eq "protein";
3599
3600     foreach $pattern ( @{ $patterns } )
3601     {
3602         $fh_out = Maasha::Common::write_open( $pat_file );
3603
3604         print $fh_out "$pattern\n";
3605
3606         close $fh_out;
3607
3608         if ( $options->{ 'genome' } ) {
3609             `scan_for_matches $arg $pat_file < $genome_file > $out_file`;
3610             # Maasha::Common::run( "scan_for_matches", "$arg $pat_file < $genome_file > $out_file" );
3611         } else {
3612             `scan_for_matches $arg $pat_file < $seq_file > $out_file`;
3613             # Maasha::Common::run( "scan_for_matches", "$arg $pat_file < $seq_file > $out_file" );
3614         }
3615
3616         $fh_in = Maasha::Common::read_open( $out_file );
3617
3618         while ( $entry = Maasha::Fasta::get_entry( $fh_in ) )
3619         {
3620             $result = Maasha::Patscan::parse_scan_result( $entry, $pattern );
3621
3622             if ( $options->{ 'genome' } )
3623             {
3624                 $result->{ "CHR" }     = $result->{ "S_ID" };
3625                 $result->{ "CHR_BEG" } = $result->{ "S_BEG" };
3626                 $result->{ "CHR_END" } = $result->{ "S_END" };
3627
3628                 delete $result->{ "S_ID" };
3629                 delete $result->{ "S_BEG" };
3630                 delete $result->{ "S_END" };
3631             }
3632             else
3633             {
3634                 $result->{ "S_ID" } = $head_hash{ $result->{ "S_ID" } };
3635             }
3636
3637             put_record( $result, $out );
3638         }
3639
3640         close $fh_in;
3641     }
3642
3643     unlink $pat_file;
3644     unlink $seq_file;
3645     unlink $out_file;
3646 }
3647
3648
3649 sub script_create_blast_db
3650 {
3651     # Martin A. Hansen, September 2007.
3652
3653     # Creates a NCBI BLAST database with formatdb
3654
3655     my ( $in,        # handle to in stream
3656          $out,       # handle to out stream
3657          $options,   # options hash
3658        ) = @_;
3659
3660     # Returns nothing.
3661
3662     my ( $fh, $seq_type, $path, $record, $entry );
3663
3664     $path = $options->{ "database" };
3665
3666     $fh = Maasha::Common::write_open( $path );
3667
3668     while ( $record = get_record( $in ) )
3669     {
3670         put_record( $record, $out ) if not $options->{ "no_stream" };
3671
3672         if ( $entry = Maasha::Fasta::biopiece2fasta( $record ) )
3673         {
3674             $seq_type = Maasha::Seq::seq_guess_type( $entry->[ SEQ ] ) if not $seq_type;
3675
3676             Maasha::Fasta::put_entry( $entry, $fh );
3677         }
3678     }
3679
3680     close $fh;
3681
3682     if ( $seq_type eq "protein" ) {
3683         Maasha::Common::run( "formatdb", "-p T -i $path -t $options->{ 'database' }" );
3684     } else {
3685         Maasha::Common::run( "formatdb", "-p F -i $path -t $options->{ 'database' }" );
3686     }
3687
3688     unlink $path;
3689 }
3690
3691
3692 sub script_blast_seq
3693 {
3694     # Martin A. Hansen, September 2007.
3695
3696     # BLASTs sequences in stream against a given database.
3697
3698     my ( $in,        # handle to in stream
3699          $out,       # handle to out stream
3700          $options,   # options hash
3701        ) = @_;
3702
3703     # Returns nothing.
3704
3705     my ( $genome, $q_type, $s_type, $tmp_in, $tmp_out, $fh_in, $fh_out, $record, $line, @fields, $entry );
3706
3707     $options->{ "e_val" }  = 10 if not defined $options->{ "e_val" };
3708     $options->{ "filter" } = "F";
3709     $options->{ "filter" } = "T" if $options->{ "filter" };
3710     $options->{ "cpus" } ||= 1;
3711
3712     $options->{ "database" } = "$ENV{ 'BP_DATA' }/genomes/$options->{ 'genome' }/blast/$options->{ 'genome' }.fna" if $options->{ 'genome' };
3713
3714     $tmp_in  = "$BP_TMP/blast_query.seq";
3715     $tmp_out = "$BP_TMP/blast.result";
3716
3717     $fh_out = Maasha::Filesys::file_write_open( $tmp_in );
3718
3719     while ( $record = get_record( $in ) )
3720     {
3721         if ( $entry = Maasha::Fasta::biopiece2fasta( $record ) )
3722         {
3723             $q_type = Maasha::Seq::seq_guess_type( $entry->[ SEQ ] ) if not $q_type;
3724
3725             Maasha::Fasta::put_entry( $entry, $fh_out );
3726         }
3727
3728         put_record( $record, $out );
3729     }
3730
3731     close $fh_out;
3732
3733     if ( -f $options->{ 'database' } . ".phr" ) {
3734         $s_type = "protein";
3735     } else {
3736         $s_type = "nucleotide";
3737     }
3738
3739     if ( not $options->{ 'program' } )
3740     {
3741         if ( $q_type ne "protein" and $s_type ne "protein" ) {
3742             $options->{ 'program' } = "blastn";
3743         } elsif ( $q_type eq "protein" and $s_type eq "protein" ) {
3744             $options->{ 'program' } = "blastp";
3745         } elsif ( $q_type ne "protein" and $s_type eq "protein" ) {
3746             $options->{ 'program' } = "blastx";
3747         } elsif ( $q_type eq "protein" and $s_type ne "protein" ) {
3748             $options->{ 'program' } = "tblastn";
3749         }
3750     }
3751
3752     if ( $options->{ 'verbose' } )
3753     {
3754         Maasha::Common::run(
3755             "blastall",
3756             join( " ",
3757                 "-p $options->{ 'program' }",
3758                 "-e $options->{ 'e_val' }",
3759                 "-a $options->{ 'cpus' }",
3760                 "-m 8",
3761                 "-i $tmp_in",
3762                 "-d $options->{ 'database' }",
3763                 "-F $options->{ 'filter' }",
3764                 "-o $tmp_out",
3765             ),
3766             1
3767         );
3768     }
3769     else
3770     {
3771         Maasha::Common::run(
3772             "blastall",
3773             join( " ",
3774                 "-p $options->{ 'program' }",
3775                 "-e $options->{ 'e_val' }",
3776                 "-a $options->{ 'cpus' }",
3777                 "-m 8",
3778                 "-i $tmp_in",
3779                 "-d $options->{ 'database' }",
3780                 "-F $options->{ 'filter' }",
3781                 "-o $tmp_out",
3782                 "> /dev/null 2>&1"
3783             ),
3784             1
3785         );
3786     }
3787
3788     unlink $tmp_in;
3789
3790     $fh_out = Maasha::Filesys::file_read_open( $tmp_out );
3791
3792     undef $record;
3793
3794     while ( $line = <$fh_out> )
3795     {
3796         chomp $line;
3797
3798         next if $line =~ /^#/;
3799
3800         @fields = split /\s+/, $line;
3801
3802         $record->{ "REC_TYPE" }   = "BLAST";
3803         $record->{ "Q_ID" }       = $fields[ 0 ];
3804         $record->{ "S_ID" }       = $fields[ 1 ];
3805         $record->{ "IDENT" }      = $fields[ 2 ];
3806         $record->{ "ALIGN_LEN" }  = $fields[ 3 ];
3807         $record->{ "MISMATCHES" } = $fields[ 4 ];
3808         $record->{ "GAPS" }       = $fields[ 5 ];
3809         $record->{ "Q_BEG" }      = $fields[ 6 ] - 1; # BLAST is 1-based
3810         $record->{ "Q_END" }      = $fields[ 7 ] - 1; # BLAST is 1-based
3811         $record->{ "S_BEG" }      = $fields[ 8 ] - 1; # BLAST is 1-based
3812         $record->{ "S_END" }      = $fields[ 9 ] - 1; # BLAST is 1-based
3813         $record->{ "E_VAL" }      = $fields[ 10 ];
3814         $record->{ "BIT_SCORE" }  = $fields[ 11 ];
3815
3816         if ( $record->{ "S_BEG" } > $record->{ "S_END" } )
3817         {
3818             $record->{ "STRAND" } = '-';
3819
3820             ( $record->{ "S_BEG" }, $record->{ "S_END" } ) = ( $record->{ "S_END" }, $record->{ "S_BEG" } );
3821         }
3822         else
3823         {
3824             $record->{ "STRAND" } = '+';
3825         }
3826
3827         put_record( $record, $out );
3828     }
3829
3830     close $fh_out;
3831
3832     unlink $tmp_out;
3833 }
3834
3835
3836 sub script_blat_seq
3837 {
3838     # Martin A. Hansen, August 2007.
3839
3840     # BLATs sequences in stream against a given genome.
3841
3842     my ( $in,        # handle to in stream
3843          $out,       # handle to out stream
3844          $options,   # options hash
3845        ) = @_;
3846
3847     # Returns nothing.
3848
3849     my ( $blat_args, $genome_file, $query_file, $fh_in, $fh_out, $type, $record, $result_file, $entries, $entry );
3850
3851     $genome_file = "$ENV{ 'BP_DATA' }/genomes/$options->{ 'genome' }/fasta/$options->{ 'genome' }.fna";
3852
3853     $options->{ 'tile_size' }    ||= 11;
3854     $options->{ 'one_off' }      ||= 0;
3855     $options->{ 'min_identity' } ||= 90;
3856     $options->{ 'min_score' }    ||= 0;
3857     $options->{ 'step_size' }    ||= $options->{ 'tile_size' };
3858
3859     $blat_args .= " -tileSize=$options->{ 'tile_size' }";
3860     $blat_args .= " -oneOff=$options->{ 'one_off' }";
3861     $blat_args .= " -minIdentity=$options->{ 'min_identity' }";
3862     $blat_args .= " -minScore=$options->{ 'min_score' }";
3863     $blat_args .= " -stepSize=$options->{ 'step_size' }";
3864 #    $blat_args .= " -ooc=" . Maasha::Config::genome_blat_ooc( $options->{ "genome" }, 11 ) if $options->{ 'ooc' };
3865
3866     $query_file = "$BP_TMP/blat.seq";
3867
3868     $fh_out = Maasha::Common::write_open( $query_file );
3869
3870     while ( $record = get_record( $in ) )
3871     {
3872         if ( $entry = Maasha::Fasta::biopiece2fasta( $record ) )
3873         {
3874             $type = Maasha::Seq::seq_guess_type( $entry->[ SEQ ] ) if not $type;
3875             Maasha::Fasta::put_entry( $entry, $fh_out, 80 );
3876         }
3877
3878         put_record( $record, $out );
3879     }
3880
3881     close $fh_out;
3882
3883     $blat_args .= " -t=dnax" if $type eq "protein";
3884     $blat_args .= " -q=$type";
3885
3886     $result_file = "$BP_TMP/blat.psl";
3887
3888     Maasha::Common::run( "blat", "$genome_file $query_file $blat_args $result_file > /dev/null 2>&1" );
3889
3890     unlink $query_file;
3891
3892     $entries = Maasha::UCSC::psl_get_entries( $result_file );
3893
3894     map { put_record( $_, $out ) } @{ $entries };
3895
3896     unlink $result_file;
3897 }
3898
3899
3900 sub script_soap_seq
3901 {
3902     # Martin A. Hansen, July 2008.
3903
3904     # soap sequences in stream against a given file or genome.
3905
3906     my ( $in,        # handle to in stream
3907          $out,       # handle to out stream
3908          $options,   # options hash
3909        ) = @_;
3910
3911     # Returns nothing.
3912
3913     my ( $genome, $tmp_in, $tmp_out, $fh_in, $fh_out, $record, $line, @fields, $entry, $count, $args );
3914
3915     $options->{ "seed_size" }  ||= 10;
3916     $options->{ "mismatches" } ||= 2;
3917     $options->{ "gap_size" }   ||= 0;
3918     $options->{ "cpus" }       ||= 1;
3919
3920     if ( $options->{ "genome" } ) {
3921         $options->{ "in_file" } = "$ENV{ 'BP_DATA' }/genomes/$options->{ 'genome' }/fasta/$options->{ 'genome' }.fna";
3922     }
3923
3924     $tmp_in  = "$BP_TMP/soap_query.seq";
3925     $tmp_out = "$BP_TMP/soap.result";
3926
3927     $fh_out = Maasha::Common::write_open( $tmp_in );
3928
3929     $count = 0;
3930
3931     while ( $record = get_record( $in ) )
3932     {
3933         if ( $entry = Maasha::Fasta::biopiece2fasta( $record ) )
3934         {
3935             Maasha::Fasta::put_entry( $entry, $fh_out );
3936
3937             $count++;
3938         }
3939
3940         put_record( $record, $out );
3941     }
3942
3943     close $fh_out;
3944
3945     if ( $count > 0 )
3946     {
3947         $args = join( " ",
3948             "-s $options->{ 'seed_size' }",
3949             "-r 2",
3950             "-a $tmp_in",
3951             "-v $options->{ 'mismatches' }",
3952             "-g $options->{ 'gap_size' }",
3953             "-p $options->{ 'cpus' }",
3954             "-d $options->{ 'in_file' }",
3955             "-o $tmp_out",
3956         );
3957
3958         $args .= " > /dev/null 2>&1" if not $options->{ 'verbose' };
3959
3960         Maasha::Common::run( "soap", $args, 1 );
3961
3962         unlink $tmp_in;
3963
3964         $fh_out = Maasha::Common::read_open( $tmp_out );
3965
3966         undef $record;
3967
3968         while ( $line = <$fh_out> )
3969         {
3970             chomp $line;
3971
3972             @fields = split /\t/, $line;
3973
3974             $record->{ "REC_TYPE" }   = "SOAP";
3975             $record->{ "Q_ID" }       = $fields[ 0 ];
3976             $record->{ "SCORE" }      = $fields[ 3 ];
3977             $record->{ "STRAND" }     = $fields[ 6 ];
3978             $record->{ "S_ID" }       = $fields[ 7 ];
3979             $record->{ "S_BEG" }      = $fields[ 8 ] - 1; # soap is 1-based
3980             $record->{ "S_END" }      = $fields[ 8 ] + $fields[ 5 ] - 2;
3981
3982             put_record( $record, $out );
3983         }
3984
3985         close $fh_out;
3986     }
3987
3988     unlink $tmp_out;
3989 }
3990
3991
3992 sub script_match_seq
3993 {
3994     # Martin A. Hansen, August 2007.
3995
3996     # BLATs sequences in stream against a given genome.
3997
3998     my ( $in,        # handle to in stream
3999          $out,       # handle to out stream
4000          $options,   # options hash
4001        ) = @_;
4002
4003     # Returns nothing.
4004
4005     my ( $record, @entries, $results );
4006
4007     $options->{ "word_size" } ||= 20;
4008     $options->{ "direction" } ||= "both";
4009
4010     while ( $record = get_record( $in ) )
4011     {
4012         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } ) {
4013             push @entries, [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ];
4014         }
4015
4016         put_record( $record, $out );
4017     }
4018
4019     if ( @entries == 1 )
4020     {
4021         $results = Maasha::Match::match_mummer( [ $entries[ 0 ] ], [ $entries[ 0 ] ], $options, $BP_TMP );
4022
4023         map { put_record( $_, $out ) } @{ $results };
4024     }
4025     elsif ( @entries == 2 )
4026     {
4027         $results = Maasha::Match::match_mummer( [ $entries[ 0 ] ], [ $entries[ 1 ] ], $options, $BP_TMP );
4028
4029         map { put_record( $_, $out ) } @{ $results };
4030     }
4031 }
4032
4033
4034 sub script_create_vmatch_index
4035 {
4036     # Martin A. Hansen, January 2008.
4037
4038     # Create a vmatch index from sequences in the stream.
4039
4040     my ( $in,        # handle to in stream
4041          $out,       # handle to out stream
4042          $options,   # options hash
4043        ) = @_;
4044
4045     # Returns nothing.
4046
4047     my ( $record, $file_tmp, $fh_tmp, $type, $entry );
4048
4049     if ( $options->{ "index_name" } )
4050     {
4051         $file_tmp = $options->{ 'index_name' };
4052         $fh_tmp   = Maasha::Common::write_open( $file_tmp );
4053     }
4054
4055     while ( $record = get_record( $in ) )
4056     {
4057         if ( $options->{ "index_name" } and $entry = Maasha::Fasta::biopiece2fasta( $record ) )
4058         {
4059             Maasha::Fasta::put_entry( $entry, $fh_tmp );
4060
4061             $type = Maasha::Seq::seq_guess_type( $entry->[ SEQ ] ) if not defined $type;
4062         }
4063
4064         put_record( $record, $out ) if not $options->{ "no_stream" };
4065     }
4066
4067     if ( $options->{ "index_name" } )
4068     {
4069         close $fh_tmp;
4070
4071         if ( $type eq "protein" ) {
4072             Maasha::Common::run( "mkvtree", "-db $file_tmp -protein -pl $options->{ 'prefix_length' } -allout -indexname $file_tmp > /dev/null 2>&1" );
4073         } else {
4074             Maasha::Common::run( "mkvtree", "-db $file_tmp -dna -pl $options->{ 'prefix_length' } -allout -indexname $file_tmp > /dev/null 2>&1" );
4075         }
4076
4077         unlink $file_tmp;
4078     }
4079 }
4080
4081
4082 sub script_vmatch_seq
4083 {
4084     # Martin A. Hansen, August 2007.
4085
4086     # Vmatches sequences in stream against a given genome.
4087
4088     my ( $in,        # handle to in stream
4089          $out,       # handle to out stream
4090          $options,   # options hash
4091        ) = @_;
4092
4093     # Returns nothing.
4094
4095     my ( @index_files, @records, $result_file, $fh_in, $record, %hash );
4096
4097     $options->{ 'count' } = 1 if $options->{ 'max_hits' };
4098
4099     if ( $options->{ "index_name" } )
4100     {
4101         @index_files = $options->{ "index_name" };
4102     }
4103     else
4104     {
4105         @index_files = Maasha::Common::ls_files( "$ENV{ 'BP_DATA' }/genomes/$options->{ 'genome' }/vmatch" );
4106
4107         map { $_ =~ /^(.+)\.[a-z1]{3,4}$/; $hash{ $1 } = 1 } @index_files;
4108
4109         @index_files = sort keys %hash;
4110     }
4111
4112     while ( $record = get_record( $in ) )
4113     {
4114         push @records, $record;
4115
4116         put_record( $record, $out );
4117     }
4118
4119     $result_file = Maasha::Match::match_vmatch( $BP_TMP, \@records, \@index_files, $options );
4120
4121     undef @records;
4122
4123     $fh_in = Maasha::Common::read_open( $result_file );
4124
4125     while ( $record = Maasha::Match::vmatch_get_entry( $fh_in ) ) {
4126         put_record( $record, $out );
4127     }
4128
4129     close $fh_in;
4130
4131     unlink $result_file;
4132 }
4133
4134
4135 sub script_write_fasta
4136 {
4137     # Martin A. Hansen, August 2007.
4138
4139     # Write FASTA entries from sequences in stream.
4140
4141     my ( $in,        # handle to in stream
4142          $out,       # handle to out stream
4143          $options,   # options hash
4144        ) = @_;
4145
4146     # Returns nothing.
4147
4148     my ( $record, $fh, $entry );
4149
4150     $fh = write_stream( $options->{ "data_out" }, $options->{ "compress" } );
4151
4152     while ( $record = get_record( $in ) )
4153     {
4154         if ( $entry = Maasha::Fasta::biopiece2fasta( $record ) ) {
4155             Maasha::Fasta::put_entry( $entry, $fh, $options->{ "wrap" } );
4156         }
4157
4158         put_record( $record, $out ) if not $options->{ "no_stream" };
4159     }
4160
4161     close $fh;
4162 }
4163
4164
4165 sub script_write_align
4166 {
4167     # Martin A. Hansen, August 2007.
4168
4169     # Write pretty alignments aligned sequences in stream.
4170
4171     my ( $in,        # handle to in stream
4172          $out,       # handle to out stream
4173          $options,   # options hash
4174        ) = @_;
4175
4176     # Returns nothing.
4177
4178     my ( $fh, $record, @entries );
4179
4180     $fh = write_stream( $options->{ "data_out" } ) ;
4181
4182     while ( $record = get_record( $in ) )
4183     {
4184         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } ) {
4185             push @entries, [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ];
4186         }
4187
4188         put_record( $record, $out ) if not $options->{ "no_stream" };
4189     }
4190
4191     if ( scalar( @entries ) == 2 ) {
4192         Maasha::Align::align_print_pairwise( $entries[ 0 ], $entries[ 1 ], $fh, $options->{ "wrap" } );
4193     } elsif ( scalar ( @entries ) > 2 ) {
4194         Maasha::Align::align_print_multi( \@entries, $fh, $options->{ "wrap" }, $options->{ "no_ruler" }, $options->{ "no_consensus" } );
4195     }
4196
4197     close $fh if $fh;
4198 }
4199
4200
4201 sub script_write_blast
4202 {
4203     # Martin A. Hansen, November 2007.
4204
4205     # Write data in blast table format (-m8 and 9).
4206
4207     my ( $in,        # handle to in stream
4208          $out,       # handle to out stream
4209          $options,   # options hash
4210        ) = @_;
4211
4212     # Returns nothing.
4213
4214     my ( $fh, $record, $first );
4215
4216     $fh = write_stream( $options->{ "data_out" }, $options->{ "compress" } ) ;
4217
4218     $first = 1;
4219
4220     while ( $record = get_record( $in ) )
4221     {
4222         if ( $record->{ "REC_TYPE" } eq "BLAST" )
4223         {
4224             if ( $options->{ "comment" } and $first )
4225             {
4226                 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";
4227
4228                 $first = 0;
4229             }
4230
4231             if ( $record->{ "STRAND" } eq "-" ) {
4232                 ( $record->{ "S_BEG" }, $record->{ "S_END" } ) = ( $record->{ "S_END" }, $record->{ "S_BEG" } );
4233             }
4234
4235             print $fh join( "\t",
4236                 $record->{ "Q_ID" },
4237                 $record->{ "S_ID" },
4238                 $record->{ "IDENT" },
4239                 $record->{ "ALIGN_LEN" },
4240                 $record->{ "MISMATCHES" },
4241                 $record->{ "GAPS" },
4242                 $record->{ "Q_BEG" } + 1,
4243                 $record->{ "Q_END" } + 1,
4244                 $record->{ "S_BEG" } + 1,
4245                 $record->{ "S_END" } + 1,
4246                 $record->{ "E_VAL" },
4247                 $record->{ "BIT_SCORE" }
4248             ), "\n";
4249         }
4250
4251         put_record( $record, $out ) if not $options->{ "no_stream" };
4252     }
4253
4254     close $fh;
4255 }
4256
4257
4258 sub script_write_tab
4259 {
4260     # Martin A. Hansen, August 2007.
4261
4262     # Write data as table.
4263
4264     my ( $in,        # handle to in stream
4265          $out,       # handle to out stream
4266          $options,   # options hash
4267        ) = @_;
4268
4269     # Returns nothing.
4270
4271     my ( $fh, $record, $key, @keys, @vals, $ok, %no_keys, $A, $B );
4272
4273     $options->{ "delimit" } ||= "\t";
4274
4275     map { $no_keys{ $_ } = 1 } @{ $options->{ "no_keys" } };
4276
4277     $fh = write_stream( $options->{ "data_out" }, $options->{ "compress" } );
4278
4279     while ( $record = get_record( $in ) )
4280     {
4281         undef @vals;
4282         $ok = 1;
4283
4284         if ( $options->{ "keys" } )
4285         {
4286             map { $ok = 0 if not exists $record->{ $_ } } @{ $options->{ "keys" } };
4287
4288             if ( $ok )
4289             {
4290                 foreach $key ( @{ $options->{ "keys" }  } )
4291                 {
4292                     if ( exists $record->{ $key } )
4293                     {
4294                         push @keys, $key if $options->{ "comment" };
4295                         push @vals, $record->{ $key };
4296                     }
4297                 }
4298              }
4299         }
4300         else
4301         {
4302             foreach $key ( sort { $A = $a; $B = $b; $A =~ s/^V(\d+)$/$1/; $B =~ s/^V(\d+)$/$1/; $A <=> $B } keys %{ $record } )
4303             {
4304                 next if exists $no_keys{ $key };
4305
4306                 push @keys, $key if $options->{ "comment" };
4307                 push @vals, $record->{ $key };
4308             }
4309         }
4310
4311         if ( @keys and $options->{ "comment" } )
4312         {
4313             print $fh "#", join( $options->{ "delimit" }, @keys ), "\n";
4314
4315             delete $options->{ "comment" };
4316         }
4317
4318         print $fh join( $options->{ "delimit" }, @vals ), "\n" if @vals;
4319
4320         put_record( $record, $out ) if not $options->{ "no_stream" };
4321     }
4322
4323     close $fh;
4324 }
4325
4326
4327 sub script_write_bed
4328 {
4329     # Martin A. Hansen, August 2007.
4330
4331     # Write BED format for the UCSC genome browser using records in stream.
4332
4333     my ( $in,        # handle to in stream
4334          $out,       # handle to out stream
4335          $options,   # options hash
4336        ) = @_;
4337
4338     # Returns nothing.
4339
4340     my ( $cols, $fh, $record, $bed_entry, $new_record );
4341
4342     $cols = $options->{ 'cols' }->[ 0 ];
4343
4344     $fh = write_stream( $options->{ 'data_out' }, $options->{ 'compress' } );
4345
4346     while ( $record = get_record( $in ) )
4347     {
4348         $record = Maasha::UCSC::psl2record( $record ) if $record->{ 'tBaseInsert' }; # Dirty addition to allow Affy data from MySQL to be dumped
4349
4350         if ( $bed_entry = Maasha::UCSC::BED::biopiece2bed( $record, $cols ) ) {
4351             Maasha::UCSC::BED::bed_entry_put( $bed_entry, $fh, $cols, $options->{ 'check' } );
4352         }
4353
4354         put_record( $record, $out ) if not $options->{ 'no_stream' };
4355     }
4356
4357     close $fh;
4358 }
4359
4360
4361 sub script_write_psl
4362 {
4363     # Martin A. Hansen, August 2007.
4364
4365     # Write PSL output from stream.
4366
4367     my ( $in,        # handle to in stream
4368          $out,       # handle to out stream
4369          $options,   # options hash
4370        ) = @_;
4371
4372     # Returns nothing.
4373
4374     my ( $fh, $record, @output, $first );
4375
4376     $first = 1;
4377
4378     $fh = write_stream( $options->{ "data_out" }, $options->{ "compress" } );
4379
4380     while ( $record = get_record( $in ) )
4381     {
4382         put_record( $record, $out ) if not $options->{ "no_stream" };
4383
4384         if ( $record->{ "REC_TYPE" } and $record->{ "REC_TYPE" } eq "PSL" )
4385         {
4386             Maasha::UCSC::psl_put_header( $fh ) if $first;
4387             Maasha::UCSC::psl_put_entry( $record, $fh );
4388             $first = 0;
4389         }
4390     }
4391
4392     close $fh;
4393 }
4394
4395
4396 sub script_write_fixedstep
4397 {
4398     # Martin A. Hansen, Juli 2008.
4399
4400     # Write fixedStep entries from recrods in the stream.
4401
4402     my ( $in,        # handle to in stream
4403          $out,       # handle to out stream
4404          $options,   # options hash
4405        ) = @_;
4406
4407     # Returns nothing.
4408
4409     my ( $fh, $record, $entry );
4410
4411     $fh = write_stream( $options->{ "data_out" }, $options->{ "compress" } );
4412
4413     while ( $record = get_record( $in ) )
4414     {
4415         if ( $entry = Maasha::UCSC::Wiggle::biopiece2fixedstep( $record ) ) {
4416             Maasha::UCSC::Wiggle::fixedstep_entry_put( $entry, $fh );
4417         }
4418
4419         put_record( $record, $out ) if not $options->{ "no_stream" };
4420     }
4421
4422     close $fh;
4423 }
4424
4425
4426 sub script_write_2bit
4427 {
4428     # Martin A. Hansen, March 2008.
4429
4430     # Write sequence entries from stream in 2bit format.
4431
4432     my ( $in,        # handle to in stream
4433          $out,       # handle to out stream
4434          $options,   # options hash
4435        ) = @_;
4436
4437     # Returns nothing.
4438
4439     my ( $record, $mask, $tmp_file, $fh_tmp, $fh_in, $fh_out, $entry );
4440
4441     $mask = 1 if not $options->{ "no_mask" };
4442
4443     $tmp_file = "$BP_TMP/write_2bit.fna";
4444     $fh_tmp   = Maasha::Common::write_open( $tmp_file );
4445
4446     $fh_out = write_stream( $options->{ "data_out" } );
4447
4448     while ( $record = get_record( $in ) )
4449     {
4450         if ( $entry = Maasha::Fasta::biopiece2fasta( $record ) ) {
4451             Maasha::Fasta::put_entry( $entry, $fh_tmp );
4452         }
4453
4454         put_record( $record, $out ) if not $options->{ "no_stream" };
4455     }
4456
4457     close $fh_tmp;
4458
4459     $fh_in = Maasha::Common::read_open( $tmp_file );
4460
4461     Maasha::TwoBit::fasta2twobit( $fh_in, $fh_out, $mask );
4462
4463     close $fh_in;
4464     close $fh_out;
4465
4466     unlink $tmp_file;
4467 }
4468
4469
4470 sub script_write_solid
4471 {
4472     # Martin A. Hansen, April 2008.
4473
4474     # Write di-base encoded Solid sequence from entries in stream.
4475
4476     my ( $in,        # handle to in stream
4477          $out,       # handle to out stream
4478          $options,   # options hash
4479        ) = @_;
4480
4481     # Returns nothing.
4482
4483     my ( $record, $fh, $entry );
4484
4485     $fh = write_stream( $options->{ "data_out" }, $options->{ "compress" } );
4486
4487     while ( $record = get_record( $in ) )
4488     {
4489         if ( $entry = Maasha::Fasta::biopiece2fasta( $record ) )
4490         {
4491             $entry->[ SEQ ] = Maasha::Solid::seq2color_space( uc $entry->[ SEQ ] );
4492
4493             Maasha::Fasta::put_entry( $entry, $fh, $options->{ "wrap" } );
4494         }
4495
4496         put_record( $record, $out ) if not $options->{ "no_stream" };
4497     }
4498
4499     close $fh;
4500 }
4501
4502
4503 sub script_write_ucsc_config
4504 {
4505     # Martin A. Hansen, November 2008.
4506
4507     # Write UCSC Genome Broser configuration (.ra file type) from
4508     # records in the stream.
4509
4510     my ( $in,        # handle to in stream
4511          $out,       # handle to out stream
4512          $options,   # options hash
4513        ) = @_;
4514
4515     # Returns nothing.
4516
4517     my ( $record, $fh );
4518
4519     $fh = write_stream( $options->{ "data_out" } );
4520
4521     while ( $record = get_record( $in ) )
4522     {
4523         Maasha::UCSC::ucsc_config_entry_put( $record, $fh ) if $record->{ "REC_TYPE" } eq "UCSC Config";
4524
4525         put_record( $record, $out ) if not $options->{ "no_stream" };
4526     }
4527
4528     close $fh;
4529 }
4530
4531
4532 sub script_plot_seqlogo
4533 {
4534     # Martin A. Hansen, August 2007.
4535
4536     # Calculates and writes a sequence logo for alignments.
4537
4538     my ( $in,        # handle to in stream
4539          $out,       # handle to out stream
4540          $options,   # options hash
4541        ) = @_;
4542
4543     # Returns nothing.
4544
4545     my ( $record, @entries, $logo, $fh );
4546
4547     while ( $record = get_record( $in ) )
4548     {
4549         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } ) {
4550             push @entries, [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ];
4551         }
4552
4553         put_record( $record, $out ) if not $options->{ "no_stream" };
4554     }
4555
4556     $logo = Maasha::Plot::seq_logo( \@entries );
4557
4558     $fh = write_stream( $options->{ "data_out" } );
4559
4560     print $fh $logo;
4561
4562     close $fh;
4563 }
4564
4565
4566 sub script_plot_phastcons_profiles
4567 {
4568     # Martin A. Hansen, January 2008.
4569
4570     # Plots PhastCons profiles.
4571
4572     my ( $in,        # handle to in stream
4573          $out,       # handle to out stream
4574          $options,   # options hash
4575        ) = @_;
4576
4577     # Returns nothing.
4578
4579     my ( $phastcons_file, $phastcons_index, $index, $fh_phastcons, $record, $scores, $AoA, $plot, $fh );
4580
4581     $options->{ "title" } ||= "PhastCons Profiles";
4582
4583     $phastcons_file  = Maasha::Config::genome_phastcons( $options->{ "genome" } );
4584     $phastcons_index = Maasha::Config::genome_phastcons_index( $options->{ "genome" } );
4585
4586     $index           = Maasha::UCSC::fixedstep_index_retrieve( $phastcons_index );
4587     $fh_phastcons    = Maasha::Common::read_open( $phastcons_file );
4588
4589     while ( $record = get_record( $in ) )
4590     {
4591         if ( $record->{ "CHR" } and $record->{ "CHR_BEG" } and $record->{ "CHR_END" } )
4592         {
4593             $scores = Maasha::UCSC::fixedstep_index_lookup( $index, $fh_phastcons, $record->{ "CHR" },
4594                                                                                    $record->{ "CHR_BEG" },
4595                                                                                    $record->{ "CHR_END" },
4596                                                                                    $options->{ "flank" } );
4597
4598             push @{ $AoA }, [ @{ $scores } ];
4599         }
4600
4601         put_record( $record, $out ) if not $options->{ "no_stream" };
4602     }
4603
4604     Maasha::UCSC::phastcons_normalize( $AoA );
4605
4606     $AoA = [ [ Maasha::UCSC::phastcons_mean( $AoA ) ] ] if $options->{ "mean" };
4607     $AoA = [ [ Maasha::UCSC::phastcons_median( $AoA ) ] ] if $options->{ "median" };
4608
4609     $AoA = Maasha::Matrix::matrix_flip( $AoA );
4610
4611     $plot = Maasha::Plot::lineplot_simple( $AoA, $options, $BP_TMP );
4612
4613     $fh = write_stream( $options->{ "data_out" } );
4614
4615     print $fh "$_\n" foreach @{ $plot };
4616
4617     close $fh;
4618 }
4619
4620
4621 sub script_analyze_bed
4622 {
4623     # Martin A. Hansen, March 2008.
4624
4625     # Analyze BED entries in stream.
4626
4627     my ( $in,        # handle to in stream
4628          $out,       # handle to out stream
4629          $options,   # options hash
4630        ) = @_;
4631
4632     # Returns nothing.
4633
4634     my ( $record );
4635
4636     while ( $record = get_record( $in ) )
4637     {
4638         $record = Maasha::UCSC::bed_analyze( $record ) if $record->{ "REC_TYPE" } eq "BED";
4639
4640         put_record( $record, $out );
4641     }
4642 }
4643
4644
4645 sub script_analyze_vals
4646 {
4647     # Martin A. Hansen, August 2007.
4648
4649     # Analyze values for given keys in stream.
4650
4651     my ( $in,        # handle to in stream
4652          $out,       # handle to out stream
4653          $options,   # options hash
4654        ) = @_;
4655
4656     # Returns nothing.
4657
4658     my ( $record, $key, @keys, %key_hash, $analysis, $len );
4659
4660     map { $key_hash{ $_ } = 1 } @{ $options->{ "keys" } };
4661
4662     while ( $record = get_record( $in ) )
4663     {
4664         foreach $key ( keys %{ $record } )
4665         {
4666             next if $options->{ "keys" } and not exists $key_hash{ $key };
4667
4668             $analysis->{ $key }->{ "COUNT" }++;
4669
4670             if ( Maasha::Calc::is_a_number( $record->{ $key } ) )
4671             {
4672                 $analysis->{ $key }->{ "TYPE" } = "num";
4673                 $analysis->{ $key }->{ "SUM" } += $record->{ $key };
4674                 $analysis->{ $key }->{ "MAX" } = $record->{ $key } if $record->{ $key } > $analysis->{ $key }->{ "MAX" } or not $analysis->{ $key }->{ "MAX" };
4675                 $analysis->{ $key }->{ "MIN" } = $record->{ $key } if $record->{ $key } < $analysis->{ $key }->{ "MIN" } or not $analysis->{ $key }->{ "MIN" };
4676             }
4677             else
4678             {
4679                 $len = length $record->{ $key };
4680
4681                 $analysis->{ $key }->{ "TYPE" } = "alph";
4682                 $analysis->{ $key }->{ "SUM" } += $len;
4683                 $analysis->{ $key }->{ "MAX" } = $len if $len > $analysis->{ $key }->{ "MAX" } or not $analysis->{ $key }->{ "MAX" };
4684                 $analysis->{ $key }->{ "MIN" } = $len if $len < $analysis->{ $key }->{ "MIM" } or not $analysis->{ $key }->{ "MIN" };
4685             }
4686         }
4687
4688         put_record( $record, $out ) if not $options->{ "no_stream" };
4689     }
4690
4691     foreach $key ( keys %{ $analysis } )
4692     {
4693         $analysis->{ $key }->{ "MEAN" } = sprintf "%.2f", $analysis->{ $key }->{ "SUM" } / $analysis->{ $key }->{ "COUNT" };
4694         $analysis->{ $key }->{ "SUM" }  = sprintf "%.2f", $analysis->{ $key }->{ "SUM" };
4695     }
4696
4697     my ( $keys, $types, $counts, $mins, $maxs, $sums, $means );
4698
4699     $keys   = "KEY  ";
4700     $types  = "TYPE ";
4701     $counts = "COUNT";
4702     $mins   = "MIN  ";
4703     $maxs   = "MAX  ";
4704     $sums   = "SUM  ";
4705     $means  = "MEAN ";
4706
4707     if ( $options->{ "keys" } ) {
4708         @keys = @{ $options->{ "keys" } };
4709     } else {
4710         @keys = keys %{ $analysis };
4711     }
4712
4713     foreach $key ( @keys )
4714     {
4715         $keys   .= sprintf "% 15s", $key;
4716         $types  .= sprintf "% 15s", $analysis->{ $key }->{ "TYPE" };
4717         $counts .= sprintf "% 15s", $analysis->{ $key }->{ "COUNT" };
4718         $mins   .= sprintf "% 15s", $analysis->{ $key }->{ "MIN" };
4719         $maxs   .= sprintf "% 15s", $analysis->{ $key }->{ "MAX" };
4720         $sums   .= sprintf "% 15s", $analysis->{ $key }->{ "SUM" };
4721         $means  .= sprintf "% 15s", $analysis->{ $key }->{ "MEAN" };
4722     }
4723
4724     print $out "$keys\n";
4725     print $out "$types\n";
4726     print $out "$counts\n";
4727     print $out "$mins\n";
4728     print $out "$maxs\n";
4729     print $out "$sums\n";
4730     print $out "$means\n";
4731 }
4732
4733
4734 sub script_head_records
4735 {
4736     # Martin A. Hansen, August 2007.
4737
4738     # Display the first sequences in stream.
4739
4740     my ( $in,        # handle to in stream
4741          $out,       # handle to out stream
4742          $options,   # options hash
4743        ) = @_;
4744
4745     # Returns nothing.
4746
4747     my ( $record, $count );
4748
4749     $options->{ "num" } ||= 10;
4750
4751     $count = 0;
4752
4753     while ( $record = get_record( $in ) )
4754     {
4755         $count++;
4756
4757         put_record( $record, $out );
4758
4759         last if $count == $options->{ "num" };
4760     }
4761 }
4762
4763
4764 sub script_remove_keys
4765 {
4766     # Martin A. Hansen, August 2007.
4767
4768     # Remove keys from stream.
4769
4770     my ( $in,        # handle to in stream
4771          $out,       # handle to out stream
4772          $options,   # options hash
4773        ) = @_;
4774
4775     # Returns nothing.
4776
4777     my ( $record, $new_record );
4778
4779     while ( $record = get_record( $in ) )
4780     {
4781         if ( $options->{ "keys" } )
4782         {
4783             map { delete $record->{ $_ } } @{ $options->{ "keys" } };
4784         }
4785         elsif ( $options->{ "save_keys" } )
4786         {
4787             map { $new_record->{ $_ } = $record->{ $_ } if exists $record->{ $_ } } @{ $options->{ "save_keys" } };
4788
4789             $record = $new_record;
4790         }
4791
4792         put_record( $record, $out ) if keys %{ $record };
4793     }
4794 }
4795
4796
4797 sub script_remove_adaptor
4798 {
4799     # Martin A. Hansen, August 2008.
4800
4801     # Find and remove adaptor from sequences in the stream.
4802
4803     my ( $in,        # handle to in stream
4804          $out,       # handle to out stream
4805          $options,   # options hash
4806        ) = @_;
4807
4808     # Returns nothing.
4809
4810     my ( $record, $adaptor, $seq, $adaptor_len, $seq_len, $offset, $max_match, $max_mismatch, $pos );
4811
4812     $options->{ "remove" } ||= "after";
4813
4814     $max_mismatch = $options->{ "mismatches" } || 0;
4815     $offset       = $options->{ "offset" };
4816
4817     if ( not defined $offset ) {
4818         $offset = 0;
4819     } else {
4820         $offset--;
4821     }
4822
4823     $adaptor      = uc $options->{ "adaptor" };
4824     $adaptor_len  = length $adaptor;
4825
4826     while ( $record = get_record( $in ) )
4827     {
4828         if ( $record->{ "SEQ" } )
4829         {
4830             $seq     = uc $record->{ "SEQ" };
4831             $seq_len = length $seq;
4832
4833             $pos = Maasha::Common::index_m( $seq, $adaptor, $seq_len, $adaptor_len, $offset, $max_mismatch );
4834
4835             $record->{ "ADAPTOR_POS" } = $pos;
4836
4837             if ( $pos >= 0 and $options->{ "remove" } ne "skip" )
4838             {
4839                 if ( $options->{ "remove" } eq "after" )
4840                 {
4841                     $record->{ "SEQ" }     = substr $record->{ "SEQ" }, 0, $pos;
4842                     $record->{ "SEQ_LEN" } = $pos;
4843                 }
4844                 else
4845                 {
4846                     $record->{ "SEQ" }     = substr $record->{ "SEQ" }, $pos + $adaptor_len;
4847                     $record->{ "SEQ_LEN" } = length $record->{ "SEQ" };
4848                 }
4849             }
4850
4851             put_record( $record, $out );
4852         }
4853         else
4854         {
4855             put_record( $record, $out );
4856         }
4857     }
4858 }
4859
4860
4861 sub script_remove_mysql_tables
4862 {
4863     # Martin A. Hansen, November 2008.
4864
4865     # Remove MySQL tables from values in stream.
4866
4867     my ( $in,        # handle to in stream
4868          $out,       # handle to out stream
4869          $options,   # options hash
4870        ) = @_;
4871
4872     # Returns nothing.
4873
4874     my ( $record, %table_hash, $dbh, $table );
4875
4876     $options->{ "user" }     ||= Maasha::UCSC::ucsc_get_user();
4877     $options->{ "password" } ||= Maasha::UCSC::ucsc_get_password();
4878
4879     map { $table_hash{ $_ } = 1 } @{ $options->{ 'tables' } };
4880
4881     while ( $record = get_record( $in ) )
4882     {
4883         map { $table_hash{ $record->{ $_ } } = 1 } @{ $options->{ 'keys' } };
4884
4885         put_record( $record, $out ) if not $options->{ 'no_stream' };
4886     }
4887
4888     $dbh = Maasha::SQL::connect( $options->{ "database" }, $options->{ "user" }, $options->{ "password" } );
4889
4890     foreach $table ( sort keys %table_hash )
4891     {
4892         if ( Maasha::SQL::table_exists( $dbh, $table ) )
4893         {
4894             print STDERR qq(Removing table "$table" from database "$options->{ 'database' }" ... ) if $options->{ 'verbose' };
4895             Maasha::SQL::delete_table( $dbh, $table );
4896             print STDERR "done.\n" if $options->{ 'verbose' };
4897         }
4898         else
4899         {
4900             print STDERR qq(WARNING: table "$table" not found in database "$options->{ 'database' }\n");
4901         }
4902     }
4903
4904     Maasha::SQL::disconnect( $dbh );
4905 }
4906
4907
4908 sub script_remove_ucsc_tracks
4909 {
4910     # Martin A. Hansen, November 2008.
4911
4912     # Remove track from MySQL tables and config file.
4913
4914     my ( $in,        # handle to in stream
4915          $out,       # handle to out stream
4916          $options,   # options hash
4917        ) = @_;
4918
4919     # Returns nothing.
4920
4921     my ( $record, %track_hash, $fh_in, $fh_out, $track, @tracks, @new_tracks, $dbh );
4922
4923     $options->{ 'user' }        ||= Maasha::UCSC::ucsc_get_user();
4924     $options->{ 'password' }    ||= Maasha::UCSC::ucsc_get_password();
4925     $options->{ 'config_file' } ||= "$ENV{ 'HOME' }/ucsc/my_tracks.ra";
4926
4927     map { $track_hash{ $_ } = 1 } @{ $options->{ 'tracks' } };
4928
4929     while ( $record = get_record( $in ) )
4930     {
4931         map { $track_hash{ $record->{ $_ } } = 1 } @{ $options->{ 'keys' } };
4932
4933         put_record( $record, $out ) if not $options->{ 'no_stream' };
4934     }
4935
4936     $fh_in = Maasha::Filesys::file_read_open( $options->{ 'config_file' } );
4937
4938     while ( $track = Maasha::UCSC::ucsc_config_entry_get( $fh_in ) ) {
4939         push @tracks, $track;
4940     }
4941
4942     close $fh_in;
4943
4944     foreach $track ( @tracks )
4945     {
4946         if ( $track->{ 'database' } eq $options->{ 'database' } and exists $track_hash{ $track->{ 'track' } } ) {
4947             print STDERR qq(Removing track "$track->{ 'track' }" from config file.\n) if $options->{ 'verbose' };
4948         } else {
4949             push @new_tracks, $track;
4950         }
4951     }
4952
4953     rename "$options->{ 'config_file' }", "$options->{ 'config_file' }~";
4954
4955     $fh_out = Maasha::Filesys::file_write_open( $options->{ 'config_file' } );
4956
4957     map { Maasha::UCSC::ucsc_config_entry_put( $_, $fh_out ) } @new_tracks;
4958
4959     close $fh_out;
4960
4961     # ---- locate track in database ----
4962
4963     $dbh = Maasha::SQL::connect( $options->{ "database" }, $options->{ "user" }, $options->{ "password" } );
4964
4965     foreach $track ( sort keys %track_hash )
4966     {
4967         if ( Maasha::SQL::table_exists( $dbh, $track ) )
4968         {
4969             print STDERR qq(Removing table "$track" from database "$options->{ 'database' }" ... ) if $options->{ 'verbose' };
4970             Maasha::SQL::delete_table( $dbh, $track );
4971             print STDERR "done.\n" if $options->{ 'verbose' };
4972         }
4973         else
4974         {
4975             print STDERR qq(WARNING: table "$track" not found in database "$options->{ 'database' }\n");
4976         }
4977     }
4978
4979     Maasha::SQL::disconnect( $dbh );
4980
4981     Maasha::Common::run( "ucscMakeTracks.pl", "-b > /dev/null 2>&1" );
4982 }
4983
4984
4985 sub script_rename_keys
4986 {
4987     # Martin A. Hansen, August 2007.
4988
4989     # Rename keys in stream.
4990
4991     my ( $in,        # handle to in stream
4992          $out,       # handle to out stream
4993          $options,   # options hash
4994        ) = @_;
4995
4996     # Returns nothing.
4997
4998     my ( $record );
4999
5000     while ( $record = get_record( $in ) )
5001     {
5002         if ( exists $record->{ $options->{ "keys" }->[ 0 ] } )
5003         {
5004             $record->{ $options->{ "keys" }->[ 1 ] } = $record->{ $options->{ "keys" }->[ 0 ] };
5005
5006             delete $record->{ $options->{ "keys" }->[ 0 ] };
5007         }
5008
5009         put_record( $record, $out );
5010     }
5011 }
5012
5013
5014 sub script_uniq_vals
5015 {
5016     # Martin A. Hansen, August 2007.
5017
5018     # Find unique values in stream.
5019
5020     my ( $in,        # handle to in stream
5021          $out,       # handle to out stream
5022          $options,   # options hash
5023        ) = @_;
5024
5025     # Returns nothing.
5026
5027     my ( %hash, $record );
5028
5029     while ( $record = get_record( $in ) )
5030     {
5031         if ( $record->{ $options->{ "key" } } )
5032         {
5033             if ( not $hash{ $record->{ $options->{ "key" } } } and not $options->{ "invert" } )
5034             {
5035                 put_record( $record, $out );
5036
5037                 $hash{ $record->{ $options->{ "key" } } } = 1;
5038             }
5039             elsif ( $hash{ $record->{ $options->{ "key" } } } and $options->{ "invert" } )
5040             {
5041                 put_record( $record, $out );
5042             }
5043             else
5044             {
5045                 $hash{ $record->{ $options->{ "key" } } } = 1;
5046             }
5047         }
5048         else
5049         {
5050             put_record( $record, $out );
5051         }
5052     }
5053 }
5054
5055
5056 sub script_merge_vals
5057 {
5058     # Martin A. Hansen, August 2007.
5059
5060     # Rename keys in stream.
5061
5062     my ( $in,        # handle to in stream
5063          $out,       # handle to out stream
5064          $options,   # options hash
5065        ) = @_;
5066
5067     # Returns nothing.
5068
5069     my ( $record, @join, $i );
5070
5071     $options->{ "delimit" } ||= '_';
5072
5073     while ( $record = get_record( $in ) )
5074     {
5075         if ( exists $record->{ $options->{ "keys" }->[ 0 ] } )
5076         {
5077             @join = $record->{ $options->{ "keys" }->[ 0 ] };
5078
5079             for ( $i = 1; $i < @{ $options->{ "keys" } }; $i++ ) {
5080                 push @join, $record->{ $options->{ "keys" }->[ $i ] } if exists $record->{ $options->{ "keys" }->[ $i ] };
5081             }
5082
5083             $record->{ $options->{ "keys" }->[ 0 ] } = join $options->{ "delimit" }, @join;
5084         }
5085
5086         put_record( $record, $out );
5087     }
5088 }
5089
5090
5091 sub script_merge_records
5092 {
5093     # Martin A. Hansen, July 2008.
5094
5095     # Merges records in the stream based on identical values of two given keys.
5096
5097     my ( $in,        # handle to in stream
5098          $out,       # handle to out stream
5099          $options,   # options hash
5100        ) = @_;
5101
5102     # Returns nothing.
5103
5104     my ( $merge, $record, $file1, $file2, $fh1, $fh2, $key1, $key2, @keys1, @keys2, @vals1, @vals2,
5105          $num1, $num2, $num, $cmp, $i );
5106
5107     $merge = $options->{ "merge" } || "AandB";
5108
5109     $file1 = "$BP_TMP/merge_records1.tmp";
5110     $file2 = "$BP_TMP/merge_records2.tmp";
5111
5112     $fh1   = Maasha::Common::write_open( $file1 );
5113     $fh2   = Maasha::Common::write_open( $file2 );
5114
5115     $key1  = $options->{ "keys" }->[ 0 ];
5116     $key2  = $options->{ "keys" }->[ 1 ];
5117
5118     $num   = $key2 =~ s/n$//;
5119     $num1  = 0;
5120     $num2  = 0;
5121
5122     while ( $record = get_record( $in ) )
5123     {
5124         if ( exists $record->{ $key1 } )
5125         {
5126             @keys1 = $key1;
5127             @vals1 = $record->{ $key1 };
5128
5129             delete $record->{ $key1 };
5130
5131             map { push @keys1, $_; push @vals1, $record->{ $_ } } keys %{ $record };
5132
5133             print $fh1 join( "\t", @vals1 ), "\n";
5134
5135             $num1++;
5136         }
5137         elsif ( exists $record->{ $key2 } )
5138         {
5139             @keys2 = $key2;
5140             @vals2 = $record->{ $key2 };
5141
5142             delete $record->{ $key2 };
5143
5144             map { push @keys2, $_; push @vals2, $record->{ $_ } } keys %{ $record };
5145
5146             print $fh2 join( "\t", @vals2 ), "\n";
5147
5148             $num2++;
5149         }
5150     }
5151
5152     close $fh1;
5153     close $fh2;
5154
5155     if ( $num )
5156     {
5157         Maasha::Common::run( "sort", "-k 1,1n $file1 > $file1.sort" ) and rename "$file1.sort", $file1;
5158         Maasha::Common::run( "sort", "-k 1,1n $file2 > $file2.sort" ) and rename "$file2.sort", $file2;
5159     }
5160     else
5161     {
5162         Maasha::Common::run( "sort", "-k 1,1 $file1 > $file1.sort" ) and rename "$file1.sort", $file1;
5163         Maasha::Common::run( "sort", "-k 1,1 $file2 > $file2.sort" ) and rename "$file2.sort", $file2;
5164     }
5165
5166     $fh1 = Maasha::Common::read_open( $file1 );
5167     $fh2 = Maasha::Common::read_open( $file2 );
5168
5169     @vals1 = Maasha::Common::get_fields( $fh1 );
5170     @vals2 = Maasha::Common::get_fields( $fh2 );
5171
5172     while ( $num1 > 0 and $num2 > 0 )
5173     {
5174         undef $record;
5175
5176         if ( $num ) {
5177             $cmp = $vals1[ 0 ] <=> $vals2[ 0 ];
5178         } else {
5179             $cmp = $vals1[ 0 ] cmp $vals2[ 0 ];
5180         }
5181
5182         if ( $cmp < 0 )
5183         {
5184             if ( $merge =~ /^(AorB|AnotB)$/ )
5185             {
5186                 for ( $i = 0; $i < @keys1; $i++ ) {
5187                     $record->{ $keys1[ $i ] } = $vals1[ $i ];
5188                 }
5189
5190                 put_record( $record, $out );
5191             }
5192
5193             @vals1 = Maasha::Common::get_fields( $fh1 );
5194             $num1--;
5195         }
5196         elsif ( $cmp > 0 )
5197         {
5198             if ( $merge =~ /^(BorA|BnotA)$/ )
5199             {
5200                 for ( $i = 0; $i < @keys2; $i++ ) {
5201                     $record->{ $keys2[ $i ] } = $vals2[ $i ];
5202                 }
5203
5204                 put_record( $record, $out );
5205             }
5206
5207             @vals2 = Maasha::Common::get_fields( $fh2 );
5208             $num2--;
5209         }
5210         else
5211         {
5212             if ( $merge =~ /^(AandB|AorB|BorA)$/ )
5213             {
5214                 for ( $i = 0; $i < @keys1; $i++ ) {
5215                     $record->{ $keys1[ $i ] } = $vals1[ $i ];
5216                 }
5217
5218                 for ( $i = 1; $i < @keys2; $i++ ) {
5219                     $record->{ $keys2[ $i ] } = $vals2[ $i ];
5220                 }
5221
5222                 put_record( $record, $out );
5223             }
5224
5225             @vals1 = Maasha::Common::get_fields( $fh1 );
5226             @vals2 = Maasha::Common::get_fields( $fh2 );
5227             $num1--;
5228             $num2--;
5229         }
5230     }
5231
5232     close $fh1;
5233     close $fh2;
5234
5235     unlink $file1;
5236     unlink $file2;
5237
5238     if ( $num1 > 0 and $merge =~ /^(AorB|AnotB)$/ )
5239     {
5240         undef $record;
5241
5242         for ( $i = 0; $i < @keys1; $i++ ) {
5243             $record->{ $keys1[ $i ] } = $vals1[ $i ];
5244         }
5245
5246         put_record( $record, $out );
5247     }
5248
5249     if ( $num2 > 0 and $merge =~ /^(BorA|BnotA)$/ )
5250     {
5251         undef $record;
5252
5253         for ( $i = 0; $i < @keys2; $i++ ) {
5254             $record->{ $keys2[ $i ] } = $vals2[ $i ];
5255         }
5256
5257         put_record( $record, $out );
5258     }
5259 }
5260
5261
5262 sub script_grab
5263 {
5264     # Martin A. Hansen, August 2007.
5265
5266     # Grab for records in stream.
5267
5268     my ( $in,        # handle to in stream
5269          $out,       # handle to out stream
5270          $options,   # options hash
5271        ) = @_;
5272
5273     # Returns nothing.
5274
5275     my ( $keys, $vals_only, $keys_only, $invert, $patterns, $pattern, $regex, $record, $key, $op, $val, %lookup_hash, $found );
5276
5277     $keys      = $options->{ 'keys' };
5278     $vals_only = $options->{ 'vals_only' };
5279     $keys_only = $options->{ 'keys_only' };
5280     $invert    = $options->{ 'invert' };
5281
5282     if ( $options->{ 'patterns' } )
5283     {
5284         $patterns = [ split ",", $options->{ 'patterns' } ];
5285     }
5286     elsif ( -f $options->{ 'patterns_in' } )
5287     {
5288         $patterns = Maasha::Patscan::read_patterns( $options->{ 'patterns_in' } );
5289     }
5290     elsif ( $options->{ 'regex' } )
5291     {
5292         if ( $options->{ 'case_insensitive' } ) {
5293             $regex = qr/$options->{ 'regex' }/i;
5294         } else {
5295             $regex = qr/$options->{ 'regex' }/;
5296         }
5297     }
5298     elsif ( -f $options->{ 'exact_in' } )
5299     {
5300         $patterns = Maasha::Patscan::read_patterns( $options->{ 'exact_in' } );
5301
5302         map { $lookup_hash{ $_ } = 1 } @{ $patterns };
5303
5304         undef $patterns;
5305     }
5306     elsif ( $options->{ 'eval' } )
5307     {
5308         if ( $options->{ 'eval' } =~ /^([^><=! ]+)\s*(>=|<=|>|<|=|!=|eq|ne)\s*(.+)$/ )
5309         {
5310             $key = $1;
5311             $op  = $2;
5312             $val = $3;
5313         }
5314     }
5315
5316     while ( $record = get_record( $in ) )
5317     {
5318         $found = 0;
5319
5320         if ( %lookup_hash ) {
5321             $found = grab_lookup( \%lookup_hash, $record, $keys, $vals_only, $keys_only );
5322         } elsif ( $patterns ) {
5323             $found = grab_patterns( $patterns, $record, $keys, $vals_only, $keys_only );
5324         } elsif ( $regex ) {
5325             $found = grab_regex( $regex, $record, $keys, $vals_only, $keys_only );
5326         } elsif ( $op ) {
5327             $found = grab_eval( $key, $op, $val, $record );
5328         }
5329
5330         if ( $found and not $invert ) {
5331             put_record( $record, $out );
5332         } elsif ( not $found and $invert ) {
5333             put_record( $record, $out );
5334         }
5335     }
5336 }
5337
5338
5339 sub script_compute
5340 {
5341     # Martin A. Hansen, August 2007.
5342
5343     # Evaluate extression for records in stream.
5344
5345     my ( $in,        # handle to in stream
5346          $out,       # handle to out stream
5347          $options,   # options hash
5348        ) = @_;
5349
5350     # Returns nothing.
5351
5352     my ( $record, $eval_key, @keys, $eval_val );
5353
5354     while ( $record = get_record( $in ) )
5355     {
5356         if ( $options->{ "eval" } )
5357         {
5358             if ( $options->{ "eval" } =~ /^(\S+)\s*=\s*(.+)$/ )
5359             {
5360                 $eval_key = $1;
5361                 $eval_val = $2;
5362
5363                 if ( not @keys )
5364                 {
5365                     @keys = split /\s+|\+|-|\*|\/|\*\*/, $eval_val;
5366
5367                     @keys = grep { exists $record->{ $_ } } @keys;
5368                 }
5369
5370                 map { $eval_val =~ s/\Q$_\E/$record->{ $_ }/g } @keys;
5371
5372                 $record->{ $eval_key } = eval "$eval_val";
5373                 Maasha::Common::error( qq(eval "$eval_key = $eval_val" failed -> $@) ) if $@;
5374             }
5375             else
5376             {
5377                 Maasha::Common::error( qq(Bad compute expression: "$options->{ 'eval' }"\n) );
5378             }
5379         }
5380
5381         put_record( $record, $out );
5382     }
5383 }
5384
5385
5386 sub script_flip_tab
5387 {
5388     # Martin A. Hansen, June 2008.
5389
5390     # Flip a table.
5391
5392     my ( $in,        # handle to in stream
5393          $out,       # handle to out stream
5394          $options,   # options hash
5395        ) = @_;
5396
5397     # Returns nothing.
5398
5399     my ( $record, $key, $A, $B, @rows, @matrix, $row, $i );
5400
5401     while ( $record = get_record( $in ) )
5402     {
5403         undef @rows;
5404
5405         foreach $key ( sort { $A = $a; $B = $b; $A =~ s/^V(\d+)$/$1/; $B =~ s/^V(\d+)$/$1/; $A <=> $B } keys %{ $record } )
5406         {
5407             push @rows, $record->{ $key };
5408
5409         }
5410
5411         push @matrix, [ @rows ];
5412     }
5413
5414     undef $record;
5415
5416     @matrix = Maasha::Matrix::matrix_flip( \@matrix );
5417
5418     foreach $row ( @matrix )
5419     {
5420         for ( $i = 0; $i < @{ $row }; $i++ ) {
5421             $record->{ "V$i" } = $row->[ $i ];
5422         }
5423
5424         put_record( $record, $out );
5425     }
5426 }
5427
5428
5429 sub script_add_ident
5430 {
5431     # Martin A. Hansen, May 2008.
5432
5433     # Add a unique identifier to each record in stream.
5434
5435     my ( $in,        # handle to in stream
5436          $out,       # handle to out stream
5437          $options,   # options hash
5438        ) = @_;
5439
5440     # Returns nothing.
5441
5442     my ( $record, $key, $prefix, $i );
5443
5444     $key    = $options->{ "key" }    || "ID";
5445     $prefix = $options->{ "prefix" } || "ID";
5446
5447     $i = 0;
5448
5449     while ( $record = get_record( $in ) )
5450     {
5451         $record->{ $key } = sprintf( "$prefix%08d", $i );
5452
5453         put_record( $record, $out );
5454
5455         $i++;
5456     }
5457 }
5458
5459
5460 sub script_count_records
5461 {
5462     # Martin A. Hansen, August 2007.
5463
5464     # Count records in stream.
5465
5466     my ( $in,        # handle to in stream
5467          $out,       # handle to out stream
5468          $options,   # options hash
5469        ) = @_;
5470
5471     # Returns nothing.
5472
5473     my ( $record, $count, $result, $fh, $line );
5474
5475     $count = 0;
5476
5477     if ( $options->{ "no_stream" } )
5478     {
5479         while ( $line = <$in> )
5480         {
5481             chomp $line;
5482
5483             $count++ if $line eq "---";
5484         }
5485     }
5486     else
5487     {
5488         while ( $record = get_record( $in ) )
5489         {
5490             put_record( $record, $out );
5491
5492             $count++;
5493         }
5494     }
5495
5496     $result = { "RECORDS_COUNT" => $count };
5497
5498     $fh = write_stream( $options->{ "data_out" } );
5499
5500     put_record( $result, $fh );
5501
5502     close $fh;
5503 }
5504
5505
5506 sub script_random_records
5507 {
5508     # Martin A. Hansen, August 2007.
5509
5510     # Pick a number or random records from stream.
5511
5512     my ( $in,        # handle to in stream
5513          $out,       # handle to out stream
5514          $options,   # options hash
5515        ) = @_;
5516
5517     # Returns nothing.
5518
5519     my ( $record, $tmp_file, $fh_out, $fh_in, $count, $i, %rand_hash, $rand, $max );
5520
5521     $options->{ "num" } ||= 10;
5522
5523     $tmp_file = "$BP_TMP/random_records.tmp";
5524
5525     $fh_out = Maasha::Common::write_open( $tmp_file );
5526
5527     $count = 0;
5528
5529     while ( $record = get_record( $in ) )
5530     {
5531         put_record( $record, $fh_out );
5532
5533         $count++;
5534     }
5535
5536     close $fh_out;
5537
5538     $max = 0;
5539     $i   = 0;
5540
5541     Maasha::Common::error( qq(Requested random records > records in stream) ) if $options->{ "num" } > $count;
5542
5543     while ( $i < $options->{ "num" } )
5544     {
5545         $rand = int( rand( $count ) );
5546
5547         if ( not exists $rand_hash{ $rand } )
5548         {
5549             $rand_hash{ $rand } = 1;
5550
5551             $max = $rand if $rand > $max;
5552
5553             $i++;
5554         }
5555     }
5556
5557     $fh_in = Maasha::Common::read_open( $tmp_file );
5558
5559     $count = 0;
5560
5561     while ( $record = get_record( $fh_in ) )
5562     {
5563         put_record( $record, $out ) if exists $rand_hash{ $count };
5564
5565         last if $count == $max;
5566
5567         $count++;
5568     }
5569
5570     close $fh_in;
5571
5572     unlink $tmp_file;
5573 }
5574
5575
5576 sub script_sort_records
5577 {
5578     # Martin A. Hansen, August 2007.
5579
5580     # Sort to sort records according to keys.
5581
5582     my ( $in,        # handle to in stream
5583          $out,       # handle to out stream
5584          $options,   # options hash
5585        ) = @_;
5586
5587     # Returns nothing.
5588
5589     my ( @keys, $key, @sort_cmd, $sort_str, $sort_sub, @records, $record, $i );
5590
5591     foreach $key ( @{ $options->{ "keys" } } )
5592     {
5593         if ( $key =~ s/n$// ) {
5594             push @sort_cmd, qq(\$a->{ "$key" } <=> \$b->{ "$key" });
5595         } else {
5596             push @sort_cmd, qq(\$a->{ "$key" } cmp \$b->{ "$key" });
5597         }
5598     }
5599
5600     $sort_str = join " or ", @sort_cmd;
5601     $sort_sub = eval "sub { $sort_str }";   # NB security issue!
5602
5603     while ( $record = get_record( $in ) ) {
5604         push @records, $record;
5605     }
5606
5607     @records = sort $sort_sub @records;
5608
5609     if ( $options->{ "reverse" } )
5610     {
5611         for ( $i = scalar @records - 1; $i >= 0; $i-- ) {
5612             put_record( $records[ $i ], $out );
5613         }
5614     }
5615     else
5616     {
5617         for ( $i = 0; $i < scalar @records; $i++ ) {
5618             put_record( $records[ $i ], $out );
5619         }
5620     }
5621 }
5622
5623
5624 sub script_count_vals
5625 {
5626     # Martin A. Hansen, August 2007.
5627
5628     # Count records in stream.
5629
5630     my ( $in,        # handle to in stream
5631          $out,       # handle to out stream
5632          $options,   # options hash
5633        ) = @_;
5634
5635     # Returns nothing.
5636
5637     my ( $num, $record, %count_hash, @records, $tmp_file, $fh_out, $fh_in, $cache );
5638
5639     $tmp_file = "$BP_TMP/count_cache.tmp";
5640
5641     $fh_out   = Maasha::Common::write_open( $tmp_file );
5642
5643     $cache    = 0;
5644     $num      = 0;
5645
5646     while ( $record = get_record( $in ) )
5647     {
5648         map { $count_hash{ $_ }{ $record->{ $_ } }++ if exists $record->{ $_ } } @{ $options->{ "keys" } };
5649
5650         push @records, $record;
5651
5652         if ( scalar @records > 5_000_000 )   # too many records to hold in memory - use disk cache
5653         {
5654             map { put_record( $_, $fh_out ) } @records;
5655
5656             undef @records;
5657
5658             $cache = 1;
5659         }
5660
5661         print STDERR "verbose: records read $num\n" if ( $options->{ 'verbose' } and ( $num % 1_000_000 ) == 0 );
5662
5663         $num++;
5664     }
5665
5666     close $fh_out;
5667
5668     if ( $cache )
5669     {
5670         $num   = 0;
5671
5672         $fh_in = Maasha::Common::read_open( $tmp_file );
5673
5674         while ( $record = get_record( $fh_in ) )
5675         {
5676             map { $record->{ $_ . "_COUNT" } = $count_hash{ $_ }{ $record->{ $_ } } if exists $record->{ $_ } } @{ $options->{ "keys" } };
5677
5678             put_record( $record, $out );
5679
5680             print STDERR "verbose: cache read $num\n" if ( $options->{ 'verbose' } and ( $num % 1_000_000 ) == 0 );
5681
5682             $num++;
5683         }
5684
5685         close $fh_in;
5686     }
5687
5688     foreach $record ( @records )
5689     {
5690         map { $record->{ $_ . "_COUNT" } = $count_hash{ $_ }{ $record->{ $_ } } if exists $record->{ $_ } } @{ $options->{ "keys" } };
5691
5692         put_record( $record, $out );
5693     }
5694
5695     unlink $tmp_file;
5696 }
5697
5698
5699 sub script_plot_histogram
5700 {
5701     # Martin A. Hansen, September 2007.
5702
5703     # Plot a simple histogram for a given key using GNU plot.
5704
5705     my ( $in,        # handle to in stream
5706          $out,       # handle to out stream
5707          $options,   # options hash
5708        ) = @_;
5709
5710     # Returns nothing.
5711
5712     my ( $record, %data_hash, $max, @data_list, $i, $result, $fh );
5713
5714     $options->{ "title" } ||= "Histogram";
5715     $options->{ "sort" }  ||= "num";
5716
5717     while ( $record = get_record( $in ) )
5718     {
5719         $data_hash{ $record->{ $options->{ "key" } } }++ if defined $record->{ $options->{ "key" } };
5720
5721         put_record( $record, $out ) if not $options->{ "no_stream" };
5722     }
5723
5724     if ( $options->{ "sort" } eq "num" ) {
5725         map { push @data_list, [ $_, $data_hash{ $_ } ] } sort { $a <=> $b } keys %data_hash;
5726     } else {
5727         map { push @data_list, [ $_, $data_hash{ $_ } ] } sort keys %data_hash;
5728     }
5729
5730     $result = Maasha::Plot::histogram_simple( \@data_list, $options );
5731
5732     $fh = write_stream( $options->{ "data_out" } );
5733
5734     print $fh "$_\n" foreach @{ $result };
5735
5736     close $fh;
5737 }
5738
5739
5740 sub script_plot_lendist
5741 {
5742     # Martin A. Hansen, August 2007.
5743
5744     # Plot length distribution using GNU plot.
5745
5746     my ( $in,        # handle to in stream
5747          $out,       # handle to out stream
5748          $options,   # options hash
5749        ) = @_;
5750
5751     # Returns nothing.
5752
5753     my ( $record, %data_hash, $max, @data_list, $i, $result, $fh );
5754
5755     $options->{ "title" } ||= "Length Distribution";
5756
5757     while ( $record = get_record( $in ) )
5758     {
5759         $data_hash{ $record->{ $options->{ "key" } } }++ if defined $record->{ $options->{ "key" } };
5760
5761         put_record( $record, $out ) if not $options->{ "no_stream" };
5762     }
5763
5764     $max = Maasha::Calc::list_max( [ keys %data_hash ] );
5765
5766     for ( $i = 0; $i < $max; $i++ ) {
5767         push @data_list, [ $i, $data_hash{ $i } || 0 ];
5768     }
5769
5770     $result = Maasha::Plot::histogram_lendist( \@data_list, $options );
5771
5772     $fh = write_stream( $options->{ "data_out" } );
5773
5774     print $fh "$_\n" foreach @{ $result };
5775
5776     close $fh;
5777 }
5778
5779
5780 sub script_plot_chrdist
5781 {
5782     # Martin A. Hansen, August 2007.
5783
5784     # Plot chromosome distribution using GNU plot.
5785
5786     my ( $in,        # handle to in stream
5787          $out,       # handle to out stream
5788          $options,   # options hash
5789        ) = @_;
5790
5791     # Returns nothing.
5792
5793     my ( $record, %data_hash, @data_list, $elem, $sort_key, $count, $result, $fh );
5794
5795     $options->{ "title" } ||= "Chromosome Distribution";
5796
5797     while ( $record = get_record( $in ) )
5798     {
5799         if ( $record->{ "CHR" } ) {                                                             # generic
5800             $data_hash{ $record->{ "CHR" } }++;
5801         } elsif ( $record->{ "REC_TYPE" } eq "PATSCAN" and $record->{ "S_ID" } =~ /^chr/i ) {   # patscan
5802             $data_hash{ $record->{ "S_ID" } }++;
5803         } elsif ( $record->{ "REC_TYPE" } eq "PSL" and $record->{ "S_ID" } =~ /^chr/i ) {       # BLAT / PSL
5804             $data_hash{ $record->{ "S_ID" } }++;
5805         } elsif ( $record->{ "REC_TYPE" } eq "BLAST" and $record->{ "S_ID" } =~ /^chr/i ) {     # BLAST
5806             $data_hash{ $record->{ "S_ID" } }++;
5807         }
5808
5809         put_record( $record, $out ) if not $options->{ "no_stream" };
5810     }
5811
5812     foreach $elem ( keys %data_hash )
5813     {
5814         $sort_key = $elem;
5815
5816         $sort_key =~ s/chr//i;
5817
5818         $sort_key =~ s/^X(.*)/99$1/;
5819         $sort_key =~ s/^Y(.*)/99$1/;
5820         $sort_key =~ s/^Z(.*)/999$1/;
5821         $sort_key =~ s/^M(.*)/9999$1/;
5822         $sort_key =~ s/^U(.*)/99999$1/;
5823
5824         $count = $sort_key =~ tr/_//;
5825
5826         $sort_key =~ s/_.*/"999999" x $count/ex;
5827
5828         push @data_list, [ $elem, $data_hash{ $elem }, $sort_key ];
5829     }
5830
5831     @data_list = sort { $a->[ 2 ] <=> $b->[ 2 ] } @data_list;
5832
5833     $result = Maasha::Plot::histogram_chrdist( \@data_list, $options );
5834
5835     $fh = write_stream( $options->{ "data_out" } );
5836
5837     print $fh "$_\n" foreach @{ $result };
5838
5839     close $fh;
5840 }
5841
5842
5843 sub script_plot_karyogram
5844 {
5845     # Martin A. Hansen, August 2007.
5846
5847     # Plot hits on karyogram.
5848
5849     my ( $in,        # handle to in stream
5850          $out,       # handle to out stream
5851          $options,   # options hash
5852        ) = @_;
5853
5854     # Returns nothing.
5855
5856     my ( %options, $record, @data, $fh, $result, %data_hash );
5857
5858     $options->{ "genome" }     ||= "human";
5859     $options->{ "feat_color" } ||= "black";
5860
5861     while ( $record = get_record( $in ) )
5862     {
5863         if ( $record->{ "CHR" } and $record->{ "CHR_BEG" } and $record->{ "CHR_END" } )
5864         {
5865             push @{ $data_hash{ $record->{ "CHR" } } }, [ $record->{ "CHR_BEG" }, $record->{ "CHR_END" }, $options->{ "feat_color" } ];
5866         }
5867
5868         put_record( $record, $out ) if not $options->{ "no_stream" };
5869     }
5870
5871     $result = Maasha::Plot::karyogram( \%data_hash, \%options );
5872
5873     $fh = write_stream( $options->{ "data_out" } );
5874
5875     print $fh $result;
5876
5877     close $fh;
5878 }
5879
5880
5881 sub script_plot_matches
5882 {
5883     # Martin A. Hansen, August 2007.
5884
5885     # Plot matches in 2D generating a dotplot.
5886
5887     my ( $in,        # handle to in stream
5888          $out,       # handle to out stream
5889          $options,   # options hash
5890        ) = @_;
5891
5892     # Returns nothing.
5893
5894     my ( $record, @data, $fh, $result, %data_hash );
5895
5896     $options->{ "direction" } ||= "both";
5897
5898     while ( $record = get_record( $in ) )
5899     {
5900         if ( defined $record->{ "Q_BEG" } and defined $record->{ "S_BEG" } and $record->{ "Q_END" } and $record->{ "S_END" } ) {
5901             push @data, $record;
5902         }
5903
5904         put_record( $record, $out ) if not $options->{ "no_stream" };
5905     }
5906
5907     $options->{ "title" }  ||= "plot_matches";
5908     $options->{ "xlabel" } ||= $data[ 0 ]->{ "Q_ID" };
5909     $options->{ "ylabel" } ||= $data[ 0 ]->{ "S_ID" };
5910
5911     $result = Maasha::Plot::dotplot_matches( \@data, $options, $BP_TMP );
5912
5913     $fh = write_stream( $options->{ "data_out" } );
5914
5915     print $fh "$_\n" foreach @{ $result };
5916
5917     close $fh;
5918 }
5919
5920
5921 sub script_length_vals
5922 {
5923     # Martin A. Hansen, August 2007.
5924
5925     # Determine the length of the value for given keys.
5926
5927     my ( $in,        # handle to in stream
5928          $out,       # handle to out stream
5929          $options,   # options hash
5930        ) = @_;
5931
5932     # Returns nothing.
5933
5934     my ( $record, $key );
5935
5936     while ( $record = get_record( $in ) )
5937     {
5938         foreach $key ( @{ $options->{ "keys" } } )
5939         {
5940             if ( $record->{ $key } ) {
5941                 $record->{ $key . "_LEN" } = length $record->{ $key };
5942             }
5943         }
5944
5945         put_record( $record, $out );
5946     }
5947 }
5948
5949
5950 sub script_sum_vals
5951 {
5952     # Martin A. Hansen, August 2007.
5953
5954     # Calculates the sums for values of given keys.
5955
5956     my ( $in,        # handle to in stream
5957          $out,       # handle to out stream
5958          $options,   # options hash
5959        ) = @_;
5960
5961     # Returns nothing.
5962
5963     my ( $record, $key, %sum_hash, $fh );
5964
5965     while ( $record = get_record( $in ) )
5966     {
5967         foreach $key ( @{ $options->{ "keys" } } )
5968         {
5969             if ( $record->{ $key } ) {
5970                 $sum_hash{ $key } += $record->{ $key };
5971             }
5972         }
5973
5974         put_record( $record, $out ) if not $options->{ "no_stream" };
5975     }
5976
5977     $fh = write_stream( $options->{ "data_out" } );
5978
5979     foreach $key ( @{ $options->{ "keys" } } ) {
5980         put_record( { $key . "_SUM" => $sum_hash{ $key } || 0 } , $fh );
5981     }
5982
5983     close $fh;
5984 }
5985
5986
5987 sub script_mean_vals
5988 {
5989     # Martin A. Hansen, August 2007.
5990
5991     # Calculate the mean of values of given keys.
5992
5993     my ( $in,        # handle to in stream
5994          $out,       # handle to out stream
5995          $options,   # options hash
5996        ) = @_;
5997
5998     # Returns nothing.
5999
6000     my ( $record, $key, %sum_hash, %count_hash, $mean, $fh );
6001
6002     while ( $record = get_record( $in ) )
6003     {
6004         foreach $key ( @{ $options->{ "keys" } } )
6005         {
6006             if ( $record->{ $key } )
6007             {
6008                 $sum_hash{ $key } += $record->{ $key };
6009                 $count_hash{ $key }++;
6010             }
6011         }
6012
6013         put_record( $record, $out ) if not $options->{ "no_stream" };
6014     }
6015
6016     $fh = write_stream( $options->{ "data_out" } );
6017
6018     foreach $key ( @{ $options->{ "keys" } } )
6019     {
6020         if ( $count_hash{ $key } ) {
6021             $mean = sprintf( "%.2f", ( $sum_hash{ $key } / $count_hash{ $key } ) );
6022         } else {
6023             $mean = "N/A";
6024         }
6025
6026         put_record( { $key . "_MEAN" => $mean } , $fh );
6027     }
6028
6029     close $fh;
6030 }
6031
6032
6033 sub script_median_vals
6034 {
6035     # Martin A. Hansen, March 2008.
6036
6037     # Calculate the median values of given keys.
6038
6039     my ( $in,        # handle to in stream
6040          $out,       # handle to out stream
6041          $options,   # options hash
6042        ) = @_;
6043
6044     # Returns nothing.
6045
6046     my ( $record, $key, %median_hash, $median, $fh );
6047
6048     while ( $record = get_record( $in ) )
6049     {
6050         foreach $key ( @{ $options->{ "keys" } } ) {
6051             push @{ $median_hash{ $key } }, $record->{ $key } if defined $record->{ $key };
6052         }
6053
6054         put_record( $record, $out ) if not $options->{ "no_stream" };
6055     }
6056
6057     $fh = write_stream( $options->{ "data_out" } );
6058
6059     foreach $key ( @{ $options->{ "keys" } } )
6060     {
6061         if ( $median_hash{ $key } ) {
6062             $median = Maasha::Calc::median( $median_hash{ $key } );
6063         } else {
6064             $median = "N/A";
6065         }
6066
6067         put_record( { $key . "_MEDIAN" => $median } , $fh );
6068     }
6069
6070     close $fh;
6071 }
6072
6073
6074 sub script_max_vals
6075 {
6076     # Martin A. Hansen, February 2008.
6077
6078     # Determine the maximum values of given keys.
6079
6080     my ( $in,        # handle to in stream
6081          $out,       # handle to out stream
6082          $options,   # options hash
6083        ) = @_;
6084
6085     # Returns nothing.
6086
6087     my ( $record, $key, $fh, %max_hash, $max_record );
6088
6089     while ( $record = get_record( $in ) )
6090     {
6091         foreach $key ( @{ $options->{ "keys" } } )
6092         {
6093             if ( $record->{ $key } )
6094             {
6095                 $max_hash{ $key } = $record->{ $key } if $record->{ $key } > $max_hash{ $key };
6096             }
6097         }
6098
6099         put_record( $record, $out ) if not $options->{ "no_stream" };
6100     }
6101
6102     $fh = write_stream( $options->{ "data_out" } );
6103
6104     foreach $key ( @{ $options->{ "keys" } } )
6105     {
6106         $max_record->{ $key . "_MAX" } = $max_hash{ $key };
6107     }
6108
6109     put_record( $max_record, $fh );
6110
6111     close $fh;
6112 }
6113
6114
6115 sub script_min_vals
6116 {
6117     # Martin A. Hansen, February 2008.
6118
6119     # Determine the minimum values of given keys.
6120
6121     my ( $in,        # handle to in stream
6122          $out,       # handle to out stream
6123          $options,   # options hash
6124        ) = @_;
6125
6126     # Returns nothing.
6127
6128     my ( $record, $key, $fh, %min_hash, $min_record );
6129
6130     while ( $record = get_record( $in ) )
6131     {
6132         foreach $key ( @{ $options->{ "keys" } } )
6133         {
6134             if ( defined $record->{ $key } )
6135             {
6136                 if ( exists $min_hash{ $key } ) {
6137                     $min_hash{ $key } = $record->{ $key } if $record->{ $key } < $min_hash{ $key };
6138                 } else {
6139                     $min_hash{ $key } = $record->{ $key };
6140                 }
6141             }
6142         }
6143
6144         put_record( $record, $out ) if not $options->{ "no_stream" };
6145     }
6146
6147     $fh = write_stream( $options->{ "data_out" } );
6148
6149     foreach $key ( @{ $options->{ "keys" } } )
6150     {
6151         $min_record->{ $key . "_MIN" } = $min_hash{ $key };
6152     }
6153
6154     put_record( $min_record, $fh );
6155
6156     close $fh;
6157 }
6158
6159
6160 sub script_upload_to_ucsc
6161 {
6162     # Martin A. Hansen, August 2007.
6163
6164     # Calculate the mean of values of given keys.
6165
6166     # This routine has developed into the most ugly hack. Do something!
6167
6168     my ( $in,        # handle to in stream
6169          $out,       # handle to out stream
6170          $options,   # options hash
6171        ) = @_;
6172
6173     # Returns nothing.
6174
6175     my ( $record, $file, $wib_file, $wig_file, $wib_dir, $fh_out, $i, $first, $format, $type, $columns, $append, $entry );
6176
6177     $options->{ "short_label" } ||= $options->{ 'table' };
6178     $options->{ "long_label" }  ||= $options->{ 'table' };
6179     $options->{ "group" }       ||= $ENV{ "LOGNAME" };
6180     $options->{ "priority" }    ||= 1;
6181     $options->{ "visibility" }  ||= "pack";
6182     $options->{ "color" }       ||= join( ",", int( rand( 255 ) ), int( rand( 255 ) ), int( rand( 255 ) ) );
6183     $options->{ "chunk_size" }  ||= 10_000_000_000;    # Due to 32-bit UCSC compilation really large tables cannot be loaded in one go.
6184
6185     $file   = "$BP_TMP/ucsc_upload.tmp";
6186     $append = 0;
6187     $first  = 1;
6188     $i      = 0;
6189
6190     $fh_out = Maasha::Common::write_open( $file );
6191
6192     while ( $record = get_record( $in ) )
6193     {
6194         put_record( $record, $out ) if not $options->{ "no_stream" };
6195
6196         if ( $record->{ "REC_TYPE" } eq "fixed_step" )
6197         {
6198             $format = "WIGGLE";
6199
6200             if ( $entry = Maasha::UCSC::Wiggle::biopiece2fixedstep( $record ) ) {
6201                 Maasha::UCSC::Wiggle::fixedstep_entry_put( $entry, $fh_out );
6202             }
6203         }
6204         elsif ( $record->{ "REC_TYPE" } eq "PSL" )
6205         {
6206             $format = "PSL";
6207
6208             Maasha::UCSC::psl_put_header( $fh_out ) if $first;
6209             Maasha::UCSC::psl_put_entry( $record, $fh_out );
6210
6211             $first = 0;
6212         }
6213         elsif ( $record->{ "REC_TYPE" } eq "BED" and $record->{ "SEC_STRUCT" } )
6214         {
6215             # chrom chromStart  chromEnd    name    score   strand  size    secStr  conf
6216
6217             $format  = "BED_SS";
6218
6219             print $fh_out join ( "\t",
6220                 $record->{ "CHR" },
6221                 $record->{ "CHR_BEG" },
6222                 $record->{ "CHR_END" } + 1,
6223                 $record->{ "Q_ID" },
6224                 $record->{ "SCORE" },
6225                 $record->{ "STRAND" },
6226                 $record->{ "SIZE" },
6227                 $record->{ "SEC_STRUCT" },
6228                 $record->{ "CONF" },
6229             ), "\n";
6230         }
6231         elsif ( $record->{ "REC_TYPE" } eq "BED" )
6232         {
6233             $format  = "BED";
6234             $columns = $record->{ "BED_COLS" };
6235
6236             if ( $entry = Maasha::UCSC::BED::biopiece2bed( $record, $columns ) ) {
6237                 Maasha::UCSC::BED::bed_entry_put( $entry, $fh_out, $columns, $options->{ 'check' } );
6238             }
6239         }
6240         elsif ( $record->{ "REC_TYPE" } eq "PATSCAN" and $record->{ "CHR" } )
6241         {
6242             $format  = "BED";
6243             $columns = 6;
6244
6245             if ( $entry = Maasha::UCSC::BED::biopiece2bed( $record, $columns ) ) {
6246                 Maasha::UCSC::BED::bed_entry_put( $entry, $fh_out, $columns, $options->{ 'check' } );
6247             }
6248         }
6249         elsif ( $record->{ "REC_TYPE" } eq "BLAST" and $record->{ "S_ID" } =~ /^chr/ )
6250         {
6251             $format  = "BED";
6252             $columns = 6;
6253
6254             $record->{ "SCORE" } = $record->{ "BIT_SCORE" } * 1000;
6255
6256             if ( $entry = Maasha::UCSC::BED::biopiece2bed( $record, $columns ) ) {
6257                 Maasha::UCSC::BED::bed_entry_put( $entry, $fh_out, $columns, $options->{ 'check' } );
6258             }
6259         }
6260         elsif ( $record->{ "REC_TYPE" } eq "VMATCH" and $record->{ "S_ID" } =~ /^chr/i )
6261         {
6262             $format  = "BED";
6263             $columns = 6;
6264
6265             if ( $entry = Maasha::UCSC::BED::biopiece2bed( $record, $columns ) ) {
6266                 Maasha::UCSC::BED::bed_entry_put( $entry, $fh_out, $columns, $options->{ 'check' } );
6267             }
6268         }
6269
6270         if ( $i == $options->{ "chunk_size" } )
6271         {
6272             close $fh_out;
6273
6274             if ( $format eq "BED" ) {
6275                 Maasha::UCSC::bed_upload_to_ucsc( $BP_TMP, $file, $options, $append );
6276             } elsif ( $format eq "PSL" ) {
6277                 Maasha::UCSC::psl_upload_to_ucsc( $file, $options, $append );
6278             }
6279
6280             unlink $file;
6281
6282             $first = 1;
6283
6284             $append = 1;
6285
6286             $fh_out = Maasha::Common::write_open( $file );
6287         }
6288
6289         $i++;
6290     }
6291
6292     close $fh_out;
6293
6294     if ( exists $options->{ "database" } and $options->{ "table" } )
6295     {
6296         if ( $format eq "BED" )
6297         {
6298             $type = "bed $columns";
6299
6300             Maasha::UCSC::bed_upload_to_ucsc( $BP_TMP, $file, $options, $append );
6301         }
6302         elsif ( $format eq "BED_SS" )
6303         {
6304             $type = "type bed 6 +";
6305
6306             Maasha::UCSC::bed_upload_to_ucsc( $BP_TMP, $file, $options, $append );
6307         }
6308         elsif ( $format eq "PSL" )
6309         {
6310             $type = "psl";
6311
6312             Maasha::UCSC::psl_upload_to_ucsc( $file, $options, $append );
6313         }
6314         elsif ( $format eq "WIGGLE" )
6315         {
6316             $options->{ "visibility" } = "full";
6317
6318             $wig_file = "$options->{ 'table' }.wig";
6319             $wib_file = "$options->{ 'table' }.wib";
6320
6321             $wib_dir  = "$ENV{ 'HOME' }/ucsc/wib";
6322
6323             Maasha::Common::dir_create_if_not_exists( $wib_dir );
6324
6325             if ( $options->{ 'verbose' } ) {
6326                 `cd $BP_TMP && wigEncode $file $wig_file $wib_file`;
6327             } else {
6328                 `cd $BP_TMP && wigEncode $file $wig_file $wib_file > /dev/null 2>&1`;
6329             }
6330
6331             Maasha::Common::run( "mv", "$BP_TMP/$wib_file $wib_dir" );
6332
6333             unlink $file;
6334
6335             $file = $wig_file;
6336
6337             $type = "wig 0";
6338
6339             Maasha::UCSC::wiggle_upload_to_ucsc( $BP_TMP, $wib_dir, $file, $options );
6340         }
6341
6342         unlink $file;
6343
6344         Maasha::UCSC::ucsc_update_config( $options, $type );
6345     }
6346 }
6347
6348
6349 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
6350
6351
6352 sub grab_lookup
6353 {
6354     # Martin A. Hansen, November 2009.
6355
6356     # Uses keys from a lookup hash to search records. Optionally, a list of
6357     # keys can be given so the lookup is limited to these, also, flags
6358     # can be given to limit lookup to keys or vals only. Returns 1 if lookup
6359     # succeeded, else 0.
6360
6361     my ( $lookup_hash,   # hashref with patterns
6362          $record,        # hashref
6363          $keys,          # list of keys   - OPTIONAL
6364          $vals_only,     # only vals flag - OPTIONAL
6365          $keys_only,     # only keys flag - OPTIONAL
6366        ) = @_;
6367
6368     # Returns boolean.
6369
6370     if ( $keys )
6371     {
6372         map { return 1 if exists $lookup_hash->{ $record->{ $_ } } } @{ $keys };
6373     }
6374     else
6375     {
6376         if ( not $vals_only ) {
6377             map { return 1 if exists $lookup_hash->{ $_ } } keys %{ $record };
6378         }
6379
6380         if ( not $keys_only ) {
6381             map { return 1 if exists $lookup_hash->{ $record->{ $_ } } } keys %{ $record };
6382         }
6383     }
6384
6385     return 0;
6386 }
6387
6388
6389 sub grab_patterns
6390 {
6391     # Martin A. Hansen, November 2009.
6392
6393     # Uses patterns to match records containing the pattern as a substring.
6394     # Returns 1 if the record is matched, else 0.
6395
6396     my ( $patterns,   # list of patterns
6397          $record,     # hashref
6398          $keys,       # list of keys   - OPTIONAL
6399          $vals_only,  # only vals flag - OPTIONAL
6400          $keys_only,  # only keys flag - OPTIONAL
6401        ) = @_;
6402
6403     # Returns boolean.
6404
6405     my ( $pattern );
6406
6407     foreach $pattern ( @{ $patterns } )
6408     {
6409         if ( $keys )
6410         {
6411             map { return 1 if index( $record->{ $_ }, $pattern ) >= 0 } @{ $keys };
6412         }
6413         else
6414         {
6415             if ( not $vals_only ) {
6416                 map { return 1 if index( $_, $pattern ) >= 0 } keys %{ $record };
6417             }
6418
6419             if ( not $keys_only ) {
6420                 map { return 1 if index( $record->{ $_ }, $pattern ) >= 0 } keys %{ $record };
6421             }
6422         }
6423     }
6424
6425     return 0;
6426 }
6427
6428
6429 sub grab_regex
6430 {
6431     # Martin A. Hansen, November 2009.
6432
6433     # Uses regex to match records.
6434     # Returns 1 if the record is matched, else 0.
6435
6436     my ( $regex,      # regex to match
6437          $record,     # hashref
6438          $keys,       # list of keys   - OPTIONAL
6439          $vals_only,  # only vals flag - OPTIONAL
6440          $keys_only,  # only keys flag - OPTIONAL
6441        ) = @_;
6442
6443     # Returns boolean.
6444
6445     if ( $keys )
6446     {
6447         map { return 1 if $record->{ $_ } =~ /$regex/ } @{ $keys };
6448     }
6449     else
6450     {
6451         if ( not $vals_only ) {
6452             map { return 1 if $_ =~ /$regex/ } keys %{ $record };
6453         }
6454
6455         if ( not $keys_only ) {
6456             map { return 1 if $record->{ $_ } =~ /$regex/ } keys %{ $record };
6457         }
6458     }
6459
6460     return 0;
6461 }
6462
6463
6464 sub grab_eval
6465 {
6466     # Martin A. Hansen, November 2009.
6467
6468     # Test if the value of a given record key evaluates according
6469     # to a given operator. Returns 1 if eval is OK, else 0.
6470
6471     my ( $key,     # record key
6472          $op,      # operator
6473          $val,     # value
6474          $record,  # hashref
6475        ) = @_;
6476
6477     # Returns boolean.
6478
6479     if ( defined $record->{ $key } )
6480     {
6481         return 1 if ( $op eq "<" and $record->{ $key } < $val );
6482         return 1 if ( $op eq ">" and $record->{ $key } > $val );
6483         return 1 if ( $op eq ">=" and $record->{ $key } >= $val );
6484         return 1 if ( $op eq "<=" and $record->{ $key } <= $val );
6485         return 1 if ( $op eq "=" and $record->{ $key } == $val );
6486         return 1 if ( $op eq "!=" and $record->{ $key } != $val );
6487         return 1 if ( $op eq "eq" and $record->{ $key } eq $val );
6488         return 1 if ( $op eq "ne" and $record->{ $key } ne $val );
6489     }
6490
6491     return 0;
6492 }
6493
6494
6495 sub read_stream
6496 {
6497     # Martin A. Hansen, July 2007.
6498
6499     # Opens a stream to STDIN or a file,
6500
6501     my ( $path,   # path - OPTIONAL
6502        ) = @_;
6503
6504     # Returns filehandle.
6505
6506     my ( $fh );
6507
6508     if ( not -t STDIN ) {
6509         $fh = Maasha::Common::read_stdin();
6510     } elsif ( not $path ) {
6511 #        Maasha::Common::error( qq(no data stream) );
6512     } else {
6513         $fh = Maasha::Common::read_open( $path );
6514     }
6515
6516 #    $fh->autoflush(1) if $fh;  # Disable file buffer for debugging.
6517
6518     return $fh;
6519 }
6520
6521
6522 sub write_stream
6523 {
6524     # Martin A. Hansen, August 2007.
6525
6526     # Opens a stream to STDOUT or a file.
6527
6528     my ( $path,   # path          - OPTIONAL
6529          $gzip,   # compress data - OPTIONAL
6530        ) = @_;
6531
6532     # Returns filehandle.
6533
6534     my ( $fh );
6535
6536     if ( $path ) {
6537         $fh = Maasha::Common::write_open( $path, $gzip );
6538     } else {
6539         $fh = Maasha::Common::write_stdout();
6540     }
6541
6542     return $fh;
6543 }
6544
6545
6546 sub get_record
6547 {
6548     # Martin A. Hansen, July 2007.
6549
6550     # Reads one record at a time and converts that record
6551     # to a Perl data structure (a hash) which is returned.
6552
6553     my ( $fh,   # handle to stream
6554        ) = @_;
6555
6556     # Returns a hash.
6557
6558     my ( $block, @lines, $line, $key, $value, %record );
6559
6560     local $/ = "\n---\n";
6561
6562     $block = <$fh>;
6563
6564     chomp $block;
6565
6566     return if not defined $block;
6567
6568     @lines = split "\n", $block;
6569
6570     foreach $line ( @lines )
6571     {
6572         ( $key, $value ) = split ": ", $line, 2;
6573
6574         $record{ $key } = $value;
6575     }
6576
6577     return wantarray ? %record : \%record;
6578 }
6579
6580
6581 sub put_record
6582 {
6583     # Martin A. Hansen, July 2007.
6584
6585     # Given a Perl datastructure (a hash ref) emits this to STDOUT or a filehandle.
6586
6587     my ( $data,   # data structure
6588          $fh,     # file handle - OPTIONAL
6589        ) = @_;
6590
6591     # Returns nothing.
6592
6593     if ( scalar keys %{ $data } )
6594     {
6595         if ( $fh )
6596         {
6597             map { print $fh "$_: $data->{ $_ }\n" } keys %{ $data };
6598             print $fh "---\n";
6599         }
6600         else
6601         {
6602             map { print "$_: $data->{ $_ }\n" } keys %{ $data };
6603             print "---\n";
6604         }
6605     }
6606
6607     undef $data;
6608 }
6609
6610
6611 sub getopt_files
6612 {
6613     # Martin A. Hansen, November 2007.
6614
6615     # Extracts files from an explicit GetOpt::Long argument
6616     # allowing for the use of glob. E.g.
6617     # --data_in=test.fna
6618     # --data_in=test.fna,test2.fna
6619     # --data_in=*.fna
6620     # --data_in=test.fna,/dir/*.fna
6621
6622     my ( $option,   # option from GetOpt::Long
6623        ) = @_;
6624
6625     # Returns a list.
6626
6627     my ( $elem, @files );
6628
6629     foreach $elem ( split ",", $option )
6630     {
6631         if ( -f $elem ) {
6632             push @files, $elem;
6633         } elsif ( $elem =~ /\*/ ) {
6634             push @files, glob( $elem );
6635         }
6636     }
6637
6638     return wantarray ? @files : \@files;
6639 }
6640
6641
6642 sub sig_handler
6643 {
6644     # Martin A. Hansen, April 2008.
6645
6646     # Removes temporary directory and exits gracefully.
6647     # This subroutine is meant to be run always as the last
6648     # thing even if a script is dies or is interrupted
6649     # or killed.
6650
6651     my ( $sig,   # signal from the %SIG
6652        ) = @_;
6653
6654     # print STDERR "signal->$sig<-\n";
6655
6656     chomp $sig;
6657
6658     sleep 1;
6659
6660     if ( -d $BP_TMP )
6661     {
6662         if ( $sig =~ /MAASHA_ERROR/ ) {
6663             print STDERR "\nProgram '$script' had an error"                     . "  -  Please wait for temporary data to be removed\n";
6664         } elsif ( $sig eq "INT" ) {
6665             print STDERR "\nProgram '$script' interrupted (ctrl-c was pressed)" . "  -  Please wait for temporary data to be removed\n";
6666         } elsif ( $sig eq "TERM" ) {
6667             print STDERR "\nProgram '$script' terminated (someone used kill?)"  . "  -  Please wait for temporary data to be removed\n";
6668         } else {
6669             print STDERR "\nProgram '$script' died->$sig"                       . "  -  Please wait for temporary data to be removed\n";
6670         }
6671
6672         clean_tmp();
6673     }
6674
6675     exit( 0 );
6676 }
6677
6678
6679 sub clean_tmp
6680 {
6681     # Martin A. Hansen, July 2008.
6682
6683     # Cleans out any unused temporary files and directories in BP_TMP.
6684
6685     # Returns nothing.
6686
6687     my ( $tmpdir, @dirs, $curr_pid, $dir, $user, $sid, $pid );
6688
6689     $tmpdir = $ENV{ 'BP_TMP' } || Maasha::Common::error( 'No BP_TMP variable in environment.' );
6690
6691     $curr_pid = Maasha::Common::get_processid();
6692
6693     @dirs = Maasha::Common::ls_dirs( $tmpdir );
6694
6695     foreach $dir ( @dirs )
6696     {
6697         if ( $dir =~ /^$tmpdir\/(.+)_(\d+)_(\d+)_bp_tmp$/ )
6698         {
6699             $user = $1;
6700             $sid  = $2;
6701             $pid  = $3;
6702
6703             if ( $user eq Maasha::Common::get_user() )
6704             {
6705                 if ( not Maasha::Common::process_running( $pid ) )
6706                 {
6707                     # print STDERR "Removing stale dir: $dir\n";
6708                     Maasha::Common::dir_remove( $dir );
6709                 }
6710                 elsif ( $pid == $curr_pid )
6711                 {
6712                     # print STDERR "Removing current dir: $dir\n";
6713                     Maasha::Common::dir_remove( $dir );
6714                 }
6715             }
6716         }
6717     }
6718 }
6719
6720
6721 END
6722 {
6723     clean_tmp();
6724 }
6725
6726
6727 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
6728
6729 1;
6730
6731 __END__