code_perl/Maasha/Biopieces.pm

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