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