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|v
 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                 REC_TYPE   => 'SOLID',
2023                 SEQ_NAME   => $seq_name,
2024                 SEQ_CS     => $seq_cs,
2025                 SEQ_QUAL   => join( ";", @scores ),
2026                 SEQ_LEN    => length $seq_cs,
2027                 SEQ        => join( "", @seqs ),
2028                 SCORE_MEAN => sprintf( "%.2f", Maasha::Calc::mean( \@scores ) ),
2029             };
2030
2031             put_record( $record, $out );
2032
2033             goto NUM if $options->{ "num" } and $num == $options->{ "num" };
2034
2035             $num++;
2036         }
2037
2038         close $data_in;
2039     }
2040
2041     NUM:
2042
2043     close $data_in if $data_in;
2044 }
2045
2046
2047 sub script_read_mysql
2048 {
2049     # Martin A. Hansen, May 2008.
2050
2051     # Read a MySQL query into stream.
2052
2053     my ( $in,        # handle to in stream
2054          $out,       # handle to out stream
2055          $options,   # options hash
2056        ) = @_;
2057
2058     # Returns nothing.
2059
2060     my ( $record, $dbh, $results );
2061
2062     $options->{ "user" }     ||= Maasha::UCSC::ucsc_get_user();
2063     $options->{ "password" } ||= Maasha::UCSC::ucsc_get_password();
2064
2065     while ( $record = get_record( $in ) ) {
2066         put_record( $record, $out );
2067     }
2068
2069     $dbh = Maasha::SQL::connect( $options->{ "database" }, $options->{ "user" }, $options->{ "password" } );
2070
2071     $results = Maasha::SQL::query_hashref_list( $dbh, $options->{ "query" } );
2072
2073     Maasha::SQL::disconnect( $dbh );
2074
2075     map { put_record( $_ ) } @{ $results };
2076 }
2077
2078
2079 sub script_format_genome
2080 {
2081     # Martin A. Hansen, Juli 2008.
2082
2083     # Format a genome to speficed formats.
2084
2085     my ( $in,        # handle to in stream
2086          $out,       # handle to out stream
2087          $options,   # options hash
2088        ) = @_;
2089
2090     # Returns nothing.
2091
2092     my ( $dir, $genome, $fasta_dir, $phastcons_dir, $vals, $fh_out, $record, $format, $index, $entry );
2093
2094     $dir    = $options->{ 'dir' } || $ENV{ 'BP_DATA' };
2095     $genome = $options->{ 'genome' };
2096
2097     Maasha::Common::error( "Directory: $dir does not exist" ) if not -d $dir;
2098     Maasha::Common::dir_create_if_not_exists( "$dir/genomes" );
2099     Maasha::Common::dir_create_if_not_exists( "$dir/genomes/$genome" );
2100
2101     if ( grep { $_ =~ /fasta|blast|vmatch/i } @{ $options->{ "formats" } } )
2102     {
2103         if ( -f "$dir/genomes/$genome/fasta/$genome.fna" )
2104         {
2105             $fasta_dir = "$dir/genomes/$genome/fasta";
2106         }
2107         else
2108         {
2109             Maasha::Common::dir_create_if_not_exists( "$dir/genomes/$genome/fasta" );
2110
2111             $fasta_dir = "$dir/genomes/$genome/fasta";
2112
2113             $fh_out = Maasha::Common::write_open( "$fasta_dir/$genome.fna" );
2114         }
2115     }
2116     elsif ( grep { $_ =~ /phastcons/i } @{ $options->{ "formats" } } )
2117     {
2118         Maasha::Common::dir_create_if_not_exists( "$dir/genomes/$genome/phastcons" );
2119
2120         $phastcons_dir = "$dir/genomes/$genome/phastcons";
2121
2122         $fh_out = Maasha::Common::write_open( "$phastcons_dir/$genome.pp" );
2123     }
2124
2125     while ( $record = get_record( $in ) )
2126     {
2127         if ( $fh_out and $entry = record2fasta( $record ) )
2128         {
2129             Maasha::Fasta::put_entry( $entry, $fh_out, $options->{ "wrap" } );
2130         }
2131         elsif ( $fh_out and $record->{ "CHR" } and $record->{ "CHR_BEG" } and $record->{ "STEP" } and $record->{ "VALS" } )
2132         {
2133             print $fh_out "fixedStep chrom=$record->{ 'CHR' } start=$record->{ 'CHR_BEG' } step=$record->{ 'STEP' }\n";
2134
2135             $vals = $record->{ 'VALS' };
2136
2137             $vals =~ tr/,/\n/;
2138
2139             print $fh_out "$vals\n";
2140         }
2141
2142         put_record( $record, $out ) if not $options->{ "no_stream" };
2143     }
2144
2145     foreach $format ( @{ $options->{ 'formats' } } )
2146     {
2147         if    ( $format =~ /^fasta$/i )     { Maasha::Fasta::fasta_index( "$fasta_dir/$genome.fna", "$dir/genomes/$genome/fasta/$genome.index" ) }
2148         elsif ( $format =~ /^blast$/i )     { Maasha::NCBI::blast_index( "$genome.fna", $fasta_dir, "$dir/genomes/$genome/blast", "dna", $genome ) }
2149         elsif ( $format =~ /^blat$/i )      { print STDERR "BLAT FORMAT NOT IMPLEMENTED" }
2150         elsif ( $format =~ /^vmatch$/i )    { Maasha::Match::vmatch_index( "$genome.fna", $fasta_dir, "$dir/genomes/$genome/vmatch", $BP_TMP ) }
2151         elsif ( $format =~ /^phastcons$/i ) { Maasha::UCSC::phastcons_index( "$genome.pp", $phastcons_dir ) }
2152     }
2153
2154     close $fh_out if $fh_out;
2155 }
2156
2157
2158 sub script_length_seq
2159 {
2160     # Martin A. Hansen, August 2007.
2161
2162     # Determine the length of sequences in stream.
2163
2164     my ( $in,        # handle to in stream
2165          $out,       # handle to out stream
2166          $options,   # options hash
2167        ) = @_;
2168
2169     # Returns nothing.
2170
2171     my ( $record, $total );
2172
2173     while ( $record = get_record( $in ) )
2174     {
2175         if ( $record->{ "SEQ" } )
2176         {
2177             $record->{ "SEQ_LEN" } = length $record->{ "SEQ" };
2178             $total += $record->{ "SEQ_LEN" };
2179         }
2180
2181         put_record( $record, $out ) if not $options->{ "no_stream" };
2182     }
2183
2184     put_record( { TOTAL_SEQ_LEN => $total }, $out );
2185 }
2186
2187
2188 sub script_uppercase_seq
2189 {
2190     # Martin A. Hansen, August 2007.
2191
2192     # Uppercases sequences in stream.
2193
2194     my ( $in,    # handle to in stream
2195          $out,   # handle to out stream
2196        ) = @_;
2197
2198     # Returns nothing.
2199
2200     my ( $record );
2201
2202     while ( $record = get_record( $in ) )
2203     {
2204         $record->{ "SEQ" } = uc $record->{ "SEQ" } if $record->{ "SEQ" };
2205
2206         put_record( $record, $out );
2207     }
2208 }
2209
2210
2211 sub script_shuffle_seq
2212 {
2213     # Martin A. Hansen, December 2007.
2214
2215     # Shuffle sequences in stream.
2216
2217     my ( $in,    # handle to in stream
2218          $out,   # handle to out stream
2219        ) = @_;
2220
2221     # Returns nothing.
2222
2223     my ( $record );
2224
2225     while ( $record = get_record( $in ) )
2226     {
2227         $record->{ "SEQ" } = Maasha::Seq::seq_shuffle( $record->{ "SEQ" } ) if $record->{ "SEQ" };
2228
2229         put_record( $record, $out );
2230     }
2231 }
2232
2233
2234 sub script_analyze_seq
2235 {
2236     # Martin A. Hansen, August 2007.
2237
2238     # Analyze sequence composition of sequences in stream.
2239
2240     my ( $in,     # handle to in stream
2241          $out,    # handle to out stream
2242        ) = @_;
2243
2244     # Returns nothing.
2245
2246     my ( $record, $analysis );
2247
2248     while ( $record = get_record( $in ) )
2249     {
2250         if ( $record->{ "SEQ" } )
2251         {
2252             $analysis = Maasha::Seq::seq_analyze( $record->{ "SEQ" } );
2253
2254             map { $record->{ $_ } = $analysis->{ $_ } } keys %{ $analysis };
2255         }
2256
2257         put_record( $record, $out );
2258     }
2259 }
2260
2261
2262 sub script_analyze_tags
2263 {
2264     # Martin A. Hansen, August 2008.
2265
2266     # Analyze sequence tags in stream.
2267
2268     my ( $in,     # handle to in stream
2269          $out,    # handle to out stream
2270        ) = @_;
2271
2272     # Returns nothing.
2273
2274     my ( $record, $analysis, %len_hash, %clone_hash, $clones, $key, $tag_record );
2275
2276     while ( $record = get_record( $in ) )
2277     {
2278         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } )
2279         {
2280             if ( $record->{ "SEQ_NAME" } =~ /_(\d+)$/ )
2281             {
2282                 $clones = $1;
2283
2284                 $len_hash{ length( $record->{ "SEQ" } ) }++;
2285                 $clone_hash{ length( $record->{ "SEQ" } ) } += $clones;
2286             }
2287         }
2288         elsif ( $record->{ "Q_ID" } and $record->{ "BED_LEN" } )
2289         {
2290             if ( $record->{ "Q_ID" } =~ /_(\d+)$/ )
2291             {
2292                 $clones = $1;
2293
2294                 $len_hash{ $record->{ "BED_LEN" } }++;
2295                 $clone_hash{ $record->{ "BED_LEN" } } += $clones;
2296             }
2297         }
2298     }
2299
2300     foreach $key ( sort { $a <=> $b } keys %len_hash )
2301     {
2302         $tag_record->{ "TAG_LEN" }    = $key;
2303         $tag_record->{ "TAG_COUNT" }  = $len_hash{ $key };
2304         $tag_record->{ "TAG_CLONES" } = $clone_hash{ $key };
2305
2306         put_record( $tag_record, $out );
2307     }
2308 }
2309
2310
2311 sub script_complexity_seq
2312 {
2313     # Martin A. Hansen, May 2008.
2314
2315     # Generates an index calculated as the most common di-residue over
2316     # the sequence length for all sequences in stream.
2317
2318     my ( $in,     # handle to in stream
2319          $out,    # handle to out stream
2320        ) = @_;
2321
2322     # Returns nothing.
2323
2324     my ( $record, $index );
2325
2326     while ( $record = get_record( $in ) )
2327     {
2328         $record->{ "SEQ_COMPLEXITY" } = sprintf( "%.2f", Maasha::Seq::seq_complexity( $record->{ "SEQ" } ) ) if $record->{ "SEQ" };
2329
2330         put_record( $record, $out );
2331     }
2332 }
2333
2334
2335 sub script_oligo_freq
2336 {
2337     # Martin A. Hansen, August 2007.
2338
2339     # Determine the length of sequences in stream.
2340
2341     my ( $in,        # handle to in stream
2342          $out,       # handle to out stream
2343          $options,   # options hash
2344        ) = @_;
2345
2346     # Returns nothing.
2347
2348     my ( $record, %oligos, @freq_table );
2349
2350     $options->{ "word_size" } ||= 7;
2351
2352     while ( $record = get_record( $in ) )
2353     {
2354         if ( $record->{ "SEQ" } )
2355         {
2356             map { $oligos{ $_ }++ } Maasha::Seq::seq2oligos( \$record->{ "SEQ" }, $options->{ "word_size" } );
2357
2358             if ( not $options->{ "all" } )
2359             {
2360                 @freq_table = Maasha::Seq::oligo_freq( \%oligos );
2361
2362                 map { put_record( $_, $out ) } @freq_table;
2363
2364                 undef %oligos;
2365             }
2366         }
2367
2368         put_record( $record, $out );
2369     }
2370
2371     if ( $options->{ "all" } )
2372     {
2373         @freq_table = Maasha::Seq::oligo_freq( \%oligos );
2374
2375         map { put_record( $_, $out ) } @freq_table;
2376     }
2377 }
2378
2379
2380 sub script_create_weight_matrix
2381 {
2382     # Martin A. Hansen, August 2007.
2383
2384     # Creates a weight matrix from an alignmnet.
2385
2386     my ( $in,        # handle to in stream
2387          $out,       # handle to out stream
2388          $options,   # options hash
2389        ) = @_;
2390
2391     # Returns nothing.
2392
2393     my ( $record, $count, $i, $res, %freq_hash, %res_hash, $freq );
2394
2395     $count = 0;
2396
2397     while ( $record = get_record( $in ) )
2398     {
2399         if ( $record->{ "SEQ" } )
2400         {
2401             for ( $i = 0; $i < length $record->{ "SEQ" }; $i++ )
2402             {
2403                 $res = substr $record->{ "SEQ" }, $i, 1;
2404
2405                 $freq_hash{ $i }{ $res }++;
2406                 $res_hash{ $res } = 1;
2407             }
2408
2409             $count++;
2410         }
2411         else
2412         {
2413             put_record( $record, $out );
2414         }
2415     }
2416
2417     foreach $res ( sort keys %res_hash )
2418     {
2419         undef $record;
2420
2421         $record->{ "V0" } = $res;
2422
2423         for ( $i = 0; $i < keys %freq_hash; $i++ )
2424         {
2425             $freq = $freq_hash{ $i }{ $res } || 0;
2426
2427             if ( $options->{ "percent" } ) {
2428                 $freq = sprintf( "%.0f", 100 * $freq / $count ) if $freq > 0;
2429             }
2430
2431             $record->{ "V" . ( $i + 1 ) } = $freq;
2432         }
2433
2434         put_record( $record, $out );
2435     }
2436 }
2437
2438
2439 sub script_calc_bit_scores
2440 {
2441     # Martin A. Hansen, March 2007.
2442
2443     # Calculates the bit scores for each position from an alignmnet in the stream.
2444
2445     my ( $in,        # handle to in stream
2446          $out,       # handle to out stream
2447        ) = @_;
2448
2449     # Returns nothing.
2450
2451     my ( $record, $type, $count, $i, $res, %freq_hash, $bit_max, $bit_height, $bit_diff );
2452
2453     $count = 0;
2454
2455     while ( $record = get_record( $in ) )
2456     {
2457         if ( $record->{ "SEQ" } )
2458         {
2459             $type = Maasha::Seq::seq_guess_type( $record->{ "SEQ" } ) if not $type;
2460
2461             for ( $i = 0; $i < length $record->{ "SEQ" }; $i++ )
2462             {
2463                 $res = substr $record->{ "SEQ" }, $i, 1;
2464
2465                 next if $res =~ /-|_|~|\./;
2466
2467                 $freq_hash{ $i }{ $res }++;
2468             }
2469
2470             $count++;
2471         }
2472         else
2473         {
2474             put_record( $record, $out );
2475         }
2476     }
2477
2478     undef $record;
2479
2480     if ( $type eq "protein" ) {
2481         $bit_max = 4;
2482     } else {
2483         $bit_max = 2;
2484     }
2485
2486     for ( $i = 0; $i < keys %freq_hash; $i++ )
2487     {
2488         $bit_height = Maasha::Seq::seqlogo_calc_bit_height( $freq_hash{ $i }, $count );
2489
2490         $bit_diff = $bit_max - $bit_height;
2491
2492         $record->{ "V" . ( $i ) } = sprintf( "%.2f", $bit_diff );
2493     }
2494
2495     put_record( $record, $out );
2496 }
2497
2498
2499 sub script_calc_fixedstep
2500 {
2501     # Martin A. Hansen, September 2008.
2502
2503     # Calculates fixedstep entries from data in the stream.
2504
2505     my ( $in,        # handle to in stream
2506          $out,       # handle to out stream
2507          $options,   # options hash
2508        ) = @_;
2509
2510     # Returns nothing.
2511
2512     my ( $record, %fh_hash, $fh_in, $fh_out, $chr, $chr, $beg, $end, $q_id, $block, $entry, $clones, $beg_block, $max, $i );
2513
2514     while ( $record = get_record( $in ) )
2515     {
2516         $record->{ "CHR" }     = $record->{ "S_ID" }  if not defined $record->{ "CHR" };
2517         $record->{ "CHR_BEG" } = $record->{ "S_BEG" } if not defined $record->{ "CHR_BEG" };
2518         $record->{ "CHR_END" } = $record->{ "S_END" } if not defined $record->{ "CHR_END" };
2519
2520         if ( $record->{ "CHR" } and defined $record->{ "CHR_BEG" } and $record->{ "CHR_END" } )
2521         {
2522             $fh_hash{ $record->{ "CHR" } } = Maasha::Common::write_open( "$BP_TMP/$record->{ 'CHR' }" ) if not exists $fh_hash{ $record->{ "CHR" } };
2523
2524             $fh_out = $fh_hash{ $record->{ "CHR" } };
2525
2526             Maasha::UCSC::bed_put_entry( $record, $fh_out, 5 );
2527         }
2528     }
2529
2530     map { close $_ } keys %fh_hash;
2531
2532     foreach $chr ( sort keys %fh_hash )
2533     {
2534         Maasha::Common::run( "bedSort", "$BP_TMP/$chr $BP_TMP/$chr" );
2535
2536         $fh_in = Maasha::Common::read_open( "$BP_TMP/$chr" );
2537
2538         undef $block;
2539
2540         while ( $entry = Maasha::UCSC::bed_get_entry( $fh_in, 5 ) )
2541         {
2542             $chr  = $entry->{ 'CHR' };
2543             $beg  = $entry->{ 'CHR_BEG' };
2544             $end  = $entry->{ 'CHR_END' };
2545             $q_id = $entry->{ 'Q_ID' };
2546
2547             if ( $options->{ "score" } ) {
2548                 $clones = $entry->{ 'SCORE' };
2549             } elsif ( $q_id =~ /_(\d+)$/ ) {
2550                 $clones = $1;
2551             } else {
2552                 $clones = 1;
2553             }
2554
2555             if ( $block )
2556             {
2557                 if ( $beg > $max )
2558                 {
2559                     map { $_ = sprintf( "%.4f", Maasha::Calc::log10( $_ ) ) } @{ $block } if $options->{ "log10" };
2560
2561                     $record->{ "CHR" }      = $chr;
2562                     $record->{ "CHR_BEG" }  = $beg_block;
2563                     $record->{ "STEP" }     = 1;
2564                     $record->{ "VALS" }     = join ";", @{ $block };
2565                     $record->{ "REC_TYPE" } = "fixed_step";
2566
2567                     put_record( $record, $out );
2568
2569                     undef $block;
2570                 }
2571                 else
2572                 {
2573                     for ( $i = $beg - $beg_block; $i < ( $beg - $beg_block ) + ( $end - $beg ); $i++ ) {
2574                         $block->[ $i ] += $clones;
2575                     }
2576
2577                     $max = Maasha::Calc::max( $max, $end );
2578                 }
2579             }
2580
2581             if ( not $block )
2582             {
2583                 $beg_block = $beg;
2584                 $max       = $end;
2585
2586                 for ( $i = 0; $i < ( $end - $beg ); $i++ ) {
2587                     $block->[ $i ] += $clones;
2588                 }
2589             }
2590         }
2591
2592         close $fh_in;
2593
2594         map { $_ = sprintf( "%.4f", Maasha::Calc::log10( $_ ) ) } @{ $block } if $options->{ "log10" };
2595
2596         $record->{ "CHR" }      = $chr;
2597         $record->{ "CHR_BEG" }  = $beg_block;
2598         $record->{ "STEP" }     = 1;
2599         $record->{ "VALS" }     = join ";", @{ $block };
2600         $record->{ "REC_TYPE" } = "fixed_step";
2601
2602         put_record( $record, $out );
2603
2604         unlink "$BP_TMP/$chr";
2605     }
2606 }
2607
2608
2609 sub script_reverse_seq
2610 {
2611     # Martin A. Hansen, August 2007.
2612
2613     # Reverse sequence in record.
2614
2615     my ( $in,    # handle to in stream
2616          $out,   # handle to out stream
2617        ) = @_;
2618
2619     # Returns nothing.
2620
2621     my ( $record );
2622
2623     while ( $record = get_record( $in ) )
2624     {
2625         if ( $record->{ "SEQ" } ) {
2626             $record->{ "SEQ" } = reverse $record->{ "SEQ" };
2627         }
2628
2629         put_record( $record, $out );
2630     }
2631 }
2632
2633
2634 sub script_complement_seq
2635 {
2636     # Martin A. Hansen, August 2007.
2637
2638     # Complement sequence in record.
2639
2640     my ( $in,     # handle to in stream
2641          $out,    # handle to out stream
2642        ) = @_;
2643
2644     # Returns nothing.
2645
2646     my ( $record, $type );
2647
2648     while ( $record = get_record( $in ) )
2649     {
2650         if ( $record->{ "SEQ" } )
2651         {
2652             if ( not $type ) {
2653                 $type = Maasha::Seq::seq_guess_type( $record->{ "SEQ" } );
2654             }
2655
2656             if ( $type eq "rna" ) {
2657                 Maasha::Seq::rna_comp( \$record->{ "SEQ" } );
2658             } elsif ( $type eq "dna" ) {
2659                 Maasha::Seq::dna_comp( \$record->{ "SEQ" } );
2660             }
2661         }
2662
2663         put_record( $record, $out );
2664     }
2665 }
2666
2667
2668 sub script_remove_indels
2669 {
2670     # Martin A. Hansen, August 2007.
2671
2672     # Remove indels from sequences in stream.
2673
2674     my ( $in,     # handle to in stream
2675          $out,    # handle to out stream
2676        ) = @_;
2677
2678     # Returns nothing.
2679
2680     my ( $record );
2681
2682     while ( $record = get_record( $in ) )
2683     {
2684         $record->{ 'SEQ' } =~ tr/-~.//d if $record->{ "SEQ" };
2685
2686         put_record( $record, $out );
2687     }
2688 }
2689
2690
2691 sub script_transliterate_seq
2692 {
2693     # Martin A. Hansen, August 2007.
2694
2695     # Transliterate chars from sequence in record.
2696
2697     my ( $in,        # handle to in stream
2698          $out,       # handle to out stream
2699          $options,   # options hash
2700        ) = @_;
2701
2702     # Returns nothing.
2703
2704     my ( $record, $search, $replace, $delete );
2705
2706     $search  = $options->{ "search" }  || "";
2707     $replace = $options->{ "replace" } || "";
2708     $delete  = $options->{ "delete" }  || "";
2709
2710     while ( $record = get_record( $in ) )
2711     {
2712         if ( $record->{ "SEQ" } )
2713         {
2714             if ( $search and $replace ) {
2715                 eval "\$record->{ 'SEQ' } =~ tr/$search/$replace/";
2716             } elsif ( $delete ) {
2717                 eval "\$record->{ 'SEQ' } =~ tr/$delete//d";
2718             }
2719         }
2720
2721         put_record( $record, $out );
2722     }
2723 }
2724
2725
2726 sub script_transliterate_vals
2727 {
2728     # Martin A. Hansen, April 2008.
2729
2730     # Transliterate chars from values in record.
2731
2732     my ( $in,        # handle to in stream
2733          $out,       # handle to out stream
2734          $options,   # options hash
2735        ) = @_;
2736
2737     # Returns nothing.
2738
2739     my ( $record, $search, $replace, $delete, $key );
2740
2741     $search  = $options->{ "search" }  || "";
2742     $replace = $options->{ "replace" } || "";
2743     $delete  = $options->{ "delete" }  || "";
2744
2745     while ( $record = get_record( $in ) )
2746     {
2747         foreach $key ( @{ $options->{ "keys" } } )
2748         {
2749             if ( exists $record->{ $key } )
2750             {
2751                 if ( $search and $replace ) {
2752                     eval "\$record->{ $key } =~ tr/$search/$replace/";
2753                 } elsif ( $delete ) {
2754                     eval "\$record->{ $key } =~ tr/$delete//d";
2755                 }
2756             }
2757         }
2758
2759         put_record( $record, $out );
2760     }
2761 }
2762
2763
2764 sub script_translate_seq
2765 {
2766     # Martin A. Hansen, February 2008.
2767
2768     # Translate DNA sequence into protein sequence.
2769
2770     my ( $in,        # handle to in stream
2771          $out,       # handle to out stream
2772          $options,   # options hash
2773        ) = @_;
2774
2775     # Returns nothing.
2776
2777     my ( $record, $frame, %new_record );
2778
2779     $options->{ "frames" } ||= [ 1, 2, 3, -1, -2, -3 ];
2780
2781     while ( $record = get_record( $in ) )
2782     {
2783         if ( $record->{ "SEQ" } )
2784         {
2785             if ( Maasha::Seq::seq_guess_type( $record->{ "SEQ" } ) eq "dna" )
2786             {
2787                 foreach $frame ( @{ $options->{ "frames" } } )
2788                 {
2789                     %new_record = %{ $record };
2790
2791                     $new_record{ "SEQ" }     = Maasha::Seq::translate( $record->{ "SEQ" }, $frame );
2792                     $new_record{ "SEQ_LEN" } = length $new_record{ "SEQ" };
2793                     $new_record{ "FRAME" }   = $frame;
2794
2795                     put_record( \%new_record, $out );
2796                 }
2797             }
2798         }
2799         else
2800         {
2801             put_record( $record, $out );
2802         }
2803     }
2804 }
2805
2806
2807 sub script_extract_seq
2808 {
2809     # Martin A. Hansen, August 2007.
2810
2811     # Extract subsequences from sequences in record.
2812
2813     my ( $in,        # handle to in stream
2814          $out,       # handle to out stream
2815          $options,   # options hash
2816        ) = @_;
2817
2818     # Returns nothing.
2819
2820     my ( $beg, $end, $len, $record );
2821
2822     if ( not defined $options->{ "beg" } or $options->{ "beg" } < 0 ) {
2823         $beg = 0;
2824     } else {
2825         $beg = $options->{ "beg" } - 1;   # correcting for start offset
2826     }
2827
2828     if ( defined $options->{ "end" } and $options->{ "end" } - 1 < $beg ) {
2829         $end = $beg - 1;
2830     } elsif ( defined $options->{ "end" } ) {
2831         $end = $options->{ "end" } - 1;   # correcting for start offset
2832     }
2833
2834     $len = $options->{ "len" };
2835
2836 #    print "beg->$beg,  end->$end,  len->$len\n";
2837
2838     while ( $record = get_record( $in ) )
2839     {
2840         if ( $record->{ "SEQ" } )
2841         {
2842             if ( defined $beg and defined $end )
2843             {
2844                 if ( $end - $beg + 1 > length $record->{ "SEQ" } ) {
2845                     $record->{ "SEQ" } = substr $record->{ "SEQ" }, $beg;
2846                 } else {
2847                     $record->{ "SEQ" } = substr $record->{ "SEQ" }, $beg, $end - $beg + 1;
2848                 }
2849             }
2850             elsif ( defined $beg and defined $len )
2851             {
2852                 if ( $len > length $record->{ "SEQ" } ) {
2853                     $record->{ "SEQ" } = substr $record->{ "SEQ" }, $beg;
2854                 } else {
2855                     $record->{ "SEQ" } = substr $record->{ "SEQ" }, $beg, $len;
2856                 }
2857             }
2858             elsif ( defined $beg )
2859             {
2860                 $record->{ "SEQ" } = substr $record->{ "SEQ" }, $beg;
2861             }
2862         }
2863
2864         $record->{ "SEQ_LEN" } = length $record->{ "SEQ" };
2865
2866         put_record( $record, $out );
2867     }
2868 }
2869
2870
2871 sub script_get_genome_seq
2872 {
2873     # Martin A. Hansen, December 2007.
2874
2875     # Gets a subsequence from a genome.
2876
2877     my ( $in,        # handle to in stream
2878          $out,       # handle to out stream
2879          $options,   # options hash
2880        ) = @_;
2881
2882     # Returns nothing.
2883
2884     my ( $record, $genome, $genome_file, $index_file, $index, $fh, $index_head, $index_beg, $index_len, $beg, $len, %lookup_hash, @begs, @lens, $i );
2885
2886     $options->{ "flank" } ||= 0;
2887
2888     if ( $options->{ "genome" } )
2889     {
2890         $genome      = $options->{ "genome" };
2891
2892         $genome_file = "$ENV{ 'BP_DATA' }/genomes/$genome/fasta/$genome.fna";
2893         $index_file  = "$ENV{ 'BP_DATA' }/genomes/$genome/fasta/$genome.index";
2894
2895         $fh          = Maasha::Common::read_open( $genome_file );
2896         $index       = Maasha::Fasta::index_retrieve( $index_file );
2897
2898         shift @{ $index }; # Get rid of the file size info
2899
2900         map { $lookup_hash{ $_->[ 0 ] } = [ $_->[ 1 ], $_->[ 2 ] ] } @{ $index };
2901
2902         if ( exists $lookup_hash{ $options->{ "chr" } } and defined $options->{ "beg" } and ( defined $options->{ "end" } or defined $options->{ "len" } ) )
2903         {
2904             ( $index_beg, $index_len ) = @{ $lookup_hash{ $options->{ "chr" } } };
2905
2906             $beg = $index_beg + $options->{ "beg" } - 1;
2907
2908             if ( $options->{ "len" } ) {
2909                 $len = $options->{ "len" };
2910             } elsif ( $options->{ "end" } ) {
2911                 $len = ( $options->{ "end" } - $options->{ "beg" } + 1 );
2912             }
2913
2914             $beg -= $options->{ "flank" };
2915             $len += 2 * $options->{ "flank" };
2916
2917             if ( $beg <= $index_beg )
2918             {
2919                 $len -= $index_beg - $beg;
2920                 $beg = $index_beg;
2921             }
2922
2923             $len = $index_beg + $index_len - $beg if $beg + $len > $index_beg + $index_len;
2924
2925             next if $beg > $index_beg + $index_len;
2926
2927             $record->{ "CHR" }     = $options->{ "chr" };
2928             $record->{ "CHR_BEG" } = $beg - $index_beg;
2929             $record->{ "CHR_END" } = $record->{ "CHR_BEG" } + $len - 1;
2930
2931             $record->{ "SEQ" }     = Maasha::Common::file_read( $fh, $beg, $len );
2932             $record->{ "SEQ_LEN" } = $len;
2933
2934             put_record( $record, $out );
2935         }
2936     }
2937
2938     while ( $record = get_record( $in ) )
2939     {
2940         if ( $options->{ "genome" } and not $record->{ "SEQ" } )
2941         {
2942             if ( $record->{ "REC_TYPE" } eq "BED" and exists $lookup_hash{ $record->{ "CHR" } } )
2943             {
2944                 ( $index_beg, $index_len ) = @{ $lookup_hash{ $record->{ "CHR" } } };
2945
2946                 $beg = $record->{ "CHR_BEG" } + $index_beg;
2947                 $len = $record->{ "CHR_END" } - $record->{ "CHR_BEG" } + 1;
2948             }
2949             elsif ( $record->{ "REC_TYPE" } eq "PSL" and exists $lookup_hash{ $record->{ "S_ID" } } )
2950             {
2951                 ( $index_beg, $index_len ) = @{ $lookup_hash{ $record->{ "S_ID" } } };
2952
2953                 $beg = $record->{ "S_BEG" } + $index_beg;
2954                 $len = $record->{ "S_END" } - $record->{ "S_BEG" } + 1;
2955             }
2956             elsif ( $record->{ "REC_TYPE" } eq "BLAST" and exists $lookup_hash{ $record->{ "S_ID" } } )
2957             {
2958                 ( $index_beg, $index_len ) = @{ $lookup_hash{ $record->{ "S_ID" } } };
2959
2960                 $beg = $record->{ "S_BEG" } + $index_beg;
2961                 $len = $record->{ "S_END" } - $record->{ "S_BEG" } + 1;
2962             }
2963
2964             $beg -= $options->{ "flank" };
2965             $len += 2 * $options->{ "flank" };
2966
2967             if ( $beg <= $index_beg )
2968             {
2969                 $len -= $index_beg - $beg;
2970                 $beg = $index_beg;
2971             }
2972
2973             $len = $index_beg + $index_len - $beg if $beg + $len > $index_beg + $index_len;
2974
2975             next if $beg > $index_beg + $index_len;
2976
2977             $record->{ "CHR_BEG" } = $beg - $index_beg;
2978             $record->{ "CHR_END" } = $record->{ "CHR_BEG" } + $len - 1;
2979
2980             $record->{ "SEQ" } = Maasha::Common::file_read( $fh, $beg, $len );
2981
2982             if ( $record->{ "STRAND" } and $record->{ "STRAND" } eq "-" )
2983             {
2984                 Maasha::Seq::dna_comp( \$record->{ "SEQ" } );
2985                 $record->{ "SEQ" } = reverse $record->{ "SEQ" };
2986             }
2987
2988             if ( $options->{ "mask" } )
2989             {
2990                 if ( $record->{ "BLOCKCOUNT" } > 1 ) # uppercase hit block segments and lowercase the rest.
2991                 {
2992                     $record->{ "SEQ" } = lc $record->{ "SEQ" };
2993
2994                     @begs = split ",", $record->{ "Q_BEGS" };
2995                     @lens = split ",", $record->{ "BLOCKSIZES" };
2996
2997                     for ( $i = 0; $i < @begs; $i++ ) {
2998                         substr $record->{ "SEQ" }, $begs[ $i ], $lens[ $i ], uc substr $record->{ "SEQ" }, $begs[ $i ], $lens[ $i ];
2999                     }
3000                 }
3001             }
3002         }
3003
3004         put_record( $record, $out );
3005     }
3006
3007     close $fh if $fh;
3008 }
3009
3010
3011 sub script_get_genome_align
3012 {
3013     # Martin A. Hansen, April 2008.
3014
3015     # Gets a subalignment from a multiple genome alignment.
3016
3017     my ( $in,        # handle to in stream
3018          $out,       # handle to out stream
3019          $options,   # options hash
3020        ) = @_;
3021
3022     # Returns nothing.
3023
3024     my ( $record, $maf_track, $align, $align_num, $beg, $end, $len, $entry );
3025
3026     $options->{ "strand" } ||= "+";
3027
3028     $align_num = 1;
3029
3030     $maf_track = Maasha::Config::maf_track( $options->{ "genome" } );
3031
3032     if ( $options->{ "chr" } and $options->{ "beg" } and ( $options->{ "end" } or $options->{ "len" } ) )
3033     {
3034         $beg = $options->{ "beg" } - 1;
3035
3036         if ( $options->{ "end" } ) {
3037             $end = $options->{ "end" };
3038         } elsif ( $options->{ "len" } ) {
3039             $end = $beg + $options->{ "len" };
3040         }
3041
3042         $align = Maasha::UCSC::maf_extract( $BP_TMP, $options->{ "genome" }, $maf_track, $options->{ "chr" }, $beg, $end, $options->{ "strand" } );
3043
3044         foreach $entry ( @{ $align } )
3045         {
3046             $entry->{ "CHR" }     = $record->{ "CHR" };
3047             $entry->{ "CHR_BEG" } = $record->{ "CHR_BEG" };
3048             $entry->{ "CHR_END" } = $record->{ "CHR_END" };
3049             $entry->{ "STRAND" }  = $record->{ "STRAND" } || '+';
3050             $entry->{ "Q_ID" }    = $record->{ "Q_ID" };
3051             $entry->{ "SCORE" }   = $record->{ "SCORE" };
3052
3053             put_record( $entry, $out );
3054         }
3055     }
3056
3057     while ( $record = get_record( $in ) )
3058     {
3059         if ( $record->{ "REC_TYPE" } eq "BED" )
3060         {
3061             $align = Maasha::UCSC::maf_extract( $BP_TMP, $options->{ "genome" }, $maf_track, $record->{ "CHR" }, $record->{ "CHR_BEG" }, $record->{ "CHR_END" }, $record->{ "STRAND" } );
3062         }
3063         elsif ( $record->{ "REC_TYPE" } eq "PSL" )
3064         {
3065             $align = Maasha::UCSC::maf_extract( $BP_TMP, $options->{ "genome" }, $maf_track, $record->{ "S_ID" }, $record->{ "S_BEG" }, $record->{ "S_END" }, $record->{ "STRAND" } );
3066         }
3067         elsif ( $record->{ "REC_TYPE" } eq "BLAST" )
3068         {
3069             $align = Maasha::UCSC::maf_extract( $BP_TMP, $options->{ "genome" }, $maf_track, $record->{ "S_ID" }, $record->{ "S_BEG" }, $record->{ "S_END" }, $record->{ "STRAND" } );
3070         }
3071
3072         foreach $entry ( @{ $align } )
3073         {
3074             $entry->{ "CHR" }     = $record->{ "CHR" };
3075             $entry->{ "CHR_BEG" } = $record->{ "CHR_BEG" };
3076             $entry->{ "CHR_END" } = $record->{ "CHR_END" };
3077             $entry->{ "STRAND" }  = $record->{ "STRAND" };
3078             $entry->{ "Q_ID" }    = $record->{ "Q_ID" };
3079             $entry->{ "SCORE" }   = $record->{ "SCORE" };
3080
3081             put_record( $entry, $out );
3082         }
3083
3084         $align_num++;
3085     }
3086 }
3087
3088
3089 sub script_get_genome_phastcons
3090 {
3091     # Martin A. Hansen, February 2008.
3092
3093     # Get phastcons scores from genome intervals.
3094
3095     my ( $in,        # handle to in stream
3096          $out,       # handle to out stream
3097          $options,   # options hash
3098        ) = @_;
3099
3100     # Returns nothing.
3101
3102     my ( $phastcons_file, $phastcons_index, $index, $fh_phastcons, $scores, $record );
3103
3104     $options->{ "flank" } ||= 0;
3105
3106     $phastcons_file  = Maasha::Config::genome_phastcons( $options->{ "genome" } );
3107     $phastcons_index = Maasha::Config::genome_phastcons_index( $options->{ "genome" } );
3108
3109     $index           = Maasha::UCSC::fixedstep_index_retrieve( $phastcons_index );
3110     $fh_phastcons    = Maasha::Common::read_open( $phastcons_file );
3111
3112     if ( defined $options->{ "chr" } and defined $options->{ "beg" } and ( defined $options->{ "end" } or defined $options->{ "len" } ) )
3113     {
3114         $options->{ "beg" } -= 1;   # request is 1-based
3115         $options->{ "end" } -= 1;   # request is 1-based
3116
3117         if ( $options->{ "len" } ) {
3118             $options->{ "end" } = $options->{ "beg" } + $options->{ "len" } - 1;
3119         }
3120
3121         $scores = Maasha::UCSC::fixedstep_index_lookup( $index, $fh_phastcons, $options->{ "chr" }, $options->{ "beg" }, $options->{ "end" }, $options->{ "flank" } );
3122
3123         $record->{ "CHR" }       = $options->{ "chr" };
3124         $record->{ "CHR_BEG" }   = $options->{ "beg" } - $options->{ "flank" };
3125         $record->{ "CHR_END" }   = $options->{ "end" } + $options->{ "flank" };
3126
3127         $record->{ "PHASTCONS" }   = join ",", @{ $scores };
3128         $record->{ "PHAST_COUNT" } = scalar @{ $scores };  # DEBUG
3129
3130         put_record( $record, $out );
3131     }
3132
3133     while ( $record = get_record( $in ) )
3134     {
3135         if ( $record->{ "REC_TYPE" } eq "BED" )
3136         {
3137             $scores = Maasha::UCSC::fixedstep_index_lookup( $index, $fh_phastcons, $record->{ "CHR" }, $record->{ "CHR_BEG" }, $record->{ "CHR_END" }, $options->{ "flank" } );
3138         }
3139         elsif ( $record->{ "REC_TYPE" } eq "PSL" )
3140         {
3141             $scores = Maasha::UCSC::fixedstep_index_lookup( $index, $fh_phastcons, $record->{ "S_ID" }, $record->{ "S_BEG" }, $record->{ "S_END" }, $options->{ "flank" } );
3142         }
3143         elsif ( $record->{ "REC_TYPE" } eq "BLAST" )
3144         {
3145             $scores = Maasha::UCSC::fixedstep_index_lookup( $index, $fh_phastcons, $record->{ "S_ID" }, $record->{ "S_BEG" }, $record->{ "S_END" }, $options->{ "flank" } );
3146         }
3147
3148         $record->{ "PHASTCONS" } = join ",", @{ $scores } if @{ $scores };
3149 #        $record->{ "PHAST_COUNT" } = @{ $scores } if @{ $scores };  # DEBUG
3150
3151         put_record( $record, $out );
3152     }
3153
3154     close $fh_phastcons if $fh_phastcons;
3155 }
3156
3157
3158 sub script_fold_seq
3159 {
3160     # Martin A. Hansen, December 2007.
3161
3162     # Folds sequences in stream into secondary structures.
3163
3164     my ( $in,     # handle to in stream
3165          $out,    # handle to out stream
3166        ) = @_;
3167
3168     # Returns nothing.
3169
3170     my ( $record, $type, $struct, $index );
3171
3172     while ( $record = get_record( $in ) )
3173     {
3174         if ( $record->{ "SEQ" } )
3175         {
3176             if ( not $type ) {
3177                 $type = Maasha::Seq::seq_guess_type( $record->{ "SEQ" } );
3178             }
3179
3180             if ( $type ne "protein" )
3181             {
3182                 ( $struct, $index ) = Maasha::Seq::fold_struct_rnafold( $record->{ "SEQ" } );
3183                 $record->{ "SEC_STRUCT" }  = $struct;
3184                 $record->{ "FREE_ENERGY" } = $index;
3185                 $record->{ "SCORE" }       = abs int $index;
3186                 $record->{ "SIZE" }        = length $struct;
3187                 $record->{ "CONF" }        = "1," x $record->{ "SIZE" };
3188             }
3189         }
3190
3191         put_record( $record, $out );
3192     }
3193 }
3194
3195
3196 sub script_split_seq
3197 {
3198     # Martin A. Hansen, August 2007.
3199
3200     # Split a sequence in stream into words.
3201
3202     my ( $in,        # handle to in stream
3203          $out,       # handle to out stream
3204          $options,   # options hash
3205        ) = @_;
3206
3207     # Returns nothing.
3208
3209     my ( $record, $new_record, $i, $subseq, %lookup );
3210
3211     $options->{ "word_size" } ||= 7;
3212
3213     while ( $record = get_record( $in ) )
3214     {
3215         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } )
3216         {
3217             for ( $i = 0; $i < length( $record->{ "SEQ" } ) - $options->{ "word_size" } + 1; $i++ )
3218             {
3219                 $subseq = substr $record->{ "SEQ" }, $i, $options->{ "word_size" };
3220
3221                 if ( $options->{ "uniq" } and not $lookup{ $subseq } )
3222                 {
3223                     $new_record->{ "SEQ_NAME" } = $record->{ "SEQ_NAME" } . "[" . ( $i + 1 ) . "-" . ( $i + $options->{ "word_size" } ) . "]";
3224                     $new_record->{ "SEQ" }      = $subseq;
3225
3226                     put_record( $new_record, $out );
3227
3228                     $lookup{ $subseq } = 1;
3229                 }
3230                 else
3231                 {
3232                     $new_record->{ "SEQ_NAME" } = $record->{ "SEQ_NAME" } . "[" . ( $i + 1 ) . "-" . ( $i + $options->{ "word_size" } ) . "]";
3233                     $new_record->{ "SEQ" }      = $subseq;
3234
3235                     put_record( $new_record, $out );
3236                 }
3237             }
3238         }
3239         else
3240         {
3241             put_record( $record, $out );
3242         }
3243     }
3244 }
3245
3246
3247 sub script_split_bed
3248 {
3249     # Martin A. Hansen, June 2008.
3250
3251     # Split a BED record into overlapping windows.
3252
3253     my ( $in,        # handle to in stream
3254          $out,       # handle to out stream
3255          $options,   # options hash
3256        ) = @_;
3257
3258     # Returns nothing.
3259
3260     my ( $record, $new_record, $i );
3261
3262     $options->{ "window_size" } ||= 20;
3263     $options->{ "step_size" }   ||= 1;
3264
3265     while ( $record = get_record( $in ) )
3266     {
3267         if ( $record->{ "CHR" } and $record->{ "CHR_BEG" } and $record->{ "CHR_END" } )
3268         {
3269             $record->{ "BED_LEN" } = $record->{ "CHR_END" } - $record->{ "CHR_BEG" } + 1;
3270
3271             for ( $i = 0; $i < $record->{ "BED_LEN" } - $options->{ "window_size" }; $i += $options->{ "step_size" } )
3272             {
3273                 $new_record->{ "REC_TYPE" } = "BED";
3274                 $new_record->{ "CHR" }      = $record->{ "CHR" };
3275                 $new_record->{ "CHR_BEG" }  = $record->{ "CHR_BEG" } + $i;
3276                 $new_record->{ "CHR_END" }  = $record->{ "CHR_BEG" } + $i + $options->{ "window_size" };
3277                 $new_record->{ "BED_LEN" }  = $options->{ "window_size" };
3278                 $new_record->{ "Q_ID" }     = $record->{ "Q_ID" } . "_$i";
3279                 $new_record->{ "SCORE" }    = $record->{ "SCORE" };
3280                 $new_record->{ "STRAND" }   = $record->{ "STRAND" };
3281
3282                 put_record( $new_record, $out );
3283             }
3284         }
3285         else
3286         {
3287             put_record( $record, $out );
3288         }
3289     }
3290 }
3291
3292
3293 sub script_align_seq
3294 {
3295     # Martin A. Hansen, August 2007.
3296
3297     # Align sequences in stream.
3298
3299     my ( $in,    # handle to in stream
3300          $out,   # handle to out stream
3301        ) = @_;
3302
3303     # Returns nothing.
3304
3305     my ( $record, @entries, $entry );
3306
3307     while ( $record = get_record( $in ) )
3308     {
3309         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } ) {
3310             push @entries, [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ];
3311         } elsif ( $record->{ "Q_ID" } and $record->{ "SEQ" } ) {
3312             push @entries, [ $record->{ "Q_ID" }, $record->{ "SEQ" } ];
3313         } else {
3314             put_record( $record, $out );
3315         }
3316     }
3317
3318     @entries = Maasha::Align::align( \@entries );
3319
3320     foreach $entry ( @entries )
3321     {
3322         if ( $entry->[ SEQ_NAME ] and $entry->[ SEQ ] )
3323         {
3324             $record = {
3325                 SEQ_NAME => $entry->[ SEQ_NAME ],
3326                 SEQ      => $entry->[ SEQ ],
3327             };
3328
3329             put_record( $record, $out );
3330         }
3331     }
3332 }
3333
3334
3335 sub script_tile_seq
3336 {
3337     # Martin A. Hansen, February 2008.
3338
3339     # Using the first sequence in stream as reference, tile
3340     # all subsequent sequences based on pairwise alignments.
3341
3342     my ( $in,        # handle to in stream
3343          $out,       # handle to out stream
3344          $options,   # options hash
3345        ) = @_;
3346
3347     # Returns nothing.
3348
3349     my ( $record, $first, $ref_entry, @entries );
3350
3351     $first = 1;
3352
3353     while ( $record = get_record( $in ) )
3354     {
3355         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } )
3356         {
3357             if ( $first )
3358             {
3359                 $ref_entry = [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ];
3360
3361                 $first = 0;
3362             }
3363             else
3364             {
3365                 push @entries, [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ];
3366             }
3367         }
3368         else
3369         {
3370             put_record( $record, $out );
3371         }
3372     }
3373
3374     @entries = Maasha::Align::align_tile( $ref_entry, \@entries, $options );
3375
3376     map { put_record( { SEQ_NAME => $_->[ SEQ_NAME ], SEQ => $_->[ SEQ ] }, $out ) } @entries;
3377 }
3378
3379
3380 sub script_invert_align
3381 {
3382     # Martin A. Hansen, February 2008.
3383
3384     # Inverts an alignment showing only non-mathing residues
3385     # using the first sequence as reference.
3386
3387     my ( $in,        # handle to in stream
3388          $out,       # handle to out stream
3389          $options,   # options hash
3390        ) = @_;
3391
3392     # Returns nothing.
3393
3394     my ( $record, @entries );
3395
3396     while ( $record = get_record( $in ) )
3397     {
3398         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } )
3399         {
3400             push @entries, [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ];
3401         }
3402         else
3403         {
3404             put_record( $record, $out );
3405         }
3406     }
3407
3408     Maasha::Align::align_invert( \@entries, $options->{ "soft" } );
3409
3410     map { put_record( { SEQ_NAME => $_->[ SEQ_NAME ], SEQ => $_->[ SEQ ] }, $out ) } @entries;
3411 }
3412
3413
3414 sub script_patscan_seq
3415 {
3416     # Martin A. Hansen, August 2007.
3417
3418     # Locates patterns in sequences using scan_for_matches.
3419
3420     my ( $in,        # handle to in stream
3421          $out,       # handle to out stream
3422          $options,   # options hash
3423        ) = @_;
3424
3425     # Returns nothing.
3426
3427     my ( $genome_file, @args, $arg, $type, $seq_file, $pat_file, $out_file, $fh_in, $fh_out, $record, $patterns, $pattern, $entry, $result, %head_hash, $i );
3428
3429     if ( $options->{ "patterns" } ) {
3430         $patterns = Maasha::Patscan::parse_patterns( $options->{ "patterns" } );
3431     } elsif ( -f $options->{ "patterns_in" } ) {
3432         $patterns = Maasha::Patscan::read_patterns( $options->{ "patterns_in" } );
3433     }
3434
3435     $genome_file = "$ENV{ 'BP_DATA' }/genomes/$options->{ 'genome' }/fasta/$options->{ 'genome' }.fna" if $options->{ 'genome' };
3436
3437     push @args, "-c"                            if $options->{ "comp" };
3438     push @args, "-m $options->{ 'max_hits' }"   if $options->{ 'max_hits' };
3439     push @args, "-n $options->{ 'max_misses' }" if $options->{ 'max_hits' };
3440
3441     $seq_file = "$BP_TMP/patscan.seq";
3442     $pat_file = "$BP_TMP/patscan.pat";
3443     $out_file = "$BP_TMP/patscan.out";
3444
3445     $fh_out = Maasha::Common::write_open( $seq_file );
3446
3447     $i = 0;
3448
3449     while ( $record = get_record( $in ) )
3450     {
3451         if ( $record->{ "SEQ" } and $record->{ "SEQ_NAME" } )
3452         {
3453             $type = Maasha::Seq::seq_guess_type( $record->{ "SEQ" } ) if not $type;
3454
3455             Maasha::Fasta::put_entry( [ $i, $record->{ "SEQ" } ], $fh_out );
3456
3457             $head_hash{ $i } = $record->{ "SEQ_NAME" };
3458
3459             $i++;
3460         }
3461     }
3462
3463     close $fh_out;
3464
3465     $arg  = join " ", @args;
3466     $arg .= " -p" if $type eq "protein";
3467
3468     foreach $pattern ( @{ $patterns } )
3469     {
3470         $fh_out = Maasha::Common::write_open( $pat_file );
3471
3472         print $fh_out "$pattern\n";
3473
3474         close $fh_out;
3475
3476         if ( $options->{ 'genome' } ) {
3477             `scan_for_matches $arg $pat_file < $genome_file > $out_file`;
3478             # Maasha::Common::run( "scan_for_matches", "$arg $pat_file < $genome_file > $out_file" );
3479         } else {
3480             `scan_for_matches $arg $pat_file < $seq_file > $out_file`;
3481             # Maasha::Common::run( "scan_for_matches", "$arg $pat_file < $seq_file > $out_file" );
3482         }
3483
3484         $fh_in = Maasha::Common::read_open( $out_file );
3485
3486         while ( $entry = Maasha::Fasta::get_entry( $fh_in ) )
3487         {
3488             $result = Maasha::Patscan::parse_scan_result( $entry, $pattern );
3489
3490             if ( $options->{ 'genome' } )
3491             {
3492                 $result->{ "CHR" }     = $result->{ "S_ID" };
3493                 $result->{ "CHR_BEG" } = $result->{ "S_BEG" };
3494                 $result->{ "CHR_END" } = $result->{ "S_END" };
3495
3496                 delete $result->{ "S_ID" };
3497                 delete $result->{ "S_BEG" };
3498                 delete $result->{ "S_END" };
3499             }
3500             else
3501             {
3502                 $result->{ "S_ID" } = $head_hash{ $result->{ "S_ID" } };
3503             }
3504
3505             put_record( $result, $out );
3506         }
3507
3508         close $fh_in;
3509     }
3510
3511     unlink $pat_file;
3512     unlink $seq_file;
3513     unlink $out_file;
3514 }
3515
3516
3517 sub script_create_blast_db
3518 {
3519     # Martin A. Hansen, September 2007.
3520
3521     # Creates a NCBI BLAST database with formatdb
3522
3523     my ( $in,        # handle to in stream
3524          $out,       # handle to out stream
3525          $options,   # options hash
3526        ) = @_;
3527
3528     # Returns nothing.
3529
3530     my ( $fh, $seq_type, $path, $record, $entry );
3531
3532     $path = $options->{ "database" };
3533
3534     $fh = Maasha::Common::write_open( $path );
3535
3536     while ( $record = get_record( $in ) )
3537     {
3538         put_record( $record, $out ) if not $options->{ "no_stream" };
3539
3540         if ( $entry = record2fasta( $record ) )
3541         {
3542             $seq_type = Maasha::Seq::seq_guess_type( $entry->[ SEQ ] ) if not $seq_type;
3543
3544             Maasha::Fasta::put_entry( $entry, $fh );
3545         }
3546     }
3547
3548     close $fh;
3549
3550     if ( $seq_type eq "protein" ) {
3551         Maasha::Common::run( "formatdb", "-p T -i $path -t $options->{ 'database' }" );
3552     } else {
3553         Maasha::Common::run( "formatdb", "-p F -i $path -t $options->{ 'database' }" );
3554     }
3555
3556     unlink $path;
3557 }
3558
3559
3560 sub script_blast_seq
3561 {
3562     # Martin A. Hansen, September 2007.
3563
3564     # BLASTs sequences in stream against a given database.
3565
3566     my ( $in,        # handle to in stream
3567          $out,       # handle to out stream
3568          $options,   # options hash
3569        ) = @_;
3570
3571     # Returns nothing.
3572
3573     my ( $genome, $q_type, $s_type, $tmp_in, $tmp_out, $fh_in, $fh_out, $record, $line, @fields, $entry );
3574
3575     $options->{ "e_val" }  = 10 if not defined $options->{ "e_val" };
3576     $options->{ "filter" } = "F";
3577     $options->{ "filter" } = "T" if $options->{ "filter" };
3578     $options->{ "cpus" } ||= 1;
3579
3580     $options->{ "database" } = "$ENV{ 'BP_DATA' }/genomes/$options->{ 'genome' }/blast/$options->{ 'genome' }.fna" if $options->{ 'genome' };
3581
3582     $tmp_in  = "$BP_TMP/blast_query.seq";
3583     $tmp_out = "$BP_TMP/blast.result";
3584
3585     $fh_out = Maasha::Common::write_open( $tmp_in );
3586
3587     while ( $record = get_record( $in ) )
3588     {
3589         if ( $entry = record2fasta( $record ) )
3590         {
3591             $q_type = Maasha::Seq::seq_guess_type( $entry->[ SEQ ] ) if not $q_type;
3592
3593             Maasha::Fasta::put_entry( $entry, $fh_out );
3594         }
3595
3596         put_record( $record, $out );
3597     }
3598
3599     close $fh_out;
3600
3601     if ( -f $options->{ 'database' } . ".phr" ) {
3602         $s_type = "protein";
3603     } else {
3604         $s_type = "nucleotide";
3605     }
3606
3607     if ( not $options->{ 'program' } )
3608     {
3609         if ( $q_type ne "protein" and $s_type ne "protein" ) {
3610             $options->{ 'program' } = "blastn";
3611         } elsif ( $q_type eq "protein" and $s_type eq "protein" ) {
3612             $options->{ 'program' } = "blastp";
3613         } elsif ( $q_type ne "protein" and $s_type eq "protein" ) {
3614             $options->{ 'program' } = "blastx";
3615         } elsif ( $q_type eq "protein" and $s_type ne "protein" ) {
3616             $options->{ 'program' } = "tblastn";
3617         }
3618     }
3619
3620     if ( $options->{ 'verbose' } )
3621     {
3622         Maasha::Common::run(
3623             "blastall",
3624             join( " ",
3625                 "-p $options->{ 'program' }",
3626                 "-e $options->{ 'e_val' }",
3627                 "-a $options->{ 'cpus' }",
3628                 "-m 8",
3629                 "-i $tmp_in",
3630                 "-d $options->{ 'database' }",
3631                 "-F $options->{ 'filter' }",
3632                 "-o $tmp_out",
3633             ),
3634             1
3635         );
3636     }
3637     else
3638     {
3639         Maasha::Common::run(
3640             "blastall",
3641             join( " ",
3642                 "-p $options->{ 'program' }",
3643                 "-e $options->{ 'e_val' }",
3644                 "-a $options->{ 'cpus' }",
3645                 "-m 8",
3646                 "-i $tmp_in",
3647                 "-d $options->{ 'database' }",
3648                 "-F $options->{ 'filter' }",
3649                 "-o $tmp_out",
3650                 "> /dev/null 2>&1"
3651             ),
3652             1
3653         );
3654     }
3655
3656     unlink $tmp_in;
3657
3658     $fh_out = Maasha::Common::read_open( $tmp_out );
3659
3660     undef $record;
3661
3662     while ( $line = <$fh_out> )
3663     {
3664         chomp $line;
3665
3666         next if $line =~ /^#/;
3667
3668         @fields = split /\s+/, $line;
3669
3670         $record->{ "REC_TYPE" }   = "BLAST";
3671         $record->{ "Q_ID" }       = $fields[ 0 ];
3672         $record->{ "S_ID" }       = $fields[ 1 ];
3673         $record->{ "IDENT" }      = $fields[ 2 ];
3674         $record->{ "ALIGN_LEN" }  = $fields[ 3 ];
3675         $record->{ "MISMATCHES" } = $fields[ 4 ];
3676         $record->{ "GAPS" }       = $fields[ 5 ];
3677         $record->{ "Q_BEG" }      = $fields[ 6 ] - 1; # BLAST is 1-based
3678         $record->{ "Q_END" }      = $fields[ 7 ] - 1; # BLAST is 1-based
3679         $record->{ "S_BEG" }      = $fields[ 8 ] - 1; # BLAST is 1-based
3680         $record->{ "S_END" }      = $fields[ 9 ] - 1; # BLAST is 1-based
3681         $record->{ "E_VAL" }      = $fields[ 10 ];
3682         $record->{ "BIT_SCORE" }  = $fields[ 11 ];
3683
3684         if ( $record->{ "S_BEG" } > $record->{ "S_END" } )
3685         {
3686             $record->{ "STRAND" } = '-';
3687
3688             ( $record->{ "S_BEG" }, $record->{ "S_END" } ) = ( $record->{ "S_END" }, $record->{ "S_BEG" } );
3689         }
3690         else
3691         {
3692             $record->{ "STRAND" } = '+';
3693         }
3694
3695         put_record( $record, $out );
3696     }
3697
3698     close $fh_out;
3699
3700     unlink $tmp_out;
3701 }
3702
3703
3704 sub script_blat_seq
3705 {
3706     # Martin A. Hansen, August 2007.
3707
3708     # BLATs sequences in stream against a given genome.
3709
3710     my ( $in,        # handle to in stream
3711          $out,       # handle to out stream
3712          $options,   # options hash
3713        ) = @_;
3714
3715     # Returns nothing.
3716
3717     my ( $blat_args, $genome_file, $query_file, $fh_in, $fh_out, $type, $record, $result_file, $entries, $entry );
3718
3719     $genome_file = "$ENV{ 'BP_DATA' }/genomes/$options->{ 'genome' }/fasta/$options->{ 'genome' }.fna";
3720
3721     $options->{ 'tile_size' }    ||= 11;
3722     $options->{ 'one_off' }      ||= 0;
3723     $options->{ 'min_identity' } ||= 90;
3724     $options->{ 'min_score' }    ||= 0;
3725     $options->{ 'step_size' }    ||= $options->{ 'tile_size' };
3726
3727     $blat_args .= " -tileSize=$options->{ 'tile_size' }";
3728     $blat_args .= " -oneOff=$options->{ 'one_off' }";
3729     $blat_args .= " -minIdentity=$options->{ 'min_identity' }";
3730     $blat_args .= " -minScore=$options->{ 'min_score' }";
3731     $blat_args .= " -stepSize=$options->{ 'step_size' }";
3732 #    $blat_args .= " -ooc=" . Maasha::Config::genome_blat_ooc( $options->{ "genome" }, 11 ) if $options->{ 'ooc' };
3733
3734     $query_file = "$BP_TMP/blat.seq";
3735
3736     $fh_out = Maasha::Common::write_open( $query_file );
3737
3738     while ( $record = get_record( $in ) )
3739     {
3740         if ( $entry = record2fasta( $record ) )
3741         {
3742             $type = Maasha::Seq::seq_guess_type( $entry->[ SEQ ] ) if not $type;
3743             Maasha::Fasta::put_entry( $entry, $fh_out, 80 );
3744         }
3745
3746         put_record( $record, $out );
3747     }
3748
3749     close $fh_out;
3750
3751     $blat_args .= " -t=dnax" if $type eq "protein";
3752     $blat_args .= " -q=$type";
3753
3754     $result_file = "$BP_TMP/blat.psl";
3755
3756     Maasha::Common::run( "blat", "$genome_file $query_file $blat_args $result_file > /dev/null 2>&1" );
3757
3758     unlink $query_file;
3759
3760     $entries = Maasha::UCSC::psl_get_entries( $result_file );
3761
3762     map { put_record( $_, $out ) } @{ $entries };
3763
3764     unlink $result_file;
3765 }
3766
3767
3768 sub script_soap_seq
3769 {
3770     # Martin A. Hansen, July 2008.
3771
3772     # soap sequences in stream against a given file or genome.
3773
3774     my ( $in,        # handle to in stream
3775          $out,       # handle to out stream
3776          $options,   # options hash
3777        ) = @_;
3778
3779     # Returns nothing.
3780
3781     my ( $genome, $tmp_in, $tmp_out, $fh_in, $fh_out, $record, $line, @fields, $entry, $count, $args );
3782
3783     $options->{ "seed_size" }  ||= 10;
3784     $options->{ "mismatches" } ||= 2;
3785     $options->{ "gap_size" }   ||= 0;
3786     $options->{ "cpus" }       ||= 1;
3787
3788     if ( $options->{ "genome" } ) {
3789         $options->{ "in_file" } = "$ENV{ 'BP_DATA' }/genomes/$options->{ 'genome' }/fasta/$options->{ 'genome' }.fna";
3790     }
3791
3792     $tmp_in  = "$BP_TMP/soap_query.seq";
3793     $tmp_out = "$BP_TMP/soap.result";
3794
3795     $fh_out = Maasha::Common::write_open( $tmp_in );
3796
3797     $count = 0;
3798
3799     while ( $record = get_record( $in ) )
3800     {
3801         if ( $entry = record2fasta( $record ) )
3802         {
3803             Maasha::Fasta::put_entry( $entry, $fh_out );
3804
3805             $count++;
3806         }
3807
3808         put_record( $record, $out );
3809     }
3810
3811     close $fh_out;
3812
3813     if ( $count > 0 )
3814     {
3815         $args = join( " ",
3816             "-s $options->{ 'seed_size' }",
3817             "-r 2",
3818             "-a $tmp_in",
3819             "-v $options->{ 'mismatches' }",
3820             "-g $options->{ 'gap_size' }",
3821             "-p $options->{ 'cpus' }",
3822             "-d $options->{ 'in_file' }",
3823             "-o $tmp_out",
3824         );
3825
3826         $args .= " > /dev/null 2>&1" if not $options->{ 'verbose' };
3827
3828         Maasha::Common::run( "soap", $args, 1 );
3829
3830         unlink $tmp_in;
3831
3832         $fh_out = Maasha::Common::read_open( $tmp_out );
3833
3834         undef $record;
3835
3836         while ( $line = <$fh_out> )
3837         {
3838             chomp $line;
3839
3840             @fields = split /\t/, $line;
3841
3842             $record->{ "REC_TYPE" }   = "SOAP";
3843             $record->{ "Q_ID" }       = $fields[ 0 ];
3844             $record->{ "SCORE" }      = $fields[ 3 ];
3845             $record->{ "STRAND" }     = $fields[ 6 ];
3846             $record->{ "S_ID" }       = $fields[ 7 ];
3847             $record->{ "S_BEG" }      = $fields[ 8 ] - 1; # soap is 1-based
3848             $record->{ "S_END" }      = $fields[ 8 ] + $fields[ 5 ] - 2;
3849
3850             put_record( $record, $out );
3851         }
3852
3853         close $fh_out;
3854     }
3855
3856     unlink $tmp_out;
3857 }
3858
3859
3860 sub script_match_seq
3861 {
3862     # Martin A. Hansen, August 2007.
3863
3864     # BLATs sequences in stream against a given genome.
3865
3866     my ( $in,        # handle to in stream
3867          $out,       # handle to out stream
3868          $options,   # options hash
3869        ) = @_;
3870
3871     # Returns nothing.
3872
3873     my ( $record, @entries, $results );
3874
3875     $options->{ "word_size" } ||= 20;
3876     $options->{ "direction" } ||= "both";
3877
3878     while ( $record = get_record( $in ) )
3879     {
3880         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } ) {
3881             push @entries, [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ];
3882         }
3883
3884         put_record( $record, $out );
3885     }
3886
3887     if ( @entries == 1 )
3888     {
3889         $results = Maasha::Match::match_mummer( [ $entries[ 0 ] ], [ $entries[ 0 ] ], $options, $BP_TMP );
3890
3891         map { put_record( $_, $out ) } @{ $results };
3892     }
3893     elsif ( @entries == 2 )
3894     {
3895         $results = Maasha::Match::match_mummer( [ $entries[ 0 ] ], [ $entries[ 1 ] ], $options, $BP_TMP );
3896
3897         map { put_record( $_, $out ) } @{ $results };
3898     }
3899 }
3900
3901
3902 sub script_create_vmatch_index
3903 {
3904     # Martin A. Hansen, January 2008.
3905
3906     # Create a vmatch index from sequences in the stream.
3907
3908     my ( $in,        # handle to in stream
3909          $out,       # handle to out stream
3910          $options,   # options hash
3911        ) = @_;
3912
3913     # Returns nothing.
3914
3915     my ( $record, $file_tmp, $fh_tmp, $type, $entry );
3916
3917     if ( $options->{ "index_name" } )
3918     {
3919         $file_tmp = $options->{ 'index_name' };
3920         $fh_tmp   = Maasha::Common::write_open( $file_tmp );
3921     }
3922
3923     while ( $record = get_record( $in ) )
3924     {
3925         if ( $options->{ "index_name" } and $entry = record2fasta( $record ) )
3926         {
3927             Maasha::Fasta::put_entry( $entry, $fh_tmp );
3928
3929             $type = Maasha::Seq::seq_guess_type( $entry->[ SEQ ] ) if not defined $type;
3930         }
3931
3932         put_record( $record, $out ) if not $options->{ "no_stream" };
3933     }
3934
3935     if ( $options->{ "index_name" } )
3936     {
3937         close $fh_tmp;
3938
3939         if ( $type eq "protein" ) {
3940             Maasha::Common::run( "mkvtree", "-db $file_tmp -protein -pl $options->{ 'prefix_length' } -allout -indexname $file_tmp > /dev/null 2>&1" );
3941         } else {
3942             Maasha::Common::run( "mkvtree", "-db $file_tmp -dna -pl $options->{ 'prefix_length' } -allout -indexname $file_tmp > /dev/null 2>&1" );
3943         }
3944
3945         unlink $file_tmp;
3946     }
3947 }
3948
3949
3950 sub script_vmatch_seq
3951 {
3952     # Martin A. Hansen, August 2007.
3953
3954     # Vmatches sequences in stream against a given genome.
3955
3956     my ( $in,        # handle to in stream
3957          $out,       # handle to out stream
3958          $options,   # options hash
3959        ) = @_;
3960
3961     # Returns nothing.
3962
3963     my ( @index_files, @records, $result_file, $fh_in, $record, %hash );
3964
3965     $options->{ 'count' } = 1 if $options->{ 'max_hits' };
3966
3967     if ( $options->{ "index_name" } )
3968     {
3969         @index_files = $options->{ "index_name" };
3970     }
3971     else
3972     {
3973         @index_files = Maasha::Common::ls_files( "$ENV{ 'BP_DATA' }/genomes/$options->{ 'genome' }/vmatch" );
3974
3975         map { $_ =~ /^(.+)\.[a-z1]{3,4}$/; $hash{ $1 } = 1 } @index_files;
3976
3977         @index_files = sort keys %hash;
3978     }
3979
3980     while ( $record = get_record( $in ) )
3981     {
3982         push @records, $record;
3983
3984         put_record( $record, $out );
3985     }
3986
3987     $result_file = Maasha::Match::match_vmatch( $BP_TMP, \@records, \@index_files, $options );
3988
3989     undef @records;
3990
3991     $fh_in = Maasha::Common::read_open( $result_file );
3992
3993     while ( $record = Maasha::Match::vmatch_get_entry( $fh_in ) ) {
3994         put_record( $record, $out );
3995     }
3996
3997     close $fh_in;
3998
3999     unlink $result_file;
4000 }
4001
4002
4003 sub script_write_fasta
4004 {
4005     # Martin A. Hansen, August 2007.
4006
4007     # Write FASTA entries from sequences in stream.
4008
4009     my ( $in,        # handle to in stream
4010          $out,       # handle to out stream
4011          $options,   # options hash
4012        ) = @_;
4013
4014     # Returns nothing.
4015
4016     my ( $record, $fh, $entry );
4017
4018     $fh = write_stream( $options->{ "data_out" }, $options->{ "compress" } );
4019
4020     while ( $record = get_record( $in ) )
4021     {
4022         if ( $entry = record2fasta( $record ) ) {
4023             Maasha::Fasta::put_entry( $entry, $fh, $options->{ "wrap" } );
4024         }
4025
4026         put_record( $record, $out ) if not $options->{ "no_stream" };
4027     }
4028
4029     close $fh;
4030 }
4031
4032
4033 sub script_write_align
4034 {
4035     # Martin A. Hansen, August 2007.
4036
4037     # Write pretty alignments aligned sequences in stream.
4038
4039     my ( $in,        # handle to in stream
4040          $out,       # handle to out stream
4041          $options,   # options hash
4042        ) = @_;
4043
4044     # Returns nothing.
4045
4046     my ( $fh, $record, @entries );
4047
4048     $fh = write_stream( $options->{ "data_out" } ) ;
4049
4050     while ( $record = get_record( $in ) )
4051     {
4052         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } ) {
4053             push @entries, [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ];
4054         }
4055
4056         put_record( $record, $out ) if not $options->{ "no_stream" };
4057     }
4058
4059     if ( scalar( @entries ) == 2 ) {
4060         Maasha::Align::align_print_pairwise( $entries[ 0 ], $entries[ 1 ], $fh, $options->{ "wrap" } );
4061     } elsif ( scalar ( @entries ) > 2 ) {
4062         Maasha::Align::align_print_multi( \@entries, $fh, $options->{ "wrap" }, $options->{ "no_ruler" }, $options->{ "no_consensus" } );
4063     }
4064
4065     close $fh if $fh;
4066 }
4067
4068
4069 sub script_write_blast
4070 {
4071     # Martin A. Hansen, November 2007.
4072
4073     # Write data in blast table format (-m8 and 9).
4074
4075     my ( $in,        # handle to in stream
4076          $out,       # handle to out stream
4077          $options,   # options hash
4078        ) = @_;
4079
4080     # Returns nothing.
4081
4082     my ( $fh, $record, $first );
4083
4084     $fh = write_stream( $options->{ "data_out" }, $options->{ "compress" } ) ;
4085
4086     $first = 1;
4087
4088     while ( $record = get_record( $in ) )
4089     {
4090         if ( $record->{ "REC_TYPE" } eq "BLAST" )
4091         {
4092             if ( $options->{ "comment" } and $first )
4093             {
4094                 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";
4095
4096                 $first = 0;
4097             }
4098
4099             if ( $record->{ "STRAND" } eq "-" ) {
4100                 ( $record->{ "S_BEG" }, $record->{ "S_END" } ) = ( $record->{ "S_END" }, $record->{ "S_BEG" } );
4101             }
4102
4103             print $fh join( "\t",
4104                 $record->{ "Q_ID" },
4105                 $record->{ "S_ID" },
4106                 $record->{ "IDENT" },
4107                 $record->{ "ALIGN_LEN" },
4108                 $record->{ "MISMATCHES" },
4109                 $record->{ "GAPS" },
4110                 $record->{ "Q_BEG" } + 1,
4111                 $record->{ "Q_END" } + 1,
4112                 $record->{ "S_BEG" } + 1,
4113                 $record->{ "S_END" } + 1,
4114                 $record->{ "E_VAL" },
4115                 $record->{ "BIT_SCORE" }
4116             ), "\n";
4117         }
4118
4119         put_record( $record, $out ) if not $options->{ "no_stream" };
4120     }
4121
4122     close $fh;
4123 }
4124
4125
4126 sub script_write_tab
4127 {
4128     # Martin A. Hansen, August 2007.
4129
4130     # Write data as table.
4131
4132     my ( $in,        # handle to in stream
4133          $out,       # handle to out stream
4134          $options,   # options hash
4135        ) = @_;
4136
4137     # Returns nothing.
4138
4139     my ( $fh, $record, $key, @keys, @vals, $ok, %no_keys, $A, $B );
4140
4141     $options->{ "delimit" } ||= "\t";
4142
4143     map { $no_keys{ $_ } = 1 } @{ $options->{ "no_keys" } };
4144
4145     $fh = write_stream( $options->{ "data_out" }, $options->{ "compress" } );
4146
4147     while ( $record = get_record( $in ) )
4148     {
4149         undef @vals;
4150         $ok = 1;
4151
4152         if ( $options->{ "keys" } )
4153         {
4154             map { $ok = 0 if not exists $record->{ $_ } } @{ $options->{ "keys" } };
4155
4156             if ( $ok )
4157             {
4158                 foreach $key ( @{ $options->{ "keys" }  } )
4159                 {
4160                     if ( exists $record->{ $key } )
4161                     {
4162                         push @keys, $key if $options->{ "comment" };
4163                         push @vals, $record->{ $key };
4164                     }
4165                 }
4166              }
4167         }
4168         else
4169         {
4170             foreach $key ( sort { $A = $a; $B = $b; $A =~ s/^V(\d+)$/$1/; $B =~ s/^V(\d+)$/$1/; $A <=> $B } keys %{ $record } )
4171             {
4172                 next if exists $no_keys{ $key };
4173
4174                 push @keys, $key if $options->{ "comment" };
4175                 push @vals, $record->{ $key };
4176             }
4177         }
4178
4179         if ( @keys and $options->{ "comment" } )
4180         {
4181             print $fh "#", join( $options->{ "delimit" }, @keys ), "\n";
4182
4183             delete $options->{ "comment" };
4184         }
4185
4186         print $fh join( $options->{ "delimit" }, @vals ), "\n" if @vals;
4187
4188         put_record( $record, $out ) if not $options->{ "no_stream" };
4189     }
4190
4191     close $fh;
4192 }
4193
4194
4195 sub script_write_bed
4196 {
4197     # Martin A. Hansen, August 2007.
4198
4199     # Write BED format for the UCSC genome browser using records in stream.
4200
4201     my ( $in,        # handle to in stream
4202          $out,       # handle to out stream
4203          $options,   # options hash
4204        ) = @_;
4205
4206     # Returns nothing.
4207
4208     my ( $fh, $record, $new_record );
4209
4210     $fh = write_stream( $options->{ "data_out" }, $options->{ "compress" } );
4211
4212     while ( $record = get_record( $in ) )
4213     {
4214         if ( $record->{ "REC_TYPE" } eq "BED" )                                             # ---- Hits from BED ----
4215         {
4216             Maasha::UCSC::bed_put_entry( $record, $fh, $record->{ "BED_COLS" } );
4217         }
4218         elsif ( $record->{ "REC_TYPE" } eq "PSL" and $record->{ "S_ID" } =~ /^chr/i )       # ---- Hits from BLAT (PSL) ----
4219         {
4220             $new_record->{ "CHR" }     = $record->{ "S_ID" };
4221             $new_record->{ "CHR_BEG" } = $record->{ "S_BEG" };
4222             $new_record->{ "CHR_END" } = $record->{ "S_END" };
4223             $new_record->{ "Q_ID" }    = $record->{ "Q_ID" };
4224             $new_record->{ "SCORE" }   = $record->{ "SCORE" } || 999;
4225             $new_record->{ "STRAND" }  = $record->{ "STRAND" };
4226
4227             Maasha::UCSC::bed_put_entry( $new_record, $fh, 6 );
4228         }
4229         elsif ( $record->{ "REC_TYPE" } eq "PATSCAN" and $record->{ "CHR" } )               # ---- Hits from patscan_seq ----
4230         {
4231             Maasha::UCSC::bed_put_entry( $record, $fh, 6 );
4232         }
4233         elsif ( $record->{ "REC_TYPE" } eq "BLAST" and $record->{ "S_ID" } =~ /^chr/i )     # ---- Hits from BLAST ----
4234         {
4235             $new_record->{ "CHR" }     = $record->{ "S_ID" };
4236             $new_record->{ "CHR_BEG" } = $record->{ "S_BEG" };
4237             $new_record->{ "CHR_END" } = $record->{ "S_END" };
4238             $new_record->{ "Q_ID" }    = $record->{ "Q_ID" };
4239             $new_record->{ "SCORE" }   = $record->{ "SCORE" } || 999; # or use E_VAL somehow
4240             $new_record->{ "STRAND" }  = $record->{ "STRAND" };
4241
4242             Maasha::UCSC::bed_put_entry( $new_record, $fh, 6 );
4243         }
4244         elsif ( $record->{ "REC_TYPE" } eq "VMATCH" and $record->{ "S_ID" } =~ /^chr/i )    # ---- Hits from Vmatch ----
4245         {
4246             $new_record->{ "CHR" }     = $record->{ "S_ID" };
4247             $new_record->{ "CHR_BEG" } = $record->{ "S_BEG" };
4248             $new_record->{ "CHR_END" } = $record->{ "S_END" };
4249             $new_record->{ "Q_ID" }    = $record->{ "Q_ID" };
4250             $new_record->{ "SCORE" }   = $record->{ "SCORE" } || 999; # or use E_VAL somehow
4251             $new_record->{ "STRAND" }  = $record->{ "STRAND" };
4252
4253             Maasha::UCSC::bed_put_entry( $new_record, $fh, 6 );
4254         }
4255         elsif ( $record->{ "REC_TYPE" } eq "SOAP" and $record->{ "S_ID" } =~ /^chr/i )    # ---- Hits from Vmatch ----
4256         {
4257             $new_record->{ "CHR" }     = $record->{ "S_ID" };
4258             $new_record->{ "CHR_BEG" } = $record->{ "S_BEG" };
4259             $new_record->{ "CHR_END" } = $record->{ "S_END" };
4260             $new_record->{ "Q_ID" }    = $record->{ "Q_ID" };
4261             $new_record->{ "SCORE" }   = $record->{ "SCORE" } || 999;
4262             $new_record->{ "STRAND" }  = $record->{ "STRAND" };
4263
4264             Maasha::UCSC::bed_put_entry( $new_record, $fh, 6 );
4265         }
4266         elsif ( $record->{ "CHR" } and defined $record->{ "CHR_BEG" } and $record->{ "CHR_END" } )  # ---- Generic data from tables ----
4267         {
4268             Maasha::UCSC::bed_put_entry( $record, $fh );
4269         }
4270
4271         put_record( $record, $out ) if not $options->{ "no_stream" };
4272     }
4273
4274     close $fh;
4275 }
4276
4277
4278 sub script_write_psl
4279 {
4280     # Martin A. Hansen, August 2007.
4281
4282     # Write PSL output from stream.
4283
4284     my ( $in,        # handle to in stream
4285          $out,       # handle to out stream
4286          $options,   # options hash
4287        ) = @_;
4288
4289     # Returns nothing.
4290
4291     my ( $fh, $record, @output, $first );
4292
4293     $first = 1;
4294
4295     $fh = write_stream( $options->{ "data_out" }, $options->{ "compress" } );
4296
4297     while ( $record = get_record( $in ) )
4298     {
4299         put_record( $record, $out ) if not $options->{ "no_stream" };
4300
4301         if ( $record->{ "REC_TYPE" } and $record->{ "REC_TYPE" } eq "PSL" )
4302         {
4303             Maasha::UCSC::psl_put_header( $fh ) if $first;
4304             Maasha::UCSC::psl_put_entry( $record, $fh );
4305             $first = 0;
4306         }
4307     }
4308
4309     close $fh;
4310 }
4311
4312
4313 sub script_write_fixedstep
4314 {
4315     # Martin A. Hansen, Juli 2008.
4316
4317     # Write fixedStep entries from recrods in the stream.
4318
4319     my ( $in,        # handle to in stream
4320          $out,       # handle to out stream
4321          $options,   # options hash
4322        ) = @_;
4323
4324     # Returns nothing.
4325
4326     my ( $fh, $record, $vals );
4327
4328     $fh = write_stream( $options->{ "data_out" }, $options->{ "compress" } );
4329
4330     while ( $record = get_record( $in ) )
4331     {
4332         if ( $record->{ "CHR" } and $record->{ "CHR_BEG" } and $record->{ "STEP" } and $record->{ "VALS" } )
4333         {
4334             print $fh "fixedStep chrom=$record->{ 'CHR' } start=$record->{ 'CHR_BEG' } step=$record->{ 'STEP' }\n";
4335
4336             $vals = $record->{ 'VALS' };
4337
4338             $vals =~ tr/;/\n/;
4339
4340             print $fh "$vals\n";
4341         }
4342
4343         put_record( $record, $out ) if not $options->{ "no_stream" };
4344     }
4345
4346     close $fh;
4347 }
4348
4349
4350 sub script_write_2bit
4351 {
4352     # Martin A. Hansen, March 2008.
4353
4354     # Write sequence entries from stream in 2bit format.
4355
4356     my ( $in,        # handle to in stream
4357          $out,       # handle to out stream
4358          $options,   # options hash
4359        ) = @_;
4360
4361     # Returns nothing.
4362
4363     my ( $record, $mask, $tmp_file, $fh_tmp, $fh_in, $fh_out, $entry );
4364
4365     $mask = 1 if not $options->{ "no_mask" };
4366
4367     $tmp_file = "$BP_TMP/write_2bit.fna";
4368     $fh_tmp   = Maasha::Common::write_open( $tmp_file );
4369
4370     $fh_out = write_stream( $options->{ "data_out" } );
4371
4372     while ( $record = get_record( $in ) )
4373     {
4374         if ( $entry = record2fasta( $record ) ) {
4375             Maasha::Fasta::put_entry( $entry, $fh_tmp );
4376         }
4377
4378         put_record( $record, $out ) if not $options->{ "no_stream" };
4379     }
4380
4381     close $fh_tmp;
4382
4383     $fh_in = Maasha::Common::read_open( $tmp_file );
4384
4385     Maasha::TwoBit::fasta2twobit( $fh_in, $fh_out, $mask );
4386
4387     close $fh_in;
4388     close $fh_out;
4389
4390     unlink $tmp_file;
4391 }
4392
4393
4394 sub script_write_solid
4395 {
4396     # Martin A. Hansen, April 2008.
4397
4398     # Write di-base encoded Solid sequence from entries in stream.
4399
4400     my ( $in,        # handle to in stream
4401          $out,       # handle to out stream
4402          $options,   # options hash
4403        ) = @_;
4404
4405     # Returns nothing.
4406
4407     my ( $record, $fh, $entry );
4408
4409     $fh = write_stream( $options->{ "data_out" }, $options->{ "compress" } );
4410
4411     while ( $record = get_record( $in ) )
4412     {
4413         if ( $entry = record2fasta( $record ) )
4414         {
4415             $entry->[ SEQ ] = Maasha::Solid::seq2color_space( uc $entry->[ SEQ ] );
4416
4417             Maasha::Fasta::put_entry( $entry, $fh, $options->{ "wrap" } );
4418         }
4419
4420         put_record( $record, $out ) if not $options->{ "no_stream" };
4421     }
4422
4423     close $fh;
4424 }
4425
4426
4427 sub script_plot_seqlogo
4428 {
4429     # Martin A. Hansen, August 2007.
4430
4431     # Calculates and writes a sequence logo for alignments.
4432
4433     my ( $in,        # handle to in stream
4434          $out,       # handle to out stream
4435          $options,   # options hash
4436        ) = @_;
4437
4438     # Returns nothing.
4439
4440     my ( $record, @entries, $logo, $fh );
4441
4442     while ( $record = get_record( $in ) )
4443     {
4444         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } ) {
4445             push @entries, [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ];
4446         }
4447
4448         put_record( $record, $out ) if not $options->{ "no_stream" };
4449     }
4450
4451     $logo = Maasha::Plot::seq_logo( \@entries );
4452
4453     $fh = write_stream( $options->{ "data_out" } );
4454
4455     print $fh $logo;
4456
4457     close $fh;
4458 }
4459
4460
4461 sub script_plot_phastcons_profiles
4462 {
4463     # Martin A. Hansen, January 2008.
4464
4465     # Plots PhastCons profiles.
4466
4467     my ( $in,        # handle to in stream
4468          $out,       # handle to out stream
4469          $options,   # options hash
4470        ) = @_;
4471
4472     # Returns nothing.
4473
4474     my ( $phastcons_file, $phastcons_index, $index, $fh_phastcons, $record, $scores, $AoA, $plot, $fh );
4475
4476     $options->{ "title" } ||= "PhastCons Profiles";
4477
4478     $phastcons_file  = Maasha::Config::genome_phastcons( $options->{ "genome" } );
4479     $phastcons_index = Maasha::Config::genome_phastcons_index( $options->{ "genome" } );
4480
4481     $index           = Maasha::UCSC::fixedstep_index_retrieve( $phastcons_index );
4482     $fh_phastcons    = Maasha::Common::read_open( $phastcons_file );
4483
4484     while ( $record = get_record( $in ) )
4485     {
4486         if ( $record->{ "CHR" } and $record->{ "CHR_BEG" } and $record->{ "CHR_END" } )
4487         {
4488             $scores = Maasha::UCSC::fixedstep_index_lookup( $index, $fh_phastcons, $record->{ "CHR" },
4489                                                                                    $record->{ "CHR_BEG" },
4490                                                                                    $record->{ "CHR_END" },
4491                                                                                    $options->{ "flank" } );
4492
4493             push @{ $AoA }, [ @{ $scores } ];
4494         }
4495
4496         put_record( $record, $out ) if not $options->{ "no_stream" };
4497     }
4498
4499     Maasha::UCSC::phastcons_normalize( $AoA );
4500
4501     $AoA = [ [ Maasha::UCSC::phastcons_mean( $AoA ) ] ] if $options->{ "mean" };
4502     $AoA = [ [ Maasha::UCSC::phastcons_median( $AoA ) ] ] if $options->{ "median" };
4503
4504     $AoA = Maasha::Matrix::matrix_flip( $AoA );
4505
4506     $plot = Maasha::Plot::lineplot_simple( $AoA, $options, $BP_TMP );
4507
4508     $fh = write_stream( $options->{ "data_out" } );
4509
4510     print $fh "$_\n" foreach @{ $plot };
4511
4512     close $fh;
4513 }
4514
4515
4516 sub script_analyze_bed
4517 {
4518     # Martin A. Hansen, March 2008.
4519
4520     # Analyze BED entries in stream.
4521
4522     my ( $in,        # handle to in stream
4523          $out,       # handle to out stream
4524          $options,   # options hash
4525        ) = @_;
4526
4527     # Returns nothing.
4528
4529     my ( $record );
4530
4531     while ( $record = get_record( $in ) )
4532     {
4533         $record = Maasha::UCSC::bed_analyze( $record ) if $record->{ "REC_TYPE" } eq "BED";
4534
4535         put_record( $record, $out );
4536     }
4537 }
4538
4539
4540 sub script_analyze_vals
4541 {
4542     # Martin A. Hansen, August 2007.
4543
4544     # Analyze values for given keys in stream.
4545
4546     my ( $in,        # handle to in stream
4547          $out,       # handle to out stream
4548          $options,   # options hash
4549        ) = @_;
4550
4551     # Returns nothing.
4552
4553     my ( $record, $key, @keys, %key_hash, $analysis, $len );
4554
4555     map { $key_hash{ $_ } = 1 } @{ $options->{ "keys" } };
4556
4557     while ( $record = get_record( $in ) )
4558     {
4559         foreach $key ( keys %{ $record } )
4560         {
4561             next if $options->{ "keys" } and not exists $key_hash{ $key };
4562
4563             $analysis->{ $key }->{ "COUNT" }++;
4564
4565             if ( Maasha::Calc::is_a_number( $record->{ $key } ) )
4566             {
4567                 $analysis->{ $key }->{ "TYPE" } = "num";
4568                 $analysis->{ $key }->{ "SUM" } += $record->{ $key };
4569                 $analysis->{ $key }->{ "MAX" } = $record->{ $key } if $record->{ $key } > $analysis->{ $key }->{ "MAX" } or not $analysis->{ $key }->{ "MAX" };
4570                 $analysis->{ $key }->{ "MIN" } = $record->{ $key } if $record->{ $key } < $analysis->{ $key }->{ "MIN" } or not $analysis->{ $key }->{ "MIN" };
4571             }
4572             else
4573             {
4574                 $len = length $record->{ $key };
4575
4576                 $analysis->{ $key }->{ "TYPE" } = "alph";
4577                 $analysis->{ $key }->{ "SUM" } += $len;
4578                 $analysis->{ $key }->{ "MAX" } = $len if $len > $analysis->{ $key }->{ "MAX" } or not $analysis->{ $key }->{ "MAX" };
4579                 $analysis->{ $key }->{ "MIN" } = $len if $len < $analysis->{ $key }->{ "MIM" } or not $analysis->{ $key }->{ "MIN" };
4580             }
4581         }
4582
4583         put_record( $record, $out ) if not $options->{ "no_stream" };
4584     }
4585
4586     foreach $key ( keys %{ $analysis } )
4587     {
4588         $analysis->{ $key }->{ "MEAN" } = sprintf "%.2f", $analysis->{ $key }->{ "SUM" } / $analysis->{ $key }->{ "COUNT" };
4589         $analysis->{ $key }->{ "SUM" }  = sprintf "%.2f", $analysis->{ $key }->{ "SUM" };
4590     }
4591
4592     my ( $keys, $types, $counts, $mins, $maxs, $sums, $means );
4593
4594     $keys   = "KEY  ";
4595     $types  = "TYPE ";
4596     $counts = "COUNT";
4597     $mins   = "MIN  ";
4598     $maxs   = "MAX  ";
4599     $sums   = "SUM  ";
4600     $means  = "MEAN ";
4601
4602     if ( $options->{ "keys" } ) {
4603         @keys = @{ $options->{ "keys" } };
4604     } else {
4605         @keys = keys %{ $analysis };
4606     }
4607
4608     foreach $key ( @keys )
4609     {
4610         $keys   .= sprintf "% 15s", $key;
4611         $types  .= sprintf "% 15s", $analysis->{ $key }->{ "TYPE" };
4612         $counts .= sprintf "% 15s", $analysis->{ $key }->{ "COUNT" };
4613         $mins   .= sprintf "% 15s", $analysis->{ $key }->{ "MIN" };
4614         $maxs   .= sprintf "% 15s", $analysis->{ $key }->{ "MAX" };
4615         $sums   .= sprintf "% 15s", $analysis->{ $key }->{ "SUM" };
4616         $means  .= sprintf "% 15s", $analysis->{ $key }->{ "MEAN" };
4617     }
4618
4619     print $out "$keys\n";
4620     print $out "$types\n";
4621     print $out "$counts\n";
4622     print $out "$mins\n";
4623     print $out "$maxs\n";
4624     print $out "$sums\n";
4625     print $out "$means\n";
4626 }
4627
4628
4629 sub script_head_records
4630 {
4631     # Martin A. Hansen, August 2007.
4632
4633     # Display the first sequences in stream.
4634
4635     my ( $in,        # handle to in stream
4636          $out,       # handle to out stream
4637          $options,   # options hash
4638        ) = @_;
4639
4640     # Returns nothing.
4641
4642     my ( $record, $count );
4643
4644     $options->{ "num" } ||= 10;
4645
4646     $count = 0;
4647
4648     while ( $record = get_record( $in ) )
4649     {
4650         $count++;
4651
4652         put_record( $record, $out );
4653
4654         last if $count == $options->{ "num" };
4655     }
4656 }
4657
4658
4659 sub script_remove_keys
4660 {
4661     # Martin A. Hansen, August 2007.
4662
4663     # Remove keys from stream.
4664
4665     my ( $in,        # handle to in stream
4666          $out,       # handle to out stream
4667          $options,   # options hash
4668        ) = @_;
4669
4670     # Returns nothing.
4671
4672     my ( $record, $new_record );
4673
4674     while ( $record = get_record( $in ) )
4675     {
4676         if ( $options->{ "keys" } )
4677         {
4678             map { delete $record->{ $_ } } @{ $options->{ "keys" } };
4679         }
4680         elsif ( $options->{ "save_keys" } )
4681         {
4682             map { $new_record->{ $_ } = $record->{ $_ } if exists $record->{ $_ } } @{ $options->{ "save_keys" } };
4683
4684             $record = $new_record;
4685         }
4686
4687         put_record( $record, $out ) if keys %{ $record };
4688     }
4689 }
4690
4691
4692 sub script_remove_adaptor
4693 {
4694     # Martin A. Hansen, August 2008.
4695
4696     # Find and remove adaptor from sequences in the stream.
4697
4698     my ( $in,        # handle to in stream
4699          $out,       # handle to out stream
4700          $options,   # options hash
4701        ) = @_;
4702
4703     # Returns nothing.
4704
4705     my ( $record, $adaptor, $seq, $adaptor_len, $seq_len, $offset, $max_match, $max_mismatch, $pos );
4706
4707     $options->{ "remove" } ||= "after";
4708
4709     $max_mismatch = $options->{ "mismatches" } || 0;
4710     $offset       = $options->{ "offset" };
4711
4712     if ( not defined $offset ) {
4713         $offset = 0;
4714     } else {
4715         $offset--;
4716     }
4717
4718     $adaptor      = uc $options->{ "adaptor" };
4719     $adaptor_len  = length $adaptor;
4720
4721     while ( $record = get_record( $in ) )
4722     {
4723         if ( $record->{ "SEQ" } )
4724         {
4725             $seq     = uc $record->{ "SEQ" };
4726             $seq_len = length $seq;
4727
4728             $pos = Maasha::Common::index_m( $seq, $adaptor, $seq_len, $adaptor_len, $offset, $max_mismatch );
4729
4730             $record->{ "ADAPTOR_POS" } = $pos;
4731
4732             if ( $pos >= 0 and $options->{ "remove" } ne "skip" )
4733             {
4734                 if ( $options->{ "remove" } eq "after" )
4735                 {
4736                     $record->{ "SEQ" }     = substr $record->{ "SEQ" }, 0, $pos;
4737                     $record->{ "SEQ_LEN" } = $pos;
4738                 }
4739                 else
4740                 {
4741                     $record->{ "SEQ" }     = substr $record->{ "SEQ" }, $pos + $adaptor_len;
4742                     $record->{ "SEQ_LEN" } = length $record->{ "SEQ" };
4743                 }
4744             }
4745
4746             put_record( $record, $out );
4747         }
4748         else
4749         {
4750             put_record( $record, $out );
4751         }
4752     }
4753 }
4754
4755
4756 sub script_rename_keys
4757 {
4758     # Martin A. Hansen, August 2007.
4759
4760     # Rename keys in stream.
4761
4762     my ( $in,        # handle to in stream
4763          $out,       # handle to out stream
4764          $options,   # options hash
4765        ) = @_;
4766
4767     # Returns nothing.
4768
4769     my ( $record );
4770
4771     while ( $record = get_record( $in ) )
4772     {
4773         if ( exists $record->{ $options->{ "keys" }->[ 0 ] } )
4774         {
4775             $record->{ $options->{ "keys" }->[ 1 ] } = $record->{ $options->{ "keys" }->[ 0 ] };
4776
4777             delete $record->{ $options->{ "keys" }->[ 0 ] };
4778         }
4779
4780         put_record( $record, $out );
4781     }
4782 }
4783
4784
4785 sub script_uniq_vals
4786 {
4787     # Martin A. Hansen, August 2007.
4788
4789     # Find unique values in stream.
4790
4791     my ( $in,        # handle to in stream
4792          $out,       # handle to out stream
4793          $options,   # options hash
4794        ) = @_;
4795
4796     # Returns nothing.
4797
4798     my ( %hash, $record );
4799
4800     while ( $record = get_record( $in ) )
4801     {
4802         if ( $record->{ $options->{ "key" } } )
4803         {
4804             if ( not $hash{ $record->{ $options->{ "key" } } } and not $options->{ "invert" } )
4805             {
4806                 put_record( $record, $out );
4807
4808                 $hash{ $record->{ $options->{ "key" } } } = 1;
4809             }
4810             elsif ( $hash{ $record->{ $options->{ "key" } } } and $options->{ "invert" } )
4811             {
4812                 put_record( $record, $out );
4813             }
4814             else
4815             {
4816                 $hash{ $record->{ $options->{ "key" } } } = 1;
4817             }
4818         }
4819         else
4820         {
4821             put_record( $record, $out );
4822         }
4823     }
4824 }
4825
4826
4827 sub script_merge_vals
4828 {
4829     # Martin A. Hansen, August 2007.
4830
4831     # Rename keys in stream.
4832
4833     my ( $in,        # handle to in stream
4834          $out,       # handle to out stream
4835          $options,   # options hash
4836        ) = @_;
4837
4838     # Returns nothing.
4839
4840     my ( $record, @join, $i );
4841
4842     $options->{ "delimit" } ||= '_';
4843
4844     while ( $record = get_record( $in ) )
4845     {
4846         if ( exists $record->{ $options->{ "keys" }->[ 0 ] } )
4847         {
4848             @join = $record->{ $options->{ "keys" }->[ 0 ] };
4849
4850             for ( $i = 1; $i < @{ $options->{ "keys" } }; $i++ ) {
4851                 push @join, $record->{ $options->{ "keys" }->[ $i ] } if exists $record->{ $options->{ "keys" }->[ $i ] };
4852             }
4853
4854             $record->{ $options->{ "keys" }->[ 0 ] } = join $options->{ "delimit" }, @join;
4855         }
4856
4857         put_record( $record, $out );
4858     }
4859 }
4860
4861
4862 sub script_merge_records
4863 {
4864     # Martin A. Hansen, July 2008.
4865
4866     # Merges records in the stream based on identical values of two given keys.
4867
4868     my ( $in,        # handle to in stream
4869          $out,       # handle to out stream
4870          $options,   # options hash
4871        ) = @_;
4872
4873     # Returns nothing.
4874
4875     my ( $merge, $record, $file1, $file2, $fh1, $fh2, $key1, $key2, @keys1, @keys2, @vals1, @vals2,
4876          $num1, $num2, $num, $cmp, $i );
4877
4878     $merge = $options->{ "merge" } || "AandB";
4879
4880     $file1 = "$BP_TMP/merge_records1.tmp";
4881     $file2 = "$BP_TMP/merge_records2.tmp";
4882
4883     $fh1   = Maasha::Common::write_open( $file1 );
4884     $fh2   = Maasha::Common::write_open( $file2 );
4885
4886     $key1  = $options->{ "keys" }->[ 0 ];
4887     $key2  = $options->{ "keys" }->[ 1 ];
4888
4889     $num   = $key2 =~ s/n$//;
4890     $num1  = 0;
4891     $num2  = 0;
4892
4893     while ( $record = get_record( $in ) )
4894     {
4895         if ( exists $record->{ $key1 } )
4896         {
4897             @keys1 = $key1;
4898             @vals1 = $record->{ $key1 };
4899
4900             delete $record->{ $key1 };
4901
4902             map { push @keys1, $_; push @vals1, $record->{ $_ } } keys %{ $record };
4903
4904             print $fh1 join( "\t", @vals1 ), "\n";
4905
4906             $num1++;
4907         }
4908         elsif ( exists $record->{ $key2 } )
4909         {
4910             @keys2 = $key2;
4911             @vals2 = $record->{ $key2 };
4912
4913             delete $record->{ $key2 };
4914
4915             map { push @keys2, $_; push @vals2, $record->{ $_ } } keys %{ $record };
4916
4917             print $fh2 join( "\t", @vals2 ), "\n";
4918
4919             $num2++;
4920         }
4921     }
4922
4923     close $fh1;
4924     close $fh2;
4925
4926     if ( $num )
4927     {
4928         Maasha::Common::run( "sort", "-k 1,1n $file1 > $file1.sort" ) and rename "$file1.sort", $file1;
4929         Maasha::Common::run( "sort", "-k 1,1n $file2 > $file2.sort" ) and rename "$file2.sort", $file2;
4930     }
4931     else
4932     {
4933         Maasha::Common::run( "sort", "-k 1,1 $file1 > $file1.sort" ) and rename "$file1.sort", $file1;
4934         Maasha::Common::run( "sort", "-k 1,1 $file2 > $file2.sort" ) and rename "$file2.sort", $file2;
4935     }
4936
4937     $fh1 = Maasha::Common::read_open( $file1 );
4938     $fh2 = Maasha::Common::read_open( $file2 );
4939
4940     @vals1 = Maasha::Common::get_fields( $fh1 );
4941     @vals2 = Maasha::Common::get_fields( $fh2 );
4942
4943     while ( $num1 > 0 and $num2 > 0 )
4944     {
4945         undef $record;
4946
4947         if ( $num ) {
4948             $cmp = $vals1[ 0 ] <=> $vals2[ 0 ];
4949         } else {
4950             $cmp = $vals1[ 0 ] cmp $vals2[ 0 ];
4951         }
4952
4953         if ( $cmp < 0 )
4954         {
4955             if ( $merge =~ /^(AorB|AnotB)$/ )
4956             {
4957                 for ( $i = 0; $i < @keys1; $i++ ) {
4958                     $record->{ $keys1[ $i ] } = $vals1[ $i ];
4959                 }
4960
4961                 put_record( $record, $out );
4962             }
4963
4964             @vals1 = Maasha::Common::get_fields( $fh1 );
4965             $num1--;
4966         }
4967         elsif ( $cmp > 0 )
4968         {
4969             if ( $merge =~ /^(BorA|BnotA)$/ )
4970             {
4971                 for ( $i = 0; $i < @keys2; $i++ ) {
4972                     $record->{ $keys2[ $i ] } = $vals2[ $i ];
4973                 }
4974
4975                 put_record( $record, $out );
4976             }
4977
4978             @vals2 = Maasha::Common::get_fields( $fh2 );
4979             $num2--;
4980         }
4981         else
4982         {
4983             if ( $merge =~ /^(AandB|AorB|BorA)$/ )
4984             {
4985                 for ( $i = 0; $i < @keys1; $i++ ) {
4986                     $record->{ $keys1[ $i ] } = $vals1[ $i ];
4987                 }
4988
4989                 for ( $i = 1; $i < @keys2; $i++ ) {
4990                     $record->{ $keys2[ $i ] } = $vals2[ $i ];
4991                 }
4992
4993                 put_record( $record, $out );
4994             }
4995
4996             @vals1 = Maasha::Common::get_fields( $fh1 );
4997             @vals2 = Maasha::Common::get_fields( $fh2 );
4998             $num1--;
4999             $num2--;
5000         }
5001     }
5002
5003     close $fh1;
5004     close $fh2;
5005
5006     unlink $file1;
5007     unlink $file2;
5008
5009     if ( $num1 > 0 and $merge =~ /^(AorB|AnotB)$/ )
5010     {
5011         undef $record;
5012
5013         for ( $i = 0; $i < @keys1; $i++ ) {
5014             $record->{ $keys1[ $i ] } = $vals1[ $i ];
5015         }
5016
5017         put_record( $record, $out );
5018     }
5019
5020     if ( $num2 > 0 and $merge =~ /^(BorA|BnotA)$/ )
5021     {
5022         undef $record;
5023
5024         for ( $i = 0; $i < @keys2; $i++ ) {
5025             $record->{ $keys2[ $i ] } = $vals2[ $i ];
5026         }
5027
5028         put_record( $record, $out );
5029     }
5030 }
5031
5032
5033 sub script_grab
5034 {
5035     # Martin A. Hansen, August 2007.
5036
5037     # Grab for records in stream.
5038
5039     my ( $in,        # handle to in stream
5040          $out,       # handle to out stream
5041          $options,   # options hash
5042        ) = @_;
5043
5044     # Returns nothing.
5045
5046     my ( $patterns, $pattern, $record, $key, $pos, $op, $val, %lookup_hash );
5047
5048     if ( $options->{ "patterns" } )
5049     {
5050         $patterns = [ split ",", $options->{ "patterns" } ];
5051     }
5052     elsif ( -f $options->{ "patterns_in" } )
5053     {
5054         $patterns = Maasha::Patscan::read_patterns( $options->{ "patterns_in" } );
5055     }
5056     elsif ( -f $options->{ "exact_in" } )
5057     {
5058         $patterns = Maasha::Patscan::read_patterns( $options->{ "exact_in" } );
5059
5060         map { $lookup_hash{ $_ } = 1 } @{ $patterns };
5061
5062         undef $patterns;
5063     }
5064
5065     if ( $options->{ "eval" } )
5066     {
5067         if ( $options->{ "eval" } =~ /^([^><=! ]+)\s*(>=|<=|>|<|=|!=|eq|ne)\s*(.+)$/ )
5068         {
5069             $key = $1;
5070             $op  = $2;
5071             $val = $3;
5072         }
5073     }
5074
5075     while ( $record = get_record( $in ) )
5076     {
5077         $pos = -1;
5078
5079         if ( %lookup_hash )
5080         {
5081             if ( $options->{ "keys" } )
5082             {
5083                 foreach $key ( @{ $options->{ "keys" } } )
5084                 {
5085                     if ( exists $lookup_hash{ $record->{ $key } } )
5086                     {
5087                         $pos = 1;
5088                         goto FOUND;
5089                     }
5090                 }
5091             }
5092             else
5093             {
5094                 foreach $key ( keys %{ $record } )
5095                 {
5096                     if ( not $options->{ "vals_only" } )
5097                     {
5098                         if ( exists $lookup_hash{ $key } )
5099                         {
5100                             $pos = 1;
5101                             goto FOUND;
5102                         }
5103                     }
5104
5105                     if ( not $options->{ "keys_only" } )
5106                     {
5107                         if ( exists $lookup_hash{ $record->{ $key } } )
5108                         {
5109                             $pos = 1;
5110                             goto FOUND;
5111                         }
5112                     }
5113                 }
5114             }
5115         }
5116         elsif ( $patterns )
5117         {
5118             foreach $pattern ( @{ $patterns } )
5119             {
5120                 if ( $options->{ "keys" } )
5121                 {
5122                     foreach $key ( @{ $options->{ "keys" } } )
5123                     {
5124                         $pos = index $record->{ $key }, $pattern;
5125
5126                         goto FOUND if $pos >= 0;
5127                     }
5128                 }
5129                 else
5130                 {
5131                     foreach $key ( keys %{ $record } )
5132                     {
5133                         if ( not $options->{ "vals_only" } )
5134                         {
5135                             $pos = index $key, $pattern;
5136
5137                             goto FOUND if $pos >= 0;
5138                         }
5139
5140                         if ( not $options->{ "keys_only" } )
5141                         {
5142                             $pos = index $record->{ $key }, $pattern;
5143
5144                             goto FOUND if $pos >= 0;
5145                         }
5146                     }
5147                 }
5148             }
5149         }
5150         elsif ( $options->{ "regex" } )
5151         {
5152             if ( $options->{ "keys" } )
5153             {
5154                 foreach $key ( @{ $options->{ "keys" } } )
5155                 {
5156                     if ( $options->{ "case_insensitive" } ) {
5157                         $pos = 1 if $record->{ $key } =~ /$options->{'regex'}/i;
5158                     } else {
5159                         $pos = 1 if $record->{ $key } =~ /$options->{'regex'}/;
5160                     }
5161
5162                     goto FOUND if $pos >= 0;
5163                 }
5164             }
5165             else
5166             {
5167                 foreach $key ( keys %{ $record } )
5168                 {
5169                     if ( not $options->{ "vals_only" } )
5170                     {
5171                         if ( $options->{ "case_insensitive" } ) {
5172                             $pos = 1 if $key =~ /$options->{'regex'}/i;
5173                         } else {
5174                             $pos = 1 if $key =~ /$options->{'regex'}/;
5175                         }
5176
5177                         goto FOUND if $pos >= 0;
5178                     }
5179
5180                     if ( not $options->{ "keys_only" } )
5181                     {
5182                         if ( $options->{ "case_insensitive" } ) {
5183                             $pos = 1 if $record->{ $key } =~ /$options->{'regex'}/i;
5184                         } else {
5185                             $pos = 1 if $record->{ $key } =~ /$options->{'regex'}/;
5186                         }
5187
5188                         goto FOUND if $pos >= 0;
5189                     }
5190                 }
5191             }
5192         }
5193         elsif ( $options->{ "eval" } )
5194         {
5195             if ( defined $record->{ $key } )
5196             {
5197                 if ( $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 "=" and $record->{ $key } == $val ) {
5206                     $pos = 1 and goto FOUND;
5207                 } elsif ( $op eq "!=" and $record->{ $key } != $val ) {
5208                     $pos = 1 and goto FOUND;
5209                 } elsif ( $op eq "eq" and $record->{ $key } eq $val ) {
5210                     $pos = 1 and goto FOUND;
5211                 } elsif ( $op eq "ne" and $record->{ $key } ne $val ) {
5212                     $pos = 1 and goto FOUND;
5213                 }
5214             }
5215         }
5216
5217         FOUND:
5218
5219         if ( $pos >= 0 and not $options->{ "invert" } ) {
5220             put_record( $record, $out );
5221         } elsif ( $pos < 0 and $options->{ "invert" } ) {
5222             put_record( $record, $out );
5223         }
5224     }
5225 }
5226
5227
5228 sub script_compute
5229 {
5230     # Martin A. Hansen, August 2007.
5231
5232     # Evaluate extression for records in stream.
5233
5234     my ( $in,        # handle to in stream
5235          $out,       # handle to out stream
5236          $options,   # options hash
5237        ) = @_;
5238
5239     # Returns nothing.
5240
5241     my ( $record, $eval_key, $eval_val, $check, @keys );
5242
5243     while ( $record = get_record( $in ) )
5244     {
5245         if ( $options->{ "eval" } )
5246         {
5247             if ( $options->{ "eval" } =~ /^(.+)\s*=\s*(.+)$/ )
5248             {
5249                 $eval_key = $1;
5250                 $eval_val = $2;
5251             }
5252
5253             if ( not $check )
5254             {
5255                 @keys = split /\W+/, $eval_val;
5256                 @keys = grep { ! /^\d+$/ } @keys;
5257
5258                 $check = 1;
5259             }
5260
5261             map { $eval_val =~ s/$_/$record->{ $_ }/g } @keys;
5262
5263             $record->{ $eval_key } = eval "$eval_val" or Maasha::Common::error( "eval failed -> $@" );
5264         }
5265
5266         put_record( $record, $out );
5267     }
5268 }
5269
5270
5271 sub script_flip_tab
5272 {
5273     # Martin A. Hansen, June 2008.
5274
5275     # Flip a table.
5276
5277     my ( $in,        # handle to in stream
5278          $out,       # handle to out stream
5279          $options,   # options hash
5280        ) = @_;
5281
5282     # Returns nothing.
5283
5284     my ( $record, $key, $A, $B, @rows, @matrix, $row, $i );
5285
5286     while ( $record = get_record( $in ) )
5287     {
5288         undef @rows;
5289
5290         foreach $key ( sort { $A = $a; $B = $b; $A =~ s/^V(\d+)$/$1/; $B =~ s/^V(\d+)$/$1/; $A <=> $B } keys %{ $record } )
5291         {
5292             push @rows, $record->{ $key };
5293
5294         }
5295
5296         push @matrix, [ @rows ];
5297     }
5298
5299     undef $record;
5300
5301     @matrix = Maasha::Matrix::matrix_flip( \@matrix );
5302
5303     foreach $row ( @matrix )
5304     {
5305         for ( $i = 0; $i < @{ $row }; $i++ ) {
5306             $record->{ "V$i" } = $row->[ $i ];
5307         }
5308
5309         put_record( $record, $out );
5310     }
5311 }
5312
5313
5314 sub script_add_ident
5315 {
5316     # Martin A. Hansen, May 2008.
5317
5318     # Add a unique identifier to each record in stream.
5319
5320     my ( $in,        # handle to in stream
5321          $out,       # handle to out stream
5322          $options,   # options hash
5323        ) = @_;
5324
5325     # Returns nothing.
5326
5327     my ( $record, $key, $prefix, $i );
5328
5329     $key    = $options->{ "key" }    || "ID";
5330     $prefix = $options->{ "prefix" } || "ID";
5331
5332     $i = 0;
5333
5334     while ( $record = get_record( $in ) )
5335     {
5336         $record->{ $key } = sprintf( "$prefix%08d", $i );
5337
5338         put_record( $record, $out );
5339
5340         $i++;
5341     }
5342 }
5343
5344
5345 sub script_count_records
5346 {
5347     # Martin A. Hansen, August 2007.
5348
5349     # Count records in stream.
5350
5351     my ( $in,        # handle to in stream
5352          $out,       # handle to out stream
5353          $options,   # options hash
5354        ) = @_;
5355
5356     # Returns nothing.
5357
5358     my ( $record, $count, $result, $fh, $line );
5359
5360     $count = 0;
5361
5362     if ( $options->{ "no_stream" } )
5363     {
5364         while ( $line = <$in> )
5365         {
5366             chomp $line;
5367
5368             $count++ if $line eq "---";
5369         }
5370     }
5371     else
5372     {
5373         while ( $record = get_record( $in ) )
5374         {
5375             put_record( $record, $out );
5376
5377             $count++;
5378         }
5379     }
5380
5381     $result = { "RECORDS_COUNT" => $count };
5382
5383     $fh = write_stream( $options->{ "data_out" } );
5384
5385     put_record( $result, $fh );
5386
5387     close $fh;
5388 }
5389
5390
5391 sub script_random_records
5392 {
5393     # Martin A. Hansen, August 2007.
5394
5395     # Pick a number or random records from stream.
5396
5397     my ( $in,        # handle to in stream
5398          $out,       # handle to out stream
5399          $options,   # options hash
5400        ) = @_;
5401
5402     # Returns nothing.
5403
5404     my ( $record, $tmp_file, $fh_out, $fh_in, $count, $i, %rand_hash, $rand, $max );
5405
5406     $options->{ "num" } ||= 10;
5407
5408     $tmp_file = "$BP_TMP/random_records.tmp";
5409
5410     $fh_out = Maasha::Common::write_open( $tmp_file );
5411
5412     $count = 0;
5413
5414     while ( $record = get_record( $in ) )
5415     {
5416         put_record( $record, $fh_out );
5417
5418         $count++;
5419     }
5420
5421     close $fh_out;
5422
5423     $max = 0;
5424     $i   = 0;
5425
5426     Maasha::Common::error( qq(Requested random records > records in stream) ) if $options->{ "num" } > $count;
5427
5428     while ( $i < $options->{ "num" } )
5429     {
5430         $rand = int( rand( $count ) );
5431
5432         if ( not exists $rand_hash{ $rand } )
5433         {
5434             $rand_hash{ $rand } = 1;
5435
5436             $max = $rand if $rand > $max;
5437
5438             $i++;
5439         }
5440     }
5441
5442     $fh_in = Maasha::Common::read_open( $tmp_file );
5443
5444     $count = 0;
5445
5446     while ( $record = get_record( $fh_in ) )
5447     {
5448         put_record( $record, $out ) if exists $rand_hash{ $count };
5449
5450         last if $count == $max;
5451
5452         $count++;
5453     }
5454
5455     close $fh_in;
5456
5457     unlink $tmp_file;
5458 }
5459
5460
5461 sub script_sort_records
5462 {
5463     # Martin A. Hansen, August 2007.
5464
5465     # Sort to sort records according to keys.
5466
5467     my ( $in,        # handle to in stream
5468          $out,       # handle to out stream
5469          $options,   # options hash
5470        ) = @_;
5471
5472     # Returns nothing.
5473
5474     my ( @keys, $key, @sort_cmd, $sort_str, $sort_sub, @records, $record, $i );
5475
5476     foreach $key ( @{ $options->{ "keys" } } )
5477     {
5478         if ( $key =~ s/n$// ) {
5479             push @sort_cmd, qq(\$a->{ "$key" } <=> \$b->{ "$key" });
5480         } else {
5481             push @sort_cmd, qq(\$a->{ "$key" } cmp \$b->{ "$key" });
5482         }
5483     }
5484
5485     $sort_str = join " or ", @sort_cmd;
5486     $sort_sub = eval "sub { $sort_str }";   # NB security issue!
5487
5488     while ( $record = get_record( $in ) ) {
5489         push @records, $record;
5490     }
5491
5492     @records = sort $sort_sub @records;
5493
5494     if ( $options->{ "reverse" } )
5495     {
5496         for ( $i = scalar @records - 1; $i >= 0; $i-- ) {
5497             put_record( $records[ $i ], $out );
5498         }
5499     }
5500     else
5501     {
5502         for ( $i = 0; $i < scalar @records; $i++ ) {
5503             put_record( $records[ $i ], $out );
5504         }
5505     }
5506 }
5507
5508
5509 sub script_count_vals
5510 {
5511     # Martin A. Hansen, August 2007.
5512
5513     # Count records in stream.
5514
5515     my ( $in,        # handle to in stream
5516          $out,       # handle to out stream
5517          $options,   # options hash
5518        ) = @_;
5519
5520     # Returns nothing.
5521
5522     my ( $num, $record, %count_hash, @records, $tmp_file, $fh_out, $fh_in, $cache );
5523
5524     $tmp_file = "$BP_TMP/count_cache.tmp";
5525
5526     $fh_out   = Maasha::Common::write_open( $tmp_file );
5527
5528     $cache    = 0;
5529     $num      = 0;
5530
5531     while ( $record = get_record( $in ) )
5532     {
5533         map { $count_hash{ $_ }{ $record->{ $_ } }++ if exists $record->{ $_ } } @{ $options->{ "keys" } };
5534
5535         push @records, $record;
5536
5537         if ( scalar @records > 5_000_000 )   # too many records to hold in memory - use disk cache
5538         {
5539             map { put_record( $_, $fh_out ) } @records;
5540
5541             undef @records;
5542
5543             $cache = 1;
5544         }
5545
5546         print STDERR "verbose: records read $num\n" if ( $options->{ 'verbose' } and ( $num % 1_000_000 ) == 0 );
5547
5548         $num++;
5549     }
5550
5551     close $fh_out;
5552
5553     if ( $cache )
5554     {
5555         $num   = 0;
5556
5557         $fh_in = Maasha::Common::read_open( $tmp_file );
5558
5559         while ( $record = get_record( $fh_in ) )
5560         {
5561             map { $record->{ $_ . "_COUNT" } = $count_hash{ $_ }{ $record->{ $_ } } if exists $record->{ $_ } } @{ $options->{ "keys" } };
5562
5563             put_record( $record, $out );
5564
5565             print STDERR "verbose: cache read $num\n" if ( $options->{ 'verbose' } and ( $num % 1_000_000 ) == 0 );
5566
5567             $num++;
5568         }
5569
5570         close $fh_in;
5571     }
5572
5573     foreach $record ( @records )
5574     {
5575         map { $record->{ $_ . "_COUNT" } = $count_hash{ $_ }{ $record->{ $_ } } if exists $record->{ $_ } } @{ $options->{ "keys" } };
5576
5577         put_record( $record, $out );
5578     }
5579
5580     unlink $tmp_file;
5581 }
5582
5583
5584 sub script_plot_histogram
5585 {
5586     # Martin A. Hansen, September 2007.
5587
5588     # Plot a simple histogram for a given key using GNU plot.
5589
5590     my ( $in,        # handle to in stream
5591          $out,       # handle to out stream
5592          $options,   # options hash
5593        ) = @_;
5594
5595     # Returns nothing.
5596
5597     my ( $record, %data_hash, $max, @data_list, $i, $result, $fh );
5598
5599     $options->{ "title" } ||= "Histogram";
5600     $options->{ "sort" }  ||= "num";
5601
5602     while ( $record = get_record( $in ) )
5603     {
5604         $data_hash{ $record->{ $options->{ "key" } } }++ if defined $record->{ $options->{ "key" } };
5605
5606         put_record( $record, $out ) if not $options->{ "no_stream" };
5607     }
5608
5609     if ( $options->{ "sort" } eq "num" ) {
5610         map { push @data_list, [ $_, $data_hash{ $_ } ] } sort { $a <=> $b } keys %data_hash;
5611     } else {
5612         map { push @data_list, [ $_, $data_hash{ $_ } ] } sort keys %data_hash;
5613     }
5614
5615     $result = Maasha::Plot::histogram_simple( \@data_list, $options );
5616
5617     $fh = write_stream( $options->{ "data_out" } );
5618
5619     print $fh "$_\n" foreach @{ $result };
5620
5621     close $fh;
5622 }
5623
5624
5625 sub script_plot_lendist
5626 {
5627     # Martin A. Hansen, August 2007.
5628
5629     # Plot length distribution using GNU plot.
5630
5631     my ( $in,        # handle to in stream
5632          $out,       # handle to out stream
5633          $options,   # options hash
5634        ) = @_;
5635
5636     # Returns nothing.
5637
5638     my ( $record, %data_hash, $max, @data_list, $i, $result, $fh );
5639
5640     $options->{ "title" } ||= "Length Distribution";
5641
5642     while ( $record = get_record( $in ) )
5643     {
5644         $data_hash{ $record->{ $options->{ "key" } } }++ if defined $record->{ $options->{ "key" } };
5645
5646         put_record( $record, $out ) if not $options->{ "no_stream" };
5647     }
5648
5649     $max = Maasha::Calc::list_max( [ keys %data_hash ] );
5650
5651     for ( $i = 0; $i < $max; $i++ ) {
5652         push @data_list, [ $i, $data_hash{ $i } || 0 ];
5653     }
5654
5655     $result = Maasha::Plot::histogram_lendist( \@data_list, $options );
5656
5657     $fh = write_stream( $options->{ "data_out" } );
5658
5659     print $fh "$_\n" foreach @{ $result };
5660
5661     close $fh;
5662 }
5663
5664
5665 sub script_plot_chrdist
5666 {
5667     # Martin A. Hansen, August 2007.
5668
5669     # Plot chromosome distribution using GNU plot.
5670
5671     my ( $in,        # handle to in stream
5672          $out,       # handle to out stream
5673          $options,   # options hash
5674        ) = @_;
5675
5676     # Returns nothing.
5677
5678     my ( $record, %data_hash, @data_list, $elem, $sort_key, $count, $result, $fh );
5679
5680     $options->{ "title" } ||= "Chromosome Distribution";
5681
5682     while ( $record = get_record( $in ) )
5683     {
5684         if ( $record->{ "CHR" } ) {                                                             # generic
5685             $data_hash{ $record->{ "CHR" } }++;
5686         } elsif ( $record->{ "REC_TYPE" } eq "PATSCAN" and $record->{ "S_ID" } =~ /^chr/i ) {   # patscan
5687             $data_hash{ $record->{ "S_ID" } }++;
5688         } elsif ( $record->{ "REC_TYPE" } eq "PSL" and $record->{ "S_ID" } =~ /^chr/i ) {       # BLAT / PSL
5689             $data_hash{ $record->{ "S_ID" } }++;
5690         } elsif ( $record->{ "REC_TYPE" } eq "BLAST" and $record->{ "S_ID" } =~ /^chr/i ) {     # BLAST
5691             $data_hash{ $record->{ "S_ID" } }++;
5692         }
5693
5694         put_record( $record, $out ) if not $options->{ "no_stream" };
5695     }
5696
5697     foreach $elem ( keys %data_hash )
5698     {
5699         $sort_key = $elem;
5700
5701         $sort_key =~ s/chr//i;
5702
5703         $sort_key =~ s/^X(.*)/99$1/;
5704         $sort_key =~ s/^Y(.*)/99$1/;
5705         $sort_key =~ s/^Z(.*)/999$1/;
5706         $sort_key =~ s/^M(.*)/9999$1/;
5707         $sort_key =~ s/^U(.*)/99999$1/;
5708
5709         $count = $sort_key =~ tr/_//;
5710
5711         $sort_key =~ s/_.*/"999999" x $count/ex;
5712
5713         push @data_list, [ $elem, $data_hash{ $elem }, $sort_key ];
5714     }
5715
5716     @data_list = sort { $a->[ 2 ] <=> $b->[ 2 ] } @data_list;
5717
5718     $result = Maasha::Plot::histogram_chrdist( \@data_list, $options );
5719
5720     $fh = write_stream( $options->{ "data_out" } );
5721
5722     print $fh "$_\n" foreach @{ $result };
5723
5724     close $fh;
5725 }
5726
5727
5728 sub script_plot_karyogram
5729 {
5730     # Martin A. Hansen, August 2007.
5731
5732     # Plot hits on karyogram.
5733
5734     my ( $in,        # handle to in stream
5735          $out,       # handle to out stream
5736          $options,   # options hash
5737        ) = @_;
5738
5739     # Returns nothing.
5740
5741     my ( %options, $record, @data, $fh, $result, %data_hash );
5742
5743     $options->{ "genome" }     ||= "human";
5744     $options->{ "feat_color" } ||= "black";
5745
5746     while ( $record = get_record( $in ) )
5747     {
5748         if ( $record->{ "CHR" } and $record->{ "CHR_BEG" } and $record->{ "CHR_END" } )
5749         {
5750             push @{ $data_hash{ $record->{ "CHR" } } }, [ $record->{ "CHR_BEG" }, $record->{ "CHR_END" }, $options->{ "feat_color" } ];
5751         }
5752
5753         put_record( $record, $out ) if not $options->{ "no_stream" };
5754     }
5755
5756     $result = Maasha::Plot::karyogram( \%data_hash, \%options );
5757
5758     $fh = write_stream( $options->{ "data_out" } );
5759
5760     print $fh $result;
5761
5762     close $fh;
5763 }
5764
5765
5766 sub script_plot_matches
5767 {
5768     # Martin A. Hansen, August 2007.
5769
5770     # Plot matches in 2D generating a dotplot.
5771
5772     my ( $in,        # handle to in stream
5773          $out,       # handle to out stream
5774          $options,   # options hash
5775        ) = @_;
5776
5777     # Returns nothing.
5778
5779     my ( $record, @data, $fh, $result, %data_hash );
5780
5781     $options->{ "direction" } ||= "both";
5782
5783     while ( $record = get_record( $in ) )
5784     {
5785         if ( defined $record->{ "Q_BEG" } and defined $record->{ "S_BEG" } and $record->{ "Q_END" } and $record->{ "S_END" } ) {
5786             push @data, $record;
5787         }
5788
5789         put_record( $record, $out ) if not $options->{ "no_stream" };
5790     }
5791
5792     $options->{ "title" }  ||= "plot_matches";
5793     $options->{ "xlabel" } ||= $data[ 0 ]->{ "Q_ID" };
5794     $options->{ "ylabel" } ||= $data[ 0 ]->{ "S_ID" };
5795
5796     $result = Maasha::Plot::dotplot_matches( \@data, $options, $BP_TMP );
5797
5798     $fh = write_stream( $options->{ "data_out" } );
5799
5800     print $fh "$_\n" foreach @{ $result };
5801
5802     close $fh;
5803 }
5804
5805
5806 sub script_length_vals
5807 {
5808     # Martin A. Hansen, August 2007.
5809
5810     # Determine the length of the value for given keys.
5811
5812     my ( $in,        # handle to in stream
5813          $out,       # handle to out stream
5814          $options,   # options hash
5815        ) = @_;
5816
5817     # Returns nothing.
5818
5819     my ( $record, $key );
5820
5821     while ( $record = get_record( $in ) )
5822     {
5823         foreach $key ( @{ $options->{ "keys" } } )
5824         {
5825             if ( $record->{ $key } ) {
5826                 $record->{ $key . "_LEN" } = length $record->{ $key };
5827             }
5828         }
5829
5830         put_record( $record, $out );
5831     }
5832 }
5833
5834
5835 sub script_sum_vals
5836 {
5837     # Martin A. Hansen, August 2007.
5838
5839     # Calculates the sums for values of given keys.
5840
5841     my ( $in,        # handle to in stream
5842          $out,       # handle to out stream
5843          $options,   # options hash
5844        ) = @_;
5845
5846     # Returns nothing.
5847
5848     my ( $record, $key, %sum_hash, $fh );
5849
5850     while ( $record = get_record( $in ) )
5851     {
5852         foreach $key ( @{ $options->{ "keys" } } )
5853         {
5854             if ( $record->{ $key } ) {
5855                 $sum_hash{ $key } += $record->{ $key };
5856             }
5857         }
5858
5859         put_record( $record, $out ) if not $options->{ "no_stream" };
5860     }
5861
5862     $fh = write_stream( $options->{ "data_out" } );
5863
5864     foreach $key ( @{ $options->{ "keys" } } ) {
5865         put_record( { $key . "_SUM" => $sum_hash{ $key } || 0 } , $fh );
5866     }
5867
5868     close $fh;
5869 }
5870
5871
5872 sub script_mean_vals
5873 {
5874     # Martin A. Hansen, August 2007.
5875
5876     # Calculate the mean of values of given keys.
5877
5878     my ( $in,        # handle to in stream
5879          $out,       # handle to out stream
5880          $options,   # options hash
5881        ) = @_;
5882
5883     # Returns nothing.
5884
5885     my ( $record, $key, %sum_hash, %count_hash, $mean, $fh );
5886
5887     while ( $record = get_record( $in ) )
5888     {
5889         foreach $key ( @{ $options->{ "keys" } } )
5890         {
5891             if ( $record->{ $key } )
5892             {
5893                 $sum_hash{ $key } += $record->{ $key };
5894                 $count_hash{ $key }++;
5895             }
5896         }
5897
5898         put_record( $record, $out ) if not $options->{ "no_stream" };
5899     }
5900
5901     $fh = write_stream( $options->{ "data_out" } );
5902
5903     foreach $key ( @{ $options->{ "keys" } } )
5904     {
5905         if ( $count_hash{ $key } ) {
5906             $mean = sprintf( "%.2f", ( $sum_hash{ $key } / $count_hash{ $key } ) );
5907         } else {
5908             $mean = "N/A";
5909         }
5910
5911         put_record( { $key . "_MEAN" => $mean } , $fh );
5912     }
5913
5914     close $fh;
5915 }
5916
5917
5918 sub script_median_vals
5919 {
5920     # Martin A. Hansen, March 2008.
5921
5922     # Calculate the median values of given keys.
5923
5924     my ( $in,        # handle to in stream
5925          $out,       # handle to out stream
5926          $options,   # options hash
5927        ) = @_;
5928
5929     # Returns nothing.
5930
5931     my ( $record, $key, %median_hash, $median, $fh );
5932
5933     while ( $record = get_record( $in ) )
5934     {
5935         foreach $key ( @{ $options->{ "keys" } } ) {
5936             push @{ $median_hash{ $key } }, $record->{ $key } if defined $record->{ $key };
5937         }
5938
5939         put_record( $record, $out ) if not $options->{ "no_stream" };
5940     }
5941
5942     $fh = write_stream( $options->{ "data_out" } );
5943
5944     foreach $key ( @{ $options->{ "keys" } } )
5945     {
5946         if ( $median_hash{ $key } ) {
5947             $median = Maasha::Calc::median( $median_hash{ $key } );
5948         } else {
5949             $median = "N/A";
5950         }
5951
5952         put_record( { $key . "_MEDIAN" => $median } , $fh );
5953     }
5954
5955     close $fh;
5956 }
5957
5958
5959 sub script_max_vals
5960 {
5961     # Martin A. Hansen, February 2008.
5962
5963     # Determine the maximum values of given keys.
5964
5965     my ( $in,        # handle to in stream
5966          $out,       # handle to out stream
5967          $options,   # options hash
5968        ) = @_;
5969
5970     # Returns nothing.
5971
5972     my ( $record, $key, $fh, %max_hash, $max_record );
5973
5974     while ( $record = get_record( $in ) )
5975     {
5976         foreach $key ( @{ $options->{ "keys" } } )
5977         {
5978             if ( $record->{ $key } )
5979             {
5980                 $max_hash{ $key } = $record->{ $key } if $record->{ $key } > $max_hash{ $key };
5981             }
5982         }
5983
5984         put_record( $record, $out ) if not $options->{ "no_stream" };
5985     }
5986
5987     $fh = write_stream( $options->{ "data_out" } );
5988
5989     foreach $key ( @{ $options->{ "keys" } } )
5990     {
5991         $max_record->{ $key . "_MAX" } = $max_hash{ $key };
5992     }
5993
5994     put_record( $max_record, $fh );
5995
5996     close $fh;
5997 }
5998
5999
6000 sub script_min_vals
6001 {
6002     # Martin A. Hansen, February 2008.
6003
6004     # Determine the minimum values of given keys.
6005
6006     my ( $in,        # handle to in stream
6007          $out,       # handle to out stream
6008          $options,   # options hash
6009        ) = @_;
6010
6011     # Returns nothing.
6012
6013     my ( $record, $key, $fh, %min_hash, $min_record );
6014
6015     while ( $record = get_record( $in ) )
6016     {
6017         foreach $key ( @{ $options->{ "keys" } } )
6018         {
6019             if ( defined $record->{ $key } )
6020             {
6021                 if ( exists $min_hash{ $key } ) {
6022                     $min_hash{ $key } = $record->{ $key } if $record->{ $key } < $min_hash{ $key };
6023                 } else {
6024                     $min_hash{ $key } = $record->{ $key };
6025                 }
6026             }
6027         }
6028
6029         put_record( $record, $out ) if not $options->{ "no_stream" };
6030     }
6031
6032     $fh = write_stream( $options->{ "data_out" } );
6033
6034     foreach $key ( @{ $options->{ "keys" } } )
6035     {
6036         $min_record->{ $key . "_MIN" } = $min_hash{ $key };
6037     }
6038
6039     put_record( $min_record, $fh );
6040
6041     close $fh;
6042 }
6043
6044
6045 sub script_upload_to_ucsc
6046 {
6047     # Martin A. Hansen, August 2007.
6048
6049     # Calculate the mean of values of given keys.
6050
6051     # This routine has developed into the most ugly hack. Do something!
6052
6053     my ( $in,        # handle to in stream
6054          $out,       # handle to out stream
6055          $options,   # options hash
6056        ) = @_;
6057
6058     # Returns nothing.
6059
6060     my ( $record, $file, $wib_file, $wig_file, $wib_dir, $fh_out, $i, $first, $format, $type, $columns, $append, $vals );
6061
6062     $options->{ "short_label" } ||= $options->{ 'table' };
6063     $options->{ "long_label" }  ||= $options->{ 'table' };
6064     $options->{ "group" }       ||= $ENV{ "LOGNAME" };
6065     $options->{ "priority" }    ||= 1;
6066     $options->{ "visibility" }  ||= "pack";
6067     $options->{ "color" }       ||= join( ",", int( rand( 255 ) ), int( rand( 255 ) ), int( rand( 255 ) ) );
6068     $options->{ "chunk_size" }  ||= 10_000_000_000;    # Due to 32-bit UCSC compilation really large tables cannot be loaded in one go.
6069
6070     $file = "$BP_TMP/ucsc_upload.tmp";
6071
6072     $append = 0;
6073
6074     $first = 1;
6075
6076     $i = 0;
6077
6078     $fh_out = Maasha::Common::write_open( $file );
6079
6080     while ( $record = get_record( $in ) )
6081     {
6082         put_record( $record, $out ) if not $options->{ "no_stream" };
6083
6084         if ( $record->{ "REC_TYPE" } eq "fixed_step" )
6085         {
6086             $vals = $record->{ "VALS" };
6087             $vals =~ tr/;/\n/;
6088
6089             print $fh_out "fixedStep chrom=$record->{ 'CHR' } start=$record->{ 'CHR_BEG' } step=$record->{ 'STEP' }\n";
6090             print $fh_out "$vals\n";
6091
6092             $format = "WIGGLE" if not $format;
6093         }
6094         elsif ( $record->{ "REC_TYPE" } eq "PSL" )
6095         {
6096             Maasha::UCSC::psl_put_header( $fh_out ) if $first;
6097             Maasha::UCSC::psl_put_entry( $record, $fh_out );
6098
6099             $first = 0;
6100
6101             $format = "PSL" if not $format;
6102         }
6103         elsif ( $record->{ "REC_TYPE" } eq "BED" and $record->{ "SEC_STRUCT" } )
6104         {
6105             # chrom chromStart  chromEnd    name    score   strand  size    secStr  conf
6106
6107             print $fh_out join ( "\t",
6108                 $record->{ "CHR" },
6109                 $record->{ "CHR_BEG" },
6110                 $record->{ "CHR_END" } + 1,
6111                 $record->{ "Q_ID" },
6112                 $record->{ "SCORE" },
6113                 $record->{ "STRAND" },
6114                 $record->{ "SIZE" },
6115                 $record->{ "SEC_STRUCT" },
6116                 $record->{ "CONF" },
6117             ), "\n";
6118
6119             $format  = "BED_SS" if not $format;
6120         }
6121         elsif ( $record->{ "REC_TYPE" } eq "BED" )
6122         {
6123             Maasha::UCSC::bed_put_entry( $record, $fh_out, $record->{ "BED_COLS" } );
6124
6125             $format  = "BED"                   if not $format;
6126             $columns = $record->{ "BED_COLS" } if not $columns;
6127         }
6128         elsif ( $record->{ "REC_TYPE" } eq "PATSCAN" and $record->{ "CHR" } )
6129         {
6130             Maasha::UCSC::bed_put_entry( $record, $fh_out, 6 );
6131
6132             $format  = "BED" if not $format;
6133             $columns = 6     if not $columns;
6134         }
6135         elsif ( $record->{ "REC_TYPE" } eq "BLAST" and $record->{ "S_ID" } =~ /^chr/ )
6136         {
6137             $record->{ "CHR" }     = $record->{ "S_ID" };
6138             $record->{ "CHR_BEG" } = $record->{ "S_BEG" };
6139             $record->{ "CHR_END" } = $record->{ "S_END" };
6140             $record->{ "SCORE" }   = $record->{ "BIT_SCORE" } * 1000;
6141
6142             $format  = "BED" if not $format;
6143             $columns = 6     if not $columns;
6144
6145             Maasha::UCSC::bed_put_entry( $record, $fh_out );
6146         }
6147         elsif ( $record->{ "REC_TYPE" } eq "VMATCH" and $record->{ "S_ID" } =~ /^chr/i )
6148         {
6149             $record->{ "CHR" }     = $record->{ "S_ID" };
6150             $record->{ "CHR_BEG" } = $record->{ "S_BEG" };
6151             $record->{ "CHR_END" } = $record->{ "S_END" };
6152             $record->{ "SCORE" }   = $record->{ "SCORE" } || 999;
6153             $record->{ "SCORE" }   = int( $record->{ "SCORE" } );
6154
6155             $format  = "BED" if not $format;
6156             $columns = 6     if not $columns;
6157
6158             Maasha::UCSC::bed_put_entry( $record, $fh_out, 6 );
6159         }
6160
6161         if ( $i == $options->{ "chunk_size" } )
6162         {
6163             close $fh_out;
6164
6165             if ( $format eq "BED" ) {
6166                 Maasha::UCSC::bed_upload_to_ucsc( $BP_TMP, $file, $options, $append );
6167             } elsif ( $format eq "PSL" ) {
6168                 Maasha::UCSC::psl_upload_to_ucsc( $file, $options, $append );
6169             }
6170
6171             unlink $file;
6172
6173             $first = 1;
6174
6175             $append = 1;
6176
6177             $fh_out = Maasha::Common::write_open( $file );
6178         }
6179
6180         $i++;
6181     }
6182
6183     close $fh_out;
6184
6185     if ( exists $options->{ "database" } and $options->{ "table" } )
6186     {
6187         if ( $format eq "BED" )
6188         {
6189             $type = "bed $columns";
6190
6191             Maasha::UCSC::bed_upload_to_ucsc( $BP_TMP, $file, $options, $append );
6192         }
6193         elsif ( $format eq "BED_SS" )
6194         {
6195             $options->{ "sec_struct" } = 1;
6196
6197             $type = "sec_struct";
6198
6199             Maasha::UCSC::bed_upload_to_ucsc( $BP_TMP, $file, $options, $append );
6200         }
6201         elsif ( $format eq "PSL" )
6202         {
6203             $type = "psl";
6204
6205             Maasha::UCSC::psl_upload_to_ucsc( $file, $options, $append );
6206         }
6207         elsif ( $format eq "WIGGLE" )
6208         {
6209             $options->{ "visibility" } = "full";
6210
6211             $wig_file = "$options->{ 'table' }.wig";
6212             $wib_file = "$options->{ 'table' }.wib";
6213
6214             $wib_dir  = "$ENV{ 'HOME' }/ucsc/wib";
6215
6216             Maasha::Common::dir_create_if_not_exists( $wib_dir );
6217
6218             if ( $options->{ 'verbose' } ) {
6219                 `cd $BP_TMP && wigEncode $file $wig_file $wib_file`;
6220             } else {
6221                 `cd $BP_TMP && wigEncode $file $wig_file $wib_file > /dev/null 2>&1`;
6222             }
6223
6224             Maasha::Common::run( "mv", "$BP_TMP/$wib_file $wib_dir" );
6225
6226             unlink $file;
6227
6228             $file = $wig_file;
6229
6230             $type = "wig 0";
6231
6232             Maasha::UCSC::wiggle_upload_to_ucsc( $BP_TMP, $wib_dir, $file, $options );
6233         }
6234
6235         unlink $file;
6236
6237         Maasha::UCSC::update_my_tracks( $options, $type );
6238     }
6239 }
6240
6241
6242 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
6243
6244
6245 sub record2fasta
6246 {
6247     # Martin A. Hansen, July 2008.
6248
6249     # Given a biopiece record converts it to a FASTA record.
6250     # If no generic SEQ or SEQ_NAME is found, the Q_* and S_* are
6251     # tried in that order.
6252
6253     my ( $record,    # record
6254        ) = @_;
6255
6256     # Returns a tuple.
6257
6258     my ( $seq_name, $seq );
6259
6260     $seq_name = $record->{ "SEQ_NAME" } || $record->{ "Q_ID" }  || $record->{ "S_ID" };
6261     $seq      = $record->{ "SEQ" }      || $record->{ "Q_SEQ" } || $record->{ "S_SEQ" };
6262
6263     if ( defined $seq_name and defined $seq ) {
6264         return wantarray ? ( $seq_name, $seq ) : [ $seq_name, $seq ];
6265     } else {
6266         return;
6267     }
6268 }
6269
6270
6271 sub read_stream
6272 {
6273     # Martin A. Hansen, July 2007.
6274
6275     # Opens a stream to STDIN or a file,
6276
6277     my ( $path,   # path - OPTIONAL
6278        ) = @_;
6279
6280     # Returns filehandle.
6281
6282     my ( $fh );
6283
6284     if ( not -t STDIN ) {
6285         $fh = Maasha::Common::read_stdin();
6286     } elsif ( not $path ) {
6287 #        Maasha::Common::error( qq(no data stream) );
6288     } else {
6289         $fh = Maasha::Common::read_open( $path );
6290     }
6291
6292 #    $fh->autoflush(1) if $fh;  # Disable file buffer for debugging.
6293
6294     return $fh;
6295 }
6296
6297
6298 sub write_stream
6299 {
6300     # Martin A. Hansen, August 2007.
6301
6302     # Opens a stream to STDOUT or a file.
6303
6304     my ( $path,   # path          - OPTIONAL
6305          $gzip,   # compress data - OPTIONAL
6306        ) = @_;
6307
6308     # Returns filehandle.
6309
6310     my ( $fh );
6311
6312     if ( $path ) {
6313         $fh = Maasha::Common::write_open( $path, $gzip );
6314     } else {
6315         $fh = Maasha::Common::write_stdout();
6316     }
6317
6318     return $fh;
6319 }
6320
6321
6322 sub get_record
6323 {
6324     # Martin A. Hansen, July 2007.
6325
6326     # Reads one record at a time and converts that record
6327     # to a Perl data structure (a hash) which is returned.
6328
6329     my ( $fh,   # handle to stream
6330        ) = @_;
6331
6332     # Returns a hash.
6333
6334     my ( $block, @lines, $line, $key, $value, %record );
6335
6336     local $/ = "\n---\n";
6337
6338     $block = <$fh>;
6339
6340     chomp $block;
6341
6342     return if not defined $block;
6343
6344     @lines = split "\n", $block;
6345
6346     foreach $line ( @lines )
6347     {
6348         ( $key, $value ) = split ": ", $line, 2;
6349
6350         $record{ $key } = $value;
6351     }
6352
6353     return wantarray ? %record : \%record;
6354 }
6355
6356
6357 sub put_record
6358 {
6359     # Martin A. Hansen, July 2007.
6360
6361     # Given a Perl datastructure (a hash ref) emits this to STDOUT or a filehandle.
6362
6363     my ( $data,   # data structure
6364          $fh,     # file handle - OPTIONAL
6365        ) = @_;
6366
6367     # Returns nothing.
6368
6369     if ( scalar keys %{ $data } )
6370     {
6371         if ( $fh )
6372         {
6373             map { print $fh "$_: $data->{ $_ }\n" } keys %{ $data };
6374             print $fh "---\n";
6375         }
6376         else
6377         {
6378             map { print "$_: $data->{ $_ }\n" } keys %{ $data };
6379             print "---\n";
6380         }
6381     }
6382
6383     undef $data;
6384 }
6385
6386
6387 sub getopt_files
6388 {
6389     # Martin A. Hansen, November 2007.
6390
6391     # Extracts files from an explicit GetOpt::Long argument
6392     # allowing for the use of glob. E.g.
6393     # --data_in=test.fna
6394     # --data_in=test.fna,test2.fna
6395     # --data_in=*.fna
6396     # --data_in=test.fna,/dir/*.fna
6397
6398     my ( $option,   # option from GetOpt::Long
6399        ) = @_;
6400
6401     # Returns a list.
6402
6403     my ( $elem, @files );
6404
6405     foreach $elem ( split ",", $option )
6406     {
6407         if ( -f $elem ) {
6408             push @files, $elem;
6409         } elsif ( $elem =~ /\*/ ) {
6410             push @files, glob( $elem );
6411         }
6412     }
6413
6414     return wantarray ? @files : \@files;
6415 }
6416
6417
6418 sub sig_handler
6419 {
6420     # Martin A. Hansen, April 2008.
6421
6422     # Removes temporary directory and exits gracefully.
6423     # This subroutine is meant to be run always as the last
6424     # thing even if a script is dies or is interrupted
6425     # or killed.
6426
6427     my ( $sig,   # signal from the %SIG
6428        ) = @_;
6429
6430     # print STDERR "signal->$sig<-\n";
6431
6432     chomp $sig;
6433
6434     sleep 1;
6435
6436     if ( -d $BP_TMP )
6437     {
6438         if ( $sig =~ /MAASHA_ERROR/ ) {
6439             print STDERR "\nProgram '$script' had an error"                     . "  -  Please wait for temporary data to be removed\n";
6440         } elsif ( $sig eq "INT" ) {
6441             print STDERR "\nProgram '$script' interrupted (ctrl-c was pressed)" . "  -  Please wait for temporary data to be removed\n";
6442         } elsif ( $sig eq "TERM" ) {
6443             print STDERR "\nProgram '$script' terminated (someone used kill?)"  . "  -  Please wait for temporary data to be removed\n";
6444         } else {
6445             print STDERR "\nProgram '$script' died->$sig"                       . "  -  Please wait for temporary data to be removed\n";
6446         }
6447
6448         clean_tmp();
6449     }
6450
6451     exit( 0 );
6452 }
6453
6454
6455 sub clean_tmp
6456 {
6457     # Martin A. Hansen, July 2008.
6458
6459     # Cleans out any unused temporary files and directories in BP_TMP.
6460
6461     # Returns nothing.
6462
6463     my ( $tmpdir, @dirs, $curr_pid, $dir, $user, $sid, $pid );
6464
6465     $tmpdir = $ENV{ 'BP_TMP' } || Maasha::Common::error( 'No BP_TMP variable in environment.' );
6466
6467     $curr_pid = Maasha::Common::get_processid();
6468
6469     @dirs = Maasha::Common::ls_dirs( $tmpdir );
6470
6471     foreach $dir ( @dirs )
6472     {
6473         if ( $dir =~ /^$tmpdir\/(.+)_(\d+)_(\d+)_bp_tmp$/ )
6474         {
6475             $user = $1;
6476             $sid  = $2;
6477             $pid  = $3;
6478
6479             if ( $user eq Maasha::Common::get_user() )
6480             {
6481                 if ( not Maasha::Common::process_running( $pid ) )
6482                 {
6483                     # print STDERR "Removing stale dir: $dir\n";
6484                     Maasha::Common::dir_remove( $dir );
6485                 }
6486                 elsif ( $pid == $curr_pid )
6487                 {
6488                     # print STDERR "Removing current dir: $dir\n";
6489                     # Maasha::Common::dir_remove( $dir );
6490                 }
6491             }
6492         }
6493     }
6494 }
6495
6496
6497 END
6498 {
6499     clean_tmp();
6500 }
6501
6502
6503 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
6504
6505 1;
6506
6507 __END__