#!/usr/bin/perl -w

##########################################################################################
## program: 		calculateGibbs.pl
## author:			Bruno Contreras-Moreira, with help from Heladia Salgado & Osbaldo Resendis
## what it does: 
## calculates Gibbs destabilization energies for DNA duplexes based on their sequences sliding windows and optionally locates
## known sequence sites in the context of the calculated energy profiles
## follows the dinucleotide nearest-neighbour algorithm proposed by Breslauer et al.(1986) PNAS, 83:3746-3750. 
##
## arguments: 	calculateGibbs.pl -i InputSeqFile (compulsory) (run ./calculateGibbs.pl for more)
##
##
## returns:			1) energy profiles for each sequence parsed at InputSeqFile
##
##
##########################################################################################

use strict;
use Getopt::Std;
$| = 1;

####################### predefined constants & variables #################################

my $K = 1.987; # cal/mol.K Boltzmann 
my $defT = 310; # default temperature(K)
my $defL = 15; # default window length
my $l = $defL; # window length
my $t = $defT; # temperature

my %thermoTable; 
# dinucleotide energies obtained from Santalucia et al(1996) Biochemistry: 35:3555­3562
# measured at NaCl 1M, 37C & pH=7
# H in kcal/mol
# S in eu
# G in kcal/mol, not used, calculated from H & S with respect to temperature t
$thermoTable{'AA/TT'}{'H'} = 8.4;
$thermoTable{'AA/TT'}{'S'} = 23.6;
$thermoTable{'AA/TT'}{'G'} = 1.02;
$thermoTable{'AT/TA'}{'H'} = 6.5;
$thermoTable{'AT/TA'}{'S'} = 18.8;
$thermoTable{'AT/TA'}{'G'} = 0.73;
$thermoTable{'TA/AT'}{'H'} = 6.3;
$thermoTable{'TA/AT'}{'S'} = 18.5;
$thermoTable{'TA/AT'}{'G'} = 0.60;
$thermoTable{'CA/GT'}{'H'} = 7.4;
$thermoTable{'CA/GT'}{'S'} = 19.3;
$thermoTable{'CA/GT'}{'G'} = 1.38;
$thermoTable{'GT/CA'}{'H'} = 8.6;
$thermoTable{'GT/CA'}{'S'} = 23.0;
$thermoTable{'GT/CA'}{'G'} = 1.43;
$thermoTable{'CT/GA'}{'H'} = 6.1;
$thermoTable{'CT/GA'}{'S'} = 16.1;
$thermoTable{'CT/GA'}{'G'} = 1.16;
$thermoTable{'GA/CT'}{'H'} = 7.7;
$thermoTable{'GA/CT'}{'S'} = 20.3;
$thermoTable{'GA/CT'}{'G'} = 1.46;
$thermoTable{'CG/GC'}{'H'} = 10.1;
$thermoTable{'CG/GC'}{'S'} = 25.5;
$thermoTable{'CG/GC'}{'G'} = 2.09;
$thermoTable{'GC/CG'}{'H'} = 11.1;
$thermoTable{'GC/CG'}{'S'} = 28.4;
$thermoTable{'GC/CG'}{'G'} = 2.28;
$thermoTable{'GG/CC'}{'H'} = 6.7;
$thermoTable{'GG/CC'}{'S'} = 15.6;
$thermoTable{'GG/CC'}{'G'} = 1.77;
my (%initGC,%initAT);
$initGC{S} = 5.9; # helix formation start costs
$initGC{G} = -1.82;
$initAT{S} = 9.0;
$initAT{G} = -2.8;

my $nullvalue = -9999; # use to pad sequence profiles (see calculateGibbsProfile)
my $min_n_of_sites_for_histogram = 5;

#######################  main function   ############################################################

## parsing command-line options ##

my $progname = "calculateGibbs.pl";
my (%opts,$parameters);
getopts('hsi:t:l:',\%opts);

if(($opts{'h'})||(scalar(keys(%opts))==0)) 
{ 
print <<HELP; # print program's help
NAME
	calculateGibbs.pl

AUTHORS
	Bruno Contreras-Moreira, Heladia Salgado, Sarath Chandra Janga & 
	Osbaldo Resendis-Antonio
	March 2005
	Computational Genomics Lab, Center for Genomic Sciences/UNAM
	Cuernavaca, México

DESCRIPTION
	Takes a DNA sequence and slides a window of length l calculating 
	the profile of Gibbs destabilization energy of a DNA duplex with 
	the same sequence. Follows the dinucleotide nearest-neighbour 
	algorithm proposed by Breslauer et al.(1986) PNAS, 83:3746-3750.

EXAMPLE
	./calculateGibbs.pl -i InputSeqFile -l 30 -t 273
        
FORMAT of InputSeqFile 
	name0 \\ AAAAAAAAAAAAAAA... \\ 
	name1 \\ CCCCCAAAAAACCCC... \\ 
	...
HELP

	print "USAGE: $progname [options]\n";
	print "\t-h \t this message\n";
	print "\t-i \t input DNA sequences file (compulsory)\n";
	print "\t-t \t temperature (K) (optional, by default t = 300)\n";
	print "\t-l \t window length (odd number, optional, by default l = 15)\nn";
	exit(-1); 
}

my ($infile) = $opts{'i'};
if	(!defined($opts{'i'})) 	{ die "# $progname : not valid input file\n"; }
if	(!-s $infile)	{ die "# $progname : input file $infile is empty\n"; }
if	(defined($opts{'l'})) { 	$l = $opts{'l'}; } 
if($l%2 == 0) { die "# $progname : -l accepts only odd numbers\n"; }
if($l <= 1) { die "# $progname : -l accepts only numbers > 1\n"; }
if	(defined($opts{'t'})) { 	$t = $opts{'t'}; }

$parameters = sprintf("-i %s -l %d -t %d",$infile,$l,$t);

print "# $progname : parameters = $parameters\n\n";

## parse input file and compute energies for each sequence ##

my ($name,$seq,$site,@siteCoords,$sbeg,$send,$act,$TF,$siteZ,$siteG,$n_of_validZ,$s,$seqlength,@gibbsProfile,%STATS,%TFSTATS);

open(SEQ, $infile) || die "# $progname : cannot open input $infile\n";
while(<SEQ>)
{
	next if(/^#/); # skip commented lines
	chomp;
	$_ =~ s/ //g;	
	
	($name,$seq,$site) = split(/\\/,$_); 
	$seq = uc($seq);
	$seqlength = length($seq);
	
	print ">$name ",length($seq)," nts "; # more things are printed at calculateGibbsProfile
	
   @gibbsProfile = &calculateGibbsProfile($seq, $l, $t, 1);
	
	for($s=0;$s<length($seq);$s++) 
	{	
		print "$gibbsProfile[$s][0] $gibbsProfile[$s][1] $gibbsProfile[$s][2] $gibbsProfile[$s][3]\n";	
	}
}
close(SEQ);

####################### other functions and procedures   ############################################
#####################################################################################################

sub calculateGibbsProfile
{
   	my ($seq,$l,$t) = @_;
	my (@sequence,@profile) = (split(//,$seq),());
	my ($halfl,$i,$n,$lasti,$lastn,$totali,$prob,$DeltaH,$DeltaS,$DeltaG,$bpKey,$motif);
	my ($totalprob,$aveDeltaG,$stddevDeltaG) = (0.0,0.0,0.0);
	
	$halfl = ($l-1)/2;
	$lasti = length($seq) - $halfl - 1;
   	for($i=$halfl;$i<=$lasti;$i++) # calculate DNA destabilization energies within sliding windows of length l
	{
		$DeltaH = $DeltaS = $DeltaG = 0.0;
		$motif = "";
		$lastn = $i + $halfl - 1;
		for ($n=$i-$halfl;$n <= $lastn;$n++) 
		{
			$bpKey = $sequence[$n] . $sequence[$n+1] . "/" . complement($sequence[$n]) . complement($sequence[$n+1]);
			if(!$thermoTable{$bpKey}{'H'}){	$bpKey = reverse($bpKey);	}	
			$DeltaH += $thermoTable{$bpKey}{'H'};
			$DeltaS += $thermoTable{$bpKey}{'S'};
			$motif .= $sequence[$n];
		}
		$motif .= $sequence[$n];
		$DeltaH *= 1000;
		$DeltaG = $DeltaH  - ($t * $DeltaS);
		#$prob = exp(($DeltaG * -1) / ($K * $t)); # prob(DeltaG)
		$prob = sprintf("%2.0f",($DeltaG * -1) / ($K * $t));		
		
		# store motif data  
      		push(@profile,[($i,$motif,$DeltaG,$prob)]);
		
		# store motif values to calculate sequence statistics
		$totalprob += $prob;
		$aveDeltaG += $DeltaG;
   }
	
	# now carry out sequence statistics
	$totali = scalar(@profile);
	$aveDeltaG /= $totali;
	for($i=0;$i<$totali;$i++) {	$stddevDeltaG += abs($profile[$i][2]-$aveDeltaG);  	}
	$stddevDeltaG /= $totali;
	for($i=0;$i<$totali;$i++) # calculate Z scores for every motif within this sequence
	{	push(@{$profile[$i]}, sprintf("%1.2f",($profile[$i][2]-$aveDeltaG)/$stddevDeltaG ));	}
	
	# pad left side of sequence profile with null values 
	@profile = reverse @profile;
	for($i=$halfl-1;$i>=0;$i--) { push(@profile,[($i,"",$nullvalue,$nullvalue)]); }
	@profile = reverse @profile;
	
	# now pad the right side
	for($i=$lasti+1;$i<length($seq);$i++) { push(@profile,[($i,"",$nullvalue,$nullvalue)]); }
	
	# print mean and standard dev before leaving
	printf(" %1.0f %1.0f\n",$aveDeltaG,$stddevDeltaG);
	
	return @profile;
}

#####################################################################################################

sub complement 
{
	my $nucl = $_[0];
	if(!$nucl){ return "";}
	$nucl =~ tr/ATGC/TACG/;
	return $nucl;
}

#####################################################################################################