#! /usr/local/bin/perl
#
################################################################################
# Requirements                                                                 #
################################################################################
#
require 5.002;
#require "GUSI.ph";   #MacOS
#
# External References
#
use GD;
use Carp;
use FileHandle;
use File::Path;
use File::Basename;
use Getopt::Std;    #UNIX

use vars qw($opt_d $opt_n $opt_h);
#
################################################################################
# density window default values                                                #
################################################################################
#
my $frame    = 100; #width of the window shifted alond the sequences
my $stepsize = 1  ; #steps between the windows
my $choose   = 0  ; #only 1 or 0! 0=no connecting lines between data points
                    #             1=grey lines between data points
#
################################################################################
# Main program global variables                                                #
################################################################################
#
$| = 1;                               #print immediately

my $abschnitt    ;
my $colvar       ;
my $file         ;
my @files        ;
my $i         = 0;
my $laenge       ;
my $maxclones = 0;
my $p         = 0;
my $seq          ;
my $segment      ;
my $signal_dens  ;
my $white        ;
my $green        ;
my $black        ;       
my $blue         ;
my $c3           ;
my $c4           ;
my $c5           ;
my $c6           ;
my $c7           ;
my $c8           ;
my $c9           ;
my $c10          ;
my $color     = 0;
my $grey         ;
my $pre_x     = 0;
my $pre_y     = 0;
#
################################################################################
# Grafic output default values                                                 #
################################################################################
#
my $offset = 50;
my $spreiz = 3;
my $connectlineclr;
#
#
# MAN page
#
my $MAN = <<MAN;
NAME

     Plot_5mCpG_density_to_png.pl: 
     Program converts Seq1 files into 5mCpG density plot

AVAILABILITY

     Requires perl 5.002

DESCRIPTION

     The methylation density, i.e. the number of methylated cytosines in 
     a given sequence length, is considered to be important for the effect 
     methylation has on gene expression. Accordingly, variations of the 
     methylation density within the investigated population of PCR products 
     (i.e. cells) are interesting. "Seq2DensPlot" displays a methylation 
     density plot. The program reads files in seq1-format carrying the 
     extension ".seq1". 
     The methylation density of individual clones is calculated as ratio 
     of methylated CpG to total CpG in a 100 bp window shifted in 1 bp steps 
     over the sequences. These are default values and they can easily be changed 
     in the source code (see: #density window default values#). 
     The number of clones with identical densities at a given position 
     is indicated by the color of the corresponding data point. 
     The color spectrum reaches from black for one clone to red for more 
     than ten clones.

SYNOPSIS

     Plot_5mCpG_density_to_png.pl [-d] directory [-n] outfile [-h]
     (replaced by file selector box under MacOS)

OPTIONS
    -d string directory with seq1-files
    -n string outfile name
    -h        Prints this message
 
AUTHOR(S)

     Christoph Grunau (cgrunau\@imb-jena.de)

CHANGE LOG

     Nov.   1st, 1998 CG   Initial version
     Okt.   1st, 1999 CG   Minor changes, MacPerl version
     Apr.  10th, 2000 CG   PNG version for UNIX
MAN
#
# Command line parsing 
#
getopts('d:n:h') || croak($MAN); #UNIX
#$opt_d = &GetFolder("Choose input folder"); #MacOS
#$opt_n = &PutFile("Choose output png file"); #MacOS

#
# Check command line arguments
#
croak($MAN) if ( $opt_h || ! defined($opt_d) || $opt_d eq "" );

chdir $opt_d;
@files = <*.seq1>;          # write all.seq1-files into a list
$maxclones = $#files + 1;        # determine number of files
#
# finish program when no files *.seq1 found
#
if ($maxclones <1) {print "Sorry. No files available! Please try option -h.\n"}
else
{
print "Processing files";
#-------------------------------------------------------
# read first sequence to determine length:
     
     $filelist = $files[0];      #initialize variables
     $seq = "";

     open (SEQUENCE, "<$filelist") || die ("Cannot open file '$filelist'.");
           while (<SEQUENCE>)
	                {
	                $file = $_;
	                chomp $file;
                 unless ($file =~ /^>/){$seq = $seq . $file}
	                }
     close SEQUENCE || die ("Cannot close open file '$filelist'.");

     $laenge = length($seq);
     
#-------------------------------------------------------
# initialize matrix:
     
     #for ($z = 0; $z<=$laenge; $z++)
     for ($z = $frame/2+1; $z<$laenge-$frame/2; $z=$z+$stepsize) 
	     {
	     for ($i = 0; $i<=$maxclones; $i++) 
	         {
	         $C_counter[$z][$i] = 0;
	         $c_counter[$z][$i] = 0;
		 $signal[$z][$i]    = 111;
	         }
	      }
$i = 0;
$z = 0;
#------------------------------------------------------
# read sequences:

foreach $filelist(@files)
     	{
     	$seq = "";
      #sequences from list
	     open (SEQUENCE, "<$filelist")|| die ("Cannot open file '$filelist'.");
	     print (".");                    #show program is running
	     while (<SEQUENCE>)
		          {
		          $file = $_;
		          chomp $file;            #read line by line, cut Newline
		          #
		          # assemble sequences from lines
		          #
		          unless ($file =~ /^>/){$seq = $seq . $file}
		          }
     	close SEQUENCE || die ("Cannot close open file '$filelist'.");
#----------------------------------------------------
# write values into matrix:

$i++;	
for ($z = $frame/2+1; $z<$laenge-$frame/2; $z=$z+$stepsize)
    {
    $segment = substr($seq,$z-$frame/2,$frame);
    unless ($segment =~ m/-+$/i || $segment =~ m/^-+/i)
           {
           while ($segment =~ /Cg/g)   {$C_counter[$z][$i]++;}
           while ($segment =~ /cg/g)   {$c_counter[$z][$i]++;}
           if (($C_counter[$z][$i]+$c_counter[$z][$i])>0)
              {
              $signal[$z][$i] =
              int(
              ($C_counter[$z][$i]/($C_counter[$z][$i]
               +$c_counter[$z][$i]))*100);
              }
	   }
    }
}  

print"\n";
print "Writing to (100 x ",$laenge,") matrix. Please wait";	

#--------------------------------------------------------------------
#transfer values into 3D-matrix:
for ($z = $frame/2+1; $z<$laenge-$frame/2; $z=$z+$stepsize)
#for ($z = 0; $z<=$laenge; $z=$z+$stepsize) 
    {
    for ($p = 0; $p<=100; $p++)
	{
	 $signal_dens[$z][$p] = 0; #initialize matrix
         }
     }
     
for ($z = $frame/2+1; $z<$laenge-$frame/2; $z=$z+$stepsize)
#for ($z = 0; $z<=$laenge; $z=$z+$stepsize) 
    {
    for ($i = 0; $i<=$maxclones; $i++) 
        {
	unless ($signal[$z][$i] == 111)
	{
	for ($p = 0; $p<=100; $p++)
	    {
	    if ($p == $signal[$z][$i]) {$signal_dens[$z][$p]++};
            #print $signal_dens[$z][$p],"\n";
	    }
	}
	}
    print".";
    }

#-----------------------------------------------------------------------
# create a new image:
$im = new GD::Image($laenge+$offset*2,$spreiz*100+$offset*4);

print "\nCreating image. Please wait";  

# allocate some colors
$white =  $im->colorAllocate(255,255,255);
$green =  $im->colorAllocate(0,180,0);
$grey  =  $im->colorAllocate(220,220,220);
$black =  $im->colorAllocate(0,0,0);       
$blue =   $im->colorAllocate(0,0,255);
$c3 = $im->colorAllocate(102,104,255);
$c4 = $im->colorAllocate(47,204,253);
$c5 = $im->colorAllocate(0,255,48);
$c6 = $im->colorAllocate(0,155,5);
$c7 = $im->colorAllocate(255,255,0);
$c8 = $im->colorAllocate(255,200,0);
$c9 = $im->colorAllocate(255,51,0);
$c10= $im->colorAllocate(255,0,0);
if ($choose==1) {$connectlineclr = $grey}
else {$connectlineclr = $white};
$im->interlaced('true');

#plot CpG density:
for ($z = $frame/2+1; $z<$laenge-$frame/2; $z=$z+$stepsize)

#for ($z = 0; $z<=$laenge; $z=$z+$stepsize) 
    {
    for ($p = 0; $p<=100; $p++)
	{
	unless ($signal_dens[$z][$p] == 0)
	       {
	       if ($signal_dens[$z][$p]==1) {$color = $black};
	       if ($signal_dens[$z][$p]==2) {$color = $blue};
	       if ($signal_dens[$z][$p]==3) {$color = $c3};
        if ($signal_dens[$z][$p]==4) {$color = $c4};
        if ($signal_dens[$z][$p]==5) {$color = $c5};
	       if ($signal_dens[$z][$p]==6) {$color = $c6};
	       if ($signal_dens[$z][$p]==7) {$color = $c7};
	       if ($signal_dens[$z][$p]==8) {$color = $c8};
	       if ($signal_dens[$z][$p]==9) {$color = $c9};
	       if ($signal_dens[$z][$p]>=10) {$color = $c10};
	       unless ($pre_x == 0)
	       {$im->line($pre_x,$pre_y,$offset+$z,$offset+($spreiz*(100-$p)),$connectlineclr)};
	       $im->arc($offset+$z,$offset+($spreiz*(100-$p)),5,5,0,360,$color);
	       $pre_x = $offset+$z;
	       $pre_y = $offset+($spreiz*(100-$p));
	       }
        }
     }

	
#frame:
$im->rectangle($offset,$offset,$offset+$laenge,$offset+$spreiz*100,$black);

#scale:
$abschnitt = int($laenge/250);
for ($z=0; $z<=$abschnitt; $z++)
{
$im->line($offset+$z*250,$offset+$spreiz*100,$offset+$z*250,$offset+$spreiz*100+10,$black);
$im->string(gdSmallFont,$offset+$z*250,$offset+$spreiz*100+15,($z*250+1)." bp",$black);
}
$im->string(gdSmallFont,$offset-30,$offset-5,"100%",$black);
$im->string(gdSmallFont,$offset-30,$offset+10,"5mCpG",$black);
$im->string(gdSmallFont,$offset-30,$offset+$spreiz*50-5," 50%",$black);
$im->string(gdSmallFont,$offset-30,$offset+$spreiz*100-5,"  0%",$black);

$im->dashedLine($offset,$offset+$spreiz*25,$offset+$laenge,$offset+$spreiz*25,$black);
$im->dashedLine($offset,$offset+$spreiz*50,$offset+$laenge,$offset+$spreiz*50,$black);
$im->dashedLine($offset,$offset+$spreiz*75,$offset+$laenge,$offset+$spreiz*75,$black);

#figure legend:

$im->filledRectangle($offset,$offset+$spreiz*100+30,$offset+400,$offset+$spreiz*100+120,$grey);


$path = `pwd`;
chomp $path;
if (!defined $opt_n)
   {
    ($base = $path) =~ s#.*[/:]##;    #from path name end to last / or :
    }

else {
     $opt_n =~ tr/a-zA-Z0-9_\.\///cd;  #UNIX
     $base  =  $opt_n;
     }

if ($base !~ /\//) {$base =  substr ($base,0,28);}        #UNIX
unless ($base =~ /\.png$/i) {$base = $base.".distr.png";}

$im->string(gdSmallFont,$offset+5,$offset-45,"Sequences analyzed in",$black);
$im->string(gdMediumBoldFont,$offset+5,$offset-25,$path,$blue);
$im->string(gdSmallFont,$offset+15,$offset+$spreiz*100+40,
"window size            = ".$frame,$green);
$im->string(gdSmallFont,$offset+15,$offset+$spreiz*100+60,
"shift                  = ".$stepsize,$green);
$im->string(gdSmallFont,$offset+15,$offset+$spreiz*100+80,
"vertical grafic spread = ".$spreiz,$green);
$im->string(gdSmallFont,$offset+15,$offset+$spreiz*100+100,
"Sequences analyzed     = ".$maxclones,$green);

$im->string(gdSmallFont,$offset+220,$offset+$spreiz*100+40,"color coding of frequencies:",$black);
$im->string(gdSmallFont,$offset+220,$offset+$spreiz*100+60,"1",$black);
$im->string(gdSmallFont,$offset+230,$offset+$spreiz*100+60,"2",$blue);
$im->string(gdSmallFont,$offset+240,$offset+$spreiz*100+60,"3",$c3);
$im->string(gdSmallFont,$offset+250,$offset+$spreiz*100+60,"4",$c4);
$im->string(gdSmallFont,$offset+260,$offset+$spreiz*100+60,"5",$c5);
$im->string(gdSmallFont,$offset+270,$offset+$spreiz*100+60,"6",$c6);
$im->string(gdSmallFont,$offset+280,$offset+$spreiz*100+60,"7",$c7);
$im->string(gdSmallFont,$offset+290,$offset+$spreiz*100+60,"8",$c8);
$im->string(gdSmallFont,$offset+300,$offset+$spreiz*100+60,"9",$c9);
$im->string(gdSmallFont,$offset+310,$offset+$spreiz*100+60,">10",$c10);

open (BILD, ">$base") || die ("Cannot open file '$base'.");
binmode STDOUT;      # make sure to write to a binary stream
print BILD $im->png; # Convert the image to png
close BILD            || die ("Cannot close open file '$base'.");;

print "\npng-image written into file '$base'.\n";
}
##############################################################################
# SUBROUTINES                                                                #
# part of the Standard File Package Utility for MacPerl 4.1.1                #
#                                                                            #
# 1994.01.05 v4.1.1 Matthias Neeracher <neeri@iis.ee.ethz.ch>                #
#  Minor changes to reflect future plans for standard file support.          #
#                                                                            #
# 1993.10.27 v1.2	wm                                                         #
#	Change the calling syntax to adopt the 4.1.0 release.                      #
#                                                                            #
# 1993.10.19 v1.1	wm                                                         #
#	convert for 4.1b6                                                          #
#                                                                            #
# 1993.8.10  V1.0                                                            #
#     Watanabe Maki (Watanabe.Maki@tko.dec.com)                              #
#                                                                            #
##############################################################################
# Name
#    PutFile/GetNewFile
# Syntax
#    $filename = &PutFile($prompt [, $default]);
#    $filename = &GetNewFile($prompt [, $default]);
# Description
#    Query a new file name to user by Standard File Dialog Box.
#    $prompt is a prompt sting on the dialog box.
#    $default is a default file name.
#
#  sub PutFile {
#      local($prompt, $default) = @_;
#      
#      &MacPerl'Choose(
#          &GUSI'AF_FILE,         # domain
#          0,                     # type
#          $prompt,               # prompt
#          "",                    # constraint
#          &GUSI'CHOOSE_NEW + ($default ?
#  &GUSI'CHOOSE_DEFAULT : 0),    
#  		  								 # flag 
#          $default               # default filename
#          );
#  }
######
# Name
#    GetFolder
# Syntax
#    $foldername = &GetFolder($prompt [, $default]);
# Description
# Query a folder name to user by Standard File Dialog Box.
# $default is the default dialog
#
#  sub GetFolder {
#  	local($prompt, $default) = @_;
#  	
#      &MacPerl'Choose(
#          &GUSI'AF_FILE,          # domain
#          0,                      # type
#          $prompt,                # prompt
#          "",                     # constraint
#          &GUSI'CHOOSE_DIR + ($default ?
#  &GUSI'CHOOSE_DEFAULT : 0),
#  		                          # flag
#  		  $default
#          );
#  }