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
  39 use vars qw( @ISA @EXPORT_OK );
  40
  41 require Exporter;
  42
  43 @ISA = qw( Exporter );
  44
  45 @EXPORT_OK = qw(
  46 );
  47
  48
  49 use constant {
  50     chrom       => 0,   # BED field names
  51     chromStart  => 1,
  52     chromEnd    => 2,
  53     name        => 3,
  54     score       => 4,
  55     strand      => 5,
  56     thickStart  => 6,
  57     thickEnd    => 7,
  58     itemRgb     => 8,
  59     blockCount  => 9,
  60     blockSizes  => 10,
  61     blockStarts => 11,
  62 };
  63
  64
  65 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> FixedStep format <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  66
  67
  68 sub fixedstep_entry_get
  69 {
  70     # Martin A. Hansen, December 2007.
  71
  72     # Given a file handle to a PhastCons file get the
  73     # next entry which is all the lines after a "fixedStep"
  74     # line and until the next "fixedStep" line or EOF.
  75
  76     my ( $fh,   # filehandle
  77        ) = @_;
  78
  79     # Returns a list of lines
  80
  81     my ( $entry, @lines );
  82
  83     local $/ = "\nfixedStep ";
  84
  85     $entry = <$fh>;
  86
  87     chomp $entry;
  88
  89     @lines = split "\n", $entry;
  90
  91     return if @lines == 0;
  92
  93     $lines[ 0 ] =~ s/fixedStep?\s*//;
  94     $lines[ 0 ] = 'fixedStep ' . $lines[ 0 ];
  95
  96     return wantarray ? @lines : \@lines;
  97 }
  98
  99
 100 sub fixedstep_entry_put
 101 {
 102     # Martin A. Hansen, April 2008.
 103
 104     # Outputs a block of fixedStep values.
 105     # Used for outputting wiggle data.
 106
 107     my ( $chr,      # chromosome
 108          $beg,      # start position
 109          $block,    # list of scores
 110          $fh,       # filehandle - OPTIONAL
 111        ) = @_;
 112
 113     # Returns nothing.
 114
 115     $beg += 1;   # fixedStep format is 1 based.
 116
 117     $fh ||= \*STDOUT;
 118
 119     print $fh "fixedStep chrom=$chr start=$beg step=1\n";
 120
 121     map { print $fh "$_\n" } @{ $block };
 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;
 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, $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             push @entry, "fixedStep chrom=$record->{ 'CHR' } start=$record->{ 'CHR_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        ) = @_;
 202
 203     # Returns a string
 204
 205     my ( $fixedstep_file, $fh_in, $fh_out, $array, $beg, $end, @block );
 206
 207     $fixedstep_file = "$bed_file.fixedstep";
 208
 209     $fh_in = Maasha::Filesys::file_read_open( $bed_file );
 210
 211     $array = fixedstep_calc_array( $fh_in, $use_score );
 212
 213     close $fh_in;
 214
 215     $fh_out = Maasha::Filesys::file_write_open( $fixedstep_file );
 216
 217     $beg = 0;
 218
 219     while ( ( $beg, $end ) = fixedstep_scan( $array, $beg ) and $beg )
 220     {
 221         @block = @{ $array }[ $beg .. $end ];
 222
 223         fixedstep_entry_put( $chr, $beg, \@block, $fh_out );
 224
 225         $beg = $end;
 226     }
 227
 228     close $fh_out;
 229
 230     return $fixedstep_file;
 231 }
 232
 233
 234 sub fixedstep_calc_array
 235 {
 236     # Martin A. Hansen, November 2008.
 237
 238     # Given a filehandle to a BED file for a single
 239     # chromosome fills an array with scores based on
 240     # clone count (or sequence read count) for the
 241     # begin position through the end position of each
 242     # BED entry.
 243
 244     my ( $fh_in,       # filehandle to BED file
 245          $use_score,   # flag indicating that the score field should be used - OPTIONAL
 246        ) = @_;
 247
 248     # Returns arrayref.
 249
 250     my ( $bed, $clones, @array );
 251
 252     while ( $bed = Maasha::UCSC::BED::bed_entry_get( $fh_in, 5 ) )
 253     {
 254         if ( $use_score ) {
 255             $clones = $bed->[ score ];
 256         } elsif ( $bed->[ name ] =~ /_(\d+)$/ ) {
 257             $clones = $1;
 258         } else {
 259             $clones = 1;
 260         }
 261
 262         map { $array[ $_ ] += $clones } $bed->[ chromStart ] .. $bed->[ chromEnd ] - 1;
 263     }
 264
 265     return wantarray ? @array : \@array;
 266 }
 267
 268
 269 sub fixedstep_scan
 270 {
 271     # Martin A. Hansen, November 2008.
 272
 273     # Given a fixedstep array locates the next block of
 274     # non-zero elements from a given begin position.
 275     # A [ begin end ] tuple of the block posittion is returned
 276     # if a block was found otherwise and empty tuple is returned.
 277
 278     my ( $array,   # array to scan
 279          $beg,     # position to start scanning from
 280        ) = @_;
 281
 282     # Returns an arrayref.
 283
 284     my ( $end );
 285
 286     while ( $beg < scalar @{ $array } and not $array->[ $beg ] ) { $beg++ }
 287
 288     $end = $beg;
 289
 290     while ( $array->[ $end ] ) { $end++ }
 291
 292     if ( $end > $beg ) {
 293         return wantarray ? ( $beg, $end ) : [ $beg, $end ];
 294     } else {
 295         return wantarray ? () : [];
 296     }
 297 }
 298
 299
 300 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 301
 302
 303 1;
 304
 305
 306 __END__