[biopieces.git] / code_perl / Maasha / UCSC / Wiggle.pm

package Maasha::UCSC::Wiggle;

# Copyright (C) 2007-2008 Martin A. Hansen.

# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 2
# of the License, or (at your option) any later version.

# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.

# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.

# http://www.gnu.org/copyleft/gpl.html


# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DESCRIPTION <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<


# Routines for manipulation of wiggle data.

# http://genome.ucsc.edu/goldenPath/help/wiggle.html


# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<


use strict;

use Data::Dumper;
use Maasha::Common;
use Maasha::Filesys;
use Maasha::Calc;

use vars qw( @ISA @EXPORT_OK );

require Exporter;

@ISA = qw( Exporter );

@EXPORT_OK = qw(
);


use constant {
    BITS        => 32,            # Number of bits in an integer
    SEQ_MAX     => 200_000_000,   # Maximum sequence size
    chrom       => 0,             # BED field names
    chromStart  => 1,
    chromEnd    => 2,
    name        => 3,
    score       => 4,
    strand      => 5,
    thickStart  => 6,
    thickEnd    => 7,
    itemRgb     => 8,
    blockCount  => 9,
    blockSizes  => 10,
    blockStarts => 11,
};


# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> FixedStep format <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<


sub fixedstep_entry_get
{
    # Martin A. Hansen, December 2007.

    # Given a file handle to a PhastCons file get the
    # next entry which is all the lines after a "fixedStep"
    # line and until the next "fixedStep" line or EOF.

    my ( $fh,   # filehandle
       ) = @_;

    # Returns a list of lines

    my ( $entry, @lines );

    local $/ = "\nfixedStep ";

    $entry = <$fh>;

    chomp $entry;

    @lines = split "\n", $entry;
    
    return if @lines == 0;

    $lines[ 0 ] =~ s/fixedStep?\s*//;
    $lines[ 0 ] = 'fixedStep ' . $lines[ 0 ];

    return wantarray ? @lines : \@lines;
}


sub fixedstep_entry_put
{
    # Martin A. Hansen, April 2008.

    # Outputs a block of fixedStep values.
    # Used for outputting wiggle data.

    my ( $entry,    # fixedStep entry
         $fh,       # filehandle - OPTIONAL
       ) = @_;

    # Returns nothing.

    $fh ||= \*STDOUT;

    print $fh join( "\n", @{ $entry } ), "\n";
}


sub fixedstep2biopiece
{
    # Martin A. Hansen, November 2008.

    # Converts a fixedstep entry to a Biopiece record.

    my ( $entry,   # fixedstep entry arrayref
       ) = @_;

    # Returns a hashref

    my ( $head, %record );

    $head = shift @{ $entry };

    if ( $head =~ /^fixedStep chrom=([^ ]+) start=(\d+) step=(\d+)$/ )
    {
        $record{ "REC_TYPE" } = "fixed_step";
        $record{ "CHR" }      = $1;
        $record{ "CHR_BEG" }  = $2 - 1;  # fixedStep is 1-based
        $record{ "STEP" }     = $3;
        $record{ "VALS" }     = join ";", @{ $entry };
    }
    else
    {
        Maasha::Common::error( qq(could not convert fixedStep entry to Biopiece record) );
    }

    return wantarray ? %record : \%record;
}


sub biopiece2fixedstep
{
    # Martin A. Hansen, November 2008.

    # Converts a Biopiece record to a fixedStep entry.
    
    my ( $record,   # Biopiece record
       ) = @_;

    # Returns a list

    my ( @entry, $beg, $vals );

    if ( exists $record->{ "REC_TYPE" } and $record->{ "REC_TYPE" } eq 'fixed_step' )
    {
        if ( exists $record->{ 'CHR' } and exists $record->{ 'CHR_BEG' } and exists $record->{ 'STEP' } )
        {
            $beg = $record->{ 'CHR_BEG' } + 1;   # fixedStep is 1-based
            push @entry, "fixedStep chrom=$record->{ 'CHR' } start=$beg step=$record->{ 'STEP' }";

            $vals = $record->{ 'VALS' };

            map { push @entry, $_ } split ";", $vals;
        }
        else
        {
            Maasha::Common::error( qq(could not convert Biopiece record to fixedStep) );
        }
    }

    return wantarray ? @entry : \@entry;
}


sub fixedstep_calc
{
    # Martin A. Hansen, November 2008.

    # Calculates fixedstep entries for a given BED file.
    # Implemented using large Perl arrays to avoid sorting.
    # Returns path to the fixedstep file.

    my ( $bed_file,   # path to BED file
         $chr,        # chromosome name
         $use_score,  # flag indicating that the score field should be used - OPTIONAL
         $use_log10,  # flag indicating that the score should be in log10   - OPTIONAL
       ) = @_;

    # Returns a string

    my ( $fixedstep_file, $fh_in, $fh_out, $array, $beg, $end, @block, $i );

    $fixedstep_file = "$bed_file.fixedstep";

    $fh_in = Maasha::Filesys::file_read_open( $bed_file );

    $array = fixedstep_calc_array( $fh_in, $use_score );

    close $fh_in;

    $fh_out = Maasha::Filesys::file_write_open( $fixedstep_file );

    $beg = 0;

    while ( ( $beg, $end ) = fixedstep_scan( $array, $beg ) and $beg != -1 )
    {
        @block = "fixedStep chrom=$chr start=" . ( $beg + 1 ) . " step=1";

        for ( $i = $beg; $i <= $end; $i++ )
        {
            if ( $use_log10 ) {
                push @block, sprintf "%.4f", Maasha::Calc::log10( $array->[ $i ] );
            } else {
                push @block, $array->[ $i ];
            }
        }

        fixedstep_entry_put( \@block, $fh_out );

        $beg = $end + 1;
    }

    close $fh_out;

    undef $array;

    return $fixedstep_file;
}


sub fixedstep_calc_array
{
    # Martin A. Hansen, November 2008.

    # Given a filehandle to a BED file for a single
    # chromosome fills an array with scores based on
    # clone count (or sequence read count) for the
    # begin position through the end position of each
    # BED entry.

    my ( $fh_in,       # filehandle to BED file
         $use_score,   # flag indicating that the score field should be used - OPTIONAL
       ) = @_;

    # Returns a bitarray.

    my ( $bed, $clones, @begs, @lens, $i, @array );

    while ( $bed = Maasha::UCSC::BED::bed_entry_get( $fh_in ) )
    {
        if ( $use_score ) {
            $clones = $bed->[ score ];
        } elsif ( $bed->[ name ] =~ /_(\d+)$/ ) {
            $clones = $1;
        } else {
            $clones = 1;
        }

        if ( $bed->[ blockCount ] and $bed->[ blockCount ] > 1 )
        {
            @begs = split ",", $bed->[ blockStarts ];
            @lens = split ",", $bed->[ blockSizes ];

            for ( $i = 0; $i < @begs; $i++ ) {
                map { $array[ $_ ] += $clones } ( $bed->[ chromStart ] + $begs[ $i ] .. $bed->[ chromStart ] + $begs[ $i ] + $lens[ $i ] - 1 );
            }
        }
        else
        {
            map { $array[ $_ ] += $clones } ( $bed->[ chromStart ] .. $bed->[ chromEnd ] - 1 );
        }
    }

    return wantarray ? @array : \@array;
}


sub fixedstep_scan
{
    # Martin A. Hansen, November 2008.

    # Given a fixedstep array locates the next block of
    # non-zero elements from a given begin position.
    # A [ begin end ] tuple of the block posittion is returned 
    # if a block was found otherwise and empty tuple is returned.

    my ( $array,   # array to scan
         $beg,     # position to start scanning from
       ) = @_;

    # Returns an arrayref.

    my ( $end );

    while ( $beg < scalar @{ $array } and not $array->[ $beg ] ) { $beg++ }

    $end = $beg;

    while ( $array->[ $end ] ) { $end++ }

    if ( $end > $beg ) {
        return wantarray ? ( $beg, $end ) : [ $beg, $end ];
    } else {
        return wantarray ? () : [];
    }
}


# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<


1;


__END__