File Coverage

lib/Bio/Graphics/Wiggle/Loader.pm

Criterion	Covered	Total	%
statement	197	377	52.2
branch	46	178	25.8
condition	25	106	23.5
subroutine	27	42	64.2
pod	5	27	18.5
total	300	730	41.1

package Bio::Graphics::Wiggle::Loader;

=head1 SYNOPSIS

  my $loader = Bio::Graphics::Wiggle::Loader->new('/base/directory/for/wigfiles','wibfilename');
  my $fh = IO::File->new('uploaded_file.txt');
  $loader->load($fh);

  my $gff3_file   = $loader->featurefile('gff3',$method,$source);
  my $featurefile = $loader->featurefile('featurefile');
  my @features    = $loader->features();

=head1 USAGE

This module loads Bio::Graphics::Wiggle files from source files that
use Jim Kent's "WIG" format:

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

Several data sets can be grouped together in a single WIG source
file. The load() method accepts the path to a WIG source file, and
will create one or more .wib ("wiggle binary") databases of
quantitative data in the directory indicated when you created the
loader. Call the featurefile() method to return a text file in either
GFF3 or Bio::Graphics::FeatureFile format, suitable for loading into a
gbrowse database.

=head2 METHODS

=over 4

=item $loader = Bio::Graphics::Wiggle::Loader->new('/base/directory' [,'my_data'])

Create a new loader. The first argument specifies the base directory
in which the loaded .wib files will be created. The second argument
specifies the base name for the created .wib files, or "track" if not
specified.

=item $loader->load($fh)

Load the data from a source WIG file opened on a filehandle.

=item $data = $loader->featurefile($type [,$method,$source])

Return the data corresponding to a GFF3 or
Bio::Graphics::FeatureFile. The returned file will have one feature
per WIG track, and a properly formatted "wigfile" attribute that
directs Bio::Graphics to the location of the quantitative data.

$type is one of "gff3" or "featurefile". In the case of "gff3", you
may specify an optional method and source for use in describing each
feature. In the case of "featurefile", the returned file will contain
GBrowse stanzas that describe a reasonable starting format to display
the data.

=item @features = $loader->features

Returns one or more Bio::Graphics::Features objects, which can be used to
create Bio::Graphics tracks with the wiggle_xyplot (and related) glyphs.

=item $loader->allow_sampling(1)

If allow_sampling() is passed a true value, then very large files
(more than 5 MB) will undergo a sampling procedure to find their
minimum and maximum values and standard deviation. Otherwise, file
will be read in its entirety to generate those statistics.

=back

=head2 EXTENSIONS

Several extensions to the WIG format "track" declaration are recognized.

=over 4

=item transform=

Specify a transform to be performed on all numeric data within this
track prior to loading into the binary wig file. Currently, the
following three declarations are recognized:

 transform=logtransform y' = 0          for  y == 0
                        y' = log(y)     for  y >  0
                        y' = -log(-y)   for  y <  0


 transform=logsquared  y' = log(y**2) for y != 0
                       y' = 0         for y == 0

 transform=none        y' = y   (no transform - the default)

=item trim=

Specify a trimming function to be performed on the data prior to
scaling. Currently, the following trim functions are recognized:

 trim=stdev1           trim to plus/minus 1 standard deviation of the mean
 trim=stdev2           trim to plus/minus 2 standard deviations of the mean (default)
 trim=stdevN           trim to plus/minus N standard deviations of the mean
 trim=none             no trimming

=back

Example entended track declaration:

 track type=wiggle_0 name="example" description="20 degrees, 2 hr"  \
       trim=stdev2 transform=logsquared

=cut

use strict;

use Carp 'croak';
use Statistics::Descriptive;
use IO::Seekable;
use File::Spec;
use Bio::Graphics::Wiggle;
use Bio::Graphics::FeatureFile;
use Text::ParseWords();
use File::stat;
use CGI 'escape';

use vars '%color_name';

# If a WIG file is very large (> 5 Mb)
use constant BIG_FILE         => 5_000_000;
use constant BIG_FILE_SAMPLES => 5_000;     # number of probes to make
use constant DEFAULT_METHOD   => 'microarray_oligo';
use constant DEFAULT_SOURCE   => '.';

sub new {
  my $class = shift;
  my $base      = shift
                 or croak "Usage: Bio::Graphics::Wiggle::Loader->new('/base/path','trackname')";
  my $trackname = shift || 'track';
  my $wigclass  = shift || 'Bio::Graphics::Wiggle';
  -d $base && -w _  or croak "$base is not a writeable directory";
  return bless {
		base            => $base,
		tracks          => {},
		trackname       => $trackname,
	        tracknum        => '000',
		track_options   => {},
		allow_sampling  => 0,
		wigclass        => $wigclass,
	       },ref $class || $class;
}
sub allow_sampling {
    my $self = shift;
    my $d    = $self->{allow_sampling};
    $self->{allow_sampling} = shift if @_;
    $d;
}
sub wigclass {
    my $self = shift;
    my $d    = $self->{wigclass};
    $self->{wigclass} = shift if @_;
    return $d;
}
sub basedir  { shift->{base}     }
sub wigfiles { shift->{wigfiles} }
sub conf_stanzas {
    my $self = shift;
    my ($method,$source) = @_;
    $method ||= DEFAULT_METHOD;
    $source ||= DEFAULT_SOURCE;

    my $tracks = $self->{tracks};
    my @lines = ();
    for my $track (sort keys %$tracks) {

	my $options    = $tracks->{$track}{display_options};
	my $name       = $options->{name} ||= $track;

	$options->{visibility} ||= 'dense';
	$options->{color}      ||= $options->{visibility} =~ /pack/i ? '255,0,0' : '0,0,0';
	$options->{altColor}   ||= $options->{visibility} =~ /pack/i ? '0,0,255' : '0,0,0';

	# stanza
	push @lines,"[$track]";
	if (my $graph_type = $options->{glyph}) {
	    if ($graph_type =~ /box/) {
		push @lines, "glyph       = wiggle_box";
	    }
	    else {
		push @lines,"glyph       = ".
		    ($graph_type =~/density/ ? 'wiggle_density' : 'wiggle_xyplot');
	    }
	}
	else {
	    push @lines,"glyph       = ".
		($options->{visibility}=~/pack/ ? 'wiggle_density' : 'wiggle_xyplot');
	}
	push @lines,"key         = $options->{name}"
	    if $options->{name};
	push @lines,"description = $options->{description}"
	    if $options->{description};
	if (my $color = $options->{color}) {
	    push @lines,"bgcolor=".format_color($color);
	}
	if (my $color = $options->{altColor}) {
	    push @lines,"fgcolor=" . format_color($color);
	}
	if (exists $options->{viewLimits} and my ($low,$hi) = split ':',$options->{viewLimits}) {
	    push @lines,"min_score   =  $low";
	    push @lines,"max_score   =  $hi";
	}
	if (exists $options->{maxHeightPixels} and my ($max,$default,$min) =
	    split ':',$options->{maxHeightPixels}) {
	    push @lines,"height  = $default";
	}
	push @lines,"smoothing        = $options->{windowingFunction}"
	    if $options->{windowingFunction};

	my $smoothing_window = $options->{smoothingWindow} || 0;

	push @lines,"smoothing window = $options->{smoothingWindow}"
	    if $options->{smoothingWindow};
	push @lines,'';
    }
    return join "\n",@lines;
}

sub featurefile {
  my $self             = shift;
  my $type             = shift;
  my ($method,$source) = @_;

  $method ||= DEFAULT_METHOD;
  $source ||= DEFAULT_SOURCE;

  $type ||= 'featurefile';
  $type =~ /^(gff3|featurefile)$/i
    or croak "featurefile type must be one of 'gff3' or 'featurefile'";

  my @lines;
  my $tracks = $self->{tracks};

  if ($type eq 'gff3') {
      push @lines,"##gff-version 3","";
  }
  else {
      push @lines,$self->conf_stanzas($method,$source),"";
  }

  for my $track (sort keys %$tracks) {
    my $options    = $tracks->{$track}{display_options};
    my $name       = $options->{name} ||= $track;
    my $seqids     = $tracks->{$track}{seqids};
    my $note       = escape($options->{description});
    my @attributes;
    push @attributes,qq(Name=$name)        if defined $name;
    push @attributes,qq(Note=$note)        if defined $note;

    # data, sorted by chromosome
    my @seqid  = sort keys %$seqids;

    for my $seqid (@seqid) {
      $seqid or next;
      $tracks->{$track}{seqids}{$seqid}{wig}->write();
      my $attributes = join ';',(@attributes,"wigfile=$seqids->{$seqid}{wigpath}");
      if ($type eq 'gff3') {
	  push @lines,join "\t",($seqid,$source,$method,
				 $seqids->{$seqid}{start},
				 $seqids->{$seqid}{end},
				 '.','.','.',
				 $attributes
	  );
      } else {
	  push @lines,'';
	  push @lines,"reference=$seqid";
	  push @lines,"$track $seqid.data $seqids->{$seqid}{start}..$seqids->{$seqid}{end} $attributes";
      }

    }

  }

  return join("\n",@lines)."\n";
}

sub features {
    my $self = shift;
    my $text = $self->featurefile('featurefile');
    my $file = Bio::Graphics::FeatureFile->new(-text=>$text);
    return $file->features;
}


sub load {
  my $self  = shift;
  my $infh  = shift;
  my $format = 'none';

  local $_;
 LINE: while (<$infh>) {
    chomp;
    next if /^#/;
    next unless /\S/;

    if (/^track/) {
      $self->process_track_line($_);
      next;
    }

    if (/^fixedStep/) {
      $self->process_fixed_step_declaration($_);
      $format = 'fixed';
    }

    if (/^variableStep/) {
      $self->process_variable_step_declaration($_);
      $format = 'variable';
    }

    if (/^\S+\s+\d+\s+\d+\s+-?[\dEe.]+/) {
      $self->process_first_bedline($_);
      $format    = 'bed';
    }

    if ($format ne 'none') {
      # remember where we are, find min and max values, return
      my $pos = tell($infh);
      $self->minmax($infh,$format eq 'bed' ? $_ : '')
	  unless $self->{track_options}{chrom} &&
	         exists $self->current_track->{seqids}{$self->{track_options}{chrom}}{min};
      seek($infh,$pos,0);

      $self->process_bed($infh,$_)        if $format eq 'bed';
      $self->process_fixedline($infh)     if $format eq 'fixed';
      $self->process_variableline($infh)  if $format eq 'variable';

      $format = 'none';
    }

    redo LINE if defined $_ && /^(track|variableStep|fixedStep)/;
  }

  return 1;
}

sub process_track_line {
  my $self      = shift;
  my $line      = shift;
  my @tokens    = shellwords($line);
  shift @tokens;
  my %options = map {split '='} @tokens;
  $options{type} eq 'wiggle_0' or croak "invalid/unknown wiggle track type $options{type}";
  delete $options{type};
  $self->{tracknum}++;
  $self->current_track->{display_options} = \%options;
}

sub process_fixed_step_declaration {
  my $self  = shift;
  my $line  = shift;
  my @tokens    = shellwords($line);
  shift @tokens;
  my %options = map {split '='} @tokens;
  exists $options{chrom}        or croak "invalid fixedStep line: need a chrom option";
  exists $options{start}        or croak "invalid fixedStep line: need a start option";
  exists $options{step}         or croak "invalid fixedStep line: need a step  option";
  $self->{track_options} = \%options;
}

sub process_variable_step_declaration {
  my $self  = shift;
  my $line  = shift;
  my @tokens    = shellwords($line);
  shift @tokens;
  my %options = map {split '='} @tokens;
  exists $options{chrom}        or croak "invalid variableStep line: need a chrom option";
  $self->{track_options} = \%options;
}

sub process_first_bedline {
  my $self  = shift;
  my $line  = shift;
  my @tokens    = shellwords($line);
  $self->{track_options} = {chrom => $tokens[0]};
}

sub current_track {
  my $self = shift;
  return $self->{tracks}{$self->{tracknum}} ||= {};
}

sub minmax {
  my $self   = shift;
  my ($infh,$bedline) = @_;
  local $_;

  my $transform  = $self->get_transform;

  my $seqids = ($self->current_track->{seqids} ||= {});
  my $chrom  = $self->{track_options}{chrom};

  if ($self->allow_sampling && (my $size = stat($infh)->size) > BIG_FILE) {
      warn "Wiggle file is very large; resorting to genome-wide sample statistics for $chrom.\n";
      $self->{FILEWIDE_STATS} ||= $self->sample_file($infh,BIG_FILE_SAMPLES);
      for (keys %{$self->{FILEWIDE_STATS}}) {
	$seqids->{$chrom}{$_} = $self->{FILEWIDE_STATS}{$_};
      }
      return;
  }

  my %stats;
  if ($bedline) {  # left-over BED line
      my @tokens = split /\s+/,$bedline;
      my $seqid  = $tokens[0];
      my $value  = $tokens[-1];
      $value = $transform->($self,$value) if $transform;
      $stats{$seqid} ||= Statistics::Descriptive::Sparse->new();
      $stats{$seqid}->add_data($value);
  }

  while (<$infh>) {
      last if /^track/;
      last if /chrom=(\S+)/ && $1 ne $chrom;
      next if /^\#|fixedStep|variableStep/;
      my @tokens = split(/\s+/,$_) or next;
      my $seqid  = @tokens > 3 ? $tokens[0] : $chrom;
      my $value  = $tokens[-1];
      $value = $transform->($self,$value) if $transform;
      $stats{$seqid} ||= Statistics::Descriptive::Sparse->new();
      $stats{$seqid}->add_data($value);
  }

  for my $seqid (keys %stats) {
      $seqids->{$seqid}{min}    = $stats{$seqid}->min();
      $seqids->{$seqid}{max}    = $stats{$seqid}->max();
      $seqids->{$seqid}{mean}   = $stats{$seqid}->mean();
      $seqids->{$seqid}{stdev}  = $stats{$seqid}->standard_deviation();
  }
}

sub sample_file {
    my $self = shift;

    my ($fh,$samples) = @_;

    my $transform  = $self->get_transform;

    my $stats = Statistics::Descriptive::Sparse->new();

    my $size = stat($fh)->size;
    my $count=0;
    while ($count < $samples) {
	seek($fh,int(rand $size),0) or die;
	scalar <$fh>; # toss first line
	my $line = <$fh>; # next full line
	$line or next;
	my @tokens = split /\s+/,$line;
	my $value  = $tokens[-1];
	next unless $value =~ /^[\d\seE.+-]+$/; # non-numeric
	$value = $transform->($self,$value) if $transform;
	$stats->add_data($value);
	$count++;
    }

    return {
	min   => $stats->min,
	max   => $stats->max,
	mean  => $stats->mean,
	stdev => $stats->standard_deviation,
    };
}

sub get_transform {
    my $self = shift;
    my $transform  = $self->current_track->{display_options}{transform};
    return $self->can($transform) if $transform;
}

# one and only transform currently defined
# Natural log of the square of the value.
# Return 0 if the value is 0
sub logsquared {
    my $self  = shift;
    my $value = shift;
    return 0 if $value == 0;
    return log($value**2);
}

sub logtransform {
  my $self  = shift;
  my $value = shift;
  return 0 if $value == 0;
  if ($value < 0) {
    return -log(-$value);
  } else {
    return log($value);
  }
}

sub process_bed {
  my $self = shift;
  my $infh = shift;
  my $oops = shift;
  my $transform  = $self->get_transform;
  $self->process_bedline($oops) if $oops;
  while (<$infh>) {
    last if /^track/;
    next if /^#/;
    chomp;
    $self->process_bedline($_);
  }
}

sub process_bedline {
  my $self = shift;
  my ($line,$transform) = @_;

  my ($seqid,$start,$end,$value) = split /\s+/,$line;
  $value = $transform->($self,$value) if $transform;
  $start++;   # to 1-based coordinates

  my $wigfile = $self->wigfile($seqid);
  $wigfile->set_range($start=>$end, $value);

  # update span
  $self->current_track->{seqids}{$seqid}{start} = $start
    unless exists $self->current_track->{seqids}{$seqid}{start}
      and $self->current_track->{seqids}{$seqid}{start} < $start;

  $self->current_track->{seqids}{$seqid}{end}   = $end
    unless exists $self->current_track->{seqids}{$seqid}{end}
      and $self->current_track->{seqids}{$seqid}{end} > $end;
}

sub process_fixedline {
  my $self  = shift;
  my $infh  = shift;
  my $seqid   = $self->{track_options}{chrom};
  my $wigfile = $self->wigfile($seqid);
  my $start   = $self->{track_options}{start};
  my $step    = $self->{track_options}{step};
  my $span    = $wigfile->span;

  # update start and end positions
  $self->{track_options}{span} ||= $wigfile->span || 1;
  my $chrom = $self->current_track->{seqids}{$seqid};
  $chrom->{start} = $start
      if !defined $chrom->{start} || $chrom->{start} > $start;
  my $end = $chrom->{start} + $span - 1;
  $chrom->{end}   = $end
      if !defined $chrom->{end} || $chrom->{end} < $end;

  my $transform  = $self->get_transform;

  # write out data in 500K chunks for efficiency
  my @buffer;
  while (<$infh>) {
    last if /^(track|variableStep|fixedStep)/;
    next if /^#/;
    chomp;
    push @buffer,$_;
    if (@buffer >= 500_000) {
	@buffer = map {$transform->($self,$_)} @buffer if $transform;
	$wigfile->set_values($start=>\@buffer);
	my $big_step = $step * @buffer;
	$start += $big_step;
	$self->current_track->{seqids}{$seqid}{end} = $start + $big_step - 1 + $span;
	@buffer = (); # reset at the end
    }

  }
  @buffer = map {$transform->($self,$_)} @buffer if $transform;
  $wigfile->set_values($start=>\@buffer)         if @buffer;
  $self->current_track->{seqids}{$seqid}{end} =
      $start + @buffer*$step - 1 + $span;
}

sub process_variableline {
  my $self  = shift;
  my $infh  = shift;
  my $seqid   = $self->{track_options}{chrom};
  my $chrom   = $self->current_track->{seqids}{$seqid};
  my $wigfile = $self->wigfile($seqid);
  my $span    = $wigfile->span;
  my $transform  = $self->get_transform;

  while (<$infh>) {
    last if /^(track|variableStep|fixedStep)/;
    next if /^#/;
    chomp;
    my ($start,$value) = split /\s+/ or next;
    $value = $transform->($self,$value) if $transform;
    eval {
	$wigfile->set_value($start=>$value);
	1;
    } or croak "Data error on line $.: $_\nDetails: $@";

    # update span
    $chrom->{start} = $start
	if !defined $chrom->{start} || $chrom->{start} > $start;
    my $end = $start + $span - 1;
    $chrom->{end}   = $end
	if !defined $chrom->{end} || $chrom->{end} < $end;

  }

  $self->current_track->{seqids}{$seqid}{end}
        ||= $self->current_track->{seqids}{$seqid}{start};
}

sub wigfile {
  my $self  = shift;
  my $seqid = shift;
  my $ts    = time();
  my $current_track = $self->{tracknum};
  my $tname         = $self->{trackname} || 'track';
  unless (exists $self->current_track->{seqids}{$seqid}{wig}) {
    my $path    = File::Spec->catfile($self->{base},"$tname\_$current_track.$seqid.$ts.wib");
    my @stats;
    foreach (qw(min max mean stdev)) {
	my $value = $self->current_track->{seqids}{$seqid}{$_} ||
	    $self->{FILEWIDE_STATS}{$_} || next;
	push @stats,($_=>$value);
   }

    my $step = $self->{track_options}{step} || 1;
    my $span = $self->{track_options}{span} ||
	$self->{track_options}{step}        ||
	1;
    my $trim      = $self->current_track->{display_options}{trim} || 'stdev10';
    my $transform = $self->current_track->{display_options}{transform};
    my $class   = $self->wigclass;
    unless ($class->can('new')) {
	warn "loading $class";
	eval "require $class";
	die $@ if $@;
    }
    my $wigfile = $class->new(
	$path,
	1,
	{
	    seqid => $seqid,
	    step  => $step,
	    span  => $span,
	    trim  => $trim,
	    @stats,
	},
	);
    $wigfile or croak "Couldn't create wigfile $wigfile: $!";
    $self->current_track->{seqids}{$seqid}{wig}     = $wigfile;
    $self->current_track->{seqids}{$seqid}{wigpath} = $path;
  }
  return $self->current_track->{seqids}{$seqid}{wig};
}

sub format_color {
  my $rgb = shift;
  return $rgb unless $rgb =~ /\d+,\d+,\d+/;
  my ($r,$g,$b) = split ',',$rgb;
  my $hex = '#'.join '',map {sprintf("%02X",$_)}($r,$g,$b);
  return translate_color($hex);
}

# use English names for the most common colors
sub translate_color {
  my $clr = shift;
  unless  (%color_name) {
    while () {
      chomp;
      my ($hex,$name) = split or next;
      $color_name{$hex} = $name;
    }
  }
  return $color_name{$clr} || $clr;
}

# work around an annoying uninit variable warning from Text::Parsewords
sub shellwords {
    my @args = @_;
    return unless @args;
    foreach(@args) {
	s/^\s+//;
	s/\s+$//;
	$_ = '' unless defined $_;
    }
    my @result = Text::ParseWords::shellwords(@args);
    return @result;
}

1;


__DATA__