code_perl/Maasha/BioRun.pm

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