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