#!/usr/bin/perl

# Codigo perl para generar los datos de la gráfica del potencial DOPE
# Archivo original descargado de http://www.salilab.org/modeller/wiki/Plotdope
# Archivo modificado por Alfonso Méndez Tenorio (2007)
# Las columnas indicadas en la subrutina profiletodat fueron modificadas a un nuevo valor
# Se elimino la parte que generaba el script para gnuplot ya que no funcionaba y
# posiblemente estaba incompleto.
# En su lugar se genera ahora un archivo tipo "Comma Separated Values (CSV)" que puede
# ser importado rapidamente por Excel u otros programas.
# En el archivo se indican la posición y los valores DOPE para el molde (template) y el modelo (model).

$save_aln_sep = $/;

#Obtain the alignment file and extract the sequences into respective variables using their code

print "Enter Alignment filename (PIR format) [Example: TvLDH-1bdmA.ali] : ";

$alnfname = <STDIN>;

unless(open(ALNFNAME, $alnfname))
{
                print "Alignment file not found. Please check if file exists or if the filename is correct\n";
                exit;
}       

print "Enter code of the Target (Model) sequence [Example: TvLDH]: ";
$targcode = <STDIN>;
chomp $targcode;

print "Enter code of the Template structure [Example: 1bdmA]: ";
$tempcode = <STDIN>; 
chomp $tempcode;

$/ = undef;

$alignment = <ALNFNAME>;

#Here we split the alignment using the >P1 found in all the PIR sequences, so that the two sequences can be obtained in individual #variables. check_aln performs this job. Once obtained, the sequence alone is extracted using the parse_aln subroutine.

@alnseqs = split(/>P1;/, $alignment);


$dtargseq = check_aln(\@alnseqs,\$targcode);

@ftargseq = parse_aln($dtargseq, $targcode);

$dtempseq = check_aln(\@alnseqs,\$tempcode);
@ftempseq = parse_aln($dtempseq, $tempcode);

#The gap positions are identified using this subroutine

@targgaps = calc_gaps(@ftargseq);

print "Target gaps are at @targgaps\n";

@tempgaps = calc_gaps(@ftempseq);

print "Template gaps are at @tempgaps\n";

#The energy profiles contain some extra information. only the energy values of the residues are required. The subroutine profiletodat() #performs this function 
$/ = $save_aln_sep;

#Check if the profile files are already present. Else ask for a filename.
$targprofile = "$targcode.profile";

if(!open(TARGPROFILE, $targprofile))
{
        print "Enter name of Target profile [Example: TvLDH.profile]: ";
        $targprofile = <STDIN>;
}

open(TARGPROFILE, $targprofile);
@targprofile = <TARGPROFILE>;


$tempprofile = "$tempcode.profile";

if(!open(TEMPPROFILE, $tempprofile))
{
        print "Enter name of Template profile [Example: 1bdmA.profile]: ";
        $tempprofile = <STDIN>;
}

open(TEMPPROFILE, $tempprofile);

@tempprofile = <TEMPPROFILE>;


@targdat = profiletodat(@targprofile);

@tempdat = profiletodat(@tempprofile);

#The gaps in the sequences do not contain energy values. Here we modify the dat file(previous step) to 
@mtargdat = modify_dat(\@targdat, \@targgaps);

@mtempdat = modify_dat(\@tempdat, \@tempgaps);

##GnuPlot does not accept arrays. Hence we save the two dat files for plotting
$gtargdat = "$targcode.dat";

open(GTARGDAT, ">$gtargdat");

print GTARGDAT @mtargdat;

$gtempdat = "$tempcode.csv";

open(GTEMPDAT, ">$gtempdat");

#print GTEMPDAT "@mtempdat";

#Este código se agregó para reunir los datos en un archivo CSV

my $nodatos = @mtempdat;

my $itemp, $itarg, $itempprev, $itargprev;

print GTEMPDAT "0,template,model\n";

$itempprev = 0;
$itargprev = 0;

for (my $i=1; $i<$nodatos; $i++){
    $itemp = $mtempdat[$i];
    $itarg = $mtargdat[$i];
    chomp $itemp;
    chomp $itarg;
    #Cuando hay un hueco, repite el valor anterior para que la gráfica no se distorcione
    #Quizas no es la estrategia adecuada.
    if ($itemp eq '?') {$itemp = $itempprev};
    if ($itarg eq '?') {$itarg = $itargprev}; 
    print GTEMPDAT "$i,$itemp,$itarg\n";
    $itempprev = $itemp;
    $itargprev = $itarg;
}

close GTEMPDAT;

#Final del codigo adicionado


#El código siguiente se eliminó porque no funciona muy bien
#Finally plot these values. plot_dope calls gnuplot and plots the data. 


#$dopeout = "dope.png";

#open(DOPEOUT, $dopeout);

#if( -e DOPEOUT)
#{
#        print "$dopeout exists\n";
#        $dopeout = check_dope();
#}



#plot_dope($gtargdat, $gtempdat, $dopeout);


###############################SUBROUTINES###################
#check alignment and send sequences for parsing
sub check_aln
{
        my($alnseqs, $code)= @_;
        
        foreach $seq(@$alnseqs)
        {
                if($seq =~ /$$code/)
                {
                        return $seq;
                }
        }
}


#Parse alignment files to obatain sequence alone
sub parse_aln
{
        
        my($seq, $code) = @_;

        my @seq = $seq;
        my @parseq = ();

        foreach $line(@seq)
        {
                $line =~ s/\n//;

                if($line =~ /^$code/)
                {
                        $line =~ s/$code.*//;
                        $line =~ s/(sequence|structureX)[:](.*)*//;
                        $line =~ s/\n//s;
                        
                        push(@parseq, $line);
                }
        }
        return @parseq;
}



#calculates gaps given the sequence alone 
sub calc_gaps
{
        my ($sequence) = @_;
        $sequence =~ s/\n//;
        $lenseq = length($sequence);
#       print "Length of sequence is $lenseq\n";
        my $position = 0;
        my $gap = 0;
        my @gaps = ();

        for($position=0;$position < length($sequence);$position++)

        {
                $gap = substr($sequence, $position-1, 1);

                if($gap eq '-')
                {
                        push(@gaps, $position);
                }
        }
        return @gaps;

}

#To convert the draft profile to working dat file
sub profiletodat
{
        my(@profile) = @_;
        my $energy='';
        my @modprofile;

        foreach my $line(@profile)
        {
                
                if($line =~ /^\s*#/)
                {
                        next;
                }
                else
                {                       
                        $energy = substr($line, 484, 10);
                        push(@modprofile,$energy);
                }
        }

        return @modprofile;
}       


#To modify the dat files. Receives reference dat file and gap file containg positions of gaps

sub modify_dat
{
        my($dat, $gaps) = @_;

        foreach $position(@$gaps)
                {
                        for($datposition=0;$datposition < scalar @$dat; ++$datposition)
                        {
                                if($datposition == $position)
                                {
                                        splice(@$dat, $datposition, 0, "?\n");
                                }
                        }
                }

        return @$dat;
}


#Check for output file dope.png. If it exists, ask for a new filename
sub check_dope
{

        print "Should I overwrite (y/n): ";
        $check = <STDIN>;
        chomp $check;

        if($check eq 'y')
        {
                exit;
                
        }
        elsif($check eq 'n')
        {
                print "Enter new filename : ";
                $dopeout=<STDIN>;
                chomp $dopeout;
                
        }
        else
        {
                print "Please enter y/n";
                exit;
        }
        
        return $dopeout;

}


#Plotting the values using gnuplot
sub plot_dope
{
        my($targdat, $tempdat, $dopeout) = @_;
        
        my $gnuplot = "/usr/bin/gnuplot";
        
        my $terminal = "set terminal png\;\n";
        my $xlabel = "set xlabel \"Residue Index\"\;\n";
        my $ylabel = "set ylabel \"DOPE Score\"\;\n";
        my $missing = "set datafile missing \"?\"\;\n";
        my $output = "set output \"$dopeout\"\;\n";
        
        my $plot = "plot \"$targdat\" with lines, \"$tempdat\" with lines\;\n";

        my @plotset = ($terminal,$xlabel,$ylabel,$missing,$output,$plot);
#       print "@plotset\n";

        my $plotscript = "plotscript.txt";
        
        open(PLOTSCRIPT, ">$plotscript");
        
        print PLOTSCRIPT "@plotset";

        close (PLOTSCRIPT);
        
        my $finalscript = "plotscript.txt";
        
#       open(FINALSCRIPT, $finalscript);
#       
#       @gnuplotscript = <FINALSCRIPT>;
        
        `$gnuplot $finalscript`;

        print "The DOPE plot is in $dopeout\n";
}
        




