code_perl/Maasha/UCSC/Wiggle.pm

   1 package Maasha::UCSC::Wiggle;
   2
   3 # Copyright (C) 2007-2010 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 # Routines for manipulation of wiggle data.
  26
  27 # http://genome.ucsc.edu/goldenPath/help/wiggle.html
  28
  29
  30 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  31
  32
  33 use warnings;
  34 use strict;
  35
  36 use Data::Dumper;
  37 use Maasha::Common;
  38 use Maasha::Filesys;
  39 use Maasha::Matrix;
  40 use Maasha::Calc;
  41
  42 use vars qw( @ISA @EXPORT_OK );
  43
  44 require Exporter;
  45
  46 @ISA = qw( Exporter );
  47
  48 @EXPORT_OK = qw(
  49 );
  50
  51
  52 use constant {
  53     BITS            => 32,            # Number of bits in an integer
  54     SEQ_MAX         => 200_000_000,   # Maximum sequence size
  55     chrom           => 0,             # BED field names
  56     chromStart      => 1,
  57     chromEnd        => 2,
  58     name            => 3,
  59     score           => 4,
  60     strand          => 5,
  61     thickStart      => 6,
  62     thickEnd        => 7,
  63     itemRgb         => 8,
  64     blockCount      => 9,
  65     blockSizes      => 10,
  66     blockStarts     => 11,
  67 };
  68
  69
  70 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> FixedStep format <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  71
  72
  73 sub fixedstep_entry_get
  74 {
  75     # Martin A. Hansen, December 2007.
  76
  77     # Given a file handle to a PhastCons file get the
  78     # next entry which is all the lines after a "fixedStep"
  79     # line and until the next "fixedStep" line or EOF.
  80
  81     my ( $fh,   # filehandle
  82        ) = @_;
  83
  84     # Returns a list of lines
  85
  86     my ( $entry, @lines );
  87
  88     local $/ = "\nfixedStep ";
  89
  90     $entry = <$fh>;
  91
  92     return if not $entry;
  93
  94     chomp $entry;
  95
  96     @lines = split "\n", $entry;
  97
  98     return if @lines == 0;
  99
 100     $lines[ 0 ] =~ s/fixedStep?\s*//;
 101     $lines[ 0 ] = 'fixedStep ' . $lines[ 0 ];
 102
 103     return wantarray ? @lines : \@lines;
 104 }
 105
 106
 107 sub fixedstep_entry_put
 108 {
 109     # Martin A. Hansen, April 2008.
 110
 111     # Outputs a block of fixedStep values.
 112     # Used for outputting wiggle data.
 113
 114     my ( $entry,    # fixedStep entry
 115          $fh,       # filehandle - OPTIONAL
 116        ) = @_;
 117
 118     # Returns nothing.
 119
 120     $fh ||= \*STDOUT;
 121
 122     print $fh join( "\n", @{ $entry } ), "\n";
 123 }
 124
 125
 126 sub fixedstep2biopiece
 127 {
 128     # Martin A. Hansen, November 2008.
 129
 130     # Converts a fixedstep entry to a Biopiece record.
 131
 132     my ( $entry,   # fixedstep entry arrayref
 133        ) = @_;
 134
 135     # Returns a hashref
 136
 137     my ( $head, %record );
 138
 139     $head = shift @{ $entry };
 140
 141     if ( $head =~ /^fixedStep chrom=([^ ]+) start=(\d+) step=(\d+)$/ )
 142     {
 143         $record{ "REC_TYPE" } = "fixed_step";
 144         $record{ "CHR" }      = $1;
 145         $record{ "CHR_BEG" }  = $2 - 1;  # fixedStep is 1-based
 146         $record{ "STEP" }     = $3;
 147         $record{ "VALS" }     = join ";", @{ $entry };
 148     }
 149     else
 150     {
 151         Maasha::Common::error( qq(could not convert fixedStep entry to Biopiece record) );
 152     }
 153
 154     return wantarray ? %record : \%record;
 155 }
 156
 157
 158 sub biopiece2fixedstep
 159 {
 160     # Martin A. Hansen, November 2008.
 161
 162     # Converts a Biopiece record to a fixedStep entry.
 163
 164     my ( $record,   # Biopiece record
 165        ) = @_;
 166
 167     # Returns a list
 168
 169     my ( @entry, $beg, $vals );
 170
 171     if ( exists $record->{ "REC_TYPE" } and $record->{ "REC_TYPE" } eq 'fixed_step' )
 172     {
 173         if ( exists $record->{ 'CHR' } and exists $record->{ 'CHR_BEG' } and exists $record->{ 'STEP' } )
 174         {
 175             $beg = $record->{ 'CHR_BEG' } + 1;   # fixedStep is 1-based
 176             push @entry, "fixedStep chrom=$record->{ 'CHR' } start=$beg step=$record->{ 'STEP' }";
 177
 178             $vals = $record->{ 'VALS' };
 179
 180             map { push @entry, $_ } split ";", $vals;
 181         }
 182         else
 183         {
 184             Maasha::Common::error( qq(could not convert Biopiece record to fixedStep) );
 185         }
 186     }
 187
 188     unless ( @entry ) {
 189         return
 190     } else {
 191         return wantarray ? @entry : \@entry;
 192     }
 193 }
 194
 195
 196 sub fixedstep_calc
 197 {
 198     # Martin A. Hansen, November 2008.
 199
 200     # Calculates fixedstep entries for a given BED file.
 201     # Implemented using large Perl arrays to avoid sorting.
 202     # Returns path to the fixedstep file.
 203
 204     my ( $bed_file,   # path to BED file
 205          $chr,        # chromosome name
 206          $log10,      # flag indicating that the score should be in log10 - OPTIONAL
 207        ) = @_;
 208
 209     # Returns a string
 210
 211     my ( $fixedstep_file, $fh_in, $fh_out, $array, $beg, $end, @block, $i );
 212
 213     $fixedstep_file = "$bed_file.fixedstep";
 214
 215     $fh_in = Maasha::Filesys::file_read_open( $bed_file );
 216
 217     $array = fixedstep_calc_array( $fh_in );
 218
 219     close $fh_in;
 220
 221     $fh_out = Maasha::Filesys::file_write_open( $fixedstep_file );
 222
 223     $beg = 0;
 224
 225     while ( ( $beg, $end ) = fixedstep_scan( $array, $beg ) and $beg != -1 )
 226     {
 227         @block = "fixedStep chrom=$chr start=" . ( $beg + 1 ) . " step=1";
 228
 229         for ( $i = $beg; $i < $end; $i++ )
 230         {
 231             if ( $log10 ) {
 232                 push @block, sprintf "%.4f", Maasha::Calc::log10( $array->[ $i ] );
 233             } else {
 234                 push @block, $array->[ $i ];
 235             }
 236         }
 237
 238         fixedstep_entry_put( \@block, $fh_out );
 239
 240         $beg = $end + 1;
 241     }
 242
 243     close $fh_out;
 244
 245     undef $array;
 246
 247     return $fixedstep_file;
 248 }
 249
 250
 251 sub fixedstep_calc_array
 252 {
 253     # Martin A. Hansen, November 2008.
 254
 255     # Given a filehandle to a BED file for a single
 256     # chromosome fills an array with SCORE for the
 257     # begin/end positions of each BED entry.
 258
 259     my ( $fh_in,   # filehandle to BED file
 260        ) = @_;
 261
 262     # Returns a bytearray.
 263
 264     my ( $bed, $score, @begs, @lens, $i, @array );
 265
 266     while ( $bed = Maasha::UCSC::BED::bed_entry_get( $fh_in ) )
 267     {
 268         if ( defined $bed->[ score ] ) {
 269             $score = $bed->[ score ];
 270         } else {
 271             $score = 1;
 272         }
 273
 274         if ( $bed->[ blockCount ] and $bed->[ blockCount ] > 1 )
 275         {
 276             @begs = split ",", $bed->[ blockStarts ];
 277             @lens = split ",", $bed->[ blockSizes ];
 278
 279             for ( $i = 0; $i < @begs; $i++ ) {
 280                 map { $array[ $_ ] += $score } ( $bed->[ chromStart ] + $begs[ $i ] .. $bed->[ chromStart ] + $begs[ $i ] + $lens[ $i ] - 1 );
 281             }
 282         }
 283         else
 284         {
 285             map { $array[ $_ ] += $score } ( $bed->[ chromStart ] .. $bed->[ chromEnd ] - 1 );
 286         }
 287     }
 288
 289     return wantarray ? @array : \@array;
 290 }
 291
 292
 293 sub fixedstep_scan
 294 {
 295     # Martin A. Hansen, November 2008.
 296
 297     # Given a fixedstep array locates the next block of
 298     # non-zero elements from a given begin position.
 299     # A [ begin end ] tuple of the block posittion is returned
 300     # if a block was found otherwise and empty tuple is returned.
 301
 302     my ( $array,   # array to scan
 303          $beg,     # position to start scanning from
 304        ) = @_;
 305
 306     # Returns an arrayref.
 307
 308     my ( $end );
 309
 310     while ( $beg < scalar @{ $array } and not $array->[ $beg ] ) { $beg++ }
 311
 312     $end = $beg;
 313
 314     while ( $array->[ $end ] ) { $end++ }
 315
 316     if ( $end > $beg ) {
 317         return wantarray ? ( $beg, $end ) : [ $beg, $end ];
 318     } else {
 319         return wantarray ? () : [];
 320     }
 321 }
 322
 323
 324 sub fixedstep_index
 325 {
 326     # Martin A. Hansen, June 2010.
 327
 328     # Indexes a specified fixedstep file and saves the index
 329     # to a specified directory and index name.
 330
 331     my ( $path,   # fixedstep file to index
 332          $dir,    # direcotory to place index
 333          $name,   # name of index
 334        ) = @_;
 335
 336     # Returns nothing.
 337
 338     my ( $index );
 339
 340     $index = fixedstep_index_create( $path );
 341     fixedstep_index_store( "$dir/$name", $index );
 342 }
 343
 344
 345 sub fixedstep_index_create
 346 {
 347     # Martin A. Hansen, June 2010.
 348
 349     # Indexes a specified fixedStep file.
 350     # The index consists of a hash with chromosomes as keys,
 351     # and a list of { chr_beg, chr_end, index_beg, index_len }
 352     # hashes as values.
 353
 354     my ( $path,   # path to fixedStep file
 355        ) = @_;
 356
 357     # Returns a hashref
 358
 359     my ( $fh, $pos, $line, $index_beg, $index_len, $chr, $step, $chr_beg, $chr_end, $len, %index );
 360
 361     $fh = Maasha::Filesys::file_read_open( $path );
 362
 363     $index_beg = 0;
 364     $index_len = 0;
 365     $chr_beg   = 0;
 366     $chr_end   = 0;
 367     $pos       = 0;
 368
 369     while ( $line = <$fh> )
 370     {
 371         if ( $line =~ /^fixedStep chrom=([^ ]+) start=(\d+) step=(\d+)/ )
 372         {
 373             if ( $chr_end > $chr_beg )
 374             {
 375                 $index_len = $pos - $index_beg;
 376
 377                 push @{ $index{ $chr } }, {
 378                     CHR_BEG   => $chr_beg,
 379                     CHR_END   => $chr_end - 1,
 380                     INDEX_BEG => $index_beg,
 381                     INDEX_LEN => $index_len,
 382                 };
 383             }
 384
 385             $chr     = $1;
 386             $chr_beg = $2 - 1;  #  fixedStep files are 1-based
 387             $step    = $3;
 388             $chr_end = $chr_beg;
 389
 390             $index_beg = $pos + length $line;
 391         }
 392         else
 393         {
 394             $chr_end++;
 395         }
 396
 397         $pos += length $line;
 398     }
 399
 400     $index_len = $pos - $index_beg;
 401
 402     push @{ $index{ $chr } }, {
 403         CHR_BEG   => $chr_beg,
 404         CHR_END   => $chr_end - 1,
 405         INDEX_BEG => $index_beg,
 406         INDEX_LEN => $index_len,
 407     };
 408
 409     close $fh;
 410
 411     return wantarray ? %index : \%index;
 412 }
 413
 414
 415 sub fixedstep_index_store
 416 {
 417     # Martin A. Hansen, January 2008.
 418
 419     # Writes a fixedStep index to binary file.
 420
 421     my ( $path,   # full path to file
 422          $index,  # list with index
 423        ) = @_;
 424
 425     # returns nothing
 426
 427     Maasha::Filesys::file_store( $path, $index );
 428 }
 429
 430
 431 sub fixedstep_index_retrieve
 432 {
 433     # Martin A. Hansen, January 2008.
 434
 435     # Retrieves a fixedStep index from binary file.
 436
 437     my ( $path,   # full path to file
 438        ) = @_;
 439
 440     # returns list
 441
 442     my $index;
 443
 444     $index = Maasha::Filesys::file_retrieve( $path );
 445
 446     return wantarray ? %{ $index } : $index;
 447 }
 448
 449
 450 sub fixedstep_index_lookup
 451 {
 452     # Martin A. Hansen, January 2008.
 453
 454     # Locate fixedStep entries from index based on
 455     # chr, chr_beg and chr_end.
 456
 457     my ( $index,     # data structure
 458          $chr,       # chromosome
 459          $chr_beg,   # chromosome beg
 460          $chr_end,   # chromosome end
 461        ) = @_;
 462
 463     # Returns a list
 464
 465     my ( @subindex );
 466
 467     if ( exists $index->{ $chr } ) {
 468         @subindex = grep { $_->{ 'CHR_END' } >= $chr_beg and $_->{ 'CHR_BEG' } <= $chr_end } @{ $index->{ $chr } };
 469     } else {
 470         Maasha::Common::error( qq(Chromosome "$chr" not found in index) );
 471     }
 472
 473     return wantarray ? @subindex : \@subindex;
 474 }
 475
 476
 477 sub wiggle_upload_to_ucsc
 478 {
 479     # Martin A. Hansen, May 2008.
 480
 481     # Upload a wiggle file to the UCSC database.
 482
 483     my ( $tmp_dir,    # temporary directory
 484          $wib_dir,    # wib directory
 485          $wig_file,   # file to upload,
 486          $options,    # argument hashref
 487        ) = @_;
 488
 489     # Returns nothing.
 490
 491     my ( $args );
 492
 493 #    $args = join " ", "-tmpDir=$tmp_dir", "-pathPrefix=$wib_dir", $options->{ "database" }, $options->{ 'table' }, $wig_file;
 494
 495 #    Maasha::Common::run( "hgLoadWiggle", "$args > /dev/null 2>&1" );
 496
 497     if ( $options->{ 'verbose' } ) {
 498         `cd $tmp_dir && hgLoadWiggle -tmpDir=$tmp_dir -pathPrefix=$wib_dir $options->{ 'database' } $options->{ 'table' } $wig_file`;
 499     } else {
 500         `cd $tmp_dir && hgLoadWiggle -tmpDir=$tmp_dir -pathPrefix=$wib_dir $options->{ 'database' } $options->{ 'table' } $wig_file > /dev/null 2>&1`;
 501     }
 502 }
 503
 504
 505 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 506
 507
 508 1;
 509
 510
 511 __END__