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