#!/usr/local/bin/perl -w
#
################################################################################
# Requirements                                                                 #
################################################################################
#
require 5.002;
#
# External References
#
use GD;
use Carp;
use FileHandle;
use File::Path;
use File::Basename;
use Getopt::Std;

use vars qw($opt_d $opt_n $opt_h);
#
################################################################################
# Main program global variables                                                #
################################################################################
#
my $next        = 0;
my $position    = -1;
my $seqname        ;
my @fcg_circles    ;
my @fcaa_circles   ;
my @fcag_circles   ;
my @fctg_circles   ;
my @fcta_circles   ;
my @files          ;
my $file           ;
my $seq            ;
my $laenge      = 0;
my $max_clones  = 0;
my $i           = 0;
my $f_circle    = 0;
my $h_circle    = 0;
my $fcag_circle = 0;
my $hcag_circle = 0;
my $fctg_circle = 0;
my $hctg_circle = 0;
my $fcta_circle = 0;
my $left_offset    ;
my $startline   = 0;
my $suche    = "-a";
my $base           ;
my $path           ;
my $maxclones   = 0;
#
################################################################################
# Grafic output default values                                                 #
################################################################################
#
my $stepdown    =   10;
my $offset_h    =   50;
my $offset_v    =  110;
my $max_x_size  = 1800;

################################################################################
# MAN page                                                                     #
################################################################################
my $MAN = <<MAN;
NAME

     Plot_CpNpGs_to_png.pl: Program converts methylation patterns from Seq1 files 
     into a png-image

AVAILABILITY

     Requires perl 5.002

DESCRIPTION

     The program serve for the visualization of methylation patterns of 
     individual sequences. 
     The program "PlSeq2png" was adapted to non-symmetric methylation in plants
     and uses a set of red, blue, green, black and yellow circles to represent 
     methylated CpGs, CpApAs, CpApGs, CpTpGs and CpTpA, respectively.
     Unmethylated Cs are not shown.
     Output is a file in png-format ("working_directory".png)

SYNOPSIS

     Plot_CpNpGs_to_png.pl [-f path_to_seq1-files] [-n outfile] [-h]

OPTIONS

    -f string Path name to seq1-files
    -n string file name of png-image for output
    -h        Prints this message
 
AUTHOR(S)

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

CHANGE LOG

     Nov.   1st, 1998 CG   Initial version
     Sep.  30th, 1999 CG   minor changes
     Apr.  10th, 2000 CG   PNG version for UNIX

MAN
################################################################################
# main program                                                                 #
################################################################################
#
# Command line parsing 
#
getopts('d:n:h') || croak($MAN);
#
# Check command line arguments
#
croak($MAN) if ( $opt_h || ! defined($opt_d) || $opt_d eq "" );

#write all .seq1-files into a list
chdir $opt_d            || croak ("Cannot find directory '$opt_d'.");
@files = <*.seq1>;

$maxclones = $#files;
#
# finish program when no files *.seq1 found
#
if ($maxclones <0) {print "Sorry. No files available! Please try option -h.\n"}
else
{
print "Creating image";
# create a new image

$im = new GD::Image($max_x_size,$maxclones*$stepdown+150);

# allocate some colors
$white = $im->colorAllocate(255,255,255);
$black = $im->colorAllocate(0,0,0);       
$red   = $im->colorAllocate(255,0,0);      
$blue  = $im->colorAllocate(0,0,255);
$green = $im->colorAllocate(0,200,0);
$yellow =$im->colorAllocate(255,200,0);

# make the background interlaced
$im->interlaced('true');

$| = 1;                          #print immediately


foreach $filelist(@files)
	{
	$seq = "";
	open (SEQUENCE, "<$filelist");  #consequtive transfer of sequence-files
	print (".");                    #show program is running
	while (<SEQUENCE>)
	    {
	    $file = $_;
	    chomp $file;                #read single line, cut Newline
	    #
	    # recognice ">" of Fasta-name, max.length 15 letters
	    # 
     if ($file =~ /^>/) {$seqname = substr($file, 1, 15)}
	    #
	    # assemble sequence from the lines
	    #
	    unless ($file =~ /^>/){$seq = $seq . $file}
	    }
	close SEQUENCE;
	
	 @fcg_circles  = ();             # set lists back to zero
	 @fcaa_circles = ();
	 @fcag_circles = ();
	 @fctg_circles = ();
	 @fcta_circles = ();
	 
	 while ($seq)
	     {
	     $next = $position + 1;
	     $position = index ($seq,"Cg",$next);
	     if ($position == -1) {last};
	     push (@fcg_circles, $position);
	     }
	 $position = -1;
	 while ($seq)
	     {
	     $next = $position + 1;
	     $position = index ($seq,"Caa",$next);
	     if ($position == -1) {last};
	     push (@fcaa_circles, $position);
	     }
	 $position = -1;
	 while ($seq)
	     {
	     $next = $position + 1;
	     $position = index ($seq,"Cag",$next);
	     if ($position == -1) {last};
	     push (@fcag_circles, $position);
	     }
	 $position = -1;
	 while ($seq)
	     {
	     $next = $position + 1;
	     $position = index ($seq,"Ctg",$next);
	     if ($position == -1) {last};
	     push (@fctg_circles, $position);
	     }
	 $position = -1;
	 while ($seq)
	     {
	     $next = $position + 1;
	     $position = index ($seq,"Cta",$next);
	     if ($position == -1) {last};
	     push (@fcta_circles, $position);
	     }
	#
	# find internal "-" and replace by dummy sign "_" 
	#
	$seq =~ s/\w{1}-+\w{1}/_/g;
        #
	# find first occurance of -g etc.
	#	
	if ($seq =~ /-g/) {$suche = "-g"};
	if ($seq =~ /-t/) {$suche = "-t"};
	if ($seq =~ /-c/) {$suche = "-c"};
	if ($seq =~ /-C/) {$suche = "-C"};
	if ($seq =~ /-n/) {$suche = "-n"};
	if ($seq =~ /-a/) {$suche = "-a"};
	$startline = index($seq,$suche);
	$startline = $startline+$offset_v;
	#
	# remove all "-" and "_" and determine sequence length 
	#
	$seq =~ s/-//g;
	$seq =~ s/_//g;
	$laenge = length($seq);
	#
        # write Fasta-name
	#
	$im->string(gdSmallFont,$startline-100,$offset_h+$stepdown-5,$seqname,$black);
	#
	# draw line representing the sequence fragment
	#
	$im->line($startline,
		  $offset_h+$stepdown,
		 ($startline+$laenge),
		  $offset_h+$stepdown,
		  $black);
	
	while (@fcaa_circles)
	{
	# draw blue circle for methylated CpApA
	$h_circle = shift(@fcaa_circles);
        $im->arc($offset_v+$h_circle,$offset_h+$stepdown,5,5,0,360,$blue);
	}
	
	while (@fcg_circles)
	{
	# draw red circle for methylated CpG
	$f_circle = shift(@fcg_circles);
        $im->arc($offset_v+$f_circle,$offset_h+$stepdown,5,5,0,360,$red);
	}
	
	while (@fcag_circles)
	{
	# draw green circle for methylated CpApG
	$fcag_circle = shift(@fcag_circles);
       
	$im->arc($offset_v+$fcag_circle,$offset_h+$stepdown,5,5,0,360,$green);
	}
	
	while (@fctg_circles)
	{
	# draw black circle for methylated CpTpG
	$fctg_circle = shift(@fctg_circles);
       
	$im->arc($offset_v+$fctg_circle,$offset_h+$stepdown,5,5,0,360,$black);
	}
	
	while (@fcta_circles)
	{
	# draw yellow circle for methylated CpTpA
	$fcta_circle = shift(@fcta_circles);
       
	$im->arc($offset_v+$fcta_circle,$offset_h+$stepdown,5,5,0,360,$yellow);
	}
      $stepdown = $stepdown+10;
  }
print "\n";

$path = `pwd`;
chomp $path;

$im->string(gdSmallFont,2,10,"Sequences analyzed in:",$black);
$im->string(gdMediumBoldFont,2,30,$path,$blue);

#legend:
$im->rectangle(490,5,762,45,$black);
$im->fill(500,35,$white);
$im->string(gdSmallFont,500,10,"5mCTA = o",$yellow);
$im->string(gdSmallFont,500,30,"5mCAA = o",$blue);
$im->string(gdSmallFont,600,10,"5mCAG = o",$green);
$im->string(gdSmallFont,600,30,"5mCTG = o",$black);
$im->string(gdSmallFont,700,10,"5mCG  = o",$red);

#x-axis:
$im->line($offset_v,$offset_h+$stepdown+20,$offset_v+1000,$offset_h+$stepdown+20,$black);
$im->line($offset_v,$offset_h+$stepdown+15,$offset_v,$offset_h+$stepdown+25,$black);
$im->line($offset_v+250,$offset_h+$stepdown+15,$offset_v+250,$offset_h+$stepdown+25,$black);
$im->line($offset_v+500,$offset_h+$stepdown+15,$offset_v+500,$offset_h+$stepdown+25,$black);
$im->line($offset_v+750,$offset_h+$stepdown+15,$offset_v+750,$offset_h+$stepdown+25,$black);
$im->line($offset_v+1000,$offset_h+$stepdown+15,$offset_v+1000,$offset_h+$stepdown+25,$black);
$im->string(gdSmallFont,$offset_v,$offset_h+$stepdown+30,"1 bp",$black);
$im->string(gdSmallFont,$offset_v+250,$offset_h+$stepdown+30,"251",$black);
$im->string(gdSmallFont,$offset_v+500,$offset_h+$stepdown+30,"501",$black);
$im->string(gdSmallFont,$offset_v+750,$offset_h+$stepdown+30,"751",$black);
$im->string(gdSmallFont,$offset_v+1000,$offset_h+$stepdown+30,"1001",$black);

$im->string(gdSmallFont,2,$offset_h+$stepdown+50,"Sequence names truncated at 15 letters. Clones in alphabetical order.",$black);
#
# write image to file:
#
if (!defined $opt_n)
   {
    ($base = $path) =~ s#.*[/:]##;    #from path name end to last /
    }

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

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

open (BILD, ">$base");
binmode STDOUT;      # make sure to write to a binary stream
print BILD $im->png; # Convert the image to png
close BILD;
print "png-image written into file '$base'.\n";
}