bindingSiteMap.pl

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

###############################################################################
###############################################################################
### Extracts from an all info noSql_parser.pl output each binding site hit  ###
### and adds the size of binding site to the output.                        ###
###############################################################################
###############################################################################

use strict;
use warnings;
use Getopt::Long;

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

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

&GetOptions(
    "help!"             => \$help,               # print this help
    "minChainBSsize:i"  => \$minChainBSsize,     # minimum size of a binding site that forms an inter chain interface. Default is 10.
    "minDomainBSsize:i" => \$minDomainBSsize,    # minimum size of a binding site that forms an inter domain interface. Default is 4.
    "zScoreCutOff:f"    => \$zScoreCutOff,       # minimum z-score for significant binding site matches. Default is 4.
) or die "\nTry \"$0 -h\" for a complete list of options\n\n";

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

# help
if ($help) {printHelp(); exit}
$| = 1;

##############################################################################
### SETTINGS #################################################################
##############################################################################
$minChainBSsize ||= 10;
$minDomainBSsize ||= 4;
$zScoreCutOff ||= 1.3;

##############################################################################
### 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\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\n", "-",$option,$row);

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

    $example  = "$0";
}
##############################################################################
sub readResultFile {
    (my $file) = @_;

    my $bs = ();
    open(F, "gzip -dc $file |") or die "\nERROR: Failed to open $file: $!\n\n";
    my %zscores = ();
    while(my $l = <F>) {
	chomp($l);
	# each line has info about Query binding site and Hit binding site
	(my $file, my $chainQ, my $resnQ, my $resiQ, my $pdbH, my $chainH, my $resnH, my $resiH, my $blosumScore, my $rmsd, my $z_score) = split(/\t/, $l);
	(my $d, my $orgQ, my $uniprotQ, my $pdbQ, my $pfamQ) = split(/\//, $file);
	my $key = "$uniprotQ\t$pdbQ\t$pfamQ\t$pdbH\t$resiQ";
	if(exists $zscores{$key}) {
	    $bs->{$uniprotQ}->{$pdbQ}->{$pfamQ}->{$pdbH}->{$resiQ} = "$resnQ\t$resiH\t$resnH\t$blosumScore\t$rmsd\t$z_score" if($z_score > $zscores{$key});
	    $zscores{$key} = $z_score;
	} else {
	    $bs->{$uniprotQ}->{$pdbQ}->{$pfamQ}->{$pdbH}->{$resiQ} = "$resnQ\t$resiH\t$resnH\t$blosumScore\t$rmsd\t$z_score";
	    $zscores{$key} = $z_score;
	}
    }
    close(F);
    return $bs;
}

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




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

if(@ARGV < 2) {
    print STDERR "\n$0 nonOverlapping_uniq_bindingSites.txt matchDir\n\n";
    exit;
}

my $matchDir = $ARGV[1];

my $bs;
while(my $l = <>) {
    chomp($l);
    #CALM_RAT 4ehqA 4.2.PF00036.120.146#M.1.PF00036.47.74 120,143,144 /mnt/nfs3/mnemo2/damian/structures/pdb_derived/all_pdb_chains/4e/4ehq.pdb1.1.A.pdb
    (my $uniprot, my $pdb, my $pfam, my $bsNos, my $file)  = split(/\t/, $l);
    my @bsNos = split(/\,/, $bsNos);
    (my $id, my $org) = split(/\_/, $uniprot);
    $pfam =~ s/\#/-/g; 
    my $chainId = substr($pdb, 4, 1);
    my $resultFile = "$matchDir/$org/$uniprot/$pdb/$pfam/$pfam"."_results.txt.gz";
    my $probisPredFile = "$matchDir/$org/$uniprot/$pdb/$pfam/$pfam"."_probisBs.txt.gz";
    if(!-e $resultFile) {
	$resultFile = "$matchDir/$org/$uniprot/$pdb/$pfam/$pfam$chainId"."_results.txt.gz";
	$probisPredFile = "$matchDir/$org/$uniprot/$pdb/$pfam/$pfam$chainId"."_probisBs.txt.gz";
    }
    if (-e $resultFile) {
	open(O, "| gzip -9 > $probisPredFile") or die "\nERROR: Failed to create probis prediction file $probisPredFile: $!\n\n";
	my $bs = &readResultFile($resultFile);
	foreach my $uniprotQ (sort keys %$bs) {
	    foreach my $pdbQ (sort keys %{$bs->{$uniprotQ}}) {
		foreach my $pfamQ (sort keys %{$bs->{$uniprotQ}->{$pdbQ}}) {
		    foreach my $pdbH (sort keys %{$bs->{$uniprotQ}->{$pdbQ}->{$pfamQ}}) {
			my $hit = "-";
			my $hitCount = 0;
			my $failCount = 0;
			my $seqQ = "-";
			my $seqH = "-";
			my $seqNoH = "-";
			my $seqNoQ = "-";
			my $blosumScore = "-";
			my $rmsd = "-";
			my $z_score = "-";
			foreach my $bsNo (@bsNos) {
			    if(exists $bs->{$uniprotQ}->{$pdbQ}->{$pfamQ}->{$pdbH}->{$bsNo}) {
				$hit = $bs->{$uniprotQ}->{$pdbQ}->{$pfamQ}->{$pdbH}->{$bsNo};
				(my $resnQ, my $resiH, my $resnH, $blosumScore, $rmsd, $z_score) = split(/\t/, $hit); 
				$seqQ = "" if($seqQ eq "-");
				$seqH = "" if($seqH eq "-");
				$seqNoH = "" if($seqNoH eq "-");
				$seqNoQ = "" if($seqNoQ eq "-");
 				$seqQ .= $resnQ;
 				$seqH .= $resnH;
 				$seqNoH .= "$resiH,";
 				$seqNoQ .= "$bsNo,";
				$hitCount++;
			    } else {
				$failCount++;
			    }
			}
			$seqNoH =~ s/,$//; 
			$seqNoQ =~ s/,$//; 
			my $total = $hitCount + $failCount;
			# first check the hit has high enough zScore
			if($z_score ne "-" and $z_score >= $zScoreCutOff) {
			    # then check that hit is large enough
			    if(($pfamQ =~ /PF/ and length($seqH) >= $minDomainBSsize) or ($pfamQ !~ /PF/ and length($seqH) >= $minChainBSsize)) {
				# header
				# queryUniprotId queryBsId queryPDBid probisPDBid queryBsNoFull queryBsNo probisBsNo queryBsSeq probisBsSeq probisBlosumScore probisRmsd probisZscore probisBsSize probisAligned probisAlignedGapRatio probisGap
				print O "$uniprotQ\t$pfamQ\t$pdbQ\t$pdbH\t$bsNos\t$seqNoQ\t$seqNoH\t$seqQ\t$seqH\t$blosumScore\t$rmsd\t$z_score\t$total\t$hitCount\t".
				($total == 0 ? 0 : sprintf("%.2f",$hitCount/$total))."\t$failCount\n";
			    }
			}
		    }
		}
	    }
	}
    } else {
	print STDERR "WARNING: Could not find result file $resultFile!\n";
    }
    close(O);
}