#!/usr/bin/perl
# wig2wig.pl

=head1 about
 convert simple .wig/.bed input table of genome location,signal
 to Bio::Graphics::Wiggle for Gbrowse
 for modenc 0903 RNA-Seq .wig data

=head sample input

track chrom=2L name=matecov39
variableStep chrom=2L
#loc   materat mates  reads
10000   1.38    111     27
10005   1.33    37      9
10010   1.42    74      21
10015   1.44    37      10

.. or default location,reads 2 columns

track chrom=2L name=matecov39
variableStep chrom=2L
#loc   reads
10000   27
10005   9
10010   21
10015   10

=head1 sample usage

 # input table.wig with several columns of read counts, pick col=3 from 0, log(count), min/max 0..10
 cat ../matecov39.merg3.wig | env id=readcov col=3 log=1 step=5 smin=0 smax=10 tn=3 \
   perl -I$dmx $dmx/wig2wig.pl > materat.tracks.gff

 # this data column is already logged materat=log(mate/read)
 cat ../matecov39.merg3.wig | env id=materat col=1 log=0 step=5 smin=-5 smax=5 tn=4 \
   perl -I$dmx $dmx/wig2wig.pl > materat.tracks.gff

 if data is stranded, it is parsed from track comment(see below), and set -sig for neg strand

=cut 

use strict;
use Bio::Graphics::Wiggle;

my $seqid= $ENV{'id'}||"DmRnaSeq0903";
my $tn=    $ENV{'tn'}||1; 
my $do_max= $ENV{'max'} || 0; 
my $do_log= $ENV{'log'} || 0; 
my $scorecol=  $ENV{'col'}||1; 

my $step= $ENV{'step'} || 1; 
my $span= $ENV{'span'} || $step; 
my $minsig= $ENV{'smin'}; 
my $maxsig= $ENV{'smax'}; 

my $sepstrand= 1; # 0 == same wig
##my ($minsig,$maxsig)= (0,250);
my ($writeable, $as_abs)= (1) x 10;
my($wig,$wigfile,$lchr,$orient,$name, $nval, $ov, @oval);
my($chrstart,$chrstop) = (undef) x 2;

$orient=0; $name= $seqid;
$sepstrand= 0 if ($do_max);
unless(defined($minsig) and defined($maxsig)) {
if($sepstrand) {
  ($minsig,$maxsig)= ($do_log) ? (0,10) : (0,250);
} else {
  ($minsig,$maxsig)= ($do_log) ? (-9,9) : (-100,100); #
}
}

sub puttrack {
  my($mn,$mx,$md, $nval);
  $nval= scalar(@oval);
  @oval= sort{$b <=> $a} @oval; 
  $mx=$oval[0]; $mn=$oval[-1]; $md=$oval[int($nval/2)]; 
  $mx=int(1000*$mx)/1000; $mn=int(1000*$mn)/1000; $md=int(1000*$md)/1000;
  print join("\t",$lchr,$seqid,"region",$chrstart,$chrstop,".",$orient,".",
    "Name=$name;wigfile=$wigfile;opt=log:$do_log,maxv:$do_max;stat=n:$nval,med:$md,min:$mn,max:$mx"),"\n"
    if($lchr);
  @oval=();
}

while(<>) {
  if(/^track/) {
    if( m/name=(\S+)/) { $name=$1; $name=~s/^["']//; $name=~s/["']$//; }
    $orient= ($name =~ /\(\-\)/) ? "-" : ($name =~ /\(\+\)/) ? "+" : "0";

  } elsif (/^variableStep/) {  # variableStep chrom=chr2L
    my($chr)= m/chrom=(\S+)/; $chr =~ s/^chr//;    

    if($chr ne $lchr or not $wig) {
      puttrack();
      undef $wig if $wig; #? release mem

      $wigfile= "track$tn.$seqid.$chr";
      $wigfile .= ".$orient" if $sepstrand;
      $wigfile .= ".wig";
      
      # this doesnt erase existing wigfile; use wig->erase to drop data
      $wig = Bio::Graphics::Wiggle->new( $wigfile, $writeable, 
      { seqid => $seqid, 
        step => $step, span => $span, 
        min => $minsig,  max => $maxsig, });
      ($chrstart,$chrstop) = (undef) x 2;
      $lchr= $chr;  $nval= 0; @oval=();
    }
    
  } elsif(/^(\d+)\s+([\d\.+-]+)/) {
    my ($start,$value) = ($1,$2);
    if($scorecol>1) { $value=(split)[$scorecol]; }

    $value= log(1 + $value) if($do_log && $value>0); # +1 makes small/ 0.1 > 0
    $value= -$value if(!$sepstrand && $orient eq "-");  ## dont need change value with diff .wig files 
    if( $do_max ) { $ov= $wig->set_valuemax( $start, $value, $as_abs);  }
    else { $ov= $wig->set_value($start,$value); }
    $nval++; push(@oval, $ov); #? for stats
    $chrstop = $start if $chrstop < $start;
    $chrstart= $start unless(defined $chrstart and $chrstart < $start);
  }
}

puttrack();