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