Bio::Graphics::Wiggle::Loader - Bio::Graphics::Wiggle::Loader documentation

  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');

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

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

use strict;

use Carp 'croak';
use Statistics::Descriptive;
use IO::Seekable;
use File::Spec;
use Bio::Graphics::Wiggle;
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;
      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 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 (<DATA>) {
      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__
#000000 black
#FFFFFF white
#0000FF blue
#00FF00 green
#FF0000 red
#FFFF00 yellow
#00FFFF cyan
#FF00FF magenta
#C0C0C0 gray


__END__