code_perl/Maasha/BioRun.pm

   1 package Maasha::BioRun;
   2
   3
   4 # Copyright (C) 2007-2009 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 that contains Biopieces which are run.
  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::Biopieces;
  38 use Maasha::Config;
  39 use Maasha::Common;
  40 use Maasha::Filesys;
  41 use Maasha::Fasta;
  42 use Maasha::Align;
  43 use Maasha::Matrix;
  44 use Maasha::Match;
  45 use Maasha::EMBL;
  46 use Maasha::Stockholm;
  47 use Maasha::Seq;
  48 use Maasha::Patscan;
  49 use Maasha::Plot;
  50 use Maasha::Calc;
  51 use Maasha::UCSC;
  52 use Maasha::UCSC::BED;
  53 use Maasha::UCSC::Wiggle;
  54 use Maasha::NCBI;
  55 use Maasha::GFF;
  56 use Maasha::TwoBit;
  57 use Maasha::Solid;
  58 use Maasha::Solexa;
  59 use Maasha::SQL;
  60 use Maasha::Gwiki;
  61
  62 use vars qw( @ISA @EXPORT_OK );
  63
  64 require Exporter;
  65
  66 @ISA = qw( Exporter );
  67
  68 use constant {
  69     SEQ_NAME => 0,
  70     SEQ      => 1,
  71 };
  72
  73
  74 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> GLOBALS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  75
  76
  77 my ( $script, $BP_TMP );
  78
  79 $script = Maasha::Common::get_scriptname();
  80 $BP_TMP = Maasha::Common::get_tmpdir();
  81
  82
  83 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> RUN SCRIPT <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  84
  85
  86 run_script( $script );
  87
  88
  89 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> SUBROUTINES <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  90
  91
  92 sub run_script
  93 {
  94     # Martin A. Hansen, August 2007.
  95
  96     # Run a specific script.
  97
  98     my ( $script,   # script name
  99        ) = @_;
 100
 101     # Returns nothing.
 102
 103     my ( $t0, $t1, $options, $in, $out );
 104
 105     Maasha::Biopieces::log_biopiece();
 106
 107     $t0 = gettimeofday();
 108
 109     $options = get_options( $script );
 110
 111     $options->{ "SCRIPT" } = $script;
 112
 113     $script = "print_usage" if ( -t STDIN and keys %{ $options } <= 1 or $options->{ 'help' } );
 114
 115     $in  = Maasha::Biopieces::read_stream( $options->{ "stream_in" } );
 116     $out = Maasha::Biopieces::write_stream( $options->{ "stream_out" } );
 117
 118     if    ( $script eq "print_usage" )              { script_print_usage(               $in, $out, $options ) }
 119     elsif ( $script eq "read_psl" )                 { script_read_psl(                  $in, $out, $options ) }
 120     elsif ( $script eq "read_bed" )                 { script_read_bed(                  $in, $out, $options ) }
 121     elsif ( $script eq "read_fixedstep" )           { script_read_fixedstep(            $in, $out, $options ) }
 122     elsif ( $script eq "read_blast_tab" )           { script_read_blast_tab(            $in, $out, $options ) }
 123     elsif ( $script eq "read_embl" )                { script_read_embl(                 $in, $out, $options ) }
 124     elsif ( $script eq "read_stockholm" )           { script_read_stockholm(            $in, $out, $options ) }
 125     elsif ( $script eq "read_phastcons" )           { script_read_phastcons(            $in, $out, $options ) }
 126     elsif ( $script eq "read_soft" )                { script_read_soft(                 $in, $out, $options ) }
 127     elsif ( $script eq "read_gff" )                 { script_read_gff(                  $in, $out, $options ) }
 128     elsif ( $script eq "read_2bit" )                { script_read_2bit(                 $in, $out, $options ) }
 129     elsif ( $script eq "read_solexa" )              { script_read_solexa(               $in, $out, $options ) }
 130     elsif ( $script eq "read_solid" )               { script_read_solid(                $in, $out, $options ) }
 131     elsif ( $script eq "read_mysql" )               { script_read_mysql(                $in, $out, $options ) }
 132     elsif ( $script eq "read_ucsc_config" )         { script_read_ucsc_config(          $in, $out, $options ) }
 133     elsif ( $script eq "length_seq" )               { script_length_seq(                $in, $out, $options ) }
 134     elsif ( $script eq "uppercase_seq" )            { script_uppercase_seq(             $in, $out, $options ) }
 135     elsif ( $script eq "complexity_seq" )           { script_complexity_seq(            $in, $out, $options ) }
 136     elsif ( $script eq "oligo_freq" )               { script_oligo_freq(                $in, $out, $options ) }
 137     elsif ( $script eq "remove_indels" )            { script_remove_indels(             $in, $out, $options ) }
 138     elsif ( $script eq "transliterate_seq" )        { script_transliterate_seq(         $in, $out, $options ) }
 139     elsif ( $script eq "transliterate_vals" )       { script_transliterate_vals(        $in, $out, $options ) }
 140     elsif ( $script eq "translate_seq" )            { script_translate_seq(             $in, $out, $options ) }
 141     elsif ( $script eq "get_genome_align" )         { script_get_genome_align(          $in, $out, $options ) }
 142     elsif ( $script eq "get_genome_phastcons" )     { script_get_genome_phastcons(      $in, $out, $options ) }
 143     elsif ( $script eq "fold_seq" )                 { script_fold_seq(                  $in, $out, $options ) }
 144     elsif ( $script eq "tile_seq" )                 { script_tile_seq(                  $in, $out, $options ) }
 145     elsif ( $script eq "invert_align" )             { script_invert_align(              $in, $out, $options ) }
 146     elsif ( $script eq "patscan_seq" )              { script_patscan_seq(               $in, $out, $options ) }
 147     elsif ( $script eq "soap_seq" )                 { script_soap_seq(                  $in, $out, $options ) }
 148     elsif ( $script eq "match_seq" )                { script_match_seq(                 $in, $out, $options ) }
 149     elsif ( $script eq "create_vmatch_index" )      { script_create_vmatch_index(       $in, $out, $options ) }
 150     elsif ( $script eq "write_bed" )                { script_write_bed(                 $in, $out, $options ) }
 151     elsif ( $script eq "write_psl" )                { script_write_psl(                 $in, $out, $options ) }
 152     elsif ( $script eq "write_fixedstep" )          { script_write_fixedstep(           $in, $out, $options ) }
 153     elsif ( $script eq "write_2bit" )               { script_write_2bit(                $in, $out, $options ) }
 154     elsif ( $script eq "write_solid" )              { script_write_solid(               $in, $out, $options ) }
 155     elsif ( $script eq "write_ucsc_config" )        { script_write_ucsc_config(         $in, $out, $options ) }
 156     elsif ( $script eq "head_records" )             { script_head_records(              $in, $out, $options ) }
 157     elsif ( $script eq "remove_keys" )              { script_remove_keys(               $in, $out, $options ) }
 158     elsif ( $script eq "remove_adaptor" )           { script_remove_adaptor(            $in, $out, $options ) }
 159     elsif ( $script eq "remove_mysql_tables" )      { script_remove_mysql_tables(       $in, $out, $options ) }
 160     elsif ( $script eq "remove_ucsc_tracks" )       { script_remove_ucsc_tracks(        $in, $out, $options ) }
 161     elsif ( $script eq "rename_keys" )              { script_rename_keys(               $in, $out, $options ) }
 162     elsif ( $script eq "uniq_vals" )                { script_uniq_vals(                 $in, $out, $options ) }
 163     elsif ( $script eq "merge_vals" )               { script_merge_vals(                $in, $out, $options ) }
 164     elsif ( $script eq "merge_records" )            { script_merge_records(             $in, $out, $options ) }
 165     elsif ( $script eq "flip_tab" )                 { script_flip_tab(                  $in, $out, $options ) }
 166     elsif ( $script eq "sort_records" )             { script_sort_records(              $in, $out, $options ) }
 167     elsif ( $script eq "plot_histogram" )           { script_plot_histogram(            $in, $out, $options ) }
 168     elsif ( $script eq "plot_lendist" )             { script_plot_lendist(              $in, $out, $options ) }
 169     elsif ( $script eq "plot_chrdist" )             { script_plot_chrdist(              $in, $out, $options ) }
 170     elsif ( $script eq "plot_karyogram" )           { script_plot_karyogram(            $in, $out, $options ) }
 171     elsif ( $script eq "plot_matches" )             { script_plot_matches(              $in, $out, $options ) }
 172     elsif ( $script eq "plot_seqlogo" )             { script_plot_seqlogo(              $in, $out, $options ) }
 173     elsif ( $script eq "plot_phastcons_profiles" )  { script_plot_phastcons_profiles(   $in, $out, $options ) }
 174     elsif ( $script eq "length_vals" )              { script_length_vals(               $in, $out, $options ) }
 175     elsif ( $script eq "sum_vals" )                 { script_sum_vals(                  $in, $out, $options ) }
 176     elsif ( $script eq "mean_vals" )                { script_mean_vals(                 $in, $out, $options ) }
 177     elsif ( $script eq "median_vals" )              { script_median_vals(               $in, $out, $options ) }
 178     elsif ( $script eq "max_vals" )                 { script_max_vals(                  $in, $out, $options ) }
 179     elsif ( $script eq "min_vals" )                 { script_min_vals(                  $in, $out, $options ) }
 180     elsif ( $script eq "upload_to_ucsc" )           { script_upload_to_ucsc(            $in, $out, $options ) }
 181
 182     close $in if defined $in;
 183     close $out;
 184
 185     $t1 = gettimeofday();
 186
 187     print STDERR "Program: $script" . ( " " x ( 25 - length( $script ) ) ) . sprintf( "Run time: %.4f\n", ( $t1 - $t0 ) ) if $options->{ 'verbose' };
 188 }
 189
 190
 191 sub get_options
 192 {
 193     # Martin A. Hansen, February 2008.
 194
 195     # Gets options from commandline and checks these vigerously.
 196
 197     my ( $script,     # name of script
 198        ) = @_;
 199
 200     # Returns hash
 201
 202     my ( %options, @options, $opt, @genomes, $real );
 203
 204     if ( $script eq "print_usage" )
 205     {
 206         @options = qw(
 207             data_in|i=s
 208         );
 209     }
 210     elsif ( $script eq "read_psl" )
 211     {
 212         @options = qw(
 213             data_in|i=s
 214             num|n=s
 215         );
 216     }
 217     elsif ( $script eq "read_bed" )
 218     {
 219         @options = qw(
 220             data_in|i=s
 221             cols|c=s
 222             num|n=s
 223             check|C
 224         );
 225     }
 226     elsif ( $script eq "read_fixedstep" )
 227     {
 228         @options = qw(
 229             data_in|i=s
 230             num|n=s
 231         );
 232     }
 233     elsif ( $script eq "read_blast_tab" )
 234     {
 235         @options = qw(
 236             data_in|i=s
 237             num|n=s
 238         );
 239     }
 240     elsif ( $script eq "read_embl" )
 241     {
 242         @options = qw(
 243             data_in|i=s
 244             num|n=s
 245             keys|k=s
 246             feats|f=s
 247             quals|q=s
 248         );
 249     }
 250     elsif ( $script eq "read_stockholm" )
 251     {
 252         @options = qw(
 253             data_in|i=s
 254             num|n=s
 255         );
 256     }
 257     elsif ( $script eq "read_phastcons" )
 258     {
 259         @options = qw(
 260             data_in|i=s
 261             num|n=s
 262             min|m=s
 263             dist|d=s
 264             threshold|t=f
 265             gap|g=s
 266         );
 267     }
 268     elsif ( $script eq "read_soft" )
 269     {
 270         @options = qw(
 271             data_in|i=s
 272             samples|s=s
 273             num|n=s
 274         );
 275     }
 276     elsif ( $script eq "read_gff" )
 277     {
 278         @options = qw(
 279             data_in|i=s
 280             num|n=s
 281         );
 282     }
 283     elsif ( $script eq "read_2bit" )
 284     {
 285         @options = qw(
 286             data_in|i=s
 287             num|n=s
 288             no_mask|N
 289         );
 290     }
 291     elsif ( $script eq "read_solexa" )
 292     {
 293         @options = qw(
 294             data_in|i=s
 295             num|n=s
 296             format|f=s
 297             quality|q=s
 298         );
 299     }
 300     elsif ( $script eq "read_solid" )
 301     {
 302         @options = qw(
 303             data_in|i=s
 304             num|n=s
 305             quality|q=s
 306         );
 307     }
 308     elsif ( $script eq "read_mysql" )
 309     {
 310         @options = qw(
 311             database|d=s
 312             query|q=s
 313             user|u=s
 314             password|p=s
 315         );
 316     }
 317     elsif ( $script eq "read_ucsc_config" )
 318     {
 319         @options = qw(
 320             data_in|i=s
 321             num|n=s
 322         );
 323     }
 324     elsif ( $script eq "length_seq" )
 325     {
 326         @options = qw(
 327             no_stream|x
 328             data_out|o=s
 329         );
 330     }
 331     elsif ( $script eq "oligo_freq" )
 332     {
 333         @options = qw(
 334             word_size|w=s
 335             all|a
 336         );
 337     }
 338     elsif ( $script eq "transliterate_seq" )
 339     {
 340         @options = qw(
 341             search|s=s
 342             replace|r=s
 343             delete|d=s
 344         );
 345     }
 346     elsif ( $script eq "transliterate_vals" )
 347     {
 348         @options = qw(
 349             keys|k=s
 350             search|s=s
 351             replace|r=s
 352             delete|d=s
 353         );
 354     }
 355     elsif ( $script eq "translate_seq" )
 356     {
 357         @options = qw(
 358             frames|f=s
 359         );
 360     }
 361     elsif ( $script eq "get_genome_align" )
 362     {
 363         @options = qw(
 364             genome|g=s
 365             chr|c=s
 366             beg|b=s
 367             end|e=s
 368             len|l=s
 369             strand|s=s
 370         );
 371     }
 372     elsif ( $script eq "get_genome_phastcons" )
 373     {
 374         @options = qw(
 375             genome|g=s
 376             chr|c=s
 377             beg|b=s
 378             end|e=s
 379             len|l=s
 380             flank|f=s
 381         );
 382     }
 383     elsif ( $script eq "tile_seq" )
 384     {
 385         @options = qw(
 386             identity|i=s
 387             supress_indels|s
 388         );
 389     }
 390     elsif ( $script eq "invert_align" )
 391     {
 392         @options = qw(
 393             soft|s
 394         );
 395     }
 396     elsif ( $script eq "patscan_seq" )
 397     {
 398         @options = qw(
 399             patterns|p=s
 400             patterns_in|P=s
 401             comp|c
 402             max_hits|h=s
 403             max_misses|m=s
 404             genome|g=s
 405         );
 406     }
 407     elsif ( $script eq "soap_seq" )
 408     {
 409         @options = qw(
 410             in_file|i=s
 411             genome|g=s
 412             seed_size|s=s
 413             mismatches|m=s
 414             gap_size|G=s
 415             cpus|c=s
 416         );
 417     }
 418     elsif ( $script eq "match_seq" )
 419     {
 420         @options = qw(
 421             word_size|w=s
 422             direction|d=s
 423         );
 424     }
 425     elsif ( $script eq "create_vmatch_index" )
 426     {
 427         @options = qw(
 428             index_name|i=s
 429             prefix_length|p=s
 430             no_stream|x
 431         );
 432     }
 433     elsif ( $script eq "write_bed" )
 434     {
 435         @options = qw(
 436             cols|c=s
 437             check|C
 438             no_stream|x
 439             data_out|o=s
 440             compress|Z
 441         );
 442     }
 443     elsif ( $script eq "write_psl" )
 444     {
 445         @options = qw(
 446             no_stream|x
 447             data_out|o=s
 448             compress|Z
 449         );
 450     }
 451     elsif ( $script eq "write_fixedstep" )
 452     {
 453         @options = qw(
 454             no_stream|x
 455             data_out|o=s
 456             compress|Z
 457         );
 458     }
 459     elsif ( $script eq "write_2bit" )
 460     {
 461         @options = qw(
 462             no_stream|x
 463             data_out|o=s
 464             no_mask|N
 465         );
 466     }
 467     elsif ( $script eq "write_solid" )
 468     {
 469         @options = qw(
 470             wrap|w=s
 471             no_stream|x
 472             data_out|o=s
 473             compress|Z
 474         );
 475     }
 476     elsif ( $script eq "write_ucsc_config" )
 477     {
 478         @options = qw(
 479             no_stream|x
 480             data_out|o=s
 481         );
 482     }
 483     elsif ( $script eq "plot_seqlogo" )
 484     {
 485         @options = qw(
 486             no_stream|x
 487             data_out|o=s
 488         );
 489     }
 490     elsif ( $script eq "plot_phastcons_profiles" )
 491     {
 492         @options = qw(
 493             no_stream|x
 494             data_out|o=s
 495             genome|g=s
 496             mean|m
 497             median|M
 498             flank|f=s
 499             terminal|t=s
 500             title|T=s
 501             xlabel|X=s
 502             ylabel|Y=s
 503         );
 504     }
 505     elsif ( $script eq "head_records" )
 506     {
 507         @options = qw(
 508             num|n=s
 509         );
 510     }
 511     elsif ( $script eq "remove_keys" )
 512     {
 513         @options = qw(
 514             keys|k=s
 515             save_keys|K=s
 516         );
 517     }
 518     elsif ( $script eq "remove_adaptor" )
 519     {
 520         @options = qw(
 521             adaptor|a=s
 522             mismatches|m=s
 523             remove|r=s
 524             offset|o=s
 525         );
 526     }
 527     elsif ( $script eq "remove_mysql_tables" )
 528     {
 529         @options = qw(
 530             database|d=s
 531             tables|t=s
 532             keys|k=s
 533             user|u=s
 534             password|p=s
 535             no_stream|x
 536         );
 537     }
 538     elsif ( $script eq "remove_ucsc_tracks" )
 539     {
 540         @options = qw(
 541             database|d=s
 542             tracks|t=s
 543             keys|k=s
 544             config_file|c=s
 545             user|u=s
 546             password|p=s
 547             no_stream|x
 548         );
 549     }
 550     elsif ( $script eq "rename_keys" )
 551     {
 552         @options = qw(
 553             keys|k=s
 554         );
 555     }
 556     elsif ( $script eq "uniq_vals" )
 557     {
 558         @options = qw(
 559             key|k=s
 560             invert|i
 561         );
 562     }
 563     elsif ( $script eq "merge_vals" )
 564     {
 565         @options = qw(
 566             keys|k=s
 567             delimit|d=s
 568         );
 569     }
 570     elsif ( $script eq "merge_records" )
 571     {
 572         @options = qw(
 573             keys|k=s
 574             merge|m=s
 575         );
 576     }
 577     elsif ( $script eq "sort_records" )
 578     {
 579         @options = qw(
 580             reverse|r
 581             keys|k=s
 582         );
 583     }
 584     elsif ( $script eq "plot_histogram" )
 585     {
 586         @options = qw(
 587             no_stream|x
 588             data_out|o=s
 589             terminal|t=s
 590             title|T=s
 591             xlabel|X=s
 592             ylabel|Y=s
 593             key|k=s
 594             sort|s=s
 595         );
 596     }
 597     elsif ( $script eq "plot_lendist" )
 598     {
 599         @options = qw(
 600             no_stream|x
 601             data_out|o=s
 602             terminal|t=s
 603             title|T=s
 604             xlabel|X=s
 605             ylabel|Y=s
 606             key|k=s
 607         );
 608     }
 609     elsif ( $script eq "plot_chrdist" )
 610     {
 611         @options = qw(
 612             no_stream|x
 613             data_out|o=s
 614             terminal|t=s
 615             title|T=s
 616             xlabel|X=s
 617             ylabel|Y=s
 618         );
 619     }
 620     elsif ( $script eq "plot_karyogram" )
 621     {
 622         @options = qw(
 623             no_stream|x
 624             data_out|o=s
 625             genome|g=s
 626             feat_color|f=s
 627         );
 628     }
 629     elsif ( $script eq "plot_matches" )
 630     {
 631         @options = qw(
 632             no_stream|x
 633             data_out|o=s
 634             terminal|t=s
 635             title|T=s
 636             xlabel|X=s
 637             ylabel|Y=s
 638             direction|d=s
 639         );
 640     }
 641     elsif ( $script eq "length_vals" )
 642     {
 643         @options = qw(
 644             keys|k=s
 645         );
 646     }
 647     elsif ( $script eq "sum_vals" )
 648     {
 649         @options = qw(
 650             no_stream|x
 651             data_out|o=s
 652             keys|k=s
 653         );
 654     }
 655     elsif ( $script eq "mean_vals" )
 656     {
 657         @options = qw(
 658             no_stream|x
 659             data_out|o=s
 660             keys|k=s
 661         );
 662     }
 663     elsif ( $script eq "median_vals" )
 664     {
 665         @options = qw(
 666             no_stream|x
 667             data_out|o=s
 668             keys|k=s
 669         );
 670     }
 671     elsif ( $script eq "max_vals" )
 672     {
 673         @options = qw(
 674             no_stream|x
 675             data_out|o=s
 676             keys|k=s
 677         );
 678     }
 679     elsif ( $script eq "min_vals" )
 680     {
 681         @options = qw(
 682             no_stream|x
 683             data_out|o=s
 684             keys|k=s
 685         );
 686     }
 687     elsif ( $script eq "upload_to_ucsc" )
 688     {
 689         @options = qw(
 690             no_stream|x
 691             database|d=s
 692             table|t=s
 693             short_label|s=s
 694             long_label|l=s
 695             group|g=s
 696             priority|p=f
 697             use_score|u
 698             visibility|V=s
 699             color|c=s
 700             check|C
 701         );
 702     }
 703
 704     push @options, qw(
 705         stream_in|I=s
 706         stream_out|O=s
 707         verbose|v
 708         help|?
 709     );
 710
 711 #    print STDERR Dumper( \@options );
 712
 713     GetOptions(
 714         \%options,
 715         @options,
 716     );
 717
 718 #    print STDERR Dumper( \%options );
 719
 720     if ( -t STDIN && scalar( keys %options ) == 0 or $options{ "help" } ) {
 721         return wantarray ? %options : \%options;
 722     }
 723
 724     $options{ "cols" }      = [ split ",", $options{ "cols" } ]      if defined $options{ "cols" };
 725     $options{ "keys" }      = [ split ",", $options{ "keys" } ]      if defined $options{ "keys" };
 726     $options{ "no_keys" }   = [ split ",", $options{ "no_keys" } ]   if defined $options{ "no_keys" };
 727     $options{ "save_keys" } = [ split ",", $options{ "save_keys" } ] if defined $options{ "save_keys" };
 728     $options{ "quals" }     = [ split ",", $options{ "quals" } ]     if defined $options{ "quals" };
 729     $options{ "feats" }     = [ split ",", $options{ "feats" } ]     if defined $options{ "feats" };
 730     $options{ "frames" }    = [ split ",", $options{ "frames" } ]    if defined $options{ "frames" };
 731     $options{ "samples" }   = [ split ",", $options{ "samples" } ]   if defined $options{ "samples" };
 732     $options{ "tables" }    = [ split ",", $options{ "tables" } ]    if defined $options{ "tables" };
 733     $options{ "tracks" }    = [ split ",", $options{ "tracks" } ]    if defined $options{ "tracks" };
 734
 735     # ---- check arguments ----
 736
 737     if ( $options{ 'data_in' } )
 738     {
 739         $options{ "files" } = Maasha::Biopieces::getopt_files( $options{ 'data_in' } );
 740
 741         Maasha::Common::error( qq(Argument to --data_in must be a valid file or fileglob expression) ) if scalar @{ $options{ "files" } } == 0;
 742     }
 743
 744     map { Maasha::Common::error( qq(Argument to --cols must be a whole numbers - not "$_") ) if $_ !~ /^\d+$/ } @{ $options{ "cols" } } if $options{ "cols" };
 745
 746     # print STDERR Dumper( \%options );
 747
 748     $real = "beg|end|word_size|wrap|tile_size|len|prefix_length|mismatches|offset|num|skip|cpus|window_size|step_size";
 749
 750     foreach $opt ( keys %options )
 751     {
 752         if ( $opt =~ /stream_in|pattern_in|exact_in/ and not -f $options{ $opt } )
 753         {
 754             Maasha::Common::error( qq(Argument to --$opt must be a valid file or fileglob expression - not "$options{ $opt }") );
 755         }
 756         elsif ( $opt =~ /$real/ and $options{ $opt } !~ /^\d+$/ )
 757         {
 758             Maasha::Common::error( qq(Argument to --$opt must be a whole number - not "$options{ $opt }") );
 759         }
 760         elsif ( $opt =~ /max_hits|max_hits|max_misses|dist|edit_dist|flank|gap|hamming_dist|priority/ and $options{ $opt } !~ /^-?\d+$/ )
 761         {
 762             Maasha::Common::error( qq(Argument to --$opt must be an integer - not "$options{ $opt }") );
 763         }
 764         elsif ( $opt =~ /identity|threshold/ and $options{ $opt } !~ /^-?(?:\d+(?:\.\d*)?|\.\d+)$/ )
 765         {
 766             Maasha::Common::error( qq(Argument to --$opt must be a decimal number - not "$options{ $opt }") );
 767         }
 768         elsif ( $opt =~ /e_val/ and $options{ $opt } !~ /^([+-]?)(?=\d|\.\d)\d*(\.\d*)?([Ee]([+-]?\d+))?$/ )
 769         {
 770             Maasha::Common::error( qq(Argument to --$opt must be a float - not "$options{ $opt }") );
 771         }
 772         elsif ( $opt =~ /strand/ and $options{ $opt } !~ /^(\+|-)$/ )
 773         {
 774             Maasha::Common::error( qq(Argument to --$opt must be "+" or "-" - not "$options{ $opt }") );
 775         }
 776         elsif ( $opt eq "genome" )
 777         {
 778             @genomes = Maasha::Common::ls_dirs( "$ENV{ 'BP_DATA' }/genomes" );
 779             map { $_ =~ s/.*\/(.+)$/$1/ } @genomes;
 780
 781             if ( not grep { $_ =~ /^$options{ $opt }$/ } @genomes ) {
 782                 Maasha::Common::error( qq(Genome $options{ $opt } not found in "$ENV{ 'BP_DATA' }/genomes/") );
 783             }
 784         }
 785         elsif ( $opt eq "terminal" and not $options{ $opt } =~ /^(svg|post|dumb|x11)/ )
 786         {
 787             Maasha::Common::error( qq(Bad --$opt argument "$options{ $opt }") );
 788         }
 789         elsif ( $opt eq "table" and $options{ $opt } =~ /(-|\.)/ )
 790         {
 791             Maasha::Common::error( qq(Character '$1' is not allowed in table name: $options{ $opt }) );
 792         }
 793         elsif ( $opt eq "merge" and $options{ $opt } !~ /^(AandB|AorB|BorA|AnotB|BnotA)$/ )
 794         {
 795             Maasha::Common::error( qq(Argument to --$opt must be AandB, AorB, BorA, AnotB, or BnotA - not "$options{ $opt }") );
 796         }
 797         elsif ( $opt eq "format" and $script eq "read_solexa" and $options{ $opt } !~ /octal|decimal/ )
 798         {
 799             Maasha::Common::error( qq(Argument to --$opt must be octal or decimal - not "$options{ $opt }") );
 800         }
 801         elsif ( $opt eq "remove" and $script eq "remove_adaptor" and $options{ $opt } !~ /before|after|skip/ )
 802         {
 803             Maasha::Common::error( qq(Argument to --$opt must be before, after, or skip - not "$options{ $opt }") );
 804         }
 805     }
 806
 807     Maasha::Common::error( qq(no --database specified) )                if $script eq "remove_ucsc_tracks"  and not $options{ "database" };
 808     Maasha::Common::error( qq(no --index_name specified) )              if $script =~ /create_vmatch_index/ and not $options{ "index_name" };
 809     Maasha::Common::error( qq(no --in_file or --genome specified) )     if $script eq "soap_seq" and not $options{ "genome" } and not $options{ "in_file" };
 810     Maasha::Common::error( qq(both --in_file and --genome specified) )  if $script eq "soap_seq" and $options{ "genome" } and $options{ "in_file" };
 811     Maasha::Common::error( qq(no --genome specified) )                  if $script =~ /get_genome_align|get_genome_phastcons|plot_phastcons_profiles|plot_karyogram/ and not $options{ "genome" };
 812     Maasha::Common::error( qq(no --key specified) )                     if $script =~ /plot_lendist|plot_histogram/ and not $options{ "key" };
 813     Maasha::Common::error( qq(no --keys speficied) )                    if $script =~ /sort_records|sum_vals|mean_vals|median_vals|length_vals/ and not $options{ "keys" };
 814
 815     if ( $script eq "upload_to_ucsc" )
 816     {
 817         Maasha::Common::error( qq(no --database specified) ) if not $options{ "database" };
 818         Maasha::Common::error( qq(no --table specified) )    if not $options{ "table" };
 819     }
 820
 821     return wantarray ? %options : \%options;
 822 }
 823
 824
 825 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> SCRIPTS  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 826
 827
 828 sub script_print_usage
 829 {
 830     # Martin A. Hansen, January 2008.
 831
 832     # Retrieves usage information from file and
 833     # prints this nicely formatted.
 834
 835     my ( $in,        # handle to in stream
 836          $out,       # handle to out stream
 837          $options,   # options hash
 838        ) = @_;
 839
 840     # Returns nothing.
 841
 842     my ( $file, $wiki, $lines );
 843
 844     if ( $options->{ 'data_in' } ) {
 845         $file = $options->{ 'data_in' };
 846     } else {
 847         $file = join "", $ENV{ 'BP_DIR' }, "/bp_usage/", $options->{ 'SCRIPT' }, ".wiki";
 848     }
 849
 850     $wiki = Maasha::Gwiki::gwiki_read( $file );
 851
 852     ( $wiki->[ 2 ], $wiki->[ 3 ], $wiki->[ 0 ], $wiki->[ 1 ] ) = ( $wiki->[ 0 ], $wiki->[ 1 ], $wiki->[ 2 ], $wiki->[ 3 ] );
 853
 854     if ( not $options->{ "help" } ) {
 855         @{ $wiki } = grep { $_->[ 0 ]->{ 'SECTION' } =~ /Biopiece|summary|Usage|Options|Help/ } @{ $wiki };
 856     }
 857
 858     $lines = Maasha::Gwiki::gwiki2ascii( $wiki );
 859
 860     print STDERR "$_\n" foreach @{ $lines };
 861
 862     exit;
 863 }
 864
 865
 866 sub script_read_psl
 867 {
 868     # Martin A. Hansen, August 2007.
 869
 870     # Read psl table from stream or file.
 871
 872     my ( $in,        # handle to in stream
 873          $out,       # handle to out stream
 874          $options,   # options hash
 875        ) = @_;
 876
 877     # Returns nothing.
 878
 879     my ( $record, $file, $data_in, $num );
 880
 881     while ( $record = Maasha::Biopieces::get_record( $in ) ) {
 882         Maasha::Biopieces::put_record( $record, $out );
 883     }
 884
 885     $num = 1;
 886
 887     foreach $file ( @{ $options->{ "files" } } )
 888     {
 889         $data_in = Maasha::Common::read_open( $file );
 890
 891         while ( $record = Maasha::UCSC::psl_get_entry( $data_in ) )
 892         {
 893             Maasha::Biopieces::put_record( $record, $out );
 894
 895             goto NUM if $options->{ "num" } and $num == $options->{ "num" };
 896
 897             $num++;
 898         }
 899     }
 900
 901     NUM:
 902 }
 903
 904
 905 sub script_read_bed
 906 {
 907     # Martin A. Hansen, August 2007.
 908
 909     # Read bed table from stream or file.
 910
 911     my ( $in,        # handle to in stream
 912          $out,       # handle to out stream
 913          $options,   # options hash
 914        ) = @_;
 915
 916     # Returns nothing.
 917
 918     my ( $cols, $file, $record, $bed_entry, $data_in, $num );
 919
 920     $cols = $options->{ 'cols' }->[ 0 ];
 921
 922     while ( $record = Maasha::Biopieces::get_record( $in ) ) {
 923         Maasha::Biopieces::put_record( $record, $out );
 924     }
 925
 926     $num = 1;
 927
 928     foreach $file ( @{ $options->{ "files" } } )
 929     {
 930         $data_in = Maasha::Common::read_open( $file );
 931
 932         while ( $bed_entry = Maasha::UCSC::BED::bed_entry_get( $data_in, $cols, $options->{ 'check' } ) )
 933         {
 934             $record = Maasha::UCSC::BED::bed2biopiece( $bed_entry );
 935
 936             Maasha::Biopieces::put_record( $record, $out );
 937
 938             goto NUM if $options->{ "num" } and $num == $options->{ "num" };
 939
 940             $num++;
 941         }
 942
 943         close $data_in;
 944     }
 945
 946     NUM:
 947
 948     close $data_in if $data_in;
 949 }
 950
 951
 952 sub script_read_fixedstep
 953 {
 954     # Martin A. Hansen, Juli 2008.
 955
 956     # Read fixedstep wiggle format from stream or file.
 957
 958     my ( $in,        # handle to in stream
 959          $out,       # handle to out stream
 960          $options,   # options hash
 961        ) = @_;
 962
 963     # Returns nothing.
 964
 965     my ( $file, $record, $entry, $data_in, $num );
 966
 967     while ( $record = Maasha::Biopieces::get_record( $in ) ) {
 968         Maasha::Biopieces::put_record( $record, $out );
 969     }
 970
 971     $num = 1;
 972
 973     foreach $file ( @{ $options->{ "files" } } )
 974     {
 975         $data_in = Maasha::Common::read_open( $file );
 976
 977         while ( $entry = Maasha::UCSC::Wiggle::fixedstep_entry_get( $data_in ) )
 978         {
 979             $record = Maasha::UCSC::Wiggle::fixedstep2biopiece( $entry );
 980
 981             Maasha::Biopieces::put_record( $record, $out );
 982
 983             goto NUM if $options->{ "num" } and $num == $options->{ "num" };
 984
 985             $num++;
 986         }
 987
 988         close $data_in;
 989     }
 990
 991     NUM:
 992
 993     close $data_in if $data_in;
 994 }
 995
 996
 997 sub script_read_blast_tab
 998 {
 999     # Martin A. Hansen, September 2007.
1000
1001     # Read tabular BLAST output from NCBI blast run with -m8 or -m9.
1002
1003     my ( $in,        # handle to in stream
1004          $out,       # handle to out stream
1005          $options,   # options hash
1006        ) = @_;
1007
1008     # Returns nothing.
1009
1010     my ( $file, $line, @fields, $strand, $record, $data_in, $num );
1011
1012     while ( $record = Maasha::Biopieces::get_record( $in ) ) {
1013         Maasha::Biopieces::put_record( $record, $out );
1014     }
1015
1016     $num = 1;
1017
1018     foreach $file ( @{ $options->{ "files" } } )
1019     {
1020         $data_in = Maasha::Common::read_open( $file );
1021
1022         while ( $line = <$data_in> )
1023         {
1024             chomp $line;
1025
1026             next if $line =~ /^#/;
1027
1028             @fields = split /\t/, $line;
1029
1030             $record->{ "REC_TYPE" }   = "BLAST";
1031             $record->{ "Q_ID" }       = $fields[ 0 ];
1032             $record->{ "S_ID" }       = $fields[ 1 ];
1033             $record->{ "IDENT" }      = $fields[ 2 ];
1034             $record->{ "ALIGN_LEN" }  = $fields[ 3 ];
1035             $record->{ "MISMATCHES" } = $fields[ 4 ];
1036             $record->{ "GAPS" }       = $fields[ 5 ];
1037             $record->{ "Q_BEG" }      = $fields[ 6 ] - 1; # BLAST is 1-based
1038             $record->{ "Q_END" }      = $fields[ 7 ] - 1; # BLAST is 1-based
1039             $record->{ "S_BEG" }      = $fields[ 8 ] - 1; # BLAST is 1-based
1040             $record->{ "S_END" }      = $fields[ 9 ] - 1; # BLAST is 1-based
1041             $record->{ "E_VAL" }      = $fields[ 10 ];
1042             $record->{ "BIT_SCORE" }  = $fields[ 11 ];
1043
1044             if ( $record->{ "S_BEG" } > $record->{ "S_END" } )
1045             {
1046                 $record->{ "STRAND" } = '-';
1047
1048                 ( $record->{ "S_BEG" }, $record->{ "S_END" } ) = ( $record->{ "S_END" }, $record->{ "S_BEG" } );
1049             }
1050             else
1051             {
1052                 $record->{ "STRAND" } = '+';
1053             }
1054
1055             Maasha::Biopieces::put_record( $record, $out );
1056
1057             goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1058
1059             $num++;
1060         }
1061
1062         close $data_in;
1063     }
1064
1065     NUM:
1066
1067     close $data_in if $data_in;
1068 }
1069
1070
1071 sub script_read_embl
1072 {
1073     # Martin A. Hansen, August 2007.
1074
1075     # Read EMBL format.
1076
1077     my ( $in,        # handle to in stream
1078          $out,       # handle to out stream
1079          $options,   # options hash
1080        ) = @_;
1081
1082     # Returns nothing.
1083
1084     my ( %options2, $file, $data_in, $num, $entry, $record );
1085
1086     map { $options2{ "keys" }{ $_ } = 1 }  @{ $options->{ "keys" } };
1087     map { $options2{ "feats" }{ $_ } = 1 } @{ $options->{ "feats" } };
1088     map { $options2{ "quals" }{ $_ } = 1 } @{ $options->{ "quals" } };
1089
1090     while ( $record = Maasha::Biopieces::get_record( $in ) ) {
1091         Maasha::Biopieces::put_record( $record, $out );
1092     }
1093
1094     $num = 1;
1095
1096     foreach $file ( @{ $options->{ "files" } } )
1097     {
1098         $data_in = Maasha::Common::read_open( $file );
1099
1100         while ( $entry = Maasha::EMBL::get_embl_entry( $data_in ) )
1101         {
1102             $record = Maasha::EMBL::parse_embl_entry( $entry, \%options2 );
1103
1104             my ( $feat, $feat2, $qual, $qual_val, $record_copy );
1105
1106             $record_copy = dclone $record;
1107
1108             delete $record_copy->{ "FT" };
1109
1110             Maasha::Biopieces::put_record( $record_copy, $out );
1111
1112             delete $record_copy->{ "SEQ" };
1113
1114             foreach $feat ( keys %{ $record->{ "FT" } } )
1115             {
1116                 $record_copy->{ "FEAT_TYPE" } = $feat;
1117
1118                 foreach $feat2 ( @{ $record->{ "FT" }->{ $feat } } )
1119                 {
1120                     foreach $qual ( keys %{ $feat2 } )
1121                     {
1122                         $qual_val = join "; ", @{ $feat2->{ $qual } };
1123
1124                         $qual =~ s/^_//;
1125                         $qual = uc $qual;
1126
1127                         $record_copy->{ $qual } = $qual_val;
1128                     }
1129
1130                     Maasha::Biopieces::put_record( $record_copy, $out );
1131                 }
1132             }
1133
1134             goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1135
1136             $num++;
1137         }
1138
1139         close $data_in;
1140     }
1141
1142     NUM:
1143
1144     close $data_in if $data_in;
1145 }
1146
1147
1148 sub script_read_stockholm
1149 {
1150     # Martin A. Hansen, August 2007.
1151
1152     # Read Stockholm format.
1153
1154     my ( $in,        # handle to in stream
1155          $out,       # handle to out stream
1156          $options,   # options hash
1157        ) = @_;
1158
1159     # Returns nothing.
1160
1161     my ( $data_in, $file, $num, $entry, $record, $record_anno, $record_align, $key, $seq );
1162
1163     while ( $record = Maasha::Biopieces::get_record( $in ) ) {
1164         Maasha::Biopieces::put_record( $record, $out );
1165     }
1166
1167     $num = 1;
1168
1169     foreach $file ( @{ $options->{ "files" } } )
1170     {
1171         $data_in = Maasha::Common::read_open( $file );
1172
1173         while ( $entry = Maasha::Stockholm::get_stockholm_entry( $data_in ) )
1174         {
1175             $record = Maasha::Stockholm::parse_stockholm_entry( $entry );
1176
1177             undef $record_anno;
1178
1179             foreach $key ( keys %{ $record->{ "GF" } } ) {
1180                 $record_anno->{ $key } = $record->{ "GF" }->{ $key };
1181             }
1182
1183             $record_anno->{ "ALIGN" } = $num;
1184
1185             Maasha::Biopieces::put_record( $record_anno, $out );
1186
1187             foreach $seq ( @{ $record->{ "ALIGN" } } )
1188             {
1189                 undef $record_align;
1190
1191                 $record_align = {
1192                     SEQ_NAME  => $seq->[ 0 ],
1193                     SEQ       => $seq->[ 1 ],
1194                 };
1195
1196                 Maasha::Biopieces::put_record( $record_align, $out );
1197             }
1198
1199             goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1200
1201             $num++;
1202         }
1203
1204         close $data_in;
1205     }
1206
1207     NUM:
1208
1209     close $data_in if $data_in;
1210 }
1211
1212
1213 sub script_read_phastcons
1214 {
1215     # Martin A. Hansen, December 2007.
1216
1217     # Read PhastCons format.
1218
1219     my ( $in,        # handle to in stream
1220          $out,       # handle to out stream
1221          $options,   # options hash
1222        ) = @_;
1223
1224     # Returns nothing.
1225
1226     my ( $data_in, $file, $num, $entry, @records, $record );
1227
1228     $options->{ "min" }       ||= 10;
1229     $options->{ "dist" }      ||= 25;
1230     $options->{ "threshold" } ||= 0.8;
1231     $options->{ "gap" }       ||= 5;
1232
1233     while ( $record = Maasha::Biopieces::get_record( $in ) ) {
1234         Maasha::Biopieces::put_record( $record, $out );
1235     }
1236
1237     $num = 1;
1238
1239     foreach $file ( @{ $options->{ "files" } } )
1240     {
1241         $data_in = Maasha::Common::read_open( $file );
1242
1243         while ( $entry = Maasha::UCSC::fixedstep_get_entry( $data_in ) )
1244         {
1245             @records = Maasha::UCSC::phastcons_parse_entry( $entry, $options );
1246
1247             foreach $record ( @records )
1248             {
1249                 $record->{ "REC_TYPE" } = "BED";
1250                 $record->{ "BED_LEN" }  = $record->{ "CHR_END" } - $record->{ "CHR_BEG" } + 1;
1251
1252                 Maasha::Biopieces::put_record( $record, $out );
1253
1254                 goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1255
1256                 $num++;
1257             }
1258         }
1259
1260         close $data_in;
1261     }
1262
1263     NUM:
1264
1265     close $data_in if $data_in;
1266 }
1267
1268
1269 sub script_read_soft
1270 {
1271     # Martin A. Hansen, December 2007.
1272
1273     # Read soft format.
1274     # http://www.ncbi.nlm.nih.gov/geo/info/soft2.html
1275
1276     my ( $in,        # handle to in stream
1277          $out,       # handle to out stream
1278          $options,   # options hash
1279        ) = @_;
1280
1281     # Returns nothing.
1282
1283     my ( $data_in, $file, $num, $records, $record, $soft_index, $fh, @platforms, $plat_table, @samples, $sample, $old_end, $skip );
1284
1285     while ( $record = Maasha::Biopieces::get_record( $in ) ) {
1286         Maasha::Biopieces::put_record( $record, $out );
1287     }
1288
1289     $num = 1;
1290
1291     foreach $file ( @{ $options->{ "files" } } )
1292     {
1293         print STDERR "Creating index for file: $file\n" if $options->{ "verbose" };
1294
1295         $soft_index = Maasha::NCBI::soft_index_file( $file );
1296
1297         $fh         = Maasha::Common::read_open( $file );
1298
1299         @platforms  = grep { $_->{ "SECTION" } =~ /PLATFORM/ } @{ $soft_index };
1300
1301         print STDERR "Getting platform tables for file: $file\n" if $options->{ "verbose" };
1302
1303         $plat_table = Maasha::NCBI::soft_get_platform( $fh, $platforms[ 0 ]->{ "LINE_BEG" }, $platforms[ -1 ]->{ "LINE_END" } );
1304
1305         @samples    = grep { $_->{ "SECTION" } =~ /SAMPLE/ } @{ $soft_index };
1306
1307         $old_end    = $platforms[ -1 ]->{ "LINE_END" };
1308
1309         foreach $sample ( @samples )
1310         {
1311             $skip = 0;
1312             $skip = 1 if ( $options->{ "samples" } and grep { $sample->{ "SECTION" } !~ /$_/ } @{ $options->{ "samples" } } );
1313
1314             print STDERR "Getting samples for dataset: $sample->{ 'SECTION' }\n" if $options->{ "verbose" } and not $skip;
1315
1316             $records = Maasha::NCBI::soft_get_sample( $fh, $plat_table, $sample->{ "LINE_BEG" } - $old_end - 1, $sample->{ "LINE_END" } - $old_end - 1, $skip );
1317
1318             foreach $record ( @{ $records } )
1319             {
1320                 Maasha::Biopieces::put_record( $record, $out );
1321
1322                 goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1323
1324                 $num++;
1325             }
1326
1327             $old_end = $sample->{ "LINE_END" };
1328         }
1329
1330         close $fh;
1331     }
1332
1333     NUM:
1334
1335     close $data_in if $data_in;
1336     close $fh if $fh;
1337 }
1338
1339
1340 sub script_read_gff
1341 {
1342     # Martin A. Hansen, February 2008.
1343
1344     # Read soft format.
1345     # http://www.ncbi.nlm.nih.gov/geo/info/soft2.html
1346
1347     my ( $in,        # handle to in stream
1348          $out,       # handle to out stream
1349          $options,   # options hash
1350        ) = @_;
1351
1352     # Returns nothing.
1353
1354     my ( $data_in, $file, $fh, $num, $record, $entry );
1355
1356     while ( $record = Maasha::Biopieces::get_record( $in ) ) {
1357         Maasha::Biopieces::put_record( $record, $out );
1358     }
1359
1360     $num = 1;
1361
1362     foreach $file ( @{ $options->{ "files" } } )
1363     {
1364         $fh = Maasha::Common::read_open( $file );
1365
1366         while ( $entry = Maasha::GFF::get_entry( $fh ) )
1367         {
1368             Maasha::Biopieces::put_record( $entry, $out );
1369
1370             goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1371
1372             $num++;
1373         }
1374
1375         close $fh;
1376     }
1377
1378     NUM:
1379
1380     close $data_in if $data_in;
1381 }
1382
1383
1384 sub script_read_2bit
1385 {
1386     # Martin A. Hansen, March 2008.
1387
1388     # Read sequences from 2bit file.
1389
1390     my ( $in,        # handle to in stream
1391          $out,       # handle to out stream
1392          $options,   # options hash
1393        ) = @_;
1394
1395     # Returns nothing.
1396
1397     my ( $record, $file, $data_in, $mask, $toc, $line, $num );
1398
1399     $mask = 1 if not $options->{ "no_mask" };
1400
1401     while ( $record = Maasha::Biopieces::get_record( $in ) ) {
1402         Maasha::Biopieces::put_record( $record, $out );
1403     }
1404
1405     $num = 1;
1406
1407     foreach $file ( @{ $options->{ "files" } } )
1408     {
1409         $data_in = Maasha::Common::read_open( $file );
1410
1411         $toc = Maasha::TwoBit::twobit_get_TOC( $data_in );
1412
1413         foreach $line ( @{ $toc } )
1414         {
1415             $record->{ "SEQ_NAME" } = $line->[ 0 ];
1416             $record->{ "SEQ" }      = Maasha::TwoBit::twobit_get_seq( $data_in, $line->[ 1 ], undef, undef, $mask );
1417             $record->{ "SEQ_LEN" }  = length $record->{ "SEQ" };
1418
1419             Maasha::Biopieces::put_record( $record, $out );
1420
1421             goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1422
1423             $num++;
1424         }
1425
1426         close $data_in;
1427     }
1428
1429     NUM:
1430
1431     close $data_in if $data_in;
1432 }
1433
1434
1435 sub script_read_solexa
1436 {
1437     # Martin A. Hansen, March 2008.
1438
1439     # Read Solexa sequence reads from file.
1440
1441     my ( $in,        # handle to in stream
1442          $out,       # handle to out stream
1443          $options,   # options hash
1444        ) = @_;
1445
1446     # Returns nothing.
1447
1448     my ( $record, $file, $data_in, $entry, $num, @seqs, @scores, $i );
1449
1450     $options->{ "format" }  ||= "octal";
1451     $options->{ "quality" } ||= 20;
1452
1453     while ( $record = Maasha::Biopieces::get_record( $in ) ) {
1454         Maasha::Biopieces::put_record( $record, $out );
1455     }
1456
1457     $num = 1;
1458
1459     foreach $file ( @{ $options->{ "files" } } )
1460     {
1461         $data_in = Maasha::Common::read_open( $file );
1462
1463         if ( $options->{ "format" } eq "octal" )
1464         {
1465             while ( $entry = Maasha::Solexa::solexa_get_entry_octal( $data_in ) )
1466             {
1467                 $record = Maasha::Solexa::solexa2biopiece( $entry, $options->{ "quality" } );
1468
1469                 Maasha::Biopieces::put_record( $record, $out );
1470
1471                 goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1472
1473                 $num++;
1474             }
1475         }
1476         else
1477         {
1478             while ( $entry = Maasha::Solexa::solexa_get_entry_decimal( $data_in ) )
1479             {
1480                 $record = Maasha::Solexa::solexa2biopiece( $entry, $options->{ "quality" } );
1481
1482                 Maasha::Biopieces::put_record( $record, $out );
1483
1484                 goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1485
1486                 $num++;
1487             }
1488         }
1489
1490         close $data_in;
1491     }
1492
1493     NUM:
1494
1495     close $data_in if $data_in;
1496 }
1497
1498
1499 sub script_read_solid
1500 {
1501     # Martin A. Hansen, April 2008.
1502
1503     # Read Solid sequence from file.
1504
1505     my ( $in,        # handle to in stream
1506          $out,       # handle to out stream
1507          $options,   # options hash
1508        ) = @_;
1509
1510     # Returns nothing.
1511
1512     my ( $record, $file, $data_in, $line, $num, $seq_name, $seq_cs, $seq_qual, @scores, @seqs, $i );
1513
1514     $options->{ "quality" } ||= 15;
1515
1516     while ( $record = Maasha::Biopieces::get_record( $in ) ) {
1517         Maasha::Biopieces::put_record( $record, $out );
1518     }
1519
1520     $num = 1;
1521
1522     foreach $file ( @{ $options->{ "files" } } )
1523     {
1524         $data_in = Maasha::Common::read_open( $file );
1525
1526         while ( $line = <$data_in> )
1527         {
1528             chomp $line;
1529
1530             ( $seq_name, $seq_cs, $seq_qual ) = split /\t/, $line;
1531
1532             @scores = split /,/, $seq_qual;
1533             @seqs   = split //, Maasha::Solid::color_space2seq( $seq_cs );
1534
1535             for ( $i = 0; $i < @seqs; $i++ ) {
1536                 $seqs[ $i ] = lc $seqs[ $i ] if $scores[ $i ] < $options->{ "quality" };
1537             }
1538
1539             $record = {
1540                 REC_TYPE   => 'SOLID',
1541                 SEQ_NAME   => $seq_name,
1542                 SEQ_CS     => $seq_cs,
1543                 SEQ_QUAL   => join( ";", @scores ),
1544                 SEQ_LEN    => length $seq_cs,
1545                 SEQ        => join( "", @seqs ),
1546                 SCORE_MEAN => sprintf( "%.2f", Maasha::Calc::mean( \@scores ) ),
1547             };
1548
1549             Maasha::Biopieces::put_record( $record, $out );
1550
1551             goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1552
1553             $num++;
1554         }
1555
1556         close $data_in;
1557     }
1558
1559     NUM:
1560
1561     close $data_in if $data_in;
1562 }
1563
1564
1565 sub script_read_mysql
1566 {
1567     # Martin A. Hansen, May 2008.
1568
1569     # Read a MySQL query into stream.
1570
1571     my ( $in,        # handle to in stream
1572          $out,       # handle to out stream
1573          $options,   # options hash
1574        ) = @_;
1575
1576     # Returns nothing.
1577
1578     my ( $record, $dbh, $results );
1579
1580     $options->{ "user" }     ||= Maasha::UCSC::ucsc_get_user();
1581     $options->{ "password" } ||= Maasha::UCSC::ucsc_get_password();
1582
1583     while ( $record = Maasha::Biopieces::get_record( $in ) ) {
1584         Maasha::Biopieces::put_record( $record, $out );
1585     }
1586
1587     $dbh = Maasha::SQL::connect( $options->{ "database" }, $options->{ "user" }, $options->{ "password" } );
1588
1589     $results = Maasha::SQL::query_hashref_list( $dbh, $options->{ "query" } );
1590
1591     Maasha::SQL::disconnect( $dbh );
1592
1593     map { Maasha::Biopieces::put_record( $_ ) } @{ $results };
1594 }
1595
1596
1597 sub script_read_ucsc_config
1598 {
1599     # Martin A. Hansen, November 2008.
1600
1601     # Read track entries from UCSC Genome Browser '.ra' files.
1602
1603     my ( $in,        # handle to in stream
1604          $out,       # handle to out stream
1605          $options,   # options hash
1606        ) = @_;
1607
1608     # Returns nothing.
1609
1610     my ( $record, $file, $data_in, $entry, $num );
1611
1612     while ( $record = Maasha::Biopieces::get_record( $in ) ) {
1613         Maasha::Biopieces::put_record( $record, $out );
1614     }
1615
1616     $num = 1;
1617
1618     foreach $file ( @{ $options->{ "files" } } )
1619     {
1620         $data_in = Maasha::Common::read_open( $file );
1621
1622         while ( $record = Maasha::UCSC::ucsc_config_entry_get( $data_in ) )
1623         {
1624             $record->{ 'REC_TYPE' } = "UCSC Config";
1625
1626             Maasha::Biopieces::put_record( $record, $out );
1627
1628             goto NUM if $options->{ "num" } and $num == $options->{ "num" };
1629
1630             $num++;
1631         }
1632
1633         close $data_in;
1634     }
1635
1636     NUM:
1637
1638     close $data_in if $data_in;
1639 }
1640
1641
1642 sub script_length_seq
1643 {
1644     # Martin A. Hansen, August 2007.
1645
1646     # Determine the length of sequences in stream.
1647
1648     my ( $in,        # handle to in stream
1649          $out,       # handle to out stream
1650          $options,   # options hash
1651        ) = @_;
1652
1653     # Returns nothing.
1654
1655     my ( $record, $total );
1656
1657     while ( $record = Maasha::Biopieces::get_record( $in ) )
1658     {
1659         if ( $record->{ "SEQ" } )
1660         {
1661             $record->{ "SEQ_LEN" } = length $record->{ "SEQ" };
1662             $total += $record->{ "SEQ_LEN" };
1663         }
1664
1665         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
1666     }
1667
1668     Maasha::Biopieces::put_record( { TOTAL_SEQ_LEN => $total }, $out );
1669 }
1670
1671
1672 sub script_complexity_seq
1673 {
1674     # Martin A. Hansen, May 2008.
1675
1676     # Generates an index calculated as the most common di-residue over
1677     # the sequence length for all sequences in stream.
1678
1679     my ( $in,     # handle to in stream
1680          $out,    # handle to out stream
1681        ) = @_;
1682
1683     # Returns nothing.
1684
1685     my ( $record, $index );
1686
1687     while ( $record = Maasha::Biopieces::get_record( $in ) )
1688     {
1689         $record->{ "SEQ_COMPLEXITY" } = sprintf( "%.2f", Maasha::Seq::seq_complexity( $record->{ "SEQ" } ) ) if $record->{ "SEQ" };
1690
1691         Maasha::Biopieces::put_record( $record, $out );
1692     }
1693 }
1694
1695
1696 sub script_oligo_freq
1697 {
1698     # Martin A. Hansen, August 2007.
1699
1700     # Determine the length of sequences in stream.
1701
1702     my ( $in,        # handle to in stream
1703          $out,       # handle to out stream
1704          $options,   # options hash
1705        ) = @_;
1706
1707     # Returns nothing.
1708
1709     my ( $record, %oligos, @freq_table );
1710
1711     $options->{ "word_size" } ||= 7;
1712
1713     while ( $record = Maasha::Biopieces::get_record( $in ) )
1714     {
1715         if ( $record->{ "SEQ" } )
1716         {
1717             map { $oligos{ $_ }++ } Maasha::Seq::seq2oligos( \$record->{ "SEQ" }, $options->{ "word_size" } );
1718
1719             if ( not $options->{ "all" } )
1720             {
1721                 @freq_table = Maasha::Seq::oligo_freq( \%oligos );
1722
1723                 map { Maasha::Biopieces::put_record( $_, $out ) } @freq_table;
1724
1725                 undef %oligos;
1726             }
1727         }
1728
1729         Maasha::Biopieces::put_record( $record, $out );
1730     }
1731
1732     if ( $options->{ "all" } )
1733     {
1734         @freq_table = Maasha::Seq::oligo_freq( \%oligos );
1735
1736         map { Maasha::Biopieces::put_record( $_, $out ) } @freq_table;
1737     }
1738 }
1739
1740
1741 sub script_remove_indels
1742 {
1743     # Martin A. Hansen, August 2007.
1744
1745     # Remove indels from sequences in stream.
1746
1747     my ( $in,     # handle to in stream
1748          $out,    # handle to out stream
1749        ) = @_;
1750
1751     # Returns nothing.
1752
1753     my ( $record );
1754
1755     while ( $record = Maasha::Biopieces::get_record( $in ) )
1756     {
1757         $record->{ 'SEQ' } =~ tr/-~.//d if $record->{ "SEQ" };
1758
1759         Maasha::Biopieces::put_record( $record, $out );
1760     }
1761 }
1762
1763
1764 sub script_transliterate_seq
1765 {
1766     # Martin A. Hansen, August 2007.
1767
1768     # Transliterate chars from sequence in record.
1769
1770     my ( $in,        # handle to in stream
1771          $out,       # handle to out stream
1772          $options,   # options hash
1773        ) = @_;
1774
1775     # Returns nothing.
1776
1777     my ( $record, $search, $replace, $delete );
1778
1779     $search  = $options->{ "search" }  || "";
1780     $replace = $options->{ "replace" } || "";
1781     $delete  = $options->{ "delete" }  || "";
1782
1783     while ( $record = Maasha::Biopieces::get_record( $in ) )
1784     {
1785         if ( $record->{ "SEQ" } )
1786         {
1787             if ( $search and $replace ) {
1788                 eval "\$record->{ 'SEQ' } =~ tr/$search/$replace/";
1789             } elsif ( $delete ) {
1790                 eval "\$record->{ 'SEQ' } =~ tr/$delete//d";
1791             }
1792         }
1793
1794         Maasha::Biopieces::put_record( $record, $out );
1795     }
1796 }
1797
1798
1799 sub script_transliterate_vals
1800 {
1801     # Martin A. Hansen, April 2008.
1802
1803     # Transliterate chars from values in record.
1804
1805     my ( $in,        # handle to in stream
1806          $out,       # handle to out stream
1807          $options,   # options hash
1808        ) = @_;
1809
1810     # Returns nothing.
1811
1812     my ( $record, $search, $replace, $delete, $key );
1813
1814     $search  = $options->{ "search" }  || "";
1815     $replace = $options->{ "replace" } || "";
1816     $delete  = $options->{ "delete" }  || "";
1817
1818     while ( $record = Maasha::Biopieces::get_record( $in ) )
1819     {
1820         foreach $key ( @{ $options->{ "keys" } } )
1821         {
1822             if ( exists $record->{ $key } )
1823             {
1824                 if ( $search and $replace ) {
1825                     eval "\$record->{ $key } =~ tr/$search/$replace/";
1826                 } elsif ( $delete ) {
1827                     eval "\$record->{ $key } =~ tr/$delete//d";
1828                 }
1829             }
1830         }
1831
1832         Maasha::Biopieces::put_record( $record, $out );
1833     }
1834 }
1835
1836
1837 sub script_translate_seq
1838 {
1839     # Martin A. Hansen, February 2008.
1840
1841     # Translate DNA sequence into protein sequence.
1842
1843     my ( $in,        # handle to in stream
1844          $out,       # handle to out stream
1845          $options,   # options hash
1846        ) = @_;
1847
1848     # Returns nothing.
1849
1850     my ( $record, $frame, %new_record );
1851
1852     $options->{ "frames" } ||= [ 1, 2, 3, -1, -2, -3 ];
1853
1854     while ( $record = Maasha::Biopieces::get_record( $in ) )
1855     {
1856         if ( $record->{ "SEQ" } )
1857         {
1858             if ( Maasha::Seq::seq_guess_type( $record->{ "SEQ" } ) eq "dna" )
1859             {
1860                 foreach $frame ( @{ $options->{ "frames" } } )
1861                 {
1862                     %new_record = %{ $record };
1863
1864                     $new_record{ "SEQ" }     = Maasha::Seq::translate( $record->{ "SEQ" }, $frame );
1865                     $new_record{ "SEQ_LEN" } = length $new_record{ "SEQ" };
1866                     $new_record{ "FRAME" }   = $frame;
1867
1868                     Maasha::Biopieces::put_record( \%new_record, $out );
1869                 }
1870             }
1871         }
1872         else
1873         {
1874             Maasha::Biopieces::put_record( $record, $out );
1875         }
1876     }
1877 }
1878
1879
1880 sub script_get_genome_align
1881 {
1882     # Martin A. Hansen, April 2008.
1883
1884     # Gets a subalignment from a multiple genome alignment.
1885
1886     my ( $in,        # handle to in stream
1887          $out,       # handle to out stream
1888          $options,   # options hash
1889        ) = @_;
1890
1891     # Returns nothing.
1892
1893     my ( $record, $maf_track, $align, $align_num, $beg, $end, $len, $entry );
1894
1895     $options->{ "strand" } ||= "+";
1896
1897     $align_num = 1;
1898
1899     $maf_track = Maasha::Config::maf_track( $options->{ "genome" } );
1900
1901     if ( $options->{ "chr" } and $options->{ "beg" } and ( $options->{ "end" } or $options->{ "len" } ) )
1902     {
1903         $beg = $options->{ "beg" } - 1;
1904
1905         if ( $options->{ "end" } ) {
1906             $end = $options->{ "end" };
1907         } elsif ( $options->{ "len" } ) {
1908             $end = $beg + $options->{ "len" };
1909         }
1910
1911         $align = Maasha::UCSC::maf_extract( $BP_TMP, $options->{ "genome" }, $maf_track, $options->{ "chr" }, $beg, $end, $options->{ "strand" } );
1912
1913         foreach $entry ( @{ $align } )
1914         {
1915             $entry->{ "CHR" }     = $record->{ "CHR" };
1916             $entry->{ "CHR_BEG" } = $record->{ "CHR_BEG" };
1917             $entry->{ "CHR_END" } = $record->{ "CHR_END" };
1918             $entry->{ "STRAND" }  = $record->{ "STRAND" } || '+';
1919             $entry->{ "Q_ID" }    = $record->{ "Q_ID" };
1920             $entry->{ "SCORE" }   = $record->{ "SCORE" };
1921
1922             Maasha::Biopieces::put_record( $entry, $out );
1923         }
1924     }
1925
1926     while ( $record = Maasha::Biopieces::get_record( $in ) )
1927     {
1928         if ( $record->{ "REC_TYPE" } eq "BED" )
1929         {
1930             $align = Maasha::UCSC::maf_extract( $BP_TMP, $options->{ "genome" }, $maf_track, $record->{ "CHR" }, $record->{ "CHR_BEG" }, $record->{ "CHR_END" }, $record->{ "STRAND" } );
1931         }
1932         elsif ( $record->{ "REC_TYPE" } eq "VMATCH" )
1933         {
1934             $align = Maasha::UCSC::maf_extract( $BP_TMP, $options->{ "genome" }, $maf_track, $record->{ "S_ID" }, $record->{ "S_BEG" }, $record->{ "S_END" } + 1, $record->{ "STRAND" } );
1935         }
1936         elsif ( $record->{ "REC_TYPE" } eq "PSL" )
1937         {
1938             $align = Maasha::UCSC::maf_extract( $BP_TMP, $options->{ "genome" }, $maf_track, $record->{ "S_ID" }, $record->{ "S_BEG" }, $record->{ "S_END" }, $record->{ "STRAND" } );
1939         }
1940         elsif ( $record->{ "REC_TYPE" } eq "BLAST" )
1941         {
1942             $align = Maasha::UCSC::maf_extract( $BP_TMP, $options->{ "genome" }, $maf_track, $record->{ "S_ID" }, $record->{ "S_BEG" }, $record->{ "S_END" }, $record->{ "STRAND" } );
1943         }
1944
1945         foreach $entry ( @{ $align } )
1946         {
1947             $entry->{ "CHR" }     = $record->{ "CHR" };
1948             $entry->{ "CHR_BEG" } = $record->{ "CHR_BEG" };
1949             $entry->{ "CHR_END" } = $record->{ "CHR_END" };
1950             $entry->{ "STRAND" }  = $record->{ "STRAND" };
1951             $entry->{ "Q_ID" }    = $record->{ "Q_ID" };
1952             $entry->{ "SCORE" }   = $record->{ "SCORE" };
1953
1954             Maasha::Biopieces::put_record( $entry, $out );
1955         }
1956
1957         $align_num++;
1958     }
1959 }
1960
1961
1962 sub script_get_genome_phastcons
1963 {
1964     # Martin A. Hansen, February 2008.
1965
1966     # Get phastcons scores from genome intervals.
1967
1968     my ( $in,        # handle to in stream
1969          $out,       # handle to out stream
1970          $options,   # options hash
1971        ) = @_;
1972
1973     # Returns nothing.
1974
1975     my ( $phastcons_file, $phastcons_index, $index, $fh_phastcons, $scores, $record );
1976
1977     $options->{ "flank" } ||= 0;
1978
1979     $phastcons_file  = Maasha::Config::genome_phastcons( $options->{ "genome" } );
1980     $phastcons_index = Maasha::Config::genome_phastcons_index( $options->{ "genome" } );
1981
1982     $index           = Maasha::UCSC::fixedstep_index_retrieve( $phastcons_index );
1983     $fh_phastcons    = Maasha::Common::read_open( $phastcons_file );
1984
1985     if ( defined $options->{ "chr" } and defined $options->{ "beg" } and ( defined $options->{ "end" } or defined $options->{ "len" } ) )
1986     {
1987         $options->{ "beg" } -= 1;   # request is 1-based
1988         $options->{ "end" } -= 1;   # request is 1-based
1989
1990         if ( $options->{ "len" } ) {
1991             $options->{ "end" } = $options->{ "beg" } + $options->{ "len" } - 1;
1992         }
1993
1994         $scores = Maasha::UCSC::fixedstep_index_lookup( $index, $fh_phastcons, $options->{ "chr" }, $options->{ "beg" }, $options->{ "end" }, $options->{ "flank" } );
1995
1996         $record->{ "CHR" }       = $options->{ "chr" };
1997         $record->{ "CHR_BEG" }   = $options->{ "beg" } - $options->{ "flank" };
1998         $record->{ "CHR_END" }   = $options->{ "end" } + $options->{ "flank" };
1999
2000         $record->{ "PHASTCONS" }   = join ",", @{ $scores };
2001         $record->{ "PHAST_COUNT" } = scalar @{ $scores };  # DEBUG
2002
2003         Maasha::Biopieces::put_record( $record, $out );
2004     }
2005
2006     while ( $record = Maasha::Biopieces::get_record( $in ) )
2007     {
2008         if ( $record->{ "REC_TYPE" } eq "BED" )
2009         {
2010             $scores = Maasha::UCSC::fixedstep_index_lookup( $index, $fh_phastcons, $record->{ "CHR" }, $record->{ "CHR_BEG" }, $record->{ "CHR_END" }, $options->{ "flank" } );
2011         }
2012         elsif ( $record->{ "REC_TYPE" } eq "PSL" )
2013         {
2014             $scores = Maasha::UCSC::fixedstep_index_lookup( $index, $fh_phastcons, $record->{ "S_ID" }, $record->{ "S_BEG" }, $record->{ "S_END" }, $options->{ "flank" } );
2015         }
2016         elsif ( $record->{ "REC_TYPE" } eq "BLAST" )
2017         {
2018             $scores = Maasha::UCSC::fixedstep_index_lookup( $index, $fh_phastcons, $record->{ "S_ID" }, $record->{ "S_BEG" }, $record->{ "S_END" }, $options->{ "flank" } );
2019         }
2020
2021         $record->{ "PHASTCONS" } = join ",", @{ $scores } if @{ $scores };
2022 #        $record->{ "PHAST_COUNT" } = @{ $scores } if @{ $scores };  # DEBUG
2023
2024         Maasha::Biopieces::put_record( $record, $out );
2025     }
2026
2027     close $fh_phastcons if $fh_phastcons;
2028 }
2029
2030
2031 sub script_fold_seq
2032 {
2033     # Martin A. Hansen, December 2007.
2034
2035     # Folds sequences in stream into secondary structures.
2036
2037     my ( $in,     # handle to in stream
2038          $out,    # handle to out stream
2039        ) = @_;
2040
2041     # Returns nothing.
2042
2043     my ( $record, $type, $struct, $index );
2044
2045     while ( $record = Maasha::Biopieces::get_record( $in ) )
2046     {
2047         if ( $record->{ "SEQ" } )
2048         {
2049             if ( not $type ) {
2050                 $type = Maasha::Seq::seq_guess_type( $record->{ "SEQ" } );
2051             }
2052
2053             if ( $type ne "protein" )
2054             {
2055                 ( $struct, $index ) = Maasha::Seq::fold_struct_rnafold( $record->{ "SEQ" } );
2056                 $record->{ "SEC_STRUCT" }  = $struct;
2057                 $record->{ "FREE_ENERGY" } = $index;
2058                 $record->{ "SCORE" }       = abs int $index;
2059                 $record->{ "SIZE" }        = length $struct;
2060                 $record->{ "CONF" }        = "1," x $record->{ "SIZE" };
2061             }
2062         }
2063
2064         Maasha::Biopieces::put_record( $record, $out );
2065     }
2066 }
2067
2068
2069 sub script_tile_seq
2070 {
2071     # Martin A. Hansen, February 2008.
2072
2073     # Using the first sequence in stream as reference, tile
2074     # all subsequent sequences based on pairwise alignments.
2075
2076     my ( $in,        # handle to in stream
2077          $out,       # handle to out stream
2078          $options,   # options hash
2079        ) = @_;
2080
2081     # Returns nothing.
2082
2083     my ( $record, $first, $ref_entry, @entries );
2084
2085     $first = 1;
2086
2087     while ( $record = Maasha::Biopieces::get_record( $in ) )
2088     {
2089         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } )
2090         {
2091             if ( $first )
2092             {
2093                 $ref_entry = [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ];
2094
2095                 $first = 0;
2096             }
2097             else
2098             {
2099                 push @entries, [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ];
2100             }
2101         }
2102         else
2103         {
2104             Maasha::Biopieces::put_record( $record, $out );
2105         }
2106     }
2107
2108     @entries = Maasha::Align::align_tile( $ref_entry, \@entries, $options );
2109
2110     map { Maasha::Biopieces::put_record( { SEQ_NAME => $_->[ SEQ_NAME ], SEQ => $_->[ SEQ ] }, $out ) } @entries;
2111 }
2112
2113
2114 sub script_invert_align
2115 {
2116     # Martin A. Hansen, February 2008.
2117
2118     # Inverts an alignment showing only non-mathing residues
2119     # using the first sequence as reference.
2120
2121     my ( $in,        # handle to in stream
2122          $out,       # handle to out stream
2123          $options,   # options hash
2124        ) = @_;
2125
2126     # Returns nothing.
2127
2128     my ( $record, @entries );
2129
2130     while ( $record = Maasha::Biopieces::get_record( $in ) )
2131     {
2132         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } )
2133         {
2134             push @entries, [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ];
2135         }
2136         else
2137         {
2138             Maasha::Biopieces::put_record( $record, $out );
2139         }
2140     }
2141
2142     Maasha::Align::align_invert( \@entries, $options->{ "soft" } );
2143
2144     map { Maasha::Biopieces::put_record( { SEQ_NAME => $_->[ SEQ_NAME ], SEQ => $_->[ SEQ ] }, $out ) } @entries;
2145 }
2146
2147
2148 sub script_patscan_seq
2149 {
2150     # Martin A. Hansen, August 2007.
2151
2152     # Locates patterns in sequences using scan_for_matches.
2153
2154     my ( $in,        # handle to in stream
2155          $out,       # handle to out stream
2156          $options,   # options hash
2157        ) = @_;
2158
2159     # Returns nothing.
2160
2161     my ( $genome_file, @args, $arg, $type, $seq_file, $pat_file, $out_file, $fh_in, $fh_out, $record, $patterns, $pattern, $entry, $result, %head_hash, $i );
2162
2163     if ( $options->{ "patterns" } ) {
2164         $patterns = Maasha::Patscan::parse_patterns( $options->{ "patterns" } );
2165     } elsif ( -f $options->{ "patterns_in" } ) {
2166         $patterns = Maasha::Patscan::read_patterns( $options->{ "patterns_in" } );
2167     }
2168
2169     $genome_file = "$ENV{ 'BP_DATA' }/genomes/$options->{ 'genome' }/fasta/$options->{ 'genome' }.fna" if $options->{ 'genome' };
2170
2171     push @args, "-c"                            if $options->{ "comp" };
2172     push @args, "-m $options->{ 'max_hits' }"   if $options->{ 'max_hits' };
2173     push @args, "-n $options->{ 'max_misses' }" if $options->{ 'max_hits' };
2174
2175     $seq_file = "$BP_TMP/patscan.seq";
2176     $pat_file = "$BP_TMP/patscan.pat";
2177     $out_file = "$BP_TMP/patscan.out";
2178
2179     $fh_out = Maasha::Common::write_open( $seq_file );
2180
2181     $i = 0;
2182
2183     while ( $record = Maasha::Biopieces::get_record( $in ) )
2184     {
2185         if ( $record->{ "SEQ" } and $record->{ "SEQ_NAME" } )
2186         {
2187             $type = Maasha::Seq::seq_guess_type( $record->{ "SEQ" } ) if not $type;
2188
2189             Maasha::Fasta::put_entry( [ $i, $record->{ "SEQ" } ], $fh_out );
2190
2191             $head_hash{ $i } = $record->{ "SEQ_NAME" };
2192
2193             $i++;
2194         }
2195     }
2196
2197     close $fh_out;
2198
2199     $arg  = join " ", @args;
2200     $arg .= " -p" if $type eq "protein";
2201
2202     foreach $pattern ( @{ $patterns } )
2203     {
2204         $fh_out = Maasha::Common::write_open( $pat_file );
2205
2206         print $fh_out "$pattern\n";
2207
2208         close $fh_out;
2209
2210         if ( $options->{ 'genome' } ) {
2211             `scan_for_matches $arg $pat_file < $genome_file > $out_file`;
2212             # Maasha::Common::run( "scan_for_matches", "$arg $pat_file < $genome_file > $out_file" );
2213         } else {
2214             `scan_for_matches $arg $pat_file < $seq_file > $out_file`;
2215             # Maasha::Common::run( "scan_for_matches", "$arg $pat_file < $seq_file > $out_file" );
2216         }
2217
2218         $fh_in = Maasha::Common::read_open( $out_file );
2219
2220         while ( $entry = Maasha::Fasta::get_entry( $fh_in ) )
2221         {
2222             $result = Maasha::Patscan::parse_scan_result( $entry, $pattern );
2223
2224             if ( $options->{ 'genome' } )
2225             {
2226                 $result->{ "CHR" }     = $result->{ "S_ID" };
2227                 $result->{ "CHR_BEG" } = $result->{ "S_BEG" };
2228                 $result->{ "CHR_END" } = $result->{ "S_END" };
2229
2230                 delete $result->{ "S_ID" };
2231                 delete $result->{ "S_BEG" };
2232                 delete $result->{ "S_END" };
2233             }
2234             else
2235             {
2236                 $result->{ "S_ID" } = $head_hash{ $result->{ "S_ID" } };
2237             }
2238
2239             Maasha::Biopieces::put_record( $result, $out );
2240         }
2241
2242         close $fh_in;
2243     }
2244
2245     unlink $pat_file;
2246     unlink $seq_file;
2247     unlink $out_file;
2248 }
2249
2250
2251 sub script_soap_seq
2252 {
2253     # Martin A. Hansen, July 2008.
2254
2255     # soap sequences in stream against a given file or genome.
2256
2257     my ( $in,        # handle to in stream
2258          $out,       # handle to out stream
2259          $options,   # options hash
2260        ) = @_;
2261
2262     # Returns nothing.
2263
2264     my ( $genome, $tmp_in, $tmp_out, $fh_in, $fh_out, $record, $line, @fields, $entry, $count, $args );
2265
2266     $options->{ "seed_size" }  ||= 10;
2267     $options->{ "mismatches" } ||= 2;
2268     $options->{ "gap_size" }   ||= 0;
2269     $options->{ "cpus" }       ||= 1;
2270
2271     if ( $options->{ "genome" } ) {
2272         $options->{ "in_file" } = "$ENV{ 'BP_DATA' }/genomes/$options->{ 'genome' }/fasta/$options->{ 'genome' }.fna";
2273     }
2274
2275     $tmp_in  = "$BP_TMP/soap_query.seq";
2276     $tmp_out = "$BP_TMP/soap.result";
2277
2278     $fh_out = Maasha::Common::write_open( $tmp_in );
2279
2280     $count = 0;
2281
2282     while ( $record = Maasha::Biopieces::get_record( $in ) )
2283     {
2284         if ( $entry = Maasha::Fasta::biopiece2fasta( $record ) )
2285         {
2286             Maasha::Fasta::put_entry( $entry, $fh_out );
2287
2288             $count++;
2289         }
2290
2291         Maasha::Biopieces::put_record( $record, $out );
2292     }
2293
2294     close $fh_out;
2295
2296     if ( $count > 0 )
2297     {
2298         $args = join( " ",
2299             "-s $options->{ 'seed_size' }",
2300             "-r 2",
2301             "-a $tmp_in",
2302             "-v $options->{ 'mismatches' }",
2303             "-g $options->{ 'gap_size' }",
2304             "-p $options->{ 'cpus' }",
2305             "-d $options->{ 'in_file' }",
2306             "-o $tmp_out",
2307         );
2308
2309         $args .= " > /dev/null 2>&1" if not $options->{ 'verbose' };
2310
2311         Maasha::Common::run( "soap", $args, 1 );
2312
2313         unlink $tmp_in;
2314
2315         $fh_out = Maasha::Common::read_open( $tmp_out );
2316
2317         undef $record;
2318
2319         while ( $line = <$fh_out> )
2320         {
2321             chomp $line;
2322
2323             @fields = split /\t/, $line;
2324
2325             $record->{ "REC_TYPE" }   = "SOAP";
2326             $record->{ "Q_ID" }       = $fields[ 0 ];
2327             $record->{ "SCORE" }      = $fields[ 3 ];
2328             $record->{ "STRAND" }     = $fields[ 6 ];
2329             $record->{ "S_ID" }       = $fields[ 7 ];
2330             $record->{ "S_BEG" }      = $fields[ 8 ] - 1; # soap is 1-based
2331             $record->{ "S_END" }      = $fields[ 8 ] + $fields[ 5 ] - 2;
2332
2333             Maasha::Biopieces::put_record( $record, $out );
2334         }
2335
2336         close $fh_out;
2337     }
2338
2339     unlink $tmp_out;
2340 }
2341
2342
2343 sub script_match_seq
2344 {
2345     # Martin A. Hansen, August 2007.
2346
2347     # BLATs sequences in stream against a given genome.
2348
2349     my ( $in,        # handle to in stream
2350          $out,       # handle to out stream
2351          $options,   # options hash
2352        ) = @_;
2353
2354     # Returns nothing.
2355
2356     my ( $record, @entries, $results );
2357
2358     $options->{ "word_size" } ||= 20;
2359     $options->{ "direction" } ||= "both";
2360
2361     while ( $record = Maasha::Biopieces::get_record( $in ) )
2362     {
2363         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } ) {
2364             push @entries, [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ];
2365         }
2366
2367         Maasha::Biopieces::put_record( $record, $out );
2368     }
2369
2370     if ( @entries == 1 )
2371     {
2372         $results = Maasha::Match::match_mummer( [ $entries[ 0 ] ], [ $entries[ 0 ] ], $options, $BP_TMP );
2373
2374         map { Maasha::Biopieces::put_record( $_, $out ) } @{ $results };
2375     }
2376     elsif ( @entries == 2 )
2377     {
2378         $results = Maasha::Match::match_mummer( [ $entries[ 0 ] ], [ $entries[ 1 ] ], $options, $BP_TMP );
2379
2380         map { Maasha::Biopieces::put_record( $_, $out ) } @{ $results };
2381     }
2382 }
2383
2384
2385 sub script_create_vmatch_index
2386 {
2387     # Martin A. Hansen, January 2008.
2388
2389     # Create a vmatch index from sequences in the stream.
2390
2391     my ( $in,        # handle to in stream
2392          $out,       # handle to out stream
2393          $options,   # options hash
2394        ) = @_;
2395
2396     # Returns nothing.
2397
2398     my ( $record, $file_tmp, $fh_tmp, $type, $entry );
2399
2400     if ( $options->{ "index_name" } )
2401     {
2402         $file_tmp = $options->{ 'index_name' };
2403         $fh_tmp   = Maasha::Common::write_open( $file_tmp );
2404     }
2405
2406     while ( $record = Maasha::Biopieces::get_record( $in ) )
2407     {
2408         if ( $options->{ "index_name" } and $entry = Maasha::Fasta::biopiece2fasta( $record ) )
2409         {
2410             Maasha::Fasta::put_entry( $entry, $fh_tmp );
2411
2412             $type = Maasha::Seq::seq_guess_type( $entry->[ SEQ ] ) if not defined $type;
2413         }
2414
2415         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
2416     }
2417
2418     if ( $options->{ "index_name" } )
2419     {
2420         close $fh_tmp;
2421
2422         if ( $type eq "protein" ) {
2423             Maasha::Common::run( "mkvtree", "-db $file_tmp -protein -pl $options->{ 'prefix_length' } -allout -indexname $file_tmp > /dev/null 2>&1" );
2424         } else {
2425             Maasha::Common::run( "mkvtree", "-db $file_tmp -dna -pl $options->{ 'prefix_length' } -allout -indexname $file_tmp > /dev/null 2>&1" );
2426         }
2427
2428         unlink $file_tmp;
2429     }
2430 }
2431
2432
2433 sub script_write_bed
2434 {
2435     # Martin A. Hansen, August 2007.
2436
2437     # Write BED format for the UCSC genome browser using records in stream.
2438
2439     my ( $in,        # handle to in stream
2440          $out,       # handle to out stream
2441          $options,   # options hash
2442        ) = @_;
2443
2444     # Returns nothing.
2445
2446     my ( $cols, $fh, $record, $bed_entry, $new_record );
2447
2448     $cols = $options->{ 'cols' }->[ 0 ];
2449
2450     $fh = Maasha::Biopieces::write_stream( $options->{ 'data_out' }, $options->{ 'compress' } );
2451
2452     while ( $record = Maasha::Biopieces::get_record( $in ) )
2453     {
2454         $record = Maasha::UCSC::psl2record( $record ) if $record->{ 'tBaseInsert' }; # Dirty addition to allow Affy data from MySQL to be dumped
2455
2456         if ( $bed_entry = Maasha::UCSC::BED::biopiece2bed( $record, $cols ) ) {
2457             Maasha::UCSC::BED::bed_entry_put( $bed_entry, $fh, $cols, $options->{ 'check' } );
2458         }
2459
2460         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ 'no_stream' };
2461     }
2462
2463     close $fh;
2464 }
2465
2466
2467 sub script_write_psl
2468 {
2469     # Martin A. Hansen, August 2007.
2470
2471     # Write PSL output from stream.
2472
2473     my ( $in,        # handle to in stream
2474          $out,       # handle to out stream
2475          $options,   # options hash
2476        ) = @_;
2477
2478     # Returns nothing.
2479
2480     my ( $fh, $record, @output, $first );
2481
2482     $first = 1;
2483
2484     $fh = write_stream( $options->{ "data_out" }, $options->{ "compress" } );
2485
2486     while ( $record = Maasha::Biopieces::get_record( $in ) )
2487     {
2488         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
2489
2490         if ( $record->{ "REC_TYPE" } and $record->{ "REC_TYPE" } eq "PSL" )
2491         {
2492             Maasha::UCSC::psl_put_header( $fh ) if $first;
2493             Maasha::UCSC::psl_put_entry( $record, $fh );
2494             $first = 0;
2495         }
2496     }
2497
2498     close $fh;
2499 }
2500
2501
2502 sub script_write_fixedstep
2503 {
2504     # Martin A. Hansen, Juli 2008.
2505
2506     # Write fixedStep entries from recrods in the stream.
2507
2508     my ( $in,        # handle to in stream
2509          $out,       # handle to out stream
2510          $options,   # options hash
2511        ) = @_;
2512
2513     # Returns nothing.
2514
2515     my ( $fh, $record, $entry );
2516
2517     $fh = write_stream( $options->{ "data_out" }, $options->{ "compress" } );
2518
2519     while ( $record = Maasha::Biopieces::get_record( $in ) )
2520     {
2521         if ( $entry = Maasha::UCSC::Wiggle::biopiece2fixedstep( $record ) ) {
2522             Maasha::UCSC::Wiggle::fixedstep_entry_put( $entry, $fh );
2523         }
2524
2525         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
2526     }
2527
2528     close $fh;
2529 }
2530
2531
2532 sub script_write_2bit
2533 {
2534     # Martin A. Hansen, March 2008.
2535
2536     # Write sequence entries from stream in 2bit format.
2537
2538     my ( $in,        # handle to in stream
2539          $out,       # handle to out stream
2540          $options,   # options hash
2541        ) = @_;
2542
2543     # Returns nothing.
2544
2545     my ( $record, $mask, $tmp_file, $fh_tmp, $fh_in, $fh_out, $entry );
2546
2547     $mask = 1 if not $options->{ "no_mask" };
2548
2549     $tmp_file = "$BP_TMP/write_2bit.fna";
2550     $fh_tmp   = Maasha::Common::write_open( $tmp_file );
2551
2552     $fh_out = write_stream( $options->{ "data_out" } );
2553
2554     while ( $record = Maasha::Biopieces::get_record( $in ) )
2555     {
2556         if ( $entry = Maasha::Fasta::biopiece2fasta( $record ) ) {
2557             Maasha::Fasta::put_entry( $entry, $fh_tmp );
2558         }
2559
2560         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
2561     }
2562
2563     close $fh_tmp;
2564
2565     $fh_in = Maasha::Common::read_open( $tmp_file );
2566
2567     Maasha::TwoBit::fasta2twobit( $fh_in, $fh_out, $mask );
2568
2569     close $fh_in;
2570     close $fh_out;
2571
2572     unlink $tmp_file;
2573 }
2574
2575
2576 sub script_write_solid
2577 {
2578     # Martin A. Hansen, April 2008.
2579
2580     # Write di-base encoded Solid sequence from entries in stream.
2581
2582     my ( $in,        # handle to in stream
2583          $out,       # handle to out stream
2584          $options,   # options hash
2585        ) = @_;
2586
2587     # Returns nothing.
2588
2589     my ( $record, $fh, $entry );
2590
2591     $fh = write_stream( $options->{ "data_out" }, $options->{ "compress" } );
2592
2593     while ( $record = Maasha::Biopieces::get_record( $in ) )
2594     {
2595         if ( $entry = Maasha::Fasta::biopiece2fasta( $record ) )
2596         {
2597             $entry->[ SEQ ] = Maasha::Solid::seq2color_space( uc $entry->[ SEQ ] );
2598
2599             Maasha::Fasta::put_entry( $entry, $fh, $options->{ "wrap" } );
2600         }
2601
2602         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
2603     }
2604
2605     close $fh;
2606 }
2607
2608
2609 sub script_write_ucsc_config
2610 {
2611     # Martin A. Hansen, November 2008.
2612
2613     # Write UCSC Genome Broser configuration (.ra file type) from
2614     # records in the stream.
2615
2616     my ( $in,        # handle to in stream
2617          $out,       # handle to out stream
2618          $options,   # options hash
2619        ) = @_;
2620
2621     # Returns nothing.
2622
2623     my ( $record, $fh );
2624
2625     $fh = write_stream( $options->{ "data_out" } );
2626
2627     while ( $record = Maasha::Biopieces::get_record( $in ) )
2628     {
2629         Maasha::UCSC::ucsc_config_entry_put( $record, $fh ) if $record->{ "REC_TYPE" } eq "UCSC Config";
2630
2631         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
2632     }
2633
2634     close $fh;
2635 }
2636
2637
2638 sub script_plot_seqlogo
2639 {
2640     # Martin A. Hansen, August 2007.
2641
2642     # Calculates and writes a sequence logo for alignments.
2643
2644     my ( $in,        # handle to in stream
2645          $out,       # handle to out stream
2646          $options,   # options hash
2647        ) = @_;
2648
2649     # Returns nothing.
2650
2651     my ( $record, @entries, $logo, $fh );
2652
2653     while ( $record = Maasha::Biopieces::get_record( $in ) )
2654     {
2655         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } ) {
2656             push @entries, [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ];
2657         }
2658
2659         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
2660     }
2661
2662     $logo = Maasha::Plot::seq_logo( \@entries );
2663
2664     $fh = write_stream( $options->{ "data_out" } );
2665
2666     print $fh $logo;
2667
2668     close $fh;
2669 }
2670
2671
2672 sub script_plot_phastcons_profiles
2673 {
2674     # Martin A. Hansen, January 2008.
2675
2676     # Plots PhastCons profiles.
2677
2678     my ( $in,        # handle to in stream
2679          $out,       # handle to out stream
2680          $options,   # options hash
2681        ) = @_;
2682
2683     # Returns nothing.
2684
2685     my ( $phastcons_file, $phastcons_index, $index, $fh_phastcons, $record, $scores, $AoA, $plot, $fh );
2686
2687     $options->{ "title" } ||= "PhastCons Profiles";
2688
2689     $phastcons_file  = Maasha::Config::genome_phastcons( $options->{ "genome" } );
2690     $phastcons_index = Maasha::Config::genome_phastcons_index( $options->{ "genome" } );
2691
2692     $index           = Maasha::UCSC::fixedstep_index_retrieve( $phastcons_index );
2693     $fh_phastcons    = Maasha::Common::read_open( $phastcons_file );
2694
2695     while ( $record = Maasha::Biopieces::get_record( $in ) )
2696     {
2697         if ( $record->{ "CHR" } and $record->{ "CHR_BEG" } and $record->{ "CHR_END" } )
2698         {
2699             $scores = Maasha::UCSC::fixedstep_index_lookup( $index, $fh_phastcons, $record->{ "CHR" },
2700                                                                                    $record->{ "CHR_BEG" },
2701                                                                                    $record->{ "CHR_END" },
2702                                                                                    $options->{ "flank" } );
2703
2704             push @{ $AoA }, [ @{ $scores } ];
2705         }
2706
2707         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
2708     }
2709
2710     Maasha::UCSC::phastcons_normalize( $AoA );
2711
2712     $AoA = [ [ Maasha::UCSC::phastcons_mean( $AoA ) ] ] if $options->{ "mean" };
2713     $AoA = [ [ Maasha::UCSC::phastcons_median( $AoA ) ] ] if $options->{ "median" };
2714
2715     $AoA = Maasha::Matrix::matrix_flip( $AoA );
2716
2717     $plot = Maasha::Plot::lineplot_simple( $AoA, $options, $BP_TMP );
2718
2719     $fh = write_stream( $options->{ "data_out" } );
2720
2721     print $fh "$_\n" foreach @{ $plot };
2722
2723     close $fh;
2724 }
2725
2726
2727 sub script_head_records
2728 {
2729     # Martin A. Hansen, August 2007.
2730
2731     # Display the first sequences in stream.
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, $count );
2741
2742     $options->{ "num" } ||= 10;
2743
2744     $count = 0;
2745
2746     while ( $record = Maasha::Biopieces::get_record( $in ) )
2747     {
2748         $count++;
2749
2750         Maasha::Biopieces::put_record( $record, $out );
2751
2752         last if $count == $options->{ "num" };
2753     }
2754 }
2755
2756
2757 sub script_remove_keys
2758 {
2759     # Martin A. Hansen, August 2007.
2760
2761     # Remove keys from stream.
2762
2763     my ( $in,        # handle to in stream
2764          $out,       # handle to out stream
2765          $options,   # options hash
2766        ) = @_;
2767
2768     # Returns nothing.
2769
2770     my ( $record, $new_record );
2771
2772     while ( $record = Maasha::Biopieces::get_record( $in ) )
2773     {
2774         if ( $options->{ "keys" } )
2775         {
2776             map { delete $record->{ $_ } } @{ $options->{ "keys" } };
2777         }
2778         elsif ( $options->{ "save_keys" } )
2779         {
2780             map { $new_record->{ $_ } = $record->{ $_ } if exists $record->{ $_ } } @{ $options->{ "save_keys" } };
2781
2782             $record = $new_record;
2783         }
2784
2785         Maasha::Biopieces::put_record( $record, $out ) if keys %{ $record };
2786     }
2787 }
2788
2789
2790 sub script_remove_adaptor
2791 {
2792     # Martin A. Hansen, August 2008.
2793
2794     # Find and remove adaptor from sequences in the stream.
2795
2796     my ( $in,        # handle to in stream
2797          $out,       # handle to out stream
2798          $options,   # options hash
2799        ) = @_;
2800
2801     # Returns nothing.
2802
2803     my ( $record, $adaptor, $seq, $adaptor_len, $seq_len, $offset, $max_match, $max_mismatch, $pos );
2804
2805     $options->{ "remove" } ||= "after";
2806
2807     $max_mismatch = $options->{ "mismatches" } || 0;
2808     $offset       = $options->{ "offset" };
2809
2810     if ( not defined $offset ) {
2811         $offset = 0;
2812     } else {
2813         $offset--;
2814     }
2815
2816     $adaptor      = uc $options->{ "adaptor" };
2817     $adaptor_len  = length $adaptor;
2818
2819     while ( $record = Maasha::Biopieces::get_record( $in ) )
2820     {
2821         if ( $record->{ "SEQ" } )
2822         {
2823             $seq     = uc $record->{ "SEQ" };
2824             $seq_len = length $seq;
2825
2826             $pos = Maasha::Common::index_m( $seq, $adaptor, $seq_len, $adaptor_len, $offset, $max_mismatch );
2827
2828             $record->{ "ADAPTOR_POS" } = $pos;
2829
2830             if ( $pos >= 0 and $options->{ "remove" } ne "skip" )
2831             {
2832                 if ( $options->{ "remove" } eq "after" )
2833                 {
2834                     $record->{ "SEQ" }     = substr $record->{ "SEQ" }, 0, $pos;
2835                     $record->{ "SEQ_LEN" } = $pos;
2836                 }
2837                 else
2838                 {
2839                     $record->{ "SEQ" }     = substr $record->{ "SEQ" }, $pos + $adaptor_len;
2840                     $record->{ "SEQ_LEN" } = length $record->{ "SEQ" };
2841                 }
2842             }
2843
2844             Maasha::Biopieces::put_record( $record, $out );
2845         }
2846         else
2847         {
2848             Maasha::Biopieces::put_record( $record, $out );
2849         }
2850     }
2851 }
2852
2853
2854 sub script_remove_mysql_tables
2855 {
2856     # Martin A. Hansen, November 2008.
2857
2858     # Remove MySQL tables from values in stream.
2859
2860     my ( $in,        # handle to in stream
2861          $out,       # handle to out stream
2862          $options,   # options hash
2863        ) = @_;
2864
2865     # Returns nothing.
2866
2867     my ( $record, %table_hash, $dbh, $table );
2868
2869     $options->{ "user" }     ||= Maasha::UCSC::ucsc_get_user();
2870     $options->{ "password" } ||= Maasha::UCSC::ucsc_get_password();
2871
2872     map { $table_hash{ $_ } = 1 } @{ $options->{ 'tables' } };
2873
2874     while ( $record = Maasha::Biopieces::get_record( $in ) )
2875     {
2876         map { $table_hash{ $record->{ $_ } } = 1 } @{ $options->{ 'keys' } };
2877
2878         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ 'no_stream' };
2879     }
2880
2881     $dbh = Maasha::SQL::connect( $options->{ "database" }, $options->{ "user" }, $options->{ "password" } );
2882
2883     foreach $table ( sort keys %table_hash )
2884     {
2885         if ( Maasha::SQL::table_exists( $dbh, $table ) )
2886         {
2887             print STDERR qq(Removing table "$table" from database "$options->{ 'database' }" ... ) if $options->{ 'verbose' };
2888             Maasha::SQL::delete_table( $dbh, $table );
2889             print STDERR "done.\n" if $options->{ 'verbose' };
2890         }
2891         else
2892         {
2893             print STDERR qq(WARNING: table "$table" not found in database "$options->{ 'database' }\n");
2894         }
2895     }
2896
2897     Maasha::SQL::disconnect( $dbh );
2898 }
2899
2900
2901 sub script_remove_ucsc_tracks
2902 {
2903     # Martin A. Hansen, November 2008.
2904
2905     # Remove track from MySQL tables and config file.
2906
2907     my ( $in,        # handle to in stream
2908          $out,       # handle to out stream
2909          $options,   # options hash
2910        ) = @_;
2911
2912     # Returns nothing.
2913
2914     my ( $record, %track_hash, $fh_in, $fh_out, $track, @tracks, @new_tracks, $dbh );
2915
2916     $options->{ 'user' }        ||= Maasha::UCSC::ucsc_get_user();
2917     $options->{ 'password' }    ||= Maasha::UCSC::ucsc_get_password();
2918     $options->{ 'config_file' } ||= "$ENV{ 'HOME' }/ucsc/my_tracks.ra";
2919
2920     map { $track_hash{ $_ } = 1 } @{ $options->{ 'tracks' } };
2921
2922     while ( $record = Maasha::Biopieces::get_record( $in ) )
2923     {
2924         map { $track_hash{ $record->{ $_ } } = 1 } @{ $options->{ 'keys' } };
2925
2926         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ 'no_stream' };
2927     }
2928
2929     $fh_in = Maasha::Filesys::file_read_open( $options->{ 'config_file' } );
2930
2931     while ( $track = Maasha::UCSC::ucsc_config_entry_get( $fh_in ) ) {
2932         push @tracks, $track;
2933     }
2934
2935     close $fh_in;
2936
2937     foreach $track ( @tracks )
2938     {
2939         if ( $track->{ 'database' } eq $options->{ 'database' } and exists $track_hash{ $track->{ 'track' } } ) {
2940             print STDERR qq(Removing track "$track->{ 'track' }" from config file.\n) if $options->{ 'verbose' };
2941         } else {
2942             push @new_tracks, $track;
2943         }
2944     }
2945
2946     rename "$options->{ 'config_file' }", "$options->{ 'config_file' }~";
2947
2948     $fh_out = Maasha::Filesys::file_write_open( $options->{ 'config_file' } );
2949
2950     map { Maasha::UCSC::ucsc_config_entry_put( $_, $fh_out ) } @new_tracks;
2951
2952     close $fh_out;
2953
2954     # ---- locate track in database ----
2955
2956     $dbh = Maasha::SQL::connect( $options->{ "database" }, $options->{ "user" }, $options->{ "password" } );
2957
2958     foreach $track ( sort keys %track_hash )
2959     {
2960         if ( Maasha::SQL::table_exists( $dbh, $track ) )
2961         {
2962             print STDERR qq(Removing table "$track" from database "$options->{ 'database' }" ... ) if $options->{ 'verbose' };
2963             Maasha::SQL::delete_table( $dbh, $track );
2964             print STDERR "done.\n" if $options->{ 'verbose' };
2965         }
2966         else
2967         {
2968             print STDERR qq(WARNING: table "$track" not found in database "$options->{ 'database' }\n");
2969         }
2970     }
2971
2972     Maasha::SQL::disconnect( $dbh );
2973
2974     Maasha::Common::run( "ucscMakeTracks.pl", "-b > /dev/null 2>&1" );
2975 }
2976
2977
2978 sub script_rename_keys
2979 {
2980     # Martin A. Hansen, August 2007.
2981
2982     # Rename keys in stream.
2983
2984     my ( $in,        # handle to in stream
2985          $out,       # handle to out stream
2986          $options,   # options hash
2987        ) = @_;
2988
2989     # Returns nothing.
2990
2991     my ( $record );
2992
2993     while ( $record = Maasha::Biopieces::get_record( $in ) )
2994     {
2995         if ( exists $record->{ $options->{ "keys" }->[ 0 ] } )
2996         {
2997             $record->{ $options->{ "keys" }->[ 1 ] } = $record->{ $options->{ "keys" }->[ 0 ] };
2998
2999             delete $record->{ $options->{ "keys" }->[ 0 ] };
3000         }
3001
3002         Maasha::Biopieces::put_record( $record, $out );
3003     }
3004 }
3005
3006
3007 sub script_uniq_vals
3008 {
3009     # Martin A. Hansen, August 2007.
3010
3011     # Find unique values in stream.
3012
3013     my ( $in,        # handle to in stream
3014          $out,       # handle to out stream
3015          $options,   # options hash
3016        ) = @_;
3017
3018     # Returns nothing.
3019
3020     my ( %hash, $record );
3021
3022     while ( $record = Maasha::Biopieces::get_record( $in ) )
3023     {
3024         if ( $record->{ $options->{ "key" } } )
3025         {
3026             if ( not $hash{ $record->{ $options->{ "key" } } } and not $options->{ "invert" } )
3027             {
3028                 Maasha::Biopieces::put_record( $record, $out );
3029
3030                 $hash{ $record->{ $options->{ "key" } } } = 1;
3031             }
3032             elsif ( $hash{ $record->{ $options->{ "key" } } } and $options->{ "invert" } )
3033             {
3034                 Maasha::Biopieces::put_record( $record, $out );
3035             }
3036             else
3037             {
3038                 $hash{ $record->{ $options->{ "key" } } } = 1;
3039             }
3040         }
3041         else
3042         {
3043             Maasha::Biopieces::put_record( $record, $out );
3044         }
3045     }
3046 }
3047
3048
3049 sub script_merge_vals
3050 {
3051     # Martin A. Hansen, August 2007.
3052
3053     # Rename keys in stream.
3054
3055     my ( $in,        # handle to in stream
3056          $out,       # handle to out stream
3057          $options,   # options hash
3058        ) = @_;
3059
3060     # Returns nothing.
3061
3062     my ( $record, @join, $i );
3063
3064     $options->{ "delimit" } ||= '_';
3065
3066     while ( $record = Maasha::Biopieces::get_record( $in ) )
3067     {
3068         if ( exists $record->{ $options->{ "keys" }->[ 0 ] } )
3069         {
3070             @join = $record->{ $options->{ "keys" }->[ 0 ] };
3071
3072             for ( $i = 1; $i < @{ $options->{ "keys" } }; $i++ ) {
3073                 push @join, $record->{ $options->{ "keys" }->[ $i ] } if exists $record->{ $options->{ "keys" }->[ $i ] };
3074             }
3075
3076             $record->{ $options->{ "keys" }->[ 0 ] } = join $options->{ "delimit" }, @join;
3077         }
3078
3079         Maasha::Biopieces::put_record( $record, $out );
3080     }
3081 }
3082
3083
3084 sub script_merge_records
3085 {
3086     # Martin A. Hansen, July 2008.
3087
3088     # Merges records in the stream based on identical values of two given keys.
3089
3090     my ( $in,        # handle to in stream
3091          $out,       # handle to out stream
3092          $options,   # options hash
3093        ) = @_;
3094
3095     # Returns nothing.
3096
3097     my ( $merge, $record, $file1, $file2, $fh1, $fh2, $key1, $key2, @keys1, @keys2, @vals1, @vals2,
3098          $num1, $num2, $num, $cmp, $i );
3099
3100     $merge = $options->{ "merge" } || "AandB";
3101
3102     $file1 = "$BP_TMP/merge_records1.tmp";
3103     $file2 = "$BP_TMP/merge_records2.tmp";
3104
3105     $fh1   = Maasha::Common::write_open( $file1 );
3106     $fh2   = Maasha::Common::write_open( $file2 );
3107
3108     $key1  = $options->{ "keys" }->[ 0 ];
3109     $key2  = $options->{ "keys" }->[ 1 ];
3110
3111     $num   = $key2 =~ s/n$//;
3112     $num1  = 0;
3113     $num2  = 0;
3114
3115     while ( $record = Maasha::Biopieces::get_record( $in ) )
3116     {
3117         if ( exists $record->{ $key1 } )
3118         {
3119             @keys1 = $key1;
3120             @vals1 = $record->{ $key1 };
3121
3122             delete $record->{ $key1 };
3123
3124             map { push @keys1, $_; push @vals1, $record->{ $_ } } keys %{ $record };
3125
3126             print $fh1 join( "\t", @vals1 ), "\n";
3127
3128             $num1++;
3129         }
3130         elsif ( exists $record->{ $key2 } )
3131         {
3132             @keys2 = $key2;
3133             @vals2 = $record->{ $key2 };
3134
3135             delete $record->{ $key2 };
3136
3137             map { push @keys2, $_; push @vals2, $record->{ $_ } } keys %{ $record };
3138
3139             print $fh2 join( "\t", @vals2 ), "\n";
3140
3141             $num2++;
3142         }
3143     }
3144
3145     close $fh1;
3146     close $fh2;
3147
3148     if ( $num )
3149     {
3150         Maasha::Common::run( "sort", "-k 1,1n $file1 > $file1.sort" ) and rename "$file1.sort", $file1;
3151         Maasha::Common::run( "sort", "-k 1,1n $file2 > $file2.sort" ) and rename "$file2.sort", $file2;
3152     }
3153     else
3154     {
3155         Maasha::Common::run( "sort", "-k 1,1 $file1 > $file1.sort" ) and rename "$file1.sort", $file1;
3156         Maasha::Common::run( "sort", "-k 1,1 $file2 > $file2.sort" ) and rename "$file2.sort", $file2;
3157     }
3158
3159     $fh1 = Maasha::Common::read_open( $file1 );
3160     $fh2 = Maasha::Common::read_open( $file2 );
3161
3162     @vals1 = Maasha::Common::get_fields( $fh1 );
3163     @vals2 = Maasha::Common::get_fields( $fh2 );
3164
3165     while ( $num1 > 0 and $num2 > 0 )
3166     {
3167         undef $record;
3168
3169         if ( $num ) {
3170             $cmp = $vals1[ 0 ] <=> $vals2[ 0 ];
3171         } else {
3172             $cmp = $vals1[ 0 ] cmp $vals2[ 0 ];
3173         }
3174
3175         if ( $cmp < 0 )
3176         {
3177             if ( $merge =~ /^(AorB|AnotB)$/ )
3178             {
3179                 for ( $i = 0; $i < @keys1; $i++ ) {
3180                     $record->{ $keys1[ $i ] } = $vals1[ $i ];
3181                 }
3182
3183                 Maasha::Biopieces::put_record( $record, $out );
3184             }
3185
3186             @vals1 = Maasha::Common::get_fields( $fh1 );
3187             $num1--;
3188         }
3189         elsif ( $cmp > 0 )
3190         {
3191             if ( $merge =~ /^(BorA|BnotA)$/ )
3192             {
3193                 for ( $i = 0; $i < @keys2; $i++ ) {
3194                     $record->{ $keys2[ $i ] } = $vals2[ $i ];
3195                 }
3196
3197                 Maasha::Biopieces::put_record( $record, $out );
3198             }
3199
3200             @vals2 = Maasha::Common::get_fields( $fh2 );
3201             $num2--;
3202         }
3203         else
3204         {
3205             if ( $merge =~ /^(AandB|AorB|BorA)$/ )
3206             {
3207                 for ( $i = 0; $i < @keys1; $i++ ) {
3208                     $record->{ $keys1[ $i ] } = $vals1[ $i ];
3209                 }
3210
3211                 for ( $i = 1; $i < @keys2; $i++ ) {
3212                     $record->{ $keys2[ $i ] } = $vals2[ $i ];
3213                 }
3214
3215                 Maasha::Biopieces::put_record( $record, $out );
3216             }
3217
3218             @vals1 = Maasha::Common::get_fields( $fh1 );
3219             @vals2 = Maasha::Common::get_fields( $fh2 );
3220             $num1--;
3221             $num2--;
3222         }
3223     }
3224
3225     close $fh1;
3226     close $fh2;
3227
3228     unlink $file1;
3229     unlink $file2;
3230
3231     if ( $num1 > 0 and $merge =~ /^(AorB|AnotB)$/ )
3232     {
3233         undef $record;
3234
3235         for ( $i = 0; $i < @keys1; $i++ ) {
3236             $record->{ $keys1[ $i ] } = $vals1[ $i ];
3237         }
3238
3239         Maasha::Biopieces::put_record( $record, $out );
3240     }
3241
3242     if ( $num2 > 0 and $merge =~ /^(BorA|BnotA)$/ )
3243     {
3244         undef $record;
3245
3246         for ( $i = 0; $i < @keys2; $i++ ) {
3247             $record->{ $keys2[ $i ] } = $vals2[ $i ];
3248         }
3249
3250         Maasha::Biopieces::put_record( $record, $out );
3251     }
3252 }
3253
3254
3255 sub script_flip_tab
3256 {
3257     # Martin A. Hansen, June 2008.
3258
3259     # Flip a table.
3260
3261     my ( $in,        # handle to in stream
3262          $out,       # handle to out stream
3263          $options,   # options hash
3264        ) = @_;
3265
3266     # Returns nothing.
3267
3268     my ( $record, $key, $A, $B, @rows, @matrix, $row, $i );
3269
3270     while ( $record = Maasha::Biopieces::get_record( $in ) )
3271     {
3272         undef @rows;
3273
3274         foreach $key ( sort { $A = $a; $B = $b; $A =~ s/^V(\d+)$/$1/; $B =~ s/^V(\d+)$/$1/; $A <=> $B } keys %{ $record } )
3275         {
3276             push @rows, $record->{ $key };
3277
3278         }
3279
3280         push @matrix, [ @rows ];
3281     }
3282
3283     undef $record;
3284
3285     @matrix = Maasha::Matrix::matrix_flip( \@matrix );
3286
3287     foreach $row ( @matrix )
3288     {
3289         for ( $i = 0; $i < @{ $row }; $i++ ) {
3290             $record->{ "V$i" } = $row->[ $i ];
3291         }
3292
3293         Maasha::Biopieces::put_record( $record, $out );
3294     }
3295 }
3296
3297
3298 sub script_sort_records
3299 {
3300     # Martin A. Hansen, August 2007.
3301
3302     # Sort to sort records according to keys.
3303
3304     my ( $in,        # handle to in stream
3305          $out,       # handle to out stream
3306          $options,   # options hash
3307        ) = @_;
3308
3309     # Returns nothing.
3310
3311     my ( @keys, $key, @sort_cmd, $sort_str, $sort_sub, @records, $record, $i );
3312
3313     foreach $key ( @{ $options->{ "keys" } } )
3314     {
3315         if ( $key =~ s/n$// ) {
3316             push @sort_cmd, qq(\$a->{ "$key" } <=> \$b->{ "$key" });
3317         } else {
3318             push @sort_cmd, qq(\$a->{ "$key" } cmp \$b->{ "$key" });
3319         }
3320     }
3321
3322     $sort_str = join " or ", @sort_cmd;
3323     $sort_sub = eval "sub { $sort_str }";   # NB security issue!
3324
3325     while ( $record = Maasha::Biopieces::get_record( $in ) ) {
3326         push @records, $record;
3327     }
3328
3329     @records = sort $sort_sub @records;
3330
3331     if ( $options->{ "reverse" } )
3332     {
3333         for ( $i = scalar @records - 1; $i >= 0; $i-- ) {
3334             Maasha::Biopieces::put_record( $records[ $i ], $out );
3335         }
3336     }
3337     else
3338     {
3339         for ( $i = 0; $i < scalar @records; $i++ ) {
3340             Maasha::Biopieces::put_record( $records[ $i ], $out );
3341         }
3342     }
3343 }
3344
3345
3346 sub script_plot_histogram
3347 {
3348     # Martin A. Hansen, September 2007.
3349
3350     # Plot a simple histogram for a given key using GNU plot.
3351
3352     my ( $in,        # handle to in stream
3353          $out,       # handle to out stream
3354          $options,   # options hash
3355        ) = @_;
3356
3357     # Returns nothing.
3358
3359     my ( $record, %data_hash, $max, @data_list, $i, $result, $fh );
3360
3361     $options->{ "title" } ||= "Histogram";
3362     $options->{ "sort" }  ||= "num";
3363
3364     while ( $record = Maasha::Biopieces::get_record( $in ) )
3365     {
3366         $data_hash{ $record->{ $options->{ "key" } } }++ if defined $record->{ $options->{ "key" } };
3367
3368         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
3369     }
3370
3371     if ( $options->{ "sort" } eq "num" ) {
3372         map { push @data_list, [ $_, $data_hash{ $_ } ] } sort { $a <=> $b } keys %data_hash;
3373     } else {
3374         map { push @data_list, [ $_, $data_hash{ $_ } ] } sort keys %data_hash;
3375     }
3376
3377     $result = Maasha::Plot::histogram_simple( \@data_list, $options );
3378
3379     $fh = write_stream( $options->{ "data_out" } );
3380
3381     print $fh "$_\n" foreach @{ $result };
3382
3383     close $fh;
3384 }
3385
3386
3387 sub script_plot_lendist
3388 {
3389     # Martin A. Hansen, August 2007.
3390
3391     # Plot length distribution using GNU plot.
3392
3393     my ( $in,        # handle to in stream
3394          $out,       # handle to out stream
3395          $options,   # options hash
3396        ) = @_;
3397
3398     # Returns nothing.
3399
3400     my ( $record, %data_hash, $max, @data_list, $i, $result, $fh );
3401
3402     $options->{ "title" } ||= "Length Distribution";
3403
3404     while ( $record = Maasha::Biopieces::get_record( $in ) )
3405     {
3406         $data_hash{ $record->{ $options->{ "key" } } }++ if defined $record->{ $options->{ "key" } };
3407
3408         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
3409     }
3410
3411     $max = Maasha::Calc::list_max( [ keys %data_hash ] );
3412
3413     for ( $i = 0; $i < $max; $i++ ) {
3414         push @data_list, [ $i, $data_hash{ $i } || 0 ];
3415     }
3416
3417     $result = Maasha::Plot::histogram_lendist( \@data_list, $options );
3418
3419     $fh = write_stream( $options->{ "data_out" } );
3420
3421     print $fh "$_\n" foreach @{ $result };
3422
3423     close $fh;
3424 }
3425
3426
3427 sub script_plot_chrdist
3428 {
3429     # Martin A. Hansen, August 2007.
3430
3431     # Plot chromosome distribution using GNU plot.
3432
3433     my ( $in,        # handle to in stream
3434          $out,       # handle to out stream
3435          $options,   # options hash
3436        ) = @_;
3437
3438     # Returns nothing.
3439
3440     my ( $record, %data_hash, @data_list, $elem, $sort_key, $count, $result, $fh );
3441
3442     $options->{ "title" } ||= "Chromosome Distribution";
3443
3444     while ( $record = Maasha::Biopieces::get_record( $in ) )
3445     {
3446         if ( $record->{ "CHR" } ) {                                                             # generic
3447             $data_hash{ $record->{ "CHR" } }++;
3448         } elsif ( $record->{ "REC_TYPE" } eq "PATSCAN" and $record->{ "S_ID" } =~ /^chr/i ) {   # patscan
3449             $data_hash{ $record->{ "S_ID" } }++;
3450         } elsif ( $record->{ "REC_TYPE" } eq "PSL" and $record->{ "S_ID" } =~ /^chr/i ) {       # BLAT / PSL
3451             $data_hash{ $record->{ "S_ID" } }++;
3452         } elsif ( $record->{ "REC_TYPE" } eq "BLAST" and $record->{ "S_ID" } =~ /^chr/i ) {     # BLAST
3453             $data_hash{ $record->{ "S_ID" } }++;
3454         }
3455
3456         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
3457     }
3458
3459     foreach $elem ( keys %data_hash )
3460     {
3461         $sort_key = $elem;
3462
3463         $sort_key =~ s/chr//i;
3464
3465         $sort_key =~ s/^X(.*)/99$1/;
3466         $sort_key =~ s/^Y(.*)/99$1/;
3467         $sort_key =~ s/^Z(.*)/999$1/;
3468         $sort_key =~ s/^M(.*)/9999$1/;
3469         $sort_key =~ s/^U(.*)/99999$1/;
3470
3471         $count = $sort_key =~ tr/_//;
3472
3473         $sort_key =~ s/_.*/"999999" x $count/ex;
3474
3475         push @data_list, [ $elem, $data_hash{ $elem }, $sort_key ];
3476     }
3477
3478     @data_list = sort { $a->[ 2 ] <=> $b->[ 2 ] } @data_list;
3479
3480     $result = Maasha::Plot::histogram_chrdist( \@data_list, $options );
3481
3482     $fh = write_stream( $options->{ "data_out" } );
3483
3484     print $fh "$_\n" foreach @{ $result };
3485
3486     close $fh;
3487 }
3488
3489
3490 sub script_plot_karyogram
3491 {
3492     # Martin A. Hansen, August 2007.
3493
3494     # Plot hits on karyogram.
3495
3496     my ( $in,        # handle to in stream
3497          $out,       # handle to out stream
3498          $options,   # options hash
3499        ) = @_;
3500
3501     # Returns nothing.
3502
3503     my ( %options, $record, @data, $fh, $result, %data_hash );
3504
3505     $options->{ "genome" }     ||= "human";
3506     $options->{ "feat_color" } ||= "black";
3507
3508     while ( $record = Maasha::Biopieces::get_record( $in ) )
3509     {
3510         if ( $record->{ "CHR" } and $record->{ "CHR_BEG" } and $record->{ "CHR_END" } )
3511         {
3512             push @{ $data_hash{ $record->{ "CHR" } } }, [ $record->{ "CHR_BEG" }, $record->{ "CHR_END" }, $options->{ "feat_color" } ];
3513         }
3514
3515         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
3516     }
3517
3518     $result = Maasha::Plot::karyogram( \%data_hash, \%options );
3519
3520     $fh = write_stream( $options->{ "data_out" } );
3521
3522     print $fh $result;
3523
3524     close $fh;
3525 }
3526
3527
3528 sub script_plot_matches
3529 {
3530     # Martin A. Hansen, August 2007.
3531
3532     # Plot matches in 2D generating a dotplot.
3533
3534     my ( $in,        # handle to in stream
3535          $out,       # handle to out stream
3536          $options,   # options hash
3537        ) = @_;
3538
3539     # Returns nothing.
3540
3541     my ( $record, @data, $fh, $result, %data_hash );
3542
3543     $options->{ "direction" } ||= "both";
3544
3545     while ( $record = Maasha::Biopieces::get_record( $in ) )
3546     {
3547         if ( defined $record->{ "Q_BEG" } and defined $record->{ "S_BEG" } and $record->{ "Q_END" } and $record->{ "S_END" } ) {
3548             push @data, $record;
3549         }
3550
3551         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
3552     }
3553
3554     $options->{ "title" }  ||= "plot_matches";
3555     $options->{ "xlabel" } ||= $data[ 0 ]->{ "Q_ID" };
3556     $options->{ "ylabel" } ||= $data[ 0 ]->{ "S_ID" };
3557
3558     $result = Maasha::Plot::dotplot_matches( \@data, $options, $BP_TMP );
3559
3560     $fh = write_stream( $options->{ "data_out" } );
3561
3562     print $fh "$_\n" foreach @{ $result };
3563
3564     close $fh;
3565 }
3566
3567
3568 sub script_length_vals
3569 {
3570     # Martin A. Hansen, August 2007.
3571
3572     # Determine the length of the value for given keys.
3573
3574     my ( $in,        # handle to in stream
3575          $out,       # handle to out stream
3576          $options,   # options hash
3577        ) = @_;
3578
3579     # Returns nothing.
3580
3581     my ( $record, $key );
3582
3583     while ( $record = Maasha::Biopieces::get_record( $in ) )
3584     {
3585         foreach $key ( @{ $options->{ "keys" } } )
3586         {
3587             if ( $record->{ $key } ) {
3588                 $record->{ $key . "_LEN" } = length $record->{ $key };
3589             }
3590         }
3591
3592         Maasha::Biopieces::put_record( $record, $out );
3593     }
3594 }
3595
3596
3597 sub script_sum_vals
3598 {
3599     # Martin A. Hansen, August 2007.
3600
3601     # Calculates the sums for values of given keys.
3602
3603     my ( $in,        # handle to in stream
3604          $out,       # handle to out stream
3605          $options,   # options hash
3606        ) = @_;
3607
3608     # Returns nothing.
3609
3610     my ( $record, $key, %sum_hash, $fh );
3611
3612     while ( $record = Maasha::Biopieces::get_record( $in ) )
3613     {
3614         foreach $key ( @{ $options->{ "keys" } } )
3615         {
3616             if ( $record->{ $key } ) {
3617                 $sum_hash{ $key } += $record->{ $key };
3618             }
3619         }
3620
3621         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
3622     }
3623
3624     $fh = write_stream( $options->{ "data_out" } );
3625
3626     foreach $key ( @{ $options->{ "keys" } } ) {
3627         Maasha::Biopieces::put_record( { $key . "_SUM" => $sum_hash{ $key } || 0 } , $fh );
3628     }
3629
3630     close $fh;
3631 }
3632
3633
3634 sub script_mean_vals
3635 {
3636     # Martin A. Hansen, August 2007.
3637
3638     # Calculate the mean of values of given keys.
3639
3640     my ( $in,        # handle to in stream
3641          $out,       # handle to out stream
3642          $options,   # options hash
3643        ) = @_;
3644
3645     # Returns nothing.
3646
3647     my ( $record, $key, %sum_hash, %count_hash, $mean, $fh );
3648
3649     while ( $record = Maasha::Biopieces::get_record( $in ) )
3650     {
3651         foreach $key ( @{ $options->{ "keys" } } )
3652         {
3653             if ( $record->{ $key } )
3654             {
3655                 $sum_hash{ $key } += $record->{ $key };
3656                 $count_hash{ $key }++;
3657             }
3658         }
3659
3660         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
3661     }
3662
3663     $fh = write_stream( $options->{ "data_out" } );
3664
3665     foreach $key ( @{ $options->{ "keys" } } )
3666     {
3667         if ( $count_hash{ $key } ) {
3668             $mean = sprintf( "%.2f", ( $sum_hash{ $key } / $count_hash{ $key } ) );
3669         } else {
3670             $mean = "N/A";
3671         }
3672
3673         Maasha::Biopieces::put_record( { $key . "_MEAN" => $mean } , $fh );
3674     }
3675
3676     close $fh;
3677 }
3678
3679
3680 sub script_median_vals
3681 {
3682     # Martin A. Hansen, March 2008.
3683
3684     # Calculate the median values of given keys.
3685
3686     my ( $in,        # handle to in stream
3687          $out,       # handle to out stream
3688          $options,   # options hash
3689        ) = @_;
3690
3691     # Returns nothing.
3692
3693     my ( $record, $key, %median_hash, $median, $fh );
3694
3695     while ( $record = Maasha::Biopieces::get_record( $in ) )
3696     {
3697         foreach $key ( @{ $options->{ "keys" } } ) {
3698             push @{ $median_hash{ $key } }, $record->{ $key } if defined $record->{ $key };
3699         }
3700
3701         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
3702     }
3703
3704     $fh = write_stream( $options->{ "data_out" } );
3705
3706     foreach $key ( @{ $options->{ "keys" } } )
3707     {
3708         if ( $median_hash{ $key } ) {
3709             $median = Maasha::Calc::median( $median_hash{ $key } );
3710         } else {
3711             $median = "N/A";
3712         }
3713
3714         Maasha::Biopieces::put_record( { $key . "_MEDIAN" => $median } , $fh );
3715     }
3716
3717     close $fh;
3718 }
3719
3720
3721 sub script_max_vals
3722 {
3723     # Martin A. Hansen, February 2008.
3724
3725     # Determine the maximum values of given keys.
3726
3727     my ( $in,        # handle to in stream
3728          $out,       # handle to out stream
3729          $options,   # options hash
3730        ) = @_;
3731
3732     # Returns nothing.
3733
3734     my ( $record, $key, $fh, %max_hash, $max_record );
3735
3736     while ( $record = Maasha::Biopieces::get_record( $in ) )
3737     {
3738         foreach $key ( @{ $options->{ "keys" } } )
3739         {
3740             if ( $record->{ $key } )
3741             {
3742                 $max_hash{ $key } = $record->{ $key } if $record->{ $key } > $max_hash{ $key };
3743             }
3744         }
3745
3746         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
3747     }
3748
3749     $fh = write_stream( $options->{ "data_out" } );
3750
3751     foreach $key ( @{ $options->{ "keys" } } )
3752     {
3753         $max_record->{ $key . "_MAX" } = $max_hash{ $key };
3754     }
3755
3756     Maasha::Biopieces::put_record( $max_record, $fh );
3757
3758     close $fh;
3759 }
3760
3761
3762 sub script_min_vals
3763 {
3764     # Martin A. Hansen, February 2008.
3765
3766     # Determine the minimum values of given keys.
3767
3768     my ( $in,        # handle to in stream
3769          $out,       # handle to out stream
3770          $options,   # options hash
3771        ) = @_;
3772
3773     # Returns nothing.
3774
3775     my ( $record, $key, $fh, %min_hash, $min_record );
3776
3777     while ( $record = Maasha::Biopieces::get_record( $in ) )
3778     {
3779         foreach $key ( @{ $options->{ "keys" } } )
3780         {
3781             if ( defined $record->{ $key } )
3782             {
3783                 if ( exists $min_hash{ $key } ) {
3784                     $min_hash{ $key } = $record->{ $key } if $record->{ $key } < $min_hash{ $key };
3785                 } else {
3786                     $min_hash{ $key } = $record->{ $key };
3787                 }
3788             }
3789         }
3790
3791         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
3792     }
3793
3794     $fh = write_stream( $options->{ "data_out" } );
3795
3796     foreach $key ( @{ $options->{ "keys" } } )
3797     {
3798         $min_record->{ $key . "_MIN" } = $min_hash{ $key };
3799     }
3800
3801     Maasha::Biopieces::put_record( $min_record, $fh );
3802
3803     close $fh;
3804 }
3805
3806
3807 sub script_upload_to_ucsc
3808 {
3809     # Martin A. Hansen, August 2007.
3810
3811     # Calculate the mean of values of given keys.
3812
3813     # This routine has developed into the most ugly hack. Do something!
3814
3815     my ( $in,        # handle to in stream
3816          $out,       # handle to out stream
3817          $options,   # options hash
3818        ) = @_;
3819
3820     # Returns nothing.
3821
3822     my ( $record, $file, $wib_file, $wig_file, $wib_dir, $fh_out, $i, $first, $format, $type, $columns, $append, $entry );
3823
3824     $options->{ "short_label" } ||= $options->{ 'table' };
3825     $options->{ "long_label" }  ||= $options->{ 'table' };
3826     $options->{ "group" }       ||= $ENV{ "LOGNAME" };
3827     $options->{ "priority" }    ||= 1;
3828     $options->{ "visibility" }  ||= "pack";
3829     $options->{ "color" }       ||= join( ",", int( rand( 255 ) ), int( rand( 255 ) ), int( rand( 255 ) ) );
3830     $options->{ "chunk_size" }  ||= 10_000_000_000;    # Due to 32-bit UCSC compilation really large tables cannot be loaded in one go.
3831
3832     $file   = "$BP_TMP/ucsc_upload.tmp";
3833     $append = 0;
3834     $first  = 1;
3835     $i      = 0;
3836
3837     $fh_out = Maasha::Common::write_open( $file );
3838
3839     while ( $record = Maasha::Biopieces::get_record( $in ) )
3840     {
3841         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
3842
3843         if ( $record->{ "REC_TYPE" } eq "fixed_step" )
3844         {
3845             $format = "WIGGLE";
3846
3847             if ( $entry = Maasha::UCSC::Wiggle::biopiece2fixedstep( $record ) ) {
3848                 Maasha::UCSC::Wiggle::fixedstep_entry_put( $entry, $fh_out );
3849             }
3850         }
3851         elsif ( $record->{ "REC_TYPE" } eq "PSL" )
3852         {
3853             $format = "PSL";
3854
3855             Maasha::UCSC::psl_put_header( $fh_out ) if $first;
3856             Maasha::UCSC::psl_put_entry( $record, $fh_out );
3857
3858             $first = 0;
3859         }
3860         elsif ( $record->{ "REC_TYPE" } eq "BED" and $record->{ "SEC_STRUCT" } )
3861         {
3862             # chrom chromStart  chromEnd    name    score   strand  size    secStr  conf
3863
3864             $format  = "BED_SS";
3865
3866             print $fh_out join ( "\t",
3867                 $record->{ "CHR" },
3868                 $record->{ "CHR_BEG" },
3869                 $record->{ "CHR_END" } + 1,
3870                 $record->{ "Q_ID" },
3871                 $record->{ "SCORE" },
3872                 $record->{ "STRAND" },
3873                 $record->{ "SIZE" },
3874                 $record->{ "SEC_STRUCT" },
3875                 $record->{ "CONF" },
3876             ), "\n";
3877         }
3878         elsif ( $record->{ "REC_TYPE" } eq "BED" )
3879         {
3880             $format  = "BED";
3881             $columns = $record->{ "BED_COLS" };
3882
3883             if ( $entry = Maasha::UCSC::BED::biopiece2bed( $record, $columns ) ) {
3884                 Maasha::UCSC::BED::bed_entry_put( $entry, $fh_out, $columns, $options->{ 'check' } );
3885             }
3886         }
3887         elsif ( $record->{ "REC_TYPE" } eq "PATSCAN" and $record->{ "CHR" } )
3888         {
3889             $format  = "BED";
3890             $columns = 6;
3891
3892             if ( $entry = Maasha::UCSC::BED::biopiece2bed( $record, $columns ) ) {
3893                 Maasha::UCSC::BED::bed_entry_put( $entry, $fh_out, $columns, $options->{ 'check' } );
3894             }
3895         }
3896         elsif ( $record->{ "REC_TYPE" } eq "BLAST" and $record->{ "S_ID" } =~ /^chr/ )
3897         {
3898             $format  = "BED";
3899             $columns = 6;
3900
3901             $record->{ "SCORE" } = $record->{ "BIT_SCORE" } * 1000;
3902
3903             if ( $entry = Maasha::UCSC::BED::biopiece2bed( $record, $columns ) ) {
3904                 Maasha::UCSC::BED::bed_entry_put( $entry, $fh_out, $columns, $options->{ 'check' } );
3905             }
3906         }
3907         elsif ( $record->{ "REC_TYPE" } eq "VMATCH" and $record->{ "S_ID" } =~ /^chr/i )
3908         {
3909             $format  = "BED";
3910             $columns = 6;
3911
3912             if ( $entry = Maasha::UCSC::BED::biopiece2bed( $record, $columns ) ) {
3913                 Maasha::UCSC::BED::bed_entry_put( $entry, $fh_out, $columns, $options->{ 'check' } );
3914             }
3915         }
3916
3917         if ( $i == $options->{ "chunk_size" } )
3918         {
3919             close $fh_out;
3920
3921             if ( $format eq "BED" ) {
3922                 Maasha::UCSC::bed_upload_to_ucsc( $BP_TMP, $file, $options, $append );
3923             } elsif ( $format eq "PSL" ) {
3924                 Maasha::UCSC::psl_upload_to_ucsc( $file, $options, $append );
3925             }
3926
3927             unlink $file;
3928
3929             $first = 1;
3930
3931             $append = 1;
3932
3933             $fh_out = Maasha::Common::write_open( $file );
3934         }
3935
3936         $i++;
3937     }
3938
3939     close $fh_out;
3940
3941     if ( exists $options->{ "database" } and $options->{ "table" } )
3942     {
3943         if ( $format eq "BED" )
3944         {
3945             $type = "bed $columns";
3946
3947             Maasha::UCSC::bed_upload_to_ucsc( $BP_TMP, $file, $options, $append );
3948         }
3949         elsif ( $format eq "BED_SS" )
3950         {
3951             $type = "type bed 6 +";
3952
3953             Maasha::UCSC::bed_upload_to_ucsc( $BP_TMP, $file, $options, $append );
3954         }
3955         elsif ( $format eq "PSL" )
3956         {
3957             $type = "psl";
3958
3959             Maasha::UCSC::psl_upload_to_ucsc( $file, $options, $append );
3960         }
3961         elsif ( $format eq "WIGGLE" )
3962         {
3963             $options->{ "visibility" } = "full";
3964
3965             $wig_file = "$options->{ 'table' }.wig";
3966             $wib_file = "$options->{ 'table' }.wib";
3967
3968             $wib_dir  = "$ENV{ 'HOME' }/ucsc/wib";
3969
3970             Maasha::Common::dir_create_if_not_exists( $wib_dir );
3971
3972             if ( $options->{ 'verbose' } ) {
3973                 `cd $BP_TMP && wigEncode $file $wig_file $wib_file`;
3974             } else {
3975                 `cd $BP_TMP && wigEncode $file $wig_file $wib_file > /dev/null 2>&1`;
3976             }
3977
3978             Maasha::Common::run( "mv", "$BP_TMP/$wib_file $wib_dir" );
3979
3980             unlink $file;
3981
3982             $file = $wig_file;
3983
3984             $type = "wig 0";
3985
3986             Maasha::UCSC::wiggle_upload_to_ucsc( $BP_TMP, $wib_dir, $file, $options );
3987         }
3988
3989         unlink $file;
3990
3991         Maasha::UCSC::ucsc_update_config( $options, $type );
3992     }
3993 }
3994
3995
3996 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
3997
3998 1;
3999
4000 __END__
4001