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