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