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 strict;
  34
  35 use Data::Dumper;
  36 use Maasha::Common;
  37 use Maasha::Filesys;
  38 use Maasha::Calc;
  39
  40 use vars qw( @ISA @EXPORT_OK );
  41
  42 require Exporter;
  43
  44 @ISA = qw( Exporter );
  45
  46 @EXPORT_OK = qw(
  47 );
  48
  49
  50 use constant {
  51     BITS        => 32,            # Number of bits in an integer
  52     SEQ_MAX     => 200_000_000,   # Maximum sequence size
  53     chrom       => 0,             # BED field names
  54     chromStart  => 1,
  55     chromEnd    => 2,
  56     name        => 3,
  57     score       => 4,
  58     strand      => 5,
  59     thickStart  => 6,
  60     thickEnd    => 7,
  61     itemRgb     => 8,
  62     blockCount  => 9,
  63     blockSizes  => 10,
  64     blockStarts => 11,
  65 };
  66
  67
  68 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> FixedStep format <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  69
  70
  71 sub fixedstep_entry_get
  72 {
  73     # Martin A. Hansen, December 2007.
  74
  75     # Given a file handle to a PhastCons file get the
  76     # next entry which is all the lines after a "fixedStep"
  77     # line and until the next "fixedStep" line or EOF.
  78
  79     my ( $fh,   # filehandle
  80        ) = @_;
  81
  82     # Returns a list of lines
  83
  84     my ( $entry, @lines );
  85
  86     local $/ = "\nfixedStep ";
  87
  88     $entry = <$fh>;
  89
  90     return if not $entry;
  91
  92     chomp $entry;
  93
  94     @lines = split "\n", $entry;
  95
  96     return if @lines == 0;
  97
  98     $lines[ 0 ] =~ s/fixedStep?\s*//;
  99     $lines[ 0 ] = 'fixedStep ' . $lines[ 0 ];
 100
 101     return wantarray ? @lines : \@lines;
 102 }
 103
 104
 105 sub fixedstep_entry_put
 106 {
 107     # Martin A. Hansen, April 2008.
 108
 109     # Outputs a block of fixedStep values.
 110     # Used for outputting wiggle data.
 111
 112     my ( $entry,    # fixedStep entry
 113          $fh,       # filehandle - OPTIONAL
 114        ) = @_;
 115
 116     # Returns nothing.
 117
 118     $fh ||= \*STDOUT;
 119
 120     print $fh join( "\n", @{ $entry } ), "\n";
 121 }
 122
 123
 124 sub fixedstep2biopiece
 125 {
 126     # Martin A. Hansen, November 2008.
 127
 128     # Converts a fixedstep entry to a Biopiece record.
 129
 130     my ( $entry,   # fixedstep entry arrayref
 131        ) = @_;
 132
 133     # Returns a hashref
 134
 135     my ( $head, %record );
 136
 137     $head = shift @{ $entry };
 138
 139     if ( $head =~ /^fixedStep chrom=([^ ]+) start=(\d+) step=(\d+)$/ )
 140     {
 141         $record{ "REC_TYPE" } = "fixed_step";
 142         $record{ "CHR" }      = $1;
 143         $record{ "CHR_BEG" }  = $2 - 1;  # fixedStep is 1-based
 144         $record{ "STEP" }     = $3;
 145         $record{ "VALS" }     = join ";", @{ $entry };
 146     }
 147     else
 148     {
 149         Maasha::Common::error( qq(could not convert fixedStep entry to Biopiece record) );
 150     }
 151
 152     return wantarray ? %record : \%record;
 153 }
 154
 155
 156 sub biopiece2fixedstep
 157 {
 158     # Martin A. Hansen, November 2008.
 159
 160     # Converts a Biopiece record to a fixedStep entry.
 161
 162     my ( $record,   # Biopiece record
 163        ) = @_;
 164
 165     # Returns a list
 166
 167     my ( @entry, $beg, $vals );
 168
 169     if ( exists $record->{ "REC_TYPE" } and $record->{ "REC_TYPE" } eq 'fixed_step' )
 170     {
 171         if ( exists $record->{ 'CHR' } and exists $record->{ 'CHR_BEG' } and exists $record->{ 'STEP' } )
 172         {
 173             $beg = $record->{ 'CHR_BEG' } + 1;   # fixedStep is 1-based
 174             push @entry, "fixedStep chrom=$record->{ 'CHR' } start=$beg step=$record->{ 'STEP' }";
 175
 176             $vals = $record->{ 'VALS' };
 177
 178             map { push @entry, $_ } split ";", $vals;
 179         }
 180         else
 181         {
 182             Maasha::Common::error( qq(could not convert Biopiece record to fixedStep) );
 183         }
 184     }
 185
 186     return wantarray ? @entry : \@entry;
 187 }
 188
 189
 190 sub fixedstep_calc
 191 {
 192     # Martin A. Hansen, November 2008.
 193
 194     # Calculates fixedstep entries for a given BED file.
 195     # Implemented using large Perl arrays to avoid sorting.
 196     # Returns path to the fixedstep file.
 197
 198     my ( $bed_file,   # path to BED file
 199          $chr,        # chromosome name
 200          $use_score,  # flag indicating that the score field should be used - OPTIONAL
 201          $use_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, $use_score );
 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 ( $use_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 scores based on
 252     # clone count (or sequence read count) for the
 253     # begin position through the end position of each
 254     # BED entry.
 255
 256     my ( $fh_in,       # filehandle to BED file
 257          $use_score,   # flag indicating that the score field should be used - OPTIONAL
 258        ) = @_;
 259
 260     # Returns a bitarray.
 261
 262     my ( $bed, $clones, @begs, @lens, $i, @array );
 263
 264     while ( $bed = Maasha::UCSC::BED::bed_entry_get( $fh_in ) )
 265     {
 266         if ( $use_score ) {
 267             $clones = $bed->[ score ];
 268         } elsif ( defined $bed->[ name ] and $bed->[ name ] =~ /_(\d+)$/ ) {
 269             $clones = $1;
 270         } else {
 271             $clones = 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[ $_ ] += $clones } ( $bed->[ chromStart ] + $begs[ $i ] .. $bed->[ chromStart ] + $begs[ $i ] + $lens[ $i ] - 1 );
 281             }
 282         }
 283         else
 284         {
 285             map { $array[ $_ ] += $clones } ( $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 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 325
 326
 327 1;
 328
 329
 330 __END__