code_perl/Maasha/Biopieces.pm

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