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