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