Plotdope

Plotdope is a perl script to plot the energy values obtained as output of the model evaluation step of MODELLER.

Often, there are gaps in the alignment of the target and template. In the model evaluation step, the energy values at these positions are not computed. To compare the energy profiles of two structures, the positions of these residues should match. This script obtains the alignment and the energy profiles as input and then introduces gaps accordingly in the parsed profile files and then plots the resulting data using Gnuplot

/!\ Please Note: This script is useful only for an alignment containing two sequences (target and a template.). Also complete error checks are not yet implemented.

Updated 11/7/07 by James Haven: add some error checking and extract correct energy values.

#!/usr/bin/perl -w


$save_aln_sep = $/;


#Open the alignment file 
print "Enter Alignment filename (PIR format) [Example: TvLDH-1bdmA.ali] : ";
$alnfname = <STDIN>;
open(ALNFNAME, $alnfname) or die "Alignment file not found. Please check if file exists or if filename is correct\n";



#Get the Target Code
print "Enter code of the Target (Model) sequence [Example: TvLDH]: ";
$targcode = <STDIN>;
chomp $targcode;



#Get the Template Code
print "Enter code of the Template structure [Example: 1bdmA]: ";
$tempcode = <STDIN>; 
chomp $tempcode;


$/ = undef;


#Store alignment data
$alignment = <ALNFNAME>;
close 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);



#extract TARGET sequence
$dtargseq = check_aln(\@alnseqs,\$targcode);
@ftargseq = parse_aln($dtargseq, $targcode);



#extract TEMPLATE sequence
$dtempseq = check_aln(\@alnseqs,\$tempcode);
@ftempseq = parse_aln($dtempseq, $tempcode);



#extract TARGET gap positions
@targgaps = calc_gaps(@ftargseq);
print "Target gaps are at @targgaps\n";



#extract TEMPLATE gap positions
@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 TARGET 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) or die "Can't open Target profile: $targprofile\n";
@targprofile = <TARGPROFILE>;
close TARGPROFILE;



#Check if the TEMPLATE profile files are already present. Else ask for a filename.
$tempprofile = "$tempcode.profile";
if(!open(TEMPPROFILE, $tempprofile))
{
        print "Enter name of Template profile [Example: 1bdmA.profile]: ";
        $tempprofile = <STDIN>;
}
open(TEMPPROFILE, $tempprofile) or die "Can't open Template profile: $tempprofile\n";
@tempprofile = <TEMPPROFILE>;
close TEMPPROFILE;



#Parse TARGET Profile file
@targdat = profiletodat(@targprofile);



#Parse TEMPLATE Profile file
@tempdat = profiletodat(@tempprofile);



#The gaps in the sequences do not contain energy values. Here we modify the dat file(previous step) to 
#Correct TARGET energy data for gaps (insert ? in gapped positions) => these sites won't be plotted
@mtargdat_old = modify_dat(\@targdat, \@targgaps);
for (my $i=0; $i<scalar(@mtargdat_old)-1; $i++) {
	push @mtargdat, $i."\t".$mtargdat_old[$i]."\n";
}



#Correct TEMPLATE energy data for gaps (insert ? in gapped positions) => these sites won't be plotted
@mtempdat_old = modify_dat(\@tempdat, \@tempgaps);
for (my $i=0; $i<scalar(@mtempdat_old)-1; $i++) {
	push @mtempdat, $i."\t".$mtempdat_old[$i]."\n";
}



##GnuPlot does not accept arrays. Hence we save the two dat files for plotting
#Write TARGET dat (see previous block) to file
$gtargdat = "$targcode.dat";
open(GTARGDAT, ">$gtargdat");
print GTARGDAT @mtargdat;



#Write TEMPLATE dat (see previous block) to file
$gtempdat = "$tempcode.dat";
open(GTEMPDAT, ">$gtempdat");
print GTEMPDAT "@mtempdat";



#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 = $seq;
                }
        }
	if($return_seq eq '') {die "Can't find code: $$code in alignment; subroutine check_aln failed!\n";}
	return $return_seq;
}



#Parse alignment files to obtain 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);
                }
        }

	if (scalar(@parseq) == 0) { die "Check your alignment (PIR) file; subroutine parse_aln failed!\n"}
        return @parseq;
}



#calculates gaps given the sequence alone 
sub calc_gaps
{
        my ($sequence) = @_;
        $sequence =~ s/\n//;
        my $position = 0;
        my $gap = 0;
        my @gaps = ();

        for($position=0;$position <= length($sequence)-1;$position++)

        {
                $gap = substr($sequence, $position, 1);

                if($gap eq '-')
                {
                        push(@gaps, $position+1);
                }
        }
	if(scalar(@gaps) == 0) { die "subroutine calc_gaps failed\n"}
        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
                {                       
			@a = split /\s+/, $line;
			$energy = $a[42];                        
			#$energy = substr($line, 436, 10);
			#print $energy,"\n";
                        push(@modprofile,$energy);
                }
        }
	if (scalar(@modprofile) == 0) { die "Something is wrong with an Energy Profile file; subroutine
                                             profiletodat failed\n"}     
	shift @modprofile;
	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, "?");
                                }
                        }
                }

        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 =~ /y/i)
        {               
		return $dopeout;
        }
        elsif($check =~ /n/i)
        {
                print "Enter new filename : ";
                $dopeout=<STDIN>;
                chomp $dopeout;
                return $dopeout;
        }
        else
        {
                print "Please enter y/n\n";
		$dopeout = check_dope();
        }
        
}



#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\" using 1:2 with lines, \"$tempdat\" using 1:2 with lines\;\n";

        my @plotset = ($terminal,$xlabel,$ylabel,$missing,$output,$plot);
#       print "@plotset\n";

        my $plotscript = "plotscript.txt";
        
        open(PLOTSCRIPT, ">$plotscript") or die "Can't open file $plotscript for writing\n";
        
        print PLOTSCRIPT @plotset;

        close (PLOTSCRIPT);
        
        my $finalscript = "plotscript.txt";
        
#       open(FINALSCRIPT, $finalscript);
#       
#       @gnuplotscript = <FINALSCRIPT>;
        
        `$gnuplot $finalscript`;

        print "The DOPE plot is in $dopeout\n";
}