atomPotDiff.pl

#!/usr/bin/perl -w
# Author: Abdullah Kahraman
# Date: 07.03.2004

###############################################################################
###############################################################################
### This script calculates the square root distance between two spherical   ###
### hormonics coefficient file.                                             ###
###############################################################################
###############################################################################

use strict;
use Getopt::Long;

my (
    # variable for parameters which are read in from commandline
    $help,
   );

##############################################################################
### read all needed parameters from commandline ##############################

&GetOptions(
    "help!"      => \$help,     # print this help
) or die "\nTry \"$0 -h\" for a complete list of options\n\n";

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

# help
if ($help) {printHelp(); exit}

##############################################################################
### SETTINGS #################################################################
##############################################################################

# Get Author Information 
my $author;
my $sg= "finger  " . `whoami`;
my @authorName = `$sg`;
@authorName = grep /life/, @authorName;
foreach (@authorName){
  if (/life:\s+(\S+\s+\S+)/){
    $author =  $1;
    last;
  }
}
my $dbsubmitter = $author;

# Get Date
my @submitDate = split(" ", `date`);
my $dbsubmitDate = $submitDate[2] . "-" . $submitDate[1] . "-" .
                    $submitDate[5];

##############################################################################
### SUBROUTINES ##############################################################
############################################################################## 

###############################################################################
sub printHelp {
###############################################################################
    # prints a help about the using and parameters of this scripts 
    # (execute if user types commandline parameter -h)
    # param:  no paramaters
    # return: no return value

    my (
	$usage,
	$sourceCode,
	@rows,
	$row,
	$option,
	$scriptInfo,
	$example,
       );

    $usage = "$0 -in=<path> -field=[string] -stop=[integer] " .
             "-output [stdout|database] -help -list -verbose\n";


    print "\nUsage: " .  $usage . "\n";

    print "Valid options are:\n\n";
    open(MYSELF, "$0") or
      die "Cannot read source code file $0: $!\n";
    $sourceCode .= join "", <MYSELF>;
    close MYSELF;
    $sourceCode =~ s/^.+?\&GetOptions\(\n//s;
    $sourceCode =~ s/\n\).+$//s;
    @rows = split /\n/, $sourceCode;
    foreach $row (@rows){
        $option = $row;
	$option =~ s/\s+\"//g;
	$option =~ s/\"\s.+\#/\t\#/g;
	$option =~ s/=./\t<value> [required]/;
	$option =~ s/:./\t<value> [optional]/;
	$option =~ s/!/\t<non value> [optional]/;

	$row =~ s/^.*//;
	print "\t";
	printf("%-1s%-30s%-30s%-20s\n", "-",$option,$row);

    } # end of foreach $row (@rows)
    print "\n";
    print "Options may be abreviated, e.g. -h for --help\n\n";

    $example  = "$0 -in=../../databases/locuslink/locuslink_031007.txt ";
}

##############################################################################
### END OF SUBROUTINES########################################################
############################################################################## 


############
### MAIN ###
############

if ($ARGV[0] && $ARGV[1]) {
  my %lmPair;
  my $dist;

  open(F1, "<$ARGV[0]") or 
    die "Can't open file \"$ARGV[0]\", $!\n";

  my @posNeg1;
  my $i=0;
  while(<F1>) {
    next if(/MAP /);
    if(/^ATOM/ or /^HETATM/){
      my $atomNo   = substr($_,  6, 5);
      $atomNo =~ s/\s//g;
      my $tempFact = substr($_, 60, 6);
      $tempFact =~ s/\s//g;
      $lmPair{$atomNo} = $tempFact;
      if($tempFact>=0){$posNeg1[$i++]=+1;}
      else{$posNeg1[$i++]=-1;}

    }
  }
  close(F1);

  open(F2, "<$ARGV[1]") or 
    die "Can't open file \"$ARGV[1]\", $!\n";

  if(defined $ARGV[2]){
    my $count;
    while(<F2>) {
      next if(/MAP /);
      if(/^ATOM/ or /^HETATM/){
	$count++;
	my $atomNo   = substr($_,6,5);
	$atomNo =~ s/\s//g;
	my $tempFact = substr($_,60,6);
	$tempFact =~ s/\s//g;
	$dist += ($tempFact - $lmPair{$atomNo})**2;
      }
    }

    $dist = sqrt($dist);
    print "$ARGV[0]\t$ARGV[1]\t$dist\n";
  }
  else{
    my @posNeg2;
    $i=0;
    while(<F2>) {
      next if(/MAP /);
      if(/^ATOM/ or /^HETATM/){
	my $atomNo   = substr($_,6,5);
	$atomNo =~ s/\s//g;
	my $tempFact = substr($_,60,6);
	$tempFact =~ s/\s//g;
	if($tempFact>=0){$posNeg2[$i++]=+1;}
	else{$posNeg2[$i++]=-1;}
      }
    }
    close(F2);
    my $no=$i;
    my $c=0;
    for(my $j=0; $j<=$#posNeg1; $j++){
	if(!defined $posNeg2[$j]){
	    print STDERR "$j\t$ARGV[0]\t$ARGV[1]\n";
	}
	else{
	    $c++ if($posNeg1[$j] != $posNeg2[$j]);
	}
    }
    print "$ARGV[0]\t$ARGV[1]\t".($c/$no)."\n";
  }
}

else {
  die "\nTry \"$0 -h\" for a complete list of options\n\n";
}

## end of main
## end of script