code_perl/Maasha/UCSC.pm

   1 package Maasha::UCSC;
   2
   3 # Copyright (C) 2007 Martin A. Hansen.
   4
   5 # This program is free software; you can redistribute it and/or
   6 # modify it under the terms of the GNU General Public License
   7 # as published by the Free Software Foundation; either version 2
   8 # of the License, or (at your option) any later version.
   9
  10 # This program is distributed in the hope that it will be useful,
  11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 # GNU General Public License for more details.
  14
  15 # You should have received a copy of the GNU General Public License
  16 # along with this program; if not, write to the Free Software
  17 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18
  19 # http://www.gnu.org/copyleft/gpl.html
  20
  21
  22 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DESCRIPTION <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  23
  24
  25 # Stuff for interacting with UCSC genome browser
  26
  27
  28 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  29
  30
  31 use warnings;
  32 use strict;
  33 use vars qw ( @ISA @EXPORT );
  34
  35 use Data::Dumper;
  36 use Maasha::Common;
  37 use Maasha::Filesys;
  38 use Maasha::Calc;
  39 use Maasha::Matrix;
  40
  41 use constant {
  42     CHR          => 0,
  43     CHR_BEG      => 1,
  44     CHR_END      => 2,
  45     Q_ID         => 3,
  46     SCORE        => 4,
  47     STRAND       => 5,
  48     THICK_BEG    => 6,
  49     THICK_END    => 7,
  50     COLOR        => 8,
  51     BLOCK_COUNT  => 9,
  52     BLOCK_LENS   => 10,
  53     Q_BEGS       => 11,
  54 };
  55
  56 @ISA = qw( Exporter );
  57
  58
  59 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> TRACK FILE <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  60
  61
  62 sub ucsc_config_entry_get
  63 {
  64     # Martin A. Hansen, November 2008.
  65
  66     # Given a filehandle to a UCSC Genome Browser
  67     # config file (.ra) get the next entry and
  68     # return as a hash. Entries are separated by blank
  69     # lines. # lines are skipped unless it is the lines:
  70     # # Track added ...
  71     # # Database ...
  72
  73     my ( $fh,   # file hanlde
  74        ) = @_;
  75
  76     # Returns a hashref.
  77
  78     my ( $line, %record );
  79
  80     while ( $line = <$fh> )
  81     {
  82         $line =~ tr/\n\r//d;
  83
  84         if ( $line =~ /Track added by 'upload_to_ucsc' (\S+) (\S+)/) {
  85              $record{ 'date' } = $1;
  86              $record{ 'time' } = $2;
  87         } elsif ( $line =~ /^# date: (.+)/ ) {
  88             $record{ 'date' } = $1;
  89         } elsif ( $line =~ /^# time: (.+)/ ) {
  90             $record{ 'time' } = $1;
  91         } elsif ( $line =~ /^# (?:database:|Database) (.+)/ ) {
  92             $record{ 'database' } = $1;
  93         } elsif ( $line =~ /^#/ ) {
  94             # skip
  95         } elsif ( $line =~ /(\S+)\s+(.+)/ ) {
  96             $record{ $1 } = $2;
  97         }
  98
  99         last if $line =~ /^$/ and exists $record{ "track" };
 100     }
 101
 102     return undef if not exists $record{ "track" };
 103
 104     return wantarray ? %record : \%record;
 105 }
 106
 107
 108 sub ucsc_config_entry_put
 109 {
 110     # Martin A. Hansen, November 2008.
 111
 112     # Outputs a Biopiece record (a hashref)
 113     # to a filehandle or STDOUT.
 114
 115     my ( $record,    # hashref
 116          $fh_out,    # file handle
 117        ) = @_;
 118
 119     # Returns nothing.
 120
 121     my ( $date, $time );
 122
 123     $fh_out ||= \*STDOUT;
 124
 125     ( $date, $time ) = split " ", Maasha::Common::time_stamp();
 126
 127     $record->{ "date" } ||= $date;
 128     $record->{ "time" } ||= $time;
 129
 130     print $fh_out "\n# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n\n";
 131
 132     map { print $fh_out "# $_: $record->{ $_ }\n" if exists $record->{ $_ } } qw( date time database );
 133     map { print $fh_out "$_ $record->{ $_ }\n" if exists $record->{ $_ } } qw( track
 134                                                                                name
 135                                                                                description
 136                                                                                itemRgb
 137                                                                                db
 138                                                                                offset
 139                                                                                url
 140                                                                                htmlUrl
 141                                                                                shortLabel
 142                                                                                longLabel
 143                                                                                group
 144                                                                                priority
 145                                                                                useScore
 146                                                                                visibility
 147                                                                                maxHeightPixels
 148                                                                                color
 149                                                                                mafTrack
 150                                                                                type );
 151 }
 152
 153
 154 sub ucsc_update_config
 155 {
 156     # Martin A. Hansen, September 2007.
 157
 158     # Update the /home/user/ucsc/my_tracks.ra file and executes makeCustomTracks.pl
 159
 160     my ( $options,   # hashref
 161          $type,      # track type
 162        ) = @_;
 163
 164     # Returns nothing.
 165
 166     my ( $file, $entry, %new_entry, $fh_in, $fh_out, $was_found );
 167
 168     $file = $ENV{ "HOME" } . "/ucsc/my_tracks.ra";
 169
 170     Maasha::Filesys::file_copy( $file, "$file~" );   # create a backup
 171
 172     %new_entry = (
 173         database   => $options->{ 'database' },
 174         track      => $options->{ 'table' },
 175         shortLabel => $options->{ 'short_label' },
 176         longLabel  => $options->{ 'long_label' },
 177         group      => $options->{ 'group' },
 178         priority   => $options->{ 'priority' },
 179         visibility => $options->{ 'visibility' },
 180         color      => $options->{ 'color' },
 181         type       => $type,
 182     );
 183
 184     $new_entry{ 'useScore' }        = 1             if $options->{ 'use_score' };
 185     $new_entry{ 'mafTrack' }        = "multiz17way" if $type eq "type bed 6 +";
 186     $new_entry{ 'maxHeightPixels' } = "50:50:11"    if $type eq "wig 0";
 187
 188     $fh_in  = Maasha::Filesys::file_read_open( "$file~" );
 189     $fh_out = Maasha::Filesys::file_write_open( $file );
 190
 191     while ( $entry = ucsc_config_entry_get( $fh_in ) )
 192     {
 193         if ( $entry->{ 'database' } eq $new_entry{ 'database' } and $entry->{ 'track' } eq $new_entry{ 'track' } )
 194         {
 195             ucsc_config_entry_put( \%new_entry, $fh_out );
 196
 197             $was_found = 1;
 198         }
 199         else
 200         {
 201             ucsc_config_entry_put( $entry, $fh_out );
 202         }
 203     }
 204
 205     ucsc_config_entry_put( \%new_entry, $fh_out ) if not $was_found;
 206
 207     close $fh_in;
 208     close $fh_out;
 209
 210     Maasha::Common::run( "ucscMakeTracks.pl", "-b > /dev/null 2>&1" );
 211 }
 212
 213
 214 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> PhastCons format <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 215
 216
 217 sub phastcons_index
 218 {
 219     # Martin A. Hansen, July 2008
 220
 221     # Create a fixedStep index for PhastCons data.
 222
 223     my ( $file,   # file to index
 224          $dir,    # dir with file
 225        ) = @_;
 226
 227     # Returns nothing.
 228
 229     my ( $index );
 230
 231     $index = fixedstep_index_create( "$dir/$file" );
 232
 233     fixedstep_index_store( "$dir/$file.index", $index );
 234 }
 235
 236
 237 sub phastcons_parse_entry
 238 {
 239     # Martin A. Hansen, December 2007.
 240
 241     # Given a PhastCons entry converts this to a
 242     # list of super blocks.
 243
 244 #    $options->{ "min" }       ||= 10;
 245 #    $options->{ "dist" }      ||= 25;
 246 #    $options->{ "threshold" } ||= 0.8;
 247 #    $options->{ "gap" }       ||= 5;
 248
 249     my ( $lines,   # list of lines
 250          $args,    # argument hash
 251        ) = @_;
 252
 253     # Returns
 254
 255     my ( $info, $chr, $beg, $step, $i, $c, $j, @blocks, @super_blocks, @entries, $super_block, $block, @lens, @begs );
 256
 257     $info = shift @{ $lines };
 258
 259     if ( $info =~ /^chrom=([^ ]+) start=(\d+) step=(\d+)$/ )
 260     {
 261         $chr  = $1;
 262         $beg  = $2;
 263         $step = $3;
 264
 265         die qq(ERROR: step size $step != 1 -> problem!\n) if $step != 1; # in an ideal world should would be fixed ...
 266     }
 267
 268     $i = 0;
 269
 270     while ( $i < @{ $lines } )
 271     {
 272         if ( $lines->[ $i ] >= $args->{ "threshold" } )
 273         {
 274             $c = $i + 1;
 275
 276             while ( $c < @{ $lines } )
 277             {
 278                 if ( $lines->[ $c ] < $args->{ "threshold" } )
 279                 {
 280                     $j = $c + 1;
 281
 282                     while ( $j < @{ $lines } and $lines->[ $j ] < $args->{ "threshold" } ) {
 283                         $j++;
 284                     }
 285
 286                     if ( $j - $c > $args->{ "gap" } )
 287                     {
 288                         if ( $c - $i >= $args->{ "min" } )
 289                         {
 290                             push @blocks, {
 291                                 CHR     => $chr,
 292                                 CHR_BEG => $beg + $i - 1,
 293                                 CHR_END => $beg + $c - 2,
 294                                 CHR_LEN => $c - $i,
 295                             };
 296                         }
 297
 298                         $i = $j;
 299
 300                         last;
 301                     }
 302
 303                     $c = $j
 304                 }
 305                 else
 306                 {
 307                     $c++;
 308                 }
 309             }
 310
 311             if ( $c - $i >= $args->{ "min" } )
 312             {
 313                 push @blocks, {
 314                     CHR     => $chr,
 315                     CHR_BEG => $beg + $i - 1,
 316                     CHR_END => $beg + $c - 2,
 317                     CHR_LEN => $c - $i,
 318                 };
 319             }
 320
 321             $i = $c;
 322         }
 323         else
 324         {
 325             $i++;
 326         }
 327     }
 328
 329     $i = 0;
 330
 331     while ( $i < @blocks )
 332     {
 333         $c = $i + 1;
 334
 335         while ( $c < @blocks and $blocks[ $c ]->{ "CHR_BEG" } - $blocks[ $c - 1 ]->{ "CHR_END" } <= $args->{ "dist" } )
 336         {
 337             $c++;
 338         }
 339
 340         push @super_blocks, [ @blocks[ $i .. $c - 1 ] ];
 341
 342         $i = $c;
 343     }
 344
 345     foreach $super_block ( @super_blocks )
 346     {
 347         foreach $block ( @{ $super_block } )
 348         {
 349             push @begs, $block->{ "CHR_BEG" } - $super_block->[ 0 ]->{ "CHR_BEG" };
 350             push @lens, $block->{ "CHR_LEN" } - 1;
 351         }
 352
 353         $lens[ -1 ]++;
 354
 355         push @entries, {
 356             CHR         => $super_block->[ 0 ]->{ "CHR" },
 357             CHR_BEG     => $super_block->[ 0 ]->{ "CHR_BEG" },
 358             CHR_END     => $super_block->[ -1 ]->{ "CHR_END" },
 359             Q_ID        => "Q_ID",
 360             SCORE       => 100,
 361             STRAND      => "+",
 362             THICK_BEG   => $super_block->[ 0 ]->{ "CHR_BEG" },
 363             THICK_END   => $super_block->[ -1 ]->{ "CHR_END" } + 1,
 364             COLOR       => 0,
 365             BLOCK_COUNT => scalar @{ $super_block },
 366             BLOCK_LENS  => join( ",", @lens ),
 367             Q_BEGS      => join( ",", @begs ),
 368         };
 369
 370         undef @begs;
 371         undef @lens;
 372     }
 373
 374     return wantarray ? @entries : \@entries;
 375 }
 376
 377
 378 sub phastcons_normalize
 379 {
 380     # Martin A. Hansen, January 2008.
 381
 382     # Normalizes a list of lists with PhastCons scores,
 383     # in such a way that each list contains the same number
 384     # of PhastCons scores.
 385
 386     my ( $AoA,    # AoA with PhastCons scores
 387        ) = @_;
 388
 389     # Returns AoA.
 390
 391     my ( $list, $max, $min, $mean, $diff );
 392
 393     $min = 99999999;
 394     $max = 0;
 395
 396     foreach $list ( @{ $AoA } )
 397     {
 398         $min = scalar @{ $list } if scalar @{ $list } < $min;
 399         $max = scalar @{ $list } if scalar @{ $list } > $max;
 400     }
 401
 402     $mean = int( ( $min + $max ) / 2 );
 403
 404 #    print STDERR "min->$min   max->$max   mean->$mean\n";
 405
 406     foreach $list ( @{ $AoA } )
 407     {
 408         $diff = scalar @{ $list } - $mean;
 409
 410         phastcons_list_inflate( $list, abs( $diff ) ) if $diff < 0;
 411         phastcons_list_deflate( $list, $diff )        if $diff > 0;
 412     }
 413
 414     return wantarray ? @{ $AoA } : $AoA;
 415 }
 416
 417
 418 sub phastcons_list_inflate
 419 {
 420     # Martin A. Hansen, January 2008.
 421
 422     # Inflates a list with a given number of elements
 423     # in such a way that the extra elements are introduced
 424     # evenly over the entire length of the list. The value
 425     # of the extra elements is based on a mean of the
 426     # adjacent elements.
 427
 428     my ( $list,   # list of elements
 429          $diff,   # number of elements to introduce
 430        ) = @_;
 431
 432     # Returns nothing
 433
 434     my ( $len, $space, $i, $pos );
 435
 436     $len = scalar @{ $list };
 437
 438     $space = $len / $diff;
 439
 440     for ( $i = 0; $i < $diff; $i++ )
 441     {
 442         $pos = int( ( $space / 2 ) + $i * $space );
 443
 444         splice @{ $list }, $pos, 0, ( $list->[ $pos - 1 ] + $list->[ $pos + 1 ] ) / 2;
 445         # splice @{ $list }, $pos, 0, "X";
 446     }
 447
 448     die qq(ERROR: Bad inflate\n) if scalar @{ $list } != $len + $diff;
 449 }
 450
 451
 452 sub phastcons_list_deflate
 453 {
 454     # Martin A. Hansen, January 2008.
 455
 456     # Deflates a list by removing a given number of elements
 457     # evenly distributed over the entire list.
 458
 459     my ( $list,   # list of elements
 460          $diff,   # number of elements to remove
 461        ) = @_;
 462
 463     # Returns nothing
 464
 465     my ( $len, $space, $i, $pos );
 466
 467     $len = scalar @{ $list };
 468
 469     $space = ( $len - $diff ) / $diff;
 470
 471     for ( $i = 0; $i < $diff; $i++ )
 472     {
 473         $pos = int( ( $space / 2 ) + $i * $space );
 474
 475         splice @{ $list }, $pos, 1;
 476     }
 477
 478     die qq(ERROR: Dad deflate\n) if scalar @{ $list } != $len - $diff;
 479 }
 480
 481
 482 sub phastcons_mean
 483 {
 484     # Martin A. Hansen, January 2008.
 485
 486     # Given a normalized PhastCons matrix in an AoA,
 487     # calculate the mean for each column and return as a list.
 488
 489     my ( $AoA,    # AoA with normalized PhastCons scores
 490        ) = @_;
 491
 492     # Returns a list
 493
 494     my ( @list );
 495
 496     $AoA = Maasha::Matrix::matrix_flip( $AoA );
 497
 498     map { push @list, Maasha::Calc::mean( $_ ) } @{ $AoA };
 499
 500     return wantarray ? @list : \@list;
 501 }
 502
 503
 504 sub phastcons_median
 505 {
 506     # Martin A. Hansen, January 2008.
 507
 508     # Given a normalized PhastCons matrix in an AoA,
 509     # calculate the median for each column and return as a list.
 510
 511     my ( $AoA,    # AoA with normalized PhastCons scores
 512        ) = @_;
 513
 514     # Returns a list
 515
 516     my ( @list );
 517
 518     $AoA = Maasha::Matrix::matrix_flip( $AoA );
 519
 520     map { push @list, Maasha::Calc::median( $_ ) } @{ $AoA };
 521
 522     return wantarray ? @list : \@list;
 523 }
 524
 525
 526 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> MULTIPLE ALIGNMENT FILES <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 527
 528
 529 sub maf_extract
 530 {
 531     # Martin A. Hansen, April 2008.
 532
 533     # Executes mafFrag to extract a subalignment from a multiz track
 534     # in the UCSC genome browser database.
 535
 536     my ( $tmp_dir,    # temporary directory
 537          $database,   # genome database
 538          $table,      # table with the multiz track
 539          $chr,        # chromosome
 540          $beg,        # begin position
 541          $end,        # end position
 542          $strand,     # strand
 543        ) = @_;
 544
 545     # Returns a list of record
 546
 547     my ( $tmp_file, $align );
 548
 549     $tmp_file = "$tmp_dir/maf_extract.maf";
 550
 551     Maasha::Common::run( "mafFrag", "$database $table $chr $beg $end $strand $tmp_file" );
 552
 553     $align = maf_parse( $tmp_file );
 554
 555     unlink $tmp_file;
 556
 557     return wantarray ? @{ $align } : $align;
 558 }
 559
 560
 561 sub maf_parse
 562 {
 563     # Martin A. Hansen, April 2008.
 564
 565
 566     my ( $path,   # full path to MAF file
 567        ) = @_;
 568
 569     # Returns a list of record.
 570
 571     my ( $fh, $line, @fields, @align );
 572
 573     $fh = Maasha::Filesys::file_read_open( $path );
 574
 575     while ( $line = <$fh> )
 576     {
 577         chomp $line;
 578
 579         if ( $line =~ /^s/ )
 580         {
 581             @fields = split / /, $line;
 582
 583             push @align, {
 584                 SEQ_NAME  => $fields[ 1 ],
 585                 SEQ       => $fields[ -1 ],
 586                 ALIGN     => 1,
 587                 ALIGN_LEN => length $fields[ -1 ],
 588             }
 589         }
 590     }
 591
 592     close $fh;
 593
 594     return wantarray ? @align : \@align;
 595 }
 596
 597
 598 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> MySQL CONF <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 599
 600
 601 sub ucsc_get_user
 602 {
 603     # Martin A. Hansen, May 2008
 604
 605     # Fetches the MySQL database user name from the
 606     # .hg.conf file in the users home directory.
 607
 608     # Returns a string.
 609
 610     my ( $file, $fh, $line, $user );
 611
 612     $file = "$ENV{ 'HOME' }/.hg.conf";
 613
 614     return if not -f $file;
 615
 616     $fh = Maasha::Filesys::file_read_open( $file );
 617
 618     while ( $line = <$fh> )
 619     {
 620         chomp $line;
 621
 622         if ( $line =~ /^db\.user=(.+)/ )
 623         {
 624             $user = $1;
 625
 626             last;
 627         }
 628     }
 629
 630     close $fh;
 631
 632     return $user;
 633 }
 634
 635
 636 sub ucsc_get_password
 637 {
 638     # Martin A. Hansen, May 2008
 639
 640     # Fetches the MySQL database password from the
 641     # .hg.conf file in the users home directory.
 642
 643     # Returns a string.
 644
 645     my ( $file, $fh, $line, $password );
 646
 647     $file = "$ENV{ 'HOME' }/.hg.conf";
 648
 649     return if not -f $file;
 650
 651     $fh = Maasha::Filesys::file_read_open( $file );
 652
 653     while ( $line = <$fh> )
 654     {
 655         chomp $line;
 656
 657         if ( $line =~ /^db\.password=(.+)/ )
 658         {
 659             $password = $1;
 660
 661             last;
 662         }
 663     }
 664
 665     close $fh;
 666
 667     return $password;
 668 }
 669
 670
 671 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 672
 673
 674 __END__
 675