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::EMBL;
  43 use Maasha::Stockholm;
  44 use Maasha::Seq;
  45 use Maasha::Calc;
  46 use Maasha::UCSC;
  47 use Maasha::UCSC::BED;
  48 use Maasha::UCSC::Wiggle;
  49 use Maasha::NCBI;
  50 use Maasha::GFF;
  51 use Maasha::TwoBit;
  52 use Maasha::Solid;
  53 use Maasha::Solexa;
  54 use Maasha::SQL;
  55 use Maasha::Gwiki;
  56
  57 use vars qw( @ISA @EXPORT_OK );
  58
  59 require Exporter;
  60
  61 @ISA = qw( Exporter );
  62
  63 use constant {
  64     SEQ_NAME => 0,
  65     SEQ      => 1,
  66 };
  67
  68
  69 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> GLOBALS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  70
  71
  72 my ( $script, $BP_TMP );
  73
  74 $script = Maasha::Common::get_scriptname();
  75 $BP_TMP = Maasha::Common::get_tmpdir();
  76
  77
  78 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> RUN SCRIPT <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  79
  80
  81 run_script( $script );
  82
  83
  84 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> SUBROUTINES <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  85
  86
  87 sub run_script
  88 {
  89     # Martin A. Hansen, August 2007.
  90
  91     # Run a specific script.
  92
  93     my ( $script,   # script name
  94        ) = @_;
  95
  96     # Returns nothing.
  97
  98     my ( $t0, $t1, $options, $in, $out );
  99
 100     Maasha::Biopieces::log_biopiece();
 101
 102     $t0 = gettimeofday();
 103
 104     $options = get_options( $script );
 105
 106     $options->{ "SCRIPT" } = $script;
 107
 108     $script = "print_usage" if ( -t STDIN and keys %{ $options } <= 1 or $options->{ 'help' } );
 109
 110     $in  = Maasha::Biopieces::read_stream( $options->{ "stream_in" } );
 111     $out = Maasha::Biopieces::write_stream( $options->{ "stream_out" } );
 112
 113     if    ( $script eq "print_usage" )              { script_print_usage(               $in, $out, $options ) }
 114     elsif ( $script eq "read_psl" )                 { script_read_psl(                  $in, $out, $options ) }
 115     elsif ( $script eq "read_stockholm" )           { script_read_stockholm(            $in, $out, $options ) }
 116     elsif ( $script eq "read_phastcons" )           { script_read_phastcons(            $in, $out, $options ) }
 117     elsif ( $script eq "read_soft" )                { script_read_soft(                 $in, $out, $options ) }
 118     elsif ( $script eq "read_solexa" )              { script_read_solexa(               $in, $out, $options ) }
 119     elsif ( $script eq "read_solid" )               { script_read_solid(                $in, $out, $options ) }
 120     elsif ( $script eq "read_mysql" )               { script_read_mysql(                $in, $out, $options ) }
 121     elsif ( $script eq "uppercase_seq" )            { script_uppercase_seq(             $in, $out, $options ) }
 122     elsif ( $script eq "complexity_seq" )           { script_complexity_seq(            $in, $out, $options ) }
 123     elsif ( $script eq "get_genome_align" )         { script_get_genome_align(          $in, $out, $options ) }
 124     elsif ( $script eq "get_genome_phastcons" )     { script_get_genome_phastcons(      $in, $out, $options ) }
 125     elsif ( $script eq "soap_seq" )                 { script_soap_seq(                  $in, $out, $options ) }
 126     elsif ( $script eq "remove_mysql_tables" )      { script_remove_mysql_tables(       $in, $out, $options ) }
 127     elsif ( $script eq "remove_ucsc_tracks" )       { script_remove_ucsc_tracks(        $in, $out, $options ) }
 128     elsif ( $script eq "upload_to_ucsc" )           { script_upload_to_ucsc(            $in, $out, $options ) }
 129
 130     close $in if defined $in;
 131     close $out;
 132
 133     $t1 = gettimeofday();
 134
 135     print STDERR "Program: $script" . ( " " x ( 25 - length( $script ) ) ) . sprintf( "Run time: %.4f\n", ( $t1 - $t0 ) ) if $options->{ 'verbose' };
 136 }
 137
 138
 139 sub get_options
 140 {
 141     # Martin A. Hansen, February 2008.
 142
 143     # Gets options from commandline and checks these vigerously.
 144
 145     my ( $script,     # name of script
 146        ) = @_;
 147
 148     # Returns hash
 149
 150     my ( %options, @options, $opt, @genomes, $real );
 151
 152     if ( $script eq "print_usage" )
 153     {
 154         @options = qw(
 155             data_in|i=s
 156         );
 157     }
 158     elsif ( $script eq "read_psl" )
 159     {
 160         @options = qw(
 161             data_in|i=s
 162             num|n=s
 163         );
 164     }
 165     elsif ( $script eq "read_stockholm" )
 166     {
 167         @options = qw(
 168             data_in|i=s
 169             num|n=s
 170         );
 171     }
 172     elsif ( $script eq "read_phastcons" )
 173     {
 174         @options = qw(
 175             data_in|i=s
 176             num|n=s
 177             min|m=s
 178             dist|d=s
 179             threshold|t=f
 180             gap|g=s
 181         );
 182     }
 183     elsif ( $script eq "read_soft" )
 184     {
 185         @options = qw(
 186             data_in|i=s
 187             samples|s=s
 188             num|n=s
 189         );
 190     }
 191     elsif ( $script eq "read_solexa" )
 192     {
 193         @options = qw(
 194             data_in|i=s
 195             num|n=s
 196             format|f=s
 197             quality|q=s
 198         );
 199     }
 200     elsif ( $script eq "read_solid" )
 201     {
 202         @options = qw(
 203             data_in|i=s
 204             num|n=s
 205             quality|q=s
 206         );
 207     }
 208     elsif ( $script eq "read_mysql" )
 209     {
 210         @options = qw(
 211             database|d=s
 212             query|q=s
 213             user|u=s
 214             password|p=s
 215         );
 216     }
 217     elsif ( $script eq "get_genome_align" )
 218     {
 219         @options = qw(
 220             genome|g=s
 221             chr|c=s
 222             beg|b=s
 223             end|e=s
 224             len|l=s
 225             strand|s=s
 226         );
 227     }
 228     elsif ( $script eq "get_genome_phastcons" )
 229     {
 230         @options = qw(
 231             genome|g=s
 232             chr|c=s
 233             beg|b=s
 234             end|e=s
 235             len|l=s
 236             flank|f=s
 237         );
 238     }
 239     elsif ( $script eq "soap_seq" )
 240     {
 241         @options = qw(
 242             in_file|i=s
 243             genome|g=s
 244             seed_size|s=s
 245             mismatches|m=s
 246             gap_size|G=s
 247             cpus|c=s
 248         );
 249     }
 250     elsif ( $script eq "remove_mysql_tables" )
 251     {
 252         @options = qw(
 253             database|d=s
 254             tables|t=s
 255             keys|k=s
 256             user|u=s
 257             password|p=s
 258             no_stream|x
 259         );
 260     }
 261     elsif ( $script eq "remove_ucsc_tracks" )
 262     {
 263         @options = qw(
 264             database|d=s
 265             tracks|t=s
 266             keys|k=s
 267             config_file|c=s
 268             user|u=s
 269             password|p=s
 270             no_stream|x
 271         );
 272     }
 273     elsif ( $script eq "upload_to_ucsc" )
 274     {
 275         @options = qw(
 276             no_stream|x
 277             database|d=s
 278             table|t=s
 279             short_label|s=s
 280             long_label|l=s
 281             group|g=s
 282             priority|p=f
 283             use_score|u
 284             visibility|V=s
 285             color|c=s
 286             check|C
 287         );
 288     }
 289
 290     push @options, qw(
 291         stream_in|I=s
 292         stream_out|O=s
 293         verbose|v
 294         help|?
 295     );
 296
 297 #    print STDERR Dumper( \@options );
 298
 299     GetOptions(
 300         \%options,
 301         @options,
 302     );
 303
 304 #    print STDERR Dumper( \%options );
 305
 306     if ( -t STDIN && scalar( keys %options ) == 0 or $options{ "help" } ) {
 307         return wantarray ? %options : \%options;
 308     }
 309
 310     $options{ "cols" }      = [ split ",", $options{ "cols" } ]      if defined $options{ "cols" };
 311     $options{ "keys" }      = [ split ",", $options{ "keys" } ]      if defined $options{ "keys" };
 312     $options{ "no_keys" }   = [ split ",", $options{ "no_keys" } ]   if defined $options{ "no_keys" };
 313     $options{ "save_keys" } = [ split ",", $options{ "save_keys" } ] if defined $options{ "save_keys" };
 314     $options{ "quals" }     = [ split ",", $options{ "quals" } ]     if defined $options{ "quals" };
 315     $options{ "feats" }     = [ split ",", $options{ "feats" } ]     if defined $options{ "feats" };
 316     $options{ "frames" }    = [ split ",", $options{ "frames" } ]    if defined $options{ "frames" };
 317     $options{ "samples" }   = [ split ",", $options{ "samples" } ]   if defined $options{ "samples" };
 318     $options{ "tables" }    = [ split ",", $options{ "tables" } ]    if defined $options{ "tables" };
 319     $options{ "tracks" }    = [ split ",", $options{ "tracks" } ]    if defined $options{ "tracks" };
 320
 321     # ---- check arguments ----
 322
 323     if ( $options{ 'data_in' } )
 324     {
 325         $options{ "files" } = Maasha::Biopieces::getopt_files( $options{ 'data_in' } );
 326
 327         Maasha::Common::error( qq(Argument to --data_in must be a valid file or fileglob expression) ) if scalar @{ $options{ "files" } } == 0;
 328     }
 329
 330     map { Maasha::Common::error( qq(Argument to --cols must be a whole numbers - not "$_") ) if $_ !~ /^\d+$/ } @{ $options{ "cols" } } if $options{ "cols" };
 331
 332     # print STDERR Dumper( \%options );
 333
 334     $real = "beg|end|word_size|wrap|len|prefix_length|mismatches|offset|num|skip|cpus|window_size|step_size";
 335
 336     foreach $opt ( keys %options )
 337     {
 338         if ( $opt =~ /stream_in|pattern_in|exact_in/ and not -f $options{ $opt } )
 339         {
 340             Maasha::Common::error( qq(Argument to --$opt must be a valid file or fileglob expression - not "$options{ $opt }") );
 341         }
 342         elsif ( $opt =~ /$real/ and $options{ $opt } !~ /^\d+$/ )
 343         {
 344             Maasha::Common::error( qq(Argument to --$opt must be a whole number - not "$options{ $opt }") );
 345         }
 346         elsif ( $opt =~ /max_hits|max_hits|max_misses|dist|edit_dist|flank|gap|hamming_dist|priority/ and $options{ $opt } !~ /^-?\d+$/ )
 347         {
 348             Maasha::Common::error( qq(Argument to --$opt must be an integer - not "$options{ $opt }") );
 349         }
 350         elsif ( $opt =~ /identity|threshold/ and $options{ $opt } !~ /^-?(?:\d+(?:\.\d*)?|\.\d+)$/ )
 351         {
 352             Maasha::Common::error( qq(Argument to --$opt must be a decimal number - not "$options{ $opt }") );
 353         }
 354         elsif ( $opt =~ /e_val/ and $options{ $opt } !~ /^([+-]?)(?=\d|\.\d)\d*(\.\d*)?([Ee]([+-]?\d+))?$/ )
 355         {
 356             Maasha::Common::error( qq(Argument to --$opt must be a float - not "$options{ $opt }") );
 357         }
 358         elsif ( $opt =~ /strand/ and $options{ $opt } !~ /^(\+|-)$/ )
 359         {
 360             Maasha::Common::error( qq(Argument to --$opt must be "+" or "-" - not "$options{ $opt }") );
 361         }
 362         elsif ( $opt eq "genome" )
 363         {
 364             @genomes = Maasha::Common::ls_dirs( "$ENV{ 'BP_DATA' }/genomes" );
 365             map { $_ =~ s/.*\/(.+)$/$1/ } @genomes;
 366
 367             if ( not grep { $_ =~ /^$options{ $opt }$/ } @genomes ) {
 368                 Maasha::Common::error( qq(Genome $options{ $opt } not found in "$ENV{ 'BP_DATA' }/genomes/") );
 369             }
 370         }
 371         elsif ( $opt eq "terminal" and not $options{ $opt } =~ /^(svg|post|dumb|x11)/ )
 372         {
 373             Maasha::Common::error( qq(Bad --$opt argument "$options{ $opt }") );
 374         }
 375         elsif ( $opt eq "table" and $options{ $opt } =~ /(-|\.)/ )
 376         {
 377             Maasha::Common::error( qq(Character '$1' is not allowed in table name: $options{ $opt }) );
 378         }
 379         elsif ( $opt eq "merge" and $options{ $opt } !~ /^(AandB|AorB|BorA|AnotB|BnotA)$/ )
 380         {
 381             Maasha::Common::error( qq(Argument to --$opt must be AandB, AorB, BorA, AnotB, or BnotA - not "$options{ $opt }") );
 382         }
 383         elsif ( $opt eq "format" and $script eq "read_solexa" and $options{ $opt } !~ /octal|decimal/ )
 384         {
 385             Maasha::Common::error( qq(Argument to --$opt must be octal or decimal - not "$options{ $opt }") );
 386         }
 387     }
 388
 389     Maasha::Common::error( qq(no --database specified) )                if $script eq "remove_ucsc_tracks"  and not $options{ "database" };
 390     Maasha::Common::error( qq(no --in_file or --genome specified) )     if $script eq "soap_seq" and not $options{ "genome" } and not $options{ "in_file" };
 391     Maasha::Common::error( qq(both --in_file and --genome specified) )  if $script eq "soap_seq" and $options{ "genome" } and $options{ "in_file" };
 392     Maasha::Common::error( qq(no --genome specified) )                  if $script =~ /get_genome_align|get_genome_phastcons/ and not $options{ "genome" };
 393
 394     if ( $script eq "upload_to_ucsc" )
 395     {
 396         Maasha::Common::error( qq(no --database specified) ) if not $options{ "database" };
 397         Maasha::Common::error( qq(no --table specified) )    if not $options{ "table" };
 398     }
 399
 400     return wantarray ? %options : \%options;
 401 }
 402
 403
 404 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> SCRIPTS  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 405
 406
 407 sub script_print_usage
 408 {
 409     # Martin A. Hansen, January 2008.
 410
 411     # Retrieves usage information from file and
 412     # prints this nicely formatted.
 413
 414     my ( $in,        # handle to in stream
 415          $out,       # handle to out stream
 416          $options,   # options hash
 417        ) = @_;
 418
 419     # Returns nothing.
 420
 421     my ( $file, $wiki, $lines );
 422
 423     if ( $options->{ 'data_in' } ) {
 424         $file = $options->{ 'data_in' };
 425     } else {
 426         $file = join "", $ENV{ 'BP_DIR' }, "/bp_usage/", $options->{ 'SCRIPT' }, ".wiki";
 427     }
 428
 429     $wiki = Maasha::Gwiki::gwiki_read( $file );
 430
 431     ( $wiki->[ 2 ], $wiki->[ 3 ], $wiki->[ 0 ], $wiki->[ 1 ] ) = ( $wiki->[ 0 ], $wiki->[ 1 ], $wiki->[ 2 ], $wiki->[ 3 ] );
 432
 433     if ( not $options->{ "help" } ) {
 434         @{ $wiki } = grep { $_->[ 0 ]->{ 'SECTION' } =~ /Biopiece|summary|Usage|Options|Help/ } @{ $wiki };
 435     }
 436
 437     $lines = Maasha::Gwiki::gwiki2ascii( $wiki );
 438
 439     print STDERR "$_\n" foreach @{ $lines };
 440
 441     exit;
 442 }
 443
 444
 445 sub script_read_psl
 446 {
 447     # Martin A. Hansen, August 2007.
 448
 449     # Read psl table from stream or file.
 450
 451     my ( $in,        # handle to in stream
 452          $out,       # handle to out stream
 453          $options,   # options hash
 454        ) = @_;
 455
 456     # Returns nothing.
 457
 458     my ( $record, $file, $data_in, $num );
 459
 460     while ( $record = Maasha::Biopieces::get_record( $in ) ) {
 461         Maasha::Biopieces::put_record( $record, $out );
 462     }
 463
 464     $num = 1;
 465
 466     foreach $file ( @{ $options->{ "files" } } )
 467     {
 468         $data_in = Maasha::Common::read_open( $file );
 469
 470         while ( $record = Maasha::UCSC::psl_get_entry( $data_in ) )
 471         {
 472             Maasha::Biopieces::put_record( $record, $out );
 473
 474             goto NUM if $options->{ "num" } and $num == $options->{ "num" };
 475
 476             $num++;
 477         }
 478     }
 479
 480     NUM:
 481 }
 482
 483
 484 sub script_read_stockholm
 485 {
 486     # Martin A. Hansen, August 2007.
 487
 488     # Read Stockholm format.
 489
 490     my ( $in,        # handle to in stream
 491          $out,       # handle to out stream
 492          $options,   # options hash
 493        ) = @_;
 494
 495     # Returns nothing.
 496
 497     my ( $data_in, $file, $num, $entry, $record, $record_anno, $record_align, $key, $seq );
 498
 499     while ( $record = Maasha::Biopieces::get_record( $in ) ) {
 500         Maasha::Biopieces::put_record( $record, $out );
 501     }
 502
 503     $num = 1;
 504
 505     foreach $file ( @{ $options->{ "files" } } )
 506     {
 507         $data_in = Maasha::Common::read_open( $file );
 508
 509         while ( $entry = Maasha::Stockholm::get_stockholm_entry( $data_in ) )
 510         {
 511             $record = Maasha::Stockholm::parse_stockholm_entry( $entry );
 512
 513             undef $record_anno;
 514
 515             foreach $key ( keys %{ $record->{ "GF" } } ) {
 516                 $record_anno->{ $key } = $record->{ "GF" }->{ $key };
 517             }
 518
 519             $record_anno->{ "ALIGN" } = $num;
 520
 521             Maasha::Biopieces::put_record( $record_anno, $out );
 522
 523             foreach $seq ( @{ $record->{ "ALIGN" } } )
 524             {
 525                 undef $record_align;
 526
 527                 $record_align = {
 528                     SEQ_NAME  => $seq->[ 0 ],
 529                     SEQ       => $seq->[ 1 ],
 530                 };
 531
 532                 Maasha::Biopieces::put_record( $record_align, $out );
 533             }
 534
 535             goto NUM if $options->{ "num" } and $num == $options->{ "num" };
 536
 537             $num++;
 538         }
 539
 540         close $data_in;
 541     }
 542
 543     NUM:
 544
 545     close $data_in if $data_in;
 546 }
 547
 548
 549 sub script_read_phastcons
 550 {
 551     # Martin A. Hansen, December 2007.
 552
 553     # Read PhastCons format.
 554
 555     my ( $in,        # handle to in stream
 556          $out,       # handle to out stream
 557          $options,   # options hash
 558        ) = @_;
 559
 560     # Returns nothing.
 561
 562     my ( $data_in, $file, $num, $entry, @records, $record );
 563
 564     $options->{ "min" }       ||= 10;
 565     $options->{ "dist" }      ||= 25;
 566     $options->{ "threshold" } ||= 0.8;
 567     $options->{ "gap" }       ||= 5;
 568
 569     while ( $record = Maasha::Biopieces::get_record( $in ) ) {
 570         Maasha::Biopieces::put_record( $record, $out );
 571     }
 572
 573     $num = 1;
 574
 575     foreach $file ( @{ $options->{ "files" } } )
 576     {
 577         $data_in = Maasha::Common::read_open( $file );
 578
 579         while ( $entry = Maasha::UCSC::fixedstep_get_entry( $data_in ) )
 580         {
 581             @records = Maasha::UCSC::phastcons_parse_entry( $entry, $options );
 582
 583             foreach $record ( @records )
 584             {
 585                 $record->{ "REC_TYPE" } = "BED";
 586                 $record->{ "BED_LEN" }  = $record->{ "CHR_END" } - $record->{ "CHR_BEG" } + 1;
 587
 588                 Maasha::Biopieces::put_record( $record, $out );
 589
 590                 goto NUM if $options->{ "num" } and $num == $options->{ "num" };
 591
 592                 $num++;
 593             }
 594         }
 595
 596         close $data_in;
 597     }
 598
 599     NUM:
 600
 601     close $data_in if $data_in;
 602 }
 603
 604
 605 sub script_read_soft
 606 {
 607     # Martin A. Hansen, December 2007.
 608
 609     # Read soft format.
 610     # http://www.ncbi.nlm.nih.gov/geo/info/soft2.html
 611
 612     my ( $in,        # handle to in stream
 613          $out,       # handle to out stream
 614          $options,   # options hash
 615        ) = @_;
 616
 617     # Returns nothing.
 618
 619     my ( $data_in, $file, $num, $records, $record, $soft_index, $fh, @platforms, $plat_table, @samples, $sample, $old_end, $skip );
 620
 621     while ( $record = Maasha::Biopieces::get_record( $in ) ) {
 622         Maasha::Biopieces::put_record( $record, $out );
 623     }
 624
 625     $num = 1;
 626
 627     foreach $file ( @{ $options->{ "files" } } )
 628     {
 629         print STDERR "Creating index for file: $file\n" if $options->{ "verbose" };
 630
 631         $soft_index = Maasha::NCBI::soft_index_file( $file );
 632
 633         $fh         = Maasha::Common::read_open( $file );
 634
 635         @platforms  = grep { $_->{ "SECTION" } =~ /PLATFORM/ } @{ $soft_index };
 636
 637         print STDERR "Getting platform tables for file: $file\n" if $options->{ "verbose" };
 638
 639         $plat_table = Maasha::NCBI::soft_get_platform( $fh, $platforms[ 0 ]->{ "LINE_BEG" }, $platforms[ -1 ]->{ "LINE_END" } );
 640
 641         @samples    = grep { $_->{ "SECTION" } =~ /SAMPLE/ } @{ $soft_index };
 642
 643         $old_end    = $platforms[ -1 ]->{ "LINE_END" };
 644
 645         foreach $sample ( @samples )
 646         {
 647             $skip = 0;
 648             $skip = 1 if ( $options->{ "samples" } and grep { $sample->{ "SECTION" } !~ /$_/ } @{ $options->{ "samples" } } );
 649
 650             print STDERR "Getting samples for dataset: $sample->{ 'SECTION' }\n" if $options->{ "verbose" } and not $skip;
 651
 652             $records = Maasha::NCBI::soft_get_sample( $fh, $plat_table, $sample->{ "LINE_BEG" } - $old_end - 1, $sample->{ "LINE_END" } - $old_end - 1, $skip );
 653
 654             foreach $record ( @{ $records } )
 655             {
 656                 Maasha::Biopieces::put_record( $record, $out );
 657
 658                 goto NUM if $options->{ "num" } and $num == $options->{ "num" };
 659
 660                 $num++;
 661             }
 662
 663             $old_end = $sample->{ "LINE_END" };
 664         }
 665
 666         close $fh;
 667     }
 668
 669     NUM:
 670
 671     close $data_in if $data_in;
 672     close $fh if $fh;
 673 }
 674
 675
 676 sub script_read_solexa
 677 {
 678     # Martin A. Hansen, March 2008.
 679
 680     # Read Solexa sequence reads from file.
 681
 682     my ( $in,        # handle to in stream
 683          $out,       # handle to out stream
 684          $options,   # options hash
 685        ) = @_;
 686
 687     # Returns nothing.
 688
 689     my ( $record, $file, $data_in, $entry, $num, @seqs, @scores, $i );
 690
 691     $options->{ "format" }  ||= "octal";
 692     $options->{ "quality" } ||= 20;
 693
 694     while ( $record = Maasha::Biopieces::get_record( $in ) ) {
 695         Maasha::Biopieces::put_record( $record, $out );
 696     }
 697
 698     $num = 1;
 699
 700     foreach $file ( @{ $options->{ "files" } } )
 701     {
 702         $data_in = Maasha::Common::read_open( $file );
 703
 704         if ( $options->{ "format" } eq "octal" )
 705         {
 706             while ( $entry = Maasha::Solexa::solexa_get_entry_octal( $data_in ) )
 707             {
 708                 $record = Maasha::Solexa::solexa2biopiece( $entry, $options->{ "quality" } );
 709
 710                 Maasha::Biopieces::put_record( $record, $out );
 711
 712                 goto NUM if $options->{ "num" } and $num == $options->{ "num" };
 713
 714                 $num++;
 715             }
 716         }
 717         else
 718         {
 719             while ( $entry = Maasha::Solexa::solexa_get_entry_decimal( $data_in ) )
 720             {
 721                 $record = Maasha::Solexa::solexa2biopiece( $entry, $options->{ "quality" } );
 722
 723                 Maasha::Biopieces::put_record( $record, $out );
 724
 725                 goto NUM if $options->{ "num" } and $num == $options->{ "num" };
 726
 727                 $num++;
 728             }
 729         }
 730
 731         close $data_in;
 732     }
 733
 734     NUM:
 735
 736     close $data_in if $data_in;
 737 }
 738
 739
 740 sub script_read_solid
 741 {
 742     # Martin A. Hansen, April 2008.
 743
 744     # Read Solid sequence from file.
 745
 746     my ( $in,        # handle to in stream
 747          $out,       # handle to out stream
 748          $options,   # options hash
 749        ) = @_;
 750
 751     # Returns nothing.
 752
 753     my ( $record, $file, $data_in, $line, $num, $seq_name, $seq_cs, $seq_qual, @scores, @seqs, $i );
 754
 755     $options->{ "quality" } ||= 15;
 756
 757     while ( $record = Maasha::Biopieces::get_record( $in ) ) {
 758         Maasha::Biopieces::put_record( $record, $out );
 759     }
 760
 761     $num = 1;
 762
 763     foreach $file ( @{ $options->{ "files" } } )
 764     {
 765         $data_in = Maasha::Common::read_open( $file );
 766
 767         while ( $line = <$data_in> )
 768         {
 769             chomp $line;
 770
 771             ( $seq_name, $seq_cs, $seq_qual ) = split /\t/, $line;
 772
 773             @scores = split /,/, $seq_qual;
 774             @seqs   = split //, Maasha::Solid::color_space2seq( $seq_cs );
 775
 776             for ( $i = 0; $i < @seqs; $i++ ) {
 777                 $seqs[ $i ] = lc $seqs[ $i ] if $scores[ $i ] < $options->{ "quality" };
 778             }
 779
 780             $record = {
 781                 REC_TYPE   => 'SOLID',
 782                 SEQ_NAME   => $seq_name,
 783                 SEQ_CS     => $seq_cs,
 784                 SEQ_QUAL   => join( ";", @scores ),
 785                 SEQ_LEN    => length $seq_cs,
 786                 SEQ        => join( "", @seqs ),
 787                 SCORE_MEAN => sprintf( "%.2f", Maasha::Calc::mean( \@scores ) ),
 788             };
 789
 790             Maasha::Biopieces::put_record( $record, $out );
 791
 792             goto NUM if $options->{ "num" } and $num == $options->{ "num" };
 793
 794             $num++;
 795         }
 796
 797         close $data_in;
 798     }
 799
 800     NUM:
 801
 802     close $data_in if $data_in;
 803 }
 804
 805
 806 sub script_read_mysql
 807 {
 808     # Martin A. Hansen, May 2008.
 809
 810     # Read a MySQL query into stream.
 811
 812     my ( $in,        # handle to in stream
 813          $out,       # handle to out stream
 814          $options,   # options hash
 815        ) = @_;
 816
 817     # Returns nothing.
 818
 819     my ( $record, $dbh, $results );
 820
 821     $options->{ "user" }     ||= Maasha::UCSC::ucsc_get_user();
 822     $options->{ "password" } ||= Maasha::UCSC::ucsc_get_password();
 823
 824     while ( $record = Maasha::Biopieces::get_record( $in ) ) {
 825         Maasha::Biopieces::put_record( $record, $out );
 826     }
 827
 828     $dbh = Maasha::SQL::connect( $options->{ "database" }, $options->{ "user" }, $options->{ "password" } );
 829
 830     $results = Maasha::SQL::query_hashref_list( $dbh, $options->{ "query" } );
 831
 832     Maasha::SQL::disconnect( $dbh );
 833
 834     map { Maasha::Biopieces::put_record( $_ ) } @{ $results };
 835 }
 836
 837
 838 sub script_complexity_seq
 839 {
 840     # Martin A. Hansen, May 2008.
 841
 842     # Generates an index calculated as the most common di-residue over
 843     # the sequence length for all sequences in stream.
 844
 845     my ( $in,     # handle to in stream
 846          $out,    # handle to out stream
 847        ) = @_;
 848
 849     # Returns nothing.
 850
 851     my ( $record, $index );
 852
 853     while ( $record = Maasha::Biopieces::get_record( $in ) )
 854     {
 855         $record->{ "SEQ_COMPLEXITY" } = sprintf( "%.2f", Maasha::Seq::seq_complexity( $record->{ "SEQ" } ) ) if $record->{ "SEQ" };
 856
 857         Maasha::Biopieces::put_record( $record, $out );
 858     }
 859 }
 860
 861
 862 sub script_get_genome_align
 863 {
 864     # Martin A. Hansen, April 2008.
 865
 866     # Gets a subalignment from a multiple genome alignment.
 867
 868     my ( $in,        # handle to in stream
 869          $out,       # handle to out stream
 870          $options,   # options hash
 871        ) = @_;
 872
 873     # Returns nothing.
 874
 875     my ( $record, $maf_track, $align, $align_num, $beg, $end, $len, $entry );
 876
 877     $options->{ "strand" } ||= "+";
 878
 879     $align_num = 1;
 880
 881     $maf_track = Maasha::Config::maf_track( $options->{ "genome" } );
 882
 883     if ( $options->{ "chr" } and $options->{ "beg" } and ( $options->{ "end" } or $options->{ "len" } ) )
 884     {
 885         $beg = $options->{ "beg" } - 1;
 886
 887         if ( $options->{ "end" } ) {
 888             $end = $options->{ "end" };
 889         } elsif ( $options->{ "len" } ) {
 890             $end = $beg + $options->{ "len" };
 891         }
 892
 893         $align = Maasha::UCSC::maf_extract( $BP_TMP, $options->{ "genome" }, $maf_track, $options->{ "chr" }, $beg, $end, $options->{ "strand" } );
 894
 895         foreach $entry ( @{ $align } )
 896         {
 897             $entry->{ "CHR" }     = $record->{ "CHR" };
 898             $entry->{ "CHR_BEG" } = $record->{ "CHR_BEG" };
 899             $entry->{ "CHR_END" } = $record->{ "CHR_END" };
 900             $entry->{ "STRAND" }  = $record->{ "STRAND" } || '+';
 901             $entry->{ "Q_ID" }    = $record->{ "Q_ID" };
 902             $entry->{ "SCORE" }   = $record->{ "SCORE" };
 903
 904             Maasha::Biopieces::put_record( $entry, $out );
 905         }
 906     }
 907
 908     while ( $record = Maasha::Biopieces::get_record( $in ) )
 909     {
 910         if ( $record->{ "REC_TYPE" } eq "BED" )
 911         {
 912             $align = Maasha::UCSC::maf_extract( $BP_TMP, $options->{ "genome" }, $maf_track, $record->{ "CHR" }, $record->{ "CHR_BEG" }, $record->{ "CHR_END" }, $record->{ "STRAND" } );
 913         }
 914         elsif ( $record->{ "REC_TYPE" } eq "VMATCH" )
 915         {
 916             $align = Maasha::UCSC::maf_extract( $BP_TMP, $options->{ "genome" }, $maf_track, $record->{ "S_ID" }, $record->{ "S_BEG" }, $record->{ "S_END" } + 1, $record->{ "STRAND" } );
 917         }
 918         elsif ( $record->{ "REC_TYPE" } eq "PSL" )
 919         {
 920             $align = Maasha::UCSC::maf_extract( $BP_TMP, $options->{ "genome" }, $maf_track, $record->{ "S_ID" }, $record->{ "S_BEG" }, $record->{ "S_END" }, $record->{ "STRAND" } );
 921         }
 922         elsif ( $record->{ "REC_TYPE" } eq "BLAST" )
 923         {
 924             $align = Maasha::UCSC::maf_extract( $BP_TMP, $options->{ "genome" }, $maf_track, $record->{ "S_ID" }, $record->{ "S_BEG" }, $record->{ "S_END" }, $record->{ "STRAND" } );
 925         }
 926
 927         foreach $entry ( @{ $align } )
 928         {
 929             $entry->{ "CHR" }     = $record->{ "CHR" };
 930             $entry->{ "CHR_BEG" } = $record->{ "CHR_BEG" };
 931             $entry->{ "CHR_END" } = $record->{ "CHR_END" };
 932             $entry->{ "STRAND" }  = $record->{ "STRAND" };
 933             $entry->{ "Q_ID" }    = $record->{ "Q_ID" };
 934             $entry->{ "SCORE" }   = $record->{ "SCORE" };
 935
 936             Maasha::Biopieces::put_record( $entry, $out );
 937         }
 938
 939         $align_num++;
 940     }
 941 }
 942
 943
 944 sub script_get_genome_phastcons
 945 {
 946     # Martin A. Hansen, February 2008.
 947
 948     # Get phastcons scores from genome intervals.
 949
 950     my ( $in,        # handle to in stream
 951          $out,       # handle to out stream
 952          $options,   # options hash
 953        ) = @_;
 954
 955     # Returns nothing.
 956
 957     my ( $phastcons_file, $phastcons_index, $index, $fh_phastcons, $scores, $record );
 958
 959     $options->{ "flank" } ||= 0;
 960
 961     $phastcons_file  = Maasha::Config::genome_phastcons( $options->{ "genome" } );
 962     $phastcons_index = Maasha::Config::genome_phastcons_index( $options->{ "genome" } );
 963
 964     $index           = Maasha::UCSC::fixedstep_index_retrieve( $phastcons_index );
 965     $fh_phastcons    = Maasha::Common::read_open( $phastcons_file );
 966
 967     if ( defined $options->{ "chr" } and defined $options->{ "beg" } and ( defined $options->{ "end" } or defined $options->{ "len" } ) )
 968     {
 969         $options->{ "beg" } -= 1;   # request is 1-based
 970         $options->{ "end" } -= 1;   # request is 1-based
 971
 972         if ( $options->{ "len" } ) {
 973             $options->{ "end" } = $options->{ "beg" } + $options->{ "len" } - 1;
 974         }
 975
 976         $scores = Maasha::UCSC::fixedstep_index_lookup( $index, $fh_phastcons, $options->{ "chr" }, $options->{ "beg" }, $options->{ "end" }, $options->{ "flank" } );
 977
 978         $record->{ "CHR" }       = $options->{ "chr" };
 979         $record->{ "CHR_BEG" }   = $options->{ "beg" } - $options->{ "flank" };
 980         $record->{ "CHR_END" }   = $options->{ "end" } + $options->{ "flank" };
 981
 982         $record->{ "PHASTCONS" }   = join ",", @{ $scores };
 983         $record->{ "PHAST_COUNT" } = scalar @{ $scores };  # DEBUG
 984
 985         Maasha::Biopieces::put_record( $record, $out );
 986     }
 987
 988     while ( $record = Maasha::Biopieces::get_record( $in ) )
 989     {
 990         if ( $record->{ "REC_TYPE" } eq "BED" )
 991         {
 992             $scores = Maasha::UCSC::fixedstep_index_lookup( $index, $fh_phastcons, $record->{ "CHR" }, $record->{ "CHR_BEG" }, $record->{ "CHR_END" }, $options->{ "flank" } );
 993         }
 994         elsif ( $record->{ "REC_TYPE" } eq "PSL" )
 995         {
 996             $scores = Maasha::UCSC::fixedstep_index_lookup( $index, $fh_phastcons, $record->{ "S_ID" }, $record->{ "S_BEG" }, $record->{ "S_END" }, $options->{ "flank" } );
 997         }
 998         elsif ( $record->{ "REC_TYPE" } eq "BLAST" )
 999         {
1000             $scores = Maasha::UCSC::fixedstep_index_lookup( $index, $fh_phastcons, $record->{ "S_ID" }, $record->{ "S_BEG" }, $record->{ "S_END" }, $options->{ "flank" } );
1001         }
1002
1003         $record->{ "PHASTCONS" } = join ",", @{ $scores } if @{ $scores };
1004 #        $record->{ "PHAST_COUNT" } = @{ $scores } if @{ $scores };  # DEBUG
1005
1006         Maasha::Biopieces::put_record( $record, $out );
1007     }
1008
1009     close $fh_phastcons if $fh_phastcons;
1010 }
1011
1012
1013 sub script_soap_seq
1014 {
1015     # Martin A. Hansen, July 2008.
1016
1017     # soap sequences in stream against a given file or genome.
1018
1019     my ( $in,        # handle to in stream
1020          $out,       # handle to out stream
1021          $options,   # options hash
1022        ) = @_;
1023
1024     # Returns nothing.
1025
1026     my ( $genome, $tmp_in, $tmp_out, $fh_in, $fh_out, $record, $line, @fields, $entry, $count, $args );
1027
1028     $options->{ "seed_size" }  ||= 10;
1029     $options->{ "mismatches" } ||= 2;
1030     $options->{ "gap_size" }   ||= 0;
1031     $options->{ "cpus" }       ||= 1;
1032
1033     if ( $options->{ "genome" } ) {
1034         $options->{ "in_file" } = "$ENV{ 'BP_DATA' }/genomes/$options->{ 'genome' }/fasta/$options->{ 'genome' }.fna";
1035     }
1036
1037     $tmp_in  = "$BP_TMP/soap_query.seq";
1038     $tmp_out = "$BP_TMP/soap.result";
1039
1040     $fh_out = Maasha::Common::write_open( $tmp_in );
1041
1042     $count = 0;
1043
1044     while ( $record = Maasha::Biopieces::get_record( $in ) )
1045     {
1046         if ( $entry = Maasha::Fasta::biopiece2fasta( $record ) )
1047         {
1048             Maasha::Fasta::put_entry( $entry, $fh_out );
1049
1050             $count++;
1051         }
1052
1053         Maasha::Biopieces::put_record( $record, $out );
1054     }
1055
1056     close $fh_out;
1057
1058     if ( $count > 0 )
1059     {
1060         $args = join( " ",
1061             "-s $options->{ 'seed_size' }",
1062             "-r 2",
1063             "-a $tmp_in",
1064             "-v $options->{ 'mismatches' }",
1065             "-g $options->{ 'gap_size' }",
1066             "-p $options->{ 'cpus' }",
1067             "-d $options->{ 'in_file' }",
1068             "-o $tmp_out",
1069         );
1070
1071         $args .= " > /dev/null 2>&1" if not $options->{ 'verbose' };
1072
1073         Maasha::Common::run( "soap", $args, 1 );
1074
1075         unlink $tmp_in;
1076
1077         $fh_out = Maasha::Common::read_open( $tmp_out );
1078
1079         undef $record;
1080
1081         while ( $line = <$fh_out> )
1082         {
1083             chomp $line;
1084
1085             @fields = split /\t/, $line;
1086
1087             $record->{ "REC_TYPE" }   = "SOAP";
1088             $record->{ "Q_ID" }       = $fields[ 0 ];
1089             $record->{ "SCORE" }      = $fields[ 3 ];
1090             $record->{ "STRAND" }     = $fields[ 6 ];
1091             $record->{ "S_ID" }       = $fields[ 7 ];
1092             $record->{ "S_BEG" }      = $fields[ 8 ] - 1; # soap is 1-based
1093             $record->{ "S_END" }      = $fields[ 8 ] + $fields[ 5 ] - 2;
1094
1095             Maasha::Biopieces::put_record( $record, $out );
1096         }
1097
1098         close $fh_out;
1099     }
1100
1101     unlink $tmp_out;
1102 }
1103
1104
1105 sub script_remove_mysql_tables
1106 {
1107     # Martin A. Hansen, November 2008.
1108
1109     # Remove MySQL tables from values in stream.
1110
1111     my ( $in,        # handle to in stream
1112          $out,       # handle to out stream
1113          $options,   # options hash
1114        ) = @_;
1115
1116     # Returns nothing.
1117
1118     my ( $record, %table_hash, $dbh, $table );
1119
1120     $options->{ "user" }     ||= Maasha::UCSC::ucsc_get_user();
1121     $options->{ "password" } ||= Maasha::UCSC::ucsc_get_password();
1122
1123     map { $table_hash{ $_ } = 1 } @{ $options->{ 'tables' } };
1124
1125     while ( $record = Maasha::Biopieces::get_record( $in ) )
1126     {
1127         map { $table_hash{ $record->{ $_ } } = 1 } @{ $options->{ 'keys' } };
1128
1129         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ 'no_stream' };
1130     }
1131
1132     $dbh = Maasha::SQL::connect( $options->{ "database" }, $options->{ "user" }, $options->{ "password" } );
1133
1134     foreach $table ( sort keys %table_hash )
1135     {
1136         if ( Maasha::SQL::table_exists( $dbh, $table ) )
1137         {
1138             print STDERR qq(Removing table "$table" from database "$options->{ 'database' }" ... ) if $options->{ 'verbose' };
1139             Maasha::SQL::delete_table( $dbh, $table );
1140             print STDERR "done.\n" if $options->{ 'verbose' };
1141         }
1142         else
1143         {
1144             print STDERR qq(WARNING: table "$table" not found in database "$options->{ 'database' }\n");
1145         }
1146     }
1147
1148     Maasha::SQL::disconnect( $dbh );
1149 }
1150
1151
1152 sub script_remove_ucsc_tracks
1153 {
1154     # Martin A. Hansen, November 2008.
1155
1156     # Remove track from MySQL tables and config file.
1157
1158     my ( $in,        # handle to in stream
1159          $out,       # handle to out stream
1160          $options,   # options hash
1161        ) = @_;
1162
1163     # Returns nothing.
1164
1165     my ( $record, %track_hash, $fh_in, $fh_out, $track, @tracks, @new_tracks, $dbh );
1166
1167     $options->{ 'user' }        ||= Maasha::UCSC::ucsc_get_user();
1168     $options->{ 'password' }    ||= Maasha::UCSC::ucsc_get_password();
1169     $options->{ 'config_file' } ||= "$ENV{ 'HOME' }/ucsc/my_tracks.ra";
1170
1171     map { $track_hash{ $_ } = 1 } @{ $options->{ 'tracks' } };
1172
1173     while ( $record = Maasha::Biopieces::get_record( $in ) )
1174     {
1175         map { $track_hash{ $record->{ $_ } } = 1 } @{ $options->{ 'keys' } };
1176
1177         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ 'no_stream' };
1178     }
1179
1180     $fh_in = Maasha::Filesys::file_read_open( $options->{ 'config_file' } );
1181
1182     while ( $track = Maasha::UCSC::ucsc_config_entry_get( $fh_in ) ) {
1183         push @tracks, $track;
1184     }
1185
1186     close $fh_in;
1187
1188     foreach $track ( @tracks )
1189     {
1190         if ( $track->{ 'database' } eq $options->{ 'database' } and exists $track_hash{ $track->{ 'track' } } ) {
1191             print STDERR qq(Removing track "$track->{ 'track' }" from config file.\n) if $options->{ 'verbose' };
1192         } else {
1193             push @new_tracks, $track;
1194         }
1195     }
1196
1197     rename "$options->{ 'config_file' }", "$options->{ 'config_file' }~";
1198
1199     $fh_out = Maasha::Filesys::file_write_open( $options->{ 'config_file' } );
1200
1201     map { Maasha::UCSC::ucsc_config_entry_put( $_, $fh_out ) } @new_tracks;
1202
1203     close $fh_out;
1204
1205     # ---- locate track in database ----
1206
1207     $dbh = Maasha::SQL::connect( $options->{ "database" }, $options->{ "user" }, $options->{ "password" } );
1208
1209     foreach $track ( sort keys %track_hash )
1210     {
1211         if ( Maasha::SQL::table_exists( $dbh, $track ) )
1212         {
1213             print STDERR qq(Removing table "$track" from database "$options->{ 'database' }" ... ) if $options->{ 'verbose' };
1214             Maasha::SQL::delete_table( $dbh, $track );
1215             print STDERR "done.\n" if $options->{ 'verbose' };
1216         }
1217         else
1218         {
1219             print STDERR qq(WARNING: table "$track" not found in database "$options->{ 'database' }\n");
1220         }
1221     }
1222
1223     Maasha::SQL::disconnect( $dbh );
1224
1225     Maasha::Common::run( "ucscMakeTracks.pl", "-b > /dev/null 2>&1" );
1226 }
1227
1228
1229 sub script_upload_to_ucsc
1230 {
1231     # Martin A. Hansen, August 2007.
1232
1233     # Calculate the mean of values of given keys.
1234
1235     # This routine has developed into the most ugly hack. Do something!
1236
1237     my ( $in,        # handle to in stream
1238          $out,       # handle to out stream
1239          $options,   # options hash
1240        ) = @_;
1241
1242     # Returns nothing.
1243
1244     my ( $record, $file, $wib_file, $wig_file, $wib_dir, $fh_out, $i, $first, $format, $type, $columns, $append, $entry );
1245
1246     $options->{ "short_label" } ||= $options->{ 'table' };
1247     $options->{ "long_label" }  ||= $options->{ 'table' };
1248     $options->{ "group" }       ||= $ENV{ "LOGNAME" };
1249     $options->{ "priority" }    ||= 1;
1250     $options->{ "visibility" }  ||= "pack";
1251     $options->{ "color" }       ||= join( ",", int( rand( 255 ) ), int( rand( 255 ) ), int( rand( 255 ) ) );
1252     $options->{ "chunk_size" }  ||= 10_000_000_000;    # Due to 32-bit UCSC compilation really large tables cannot be loaded in one go.
1253
1254     $file   = "$BP_TMP/ucsc_upload.tmp";
1255     $append = 0;
1256     $first  = 1;
1257     $i      = 0;
1258
1259     $fh_out = Maasha::Common::write_open( $file );
1260
1261     while ( $record = Maasha::Biopieces::get_record( $in ) )
1262     {
1263         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ "no_stream" };
1264
1265         if ( $record->{ "REC_TYPE" } eq "fixed_step" )
1266         {
1267             $format = "WIGGLE";
1268
1269             if ( $entry = Maasha::UCSC::Wiggle::biopiece2fixedstep( $record ) ) {
1270                 Maasha::UCSC::Wiggle::fixedstep_entry_put( $entry, $fh_out );
1271             }
1272         }
1273         elsif ( $record->{ "REC_TYPE" } eq "PSL" )
1274         {
1275             $format = "PSL";
1276
1277             Maasha::UCSC::psl_put_header( $fh_out ) if $first;
1278             Maasha::UCSC::psl_put_entry( $record, $fh_out );
1279
1280             $first = 0;
1281         }
1282         elsif ( $record->{ "REC_TYPE" } eq "BED" and $record->{ "SEC_STRUCT" } )
1283         {
1284             # chrom chromStart  chromEnd    name    score   strand  size    secStr  conf
1285
1286             $format  = "BED_SS";
1287
1288             print $fh_out join ( "\t",
1289                 $record->{ "CHR" },
1290                 $record->{ "CHR_BEG" },
1291                 $record->{ "CHR_END" } + 1,
1292                 $record->{ "Q_ID" },
1293                 $record->{ "SCORE" },
1294                 $record->{ "STRAND" },
1295                 $record->{ "SIZE" },
1296                 $record->{ "SEC_STRUCT" },
1297                 $record->{ "CONF" },
1298             ), "\n";
1299         }
1300         elsif ( $record->{ "REC_TYPE" } eq "BED" )
1301         {
1302             $format  = "BED";
1303             $columns = $record->{ "BED_COLS" };
1304
1305             if ( $entry = Maasha::UCSC::BED::biopiece2bed( $record, $columns ) ) {
1306                 Maasha::UCSC::BED::bed_entry_put( $entry, $fh_out, $columns, $options->{ 'check' } );
1307             }
1308         }
1309         elsif ( $record->{ "REC_TYPE" } eq "PATSCAN" and $record->{ "CHR" } )
1310         {
1311             $format  = "BED";
1312             $columns = 6;
1313
1314             if ( $entry = Maasha::UCSC::BED::biopiece2bed( $record, $columns ) ) {
1315                 Maasha::UCSC::BED::bed_entry_put( $entry, $fh_out, $columns, $options->{ 'check' } );
1316             }
1317         }
1318         elsif ( $record->{ "REC_TYPE" } eq "BLAST" and $record->{ "S_ID" } =~ /^chr/ )
1319         {
1320             $format  = "BED";
1321             $columns = 6;
1322
1323             $record->{ "SCORE" } = $record->{ "BIT_SCORE" } * 1000;
1324
1325             if ( $entry = Maasha::UCSC::BED::biopiece2bed( $record, $columns ) ) {
1326                 Maasha::UCSC::BED::bed_entry_put( $entry, $fh_out, $columns, $options->{ 'check' } );
1327             }
1328         }
1329         elsif ( $record->{ "REC_TYPE" } eq "VMATCH" and $record->{ "S_ID" } =~ /^chr/i )
1330         {
1331             $format  = "BED";
1332             $columns = 6;
1333
1334             if ( $entry = Maasha::UCSC::BED::biopiece2bed( $record, $columns ) ) {
1335                 Maasha::UCSC::BED::bed_entry_put( $entry, $fh_out, $columns, $options->{ 'check' } );
1336             }
1337         }
1338
1339         if ( $i == $options->{ "chunk_size" } )
1340         {
1341             close $fh_out;
1342
1343             if ( $format eq "BED" ) {
1344                 Maasha::UCSC::bed_upload_to_ucsc( $BP_TMP, $file, $options, $append );
1345             } elsif ( $format eq "PSL" ) {
1346                 Maasha::UCSC::psl_upload_to_ucsc( $file, $options, $append );
1347             }
1348
1349             unlink $file;
1350
1351             $first = 1;
1352
1353             $append = 1;
1354
1355             $fh_out = Maasha::Common::write_open( $file );
1356         }
1357
1358         $i++;
1359     }
1360
1361     close $fh_out;
1362
1363     if ( exists $options->{ "database" } and $options->{ "table" } )
1364     {
1365         if ( $format eq "BED" )
1366         {
1367             $type = "bed $columns";
1368
1369             Maasha::UCSC::bed_upload_to_ucsc( $BP_TMP, $file, $options, $append );
1370         }
1371         elsif ( $format eq "BED_SS" )
1372         {
1373             $type = "type bed 6 +";
1374
1375             Maasha::UCSC::bed_upload_to_ucsc( $BP_TMP, $file, $options, $append );
1376         }
1377         elsif ( $format eq "PSL" )
1378         {
1379             $type = "psl";
1380
1381             Maasha::UCSC::psl_upload_to_ucsc( $file, $options, $append );
1382         }
1383         elsif ( $format eq "WIGGLE" )
1384         {
1385             $options->{ "visibility" } = "full";
1386
1387             $wig_file = "$options->{ 'table' }.wig";
1388             $wib_file = "$options->{ 'table' }.wib";
1389
1390             $wib_dir  = "$ENV{ 'HOME' }/ucsc/wib";
1391
1392             Maasha::Common::dir_create_if_not_exists( $wib_dir );
1393
1394             if ( $options->{ 'verbose' } ) {
1395                 `cd $BP_TMP && wigEncode $file $wig_file $wib_file`;
1396             } else {
1397                 `cd $BP_TMP && wigEncode $file $wig_file $wib_file > /dev/null 2>&1`;
1398             }
1399
1400             Maasha::Common::run( "mv", "$BP_TMP/$wib_file $wib_dir" );
1401
1402             unlink $file;
1403
1404             $file = $wig_file;
1405
1406             $type = "wig 0";
1407
1408             Maasha::UCSC::wiggle_upload_to_ucsc( $BP_TMP, $wib_dir, $file, $options );
1409         }
1410
1411         unlink $file;
1412
1413         Maasha::UCSC::ucsc_update_config( $options, $type );
1414     }
1415 }
1416
1417
1418 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1419
1420 1;
1421
1422 __END__
1423