murasaki-mpi 1.68.6-4build1 / usr / lib / murasaki / cgr-image.pl-mpi

#!/usr/bin/perl

#Copyright (C) 2006-2008 Keio University
#(Kris Popendorf) <comp@bio.keio.ac.jp> (2006)
#
#This file is part of Murasaki.
#
#Murasaki 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 3 of the License, or
#(at your option) any later version.
#
#Murasaki 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 Murasaki.  If not, see <http://www.gnu.org/licenses/>.

###############
## view histogram details as a CGR graph
###############

use strict;
use Image::Magick;
use Getopt::Long;
use Getopt::Std;
use Pod::Usage;
$Getopt::Std::STANDARD_HELP_VERSION=1;

our ($highest,$lowest);
our $res;
our $lasttick;
our $samples;
our $detailfile;
our $seedcount;
our $noIntermediates;

my ($help,$man);
my $err=GetOptions('help|?' => \$help, man => \$man,"samples=i"=>\$samples,"noIntermediates|1"=>\$noIntermediates);
pod2usage(1) if $help or $#ARGV<0 or !$err;
pod2usage(-exitstatus => 0, -verbose => 2) if $man;

foreach $detailfile (@ARGV){
  print "Loading seedlist from $detailfile...\n";
  open(INF,"<$detailfile");
  my @pixels;
  my $sample;
  my @keyl;
  while(<INF>){
    my ($keyn,$keya,$count)=split(/\s/,$_);
    $keya=~s/\.//g;
    push(@keyl,[$keya,$count]);
    $seedcount+=$count;
  }
  print "$seedcount seeds total. Mode switch predicted at sample ".log(sqrt($seedcount))/log(2)."\n";
  close(INF);
  die "No keys!" unless $#keyl>=0;
  #  print join("\n",@keyl),"\n";
  my $sample=$keyl[0][0];
  print "Available samples: ".(length($sample))."\n";
  $samples=(length($sample)) unless $samples and $samples<=(length($sample));
  $res=2<<($samples-1);
  my @graph=growGraph(1,\@keyl,0);
  print "Writing $res x $res image to file $detailfile.cgr.png...\n";
  sneakyImage("$detailfile.cgr.png",\@graph);
}

sub grabCount {
  my ($key,$keylr)=@_;
  my $count;
  foreach(grep {$_[0]=~m/$key\$/} @$keylr){
    $count+=$_[1];
  }
  return $count;
}

sub pngGraph {
  print "(using imagemagick pixels may be slow)\n";
  my ($filename,$graph)=@_;
  my $image=Image::Magick->new();
  $image->Set(size=>join("x",$res,$res));
  $image->Read('xc:white');
  my $err;
  #  $err=$image->set(type=>'Grayscale');
  warn $err if $err;
  my ($i,$max)=(0,$res*$res);
  $lasttick=time;
  foreach my $y (0..($res-1)){
    foreach my $x (0..($res-1)){
      my $v=255-int(255*(($$graph[$y][$x])/$highest));
      $v=sprintf('%02x',$v);
      $v="#".$v.$v.$v;
      $err=$image->Set("pixel[$x,$y]"=>$v);
      warn $err if $err;
      ticker($i,$max);
      $i++;
    }
  }
  $err=$image->Write(filename=>$filename);
  warn $err if $err;
}

sub ticker {
  if(time>$lasttick+5){
    $lasttick=time;
    print int($_[0]/$_[1]*100)."%...\n";
  }
}

sub growGraph {
  my ($level,$keylr,$old_graphr)=@_;
  print "Subsample level $level\n";
  my @graph=doubleGraph($old_graphr);
  $highest=0;
  foreach my $keyd (@$keylr){
    my ($key,$count)=(substr($$keyd[0],-$level,$level),$$keyd[1]);
    my @coords=key2coords($key);
    $graph[$coords[0]][$coords[1]]+=$count;
    $highest=$graph[$coords[0]][$coords[1]] if $graph[$coords[0]][$coords[1]]>$highest;
  }

  sneakyImage("$detailfile.cgr.$level.png",\@graph) unless $noIntermediates;

  return (length($$keylr[0][0])>$level and $level<$samples) ?
    growGraph($level+1,$keylr,\@graph):@graph;
}

sub doubleGraph {
  my $src=pop;
  return ([0,0],[0,0]) unless ref($src) eq "ARRAY";
  my @ret;
  my $buckets=@$src;
  $buckets*=$buckets;
  my $redist=$buckets<$seedcount;
  print $redist ? "Doubling graph and distributing...\n":"Doubling graph...\n";
  foreach my $row (@$src){
    if($row){
      my @drow;
      foreach(@$row){
	$_>>=2 if $redist;
	push(@drow,$_,$_);
      }
      push(@ret,\@drow,[@drow]);
    }else{
      push(@ret,[],[]);
    }
  }
  @{$src}=(); #we can clear the old map after that.
  return @ret;
}

sub selectnext {
  my @ret;
  foreach(qw{G C A T}){
    push(@ret,[map {substr($_,0,length($_)-1)} grep(/$_\$/,@_)]);
  }
  return @ret;
}

sub findLow {
  my $graph=shift;
  $lowest=$highest;
  foreach my $row (@$graph){
    if($row){
      foreach my $val (@$row){
	$lowest=$val if $val<$lowest;
      }
    }
  }
  return $lowest;
}

sub firstOver {
  my $limit=shift;
  for(0..$#_){
    return $_ if $limit<$_[$_];
  }
  return -1;
}

sub firstUnder {
  my $limit=shift;
  for(0..$#_){
    return $#_-$_ if $limit>$_[$#_-$_];
  }
  return -1;
}

sub findNiceRange {
  my $graph=shift;
  my @uberlist;
  foreach my $row (@$graph){
    push(@uberlist,@$row) if $row;
  }
  @uberlist=sort {$a<=>$b} @uberlist;
  $highest=$uberlist[$#uberlist];
  $lowest=$uberlist[0];
  return ($uberlist[int($#uberlist*.01)],$uberlist[int($#uberlist*.99)+1]);
}

sub writePGM { #holy crap PGMs are easy
  my ($filename,$graph)=@_;
  open(PGM,">$filename") or die "Couldnt write to $filename";
  my $localres=@$graph;
  my $maxout=65535;
  print PGM <<ENDHEADER;
P2
#$filename
$localres $localres
$maxout
ENDHEADER

  my ($bottom,$top)=findNiceRange($graph);
  my $range=$top-$bottom;
  print "Low: $lowest High: $highest\n";
  print "White: $bottom. Black: $top Range: $range\n";
  $range=1 unless $range; #no divide by 0 screw ups
  my ($i,$max)=(0,$localres*$localres);
  $lasttick=time;
  foreach my $y (0..($localres-1)){
    foreach my $x (0..($localres-1)){
      my $v=$maxout-int($maxout*(($$graph[$y][$x]-$bottom)/$range));
      $v=0 if $v<0;
      $v=$maxout if $v>$maxout;
      print PGM $v." ";
    }
    print PGM "\n";
    ticker($y,$localres-1);
  }
  close PGM;
}

sub sneakyImage {
  my ($filename,$graph)=@_;
  #make PGM version
  print "Writing uncompressed PGM version...\n";
  writePGM("$filename.pgm",$graph);
  print "Converting PGM to $filename with ImageMagick.\n";
  #make PNG version
  my $image=Image::Magick->new();
  my $err=$image->Set(depth=>16); #this seems to be ignored anyway
  warn $err if $err;
#  $image->Set(size=>join("x",$res,$res));
  $image->Read("$filename.pgm");
  $image->Write($filename);
  print "Erasing temporary PGM file.\n";
  system("rm $filename.pgm");
}

sub showmap {
  return map {join(" ",@$_)."\n"} @_;
}

sub key2coords {
  my $key=shift;
  $key=reverse $key;
  my @coord=(0,0);
  my %field=(C => [0,0],
	     G => [0,1],
	     A => [1,0],
	     T => [1,1]);
  foreach my $bit (split(//,$key)){
    next if $bit eq '.';
    $coord[0]<<=1;
    $coord[1]<<=1;
    my @add=@{$field{$bit}};
    $coord[0]+=$add[0];
    $coord[1]+=$add[1];
  }
  return @coord;
}

sub HELP_MESSAGE(){
  pod2usage(1)
}

__END__

=head1 NAME

cgr-image.pl -- Chaos Graph Representation histogram viewer

=head1 SYNOPSIS

cgr-iamge.pl [options] <alignment1.histogram.details> [alignment2.histogram.details> ...]

=head1 OPTIONS

   --samples=n -> limits to n subsamples (remember resolution is 2^n)
   --nointermediates|1 => skip all the intermediate outputs
=cut