code_perl/Maasha/UCSC/Wiggle.pm

   1 package Maasha::UCSC::Wiggle;
   2
   3 # Copyright (C) 2007-2008 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::Calc;
  40
  41 use vars qw( @ISA @EXPORT_OK );
  42
  43 require Exporter;
  44
  45 @ISA = qw( Exporter );
  46
  47 @EXPORT_OK = qw(
  48 );
  49
  50
  51 use constant {
  52     BITS        => 32,            # Number of bits in an integer
  53     SEQ_MAX     => 200_000_000,   # Maximum sequence size
  54     chrom       => 0,             # BED field names
  55     chromStart  => 1,
  56     chromEnd    => 2,
  57     name        => 3,
  58     score       => 4,
  59     strand      => 5,
  60     thickStart  => 6,
  61     thickEnd    => 7,
  62     itemRgb     => 8,
  63     blockCount  => 9,
  64     blockSizes  => 10,
  65     blockStarts => 11,
  66 };
  67
  68
  69 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> FixedStep format <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  70
  71
  72 sub fixedstep_entry_get
  73 {
  74     # Martin A. Hansen, December 2007.
  75
  76     # Given a file handle to a PhastCons file get the
  77     # next entry which is all the lines after a "fixedStep"
  78     # line and until the next "fixedStep" line or EOF.
  79
  80     my ( $fh,   # filehandle
  81        ) = @_;
  82
  83     # Returns a list of lines
  84
  85     my ( $entry, @lines );
  86
  87     local $/ = "\nfixedStep ";
  88
  89     $entry = <$fh>;
  90
  91     return if not $entry;
  92
  93     chomp $entry;
  94
  95     @lines = split "\n", $entry;
  96
  97     return if @lines == 0;
  98
  99     $lines[ 0 ] =~ s/fixedStep?\s*//;
 100     $lines[ 0 ] = 'fixedStep ' . $lines[ 0 ];
 101
 102     return wantarray ? @lines : \@lines;
 103 }
 104
 105
 106 sub fixedstep_entry_put
 107 {
 108     # Martin A. Hansen, April 2008.
 109
 110     # Outputs a block of fixedStep values.
 111     # Used for outputting wiggle data.
 112
 113     my ( $entry,    # fixedStep entry
 114          $fh,       # filehandle - OPTIONAL
 115        ) = @_;
 116
 117     # Returns nothing.
 118
 119     $fh ||= \*STDOUT;
 120
 121     print $fh join( "\n", @{ $entry } ), "\n";
 122 }
 123
 124
 125 sub fixedstep2biopiece
 126 {
 127     # Martin A. Hansen, November 2008.
 128
 129     # Converts a fixedstep entry to a Biopiece record.
 130
 131     my ( $entry,   # fixedstep entry arrayref
 132        ) = @_;
 133
 134     # Returns a hashref
 135
 136     my ( $head, %record );
 137
 138     $head = shift @{ $entry };
 139
 140     if ( $head =~ /^fixedStep chrom=([^ ]+) start=(\d+) step=(\d+)$/ )
 141     {
 142         $record{ "REC_TYPE" } = "fixed_step";
 143         $record{ "CHR" }      = $1;
 144         $record{ "CHR_BEG" }  = $2 - 1;  # fixedStep is 1-based
 145         $record{ "STEP" }     = $3;
 146         $record{ "VALS" }     = join ";", @{ $entry };
 147     }
 148     else
 149     {
 150         Maasha::Common::error( qq(could not convert fixedStep entry to Biopiece record) );
 151     }
 152
 153     return wantarray ? %record : \%record;
 154 }
 155
 156
 157 sub biopiece2fixedstep
 158 {
 159     # Martin A. Hansen, November 2008.
 160
 161     # Converts a Biopiece record to a fixedStep entry.
 162
 163     my ( $record,   # Biopiece record
 164        ) = @_;
 165
 166     # Returns a list
 167
 168     my ( @entry, $beg, $vals );
 169
 170     if ( exists $record->{ "REC_TYPE" } and $record->{ "REC_TYPE" } eq 'fixed_step' )
 171     {
 172         if ( exists $record->{ 'CHR' } and exists $record->{ 'CHR_BEG' } and exists $record->{ 'STEP' } )
 173         {
 174             $beg = $record->{ 'CHR_BEG' } + 1;   # fixedStep is 1-based
 175             push @entry, "fixedStep chrom=$record->{ 'CHR' } start=$beg step=$record->{ 'STEP' }";
 176
 177             $vals = $record->{ 'VALS' };
 178
 179             map { push @entry, $_ } split ";", $vals;
 180         }
 181         else
 182         {
 183             Maasha::Common::error( qq(could not convert Biopiece record to fixedStep) );
 184         }
 185     }
 186
 187     return wantarray ? @entry : \@entry;
 188 }
 189
 190
 191 sub fixedstep_calc
 192 {
 193     # Martin A. Hansen, November 2008.
 194
 195     # Calculates fixedstep entries for a given BED file.
 196     # Implemented using large Perl arrays to avoid sorting.
 197     # Returns path to the fixedstep file.
 198
 199     my ( $bed_file,   # path to BED file
 200          $chr,        # chromosome name
 201          $log10,      # flag indicating that the score should be in log10 - OPTIONAL
 202        ) = @_;
 203
 204     # Returns a string
 205
 206     my ( $fixedstep_file, $fh_in, $fh_out, $array, $beg, $end, @block, $i );
 207
 208     $fixedstep_file = "$bed_file.fixedstep";
 209
 210     $fh_in = Maasha::Filesys::file_read_open( $bed_file );
 211
 212     $array = fixedstep_calc_array( $fh_in );
 213
 214     close $fh_in;
 215
 216     $fh_out = Maasha::Filesys::file_write_open( $fixedstep_file );
 217
 218     $beg = 0;
 219
 220     while ( ( $beg, $end ) = fixedstep_scan( $array, $beg ) and $beg != -1 )
 221     {
 222         @block = "fixedStep chrom=$chr start=" . ( $beg + 1 ) . " step=1";
 223
 224         for ( $i = $beg; $i < $end; $i++ )
 225         {
 226             if ( $log10 ) {
 227                 push @block, sprintf "%.4f", Maasha::Calc::log10( $array->[ $i ] );
 228             } else {
 229                 push @block, $array->[ $i ];
 230             }
 231         }
 232
 233         fixedstep_entry_put( \@block, $fh_out );
 234
 235         $beg = $end + 1;
 236     }
 237
 238     close $fh_out;
 239
 240     undef $array;
 241
 242     return $fixedstep_file;
 243 }
 244
 245
 246 sub fixedstep_calc_array
 247 {
 248     # Martin A. Hansen, November 2008.
 249
 250     # Given a filehandle to a BED file for a single
 251     # chromosome fills an array with SCORE for the
 252     # begin/end positions of each BED entry.
 253
 254     my ( $fh_in,   # filehandle to BED file
 255        ) = @_;
 256
 257     # Returns a bytearray.
 258
 259     my ( $bed, $score, @begs, @lens, $i, @array );
 260
 261     while ( $bed = Maasha::UCSC::BED::bed_entry_get( $fh_in ) )
 262     {
 263         if ( defined $bed->[ score ] ) {
 264             $score = $bed->[ score ];
 265         } else {
 266             $score = 1;
 267         }
 268
 269         if ( $bed->[ blockCount ] and $bed->[ blockCount ] > 1 )
 270         {
 271             @begs = split ",", $bed->[ blockStarts ];
 272             @lens = split ",", $bed->[ blockSizes ];
 273
 274             for ( $i = 0; $i < @begs; $i++ ) {
 275                 map { $array[ $_ ] += $score } ( $bed->[ chromStart ] + $begs[ $i ] .. $bed->[ chromStart ] + $begs[ $i ] + $lens[ $i ] - 1 );
 276             }
 277         }
 278         else
 279         {
 280             map { $array[ $_ ] += $score } ( $bed->[ chromStart ] .. $bed->[ chromEnd ] - 1 );
 281         }
 282     }
 283
 284     return wantarray ? @array : \@array;
 285 }
 286
 287
 288 sub fixedstep_scan
 289 {
 290     # Martin A. Hansen, November 2008.
 291
 292     # Given a fixedstep array locates the next block of
 293     # non-zero elements from a given begin position.
 294     # A [ begin end ] tuple of the block posittion is returned
 295     # if a block was found otherwise and empty tuple is returned.
 296
 297     my ( $array,   # array to scan
 298          $beg,     # position to start scanning from
 299        ) = @_;
 300
 301     # Returns an arrayref.
 302
 303     my ( $end );
 304
 305     while ( $beg < scalar @{ $array } and not $array->[ $beg ] ) { $beg++ }
 306
 307     $end = $beg;
 308
 309     while ( $array->[ $end ] ) { $end++ }
 310
 311     if ( $end > $beg ) {
 312         return wantarray ? ( $beg, $end ) : [ $beg, $end ];
 313     } else {
 314         return wantarray ? () : [];
 315     }
 316 }
 317
 318
 319 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 320
 321
 322 1;
 323
 324
 325 __END__