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::Seq;
  44 use Maasha::Calc;
  45 use Maasha::UCSC;
  46 use Maasha::UCSC::BED;
  47 use Maasha::UCSC::Wiggle;
  48 use Maasha::NCBI;
  49 use Maasha::GFF;
  50 use Maasha::TwoBit;
  51 use Maasha::Solid;
  52 use Maasha::Solexa;
  53 use Maasha::SQL;
  54 use Maasha::Gwiki;
  55
  56 use vars qw( @ISA @EXPORT_OK );
  57
  58 require Exporter;
  59
  60 @ISA = qw( Exporter );
  61
  62 use constant {
  63     SEQ_NAME => 0,
  64     SEQ      => 1,
  65 };
  66
  67
  68 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> GLOBALS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  69
  70
  71 my ( $script, $BP_TMP );
  72
  73 $script = Maasha::Common::get_scriptname();
  74 $BP_TMP = Maasha::Common::get_tmpdir();
  75
  76
  77 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> RUN SCRIPT <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  78
  79
  80 run_script( $script );
  81
  82
  83 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> SUBROUTINES <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  84
  85
  86 sub run_script
  87 {
  88     # Martin A. Hansen, August 2007.
  89
  90     # Run a specific script.
  91
  92     my ( $script,   # script name
  93        ) = @_;
  94
  95     # Returns nothing.
  96
  97     my ( $t0, $t1, $options, $in, $out );
  98
  99     Maasha::Biopieces::log_biopiece();
 100
 101     $t0 = gettimeofday();
 102
 103     $options = get_options( $script );
 104
 105     $options->{ "SCRIPT" } = $script;
 106
 107     $script = "print_usage" if ( -t STDIN and keys %{ $options } <= 1 or $options->{ 'help' } );
 108
 109     $in  = Maasha::Biopieces::read_stream( $options->{ "stream_in" } );
 110     $out = Maasha::Biopieces::write_stream( $options->{ "stream_out" } );
 111
 112     if    ( $script eq "print_usage" )              { script_print_usage(               $in, $out, $options ) }
 113     elsif ( $script eq "complexity_seq" )           { script_complexity_seq(            $in, $out, $options ) }
 114     elsif ( $script eq "get_genome_align" )         { script_get_genome_align(          $in, $out, $options ) }
 115     elsif ( $script eq "get_genome_phastcons" )     { script_get_genome_phastcons(      $in, $out, $options ) }
 116     elsif ( $script eq "remove_mysql_tables" )      { script_remove_mysql_tables(       $in, $out, $options ) }
 117     elsif ( $script eq "remove_ucsc_tracks" )       { script_remove_ucsc_tracks(        $in, $out, $options ) }
 118
 119     close $in if defined $in;
 120     close $out;
 121
 122     $t1 = gettimeofday();
 123
 124     print STDERR "Program: $script" . ( " " x ( 25 - length( $script ) ) ) . sprintf( "Run time: %.4f\n", ( $t1 - $t0 ) ) if $options->{ 'verbose' };
 125 }
 126
 127
 128 sub get_options
 129 {
 130     # Martin A. Hansen, February 2008.
 131
 132     # Gets options from commandline and checks these vigerously.
 133
 134     my ( $script,     # name of script
 135        ) = @_;
 136
 137     # Returns hash
 138
 139     my ( %options, @options, $opt, @genomes, $real );
 140
 141     if ( $script eq "print_usage" )
 142     {
 143         @options = qw(
 144             data_in|i=s
 145         );
 146     }
 147     elsif ( $script eq "get_genome_align" )
 148     {
 149         @options = qw(
 150             genome|g=s
 151             chr|c=s
 152             beg|b=s
 153             end|e=s
 154             len|l=s
 155             strand|s=s
 156         );
 157     }
 158     elsif ( $script eq "get_genome_phastcons" )
 159     {
 160         @options = qw(
 161             genome|g=s
 162             chr|c=s
 163             beg|b=s
 164             end|e=s
 165             len|l=s
 166             flank|f=s
 167         );
 168     }
 169     elsif ( $script eq "remove_mysql_tables" )
 170     {
 171         @options = qw(
 172             database|d=s
 173             tables|t=s
 174             keys|k=s
 175             user|u=s
 176             password|p=s
 177             no_stream|x
 178         );
 179     }
 180     elsif ( $script eq "remove_ucsc_tracks" )
 181     {
 182         @options = qw(
 183             database|d=s
 184             tracks|t=s
 185             keys|k=s
 186             config_file|c=s
 187             user|u=s
 188             password|p=s
 189             no_stream|x
 190         );
 191     }
 192
 193     push @options, qw(
 194         stream_in|I=s
 195         stream_out|O=s
 196         verbose|v
 197         help|?
 198     );
 199
 200 #    print STDERR Dumper( \@options );
 201
 202     GetOptions(
 203         \%options,
 204         @options,
 205     );
 206
 207 #    print STDERR Dumper( \%options );
 208
 209     if ( -t STDIN && scalar( keys %options ) == 0 or $options{ "help" } ) {
 210         return wantarray ? %options : \%options;
 211     }
 212
 213     $options{ "cols" }      = [ split ",", $options{ "cols" } ]      if defined $options{ "cols" };
 214     $options{ "keys" }      = [ split ",", $options{ "keys" } ]      if defined $options{ "keys" };
 215     $options{ "no_keys" }   = [ split ",", $options{ "no_keys" } ]   if defined $options{ "no_keys" };
 216     $options{ "save_keys" } = [ split ",", $options{ "save_keys" } ] if defined $options{ "save_keys" };
 217     $options{ "quals" }     = [ split ",", $options{ "quals" } ]     if defined $options{ "quals" };
 218     $options{ "feats" }     = [ split ",", $options{ "feats" } ]     if defined $options{ "feats" };
 219     $options{ "frames" }    = [ split ",", $options{ "frames" } ]    if defined $options{ "frames" };
 220     $options{ "samples" }   = [ split ",", $options{ "samples" } ]   if defined $options{ "samples" };
 221     $options{ "tables" }    = [ split ",", $options{ "tables" } ]    if defined $options{ "tables" };
 222     $options{ "tracks" }    = [ split ",", $options{ "tracks" } ]    if defined $options{ "tracks" };
 223
 224     # ---- check arguments ----
 225
 226     if ( $options{ 'data_in' } )
 227     {
 228         $options{ "files" } = Maasha::Biopieces::getopt_files( $options{ 'data_in' } );
 229
 230         Maasha::Common::error( qq(Argument to --data_in must be a valid file or fileglob expression) ) if scalar @{ $options{ "files" } } == 0;
 231     }
 232
 233     map { Maasha::Common::error( qq(Argument to --cols must be a whole numbers - not "$_") ) if $_ !~ /^\d+$/ } @{ $options{ "cols" } } if $options{ "cols" };
 234
 235     # print STDERR Dumper( \%options );
 236
 237     $real = "beg|end|word_size|wrap|len|prefix_length|mismatches|offset|num|skip|cpus|window_size|step_size";
 238
 239     foreach $opt ( keys %options )
 240     {
 241         if ( $opt =~ /stream_in|pattern_in|exact_in/ and not -f $options{ $opt } )
 242         {
 243             Maasha::Common::error( qq(Argument to --$opt must be a valid file or fileglob expression - not "$options{ $opt }") );
 244         }
 245         elsif ( $opt =~ /$real/ and $options{ $opt } !~ /^\d+$/ )
 246         {
 247             Maasha::Common::error( qq(Argument to --$opt must be a whole number - not "$options{ $opt }") );
 248         }
 249         elsif ( $opt =~ /max_hits|max_hits|max_misses|dist|edit_dist|flank|gap|hamming_dist|priority/ and $options{ $opt } !~ /^-?\d+$/ )
 250         {
 251             Maasha::Common::error( qq(Argument to --$opt must be an integer - not "$options{ $opt }") );
 252         }
 253         elsif ( $opt =~ /identity|threshold/ and $options{ $opt } !~ /^-?(?:\d+(?:\.\d*)?|\.\d+)$/ )
 254         {
 255             Maasha::Common::error( qq(Argument to --$opt must be a decimal number - not "$options{ $opt }") );
 256         }
 257         elsif ( $opt =~ /e_val/ and $options{ $opt } !~ /^([+-]?)(?=\d|\.\d)\d*(\.\d*)?([Ee]([+-]?\d+))?$/ )
 258         {
 259             Maasha::Common::error( qq(Argument to --$opt must be a float - not "$options{ $opt }") );
 260         }
 261         elsif ( $opt =~ /strand/ and $options{ $opt } !~ /^(\+|-)$/ )
 262         {
 263             Maasha::Common::error( qq(Argument to --$opt must be "+" or "-" - not "$options{ $opt }") );
 264         }
 265         elsif ( $opt eq "genome" )
 266         {
 267             @genomes = Maasha::Common::ls_dirs( "$ENV{ 'BP_DATA' }/genomes" );
 268             map { $_ =~ s/.*\/(.+)$/$1/ } @genomes;
 269
 270             if ( not grep { $_ =~ /^$options{ $opt }$/ } @genomes ) {
 271                 Maasha::Common::error( qq(Genome $options{ $opt } not found in "$ENV{ 'BP_DATA' }/genomes/") );
 272             }
 273         }
 274         elsif ( $opt eq "terminal" and not $options{ $opt } =~ /^(svg|post|dumb|x11)/ )
 275         {
 276             Maasha::Common::error( qq(Bad --$opt argument "$options{ $opt }") );
 277         }
 278         elsif ( $opt eq "table" and $options{ $opt } =~ /(-|\.)/ )
 279         {
 280             Maasha::Common::error( qq(Character '$1' is not allowed in table name: $options{ $opt }) );
 281         }
 282         elsif ( $opt eq "merge" and $options{ $opt } !~ /^(AandB|AorB|BorA|AnotB|BnotA)$/ )
 283         {
 284             Maasha::Common::error( qq(Argument to --$opt must be AandB, AorB, BorA, AnotB, or BnotA - not "$options{ $opt }") );
 285         }
 286     }
 287
 288     Maasha::Common::error( qq(no --database specified) )                if $script eq "remove_ucsc_tracks"  and not $options{ "database" };
 289     Maasha::Common::error( qq(no --genome specified) )                  if $script =~ /get_genome_align|get_genome_phastcons/ and not $options{ "genome" };
 290
 291     return wantarray ? %options : \%options;
 292 }
 293
 294
 295 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> SCRIPTS  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 296
 297
 298 sub script_print_usage
 299 {
 300     # Martin A. Hansen, January 2008.
 301
 302     # Retrieves usage information from file and
 303     # prints this nicely formatted.
 304
 305     my ( $in,        # handle to in stream
 306          $out,       # handle to out stream
 307          $options,   # options hash
 308        ) = @_;
 309
 310     # Returns nothing.
 311
 312     my ( $file, $wiki, $lines );
 313
 314     if ( $options->{ 'data_in' } ) {
 315         $file = $options->{ 'data_in' };
 316     } else {
 317         $file = join "", $ENV{ 'BP_DIR' }, "/bp_usage/", $options->{ 'SCRIPT' }, ".wiki";
 318     }
 319
 320     $wiki = Maasha::Gwiki::gwiki_read( $file );
 321
 322     ( $wiki->[ 2 ], $wiki->[ 3 ], $wiki->[ 0 ], $wiki->[ 1 ] ) = ( $wiki->[ 0 ], $wiki->[ 1 ], $wiki->[ 2 ], $wiki->[ 3 ] );
 323
 324     if ( not $options->{ "help" } ) {
 325         @{ $wiki } = grep { $_->[ 0 ]->{ 'SECTION' } =~ /Biopiece|summary|Usage|Options|Help/ } @{ $wiki };
 326     }
 327
 328     $lines = Maasha::Gwiki::gwiki2ascii( $wiki );
 329
 330     print STDERR "$_\n" foreach @{ $lines };
 331
 332     exit;
 333 }
 334
 335
 336 sub script_complexity_seq
 337 {
 338     # Martin A. Hansen, May 2008.
 339
 340     # Generates an index calculated as the most common di-residue over
 341     # the sequence length for all sequences in stream.
 342
 343     my ( $in,     # handle to in stream
 344          $out,    # handle to out stream
 345        ) = @_;
 346
 347     # Returns nothing.
 348
 349     my ( $record, $index );
 350
 351     while ( $record = Maasha::Biopieces::get_record( $in ) )
 352     {
 353         $record->{ "SEQ_COMPLEXITY" } = sprintf( "%.2f", Maasha::Seq::seq_complexity( $record->{ "SEQ" } ) ) if $record->{ "SEQ" };
 354
 355         Maasha::Biopieces::put_record( $record, $out );
 356     }
 357 }
 358
 359
 360 sub script_get_genome_align
 361 {
 362     # Martin A. Hansen, April 2008.
 363
 364     # Gets a subalignment from a multiple genome alignment.
 365
 366     my ( $in,        # handle to in stream
 367          $out,       # handle to out stream
 368          $options,   # options hash
 369        ) = @_;
 370
 371     # Returns nothing.
 372
 373     my ( $record, $maf_track, $align, $align_num, $beg, $end, $len, $entry );
 374
 375     $options->{ "strand" } ||= "+";
 376
 377     $align_num = 1;
 378
 379     $maf_track = Maasha::Config::maf_track( $options->{ "genome" } );
 380
 381     if ( $options->{ "chr" } and $options->{ "beg" } and ( $options->{ "end" } or $options->{ "len" } ) )
 382     {
 383         $beg = $options->{ "beg" } - 1;
 384
 385         if ( $options->{ "end" } ) {
 386             $end = $options->{ "end" };
 387         } elsif ( $options->{ "len" } ) {
 388             $end = $beg + $options->{ "len" };
 389         }
 390
 391         $align = Maasha::UCSC::maf_extract( $BP_TMP, $options->{ "genome" }, $maf_track, $options->{ "chr" }, $beg, $end, $options->{ "strand" } );
 392
 393         foreach $entry ( @{ $align } )
 394         {
 395             $entry->{ "CHR" }     = $record->{ "CHR" };
 396             $entry->{ "CHR_BEG" } = $record->{ "CHR_BEG" };
 397             $entry->{ "CHR_END" } = $record->{ "CHR_END" };
 398             $entry->{ "STRAND" }  = $record->{ "STRAND" } || '+';
 399             $entry->{ "Q_ID" }    = $record->{ "Q_ID" };
 400             $entry->{ "SCORE" }   = $record->{ "SCORE" };
 401
 402             Maasha::Biopieces::put_record( $entry, $out );
 403         }
 404     }
 405
 406     while ( $record = Maasha::Biopieces::get_record( $in ) )
 407     {
 408         if ( $record->{ "REC_TYPE" } eq "BED" )
 409         {
 410             $align = Maasha::UCSC::maf_extract( $BP_TMP, $options->{ "genome" }, $maf_track, $record->{ "CHR" }, $record->{ "CHR_BEG" }, $record->{ "CHR_END" }, $record->{ "STRAND" } );
 411         }
 412         elsif ( $record->{ "REC_TYPE" } eq "VMATCH" )
 413         {
 414             $align = Maasha::UCSC::maf_extract( $BP_TMP, $options->{ "genome" }, $maf_track, $record->{ "S_ID" }, $record->{ "S_BEG" }, $record->{ "S_END" } + 1, $record->{ "STRAND" } );
 415         }
 416         elsif ( $record->{ "REC_TYPE" } eq "PSL" )
 417         {
 418             $align = Maasha::UCSC::maf_extract( $BP_TMP, $options->{ "genome" }, $maf_track, $record->{ "S_ID" }, $record->{ "S_BEG" }, $record->{ "S_END" }, $record->{ "STRAND" } );
 419         }
 420         elsif ( $record->{ "REC_TYPE" } eq "BLAST" )
 421         {
 422             $align = Maasha::UCSC::maf_extract( $BP_TMP, $options->{ "genome" }, $maf_track, $record->{ "S_ID" }, $record->{ "S_BEG" }, $record->{ "S_END" }, $record->{ "STRAND" } );
 423         }
 424
 425         foreach $entry ( @{ $align } )
 426         {
 427             $entry->{ "CHR" }     = $record->{ "CHR" };
 428             $entry->{ "CHR_BEG" } = $record->{ "CHR_BEG" };
 429             $entry->{ "CHR_END" } = $record->{ "CHR_END" };
 430             $entry->{ "STRAND" }  = $record->{ "STRAND" };
 431             $entry->{ "Q_ID" }    = $record->{ "Q_ID" };
 432             $entry->{ "SCORE" }   = $record->{ "SCORE" };
 433
 434             Maasha::Biopieces::put_record( $entry, $out );
 435         }
 436
 437         $align_num++;
 438     }
 439 }
 440
 441
 442 sub script_get_genome_phastcons
 443 {
 444     # Martin A. Hansen, February 2008.
 445
 446     # Get phastcons scores from genome intervals.
 447
 448     my ( $in,        # handle to in stream
 449          $out,       # handle to out stream
 450          $options,   # options hash
 451        ) = @_;
 452
 453     # Returns nothing.
 454
 455     my ( $phastcons_file, $phastcons_index, $index, $fh_phastcons, $scores, $record );
 456
 457     $options->{ "flank" } ||= 0;
 458
 459     $phastcons_file  = Maasha::Config::genome_phastcons( $options->{ "genome" } );
 460     $phastcons_index = Maasha::Config::genome_phastcons_index( $options->{ "genome" } );
 461
 462     $index           = Maasha::UCSC::fixedstep_index_retrieve( $phastcons_index );
 463     $fh_phastcons    = Maasha::Common::read_open( $phastcons_file );
 464
 465     if ( defined $options->{ "chr" } and defined $options->{ "beg" } and ( defined $options->{ "end" } or defined $options->{ "len" } ) )
 466     {
 467         $options->{ "beg" } -= 1;   # request is 1-based
 468         $options->{ "end" } -= 1;   # request is 1-based
 469
 470         if ( $options->{ "len" } ) {
 471             $options->{ "end" } = $options->{ "beg" } + $options->{ "len" } - 1;
 472         }
 473
 474         $scores = Maasha::UCSC::fixedstep_index_lookup( $index, $fh_phastcons, $options->{ "chr" }, $options->{ "beg" }, $options->{ "end" }, $options->{ "flank" } );
 475
 476         $record->{ "CHR" }       = $options->{ "chr" };
 477         $record->{ "CHR_BEG" }   = $options->{ "beg" } - $options->{ "flank" };
 478         $record->{ "CHR_END" }   = $options->{ "end" } + $options->{ "flank" };
 479
 480         $record->{ "PHASTCONS" }   = join ",", @{ $scores };
 481         $record->{ "PHAST_COUNT" } = scalar @{ $scores };  # DEBUG
 482
 483         Maasha::Biopieces::put_record( $record, $out );
 484     }
 485
 486     while ( $record = Maasha::Biopieces::get_record( $in ) )
 487     {
 488         if ( $record->{ "REC_TYPE" } eq "BED" )
 489         {
 490             $scores = Maasha::UCSC::fixedstep_index_lookup( $index, $fh_phastcons, $record->{ "CHR" }, $record->{ "CHR_BEG" }, $record->{ "CHR_END" }, $options->{ "flank" } );
 491         }
 492         elsif ( $record->{ "REC_TYPE" } eq "PSL" )
 493         {
 494             $scores = Maasha::UCSC::fixedstep_index_lookup( $index, $fh_phastcons, $record->{ "S_ID" }, $record->{ "S_BEG" }, $record->{ "S_END" }, $options->{ "flank" } );
 495         }
 496         elsif ( $record->{ "REC_TYPE" } eq "BLAST" )
 497         {
 498             $scores = Maasha::UCSC::fixedstep_index_lookup( $index, $fh_phastcons, $record->{ "S_ID" }, $record->{ "S_BEG" }, $record->{ "S_END" }, $options->{ "flank" } );
 499         }
 500
 501         $record->{ "PHASTCONS" } = join ",", @{ $scores } if @{ $scores };
 502 #        $record->{ "PHAST_COUNT" } = @{ $scores } if @{ $scores };  # DEBUG
 503
 504         Maasha::Biopieces::put_record( $record, $out );
 505     }
 506
 507     close $fh_phastcons if $fh_phastcons;
 508 }
 509
 510
 511 sub script_remove_mysql_tables
 512 {
 513     # Martin A. Hansen, November 2008.
 514
 515     # Remove MySQL tables from values in stream.
 516
 517     my ( $in,        # handle to in stream
 518          $out,       # handle to out stream
 519          $options,   # options hash
 520        ) = @_;
 521
 522     # Returns nothing.
 523
 524     my ( $record, %table_hash, $dbh, $table );
 525
 526     $options->{ "user" }     ||= Maasha::UCSC::ucsc_get_user();
 527     $options->{ "password" } ||= Maasha::UCSC::ucsc_get_password();
 528
 529     map { $table_hash{ $_ } = 1 } @{ $options->{ 'tables' } };
 530
 531     while ( $record = Maasha::Biopieces::get_record( $in ) )
 532     {
 533         map { $table_hash{ $record->{ $_ } } = 1 } @{ $options->{ 'keys' } };
 534
 535         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ 'no_stream' };
 536     }
 537
 538     $dbh = Maasha::SQL::connect( $options->{ "database" }, $options->{ "user" }, $options->{ "password" } );
 539
 540     foreach $table ( sort keys %table_hash )
 541     {
 542         if ( Maasha::SQL::table_exists( $dbh, $table ) )
 543         {
 544             print STDERR qq(Removing table "$table" from database "$options->{ 'database' }" ... ) if $options->{ 'verbose' };
 545             Maasha::SQL::delete_table( $dbh, $table );
 546             print STDERR "done.\n" if $options->{ 'verbose' };
 547         }
 548         else
 549         {
 550             print STDERR qq(WARNING: table "$table" not found in database "$options->{ 'database' }\n");
 551         }
 552     }
 553
 554     Maasha::SQL::disconnect( $dbh );
 555 }
 556
 557
 558 sub script_remove_ucsc_tracks
 559 {
 560     # Martin A. Hansen, November 2008.
 561
 562     # Remove track from MySQL tables and config file.
 563
 564     my ( $in,        # handle to in stream
 565          $out,       # handle to out stream
 566          $options,   # options hash
 567        ) = @_;
 568
 569     # Returns nothing.
 570
 571     my ( $record, %track_hash, $fh_in, $fh_out, $track, @tracks, @new_tracks, $dbh );
 572
 573     $options->{ 'user' }        ||= Maasha::UCSC::ucsc_get_user();
 574     $options->{ 'password' }    ||= Maasha::UCSC::ucsc_get_password();
 575     $options->{ 'config_file' } ||= "$ENV{ 'HOME' }/ucsc/my_tracks.ra";
 576
 577     map { $track_hash{ $_ } = 1 } @{ $options->{ 'tracks' } };
 578
 579     while ( $record = Maasha::Biopieces::get_record( $in ) )
 580     {
 581         map { $track_hash{ $record->{ $_ } } = 1 } @{ $options->{ 'keys' } };
 582
 583         Maasha::Biopieces::put_record( $record, $out ) if not $options->{ 'no_stream' };
 584     }
 585
 586     $fh_in = Maasha::Filesys::file_read_open( $options->{ 'config_file' } );
 587
 588     while ( $track = Maasha::UCSC::ucsc_config_entry_get( $fh_in ) ) {
 589         push @tracks, $track;
 590     }
 591
 592     close $fh_in;
 593
 594     foreach $track ( @tracks )
 595     {
 596         if ( $track->{ 'database' } eq $options->{ 'database' } and exists $track_hash{ $track->{ 'track' } } ) {
 597             print STDERR qq(Removing track "$track->{ 'track' }" from config file.\n) if $options->{ 'verbose' };
 598         } else {
 599             push @new_tracks, $track;
 600         }
 601     }
 602
 603     rename "$options->{ 'config_file' }", "$options->{ 'config_file' }~";
 604
 605     $fh_out = Maasha::Filesys::file_write_open( $options->{ 'config_file' } );
 606
 607     map { Maasha::UCSC::ucsc_config_entry_put( $_, $fh_out ) } @new_tracks;
 608
 609     close $fh_out;
 610
 611     # ---- locate track in database ----
 612
 613     $dbh = Maasha::SQL::connect( $options->{ "database" }, $options->{ "user" }, $options->{ "password" } );
 614
 615     foreach $track ( sort keys %track_hash )
 616     {
 617         if ( Maasha::SQL::table_exists( $dbh, $track ) )
 618         {
 619             print STDERR qq(Removing table "$track" from database "$options->{ 'database' }" ... ) if $options->{ 'verbose' };
 620             Maasha::SQL::delete_table( $dbh, $track );
 621             print STDERR "done.\n" if $options->{ 'verbose' };
 622         }
 623         else
 624         {
 625             print STDERR qq(WARNING: table "$track" not found in database "$options->{ 'database' }\n");
 626         }
 627     }
 628
 629     Maasha::SQL::disconnect( $dbh );
 630
 631     Maasha::Common::run( "ucscMakeTracks.pl", "-b > /dev/null 2>&1" );
 632 }
 633
 634
 635 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 636
 637 1;
 638
 639 __END__
 640