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     $script = "print_usage" if ( -t STDIN and keys %{ $options } <= 1 or $options->{ 'help' } );
 114
 115     $in  = Maasha::Biopieces::read_stream( $options->{ "stream_in" } );
 116     $out = Maasha::Biopieces::write_stream( $options->{ "stream_out" } );
 117
 118     if    ( $script eq "print_usage" )              { script_print_usage(               $in, $out, $options ) }
 119     elsif ( $script eq "read_psl" )                 { script_read_psl(                  $in, $out, $options ) }
 120     elsif ( $script eq "read_bed" )                 { script_read_bed(                  $in, $out, $options ) }
 121     elsif ( $script eq "read_fixedstep" )           { script_read_fixedstep(            $in, $out, $options ) }
 122     elsif ( $script eq "read_blast_tab" )           { script_read_blast_tab(            $in, $out, $options ) }
 123     elsif ( $script eq "read_embl" )                { script_read_embl(                 $in, $out, $options ) }
 124     elsif ( $script eq "read_stockholm" )           { script_read_stockholm(            $in, $out, $options ) }
 125     elsif ( $script eq "read_phastcons" )           { script_read_phastcons(            $in, $out, $options ) }
 126     elsif ( $script eq "read_soft" )                { script_read_soft(                 $in, $out, $options ) }
 127     elsif ( $script eq "read_gff" )                 { script_read_gff(                  $in, $out, $options ) }
 128     elsif ( $script eq "read_2bit" )                { script_read_2bit(                 $in, $out, $options ) }
 129     elsif ( $script eq "read_solexa" )              { script_read_solexa(               $in, $out, $options ) }
 130     elsif ( $script eq "read_solid" )               { script_read_solid(                $in, $out, $options ) }
 131     elsif ( $script eq "read_mysql" )               { script_read_mysql(                $in, $out, $options ) }
 132     elsif ( $script eq "read_ucsc_config" )         { script_read_ucsc_config(          $in, $out, $options ) }
 133     elsif ( $script eq "assemble_tag_contigs" )     { script_assemble_tag_contigs(      $in, $out, $options ) }
 134     elsif ( $script eq "length_seq" )               { script_length_seq(                $in, $out, $options ) }
 135     elsif ( $script eq "uppercase_seq" )            { script_uppercase_seq(             $in, $out, $options ) }
 136     elsif ( $script eq "shuffle_seq" )              { script_shuffle_seq(               $in, $out, $options ) }
 137     elsif ( $script eq "analyze_seq" )              { script_analyze_seq(               $in, $out, $options ) }
 138     elsif ( $script eq "analyze_tags" )             { script_analyze_tags(              $in, $out, $options ) }
 139     elsif ( $script eq "complexity_seq" )           { script_complexity_seq(            $in, $out, $options ) }
 140     elsif ( $script eq "oligo_freq" )               { script_oligo_freq(                $in, $out, $options ) }
 141     elsif ( $script eq "create_weight_matrix" )     { script_create_weight_matrix(      $in, $out, $options ) }
 142     elsif ( $script eq "calc_bit_scores" )          { script_calc_bit_scores(           $in, $out, $options ) }
 143     elsif ( $script eq "calc_fixedstep" )           { script_calc_fixedstep(            $in, $out, $options ) }
 144     elsif ( $script eq "remove_indels" )            { script_remove_indels(             $in, $out, $options ) }
 145     elsif ( $script eq "transliterate_seq" )        { script_transliterate_seq(         $in, $out, $options ) }
 146     elsif ( $script eq "transliterate_vals" )       { script_transliterate_vals(        $in, $out, $options ) }
 147     elsif ( $script eq "translate_seq" )            { script_translate_seq(             $in, $out, $options ) }
 148     elsif ( $script eq "get_genome_align" )         { script_get_genome_align(          $in, $out, $options ) }
 149     elsif ( $script eq "get_genome_phastcons" )     { script_get_genome_phastcons(      $in, $out, $options ) }
 150     elsif ( $script eq "fold_seq" )                 { script_fold_seq(                  $in, $out, $options ) }
 151     elsif ( $script eq "split_seq" )                { script_split_seq(                 $in, $out, $options ) }
 152     elsif ( $script eq "split_bed" )                { script_split_bed(                 $in, $out, $options ) }
 153     elsif ( $script eq "align_seq" )                { script_align_seq(                 $in, $out, $options ) }
 154     elsif ( $script eq "tile_seq" )                 { script_tile_seq(                  $in, $out, $options ) }
 155     elsif ( $script eq "invert_align" )             { script_invert_align(              $in, $out, $options ) }
 156     elsif ( $script eq "patscan_seq" )              { script_patscan_seq(               $in, $out, $options ) }
 157     elsif ( $script eq "blat_seq" )                 { script_blat_seq(                  $in, $out, $options ) }
 158     elsif ( $script eq "soap_seq" )                 { script_soap_seq(                  $in, $out, $options ) }
 159     elsif ( $script eq "match_seq" )                { script_match_seq(                 $in, $out, $options ) }
 160     elsif ( $script eq "create_vmatch_index" )      { script_create_vmatch_index(       $in, $out, $options ) }
 161     elsif ( $script eq "write_align" )              { script_write_align(               $in, $out, $options ) }
 162     elsif ( $script eq "write_bed" )                { script_write_bed(                 $in, $out, $options ) }
 163     elsif ( $script eq "write_psl" )                { script_write_psl(                 $in, $out, $options ) }
 164     elsif ( $script eq "write_fixedstep" )          { script_write_fixedstep(           $in, $out, $options ) }
 165     elsif ( $script eq "write_2bit" )               { script_write_2bit(                $in, $out, $options ) }
 166     elsif ( $script eq "write_solid" )              { script_write_solid(               $in, $out, $options ) }
 167     elsif ( $script eq "write_ucsc_config" )        { script_write_ucsc_config(         $in, $out, $options ) }
 168     elsif ( $script eq "head_records" )             { script_head_records(              $in, $out, $options ) }
 169     elsif ( $script eq "remove_keys" )              { script_remove_keys(               $in, $out, $options ) }
 170     elsif ( $script eq "remove_adaptor" )           { script_remove_adaptor(            $in, $out, $options ) }
 171     elsif ( $script eq "remove_mysql_tables" )      { script_remove_mysql_tables(       $in, $out, $options ) }
 172     elsif ( $script eq "remove_ucsc_tracks" )       { script_remove_ucsc_tracks(        $in, $out, $options ) }
 173     elsif ( $script eq "rename_keys" )              { script_rename_keys(               $in, $out, $options ) }
 174     elsif ( $script eq "uniq_vals" )                { script_uniq_vals(                 $in, $out, $options ) }
 175     elsif ( $script eq "merge_vals" )               { script_merge_vals(                $in, $out, $options ) }
 176     elsif ( $script eq "merge_records" )            { script_merge_records(             $in, $out, $options ) }
 177     elsif ( $script eq "compute" )                  { script_compute(                   $in, $out, $options ) }
 178     elsif ( $script eq "flip_tab" )                 { script_flip_tab(                  $in, $out, $options ) }
 179     elsif ( $script eq "add_ident" )                { script_add_ident(                 $in, $out, $options ) }
 180     elsif ( $script eq "count_records" )            { script_count_records(             $in, $out, $options ) }
 181     elsif ( $script eq "random_records" )           { script_random_records(            $in, $out, $options ) }
 182     elsif ( $script eq "sort_records" )             { script_sort_records(              $in, $out, $options ) }
 183     elsif ( $script eq "count_vals" )               { script_count_vals(                $in, $out, $options ) }
 184     elsif ( $script eq "plot_histogram" )           { script_plot_histogram(            $in, $out, $options ) }
 185     elsif ( $script eq "plot_lendist" )             { script_plot_lendist(              $in, $out, $options ) }
 186     elsif ( $script eq "plot_chrdist" )             { script_plot_chrdist(              $in, $out, $options ) }
 187     elsif ( $script eq "plot_karyogram" )           { script_plot_karyogram(            $in, $out, $options ) }
 188     elsif ( $script eq "plot_matches" )             { script_plot_matches(              $in, $out, $options ) }
 189     elsif ( $script eq "plot_seqlogo" )             { script_plot_seqlogo(              $in, $out, $options ) }
 190     elsif ( $script eq "plot_phastcons_profiles" )  { script_plot_phastcons_profiles(   $in, $out, $options ) }
 191     elsif ( $script eq "analyze_bed" )              { script_analyze_bed(               $in, $out, $options ) }
 192     elsif ( $script eq "analyze_vals" )             { script_analyze_vals(              $in, $out, $options ) }
 193     elsif ( $script eq "length_vals" )              { script_length_vals(               $in, $out, $options ) }
 194     elsif ( $script eq "sum_vals" )                 { script_sum_vals(                  $in, $out, $options ) }
 195     elsif ( $script eq "mean_vals" )                { script_mean_vals(                 $in, $out, $options ) }
 196     elsif ( $script eq "median_vals" )              { script_median_vals(               $in, $out, $options ) }
 197     elsif ( $script eq "max_vals" )                 { script_max_vals(                  $in, $out, $options ) }
 198     elsif ( $script eq "min_vals" )                 { script_min_vals(                  $in, $out, $options ) }
 199     elsif ( $script eq "upload_to_ucsc" )           { script_upload_to_ucsc(            $in, $out, $options ) }
 200
 201     close $in if defined $in;
 202     close $out;
 203
 204     $t1 = gettimeofday();
 205
 206     print STDERR "Program: $script" . ( " " x ( 25 - length( $script ) ) ) . sprintf( "Run time: %.4f\n", ( $t1 - $t0 ) ) if $options->{ 'verbose' };
 207 }
 208
 209
 210 sub get_options
 211 {
 212     # Martin A. Hansen, February 2008.
 213
 214     # Gets options from commandline and checks these vigerously.
 215
 216     my ( $script,     # name of script
 217        ) = @_;
 218
 219     # Returns hash
 220
 221     my ( %options, @options, $opt, @genomes, $real );
 222
 223     if ( $script eq "print_usage" )
 224     {
 225         @options = qw(
 226             data_in|i=s
 227         );
 228     }
 229     elsif ( $script eq "read_psl" )
 230     {
 231         @options = qw(
 232             data_in|i=s
 233             num|n=s
 234         );
 235     }
 236     elsif ( $script eq "read_bed" )
 237     {
 238         @options = qw(
 239             data_in|i=s
 240             cols|c=s
 241             num|n=s
 242             check|C
 243         );
 244     }
 245     elsif ( $script eq "read_fixedstep" )
 246     {
 247         @options = qw(
 248             data_in|i=s
 249             num|n=s
 250         );
 251     }
 252     elsif ( $script eq "read_blast_tab" )
 253     {
 254         @options = qw(
 255             data_in|i=s
 256             num|n=s
 257         );
 258     }
 259     elsif ( $script eq "read_embl" )
 260     {
 261         @options = qw(
 262             data_in|i=s
 263             num|n=s
 264             keys|k=s
 265             feats|f=s
 266             quals|q=s
 267         );
 268     }
 269     elsif ( $script eq "read_stockholm" )
 270     {
 271         @options = qw(
 272             data_in|i=s
 273             num|n=s
 274         );
 275     }
 276     elsif ( $script eq "read_phastcons" )
 277     {
 278         @options = qw(
 279             data_in|i=s
 280             num|n=s
 281             min|m=s
 282             dist|d=s
 283             threshold|t=f
 284             gap|g=s
 285         );
 286     }
 287     elsif ( $script eq "read_soft" )
 288     {
 289         @options = qw(
 290             data_in|i=s
 291             samples|s=s
 292             num|n=s
 293         );
 294     }
 295     elsif ( $script eq "read_gff" )
 296     {
 297         @options = qw(
 298             data_in|i=s
 299             num|n=s
 300         );
 301     }
 302     elsif ( $script eq "read_2bit" )
 303     {
 304         @options = qw(
 305             data_in|i=s
 306             num|n=s
 307             no_mask|N
 308         );
 309     }
 310     elsif ( $script eq "read_solexa" )
 311     {
 312         @options = qw(
 313             data_in|i=s
 314             num|n=s
 315             format|f=s
 316             quality|q=s
 317         );
 318     }
 319     elsif ( $script eq "read_solid" )
 320     {
 321         @options = qw(
 322             data_in|i=s
 323             num|n=s
 324             quality|q=s
 325         );
 326     }
 327     elsif ( $script eq "read_mysql" )
 328     {
 329         @options = qw(
 330             database|d=s
 331             query|q=s
 332             user|u=s
 333             password|p=s
 334         );
 335     }
 336     elsif ( $script eq "read_ucsc_config" )
 337     {
 338         @options = qw(
 339             data_in|i=s
 340             num|n=s
 341         );
 342     }
 343     elsif ( $script eq "assemble_tag_contigs" )
 344     {
 345         @options = qw(
 346             check|C
 347         );
 348     }
 349     elsif ( $script eq "calc_fixedstep" )
 350     {
 351         @options = qw(
 352             check|C
 353         );
 354     }
 355     elsif ( $script eq "length_seq" )
 356     {
 357         @options = qw(
 358             no_stream|x
 359             data_out|o=s
 360         );
 361     }
 362     elsif ( $script eq "oligo_freq" )
 363     {
 364         @options = qw(
 365             word_size|w=s
 366             all|a
 367         );
 368     }
 369     elsif ( $script eq "create_weight_matrix" )
 370     {
 371         @options = qw(
 372             percent|p
 373         );
 374     }
 375     elsif ( $script eq "transliterate_seq" )
 376     {
 377         @options = qw(
 378             search|s=s
 379             replace|r=s
 380             delete|d=s
 381         );
 382     }
 383     elsif ( $script eq "transliterate_vals" )
 384     {
 385         @options = qw(
 386             keys|k=s
 387             search|s=s
 388             replace|r=s
 389             delete|d=s
 390         );
 391     }
 392     elsif ( $script eq "translate_seq" )
 393     {
 394         @options = qw(
 395             frames|f=s
 396         );
 397     }
 398     elsif ( $script eq "get_genome_align" )
 399     {
 400         @options = qw(
 401             genome|g=s
 402             chr|c=s
 403             beg|b=s
 404             end|e=s
 405             len|l=s
 406             strand|s=s
 407         );
 408     }
 409     elsif ( $script eq "get_genome_phastcons" )
 410     {
 411         @options = qw(
 412             genome|g=s
 413             chr|c=s
 414             beg|b=s
 415             end|e=s
 416             len|l=s
 417             flank|f=s
 418         );
 419     }
 420     elsif ( $script eq "split_seq" )
 421     {
 422         @options = qw(
 423             word_size|w=s
 424             uniq|u
 425         );
 426     }
 427     elsif ( $script eq "split_bed" )
 428     {
 429         @options = qw(
 430             window_size|w=s
 431             step_size|s=s
 432         );
 433     }
 434     elsif ( $script eq "tile_seq" )
 435     {
 436         @options = qw(
 437             identity|i=s
 438             supress_indels|s
 439         );
 440     }
 441     elsif ( $script eq "invert_align" )
 442     {
 443         @options = qw(
 444             soft|s
 445         );
 446     }
 447     elsif ( $script eq "patscan_seq" )
 448     {
 449         @options = qw(
 450             patterns|p=s
 451             patterns_in|P=s
 452             comp|c
 453             max_hits|h=s
 454             max_misses|m=s
 455             genome|g=s
 456         );
 457     }
 458     elsif ( $script eq "blat_seq" )
 459     {
 460         @options = qw(
 461             genome|g=s
 462             tile_size|t=s
 463             step_size|s=s
 464             min_identity|m=s
 465             min_score|M=s
 466             one_off|o=s
 467             ooc|c
 468         );
 469     }
 470     elsif ( $script eq "soap_seq" )
 471     {
 472         @options = qw(
 473             in_file|i=s
 474             genome|g=s
 475             seed_size|s=s
 476             mismatches|m=s
 477             gap_size|G=s
 478             cpus|c=s
 479         );
 480     }
 481     elsif ( $script eq "match_seq" )
 482     {
 483         @options = qw(
 484             word_size|w=s
 485             direction|d=s
 486         );
 487     }
 488     elsif ( $script eq "create_vmatch_index" )
 489     {
 490         @options = qw(
 491             index_name|i=s
 492             prefix_length|p=s
 493             no_stream|x
 494         );
 495     }
 496     elsif ( $script eq "write_align" )
 497     {
 498         @options = qw(
 499             wrap|w=s
 500             no_stream|x
 501             no_ruler|R
 502             no_consensus|C
 503             data_out|o=s
 504         );
 505     }
 506     elsif ( $script eq "write_bed" )
 507     {
 508         @options = qw(
 509             cols|c=s
 510             check|C
 511             no_stream|x
 512             data_out|o=s
 513             compress|Z
 514         );
 515     }
 516     elsif ( $script eq "write_psl" )
 517     {
 518         @options = qw(
 519             no_stream|x
 520             data_out|o=s
 521             compress|Z
 522         );
 523     }
 524     elsif ( $script eq "write_fixedstep" )
 525     {
 526         @options = qw(
 527             no_stream|x
 528             data_out|o=s
 529             compress|Z
 530         );
 531     }
 532     elsif ( $script eq "write_2bit" )
 533     {
 534         @options = qw(
 535             no_stream|x
 536             data_out|o=s
 537             no_mask|N
 538         );
 539     }
 540     elsif ( $script eq "write_solid" )
 541     {
 542         @options = qw(
 543             wrap|w=s
 544             no_stream|x
 545             data_out|o=s
 546             compress|Z
 547         );
 548     }
 549     elsif ( $script eq "write_ucsc_config" )
 550     {
 551         @options = qw(
 552             no_stream|x
 553             data_out|o=s
 554         );
 555     }
 556     elsif ( $script eq "plot_seqlogo" )
 557     {
 558         @options = qw(
 559             no_stream|x
 560             data_out|o=s
 561         );
 562     }
 563     elsif ( $script eq "plot_phastcons_profiles" )
 564     {
 565         @options = qw(
 566             no_stream|x
 567             data_out|o=s
 568             genome|g=s
 569             mean|m
 570             median|M
 571             flank|f=s
 572             terminal|t=s
 573             title|T=s
 574             xlabel|X=s
 575             ylabel|Y=s
 576         );
 577     }
 578     elsif ( $script eq "analyze_vals" )
 579     {
 580         @options = qw(
 581             no_stream|x
 582             keys|k=s
 583         );
 584     }
 585     elsif ( $script eq "head_records" )
 586     {
 587         @options = qw(
 588             num|n=s
 589         );
 590     }
 591     elsif ( $script eq "remove_keys" )
 592     {
 593         @options = qw(
 594             keys|k=s
 595             save_keys|K=s
 596         );
 597     }
 598     elsif ( $script eq "remove_adaptor" )
 599     {
 600         @options = qw(
 601             adaptor|a=s
 602             mismatches|m=s
 603             remove|r=s
 604             offset|o=s
 605         );
 606     }
 607     elsif ( $script eq "remove_mysql_tables" )
 608     {
 609         @options = qw(
 610             database|d=s
 611             tables|t=s
 612             keys|k=s
 613             user|u=s
 614             password|p=s
 615             no_stream|x
 616         );
 617     }
 618     elsif ( $script eq "remove_ucsc_tracks" )
 619     {
 620         @options = qw(
 621             database|d=s
 622             tracks|t=s
 623             keys|k=s
 624             config_file|c=s
 625             user|u=s
 626             password|p=s
 627             no_stream|x
 628         );
 629     }
 630     elsif ( $script eq "rename_keys" )
 631     {
 632         @options = qw(
 633             keys|k=s
 634         );
 635     }
 636     elsif ( $script eq "uniq_vals" )
 637     {
 638         @options = qw(
 639             key|k=s
 640             invert|i
 641         );
 642     }
 643     elsif ( $script eq "merge_vals" )
 644     {
 645         @options = qw(
 646             keys|k=s
 647             delimit|d=s
 648         );
 649     }
 650     elsif ( $script eq "merge_records" )
 651     {
 652         @options = qw(
 653             keys|k=s
 654             merge|m=s
 655         );
 656     }
 657     elsif ( $script eq "compute" )
 658     {
 659         @options = qw(
 660             eval|e=s
 661         );
 662     }
 663     elsif ( $script eq "add_ident" )
 664     {
 665         @options = qw(
 666             prefix|p=s
 667             key|k=s
 668         );
 669     }
 670     elsif ( $script eq "count_records" )
 671     {
 672         @options = qw(
 673             no_stream|x
 674             data_out|o=s
 675         );
 676     }
 677     elsif ( $script eq "random_records" )
 678     {
 679         @options = qw(
 680             num|n=s
 681         );
 682     }
 683     elsif ( $script eq "sort_records" )
 684     {
 685         @options = qw(
 686             reverse|r
 687             keys|k=s
 688         );
 689     }
 690     elsif ( $script eq "count_vals" )
 691     {
 692         @options = qw(
 693             keys|k=s
 694         );
 695     }
 696     elsif ( $script eq "plot_histogram" )
 697     {
 698         @options = qw(
 699             no_stream|x
 700             data_out|o=s
 701             terminal|t=s
 702             title|T=s
 703             xlabel|X=s
 704             ylabel|Y=s
 705             key|k=s
 706             sort|s=s
 707         );
 708     }
 709     elsif ( $script eq "plot_lendist" )
 710     {
 711         @options = qw(
 712             no_stream|x
 713             data_out|o=s
 714             terminal|t=s
 715             title|T=s
 716             xlabel|X=s
 717             ylabel|Y=s
 718             key|k=s
 719         );
 720     }
 721     elsif ( $script eq "plot_chrdist" )
 722     {
 723         @options = qw(
 724             no_stream|x
 725             data_out|o=s
 726             terminal|t=s
 727             title|T=s
 728             xlabel|X=s
 729             ylabel|Y=s
 730         );
 731     }
 732     elsif ( $script eq "plot_karyogram" )
 733     {
 734         @options = qw(
 735             no_stream|x
 736             data_out|o=s
 737             genome|g=s
 738             feat_color|f=s
 739         );
 740     }
 741     elsif ( $script eq "plot_matches" )
 742     {
 743         @options = qw(
 744             no_stream|x
 745             data_out|o=s
 746             terminal|t=s
 747             title|T=s
 748             xlabel|X=s
 749             ylabel|Y=s
 750             direction|d=s
 751         );
 752     }
 753     elsif ( $script eq "length_vals" )
 754     {
 755         @options = qw(
 756             keys|k=s
 757         );
 758     }
 759     elsif ( $script eq "sum_vals" )
 760     {
 761         @options = qw(
 762             no_stream|x
 763             data_out|o=s
 764             keys|k=s
 765         );
 766     }
 767     elsif ( $script eq "mean_vals" )
 768     {
 769         @options = qw(
 770             no_stream|x
 771             data_out|o=s
 772             keys|k=s
 773         );
 774     }
 775     elsif ( $script eq "median_vals" )
 776     {
 777         @options = qw(
 778             no_stream|x
 779             data_out|o=s
 780             keys|k=s
 781         );
 782     }
 783     elsif ( $script eq "max_vals" )
 784     {
 785         @options = qw(
 786             no_stream|x
 787             data_out|o=s
 788             keys|k=s
 789         );
 790     }
 791     elsif ( $script eq "min_vals" )
 792     {
 793         @options = qw(
 794             no_stream|x
 795             data_out|o=s
 796             keys|k=s
 797         );
 798     }
 799     elsif ( $script eq "upload_to_ucsc" )
 800     {
 801         @options = qw(
 802             no_stream|x
 803             database|d=s
 804             table|t=s
 805             short_label|s=s
 806             long_label|l=s
 807             group|g=s
 808             priority|p=f
 809             use_score|u
 810             visibility|V=s
 811             color|c=s
 812             check|C
 813         );
 814     }
 815
 816     push @options, qw(
 817         stream_in|I=s
 818         stream_out|O=s
 819         verbose|v
 820         help|?
 821     );
 822
 823 #    print STDERR Dumper( \@options );
 824
 825     GetOptions(
 826         \%options,
 827         @options,
 828     );
 829
 830 #    print STDERR Dumper( \%options );
 831
 832     if ( -t STDIN && scalar( keys %options ) == 0 or $options{ "help" } ) {
 833         return wantarray ? %options : \%options;
 834     }
 835
 836     $options{ "cols" }      = [ split ",", $options{ "cols" } ]      if defined $options{ "cols" };
 837     $options{ "keys" }      = [ split ",", $options{ "keys" } ]      if defined $options{ "keys" };
 838     $options{ "no_keys" }   = [ split ",", $options{ "no_keys" } ]   if defined $options{ "no_keys" };
 839     $options{ "save_keys" } = [ split ",", $options{ "save_keys" } ] if defined $options{ "save_keys" };
 840     $options{ "quals" }     = [ split ",", $options{ "quals" } ]     if defined $options{ "quals" };
 841     $options{ "feats" }     = [ split ",", $options{ "feats" } ]     if defined $options{ "feats" };
 842     $options{ "frames" }    = [ split ",", $options{ "frames" } ]    if defined $options{ "frames" };
 843     $options{ "samples" }   = [ split ",", $options{ "samples" } ]   if defined $options{ "samples" };
 844     $options{ "tables" }    = [ split ",", $options{ "tables" } ]    if defined $options{ "tables" };
 845     $options{ "tracks" }    = [ split ",", $options{ "tracks" } ]    if defined $options{ "tracks" };
 846
 847     # ---- check arguments ----
 848
 849     if ( $options{ 'data_in' } )
 850     {
 851         $options{ "files" } = Maasha::Biopieces::getopt_files( $options{ 'data_in' } );
 852
 853         Maasha::Common::error( qq(Argument to --data_in must be a valid file or fileglob expression) ) if scalar @{ $options{ "files" } } == 0;
 854     }
 855
 856     map { Maasha::Common::error( qq(Argument to --cols must be a whole numbers - not "$_") ) if $_ !~ /^\d+$/ } @{ $options{ "cols" } } if $options{ "cols" };
 857
 858     # print STDERR Dumper( \%options );
 859
 860     $real = "beg|end|word_size|wrap|tile_size|len|prefix_length|mismatches|offset|num|skip|cpus|window_size|step_size";
 861
 862     foreach $opt ( keys %options )
 863     {
 864         if ( $opt =~ /stream_in|pattern_in|exact_in/ and not -f $options{ $opt } )
 865         {
 866             Maasha::Common::error( qq(Argument to --$opt must be a valid file or fileglob expression - not "$options{ $opt }") );
 867         }
 868         elsif ( $opt =~ /$real/ and $options{ $opt } !~ /^\d+$/ )
 869         {
 870             Maasha::Common::error( qq(Argument to --$opt must be a whole number - not "$options{ $opt }") );
 871         }
 872         elsif ( $opt =~ /max_hits|max_hits|max_misses|dist|edit_dist|flank|gap|hamming_dist|priority/ and $options{ $opt } !~ /^-?\d+$/ )
 873         {
 874             Maasha::Common::error( qq(Argument to --$opt must be an integer - not "$options{ $opt }") );
 875         }
 876         elsif ( $opt =~ /identity|threshold/ and $options{ $opt } !~ /^-?(?:\d+(?:\.\d*)?|\.\d+)$/ )
 877         {
 878             Maasha::Common::error( qq(Argument to --$opt must be a decimal number - not "$options{ $opt }") );
 879         }
 880         elsif ( $opt =~ /e_val/ and $options{ $opt } !~ /^([+-]?)(?=\d|\.\d)\d*(\.\d*)?([Ee]([+-]?\d+))?$/ )
 881         {
 882             Maasha::Common::error( qq(Argument to --$opt must be a float - not "$options{ $opt }") );
 883         }
 884         elsif ( $opt =~ /strand/ and $options{ $opt } !~ /^(\+|-)$/ )
 885         {
 886             Maasha::Common::error( qq(Argument to --$opt must be "+" or "-" - not "$options{ $opt }") );
 887         }
 888         elsif ( $opt eq "genome" )
 889         {
 890             @genomes = Maasha::Common::ls_dirs( "$ENV{ 'BP_DATA' }/genomes" );
 891             map { $_ =~ s/.*\/(.+)$/$1/ } @genomes;
 892
 893             if ( not grep { $_ =~ /^$options{ $opt }$/ } @genomes ) {
 894                 Maasha::Common::error( qq(Genome $options{ $opt } not found in "$ENV{ 'BP_DATA' }/genomes/") );
 895             }
 896         }
 897         elsif ( $opt eq "terminal" and not $options{ $opt } =~ /^(svg|post|dumb|x11)/ )
 898         {
 899             Maasha::Common::error( qq(Bad --$opt argument "$options{ $opt }") );
 900         }
 901         elsif ( $opt eq "table" and $options{ $opt } =~ /(-|\.)/ )
 902         {
 903             Maasha::Common::error( qq(Character '$1' is not allowed in table name: $options{ $opt }) );
 904         }
 905         elsif ( $opt eq "merge" and $options{ $opt } !~ /^(AandB|AorB|BorA|AnotB|BnotA)$/ )
 906         {
 907             Maasha::Common::error( qq(Argument to --$opt must be AandB, AorB, BorA, AnotB, or BnotA - not "$options{ $opt }") );
 908         }
 909         elsif ( $opt eq "format" and $script eq "read_solexa" and $options{ $opt } !~ /octal|decimal/ )
 910         {
 911             Maasha::Common::error( qq(Argument to --$opt must be octal or decimal - not "$options{ $opt }") );
 912         }
 913         elsif ( $opt eq "remove" and $script eq "remove_adaptor" and $options{ $opt } !~ /before|after|skip/ )
 914         {
 915             Maasha::Common::error( qq(Argument to --$opt must be before, after, or skip - not "$options{ $opt }") );
 916         }
 917     }
 918
 919     Maasha::Common::error( qq(no --database specified) )                if $script eq "remove_ucsc_tracks"  and not $options{ "database" };
 920     Maasha::Common::error( qq(no --index_name specified) )              if $script =~ /create_vmatch_index/ and not $options{ "index_name" };
 921     Maasha::Common::error( qq(no --in_file or --genome specified) )     if $script eq "soap_seq" and not $options{ "genome" } and not $options{ "in_file" };
 922     Maasha::Common::error( qq(both --in_file and --genome specified) )  if $script eq "soap_seq" and $options{ "genome" } and $options{ "in_file" };
 923     Maasha::Common::error( qq(no --genome specified) )                  if $script =~ /get_genome_align|get_genome_phastcons|blat_seq|plot_phastcons_profiles|plot_karyogram/ and not $options{ "genome" };
 924     Maasha::Common::error( qq(no --key specified) )                     if $script =~ /plot_lendist|plot_histogram/ and not $options{ "key" };
 925     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" };
 926
 927     if ( $script eq "upload_to_ucsc" )
 928     {
 929         Maasha::Common::error( qq(no --database specified) ) if not $options{ "database" };
 930         Maasha::Common::error( qq(no --table specified) )    if not $options{ "table" };
 931     }
 932
 933     return wantarray ? %options : \%options;
 934 }
 935
 936
 937 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> SCRIPTS  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 938
 939
 940 sub script_print_usage
 941 {
 942     # Martin A. Hansen, January 2008.
 943
 944     # Retrieves usage information from file and
 945     # prints this nicely formatted.
 946
 947     my ( $in,        # handle to in stream
 948          $out,       # handle to out stream
 949          $options,   # options hash
 950        ) = @_;
 951
 952     # Returns nothing.
 953
 954     my ( $file, $wiki, $lines );
 955
 956     if ( $options->{ 'data_in' } ) {
 957         $file = $options->{ 'data_in' };
 958     } else {
 959         $file = join "", $ENV{ 'BP_DIR' }, "/bp_usage/", $options->{ 'SCRIPT' }, ".wiki";
 960     }
 961
 962     $wiki = Maasha::Gwiki::gwiki_read( $file );
 963
 964     ( $wiki->[ 2 ], $wiki->[ 3 ], $wiki->[ 0 ], $wiki->[ 1 ] ) = ( $wiki->[ 0 ], $wiki->[ 1 ], $wiki->[ 2 ], $wiki->[ 3 ] );
 965
 966     if ( not $options->{ "help" } ) {
 967         @{ $wiki } = grep { $_->[ 0 ]->{ 'SECTION' } =~ /Biopiece|summary|Usage|Options|Help/ } @{ $wiki };
 968     }
 969
 970     $lines = Maasha::Gwiki::gwiki2ascii( $wiki );
 971
 972     print STDERR "$_\n" foreach @{ $lines };
 973
 974     exit;
 975 }
 976
 977
 978 sub script_read_psl
 979 {
 980     # Martin A. Hansen, August 2007.
 981
 982     # Read psl table from stream or file.
 983
 984     my ( $in,        # handle to in stream
 985          $out,       # handle to out stream
 986          $options,   # options hash
 987        ) = @_;
 988
 989     # Returns nothing.
 990
 991     my ( $record, $file, $data_in, $num );
 992
 993     while ( $record = Maasha::Biopieces::get_record( $in ) ) {
 994         Maasha::Biopieces::put_record( $record, $out );
 995     }
 996
 997     $num = 1;
 998
 999     foreach $file ( @{ $options->{ "files" } } )
1000     {
1001         $data_in = Maasha::Common::read_open( $file );
1002
1003         while ( $record = Maasha::UCSC::psl_get_entry( $data_in ) )
1004         {
1005             Maasha::Biopieces::put_record( $record, $out );
1006
1007             goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1008
1009             $num++;
1010         }
1011     }
1012
1013     NUM:
1014 }
1015
1016
1017 sub script_read_bed
1018 {
1019     # Martin A. Hansen, August 2007.
1020
1021     # Read bed table from stream or file.
1022
1023     my ( $in,        # handle to in stream
1024          $out,       # handle to out stream
1025          $options,   # options hash
1026        ) = @_;
1027
1028     # Returns nothing.
1029
1030     my ( $cols, $file, $record, $bed_entry, $data_in, $num );
1031
1032     $cols = $options->{ 'cols' }->[ 0 ];
1033
1034     while ( $record = Maasha::Biopieces::get_record( $in ) ) {
1035         Maasha::Biopieces::put_record( $record, $out );
1036     }
1037
1038     $num = 1;
1039
1040     foreach $file ( @{ $options->{ "files" } } )
1041     {
1042         $data_in = Maasha::Common::read_open( $file );
1043
1044         while ( $bed_entry = Maasha::UCSC::BED::bed_entry_get( $data_in, $cols, $options->{ 'check' } ) )
1045         {
1046             $record = Maasha::UCSC::BED::bed2biopiece( $bed_entry );
1047
1048             Maasha::Biopieces::put_record( $record, $out );
1049
1050             goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1051
1052             $num++;
1053         }
1054
1055         close $data_in;
1056     }
1057
1058     NUM:
1059
1060     close $data_in if $data_in;
1061 }
1062
1063
1064 sub script_read_fixedstep
1065 {
1066     # Martin A. Hansen, Juli 2008.
1067
1068     # Read fixedstep wiggle format from stream or file.
1069
1070     my ( $in,        # handle to in stream
1071          $out,       # handle to out stream
1072          $options,   # options hash
1073        ) = @_;
1074
1075     # Returns nothing.
1076
1077     my ( $file, $record, $entry, $data_in, $num );
1078
1079     while ( $record = Maasha::Biopieces::get_record( $in ) ) {
1080         Maasha::Biopieces::put_record( $record, $out );
1081     }
1082
1083     $num = 1;
1084
1085     foreach $file ( @{ $options->{ "files" } } )
1086     {
1087         $data_in = Maasha::Common::read_open( $file );
1088
1089         while ( $entry = Maasha::UCSC::Wiggle::fixedstep_entry_get( $data_in ) )
1090         {
1091             $record = Maasha::UCSC::Wiggle::fixedstep2biopiece( $entry );
1092
1093             Maasha::Biopieces::put_record( $record, $out );
1094
1095             goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1096
1097             $num++;
1098         }
1099
1100         close $data_in;
1101     }
1102
1103     NUM:
1104
1105     close $data_in if $data_in;
1106 }
1107
1108
1109 sub script_read_blast_tab
1110 {
1111     # Martin A. Hansen, September 2007.
1112
1113     # Read tabular BLAST output from NCBI blast run with -m8 or -m9.
1114
1115     my ( $in,        # handle to in stream
1116          $out,       # handle to out stream
1117          $options,   # options hash
1118        ) = @_;
1119
1120     # Returns nothing.
1121
1122     my ( $file, $line, @fields, $strand, $record, $data_in, $num );
1123
1124     while ( $record = Maasha::Biopieces::get_record( $in ) ) {
1125         Maasha::Biopieces::put_record( $record, $out );
1126     }
1127
1128     $num = 1;
1129
1130     foreach $file ( @{ $options->{ "files" } } )
1131     {
1132         $data_in = Maasha::Common::read_open( $file );
1133
1134         while ( $line = <$data_in> )
1135         {
1136             chomp $line;
1137
1138             next if $line =~ /^#/;
1139
1140             @fields = split /\t/, $line;
1141
1142             $record->{ "REC_TYPE" }   = "BLAST";
1143             $record->{ "Q_ID" }       = $fields[ 0 ];
1144             $record->{ "S_ID" }       = $fields[ 1 ];
1145             $record->{ "IDENT" }      = $fields[ 2 ];
1146             $record->{ "ALIGN_LEN" }  = $fields[ 3 ];
1147             $record->{ "MISMATCHES" } = $fields[ 4 ];
1148             $record->{ "GAPS" }       = $fields[ 5 ];
1149             $record->{ "Q_BEG" }      = $fields[ 6 ] - 1; # BLAST is 1-based
1150             $record->{ "Q_END" }      = $fields[ 7 ] - 1; # BLAST is 1-based
1151             $record->{ "S_BEG" }      = $fields[ 8 ] - 1; # BLAST is 1-based
1152             $record->{ "S_END" }      = $fields[ 9 ] - 1; # BLAST is 1-based
1153             $record->{ "E_VAL" }      = $fields[ 10 ];
1154             $record->{ "BIT_SCORE" }  = $fields[ 11 ];
1155
1156             if ( $record->{ "S_BEG" } > $record->{ "S_END" } )
1157             {
1158                 $record->{ "STRAND" } = '-';
1159
1160                 ( $record->{ "S_BEG" }, $record->{ "S_END" } ) = ( $record->{ "S_END" }, $record->{ "S_BEG" } );
1161             }
1162             else
1163             {
1164                 $record->{ "STRAND" } = '+';
1165             }
1166
1167             Maasha::Biopieces::put_record( $record, $out );
1168
1169             goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1170
1171             $num++;
1172         }
1173
1174         close $data_in;
1175     }
1176
1177     NUM:
1178
1179     close $data_in if $data_in;
1180 }
1181
1182
1183 sub script_read_embl
1184 {
1185     # Martin A. Hansen, August 2007.
1186
1187     # Read EMBL format.
1188
1189     my ( $in,        # handle to in stream
1190          $out,       # handle to out stream
1191          $options,   # options hash
1192        ) = @_;
1193
1194     # Returns nothing.
1195
1196     my ( %options2, $file, $data_in, $num, $entry, $record );
1197
1198     map { $options2{ "keys" }{ $_ } = 1 }  @{ $options->{ "keys" } };
1199     map { $options2{ "feats" }{ $_ } = 1 } @{ $options->{ "feats" } };
1200     map { $options2{ "quals" }{ $_ } = 1 } @{ $options->{ "quals" } };
1201
1202     while ( $record = Maasha::Biopieces::get_record( $in ) ) {
1203         Maasha::Biopieces::put_record( $record, $out );
1204     }
1205
1206     $num = 1;
1207
1208     foreach $file ( @{ $options->{ "files" } } )
1209     {
1210         $data_in = Maasha::Common::read_open( $file );
1211
1212         while ( $entry = Maasha::EMBL::get_embl_entry( $data_in ) )
1213         {
1214             $record = Maasha::EMBL::parse_embl_entry( $entry, \%options2 );
1215
1216             my ( $feat, $feat2, $qual, $qual_val, $record_copy );
1217
1218             $record_copy = dclone $record;
1219
1220             delete $record_copy->{ "FT" };
1221
1222             Maasha::Biopieces::put_record( $record_copy, $out );
1223
1224             delete $record_copy->{ "SEQ" };
1225
1226             foreach $feat ( keys %{ $record->{ "FT" } } )
1227             {
1228                 $record_copy->{ "FEAT_TYPE" } = $feat;
1229
1230                 foreach $feat2 ( @{ $record->{ "FT" }->{ $feat } } )
1231                 {
1232                     foreach $qual ( keys %{ $feat2 } )
1233                     {
1234                         $qual_val = join "; ", @{ $feat2->{ $qual } };
1235
1236                         $qual =~ s/^_//;
1237                         $qual = uc $qual;
1238
1239                         $record_copy->{ $qual } = $qual_val;
1240                     }
1241
1242                     Maasha::Biopieces::put_record( $record_copy, $out );
1243                 }
1244             }
1245
1246             goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1247
1248             $num++;
1249         }
1250
1251         close $data_in;
1252     }
1253
1254     NUM:
1255
1256     close $data_in if $data_in;
1257 }
1258
1259
1260 sub script_read_stockholm
1261 {
1262     # Martin A. Hansen, August 2007.
1263
1264     # Read Stockholm format.
1265
1266     my ( $in,        # handle to in stream
1267          $out,       # handle to out stream
1268          $options,   # options hash
1269        ) = @_;
1270
1271     # Returns nothing.
1272
1273     my ( $data_in, $file, $num, $entry, $record, $record_anno, $record_align, $key, $seq );
1274
1275     while ( $record = Maasha::Biopieces::get_record( $in ) ) {
1276         Maasha::Biopieces::put_record( $record, $out );
1277     }
1278
1279     $num = 1;
1280
1281     foreach $file ( @{ $options->{ "files" } } )
1282     {
1283         $data_in = Maasha::Common::read_open( $file );
1284
1285         while ( $entry = Maasha::Stockholm::get_stockholm_entry( $data_in ) )
1286         {
1287             $record = Maasha::Stockholm::parse_stockholm_entry( $entry );
1288
1289             undef $record_anno;
1290
1291             foreach $key ( keys %{ $record->{ "GF" } } ) {
1292                 $record_anno->{ $key } = $record->{ "GF" }->{ $key };
1293             }
1294
1295             $record_anno->{ "ALIGN" } = $num;
1296
1297             Maasha::Biopieces::put_record( $record_anno, $out );
1298
1299             foreach $seq ( @{ $record->{ "ALIGN" } } )
1300             {
1301                 undef $record_align;
1302
1303                 $record_align = {
1304                     SEQ_NAME  => $seq->[ 0 ],
1305                     SEQ       => $seq->[ 1 ],
1306                 };
1307
1308                 Maasha::Biopieces::put_record( $record_align, $out );
1309             }
1310
1311             goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1312
1313             $num++;
1314         }
1315
1316         close $data_in;
1317     }
1318
1319     NUM:
1320
1321     close $data_in if $data_in;
1322 }
1323
1324
1325 sub script_read_phastcons
1326 {
1327     # Martin A. Hansen, December 2007.
1328
1329     # Read PhastCons format.
1330
1331     my ( $in,        # handle to in stream
1332          $out,       # handle to out stream
1333          $options,   # options hash
1334        ) = @_;
1335
1336     # Returns nothing.
1337
1338     my ( $data_in, $file, $num, $entry, @records, $record );
1339
1340     $options->{ "min" }       ||= 10;
1341     $options->{ "dist" }      ||= 25;
1342     $options->{ "threshold" } ||= 0.8;
1343     $options->{ "gap" }       ||= 5;
1344
1345     while ( $record = Maasha::Biopieces::get_record( $in ) ) {
1346         Maasha::Biopieces::put_record( $record, $out );
1347     }
1348
1349     $num = 1;
1350
1351     foreach $file ( @{ $options->{ "files" } } )
1352     {
1353         $data_in = Maasha::Common::read_open( $file );
1354
1355         while ( $entry = Maasha::UCSC::fixedstep_get_entry( $data_in ) )
1356         {
1357             @records = Maasha::UCSC::phastcons_parse_entry( $entry, $options );
1358
1359             foreach $record ( @records )
1360             {
1361                 $record->{ "REC_TYPE" } = "BED";
1362                 $record->{ "BED_LEN" }  = $record->{ "CHR_END" } - $record->{ "CHR_BEG" } + 1;
1363
1364                 Maasha::Biopieces::put_record( $record, $out );
1365
1366                 goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1367
1368                 $num++;
1369             }
1370         }
1371
1372         close $data_in;
1373     }
1374
1375     NUM:
1376
1377     close $data_in if $data_in;
1378 }
1379
1380
1381 sub script_read_soft
1382 {
1383     # Martin A. Hansen, December 2007.
1384
1385     # Read soft format.
1386     # http://www.ncbi.nlm.nih.gov/geo/info/soft2.html
1387
1388     my ( $in,        # handle to in stream
1389          $out,       # handle to out stream
1390          $options,   # options hash
1391        ) = @_;
1392
1393     # Returns nothing.
1394
1395     my ( $data_in, $file, $num, $records, $record, $soft_index, $fh, @platforms, $plat_table, @samples, $sample, $old_end, $skip );
1396
1397     while ( $record = Maasha::Biopieces::get_record( $in ) ) {
1398         Maasha::Biopieces::put_record( $record, $out );
1399     }
1400
1401     $num = 1;
1402
1403     foreach $file ( @{ $options->{ "files" } } )
1404     {
1405         print STDERR "Creating index for file: $file\n" if $options->{ "verbose" };
1406
1407         $soft_index = Maasha::NCBI::soft_index_file( $file );
1408
1409         $fh         = Maasha::Common::read_open( $file );
1410
1411         @platforms  = grep { $_->{ "SECTION" } =~ /PLATFORM/ } @{ $soft_index };
1412
1413         print STDERR "Getting platform tables for file: $file\n" if $options->{ "verbose" };
1414
1415         $plat_table = Maasha::NCBI::soft_get_platform( $fh, $platforms[ 0 ]->{ "LINE_BEG" }, $platforms[ -1 ]->{ "LINE_END" } );
1416
1417         @samples    = grep { $_->{ "SECTION" } =~ /SAMPLE/ } @{ $soft_index };
1418
1419         $old_end    = $platforms[ -1 ]->{ "LINE_END" };
1420
1421         foreach $sample ( @samples )
1422         {
1423             $skip = 0;
1424             $skip = 1 if ( $options->{ "samples" } and grep { $sample->{ "SECTION" } !~ /$_/ } @{ $options->{ "samples" } } );
1425
1426             print STDERR "Getting samples for dataset: $sample->{ 'SECTION' }\n" if $options->{ "verbose" } and not $skip;
1427
1428             $records = Maasha::NCBI::soft_get_sample( $fh, $plat_table, $sample->{ "LINE_BEG" } - $old_end - 1, $sample->{ "LINE_END" } - $old_end - 1, $skip );
1429
1430             foreach $record ( @{ $records } )
1431             {
1432                 Maasha::Biopieces::put_record( $record, $out );
1433
1434                 goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1435
1436                 $num++;
1437             }
1438
1439             $old_end = $sample->{ "LINE_END" };
1440         }
1441
1442         close $fh;
1443     }
1444
1445     NUM:
1446
1447     close $data_in if $data_in;
1448     close $fh if $fh;
1449 }
1450
1451
1452 sub script_read_gff
1453 {
1454     # Martin A. Hansen, February 2008.
1455
1456     # Read soft format.
1457     # http://www.ncbi.nlm.nih.gov/geo/info/soft2.html
1458
1459     my ( $in,        # handle to in stream
1460          $out,       # handle to out stream
1461          $options,   # options hash
1462        ) = @_;
1463
1464     # Returns nothing.
1465
1466     my ( $data_in, $file, $fh, $num, $record, $entry );
1467
1468     while ( $record = Maasha::Biopieces::get_record( $in ) ) {
1469         Maasha::Biopieces::put_record( $record, $out );
1470     }
1471
1472     $num = 1;
1473
1474     foreach $file ( @{ $options->{ "files" } } )
1475     {
1476         $fh = Maasha::Common::read_open( $file );
1477
1478         while ( $entry = Maasha::GFF::get_entry( $fh ) )
1479         {
1480             Maasha::Biopieces::put_record( $entry, $out );
1481
1482             goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1483
1484             $num++;
1485         }
1486
1487         close $fh;
1488     }
1489
1490     NUM:
1491
1492     close $data_in if $data_in;
1493 }
1494
1495
1496 sub script_read_2bit
1497 {
1498     # Martin A. Hansen, March 2008.
1499
1500     # Read sequences from 2bit file.
1501
1502     my ( $in,        # handle to in stream
1503          $out,       # handle to out stream
1504          $options,   # options hash
1505        ) = @_;
1506
1507     # Returns nothing.
1508
1509     my ( $record, $file, $data_in, $mask, $toc, $line, $num );
1510
1511     $mask = 1 if not $options->{ "no_mask" };
1512
1513     while ( $record = Maasha::Biopieces::get_record( $in ) ) {
1514         Maasha::Biopieces::put_record( $record, $out );
1515     }
1516
1517     $num = 1;
1518
1519     foreach $file ( @{ $options->{ "files" } } )
1520     {
1521         $data_in = Maasha::Common::read_open( $file );
1522
1523         $toc = Maasha::TwoBit::twobit_get_TOC( $data_in );
1524
1525         foreach $line ( @{ $toc } )
1526         {
1527             $record->{ "SEQ_NAME" } = $line->[ 0 ];
1528             $record->{ "SEQ" }      = Maasha::TwoBit::twobit_get_seq( $data_in, $line->[ 1 ], undef, undef, $mask );
1529             $record->{ "SEQ_LEN" }  = length $record->{ "SEQ" };
1530
1531             Maasha::Biopieces::put_record( $record, $out );
1532
1533             goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1534
1535             $num++;
1536         }
1537
1538         close $data_in;
1539     }
1540
1541     NUM:
1542
1543     close $data_in if $data_in;
1544 }
1545
1546
1547 sub script_read_solexa
1548 {
1549     # Martin A. Hansen, March 2008.
1550
1551     # Read Solexa sequence reads from file.
1552
1553     my ( $in,        # handle to in stream
1554          $out,       # handle to out stream
1555          $options,   # options hash
1556        ) = @_;
1557
1558     # Returns nothing.
1559
1560     my ( $record, $file, $data_in, $entry, $num, @seqs, @scores, $i );
1561
1562     $options->{ "format" }  ||= "octal";
1563     $options->{ "quality" } ||= 20;
1564
1565     while ( $record = Maasha::Biopieces::get_record( $in ) ) {
1566         Maasha::Biopieces::put_record( $record, $out );
1567     }
1568
1569     $num = 1;
1570
1571     foreach $file ( @{ $options->{ "files" } } )
1572     {
1573         $data_in = Maasha::Common::read_open( $file );
1574
1575         if ( $options->{ "format" } eq "octal" )
1576         {
1577             while ( $entry = Maasha::Solexa::solexa_get_entry_octal( $data_in ) )
1578             {
1579                 $record = Maasha::Solexa::solexa2biopiece( $entry, $options->{ "quality" } );
1580
1581                 Maasha::Biopieces::put_record( $record, $out );
1582
1583                 goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1584
1585                 $num++;
1586             }
1587         }
1588         else
1589         {
1590             while ( $entry = Maasha::Solexa::solexa_get_entry_decimal( $data_in ) )
1591             {
1592                 $record = Maasha::Solexa::solexa2biopiece( $entry, $options->{ "quality" } );
1593
1594                 Maasha::Biopieces::put_record( $record, $out );
1595
1596                 goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1597
1598                 $num++;
1599             }
1600         }
1601
1602         close $data_in;
1603     }
1604
1605     NUM:
1606
1607     close $data_in if $data_in;
1608 }
1609
1610
1611 sub script_read_solid
1612 {
1613     # Martin A. Hansen, April 2008.
1614
1615     # Read Solid sequence from file.
1616
1617     my ( $in,        # handle to in stream
1618          $out,       # handle to out stream
1619          $options,   # options hash
1620        ) = @_;
1621
1622     # Returns nothing.
1623
1624     my ( $record, $file, $data_in, $line, $num, $seq_name, $seq_cs, $seq_qual, @scores, @seqs, $i );
1625
1626     $options->{ "quality" } ||= 15;
1627
1628     while ( $record = Maasha::Biopieces::get_record( $in ) ) {
1629         Maasha::Biopieces::put_record( $record, $out );
1630     }
1631
1632     $num = 1;
1633
1634     foreach $file ( @{ $options->{ "files" } } )
1635     {
1636         $data_in = Maasha::Common::read_open( $file );
1637
1638         while ( $line = <$data_in> )
1639         {
1640             chomp $line;
1641
1642             ( $seq_name, $seq_cs, $seq_qual ) = split /\t/, $line;
1643
1644             @scores = split /,/, $seq_qual;
1645             @seqs   = split //, Maasha::Solid::color_space2seq( $seq_cs );
1646
1647             for ( $i = 0; $i < @seqs; $i++ ) {
1648                 $seqs[ $i ] = lc $seqs[ $i ] if $scores[ $i ] < $options->{ "quality" };
1649             }
1650
1651             $record = {
1652                 REC_TYPE   => 'SOLID',
1653                 SEQ_NAME   => $seq_name,
1654                 SEQ_CS     => $seq_cs,
1655                 SEQ_QUAL   => join( ";", @scores ),
1656                 SEQ_LEN    => length $seq_cs,
1657                 SEQ        => join( "", @seqs ),
1658                 SCORE_MEAN => sprintf( "%.2f", Maasha::Calc::mean( \@scores ) ),
1659             };
1660
1661             Maasha::Biopieces::put_record( $record, $out );
1662
1663             goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1664
1665             $num++;
1666         }
1667
1668         close $data_in;
1669     }
1670
1671     NUM:
1672
1673     close $data_in if $data_in;
1674 }
1675
1676
1677 sub script_read_mysql
1678 {
1679     # Martin A. Hansen, May 2008.
1680
1681     # Read a MySQL query into stream.
1682
1683     my ( $in,        # handle to in stream
1684          $out,       # handle to out stream
1685          $options,   # options hash
1686        ) = @_;
1687
1688     # Returns nothing.
1689
1690     my ( $record, $dbh, $results );
1691
1692     $options->{ "user" }     ||= Maasha::UCSC::ucsc_get_user();
1693     $options->{ "password" } ||= Maasha::UCSC::ucsc_get_password();
1694
1695     while ( $record = Maasha::Biopieces::get_record( $in ) ) {
1696         Maasha::Biopieces::put_record( $record, $out );
1697     }
1698
1699     $dbh = Maasha::SQL::connect( $options->{ "database" }, $options->{ "user" }, $options->{ "password" } );
1700
1701     $results = Maasha::SQL::query_hashref_list( $dbh, $options->{ "query" } );
1702
1703     Maasha::SQL::disconnect( $dbh );
1704
1705     map { Maasha::Biopieces::put_record( $_ ) } @{ $results };
1706 }
1707
1708
1709 sub script_read_ucsc_config
1710 {
1711     # Martin A. Hansen, November 2008.
1712
1713     # Read track entries from UCSC Genome Browser '.ra' files.
1714
1715     my ( $in,        # handle to in stream
1716          $out,       # handle to out stream
1717          $options,   # options hash
1718        ) = @_;
1719
1720     # Returns nothing.
1721
1722     my ( $record, $file, $data_in, $entry, $num );
1723
1724     while ( $record = Maasha::Biopieces::get_record( $in ) ) {
1725         Maasha::Biopieces::put_record( $record, $out );
1726     }
1727
1728     $num = 1;
1729
1730     foreach $file ( @{ $options->{ "files" } } )
1731     {
1732         $data_in = Maasha::Common::read_open( $file );
1733
1734         while ( $record = Maasha::UCSC::ucsc_config_entry_get( $data_in ) )
1735         {
1736             $record->{ 'REC_TYPE' } = "UCSC Config";
1737
1738             Maasha::Biopieces::put_record( $record, $out );
1739
1740             goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1741
1742             $num++;
1743         }
1744
1745         close $data_in;
1746     }
1747
1748     NUM:
1749
1750     close $data_in if $data_in;
1751 }
1752
1753
1754 sub script_assemble_tag_contigs
1755 {
1756     # Martin A. Hansen, November 2008.
1757
1758     # Assemble tags from the stream into
1759     # tag contigs.
1760
1761     # The current implementation is quite
1762     # slow because of heavy use of temporary
1763     # files.
1764
1765     my ( $in,        # handle to in stream
1766          $out,       # handle to out stream
1767          $options,   # options hash
1768        ) = @_;
1769
1770     # Returns nothing.
1771
1772     my ( $bed_file, $tag_file, $fh_in, $fh_out, $cols, $record, $bed_entry, $file_hash, $chr, $strand );
1773
1774     $bed_file = "$BP_TMP/assemble_tag_contigs.bed";
1775     $fh_out   = Maasha::Filesys::file_write_open( $bed_file );
1776     $cols     = 6; # we only need the first 6 BED columns
1777
1778     while ( $record = Maasha::Biopieces::get_record( $in ) )
1779     {
1780         if ( $bed_entry = Maasha::UCSC::BED::biopiece2bed( $record, $cols ) )
1781         {
1782             $strand = $record->{ 'STRAND' } || '+';
1783
1784             Maasha::UCSC::BED::bed_entry_put( $bed_entry, $fh_out, $cols, $options->{ 'check' } );
1785         }
1786     }
1787
1788     close $fh_out;
1789
1790     $file_hash = Maasha::UCSC::BED::bed_file_split_on_chr( $bed_file, $BP_TMP, $cols );
1791
1792     unlink $bed_file;
1793
1794     foreach $chr ( sort keys %{ $file_hash } )
1795     {
1796         $bed_file = $file_hash->{ $chr };
1797         $tag_file = "$bed_file.tc";
1798
1799         Maasha::Common::run( "bed2tag_contigs", "< $bed_file > $tag_file" );
1800
1801         $fh_in = Maasha::Filesys::file_read_open( $tag_file );
1802
1803         while ( $bed_entry = Maasha::UCSC::BED::bed_entry_get( $fh_in, $cols, $options->{ 'check' } ) )
1804         {
1805             if ( $record = Maasha::UCSC::BED::bed2biopiece( $bed_entry ) ) {
1806                 Maasha::Biopieces::put_record( $record, $out );
1807             }
1808         }
1809
1810         close $fh_in;
1811
1812         unlink $bed_file;
1813         unlink $tag_file;
1814     }
1815 }
1816
1817
1818 sub script_length_seq
1819 {
1820     # Martin A. Hansen, August 2007.
1821
1822     # Determine the length of sequences in stream.
1823
1824     my ( $in,        # handle to in stream
1825          $out,       # handle to out stream
1826          $options,   # options hash
1827        ) = @_;
1828
1829     # Returns nothing.
1830
1831     my ( $record, $total );
1832
1833     while ( $record = Maasha::Biopieces::get_record( $in ) )
1834     {
1835         if ( $record->{ "SEQ" } )
1836         {
1837             $record->{ "SEQ_LEN" } = length $record->{ "SEQ" };
1838             $total += $record->{ "SEQ_LEN" };
1839         }
1840
1841         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
1842     }
1843
1844     Maasha::Biopieces::put_record( { TOTAL_SEQ_LEN => $total }, $out );
1845 }
1846
1847
1848 sub script_shuffle_seq
1849 {
1850     # Martin A. Hansen, December 2007.
1851
1852     # Shuffle sequences in stream.
1853
1854     my ( $in,    # handle to in stream
1855          $out,   # handle to out stream
1856        ) = @_;
1857
1858     # Returns nothing.
1859
1860     my ( $record );
1861
1862     while ( $record = Maasha::Biopieces::get_record( $in ) )
1863     {
1864         $record->{ "SEQ" } = Maasha::Seq::seq_shuffle( $record->{ "SEQ" } ) if $record->{ "SEQ" };
1865
1866         Maasha::Biopieces::put_record( $record, $out );
1867     }
1868 }
1869
1870
1871 sub script_analyze_seq
1872 {
1873     # Martin A. Hansen, August 2007.
1874
1875     # Analyze sequence composition of sequences in stream.
1876
1877     my ( $in,     # handle to in stream
1878          $out,    # handle to out stream
1879        ) = @_;
1880
1881     # Returns nothing.
1882
1883     my ( $record, $analysis );
1884
1885     while ( $record = Maasha::Biopieces::get_record( $in ) )
1886     {
1887         if ( $record->{ "SEQ" } )
1888         {
1889             $analysis = Maasha::Seq::seq_analyze( $record->{ "SEQ" } );
1890
1891             map { $record->{ $_ } = $analysis->{ $_ } } keys %{ $analysis };
1892         }
1893
1894         Maasha::Biopieces::put_record( $record, $out );
1895     }
1896 }
1897
1898
1899 sub script_analyze_tags
1900 {
1901     # Martin A. Hansen, August 2008.
1902
1903     # Analyze sequence tags in stream.
1904
1905     my ( $in,     # handle to in stream
1906          $out,    # handle to out stream
1907        ) = @_;
1908
1909     # Returns nothing.
1910
1911     my ( $record, $analysis, %len_hash, %clone_hash, $clones, $key, $tag_record );
1912
1913     while ( $record = Maasha::Biopieces::get_record( $in ) )
1914     {
1915         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } )
1916         {
1917             if ( $record->{ "SEQ_NAME" } =~ /_(\d+)$/ )
1918             {
1919                 $clones = $1;
1920
1921                 $len_hash{ length( $record->{ "SEQ" } ) }++;
1922                 $clone_hash{ length( $record->{ "SEQ" } ) } += $clones;
1923             }
1924         }
1925         elsif ( $record->{ "Q_ID" } and $record->{ "BED_LEN" } )
1926         {
1927             if ( $record->{ "Q_ID" } =~ /_(\d+)$/ )
1928             {
1929                 $clones = $1;
1930
1931                 $len_hash{ $record->{ "BED_LEN" } }++;
1932                 $clone_hash{ $record->{ "BED_LEN" } } += $clones;
1933             }
1934         }
1935     }
1936
1937     foreach $key ( sort { $a <=> $b } keys %len_hash )
1938     {
1939         $tag_record->{ "TAG_LEN" }    = $key;
1940         $tag_record->{ "TAG_COUNT" }  = $len_hash{ $key };
1941         $tag_record->{ "TAG_CLONES" } = $clone_hash{ $key };
1942
1943         Maasha::Biopieces::put_record( $tag_record, $out );
1944     }
1945 }
1946
1947
1948 sub script_complexity_seq
1949 {
1950     # Martin A. Hansen, May 2008.
1951
1952     # Generates an index calculated as the most common di-residue over
1953     # the sequence length for all sequences in stream.
1954
1955     my ( $in,     # handle to in stream
1956          $out,    # handle to out stream
1957        ) = @_;
1958
1959     # Returns nothing.
1960
1961     my ( $record, $index );
1962
1963     while ( $record = Maasha::Biopieces::get_record( $in ) )
1964     {
1965         $record->{ "SEQ_COMPLEXITY" } = sprintf( "%.2f", Maasha::Seq::seq_complexity( $record->{ "SEQ" } ) ) if $record->{ "SEQ" };
1966
1967         Maasha::Biopieces::put_record( $record, $out );
1968     }
1969 }
1970
1971
1972 sub script_oligo_freq
1973 {
1974     # Martin A. Hansen, August 2007.
1975
1976     # Determine the length of sequences in stream.
1977
1978     my ( $in,        # handle to in stream
1979          $out,       # handle to out stream
1980          $options,   # options hash
1981        ) = @_;
1982
1983     # Returns nothing.
1984
1985     my ( $record, %oligos, @freq_table );
1986
1987     $options->{ "word_size" } ||= 7;
1988
1989     while ( $record = Maasha::Biopieces::get_record( $in ) )
1990     {
1991         if ( $record->{ "SEQ" } )
1992         {
1993             map { $oligos{ $_ }++ } Maasha::Seq::seq2oligos( \$record->{ "SEQ" }, $options->{ "word_size" } );
1994
1995             if ( not $options->{ "all" } )
1996             {
1997                 @freq_table = Maasha::Seq::oligo_freq( \%oligos );
1998
1999                 map { Maasha::Biopieces::put_record( $_, $out ) } @freq_table;
2000
2001                 undef %oligos;
2002             }
2003         }
2004
2005         Maasha::Biopieces::put_record( $record, $out );
2006     }
2007
2008     if ( $options->{ "all" } )
2009     {
2010         @freq_table = Maasha::Seq::oligo_freq( \%oligos );
2011
2012         map { Maasha::Biopieces::put_record( $_, $out ) } @freq_table;
2013     }
2014 }
2015
2016
2017 sub script_create_weight_matrix
2018 {
2019     # Martin A. Hansen, August 2007.
2020
2021     # Creates a weight matrix from an alignmnet.
2022
2023     my ( $in,        # handle to in stream
2024          $out,       # handle to out stream
2025          $options,   # options hash
2026        ) = @_;
2027
2028     # Returns nothing.
2029
2030     my ( $record, $count, $i, $res, %freq_hash, %res_hash, $freq );
2031
2032     $count = 0;
2033
2034     while ( $record = Maasha::Biopieces::get_record( $in ) )
2035     {
2036         if ( $record->{ "SEQ" } )
2037         {
2038             for ( $i = 0; $i < length $record->{ "SEQ" }; $i++ )
2039             {
2040                 $res = substr $record->{ "SEQ" }, $i, 1;
2041
2042                 $freq_hash{ $i }{ $res }++;
2043                 $res_hash{ $res } = 1;
2044             }
2045
2046             $count++;
2047         }
2048         else
2049         {
2050             Maasha::Biopieces::put_record( $record, $out );
2051         }
2052     }
2053
2054     foreach $res ( sort keys %res_hash )
2055     {
2056         undef $record;
2057
2058         $record->{ "V0" } = $res;
2059
2060         for ( $i = 0; $i < keys %freq_hash; $i++ )
2061         {
2062             $freq = $freq_hash{ $i }{ $res } || 0;
2063
2064             if ( $options->{ "percent" } ) {
2065                 $freq = sprintf( "%.0f", 100 * $freq / $count ) if $freq > 0;
2066             }
2067
2068             $record->{ "V" . ( $i + 1 ) } = $freq;
2069         }
2070
2071         Maasha::Biopieces::put_record( $record, $out );
2072     }
2073 }
2074
2075
2076 sub script_calc_bit_scores
2077 {
2078     # Martin A. Hansen, March 2007.
2079
2080     # Calculates the bit scores for each position from an alignmnet in the stream.
2081
2082     my ( $in,        # handle to in stream
2083          $out,       # handle to out stream
2084        ) = @_;
2085
2086     # Returns nothing.
2087
2088     my ( $record, $type, $count, $i, $res, %freq_hash, $bit_max, $bit_height, $bit_diff );
2089
2090     $count = 0;
2091
2092     while ( $record = Maasha::Biopieces::get_record( $in ) )
2093     {
2094         if ( $record->{ "SEQ" } )
2095         {
2096             $type = Maasha::Seq::seq_guess_type( $record->{ "SEQ" } ) if not $type;
2097
2098             for ( $i = 0; $i < length $record->{ "SEQ" }; $i++ )
2099             {
2100                 $res = substr $record->{ "SEQ" }, $i, 1;
2101
2102                 next if $res =~ /-|_|~|\./;
2103
2104                 $freq_hash{ $i }{ $res }++;
2105             }
2106
2107             $count++;
2108         }
2109         else
2110         {
2111             Maasha::Biopieces::put_record( $record, $out );
2112         }
2113     }
2114
2115     undef $record;
2116
2117     if ( $type eq "protein" ) {
2118         $bit_max = 4;
2119     } else {
2120         $bit_max = 2;
2121     }
2122
2123     for ( $i = 0; $i < keys %freq_hash; $i++ )
2124     {
2125         $bit_height = Maasha::Seq::seqlogo_calc_bit_height( $freq_hash{ $i }, $count );
2126
2127         $bit_diff = $bit_max - $bit_height;
2128
2129         $record->{ "V" . ( $i ) } = sprintf( "%.2f", $bit_diff );
2130     }
2131
2132     Maasha::Biopieces::put_record( $record, $out );
2133 }
2134
2135
2136 sub script_calc_fixedstep
2137 {
2138     # Martin A. Hansen, September 2008.
2139
2140     # Calculates fixedstep entries from data in the stream.
2141
2142     my ( $in,        # handle to in stream
2143          $out,       # handle to out stream
2144          $options,   # options hash
2145        ) = @_;
2146
2147     # Returns nothing.
2148
2149     my ( $bed_file, $fh_in, $fh_out, $record, $file_hash, $chr, $bed_entry, $fixedstep_file, $fixedstep_entry );
2150
2151     $bed_file = "$BP_TMP/calc_fixedstep.bed";
2152     $fh_out   = Maasha::Filesys::file_write_open( $bed_file );
2153
2154     while ( $record = Maasha::Biopieces::get_record( $in ) )
2155     {
2156         if ( $bed_entry = Maasha::UCSC::BED::biopiece2bed( $record ) ) {
2157             Maasha::UCSC::BED::bed_entry_put( $bed_entry, $fh_out, undef, $options->{ 'check' } );
2158         }
2159     }
2160
2161     close $fh_out;
2162
2163     $file_hash = Maasha::UCSC::BED::bed_file_split_on_chr( $bed_file, $BP_TMP );
2164
2165     unlink $bed_file;
2166
2167     foreach $chr ( sort keys %{ $file_hash } )
2168     {
2169         $bed_file       = $file_hash->{ $chr };
2170
2171         $fixedstep_file = Maasha::UCSC::Wiggle::fixedstep_calc( $bed_file, $chr, $options->{ 'use_score' }, $options->{ 'use_log10' } );
2172
2173         #$fixedstep_file = "$bed_file.fixedstep";
2174
2175         # Maasha::Common::run( "bed2fixedstep", "< $bed_file > $fixedstep_file" );
2176
2177         $fh_in = Maasha::Filesys::file_read_open( $fixedstep_file );
2178
2179         while ( $fixedstep_entry = Maasha::UCSC::Wiggle::fixedstep_entry_get( $fh_in ) )
2180         {
2181             if ( $record = Maasha::UCSC::Wiggle::fixedstep2biopiece( $fixedstep_entry ) ) {
2182                 Maasha::Biopieces::put_record( $record, $out );
2183             }
2184         }
2185
2186         close $fh_in;
2187
2188         unlink $bed_file;
2189         unlink $fixedstep_file;
2190     }
2191 }
2192
2193
2194 sub script_remove_indels
2195 {
2196     # Martin A. Hansen, August 2007.
2197
2198     # Remove indels from sequences in stream.
2199
2200     my ( $in,     # handle to in stream
2201          $out,    # handle to out stream
2202        ) = @_;
2203
2204     # Returns nothing.
2205
2206     my ( $record );
2207
2208     while ( $record = Maasha::Biopieces::get_record( $in ) )
2209     {
2210         $record->{ 'SEQ' } =~ tr/-~.//d if $record->{ "SEQ" };
2211
2212         Maasha::Biopieces::put_record( $record, $out );
2213     }
2214 }
2215
2216
2217 sub script_transliterate_seq
2218 {
2219     # Martin A. Hansen, August 2007.
2220
2221     # Transliterate chars from sequence in record.
2222
2223     my ( $in,        # handle to in stream
2224          $out,       # handle to out stream
2225          $options,   # options hash
2226        ) = @_;
2227
2228     # Returns nothing.
2229
2230     my ( $record, $search, $replace, $delete );
2231
2232     $search  = $options->{ "search" }  || "";
2233     $replace = $options->{ "replace" } || "";
2234     $delete  = $options->{ "delete" }  || "";
2235
2236     while ( $record = Maasha::Biopieces::get_record( $in ) )
2237     {
2238         if ( $record->{ "SEQ" } )
2239         {
2240             if ( $search and $replace ) {
2241                 eval "\$record->{ 'SEQ' } =~ tr/$search/$replace/";
2242             } elsif ( $delete ) {
2243                 eval "\$record->{ 'SEQ' } =~ tr/$delete//d";
2244             }
2245         }
2246
2247         Maasha::Biopieces::put_record( $record, $out );
2248     }
2249 }
2250
2251
2252 sub script_transliterate_vals
2253 {
2254     # Martin A. Hansen, April 2008.
2255
2256     # Transliterate chars from values in record.
2257
2258     my ( $in,        # handle to in stream
2259          $out,       # handle to out stream
2260          $options,   # options hash
2261        ) = @_;
2262
2263     # Returns nothing.
2264
2265     my ( $record, $search, $replace, $delete, $key );
2266
2267     $search  = $options->{ "search" }  || "";
2268     $replace = $options->{ "replace" } || "";
2269     $delete  = $options->{ "delete" }  || "";
2270
2271     while ( $record = Maasha::Biopieces::get_record( $in ) )
2272     {
2273         foreach $key ( @{ $options->{ "keys" } } )
2274         {
2275             if ( exists $record->{ $key } )
2276             {
2277                 if ( $search and $replace ) {
2278                     eval "\$record->{ $key } =~ tr/$search/$replace/";
2279                 } elsif ( $delete ) {
2280                     eval "\$record->{ $key } =~ tr/$delete//d";
2281                 }
2282             }
2283         }
2284
2285         Maasha::Biopieces::put_record( $record, $out );
2286     }
2287 }
2288
2289
2290 sub script_translate_seq
2291 {
2292     # Martin A. Hansen, February 2008.
2293
2294     # Translate DNA sequence into protein sequence.
2295
2296     my ( $in,        # handle to in stream
2297          $out,       # handle to out stream
2298          $options,   # options hash
2299        ) = @_;
2300
2301     # Returns nothing.
2302
2303     my ( $record, $frame, %new_record );
2304
2305     $options->{ "frames" } ||= [ 1, 2, 3, -1, -2, -3 ];
2306
2307     while ( $record = Maasha::Biopieces::get_record( $in ) )
2308     {
2309         if ( $record->{ "SEQ" } )
2310         {
2311             if ( Maasha::Seq::seq_guess_type( $record->{ "SEQ" } ) eq "dna" )
2312             {
2313                 foreach $frame ( @{ $options->{ "frames" } } )
2314                 {
2315                     %new_record = %{ $record };
2316
2317                     $new_record{ "SEQ" }     = Maasha::Seq::translate( $record->{ "SEQ" }, $frame );
2318                     $new_record{ "SEQ_LEN" } = length $new_record{ "SEQ" };
2319                     $new_record{ "FRAME" }   = $frame;
2320
2321                     Maasha::Biopieces::put_record( \%new_record, $out );
2322                 }
2323             }
2324         }
2325         else
2326         {
2327             Maasha::Biopieces::put_record( $record, $out );
2328         }
2329     }
2330 }
2331
2332
2333 sub script_get_genome_align
2334 {
2335     # Martin A. Hansen, April 2008.
2336
2337     # Gets a subalignment from a multiple genome alignment.
2338
2339     my ( $in,        # handle to in stream
2340          $out,       # handle to out stream
2341          $options,   # options hash
2342        ) = @_;
2343
2344     # Returns nothing.
2345
2346     my ( $record, $maf_track, $align, $align_num, $beg, $end, $len, $entry );
2347
2348     $options->{ "strand" } ||= "+";
2349
2350     $align_num = 1;
2351
2352     $maf_track = Maasha::Config::maf_track( $options->{ "genome" } );
2353
2354     if ( $options->{ "chr" } and $options->{ "beg" } and ( $options->{ "end" } or $options->{ "len" } ) )
2355     {
2356         $beg = $options->{ "beg" } - 1;
2357
2358         if ( $options->{ "end" } ) {
2359             $end = $options->{ "end" };
2360         } elsif ( $options->{ "len" } ) {
2361             $end = $beg + $options->{ "len" };
2362         }
2363
2364         $align = Maasha::UCSC::maf_extract( $BP_TMP, $options->{ "genome" }, $maf_track, $options->{ "chr" }, $beg, $end, $options->{ "strand" } );
2365
2366         foreach $entry ( @{ $align } )
2367         {
2368             $entry->{ "CHR" }     = $record->{ "CHR" };
2369             $entry->{ "CHR_BEG" } = $record->{ "CHR_BEG" };
2370             $entry->{ "CHR_END" } = $record->{ "CHR_END" };
2371             $entry->{ "STRAND" }  = $record->{ "STRAND" } || '+';
2372             $entry->{ "Q_ID" }    = $record->{ "Q_ID" };
2373             $entry->{ "SCORE" }   = $record->{ "SCORE" };
2374
2375             Maasha::Biopieces::put_record( $entry, $out );
2376         }
2377     }
2378
2379     while ( $record = Maasha::Biopieces::get_record( $in ) )
2380     {
2381         if ( $record->{ "REC_TYPE" } eq "BED" )
2382         {
2383             $align = Maasha::UCSC::maf_extract( $BP_TMP, $options->{ "genome" }, $maf_track, $record->{ "CHR" }, $record->{ "CHR_BEG" }, $record->{ "CHR_END" }, $record->{ "STRAND" } );
2384         }
2385         elsif ( $record->{ "REC_TYPE" } eq "VMATCH" )
2386         {
2387             $align = Maasha::UCSC::maf_extract( $BP_TMP, $options->{ "genome" }, $maf_track, $record->{ "S_ID" }, $record->{ "S_BEG" }, $record->{ "S_END" } + 1, $record->{ "STRAND" } );
2388         }
2389         elsif ( $record->{ "REC_TYPE" } eq "PSL" )
2390         {
2391             $align = Maasha::UCSC::maf_extract( $BP_TMP, $options->{ "genome" }, $maf_track, $record->{ "S_ID" }, $record->{ "S_BEG" }, $record->{ "S_END" }, $record->{ "STRAND" } );
2392         }
2393         elsif ( $record->{ "REC_TYPE" } eq "BLAST" )
2394         {
2395             $align = Maasha::UCSC::maf_extract( $BP_TMP, $options->{ "genome" }, $maf_track, $record->{ "S_ID" }, $record->{ "S_BEG" }, $record->{ "S_END" }, $record->{ "STRAND" } );
2396         }
2397
2398         foreach $entry ( @{ $align } )
2399         {
2400             $entry->{ "CHR" }     = $record->{ "CHR" };
2401             $entry->{ "CHR_BEG" } = $record->{ "CHR_BEG" };
2402             $entry->{ "CHR_END" } = $record->{ "CHR_END" };
2403             $entry->{ "STRAND" }  = $record->{ "STRAND" };
2404             $entry->{ "Q_ID" }    = $record->{ "Q_ID" };
2405             $entry->{ "SCORE" }   = $record->{ "SCORE" };
2406
2407             Maasha::Biopieces::put_record( $entry, $out );
2408         }
2409
2410         $align_num++;
2411     }
2412 }
2413
2414
2415 sub script_get_genome_phastcons
2416 {
2417     # Martin A. Hansen, February 2008.
2418
2419     # Get phastcons scores from genome intervals.
2420
2421     my ( $in,        # handle to in stream
2422          $out,       # handle to out stream
2423          $options,   # options hash
2424        ) = @_;
2425
2426     # Returns nothing.
2427
2428     my ( $phastcons_file, $phastcons_index, $index, $fh_phastcons, $scores, $record );
2429
2430     $options->{ "flank" } ||= 0;
2431
2432     $phastcons_file  = Maasha::Config::genome_phastcons( $options->{ "genome" } );
2433     $phastcons_index = Maasha::Config::genome_phastcons_index( $options->{ "genome" } );
2434
2435     $index           = Maasha::UCSC::fixedstep_index_retrieve( $phastcons_index );
2436     $fh_phastcons    = Maasha::Common::read_open( $phastcons_file );
2437
2438     if ( defined $options->{ "chr" } and defined $options->{ "beg" } and ( defined $options->{ "end" } or defined $options->{ "len" } ) )
2439     {
2440         $options->{ "beg" } -= 1;   # request is 1-based
2441         $options->{ "end" } -= 1;   # request is 1-based
2442
2443         if ( $options->{ "len" } ) {
2444             $options->{ "end" } = $options->{ "beg" } + $options->{ "len" } - 1;
2445         }
2446
2447         $scores = Maasha::UCSC::fixedstep_index_lookup( $index, $fh_phastcons, $options->{ "chr" }, $options->{ "beg" }, $options->{ "end" }, $options->{ "flank" } );
2448
2449         $record->{ "CHR" }       = $options->{ "chr" };
2450         $record->{ "CHR_BEG" }   = $options->{ "beg" } - $options->{ "flank" };
2451         $record->{ "CHR_END" }   = $options->{ "end" } + $options->{ "flank" };
2452
2453         $record->{ "PHASTCONS" }   = join ",", @{ $scores };
2454         $record->{ "PHAST_COUNT" } = scalar @{ $scores };  # DEBUG
2455
2456         Maasha::Biopieces::put_record( $record, $out );
2457     }
2458
2459     while ( $record = Maasha::Biopieces::get_record( $in ) )
2460     {
2461         if ( $record->{ "REC_TYPE" } eq "BED" )
2462         {
2463             $scores = Maasha::UCSC::fixedstep_index_lookup( $index, $fh_phastcons, $record->{ "CHR" }, $record->{ "CHR_BEG" }, $record->{ "CHR_END" }, $options->{ "flank" } );
2464         }
2465         elsif ( $record->{ "REC_TYPE" } eq "PSL" )
2466         {
2467             $scores = Maasha::UCSC::fixedstep_index_lookup( $index, $fh_phastcons, $record->{ "S_ID" }, $record->{ "S_BEG" }, $record->{ "S_END" }, $options->{ "flank" } );
2468         }
2469         elsif ( $record->{ "REC_TYPE" } eq "BLAST" )
2470         {
2471             $scores = Maasha::UCSC::fixedstep_index_lookup( $index, $fh_phastcons, $record->{ "S_ID" }, $record->{ "S_BEG" }, $record->{ "S_END" }, $options->{ "flank" } );
2472         }
2473
2474         $record->{ "PHASTCONS" } = join ",", @{ $scores } if @{ $scores };
2475 #        $record->{ "PHAST_COUNT" } = @{ $scores } if @{ $scores };  # DEBUG
2476
2477         Maasha::Biopieces::put_record( $record, $out );
2478     }
2479
2480     close $fh_phastcons if $fh_phastcons;
2481 }
2482
2483
2484 sub script_fold_seq
2485 {
2486     # Martin A. Hansen, December 2007.
2487
2488     # Folds sequences in stream into secondary structures.
2489
2490     my ( $in,     # handle to in stream
2491          $out,    # handle to out stream
2492        ) = @_;
2493
2494     # Returns nothing.
2495
2496     my ( $record, $type, $struct, $index );
2497
2498     while ( $record = Maasha::Biopieces::get_record( $in ) )
2499     {
2500         if ( $record->{ "SEQ" } )
2501         {
2502             if ( not $type ) {
2503                 $type = Maasha::Seq::seq_guess_type( $record->{ "SEQ" } );
2504             }
2505
2506             if ( $type ne "protein" )
2507             {
2508                 ( $struct, $index ) = Maasha::Seq::fold_struct_rnafold( $record->{ "SEQ" } );
2509                 $record->{ "SEC_STRUCT" }  = $struct;
2510                 $record->{ "FREE_ENERGY" } = $index;
2511                 $record->{ "SCORE" }       = abs int $index;
2512                 $record->{ "SIZE" }        = length $struct;
2513                 $record->{ "CONF" }        = "1," x $record->{ "SIZE" };
2514             }
2515         }
2516
2517         Maasha::Biopieces::put_record( $record, $out );
2518     }
2519 }
2520
2521
2522 sub script_split_seq
2523 {
2524     # Martin A. Hansen, August 2007.
2525
2526     # Split a sequence in stream into words.
2527
2528     my ( $in,        # handle to in stream
2529          $out,       # handle to out stream
2530          $options,   # options hash
2531        ) = @_;
2532
2533     # Returns nothing.
2534
2535     my ( $record, $new_record, $i, $subseq, %lookup );
2536
2537     $options->{ "word_size" } ||= 7;
2538
2539     while ( $record = Maasha::Biopieces::get_record( $in ) )
2540     {
2541         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } )
2542         {
2543             for ( $i = 0; $i < length( $record->{ "SEQ" } ) - $options->{ "word_size" } + 1; $i++ )
2544             {
2545                 $subseq = substr $record->{ "SEQ" }, $i, $options->{ "word_size" };
2546
2547                 if ( $options->{ "uniq" } and not $lookup{ $subseq } )
2548                 {
2549                     $new_record->{ "SEQ_NAME" } = $record->{ "SEQ_NAME" } . "[" . ( $i + 1 ) . "-" . ( $i + $options->{ "word_size" } ) . "]";
2550                     $new_record->{ "SEQ" }      = $subseq;
2551
2552                     Maasha::Biopieces::put_record( $new_record, $out );
2553
2554                     $lookup{ $subseq } = 1;
2555                 }
2556                 else
2557                 {
2558                     $new_record->{ "SEQ_NAME" } = $record->{ "SEQ_NAME" } . "[" . ( $i + 1 ) . "-" . ( $i + $options->{ "word_size" } ) . "]";
2559                     $new_record->{ "SEQ" }      = $subseq;
2560
2561                     Maasha::Biopieces::put_record( $new_record, $out );
2562                 }
2563             }
2564         }
2565         else
2566         {
2567             Maasha::Biopieces::put_record( $record, $out );
2568         }
2569     }
2570 }
2571
2572
2573 sub script_split_bed
2574 {
2575     # Martin A. Hansen, June 2008.
2576
2577     # Split a BED record into overlapping windows.
2578
2579     my ( $in,        # handle to in stream
2580          $out,       # handle to out stream
2581          $options,   # options hash
2582        ) = @_;
2583
2584     # Returns nothing.
2585
2586     my ( $record, $new_record, $i );
2587
2588     $options->{ "window_size" } ||= 20;
2589     $options->{ "step_size" }   ||= 1;
2590
2591     while ( $record = Maasha::Biopieces::get_record( $in ) )
2592     {
2593         if ( $record->{ "CHR" } and $record->{ "CHR_BEG" } and $record->{ "CHR_END" } )
2594         {
2595             $record->{ "BED_LEN" } = $record->{ "CHR_END" } - $record->{ "CHR_BEG" } + 1;
2596
2597             for ( $i = 0; $i < $record->{ "BED_LEN" } - $options->{ "window_size" }; $i += $options->{ "step_size" } )
2598             {
2599                 $new_record->{ "REC_TYPE" } = "BED";
2600                 $new_record->{ "CHR" }      = $record->{ "CHR" };
2601                 $new_record->{ "CHR_BEG" }  = $record->{ "CHR_BEG" } + $i;
2602                 $new_record->{ "CHR_END" }  = $record->{ "CHR_BEG" } + $i + $options->{ "window_size" };
2603                 $new_record->{ "BED_LEN" }  = $options->{ "window_size" };
2604                 $new_record->{ "Q_ID" }     = $record->{ "Q_ID" } . "_$i";
2605                 $new_record->{ "SCORE" }    = $record->{ "SCORE" };
2606                 $new_record->{ "STRAND" }   = $record->{ "STRAND" };
2607
2608                 Maasha::Biopieces::put_record( $new_record, $out );
2609             }
2610         }
2611         else
2612         {
2613             Maasha::Biopieces::put_record( $record, $out );
2614         }
2615     }
2616 }
2617
2618
2619 sub script_align_seq
2620 {
2621     # Martin A. Hansen, August 2007.
2622
2623     # Align sequences in stream.
2624
2625     my ( $in,    # handle to in stream
2626          $out,   # handle to out stream
2627        ) = @_;
2628
2629     # Returns nothing.
2630
2631     my ( $record, @entries, $entry );
2632
2633     while ( $record = Maasha::Biopieces::get_record( $in ) )
2634     {
2635         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } ) {
2636             push @entries, [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ];
2637         } elsif ( $record->{ "Q_ID" } and $record->{ "SEQ" } ) {
2638             push @entries, [ $record->{ "Q_ID" }, $record->{ "SEQ" } ];
2639         } else {
2640             Maasha::Biopieces::put_record( $record, $out );
2641         }
2642     }
2643
2644     @entries = Maasha::Align::align( \@entries );
2645
2646     foreach $entry ( @entries )
2647     {
2648         if ( $entry->[ SEQ_NAME ] and $entry->[ SEQ ] )
2649         {
2650             $record = {
2651                 SEQ_NAME => $entry->[ SEQ_NAME ],
2652                 SEQ      => $entry->[ SEQ ],
2653             };
2654
2655             Maasha::Biopieces::put_record( $record, $out );
2656         }
2657     }
2658 }
2659
2660
2661 sub script_tile_seq
2662 {
2663     # Martin A. Hansen, February 2008.
2664
2665     # Using the first sequence in stream as reference, tile
2666     # all subsequent sequences based on pairwise alignments.
2667
2668     my ( $in,        # handle to in stream
2669          $out,       # handle to out stream
2670          $options,   # options hash
2671        ) = @_;
2672
2673     # Returns nothing.
2674
2675     my ( $record, $first, $ref_entry, @entries );
2676
2677     $first = 1;
2678
2679     while ( $record = Maasha::Biopieces::get_record( $in ) )
2680     {
2681         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } )
2682         {
2683             if ( $first )
2684             {
2685                 $ref_entry = [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ];
2686
2687                 $first = 0;
2688             }
2689             else
2690             {
2691                 push @entries, [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ];
2692             }
2693         }
2694         else
2695         {
2696             Maasha::Biopieces::put_record( $record, $out );
2697         }
2698     }
2699
2700     @entries = Maasha::Align::align_tile( $ref_entry, \@entries, $options );
2701
2702     map { Maasha::Biopieces::put_record( { SEQ_NAME => $_->[ SEQ_NAME ], SEQ => $_->[ SEQ ] }, $out ) } @entries;
2703 }
2704
2705
2706 sub script_invert_align
2707 {
2708     # Martin A. Hansen, February 2008.
2709
2710     # Inverts an alignment showing only non-mathing residues
2711     # using the first sequence as reference.
2712
2713     my ( $in,        # handle to in stream
2714          $out,       # handle to out stream
2715          $options,   # options hash
2716        ) = @_;
2717
2718     # Returns nothing.
2719
2720     my ( $record, @entries );
2721
2722     while ( $record = Maasha::Biopieces::get_record( $in ) )
2723     {
2724         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } )
2725         {
2726             push @entries, [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ];
2727         }
2728         else
2729         {
2730             Maasha::Biopieces::put_record( $record, $out );
2731         }
2732     }
2733
2734     Maasha::Align::align_invert( \@entries, $options->{ "soft" } );
2735
2736     map { Maasha::Biopieces::put_record( { SEQ_NAME => $_->[ SEQ_NAME ], SEQ => $_->[ SEQ ] }, $out ) } @entries;
2737 }
2738
2739
2740 sub script_patscan_seq
2741 {
2742     # Martin A. Hansen, August 2007.
2743
2744     # Locates patterns in sequences using scan_for_matches.
2745
2746     my ( $in,        # handle to in stream
2747          $out,       # handle to out stream
2748          $options,   # options hash
2749        ) = @_;
2750
2751     # Returns nothing.
2752
2753     my ( $genome_file, @args, $arg, $type, $seq_file, $pat_file, $out_file, $fh_in, $fh_out, $record, $patterns, $pattern, $entry, $result, %head_hash, $i );
2754
2755     if ( $options->{ "patterns" } ) {
2756         $patterns = Maasha::Patscan::parse_patterns( $options->{ "patterns" } );
2757     } elsif ( -f $options->{ "patterns_in" } ) {
2758         $patterns = Maasha::Patscan::read_patterns( $options->{ "patterns_in" } );
2759     }
2760
2761     $genome_file = "$ENV{ 'BP_DATA' }/genomes/$options->{ 'genome' }/fasta/$options->{ 'genome' }.fna" if $options->{ 'genome' };
2762
2763     push @args, "-c"                            if $options->{ "comp" };
2764     push @args, "-m $options->{ 'max_hits' }"   if $options->{ 'max_hits' };
2765     push @args, "-n $options->{ 'max_misses' }" if $options->{ 'max_hits' };
2766
2767     $seq_file = "$BP_TMP/patscan.seq";
2768     $pat_file = "$BP_TMP/patscan.pat";
2769     $out_file = "$BP_TMP/patscan.out";
2770
2771     $fh_out = Maasha::Common::write_open( $seq_file );
2772
2773     $i = 0;
2774
2775     while ( $record = Maasha::Biopieces::get_record( $in ) )
2776     {
2777         if ( $record->{ "SEQ" } and $record->{ "SEQ_NAME" } )
2778         {
2779             $type = Maasha::Seq::seq_guess_type( $record->{ "SEQ" } ) if not $type;
2780
2781             Maasha::Fasta::put_entry( [ $i, $record->{ "SEQ" } ], $fh_out );
2782
2783             $head_hash{ $i } = $record->{ "SEQ_NAME" };
2784
2785             $i++;
2786         }
2787     }
2788
2789     close $fh_out;
2790
2791     $arg  = join " ", @args;
2792     $arg .= " -p" if $type eq "protein";
2793
2794     foreach $pattern ( @{ $patterns } )
2795     {
2796         $fh_out = Maasha::Common::write_open( $pat_file );
2797
2798         print $fh_out "$pattern\n";
2799
2800         close $fh_out;
2801
2802         if ( $options->{ 'genome' } ) {
2803             `scan_for_matches $arg $pat_file < $genome_file > $out_file`;
2804             # Maasha::Common::run( "scan_for_matches", "$arg $pat_file < $genome_file > $out_file" );
2805         } else {
2806             `scan_for_matches $arg $pat_file < $seq_file > $out_file`;
2807             # Maasha::Common::run( "scan_for_matches", "$arg $pat_file < $seq_file > $out_file" );
2808         }
2809
2810         $fh_in = Maasha::Common::read_open( $out_file );
2811
2812         while ( $entry = Maasha::Fasta::get_entry( $fh_in ) )
2813         {
2814             $result = Maasha::Patscan::parse_scan_result( $entry, $pattern );
2815
2816             if ( $options->{ 'genome' } )
2817             {
2818                 $result->{ "CHR" }     = $result->{ "S_ID" };
2819                 $result->{ "CHR_BEG" } = $result->{ "S_BEG" };
2820                 $result->{ "CHR_END" } = $result->{ "S_END" };
2821
2822                 delete $result->{ "S_ID" };
2823                 delete $result->{ "S_BEG" };
2824                 delete $result->{ "S_END" };
2825             }
2826             else
2827             {
2828                 $result->{ "S_ID" } = $head_hash{ $result->{ "S_ID" } };
2829             }
2830
2831             Maasha::Biopieces::put_record( $result, $out );
2832         }
2833
2834         close $fh_in;
2835     }
2836
2837     unlink $pat_file;
2838     unlink $seq_file;
2839     unlink $out_file;
2840 }
2841
2842
2843 sub script_blat_seq
2844 {
2845     # Martin A. Hansen, August 2007.
2846
2847     # BLATs sequences in stream against a given genome.
2848
2849     my ( $in,        # handle to in stream
2850          $out,       # handle to out stream
2851          $options,   # options hash
2852        ) = @_;
2853
2854     # Returns nothing.
2855
2856     my ( $blat_args, $genome_file, $query_file, $fh_in, $fh_out, $type, $record, $result_file, $entries, $entry );
2857
2858     $genome_file = "$ENV{ 'BP_DATA' }/genomes/$options->{ 'genome' }/fasta/$options->{ 'genome' }.fna";
2859
2860     $options->{ 'tile_size' }    ||= 11;
2861     $options->{ 'one_off' }      ||= 0;
2862     $options->{ 'min_identity' } ||= 90;
2863     $options->{ 'min_score' }    ||= 0;
2864     $options->{ 'step_size' }    ||= $options->{ 'tile_size' };
2865
2866     $blat_args .= " -tileSize=$options->{ 'tile_size' }";
2867     $blat_args .= " -oneOff=$options->{ 'one_off' }";
2868     $blat_args .= " -minIdentity=$options->{ 'min_identity' }";
2869     $blat_args .= " -minScore=$options->{ 'min_score' }";
2870     $blat_args .= " -stepSize=$options->{ 'step_size' }";
2871 #    $blat_args .= " -ooc=" . Maasha::Config::genome_blat_ooc( $options->{ "genome" }, 11 ) if $options->{ 'ooc' };
2872
2873     $query_file = "$BP_TMP/blat.seq";
2874
2875     $fh_out = Maasha::Common::write_open( $query_file );
2876
2877     while ( $record = Maasha::Biopieces::get_record( $in ) )
2878     {
2879         if ( $entry = Maasha::Fasta::biopiece2fasta( $record ) )
2880         {
2881             $type = Maasha::Seq::seq_guess_type( $entry->[ SEQ ] ) if not $type;
2882             Maasha::Fasta::put_entry( $entry, $fh_out, 80 );
2883         }
2884
2885         Maasha::Biopieces::put_record( $record, $out );
2886     }
2887
2888     close $fh_out;
2889
2890     $blat_args .= " -t=dnax" if $type eq "protein";
2891     $blat_args .= " -q=$type";
2892
2893     $result_file = "$BP_TMP/blat.psl";
2894
2895     Maasha::Common::run( "blat", "$genome_file $query_file $blat_args $result_file > /dev/null 2>&1" );
2896
2897     unlink $query_file;
2898
2899     $entries = Maasha::UCSC::psl_get_entries( $result_file );
2900
2901     map { Maasha::Biopieces::put_record( $_, $out ) } @{ $entries };
2902
2903     unlink $result_file;
2904 }
2905
2906
2907 sub script_soap_seq
2908 {
2909     # Martin A. Hansen, July 2008.
2910
2911     # soap sequences in stream against a given file or genome.
2912
2913     my ( $in,        # handle to in stream
2914          $out,       # handle to out stream
2915          $options,   # options hash
2916        ) = @_;
2917
2918     # Returns nothing.
2919
2920     my ( $genome, $tmp_in, $tmp_out, $fh_in, $fh_out, $record, $line, @fields, $entry, $count, $args );
2921
2922     $options->{ "seed_size" }  ||= 10;
2923     $options->{ "mismatches" } ||= 2;
2924     $options->{ "gap_size" }   ||= 0;
2925     $options->{ "cpus" }       ||= 1;
2926
2927     if ( $options->{ "genome" } ) {
2928         $options->{ "in_file" } = "$ENV{ 'BP_DATA' }/genomes/$options->{ 'genome' }/fasta/$options->{ 'genome' }.fna";
2929     }
2930
2931     $tmp_in  = "$BP_TMP/soap_query.seq";
2932     $tmp_out = "$BP_TMP/soap.result";
2933
2934     $fh_out = Maasha::Common::write_open( $tmp_in );
2935
2936     $count = 0;
2937
2938     while ( $record = Maasha::Biopieces::get_record( $in ) )
2939     {
2940         if ( $entry = Maasha::Fasta::biopiece2fasta( $record ) )
2941         {
2942             Maasha::Fasta::put_entry( $entry, $fh_out );
2943
2944             $count++;
2945         }
2946
2947         Maasha::Biopieces::put_record( $record, $out );
2948     }
2949
2950     close $fh_out;
2951
2952     if ( $count > 0 )
2953     {
2954         $args = join( " ",
2955             "-s $options->{ 'seed_size' }",
2956             "-r 2",
2957             "-a $tmp_in",
2958             "-v $options->{ 'mismatches' }",
2959             "-g $options->{ 'gap_size' }",
2960             "-p $options->{ 'cpus' }",
2961             "-d $options->{ 'in_file' }",
2962             "-o $tmp_out",
2963         );
2964
2965         $args .= " > /dev/null 2>&1" if not $options->{ 'verbose' };
2966
2967         Maasha::Common::run( "soap", $args, 1 );
2968
2969         unlink $tmp_in;
2970
2971         $fh_out = Maasha::Common::read_open( $tmp_out );
2972
2973         undef $record;
2974
2975         while ( $line = <$fh_out> )
2976         {
2977             chomp $line;
2978
2979             @fields = split /\t/, $line;
2980
2981             $record->{ "REC_TYPE" }   = "SOAP";
2982             $record->{ "Q_ID" }       = $fields[ 0 ];
2983             $record->{ "SCORE" }      = $fields[ 3 ];
2984             $record->{ "STRAND" }     = $fields[ 6 ];
2985             $record->{ "S_ID" }       = $fields[ 7 ];
2986             $record->{ "S_BEG" }      = $fields[ 8 ] - 1; # soap is 1-based
2987             $record->{ "S_END" }      = $fields[ 8 ] + $fields[ 5 ] - 2;
2988
2989             Maasha::Biopieces::put_record( $record, $out );
2990         }
2991
2992         close $fh_out;
2993     }
2994
2995     unlink $tmp_out;
2996 }
2997
2998
2999 sub script_match_seq
3000 {
3001     # Martin A. Hansen, August 2007.
3002
3003     # BLATs sequences in stream against a given genome.
3004
3005     my ( $in,        # handle to in stream
3006          $out,       # handle to out stream
3007          $options,   # options hash
3008        ) = @_;
3009
3010     # Returns nothing.
3011
3012     my ( $record, @entries, $results );
3013
3014     $options->{ "word_size" } ||= 20;
3015     $options->{ "direction" } ||= "both";
3016
3017     while ( $record = Maasha::Biopieces::get_record( $in ) )
3018     {
3019         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } ) {
3020             push @entries, [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ];
3021         }
3022
3023         Maasha::Biopieces::put_record( $record, $out );
3024     }
3025
3026     if ( @entries == 1 )
3027     {
3028         $results = Maasha::Match::match_mummer( [ $entries[ 0 ] ], [ $entries[ 0 ] ], $options, $BP_TMP );
3029
3030         map { Maasha::Biopieces::put_record( $_, $out ) } @{ $results };
3031     }
3032     elsif ( @entries == 2 )
3033     {
3034         $results = Maasha::Match::match_mummer( [ $entries[ 0 ] ], [ $entries[ 1 ] ], $options, $BP_TMP );
3035
3036         map { Maasha::Biopieces::put_record( $_, $out ) } @{ $results };
3037     }
3038 }
3039
3040
3041 sub script_create_vmatch_index
3042 {
3043     # Martin A. Hansen, January 2008.
3044
3045     # Create a vmatch index from sequences in the stream.
3046
3047     my ( $in,        # handle to in stream
3048          $out,       # handle to out stream
3049          $options,   # options hash
3050        ) = @_;
3051
3052     # Returns nothing.
3053
3054     my ( $record, $file_tmp, $fh_tmp, $type, $entry );
3055
3056     if ( $options->{ "index_name" } )
3057     {
3058         $file_tmp = $options->{ 'index_name' };
3059         $fh_tmp   = Maasha::Common::write_open( $file_tmp );
3060     }
3061
3062     while ( $record = Maasha::Biopieces::get_record( $in ) )
3063     {
3064         if ( $options->{ "index_name" } and $entry = Maasha::Fasta::biopiece2fasta( $record ) )
3065         {
3066             Maasha::Fasta::put_entry( $entry, $fh_tmp );
3067
3068             $type = Maasha::Seq::seq_guess_type( $entry->[ SEQ ] ) if not defined $type;
3069         }
3070
3071         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
3072     }
3073
3074     if ( $options->{ "index_name" } )
3075     {
3076         close $fh_tmp;
3077
3078         if ( $type eq "protein" ) {
3079             Maasha::Common::run( "mkvtree", "-db $file_tmp -protein -pl $options->{ 'prefix_length' } -allout -indexname $file_tmp > /dev/null 2>&1" );
3080         } else {
3081             Maasha::Common::run( "mkvtree", "-db $file_tmp -dna -pl $options->{ 'prefix_length' } -allout -indexname $file_tmp > /dev/null 2>&1" );
3082         }
3083
3084         unlink $file_tmp;
3085     }
3086 }
3087
3088
3089 sub script_write_align
3090 {
3091     # Martin A. Hansen, August 2007.
3092
3093     # Write pretty alignments aligned sequences in stream.
3094
3095     my ( $in,        # handle to in stream
3096          $out,       # handle to out stream
3097          $options,   # options hash
3098        ) = @_;
3099
3100     # Returns nothing.
3101
3102     my ( $fh, $record, @entries );
3103
3104     $fh = write_stream( $options->{ "data_out" } ) ;
3105
3106     while ( $record = Maasha::Biopieces::get_record( $in ) )
3107     {
3108         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } ) {
3109             push @entries, [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ];
3110         }
3111
3112         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
3113     }
3114
3115     if ( scalar( @entries ) == 2 ) {
3116         Maasha::Align::align_print_pairwise( $entries[ 0 ], $entries[ 1 ], $fh, $options->{ "wrap" } );
3117     } elsif ( scalar ( @entries ) > 2 ) {
3118         Maasha::Align::align_print_multi( \@entries, $fh, $options->{ "wrap" }, $options->{ "no_ruler" }, $options->{ "no_consensus" } );
3119     }
3120
3121     close $fh if $fh;
3122 }
3123
3124
3125 sub script_write_bed
3126 {
3127     # Martin A. Hansen, August 2007.
3128
3129     # Write BED format for the UCSC genome browser using records in stream.
3130
3131     my ( $in,        # handle to in stream
3132          $out,       # handle to out stream
3133          $options,   # options hash
3134        ) = @_;
3135
3136     # Returns nothing.
3137
3138     my ( $cols, $fh, $record, $bed_entry, $new_record );
3139
3140     $cols = $options->{ 'cols' }->[ 0 ];
3141
3142     $fh = Maasha::Biopieces::write_stream( $options->{ 'data_out' }, $options->{ 'compress' } );
3143
3144     while ( $record = Maasha::Biopieces::get_record( $in ) )
3145     {
3146         $record = Maasha::UCSC::psl2record( $record ) if $record->{ 'tBaseInsert' }; # Dirty addition to allow Affy data from MySQL to be dumped
3147
3148         if ( $bed_entry = Maasha::UCSC::BED::biopiece2bed( $record, $cols ) ) {
3149             Maasha::UCSC::BED::bed_entry_put( $bed_entry, $fh, $cols, $options->{ 'check' } );
3150         }
3151
3152         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ 'no_stream' };
3153     }
3154
3155     close $fh;
3156 }
3157
3158
3159 sub script_write_psl
3160 {
3161     # Martin A. Hansen, August 2007.
3162
3163     # Write PSL output from stream.
3164
3165     my ( $in,        # handle to in stream
3166          $out,       # handle to out stream
3167          $options,   # options hash
3168        ) = @_;
3169
3170     # Returns nothing.
3171
3172     my ( $fh, $record, @output, $first );
3173
3174     $first = 1;
3175
3176     $fh = write_stream( $options->{ "data_out" }, $options->{ "compress" } );
3177
3178     while ( $record = Maasha::Biopieces::get_record( $in ) )
3179     {
3180         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
3181
3182         if ( $record->{ "REC_TYPE" } and $record->{ "REC_TYPE" } eq "PSL" )
3183         {
3184             Maasha::UCSC::psl_put_header( $fh ) if $first;
3185             Maasha::UCSC::psl_put_entry( $record, $fh );
3186             $first = 0;
3187         }
3188     }
3189
3190     close $fh;
3191 }
3192
3193
3194 sub script_write_fixedstep
3195 {
3196     # Martin A. Hansen, Juli 2008.
3197
3198     # Write fixedStep entries from recrods in the stream.
3199
3200     my ( $in,        # handle to in stream
3201          $out,       # handle to out stream
3202          $options,   # options hash
3203        ) = @_;
3204
3205     # Returns nothing.
3206
3207     my ( $fh, $record, $entry );
3208
3209     $fh = write_stream( $options->{ "data_out" }, $options->{ "compress" } );
3210
3211     while ( $record = Maasha::Biopieces::get_record( $in ) )
3212     {
3213         if ( $entry = Maasha::UCSC::Wiggle::biopiece2fixedstep( $record ) ) {
3214             Maasha::UCSC::Wiggle::fixedstep_entry_put( $entry, $fh );
3215         }
3216
3217         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
3218     }
3219
3220     close $fh;
3221 }
3222
3223
3224 sub script_write_2bit
3225 {
3226     # Martin A. Hansen, March 2008.
3227
3228     # Write sequence entries from stream in 2bit format.
3229
3230     my ( $in,        # handle to in stream
3231          $out,       # handle to out stream
3232          $options,   # options hash
3233        ) = @_;
3234
3235     # Returns nothing.
3236
3237     my ( $record, $mask, $tmp_file, $fh_tmp, $fh_in, $fh_out, $entry );
3238
3239     $mask = 1 if not $options->{ "no_mask" };
3240
3241     $tmp_file = "$BP_TMP/write_2bit.fna";
3242     $fh_tmp   = Maasha::Common::write_open( $tmp_file );
3243
3244     $fh_out = write_stream( $options->{ "data_out" } );
3245
3246     while ( $record = Maasha::Biopieces::get_record( $in ) )
3247     {
3248         if ( $entry = Maasha::Fasta::biopiece2fasta( $record ) ) {
3249             Maasha::Fasta::put_entry( $entry, $fh_tmp );
3250         }
3251
3252         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
3253     }
3254
3255     close $fh_tmp;
3256
3257     $fh_in = Maasha::Common::read_open( $tmp_file );
3258
3259     Maasha::TwoBit::fasta2twobit( $fh_in, $fh_out, $mask );
3260
3261     close $fh_in;
3262     close $fh_out;
3263
3264     unlink $tmp_file;
3265 }
3266
3267
3268 sub script_write_solid
3269 {
3270     # Martin A. Hansen, April 2008.
3271
3272     # Write di-base encoded Solid sequence from entries in stream.
3273
3274     my ( $in,        # handle to in stream
3275          $out,       # handle to out stream
3276          $options,   # options hash
3277        ) = @_;
3278
3279     # Returns nothing.
3280
3281     my ( $record, $fh, $entry );
3282
3283     $fh = write_stream( $options->{ "data_out" }, $options->{ "compress" } );
3284
3285     while ( $record = Maasha::Biopieces::get_record( $in ) )
3286     {
3287         if ( $entry = Maasha::Fasta::biopiece2fasta( $record ) )
3288         {
3289             $entry->[ SEQ ] = Maasha::Solid::seq2color_space( uc $entry->[ SEQ ] );
3290
3291             Maasha::Fasta::put_entry( $entry, $fh, $options->{ "wrap" } );
3292         }
3293
3294         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
3295     }
3296
3297     close $fh;
3298 }
3299
3300
3301 sub script_write_ucsc_config
3302 {
3303     # Martin A. Hansen, November 2008.
3304
3305     # Write UCSC Genome Broser configuration (.ra file type) from
3306     # records in the stream.
3307
3308     my ( $in,        # handle to in stream
3309          $out,       # handle to out stream
3310          $options,   # options hash
3311        ) = @_;
3312
3313     # Returns nothing.
3314
3315     my ( $record, $fh );
3316
3317     $fh = write_stream( $options->{ "data_out" } );
3318
3319     while ( $record = Maasha::Biopieces::get_record( $in ) )
3320     {
3321         Maasha::UCSC::ucsc_config_entry_put( $record, $fh ) if $record->{ "REC_TYPE" } eq "UCSC Config";
3322
3323         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
3324     }
3325
3326     close $fh;
3327 }
3328
3329
3330 sub script_plot_seqlogo
3331 {
3332     # Martin A. Hansen, August 2007.
3333
3334     # Calculates and writes a sequence logo for alignments.
3335
3336     my ( $in,        # handle to in stream
3337          $out,       # handle to out stream
3338          $options,   # options hash
3339        ) = @_;
3340
3341     # Returns nothing.
3342
3343     my ( $record, @entries, $logo, $fh );
3344
3345     while ( $record = Maasha::Biopieces::get_record( $in ) )
3346     {
3347         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } ) {
3348             push @entries, [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ];
3349         }
3350
3351         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
3352     }
3353
3354     $logo = Maasha::Plot::seq_logo( \@entries );
3355
3356     $fh = write_stream( $options->{ "data_out" } );
3357
3358     print $fh $logo;
3359
3360     close $fh;
3361 }
3362
3363
3364 sub script_plot_phastcons_profiles
3365 {
3366     # Martin A. Hansen, January 2008.
3367
3368     # Plots PhastCons profiles.
3369
3370     my ( $in,        # handle to in stream
3371          $out,       # handle to out stream
3372          $options,   # options hash
3373        ) = @_;
3374
3375     # Returns nothing.
3376
3377     my ( $phastcons_file, $phastcons_index, $index, $fh_phastcons, $record, $scores, $AoA, $plot, $fh );
3378
3379     $options->{ "title" } ||= "PhastCons Profiles";
3380
3381     $phastcons_file  = Maasha::Config::genome_phastcons( $options->{ "genome" } );
3382     $phastcons_index = Maasha::Config::genome_phastcons_index( $options->{ "genome" } );
3383
3384     $index           = Maasha::UCSC::fixedstep_index_retrieve( $phastcons_index );
3385     $fh_phastcons    = Maasha::Common::read_open( $phastcons_file );
3386
3387     while ( $record = Maasha::Biopieces::get_record( $in ) )
3388     {
3389         if ( $record->{ "CHR" } and $record->{ "CHR_BEG" } and $record->{ "CHR_END" } )
3390         {
3391             $scores = Maasha::UCSC::fixedstep_index_lookup( $index, $fh_phastcons, $record->{ "CHR" },
3392                                                                                    $record->{ "CHR_BEG" },
3393                                                                                    $record->{ "CHR_END" },
3394                                                                                    $options->{ "flank" } );
3395
3396             push @{ $AoA }, [ @{ $scores } ];
3397         }
3398
3399         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
3400     }
3401
3402     Maasha::UCSC::phastcons_normalize( $AoA );
3403
3404     $AoA = [ [ Maasha::UCSC::phastcons_mean( $AoA ) ] ] if $options->{ "mean" };
3405     $AoA = [ [ Maasha::UCSC::phastcons_median( $AoA ) ] ] if $options->{ "median" };
3406
3407     $AoA = Maasha::Matrix::matrix_flip( $AoA );
3408
3409     $plot = Maasha::Plot::lineplot_simple( $AoA, $options, $BP_TMP );
3410
3411     $fh = write_stream( $options->{ "data_out" } );
3412
3413     print $fh "$_\n" foreach @{ $plot };
3414
3415     close $fh;
3416 }
3417
3418
3419 sub script_analyze_bed
3420 {
3421     # Martin A. Hansen, March 2008.
3422
3423     # Analyze BED entries in stream.
3424
3425     my ( $in,        # handle to in stream
3426          $out,       # handle to out stream
3427          $options,   # options hash
3428        ) = @_;
3429
3430     # Returns nothing.
3431
3432     my ( $record );
3433
3434     while ( $record = Maasha::Biopieces::get_record( $in ) )
3435     {
3436         $record = Maasha::UCSC::bed_analyze( $record ) if $record->{ "REC_TYPE" } eq "BED";
3437
3438         Maasha::Biopieces::put_record( $record, $out );
3439     }
3440 }
3441
3442
3443 sub script_analyze_vals
3444 {
3445     # Martin A. Hansen, August 2007.
3446
3447     # Analyze values for given keys in stream.
3448
3449     my ( $in,        # handle to in stream
3450          $out,       # handle to out stream
3451          $options,   # options hash
3452        ) = @_;
3453
3454     # Returns nothing.
3455
3456     my ( $record, $key, @keys, %key_hash, $analysis, $len );
3457
3458     map { $key_hash{ $_ } = 1 } @{ $options->{ "keys" } };
3459
3460     while ( $record = Maasha::Biopieces::get_record( $in ) )
3461     {
3462         foreach $key ( keys %{ $record } )
3463         {
3464             next if $options->{ "keys" } and not exists $key_hash{ $key };
3465
3466             $analysis->{ $key }->{ "COUNT" }++;
3467
3468             if ( Maasha::Calc::is_a_number( $record->{ $key } ) )
3469             {
3470                 $analysis->{ $key }->{ "TYPE" } = "num";
3471                 $analysis->{ $key }->{ "SUM" } += $record->{ $key };
3472                 $analysis->{ $key }->{ "MAX" } = $record->{ $key } if $record->{ $key } > $analysis->{ $key }->{ "MAX" } or not $analysis->{ $key }->{ "MAX" };
3473                 $analysis->{ $key }->{ "MIN" } = $record->{ $key } if $record->{ $key } < $analysis->{ $key }->{ "MIN" } or not $analysis->{ $key }->{ "MIN" };
3474             }
3475             else
3476             {
3477                 $len = length $record->{ $key };
3478
3479                 $analysis->{ $key }->{ "TYPE" } = "alph";
3480                 $analysis->{ $key }->{ "SUM" } += $len;
3481                 $analysis->{ $key }->{ "MAX" } = $len if $len > $analysis->{ $key }->{ "MAX" } or not $analysis->{ $key }->{ "MAX" };
3482                 $analysis->{ $key }->{ "MIN" } = $len if $len < $analysis->{ $key }->{ "MIM" } or not $analysis->{ $key }->{ "MIN" };
3483             }
3484         }
3485
3486         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
3487     }
3488
3489     foreach $key ( keys %{ $analysis } )
3490     {
3491         $analysis->{ $key }->{ "MEAN" } = sprintf "%.2f", $analysis->{ $key }->{ "SUM" } / $analysis->{ $key }->{ "COUNT" };
3492         $analysis->{ $key }->{ "SUM" }  = sprintf "%.2f", $analysis->{ $key }->{ "SUM" };
3493     }
3494
3495     my ( $keys, $types, $counts, $mins, $maxs, $sums, $means );
3496
3497     $keys   = "KEY  ";
3498     $types  = "TYPE ";
3499     $counts = "COUNT";
3500     $mins   = "MIN  ";
3501     $maxs   = "MAX  ";
3502     $sums   = "SUM  ";
3503     $means  = "MEAN ";
3504
3505     if ( $options->{ "keys" } ) {
3506         @keys = @{ $options->{ "keys" } };
3507     } else {
3508         @keys = keys %{ $analysis };
3509     }
3510
3511     foreach $key ( @keys )
3512     {
3513         $keys   .= sprintf "% 15s", $key;
3514         $types  .= sprintf "% 15s", $analysis->{ $key }->{ "TYPE" };
3515         $counts .= sprintf "% 15s", $analysis->{ $key }->{ "COUNT" };
3516         $mins   .= sprintf "% 15s", $analysis->{ $key }->{ "MIN" };
3517         $maxs   .= sprintf "% 15s", $analysis->{ $key }->{ "MAX" };
3518         $sums   .= sprintf "% 15s", $analysis->{ $key }->{ "SUM" };
3519         $means  .= sprintf "% 15s", $analysis->{ $key }->{ "MEAN" };
3520     }
3521
3522     print $out "$keys\n";
3523     print $out "$types\n";
3524     print $out "$counts\n";
3525     print $out "$mins\n";
3526     print $out "$maxs\n";
3527     print $out "$sums\n";
3528     print $out "$means\n";
3529 }
3530
3531
3532 sub script_head_records
3533 {
3534     # Martin A. Hansen, August 2007.
3535
3536     # Display the first sequences in stream.
3537
3538     my ( $in,        # handle to in stream
3539          $out,       # handle to out stream
3540          $options,   # options hash
3541        ) = @_;
3542
3543     # Returns nothing.
3544
3545     my ( $record, $count );
3546
3547     $options->{ "num" } ||= 10;
3548
3549     $count = 0;
3550
3551     while ( $record = Maasha::Biopieces::get_record( $in ) )
3552     {
3553         $count++;
3554
3555         Maasha::Biopieces::put_record( $record, $out );
3556
3557         last if $count == $options->{ "num" };
3558     }
3559 }
3560
3561
3562 sub script_remove_keys
3563 {
3564     # Martin A. Hansen, August 2007.
3565
3566     # Remove keys from stream.
3567
3568     my ( $in,        # handle to in stream
3569          $out,       # handle to out stream
3570          $options,   # options hash
3571        ) = @_;
3572
3573     # Returns nothing.
3574
3575     my ( $record, $new_record );
3576
3577     while ( $record = Maasha::Biopieces::get_record( $in ) )
3578     {
3579         if ( $options->{ "keys" } )
3580         {
3581             map { delete $record->{ $_ } } @{ $options->{ "keys" } };
3582         }
3583         elsif ( $options->{ "save_keys" } )
3584         {
3585             map { $new_record->{ $_ } = $record->{ $_ } if exists $record->{ $_ } } @{ $options->{ "save_keys" } };
3586
3587             $record = $new_record;
3588         }
3589
3590         Maasha::Biopieces::put_record( $record, $out ) if keys %{ $record };
3591     }
3592 }
3593
3594
3595 sub script_remove_adaptor
3596 {
3597     # Martin A. Hansen, August 2008.
3598
3599     # Find and remove adaptor from sequences in the stream.
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 ( $record, $adaptor, $seq, $adaptor_len, $seq_len, $offset, $max_match, $max_mismatch, $pos );
3609
3610     $options->{ "remove" } ||= "after";
3611
3612     $max_mismatch = $options->{ "mismatches" } || 0;
3613     $offset       = $options->{ "offset" };
3614
3615     if ( not defined $offset ) {
3616         $offset = 0;
3617     } else {
3618         $offset--;
3619     }
3620
3621     $adaptor      = uc $options->{ "adaptor" };
3622     $adaptor_len  = length $adaptor;
3623
3624     while ( $record = Maasha::Biopieces::get_record( $in ) )
3625     {
3626         if ( $record->{ "SEQ" } )
3627         {
3628             $seq     = uc $record->{ "SEQ" };
3629             $seq_len = length $seq;
3630
3631             $pos = Maasha::Common::index_m( $seq, $adaptor, $seq_len, $adaptor_len, $offset, $max_mismatch );
3632
3633             $record->{ "ADAPTOR_POS" } = $pos;
3634
3635             if ( $pos >= 0 and $options->{ "remove" } ne "skip" )
3636             {
3637                 if ( $options->{ "remove" } eq "after" )
3638                 {
3639                     $record->{ "SEQ" }     = substr $record->{ "SEQ" }, 0, $pos;
3640                     $record->{ "SEQ_LEN" } = $pos;
3641                 }
3642                 else
3643                 {
3644                     $record->{ "SEQ" }     = substr $record->{ "SEQ" }, $pos + $adaptor_len;
3645                     $record->{ "SEQ_LEN" } = length $record->{ "SEQ" };
3646                 }
3647             }
3648
3649             Maasha::Biopieces::put_record( $record, $out );
3650         }
3651         else
3652         {
3653             Maasha::Biopieces::put_record( $record, $out );
3654         }
3655     }
3656 }
3657
3658
3659 sub script_remove_mysql_tables
3660 {
3661     # Martin A. Hansen, November 2008.
3662
3663     # Remove MySQL tables from values in stream.
3664
3665     my ( $in,        # handle to in stream
3666          $out,       # handle to out stream
3667          $options,   # options hash
3668        ) = @_;
3669
3670     # Returns nothing.
3671
3672     my ( $record, %table_hash, $dbh, $table );
3673
3674     $options->{ "user" }     ||= Maasha::UCSC::ucsc_get_user();
3675     $options->{ "password" } ||= Maasha::UCSC::ucsc_get_password();
3676
3677     map { $table_hash{ $_ } = 1 } @{ $options->{ 'tables' } };
3678
3679     while ( $record = Maasha::Biopieces::get_record( $in ) )
3680     {
3681         map { $table_hash{ $record->{ $_ } } = 1 } @{ $options->{ 'keys' } };
3682
3683         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ 'no_stream' };
3684     }
3685
3686     $dbh = Maasha::SQL::connect( $options->{ "database" }, $options->{ "user" }, $options->{ "password" } );
3687
3688     foreach $table ( sort keys %table_hash )
3689     {
3690         if ( Maasha::SQL::table_exists( $dbh, $table ) )
3691         {
3692             print STDERR qq(Removing table "$table" from database "$options->{ 'database' }" ... ) if $options->{ 'verbose' };
3693             Maasha::SQL::delete_table( $dbh, $table );
3694             print STDERR "done.\n" if $options->{ 'verbose' };
3695         }
3696         else
3697         {
3698             print STDERR qq(WARNING: table "$table" not found in database "$options->{ 'database' }\n");
3699         }
3700     }
3701
3702     Maasha::SQL::disconnect( $dbh );
3703 }
3704
3705
3706 sub script_remove_ucsc_tracks
3707 {
3708     # Martin A. Hansen, November 2008.
3709
3710     # Remove track from MySQL tables and config file.
3711
3712     my ( $in,        # handle to in stream
3713          $out,       # handle to out stream
3714          $options,   # options hash
3715        ) = @_;
3716
3717     # Returns nothing.
3718
3719     my ( $record, %track_hash, $fh_in, $fh_out, $track, @tracks, @new_tracks, $dbh );
3720
3721     $options->{ 'user' }        ||= Maasha::UCSC::ucsc_get_user();
3722     $options->{ 'password' }    ||= Maasha::UCSC::ucsc_get_password();
3723     $options->{ 'config_file' } ||= "$ENV{ 'HOME' }/ucsc/my_tracks.ra";
3724
3725     map { $track_hash{ $_ } = 1 } @{ $options->{ 'tracks' } };
3726
3727     while ( $record = Maasha::Biopieces::get_record( $in ) )
3728     {
3729         map { $track_hash{ $record->{ $_ } } = 1 } @{ $options->{ 'keys' } };
3730
3731         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ 'no_stream' };
3732     }
3733
3734     $fh_in = Maasha::Filesys::file_read_open( $options->{ 'config_file' } );
3735
3736     while ( $track = Maasha::UCSC::ucsc_config_entry_get( $fh_in ) ) {
3737         push @tracks, $track;
3738     }
3739
3740     close $fh_in;
3741
3742     foreach $track ( @tracks )
3743     {
3744         if ( $track->{ 'database' } eq $options->{ 'database' } and exists $track_hash{ $track->{ 'track' } } ) {
3745             print STDERR qq(Removing track "$track->{ 'track' }" from config file.\n) if $options->{ 'verbose' };
3746         } else {
3747             push @new_tracks, $track;
3748         }
3749     }
3750
3751     rename "$options->{ 'config_file' }", "$options->{ 'config_file' }~";
3752
3753     $fh_out = Maasha::Filesys::file_write_open( $options->{ 'config_file' } );
3754
3755     map { Maasha::UCSC::ucsc_config_entry_put( $_, $fh_out ) } @new_tracks;
3756
3757     close $fh_out;
3758
3759     # ---- locate track in database ----
3760
3761     $dbh = Maasha::SQL::connect( $options->{ "database" }, $options->{ "user" }, $options->{ "password" } );
3762
3763     foreach $track ( sort keys %track_hash )
3764     {
3765         if ( Maasha::SQL::table_exists( $dbh, $track ) )
3766         {
3767             print STDERR qq(Removing table "$track" from database "$options->{ 'database' }" ... ) if $options->{ 'verbose' };
3768             Maasha::SQL::delete_table( $dbh, $track );
3769             print STDERR "done.\n" if $options->{ 'verbose' };
3770         }
3771         else
3772         {
3773             print STDERR qq(WARNING: table "$track" not found in database "$options->{ 'database' }\n");
3774         }
3775     }
3776
3777     Maasha::SQL::disconnect( $dbh );
3778
3779     Maasha::Common::run( "ucscMakeTracks.pl", "-b > /dev/null 2>&1" );
3780 }
3781
3782
3783 sub script_rename_keys
3784 {
3785     # Martin A. Hansen, August 2007.
3786
3787     # Rename keys in stream.
3788
3789     my ( $in,        # handle to in stream
3790          $out,       # handle to out stream
3791          $options,   # options hash
3792        ) = @_;
3793
3794     # Returns nothing.
3795
3796     my ( $record );
3797
3798     while ( $record = Maasha::Biopieces::get_record( $in ) )
3799     {
3800         if ( exists $record->{ $options->{ "keys" }->[ 0 ] } )
3801         {
3802             $record->{ $options->{ "keys" }->[ 1 ] } = $record->{ $options->{ "keys" }->[ 0 ] };
3803
3804             delete $record->{ $options->{ "keys" }->[ 0 ] };
3805         }
3806
3807         Maasha::Biopieces::put_record( $record, $out );
3808     }
3809 }
3810
3811
3812 sub script_uniq_vals
3813 {
3814     # Martin A. Hansen, August 2007.
3815
3816     # Find unique values in stream.
3817
3818     my ( $in,        # handle to in stream
3819          $out,       # handle to out stream
3820          $options,   # options hash
3821        ) = @_;
3822
3823     # Returns nothing.
3824
3825     my ( %hash, $record );
3826
3827     while ( $record = Maasha::Biopieces::get_record( $in ) )
3828     {
3829         if ( $record->{ $options->{ "key" } } )
3830         {
3831             if ( not $hash{ $record->{ $options->{ "key" } } } and not $options->{ "invert" } )
3832             {
3833                 Maasha::Biopieces::put_record( $record, $out );
3834
3835                 $hash{ $record->{ $options->{ "key" } } } = 1;
3836             }
3837             elsif ( $hash{ $record->{ $options->{ "key" } } } and $options->{ "invert" } )
3838             {
3839                 Maasha::Biopieces::put_record( $record, $out );
3840             }
3841             else
3842             {
3843                 $hash{ $record->{ $options->{ "key" } } } = 1;
3844             }
3845         }
3846         else
3847         {
3848             Maasha::Biopieces::put_record( $record, $out );
3849         }
3850     }
3851 }
3852
3853
3854 sub script_merge_vals
3855 {
3856     # Martin A. Hansen, August 2007.
3857
3858     # Rename keys in stream.
3859
3860     my ( $in,        # handle to in stream
3861          $out,       # handle to out stream
3862          $options,   # options hash
3863        ) = @_;
3864
3865     # Returns nothing.
3866
3867     my ( $record, @join, $i );
3868
3869     $options->{ "delimit" } ||= '_';
3870
3871     while ( $record = Maasha::Biopieces::get_record( $in ) )
3872     {
3873         if ( exists $record->{ $options->{ "keys" }->[ 0 ] } )
3874         {
3875             @join = $record->{ $options->{ "keys" }->[ 0 ] };
3876
3877             for ( $i = 1; $i < @{ $options->{ "keys" } }; $i++ ) {
3878                 push @join, $record->{ $options->{ "keys" }->[ $i ] } if exists $record->{ $options->{ "keys" }->[ $i ] };
3879             }
3880
3881             $record->{ $options->{ "keys" }->[ 0 ] } = join $options->{ "delimit" }, @join;
3882         }
3883
3884         Maasha::Biopieces::put_record( $record, $out );
3885     }
3886 }
3887
3888
3889 sub script_merge_records
3890 {
3891     # Martin A. Hansen, July 2008.
3892
3893     # Merges records in the stream based on identical values of two given keys.
3894
3895     my ( $in,        # handle to in stream
3896          $out,       # handle to out stream
3897          $options,   # options hash
3898        ) = @_;
3899
3900     # Returns nothing.
3901
3902     my ( $merge, $record, $file1, $file2, $fh1, $fh2, $key1, $key2, @keys1, @keys2, @vals1, @vals2,
3903          $num1, $num2, $num, $cmp, $i );
3904
3905     $merge = $options->{ "merge" } || "AandB";
3906
3907     $file1 = "$BP_TMP/merge_records1.tmp";
3908     $file2 = "$BP_TMP/merge_records2.tmp";
3909
3910     $fh1   = Maasha::Common::write_open( $file1 );
3911     $fh2   = Maasha::Common::write_open( $file2 );
3912
3913     $key1  = $options->{ "keys" }->[ 0 ];
3914     $key2  = $options->{ "keys" }->[ 1 ];
3915
3916     $num   = $key2 =~ s/n$//;
3917     $num1  = 0;
3918     $num2  = 0;
3919
3920     while ( $record = Maasha::Biopieces::get_record( $in ) )
3921     {
3922         if ( exists $record->{ $key1 } )
3923         {
3924             @keys1 = $key1;
3925             @vals1 = $record->{ $key1 };
3926
3927             delete $record->{ $key1 };
3928
3929             map { push @keys1, $_; push @vals1, $record->{ $_ } } keys %{ $record };
3930
3931             print $fh1 join( "\t", @vals1 ), "\n";
3932
3933             $num1++;
3934         }
3935         elsif ( exists $record->{ $key2 } )
3936         {
3937             @keys2 = $key2;
3938             @vals2 = $record->{ $key2 };
3939
3940             delete $record->{ $key2 };
3941
3942             map { push @keys2, $_; push @vals2, $record->{ $_ } } keys %{ $record };
3943
3944             print $fh2 join( "\t", @vals2 ), "\n";
3945
3946             $num2++;
3947         }
3948     }
3949
3950     close $fh1;
3951     close $fh2;
3952
3953     if ( $num )
3954     {
3955         Maasha::Common::run( "sort", "-k 1,1n $file1 > $file1.sort" ) and rename "$file1.sort", $file1;
3956         Maasha::Common::run( "sort", "-k 1,1n $file2 > $file2.sort" ) and rename "$file2.sort", $file2;
3957     }
3958     else
3959     {
3960         Maasha::Common::run( "sort", "-k 1,1 $file1 > $file1.sort" ) and rename "$file1.sort", $file1;
3961         Maasha::Common::run( "sort", "-k 1,1 $file2 > $file2.sort" ) and rename "$file2.sort", $file2;
3962     }
3963
3964     $fh1 = Maasha::Common::read_open( $file1 );
3965     $fh2 = Maasha::Common::read_open( $file2 );
3966
3967     @vals1 = Maasha::Common::get_fields( $fh1 );
3968     @vals2 = Maasha::Common::get_fields( $fh2 );
3969
3970     while ( $num1 > 0 and $num2 > 0 )
3971     {
3972         undef $record;
3973
3974         if ( $num ) {
3975             $cmp = $vals1[ 0 ] <=> $vals2[ 0 ];
3976         } else {
3977             $cmp = $vals1[ 0 ] cmp $vals2[ 0 ];
3978         }
3979
3980         if ( $cmp < 0 )
3981         {
3982             if ( $merge =~ /^(AorB|AnotB)$/ )
3983             {
3984                 for ( $i = 0; $i < @keys1; $i++ ) {
3985                     $record->{ $keys1[ $i ] } = $vals1[ $i ];
3986                 }
3987
3988                 Maasha::Biopieces::put_record( $record, $out );
3989             }
3990
3991             @vals1 = Maasha::Common::get_fields( $fh1 );
3992             $num1--;
3993         }
3994         elsif ( $cmp > 0 )
3995         {
3996             if ( $merge =~ /^(BorA|BnotA)$/ )
3997             {
3998                 for ( $i = 0; $i < @keys2; $i++ ) {
3999                     $record->{ $keys2[ $i ] } = $vals2[ $i ];
4000                 }
4001
4002                 Maasha::Biopieces::put_record( $record, $out );
4003             }
4004
4005             @vals2 = Maasha::Common::get_fields( $fh2 );
4006             $num2--;
4007         }
4008         else
4009         {
4010             if ( $merge =~ /^(AandB|AorB|BorA)$/ )
4011             {
4012                 for ( $i = 0; $i < @keys1; $i++ ) {
4013                     $record->{ $keys1[ $i ] } = $vals1[ $i ];
4014                 }
4015
4016                 for ( $i = 1; $i < @keys2; $i++ ) {
4017                     $record->{ $keys2[ $i ] } = $vals2[ $i ];
4018                 }
4019
4020                 Maasha::Biopieces::put_record( $record, $out );
4021             }
4022
4023             @vals1 = Maasha::Common::get_fields( $fh1 );
4024             @vals2 = Maasha::Common::get_fields( $fh2 );
4025             $num1--;
4026             $num2--;
4027         }
4028     }
4029
4030     close $fh1;
4031     close $fh2;
4032
4033     unlink $file1;
4034     unlink $file2;
4035
4036     if ( $num1 > 0 and $merge =~ /^(AorB|AnotB)$/ )
4037     {
4038         undef $record;
4039
4040         for ( $i = 0; $i < @keys1; $i++ ) {
4041             $record->{ $keys1[ $i ] } = $vals1[ $i ];
4042         }
4043
4044         Maasha::Biopieces::put_record( $record, $out );
4045     }
4046
4047     if ( $num2 > 0 and $merge =~ /^(BorA|BnotA)$/ )
4048     {
4049         undef $record;
4050
4051         for ( $i = 0; $i < @keys2; $i++ ) {
4052             $record->{ $keys2[ $i ] } = $vals2[ $i ];
4053         }
4054
4055         Maasha::Biopieces::put_record( $record, $out );
4056     }
4057 }
4058
4059
4060 sub script_compute
4061 {
4062     # Martin A. Hansen, August 2007.
4063
4064     # Evaluate extression for records in stream.
4065
4066     my ( $in,        # handle to in stream
4067          $out,       # handle to out stream
4068          $options,   # options hash
4069        ) = @_;
4070
4071     # Returns nothing.
4072
4073     my ( $record, $eval_key, @keys, $eval_val );
4074
4075     while ( $record = Maasha::Biopieces::get_record( $in ) )
4076     {
4077         if ( $options->{ "eval" } )
4078         {
4079             if ( $options->{ "eval" } =~ /^(\S+)\s*=\s*(.+)$/ )
4080             {
4081                 $eval_key = $1;
4082                 $eval_val = $2;
4083
4084                 if ( not @keys )
4085                 {
4086                     @keys = split /\s+|\+|-|\*|\/|\*\*/, $eval_val;
4087
4088                     @keys = grep { exists $record->{ $_ } } @keys;
4089                 }
4090
4091                 map { $eval_val =~ s/\Q$_\E/$record->{ $_ }/g } @keys;
4092
4093                 $record->{ $eval_key } = eval "$eval_val";
4094                 Maasha::Common::error( qq(eval "$eval_key = $eval_val" failed -> $@) ) if $@;
4095             }
4096             else
4097             {
4098                 Maasha::Common::error( qq(Bad compute expression: "$options->{ 'eval' }"\n) );
4099             }
4100         }
4101
4102         Maasha::Biopieces::put_record( $record, $out );
4103     }
4104 }
4105
4106
4107 sub script_flip_tab
4108 {
4109     # Martin A. Hansen, June 2008.
4110
4111     # Flip a table.
4112
4113     my ( $in,        # handle to in stream
4114          $out,       # handle to out stream
4115          $options,   # options hash
4116        ) = @_;
4117
4118     # Returns nothing.
4119
4120     my ( $record, $key, $A, $B, @rows, @matrix, $row, $i );
4121
4122     while ( $record = Maasha::Biopieces::get_record( $in ) )
4123     {
4124         undef @rows;
4125
4126         foreach $key ( sort { $A = $a; $B = $b; $A =~ s/^V(\d+)$/$1/; $B =~ s/^V(\d+)$/$1/; $A <=> $B } keys %{ $record } )
4127         {
4128             push @rows, $record->{ $key };
4129
4130         }
4131
4132         push @matrix, [ @rows ];
4133     }
4134
4135     undef $record;
4136
4137     @matrix = Maasha::Matrix::matrix_flip( \@matrix );
4138
4139     foreach $row ( @matrix )
4140     {
4141         for ( $i = 0; $i < @{ $row }; $i++ ) {
4142             $record->{ "V$i" } = $row->[ $i ];
4143         }
4144
4145         Maasha::Biopieces::put_record( $record, $out );
4146     }
4147 }
4148
4149
4150 sub script_add_ident
4151 {
4152     # Martin A. Hansen, May 2008.
4153
4154     # Add a unique identifier to each record in stream.
4155
4156     my ( $in,        # handle to in stream
4157          $out,       # handle to out stream
4158          $options,   # options hash
4159        ) = @_;
4160
4161     # Returns nothing.
4162
4163     my ( $record, $key, $prefix, $i );
4164
4165     $key    = $options->{ "key" }    || "ID";
4166     $prefix = $options->{ "prefix" } || "ID";
4167
4168     $i = 0;
4169
4170     while ( $record = Maasha::Biopieces::get_record( $in ) )
4171     {
4172         $record->{ $key } = sprintf( "$prefix%08d", $i );
4173
4174         Maasha::Biopieces::put_record( $record, $out );
4175
4176         $i++;
4177     }
4178 }
4179
4180
4181 sub script_count_records
4182 {
4183     # Martin A. Hansen, August 2007.
4184
4185     # Count records in stream.
4186
4187     my ( $in,        # handle to in stream
4188          $out,       # handle to out stream
4189          $options,   # options hash
4190        ) = @_;
4191
4192     # Returns nothing.
4193
4194     my ( $record, $count, $result, $fh, $line );
4195
4196     $count = 0;
4197
4198     if ( $options->{ "no_stream" } )
4199     {
4200         while ( $line = <$in> )
4201         {
4202             chomp $line;
4203
4204             $count++ if $line eq "---";
4205         }
4206     }
4207     else
4208     {
4209         while ( $record = Maasha::Biopieces::get_record( $in ) )
4210         {
4211             Maasha::Biopieces::put_record( $record, $out );
4212
4213             $count++;
4214         }
4215     }
4216
4217     $result = { "RECORDS_COUNT" => $count };
4218
4219     $fh = write_stream( $options->{ "data_out" } );
4220
4221     Maasha::Biopieces::put_record( $result, $fh );
4222
4223     close $fh;
4224 }
4225
4226
4227 sub script_random_records
4228 {
4229     # Martin A. Hansen, August 2007.
4230
4231     # Pick a number or random records from stream.
4232
4233     my ( $in,        # handle to in stream
4234          $out,       # handle to out stream
4235          $options,   # options hash
4236        ) = @_;
4237
4238     # Returns nothing.
4239
4240     my ( $record, $tmp_file, $fh_out, $fh_in, $count, $i, %rand_hash, $rand, $max );
4241
4242     $options->{ "num" } ||= 10;
4243
4244     $tmp_file = "$BP_TMP/random_records.tmp";
4245
4246     $fh_out = Maasha::Common::write_open( $tmp_file );
4247
4248     $count = 0;
4249
4250     while ( $record = Maasha::Biopieces::get_record( $in ) )
4251     {
4252         Maasha::Biopieces::put_record( $record, $fh_out );
4253
4254         $count++;
4255     }
4256
4257     close $fh_out;
4258
4259     $max = 0;
4260     $i   = 0;
4261
4262     Maasha::Common::error( qq(Requested random records > records in stream) ) if $options->{ "num" } > $count;
4263
4264     while ( $i < $options->{ "num" } )
4265     {
4266         $rand = int( rand( $count ) );
4267
4268         if ( not exists $rand_hash{ $rand } )
4269         {
4270             $rand_hash{ $rand } = 1;
4271
4272             $max = $rand if $rand > $max;
4273
4274             $i++;
4275         }
4276     }
4277
4278     $fh_in = Maasha::Common::read_open( $tmp_file );
4279
4280     $count = 0;
4281
4282     while ( $record = Maasha::Biopieces::get_record( $fh_in ) )
4283     {
4284         Maasha::Biopieces::put_record( $record, $out ) if exists $rand_hash{ $count };
4285
4286         last if $count == $max;
4287
4288         $count++;
4289     }
4290
4291     close $fh_in;
4292
4293     unlink $tmp_file;
4294 }
4295
4296
4297 sub script_sort_records
4298 {
4299     # Martin A. Hansen, August 2007.
4300
4301     # Sort to sort records according to keys.
4302
4303     my ( $in,        # handle to in stream
4304          $out,       # handle to out stream
4305          $options,   # options hash
4306        ) = @_;
4307
4308     # Returns nothing.
4309
4310     my ( @keys, $key, @sort_cmd, $sort_str, $sort_sub, @records, $record, $i );
4311
4312     foreach $key ( @{ $options->{ "keys" } } )
4313     {
4314         if ( $key =~ s/n$// ) {
4315             push @sort_cmd, qq(\$a->{ "$key" } <=> \$b->{ "$key" });
4316         } else {
4317             push @sort_cmd, qq(\$a->{ "$key" } cmp \$b->{ "$key" });
4318         }
4319     }
4320
4321     $sort_str = join " or ", @sort_cmd;
4322     $sort_sub = eval "sub { $sort_str }";   # NB security issue!
4323
4324     while ( $record = Maasha::Biopieces::get_record( $in ) ) {
4325         push @records, $record;
4326     }
4327
4328     @records = sort $sort_sub @records;
4329
4330     if ( $options->{ "reverse" } )
4331     {
4332         for ( $i = scalar @records - 1; $i >= 0; $i-- ) {
4333             Maasha::Biopieces::put_record( $records[ $i ], $out );
4334         }
4335     }
4336     else
4337     {
4338         for ( $i = 0; $i < scalar @records; $i++ ) {
4339             Maasha::Biopieces::put_record( $records[ $i ], $out );
4340         }
4341     }
4342 }
4343
4344
4345 sub script_count_vals
4346 {
4347     # Martin A. Hansen, August 2007.
4348
4349     # Count records in stream.
4350
4351     my ( $in,        # handle to in stream
4352          $out,       # handle to out stream
4353          $options,   # options hash
4354        ) = @_;
4355
4356     # Returns nothing.
4357
4358     my ( $num, $record, %count_hash, @records, $tmp_file, $fh_out, $fh_in, $cache );
4359
4360     $tmp_file = "$BP_TMP/count_cache.tmp";
4361
4362     $fh_out   = Maasha::Common::write_open( $tmp_file );
4363
4364     $cache    = 0;
4365     $num      = 0;
4366
4367     while ( $record = Maasha::Biopieces::get_record( $in ) )
4368     {
4369         map { $count_hash{ $_ }{ $record->{ $_ } }++ if exists $record->{ $_ } } @{ $options->{ "keys" } };
4370
4371         push @records, $record;
4372
4373         if ( scalar @records > 5_000_000 )   # too many records to hold in memory - use disk cache
4374         {
4375             map { Maasha::Biopieces::put_record( $_, $fh_out ) } @records;
4376
4377             undef @records;
4378
4379             $cache = 1;
4380         }
4381
4382         print STDERR "verbose: records read $num\n" if ( $options->{ 'verbose' } and ( $num % 1_000_000 ) == 0 );
4383
4384         $num++;
4385     }
4386
4387     close $fh_out;
4388
4389     if ( $cache )
4390     {
4391         $num   = 0;
4392
4393         $fh_in = Maasha::Common::read_open( $tmp_file );
4394
4395         while ( $record = Maasha::Biopieces::get_record( $fh_in ) )
4396         {
4397             map { $record->{ $_ . "_COUNT" } = $count_hash{ $_ }{ $record->{ $_ } } if exists $record->{ $_ } } @{ $options->{ "keys" } };
4398
4399             Maasha::Biopieces::put_record( $record, $out );
4400
4401             print STDERR "verbose: cache read $num\n" if ( $options->{ 'verbose' } and ( $num % 1_000_000 ) == 0 );
4402
4403             $num++;
4404         }
4405
4406         close $fh_in;
4407     }
4408
4409     foreach $record ( @records )
4410     {
4411         map { $record->{ $_ . "_COUNT" } = $count_hash{ $_ }{ $record->{ $_ } } if exists $record->{ $_ } } @{ $options->{ "keys" } };
4412
4413         Maasha::Biopieces::put_record( $record, $out );
4414     }
4415
4416     unlink $tmp_file;
4417 }
4418
4419
4420 sub script_plot_histogram
4421 {
4422     # Martin A. Hansen, September 2007.
4423
4424     # Plot a simple histogram for a given key using GNU plot.
4425
4426     my ( $in,        # handle to in stream
4427          $out,       # handle to out stream
4428          $options,   # options hash
4429        ) = @_;
4430
4431     # Returns nothing.
4432
4433     my ( $record, %data_hash, $max, @data_list, $i, $result, $fh );
4434
4435     $options->{ "title" } ||= "Histogram";
4436     $options->{ "sort" }  ||= "num";
4437
4438     while ( $record = Maasha::Biopieces::get_record( $in ) )
4439     {
4440         $data_hash{ $record->{ $options->{ "key" } } }++ if defined $record->{ $options->{ "key" } };
4441
4442         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
4443     }
4444
4445     if ( $options->{ "sort" } eq "num" ) {
4446         map { push @data_list, [ $_, $data_hash{ $_ } ] } sort { $a <=> $b } keys %data_hash;
4447     } else {
4448         map { push @data_list, [ $_, $data_hash{ $_ } ] } sort keys %data_hash;
4449     }
4450
4451     $result = Maasha::Plot::histogram_simple( \@data_list, $options );
4452
4453     $fh = write_stream( $options->{ "data_out" } );
4454
4455     print $fh "$_\n" foreach @{ $result };
4456
4457     close $fh;
4458 }
4459
4460
4461 sub script_plot_lendist
4462 {
4463     # Martin A. Hansen, August 2007.
4464
4465     # Plot length distribution using GNU plot.
4466
4467     my ( $in,        # handle to in stream
4468          $out,       # handle to out stream
4469          $options,   # options hash
4470        ) = @_;
4471
4472     # Returns nothing.
4473
4474     my ( $record, %data_hash, $max, @data_list, $i, $result, $fh );
4475
4476     $options->{ "title" } ||= "Length Distribution";
4477
4478     while ( $record = Maasha::Biopieces::get_record( $in ) )
4479     {
4480         $data_hash{ $record->{ $options->{ "key" } } }++ if defined $record->{ $options->{ "key" } };
4481
4482         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
4483     }
4484
4485     $max = Maasha::Calc::list_max( [ keys %data_hash ] );
4486
4487     for ( $i = 0; $i < $max; $i++ ) {
4488         push @data_list, [ $i, $data_hash{ $i } || 0 ];
4489     }
4490
4491     $result = Maasha::Plot::histogram_lendist( \@data_list, $options );
4492
4493     $fh = write_stream( $options->{ "data_out" } );
4494
4495     print $fh "$_\n" foreach @{ $result };
4496
4497     close $fh;
4498 }
4499
4500
4501 sub script_plot_chrdist
4502 {
4503     # Martin A. Hansen, August 2007.
4504
4505     # Plot chromosome distribution using GNU plot.
4506
4507     my ( $in,        # handle to in stream
4508          $out,       # handle to out stream
4509          $options,   # options hash
4510        ) = @_;
4511
4512     # Returns nothing.
4513
4514     my ( $record, %data_hash, @data_list, $elem, $sort_key, $count, $result, $fh );
4515
4516     $options->{ "title" } ||= "Chromosome Distribution";
4517
4518     while ( $record = Maasha::Biopieces::get_record( $in ) )
4519     {
4520         if ( $record->{ "CHR" } ) {                                                             # generic
4521             $data_hash{ $record->{ "CHR" } }++;
4522         } elsif ( $record->{ "REC_TYPE" } eq "PATSCAN" and $record->{ "S_ID" } =~ /^chr/i ) {   # patscan
4523             $data_hash{ $record->{ "S_ID" } }++;
4524         } elsif ( $record->{ "REC_TYPE" } eq "PSL" and $record->{ "S_ID" } =~ /^chr/i ) {       # BLAT / PSL
4525             $data_hash{ $record->{ "S_ID" } }++;
4526         } elsif ( $record->{ "REC_TYPE" } eq "BLAST" and $record->{ "S_ID" } =~ /^chr/i ) {     # BLAST
4527             $data_hash{ $record->{ "S_ID" } }++;
4528         }
4529
4530         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
4531     }
4532
4533     foreach $elem ( keys %data_hash )
4534     {
4535         $sort_key = $elem;
4536
4537         $sort_key =~ s/chr//i;
4538
4539         $sort_key =~ s/^X(.*)/99$1/;
4540         $sort_key =~ s/^Y(.*)/99$1/;
4541         $sort_key =~ s/^Z(.*)/999$1/;
4542         $sort_key =~ s/^M(.*)/9999$1/;
4543         $sort_key =~ s/^U(.*)/99999$1/;
4544
4545         $count = $sort_key =~ tr/_//;
4546
4547         $sort_key =~ s/_.*/"999999" x $count/ex;
4548
4549         push @data_list, [ $elem, $data_hash{ $elem }, $sort_key ];
4550     }
4551
4552     @data_list = sort { $a->[ 2 ] <=> $b->[ 2 ] } @data_list;
4553
4554     $result = Maasha::Plot::histogram_chrdist( \@data_list, $options );
4555
4556     $fh = write_stream( $options->{ "data_out" } );
4557
4558     print $fh "$_\n" foreach @{ $result };
4559
4560     close $fh;
4561 }
4562
4563
4564 sub script_plot_karyogram
4565 {
4566     # Martin A. Hansen, August 2007.
4567
4568     # Plot hits on karyogram.
4569
4570     my ( $in,        # handle to in stream
4571          $out,       # handle to out stream
4572          $options,   # options hash
4573        ) = @_;
4574
4575     # Returns nothing.
4576
4577     my ( %options, $record, @data, $fh, $result, %data_hash );
4578
4579     $options->{ "genome" }     ||= "human";
4580     $options->{ "feat_color" } ||= "black";
4581
4582     while ( $record = Maasha::Biopieces::get_record( $in ) )
4583     {
4584         if ( $record->{ "CHR" } and $record->{ "CHR_BEG" } and $record->{ "CHR_END" } )
4585         {
4586             push @{ $data_hash{ $record->{ "CHR" } } }, [ $record->{ "CHR_BEG" }, $record->{ "CHR_END" }, $options->{ "feat_color" } ];
4587         }
4588
4589         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
4590     }
4591
4592     $result = Maasha::Plot::karyogram( \%data_hash, \%options );
4593
4594     $fh = write_stream( $options->{ "data_out" } );
4595
4596     print $fh $result;
4597
4598     close $fh;
4599 }
4600
4601
4602 sub script_plot_matches
4603 {
4604     # Martin A. Hansen, August 2007.
4605
4606     # Plot matches in 2D generating a dotplot.
4607
4608     my ( $in,        # handle to in stream
4609          $out,       # handle to out stream
4610          $options,   # options hash
4611        ) = @_;
4612
4613     # Returns nothing.
4614
4615     my ( $record, @data, $fh, $result, %data_hash );
4616
4617     $options->{ "direction" } ||= "both";
4618
4619     while ( $record = Maasha::Biopieces::get_record( $in ) )
4620     {
4621         if ( defined $record->{ "Q_BEG" } and defined $record->{ "S_BEG" } and $record->{ "Q_END" } and $record->{ "S_END" } ) {
4622             push @data, $record;
4623         }
4624
4625         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
4626     }
4627
4628     $options->{ "title" }  ||= "plot_matches";
4629     $options->{ "xlabel" } ||= $data[ 0 ]->{ "Q_ID" };
4630     $options->{ "ylabel" } ||= $data[ 0 ]->{ "S_ID" };
4631
4632     $result = Maasha::Plot::dotplot_matches( \@data, $options, $BP_TMP );
4633
4634     $fh = write_stream( $options->{ "data_out" } );
4635
4636     print $fh "$_\n" foreach @{ $result };
4637
4638     close $fh;
4639 }
4640
4641
4642 sub script_length_vals
4643 {
4644     # Martin A. Hansen, August 2007.
4645
4646     # Determine the length of the value for given keys.
4647
4648     my ( $in,        # handle to in stream
4649          $out,       # handle to out stream
4650          $options,   # options hash
4651        ) = @_;
4652
4653     # Returns nothing.
4654
4655     my ( $record, $key );
4656
4657     while ( $record = Maasha::Biopieces::get_record( $in ) )
4658     {
4659         foreach $key ( @{ $options->{ "keys" } } )
4660         {
4661             if ( $record->{ $key } ) {
4662                 $record->{ $key . "_LEN" } = length $record->{ $key };
4663             }
4664         }
4665
4666         Maasha::Biopieces::put_record( $record, $out );
4667     }
4668 }
4669
4670
4671 sub script_sum_vals
4672 {
4673     # Martin A. Hansen, August 2007.
4674
4675     # Calculates the sums for values of given keys.
4676
4677     my ( $in,        # handle to in stream
4678          $out,       # handle to out stream
4679          $options,   # options hash
4680        ) = @_;
4681
4682     # Returns nothing.
4683
4684     my ( $record, $key, %sum_hash, $fh );
4685
4686     while ( $record = Maasha::Biopieces::get_record( $in ) )
4687     {
4688         foreach $key ( @{ $options->{ "keys" } } )
4689         {
4690             if ( $record->{ $key } ) {
4691                 $sum_hash{ $key } += $record->{ $key };
4692             }
4693         }
4694
4695         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
4696     }
4697
4698     $fh = write_stream( $options->{ "data_out" } );
4699
4700     foreach $key ( @{ $options->{ "keys" } } ) {
4701         Maasha::Biopieces::put_record( { $key . "_SUM" => $sum_hash{ $key } || 0 } , $fh );
4702     }
4703
4704     close $fh;
4705 }
4706
4707
4708 sub script_mean_vals
4709 {
4710     # Martin A. Hansen, August 2007.
4711
4712     # Calculate the mean of values of given keys.
4713
4714     my ( $in,        # handle to in stream
4715          $out,       # handle to out stream
4716          $options,   # options hash
4717        ) = @_;
4718
4719     # Returns nothing.
4720
4721     my ( $record, $key, %sum_hash, %count_hash, $mean, $fh );
4722
4723     while ( $record = Maasha::Biopieces::get_record( $in ) )
4724     {
4725         foreach $key ( @{ $options->{ "keys" } } )
4726         {
4727             if ( $record->{ $key } )
4728             {
4729                 $sum_hash{ $key } += $record->{ $key };
4730                 $count_hash{ $key }++;
4731             }
4732         }
4733
4734         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
4735     }
4736
4737     $fh = write_stream( $options->{ "data_out" } );
4738
4739     foreach $key ( @{ $options->{ "keys" } } )
4740     {
4741         if ( $count_hash{ $key } ) {
4742             $mean = sprintf( "%.2f", ( $sum_hash{ $key } / $count_hash{ $key } ) );
4743         } else {
4744             $mean = "N/A";
4745         }
4746
4747         Maasha::Biopieces::put_record( { $key . "_MEAN" => $mean } , $fh );
4748     }
4749
4750     close $fh;
4751 }
4752
4753
4754 sub script_median_vals
4755 {
4756     # Martin A. Hansen, March 2008.
4757
4758     # Calculate the median values of given keys.
4759
4760     my ( $in,        # handle to in stream
4761          $out,       # handle to out stream
4762          $options,   # options hash
4763        ) = @_;
4764
4765     # Returns nothing.
4766
4767     my ( $record, $key, %median_hash, $median, $fh );
4768
4769     while ( $record = Maasha::Biopieces::get_record( $in ) )
4770     {
4771         foreach $key ( @{ $options->{ "keys" } } ) {
4772             push @{ $median_hash{ $key } }, $record->{ $key } if defined $record->{ $key };
4773         }
4774
4775         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
4776     }
4777
4778     $fh = write_stream( $options->{ "data_out" } );
4779
4780     foreach $key ( @{ $options->{ "keys" } } )
4781     {
4782         if ( $median_hash{ $key } ) {
4783             $median = Maasha::Calc::median( $median_hash{ $key } );
4784         } else {
4785             $median = "N/A";
4786         }
4787
4788         Maasha::Biopieces::put_record( { $key . "_MEDIAN" => $median } , $fh );
4789     }
4790
4791     close $fh;
4792 }
4793
4794
4795 sub script_max_vals
4796 {
4797     # Martin A. Hansen, February 2008.
4798
4799     # Determine the maximum values of given keys.
4800
4801     my ( $in,        # handle to in stream
4802          $out,       # handle to out stream
4803          $options,   # options hash
4804        ) = @_;
4805
4806     # Returns nothing.
4807
4808     my ( $record, $key, $fh, %max_hash, $max_record );
4809
4810     while ( $record = Maasha::Biopieces::get_record( $in ) )
4811     {
4812         foreach $key ( @{ $options->{ "keys" } } )
4813         {
4814             if ( $record->{ $key } )
4815             {
4816                 $max_hash{ $key } = $record->{ $key } if $record->{ $key } > $max_hash{ $key };
4817             }
4818         }
4819
4820         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
4821     }
4822
4823     $fh = write_stream( $options->{ "data_out" } );
4824
4825     foreach $key ( @{ $options->{ "keys" } } )
4826     {
4827         $max_record->{ $key . "_MAX" } = $max_hash{ $key };
4828     }
4829
4830     Maasha::Biopieces::put_record( $max_record, $fh );
4831
4832     close $fh;
4833 }
4834
4835
4836 sub script_min_vals
4837 {
4838     # Martin A. Hansen, February 2008.
4839
4840     # Determine the minimum values of given keys.
4841
4842     my ( $in,        # handle to in stream
4843          $out,       # handle to out stream
4844          $options,   # options hash
4845        ) = @_;
4846
4847     # Returns nothing.
4848
4849     my ( $record, $key, $fh, %min_hash, $min_record );
4850
4851     while ( $record = Maasha::Biopieces::get_record( $in ) )
4852     {
4853         foreach $key ( @{ $options->{ "keys" } } )
4854         {
4855             if ( defined $record->{ $key } )
4856             {
4857                 if ( exists $min_hash{ $key } ) {
4858                     $min_hash{ $key } = $record->{ $key } if $record->{ $key } < $min_hash{ $key };
4859                 } else {
4860                     $min_hash{ $key } = $record->{ $key };
4861                 }
4862             }
4863         }
4864
4865         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
4866     }
4867
4868     $fh = write_stream( $options->{ "data_out" } );
4869
4870     foreach $key ( @{ $options->{ "keys" } } )
4871     {
4872         $min_record->{ $key . "_MIN" } = $min_hash{ $key };
4873     }
4874
4875     Maasha::Biopieces::put_record( $min_record, $fh );
4876
4877     close $fh;
4878 }
4879
4880
4881 sub script_upload_to_ucsc
4882 {
4883     # Martin A. Hansen, August 2007.
4884
4885     # Calculate the mean of values of given keys.
4886
4887     # This routine has developed into the most ugly hack. Do something!
4888
4889     my ( $in,        # handle to in stream
4890          $out,       # handle to out stream
4891          $options,   # options hash
4892        ) = @_;
4893
4894     # Returns nothing.
4895
4896     my ( $record, $file, $wib_file, $wig_file, $wib_dir, $fh_out, $i, $first, $format, $type, $columns, $append, $entry );
4897
4898     $options->{ "short_label" } ||= $options->{ 'table' };
4899     $options->{ "long_label" }  ||= $options->{ 'table' };
4900     $options->{ "group" }       ||= $ENV{ "LOGNAME" };
4901     $options->{ "priority" }    ||= 1;
4902     $options->{ "visibility" }  ||= "pack";
4903     $options->{ "color" }       ||= join( ",", int( rand( 255 ) ), int( rand( 255 ) ), int( rand( 255 ) ) );
4904     $options->{ "chunk_size" }  ||= 10_000_000_000;    # Due to 32-bit UCSC compilation really large tables cannot be loaded in one go.
4905
4906     $file   = "$BP_TMP/ucsc_upload.tmp";
4907     $append = 0;
4908     $first  = 1;
4909     $i      = 0;
4910
4911     $fh_out = Maasha::Common::write_open( $file );
4912
4913     while ( $record = Maasha::Biopieces::get_record( $in ) )
4914     {
4915         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
4916
4917         if ( $record->{ "REC_TYPE" } eq "fixed_step" )
4918         {
4919             $format = "WIGGLE";
4920
4921             if ( $entry = Maasha::UCSC::Wiggle::biopiece2fixedstep( $record ) ) {
4922                 Maasha::UCSC::Wiggle::fixedstep_entry_put( $entry, $fh_out );
4923             }
4924         }
4925         elsif ( $record->{ "REC_TYPE" } eq "PSL" )
4926         {
4927             $format = "PSL";
4928
4929             Maasha::UCSC::psl_put_header( $fh_out ) if $first;
4930             Maasha::UCSC::psl_put_entry( $record, $fh_out );
4931
4932             $first = 0;
4933         }
4934         elsif ( $record->{ "REC_TYPE" } eq "BED" and $record->{ "SEC_STRUCT" } )
4935         {
4936             # chrom chromStart  chromEnd    name    score   strand  size    secStr  conf
4937
4938             $format  = "BED_SS";
4939
4940             print $fh_out join ( "\t",
4941                 $record->{ "CHR" },
4942                 $record->{ "CHR_BEG" },
4943                 $record->{ "CHR_END" } + 1,
4944                 $record->{ "Q_ID" },
4945                 $record->{ "SCORE" },
4946                 $record->{ "STRAND" },
4947                 $record->{ "SIZE" },
4948                 $record->{ "SEC_STRUCT" },
4949                 $record->{ "CONF" },
4950             ), "\n";
4951         }
4952         elsif ( $record->{ "REC_TYPE" } eq "BED" )
4953         {
4954             $format  = "BED";
4955             $columns = $record->{ "BED_COLS" };
4956
4957             if ( $entry = Maasha::UCSC::BED::biopiece2bed( $record, $columns ) ) {
4958                 Maasha::UCSC::BED::bed_entry_put( $entry, $fh_out, $columns, $options->{ 'check' } );
4959             }
4960         }
4961         elsif ( $record->{ "REC_TYPE" } eq "PATSCAN" and $record->{ "CHR" } )
4962         {
4963             $format  = "BED";
4964             $columns = 6;
4965
4966             if ( $entry = Maasha::UCSC::BED::biopiece2bed( $record, $columns ) ) {
4967                 Maasha::UCSC::BED::bed_entry_put( $entry, $fh_out, $columns, $options->{ 'check' } );
4968             }
4969         }
4970         elsif ( $record->{ "REC_TYPE" } eq "BLAST" and $record->{ "S_ID" } =~ /^chr/ )
4971         {
4972             $format  = "BED";
4973             $columns = 6;
4974
4975             $record->{ "SCORE" } = $record->{ "BIT_SCORE" } * 1000;
4976
4977             if ( $entry = Maasha::UCSC::BED::biopiece2bed( $record, $columns ) ) {
4978                 Maasha::UCSC::BED::bed_entry_put( $entry, $fh_out, $columns, $options->{ 'check' } );
4979             }
4980         }
4981         elsif ( $record->{ "REC_TYPE" } eq "VMATCH" and $record->{ "S_ID" } =~ /^chr/i )
4982         {
4983             $format  = "BED";
4984             $columns = 6;
4985
4986             if ( $entry = Maasha::UCSC::BED::biopiece2bed( $record, $columns ) ) {
4987                 Maasha::UCSC::BED::bed_entry_put( $entry, $fh_out, $columns, $options->{ 'check' } );
4988             }
4989         }
4990
4991         if ( $i == $options->{ "chunk_size" } )
4992         {
4993             close $fh_out;
4994
4995             if ( $format eq "BED" ) {
4996                 Maasha::UCSC::bed_upload_to_ucsc( $BP_TMP, $file, $options, $append );
4997             } elsif ( $format eq "PSL" ) {
4998                 Maasha::UCSC::psl_upload_to_ucsc( $file, $options, $append );
4999             }
5000
5001             unlink $file;
5002
5003             $first = 1;
5004
5005             $append = 1;
5006
5007             $fh_out = Maasha::Common::write_open( $file );
5008         }
5009
5010         $i++;
5011     }
5012
5013     close $fh_out;
5014
5015     if ( exists $options->{ "database" } and $options->{ "table" } )
5016     {
5017         if ( $format eq "BED" )
5018         {
5019             $type = "bed $columns";
5020
5021             Maasha::UCSC::bed_upload_to_ucsc( $BP_TMP, $file, $options, $append );
5022         }
5023         elsif ( $format eq "BED_SS" )
5024         {
5025             $type = "type bed 6 +";
5026
5027             Maasha::UCSC::bed_upload_to_ucsc( $BP_TMP, $file, $options, $append );
5028         }
5029         elsif ( $format eq "PSL" )
5030         {
5031             $type = "psl";
5032
5033             Maasha::UCSC::psl_upload_to_ucsc( $file, $options, $append );
5034         }
5035         elsif ( $format eq "WIGGLE" )
5036         {
5037             $options->{ "visibility" } = "full";
5038
5039             $wig_file = "$options->{ 'table' }.wig";
5040             $wib_file = "$options->{ 'table' }.wib";
5041
5042             $wib_dir  = "$ENV{ 'HOME' }/ucsc/wib";
5043
5044             Maasha::Common::dir_create_if_not_exists( $wib_dir );
5045
5046             if ( $options->{ 'verbose' } ) {
5047                 `cd $BP_TMP && wigEncode $file $wig_file $wib_file`;
5048             } else {
5049                 `cd $BP_TMP && wigEncode $file $wig_file $wib_file > /dev/null 2>&1`;
5050             }
5051
5052             Maasha::Common::run( "mv", "$BP_TMP/$wib_file $wib_dir" );
5053
5054             unlink $file;
5055
5056             $file = $wig_file;
5057
5058             $type = "wig 0";
5059
5060             Maasha::UCSC::wiggle_upload_to_ucsc( $BP_TMP, $wib_dir, $file, $options );
5061         }
5062
5063         unlink $file;
5064
5065         Maasha::UCSC::ucsc_update_config( $options, $type );
5066     }
5067 }
5068
5069
5070 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
5071
5072 1;
5073
5074 __END__
5075