createInteractions.pl

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

###############################################################################
###############################################################################
### creates interactions from binding site hits e.g. BLAST or PROBIS.       ###
###############################################################################
###############################################################################

use strict;
use warnings;
use Getopt::Long;
no autovivification;

my (
    # variable for parameters which are read in from commandline
    $help,
    $bindingSiteFile,
    $hIntFile,
    $nIntFile,
    $fusedFile,
    $probisFile,
    $isBlast,
    $isProbis,
    $isHHsearch,
    $isNeedle,
    $withParalog,
    $verbose,
   );

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

&GetOptions(
    "help!"          => \$help,             # print this help
    "bindingSite=s"  => \$bindingSiteFile,  # e.g. npdbVsEnsg_blast.tsv.gz
    "int=s"          => \$hIntFile,         # e.g. human_chain_interfaces_enspId_uniprotId_pdbId_ensgId.txt
    "nonInt=s"       => \$nIntFile,         # e.g. chain_interfaces_pdbID_uniprotAC_uniprotID_pfam_local_nonhuman.txt
    "fused:s"        => \$fusedFile,        # ignores PDB chains that are fused. e.g. fused_pdb.tsv
    "probis!"        => \$isProbis,         # use when binding site file was created by ProBIS e.g. all_probisBsSignif.txt.gz
    "probisFile:s"   => \$probisFile,       # use in combination with BLAST file on ProBIS binding sites e.g. all_probisBsSignif.txt.gz
    "blast!"         => \$isBlast,          # use when binding site file was created by BLAST e.g. npdbVsEnsg_blast.tsv.gz
    "hhsearch!"      => \$isHHsearch,       # use when binding site file was created by HHsearch e.g. human_ensembl_combined_hhsearchMatch.txt.gz
    "needle!"        => \$isNeedle,         # use when binding site file was created by needle e.g. probisVsEnsp_pairwise_needle.tsv.gz
    "paralog!"       => \$withParalog,      # predict also paraolog interactions using the target organism.
    "verbose!"       => \$verbose,          # output information on calculation progress on STDERR
) or die "\nTry \"$0 -h\" for a complete list of options\n\n";

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

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

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

##############################################################################
### 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 getInteractions {

    my %hIntsE = ();
    my %hIntsP = ();
    my %hIntsPE = ();
    my %hIntsG = ();
    open(F3, "zcat $hIntFile |") or die "\nERROR: Failed to unzip complex file $hIntFile: $!\n\n";
    while(my $l = <F3>) {
	next if($l =~ /^#/);
	chomp($l);
	my ($ensp1, $ensp2,
	    $ensg1, $ensg2,
	    $uniprot1, $uniprot2,
	    $pdb1, $pdb2,
	    $geneName1, $geneName2) = split(/\t/, $l);

	$geneName1 =~ s/\;$//;
	$geneName2 =~ s/\;$//;

	$geneName1 = "-" if($geneName1 eq "");
	$geneName2 = "-" if($geneName2 eq "");

	$hIntsE{$ensg1}->{$ensg2} = 1;
	$hIntsE{$ensg2}->{$ensg1} = 1;
	$hIntsP{$pdb1}->{$pdb2} = 1;
	$hIntsP{$pdb2}->{$pdb1} = 1;
	$hIntsPE{$pdb1}->{$ensg1} = 1;
	$hIntsPE{$pdb2}->{$ensg2} = 1;
	$hIntsG{$ensg1} = $geneName1;
	$hIntsG{$ensg2} = $geneName2;
    }
    close(F3);
    return (\%hIntsE, \%hIntsP, \%hIntsPE, \%hIntsG);
}

##############################################################################
sub getNonInteractions {

    my %nIntsP = ();
    open(F4, "zcat $nIntFile |") or die "\nERROR: Failed to unzip complex file $nIntFile: $!\n\n";
    while(my $l = <F4>) {
	next if($l =~ /^#/);
	chomp($l);
	my ($uniprotAc1, $uniprotAc2, $uniprotId1, $uniprotId2, $pdb1, $pdb2, $res1, $res2, $pfam1, $pfam2) = split(/\t/, $l);
	
	$nIntsP{$pdb1}->{$pdb2} = 1;
	$nIntsP{$pdb2}->{$pdb1} = 1;
    }
    close(F4);
    return \%nIntsP;
}

##############################################################################
sub getFusedProteins {

    my %fused = ();
    open(F5, "zcat $fusedFile |") or die "\nERROR: Failed to unzip fused file $fusedFile: $!\n\n";
    while(my $l = <F5>) {
	next if($l =~ /^#/);
	chomp($l);
	my ($uniprotIds, $pdb, $path) = split(/\t/, $l);
	$fused{$pdb}->{$uniprotIds} = 1;
    }
    close(F5);
    return \%fused;
}

##############################################################################
sub getBlastFeatures {

    my ($hIntsE, $nIntsP) = @_;

    my %bitscoreP = ();
    my %eValueP = ();
    my %bitscoreE = ();
    my %eValueE = ();
    my %n = (); # to assess how many PDBs have been found and to reduce all vs all comparison time
    my $pdbN = 0;

    open(F2, "zcat $bindingSiteFile |") or die "\nERROR: Failed to unzip binding site input file $bindingSiteFile: $!\n\n";
    my $i = 0;
    while(my $l = <F2>) {
	$i++;
	next if($l =~ /^#/);
	chomp($l);
	
	my ($ensgHit, $pdbQuery, $ident, $alignLength, $tmp1, $tmp2, $queryStart, $queryEnd, $hitStart, $hitEnd, $eValue, $bitscore)  = split(/\t/, $l);
	
	if($pdbQuery =~ /.pdb$/) {
	    $pdbQuery =~ s/(....)\..*(.)\.pdb/$1$2/;
	} else {
	    $pdbQuery =~ s/.*_.*_(.....)_.*/$1/;
	}

	next if(!defined $hIntsE->{$ensgHit} or !defined $nIntsP->{$pdbQuery});

	if(!defined $bitscoreP{$pdbQuery}){
	    $n{$pdbN++} = $pdbQuery;
	}
	if(defined $bitscoreP{$pdbQuery}->{$ensgHit}) {
	    if($bitscore > $bitscoreP{$pdbQuery}->{$ensgHit}) {
		$bitscoreP{$pdbQuery}->{$ensgHit} = $bitscore;
		$eValueP{$pdbQuery}->{$ensgHit} = $eValue;
		$bitscoreE{$ensgHit}->{$pdbQuery} = $bitscore;
		$eValueE{$ensgHit}->{$pdbQuery} = $eValue;
	    }
	} else {
	    $bitscoreP{$pdbQuery}->{$ensgHit} = $bitscore;
	    $eValueP{$pdbQuery}->{$ensgHit} = $eValue;
	    $bitscoreE{$ensgHit}->{$pdbQuery} = $bitscore;
	    $eValueE{$ensgHit}->{$pdbQuery} = $eValue;
	}
    }
    close(F2);
    return (\%eValueP, \%bitscoreP, \%eValueE, \%bitscoreE, \%n);
}
##############################################################################
sub getProbisFeatures {

    my ($hIntsPE, $nIntsP) = @_;

    my %zScoreP = ();
    my %rmsdP = ();
    my %ratioP = ();
    my %zScoreE = ();
    my %rmsdE = ();
    my %ratioE = ();
    my %n = ();
    my $pdbN = 0;
    my %pdbHit = ();

    open(F2, "zcat $bindingSiteFile |") or die "\nERROR: Failed to unzip binding site input file $bindingSiteFile: $!\n\n";
    my $i = 0;
    while(my $l = <F2>) {
	$i++;
	next if($l =~ /^#/);
	chomp($l);
	
	my ($uniprotQuery, $bsId, $pdbQuery, $pdbHit, $qBsFull, $qBs, $bs, $qBsSeq, $hBsSeq, 
	    $eValue, $rmsd, $zScore, $total, $hit, $ratio, $false) = split(/\t/, $l);

	next if(!defined $hIntsPE->{$pdbHit} or !defined $nIntsP->{$pdbQuery});

	if(!defined $zScoreP{$pdbQuery}){
	    $n{$pdbN++} = $pdbQuery;
	}

	foreach my $ensgHit (sort keys %{$hIntsPE->{$pdbHit}}) {
	    if(defined $zScoreP{$pdbQuery}->{$ensgHit}) {
		push(@{$zScoreP{$pdbQuery}->{$ensgHit}}, $zScore);
		push(@{$rmsdP{$pdbQuery}->{$ensgHit}}, $rmsd);
		push(@{$ratioP{$pdbQuery}->{$ensgHit}}, $ratio);
		push(@{$pdbHit{$pdbQuery}->{$ensgHit}}, $pdbHit);
		push(@{$zScoreE{$ensgHit}->{$pdbQuery}}, $zScore);
		push(@{$rmsdE{$ensgHit}->{$pdbQuery}}, $rmsd);
		push(@{$ratioE{$ensgHit}->{$pdbQuery}}, $ratio);
	    } else {
		push(@{$zScoreP{$pdbQuery}->{$ensgHit}}, $zScore);
		push(@{$rmsdP{$pdbQuery}->{$ensgHit}}, $rmsd);
		push(@{$ratioP{$pdbQuery}->{$ensgHit}}, $ratio);
		push(@{$pdbHit{$pdbQuery}->{$ensgHit}}, $pdbHit);
		push(@{$zScoreE{$ensgHit}->{$pdbQuery}}, $zScore);
		push(@{$rmsdE{$ensgHit}->{$pdbQuery}}, $rmsd);
		push(@{$ratioE{$ensgHit}->{$pdbQuery}}, $ratio);
	    }
	}
    }
    close(F2);
    return (\%zScoreP, \%rmsdP, \%ratioP, \%zScoreE, \%rmsdE, \%ratioE, \%n, \%pdbHit);
}
##############################################################################
sub getNeedleFeatures {

    my ($hIntsE, $nIntsP) = @_;

    my %identP = ();
    my %scoreP = ();
    my %identE = ();
    my %scoreE = ();

    my %n = ();
    my $pdbN = 0;
    my %pdbHit = ();

    open(F4, "zcat $bindingSiteFile |") or die "\nERROR: Failed to unzip binding site input file $bindingSiteFile: $!\n\n";
    my $i = 0;
    while(my $l = <F4>) {
	$i++;
	next if($l =~ /^#/);
	chomp($l);

	my ($ensgHit, $pdbQuery, $length, $ident, $sim, $gap, $identPerc, $simPerc, $gapPerc, $score, $seq1, $seq2, $alignment) = split(/\t/, $l);
	$pdbQuery =~ s/(....).*(.)\.pdb/$1$2/; 

	next if(!defined $hIntsE->{$ensgHit} or !defined $nIntsP->{$pdbQuery});

	if(!defined $identP{$pdbQuery}){
	    $n{$pdbN++} = $pdbQuery;
	}
	if(defined $identP{$pdbQuery}->{$ensgHit}) {
	    if($identPerc > $identP{$pdbQuery}->{$ensgHit}) {
		$identP{$pdbQuery}->{$ensgHit} = $identPerc;
		$scoreP{$pdbQuery}->{$ensgHit} = $score;
		$identE{$ensgHit}->{$pdbQuery} = $identPerc;
		$scoreE{$ensgHit}->{$pdbQuery} = $score;
	    }
	} else {
	    $identP{$pdbQuery}->{$ensgHit} = $identPerc;
	    $scoreP{$pdbQuery}->{$ensgHit} = $score;
	    $identE{$ensgHit}->{$pdbQuery} = $identPerc;
	    $scoreE{$ensgHit}->{$pdbQuery} = $score;
	}
    }
    close(F4);
    return (\%identP, \%scoreP, \%identE, \%scoreE, \%n);
}
##############################################################################
sub getHHsearchFeatures {

    my ($hIntsE, $nIntsP) = @_;

    my %identP = ();
    my %bitscoreP = ();
    my %probP = ();
    my %probE = ();
    my %identE = ();
    my %bitscoreE = ();

    my %n = ();
    my $pdbN = 0;
    my %pdbHit = ();

    open(F4, "zcat $bindingSiteFile |") or die "\nERROR: Failed to unzip binding site input file $bindingSiteFile: $!\n\n";
    my $i = 0;
    while(my $l = <F4>) {
	$i++;
	next if($l =~ /^#/);
	chomp($l);

	my @a = split(/\t/, $l);
	my $ensgHit = $a[1]; 
	my $pdbQuery = $a[4];
	my $prob = $a[28];
	my $ident = $a[32];
	my $bitscore = $a[38];
	my $evalue = $a[29];

	next if(!defined $hIntsE->{$ensgHit} or !defined $nIntsP->{$pdbQuery});

	if(!defined $probP{$pdbQuery}){
	    $n{$pdbN++} = $pdbQuery;
	}
	if(defined $bitscoreP{$pdbQuery}->{$ensgHit}) {
	    if($bitscore > $bitscoreP{$pdbQuery}->{$ensgHit}) {
		$probP{$pdbQuery}->{$ensgHit} = $prob;
		$identP{$pdbQuery}->{$ensgHit} = $ident;
		$bitscoreP{$pdbQuery}->{$ensgHit} = $bitscore;
		$probE{$ensgHit}->{$pdbQuery} = $prob;
		$identE{$ensgHit}->{$pdbQuery} = $ident;
		$bitscoreE{$ensgHit}->{$pdbQuery} = $bitscore;
	    }
	} else {
	    $probP{$pdbQuery}->{$ensgHit} = $prob;
	    $identP{$pdbQuery}->{$ensgHit} = $ident;
	    $bitscoreP{$pdbQuery}->{$ensgHit} = $bitscore;
	    $probE{$ensgHit}->{$pdbQuery} = $prob;
	    $identE{$ensgHit}->{$pdbQuery} = $ident;
	    $bitscoreE{$ensgHit}->{$pdbQuery} = $bitscore;
	}
    }
    close(F4);
    return (\%probP, \%identP, \%bitscoreP, \%probE, \%identE, \%bitscoreE, \%n);
}
##############################################################################
sub getBlastAndProbisFeatures {

    my ($hIntsPE, $nIntsP, $hIntsE) = @_;

    my %zScore = ();

    open(F5, "zcat $probisFile |") or die "\nERROR: Failed to unzip binding site input file $probisFile: $!\n\n";
    my $i = 0;
    while(my $l = <F5>) {
	$i++;
	next if($l =~ /^#/);
	chomp($l);
	
	my ($uniprotQuery, $bsId, $pdbQuery, $pdbHit, $qBsFull, $qBs, $bs, $qBsSeq, $hBsSeq, 
	    $eValue, $rmsd, $zScore, $total, $hit, $ratio, $false) = split(/\t/, $l);

	$zScore{$pdbHit}->{$bsId} = $zScore;
    }
    close(F5);

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

    my %bitscoreP = ();
    my %eValueP = ();
    my %zScoreP = ();
    my %bitscoreE = ();
    my %eValueE = ();
    my %zScoreE = ();
    my %bsIdP = ();
    my %bsIdE = ();
    my %n = (); # to assess how many PDBs have been found and to reduce all vs all comparison time
    my $pdbN = 0;

    open(F6, "zcat $bindingSiteFile |") or die "\nERROR: Failed to unzip binding site input file $bindingSiteFile: $!\n\n";
    $i = 0;
    while(my $l = <F6>) {
	$i++;
	next if($l =~ /^#/);
	chomp($l);
	
	my ($ensgHit, $pdbQuery, $ident, $alignLength, $tmp1, $tmp2, $queryStart, $queryEnd, $hitStart, $hitEnd, $eValue, $bitscore)  = split(/\t/, $l);
	my $bsId = "";
	if($pdbQuery =~ /^(.....)-(.*)/) {
	    $pdbQuery = $1;
	    $bsId = $2;
	} else {
	    print STDERR "ERROR: unknown format!: $pdbQuery\n";
	    exit -1;
	}
	if(!exists $zScore{$pdbQuery}->{$bsId}) {
	    print STDERR "ERROR: Not found in ProBIS file $probisFile: $pdbQuery $bsId\n";
	    exit -2;
	}
	my $zScore = $zScore{$pdbQuery}->{$bsId};

	next if(!defined $hIntsE->{$ensgHit} or !defined $nIntsP->{$pdbQuery});

	if(!defined $bitscoreP{$pdbQuery}){
	    $n{$pdbN++} = $pdbQuery;
	}
	
	if(defined $zScoreP{$pdbQuery}->{$ensgHit}) {
	    if($zScore > $zScoreP{$pdbQuery}->{$ensgHit}) {
		if($bitscore > $bitscoreP{$pdbQuery}->{$ensgHit}) {
		    $bitscoreP{$pdbQuery}->{$ensgHit} = $bitscore;
		    $eValueP{$pdbQuery}->{$ensgHit} = $eValue;
		    $zScoreP{$pdbQuery}->{$ensgHit} = $zScore;
		    $bsIdP{$pdbQuery}->{$ensgHit} = $bsId;
		    $bitscoreE{$ensgHit}->{$pdbQuery} = $bitscore;
		    $eValueE{$ensgHit}->{$pdbQuery} = $eValue;
		    $zScoreE{$ensgHit}->{$pdbQuery} = $zScore;
		    $bsIdE{$ensgHit}->{$pdbQuery} = $bsId;
		}
	    }
	} else {
	    $bitscoreP{$pdbQuery}->{$ensgHit} = $bitscore;
	    $eValueP{$pdbQuery}->{$ensgHit} = $eValue;
	    $zScoreP{$pdbQuery}->{$ensgHit} = $zScore;
	    $bsIdP{$pdbQuery}->{$ensgHit} = $bsId;
	    $bitscoreE{$ensgHit}->{$pdbQuery} = $bitscore;
	    $eValueE{$ensgHit}->{$pdbQuery} = $eValue;
	    $zScoreE{$ensgHit}->{$pdbQuery} = $zScore;
	    $bsIdE{$ensgHit}->{$pdbQuery} = $bsId;
	}
    }
    close(F6);

    return (\%eValueP, \%bitscoreP, \%zScoreP, \%bsIdP, \%eValueE, \%bitscoreE, \%zScoreE, \%bsIdE, \%n);
}
##############################################################################
sub getBlastStats {

    my ($eValue1, $eValue2, $bitscore1, $bitscore2) = @_;

    my $minBitscore = $bitscore1;
    $minBitscore = $bitscore2 if($bitscore2 < $minBitscore);
    my $maxEvalue = $eValue1;
    $maxEvalue = $eValue2 if($eValue2 > $maxEvalue);
    my $avgBitscore = ($bitscore1 + $bitscore2) / 2;
    my $avgEvalue = ($eValue1 + $eValue2) / 2;

    return ($minBitscore, $maxEvalue, $avgBitscore, $avgEvalue);
}
##############################################################################
sub getBlastProbisStats {

    my ($eValue1, $eValue2, $bitscore1, $bitscore2, $zScore1, $zScore2) = @_;

    my $minBitscore = $bitscore1;
    $minBitscore = $bitscore2 if($bitscore2 < $minBitscore);
    my $maxEvalue = $eValue1;
    $maxEvalue = $eValue2 if($eValue2 > $maxEvalue);
    my $minZscore = $zScore1;
    $minZscore = $zScore2 if($zScore2 < $minZscore);
    my $avgBitscore = ($bitscore1 + $bitscore2) / 2;
    my $avgEvalue = ($eValue1 + $eValue2) / 2;
    my $avgZscore = ($zScore1 + $zScore2) / 2;

    return ($minBitscore, $maxEvalue, $minZscore, $avgBitscore, $avgEvalue, $avgZscore);
}
##############################################################################
sub getNeedleStats {

    my ($ident1, $ident2, $score1, $score2) = @_;

    my $minIdent = $ident1;
    $minIdent = $ident2 if($ident2 < $minIdent);
    my $minScore = $score1;
    $minScore = $score2 if($score2 < $minScore);
    my $avgIdent = ($ident1 + $ident2) / 2;
    my $avgScore = ($score1 + $score2) / 2;

    return ($minScore, $minIdent, $avgScore, $avgIdent);
}
##############################################################################
sub getProbisStats {

    my ($zScore1, $zScore2, $rmsd1, $rmsd2, $ratio1, $ratio2) = @_;

    my $minZscore = $zScore1;
    $minZscore = $zScore2 if($zScore2 < $minZscore);
    my $maxRmsd = $rmsd1;
    $maxRmsd = $rmsd2 if($rmsd2 > $maxRmsd);
    my $minPvsN = $ratio1;
    $minPvsN = $ratio2 if($ratio2 < $minPvsN);
    my $avgZscore = ($zScore1 + $zScore2) / 2;
    my $avgRmsd = ($rmsd1 + $rmsd2) / 2;
    my $avgPvsN = ($ratio1 + $ratio2) / 2;

    return ($minZscore, $maxRmsd, $minPvsN, $avgZscore, $avgRmsd, $avgPvsN);
}
##############################################################################
sub getHHsearchStats {

    my ($prob1, $prob2, $ident1, $ident2, $bitscore1, $bitscore2) = @_;

    my $minProb = $prob1;
    $minProb = $prob2 if($prob2 < $minProb);
    my $minIdent = $ident1;
    $minIdent = $ident2 if($ident2 < $minIdent);
    my $minBitscore = $bitscore1;
    $minBitscore = $bitscore2 if($bitscore2 < $minBitscore);
    my $avgProb = ($prob1 + $prob2) / 2;
    my $avgIdent = ($ident1 + $ident2) / 2;
    my $avgBitscore = ($bitscore1 + $bitscore2) / 2;

    return ($minProb, $minIdent, $minBitscore, $avgProb, $avgIdent, $avgBitscore);
}

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


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

if(!defined $bindingSiteFile) {
    print STDERR "\nPlease provide a binding site input file. Try \"$0 -h\" for a complete list of options\n\n";
    exit;
}
if(!defined $hIntFile) {
    print STDERR "\nPlease provide a human interaction file. Try \"$0 -h\" for a complete list of options\n\n";
    exit;
}
if(!defined $nIntFile) {
    print STDERR "\nPlease provide a non-human interaction file. Try \"$0 -h\" for a complete list of options\n\n";
    exit;
}

print STDERR "Reading in $hIntFile ... " if(defined $verbose);
my ($hIntsE, $hIntsP, $hIntsPE, $hIntsG) = &getInteractions();
print STDERR "Done.\n" if(defined $verbose);

print STDERR "Reading in $nIntFile ... " if(defined $verbose);
my $nIntsP = &getNonInteractions();
print STDERR "Done.\n" if(defined $verbose);

my $feats = {};
my ($n, $pdbHit);
my ($identP, $identE,
    $scoreP, $scoreE,
    $eValueP, $eValueE,
    $zScoreP, $zScoreE,
    $rmsdP, $rmsdE,
    $ratioP, $ratioE,
    $probP, $probE,
    $bitscoreP, $bitscoreE,
    $bsIdP, $bsIdE);

my $fused = {};
if(defined $fusedFile) {
    print STDERR "Reading in $fusedFile ... " if(defined $verbose);
    $fused = &getFusedProteins();
    print STDERR "Done.\n" if(defined $verbose);
}
if(defined $isHHsearch) {
    print STDERR "Reading HHsearch $bindingSiteFile ... " if(defined $verbose);
    ($probP, $identP, $bitscoreP, $probE, $identE, $bitscoreE, $n) = &getHHsearchFeatures($hIntsE, $nIntsP);
    print STDERR "Done.\n" if(defined $verbose);
    $feats = $probP;
    print "#ENSG1\tENSG2\tGeneName1\tGeneName2\tPDBquery1\tPDBquery2\tMinProb\tMinIdent\tMinBitscore\tMeanProb\tMeanIdent\tMeanBitscore\ty\n";
}
elsif(defined $isBlast and defined $probisFile) {
    print STDERR "Reading ProBIS $probisFile and BLAST file $bindingSiteFile ... " if(defined $verbose);
    ($eValueP, $bitscoreP, $zScoreP, $bsIdP, $eValueE, $bitscoreE, $zScoreE, $bsIdE, $n) = &getBlastAndProbisFeatures($hIntsPE, $nIntsP, $hIntsE);
    print STDERR "Done.\n" if(defined $verbose);
    $feats = $bitscoreP;
    print "#ENSG1\tENSG2\tGeneName1\tGeneName2\tPDB1\tPDB2\tMinBitscore\tMaxEvalue\tMinZscore\tMeanBitscore\tMeanEvalue\tMeanZscore\ty\n";
} elsif(defined $isProbis) {
    print STDERR "Reading ProBIS $bindingSiteFile ... " if(defined $verbose);
    ($zScoreP, $rmsdP, $ratioP, $zScoreE, $rmsdE, $ratioE, $n, $pdbHit) = &getProbisFeatures($hIntsPE, $nIntsP);
    print STDERR "Done.\n" if(defined $verbose);
    $feats = $zScoreP;
    print "#ENSG1\tENSG2\tGeneName1\tGeneName2\tPDBhit1\tPDBhit2\tPDBquery1\tPDBquery2\tMinZscore\tMaxRmsd\tMinPvsN\tMeanZscore\tMeanRmsd\tMeanPvsN\ty\n";
}
elsif(defined $isNeedle) {
    print STDERR "Reading Needle binding sites $bindingSiteFile ... " if(defined $verbose);
    ($identP, $scoreP, $identE, $scoreE, $n) = &getNeedleFeatures($hIntsE, $nIntsP);
    print STDERR "Done.\n" if(defined $verbose);
    $feats = $identP;
    print "#ENSG1\tENSG2\tGeneName1\tGeneName2\tPDB1\tPDB2\tMinScore\tMinIdent\tMeanScore\tMeanIdent\ty\n";
} 
elsif(defined $isBlast) {
    print STDERR "Reading BLAST binding sites $bindingSiteFile ... " if(defined $verbose);
    ($eValueP, $bitscoreP, $eValueE, $bitscoreE, $n) = &getBlastFeatures($hIntsE, $nIntsP);
    print STDERR "Done.\n" if(defined $verbose);
    $feats = $bitscoreP;
    print "#ENSG1\tENSG2\tGeneName1\tGeneName2\tPDB1\tPDB2\tMinBitscore\tMaxEvalue\tMeanBitscore\tMeanEvalue\ty\n";
}


my $pdbN = (sort {$b <=> $a} keys %$n)[0];
print STDERR "Checking all vs all interactions. In total $pdbN vs $pdbN PDB structures to test ...\n" if(defined $verbose);

my %result;
my %best;
for(my $i = 0; $i <= $pdbN; $i++) {
    print STDERR int($i/$pdbN*100)."%, $i of $pdbN done.\n" if(defined $verbose and ($i % ($pdbN/10) == 0));

    for(my $j = $i + 1; $j <= $pdbN; $j++) {
	my $pdbQuery1 = $n->{$i};
	my $pdbQuery2 = $n->{$j};

#	print "1:$pdbQuery1\t$pdbQuery2\n";
	next if(!defined $nIntsP->{$pdbQuery1}->{$pdbQuery2} and !defined $withParalog);
#	print "2:$pdbQuery1\t$pdbQuery2\n";
	
	foreach my $ensgHit1 (sort keys %{$feats->{$pdbQuery1}}) {
#	    print "3:$pdbQuery1\t$pdbQuery2\t$ensgHit1\n";
	    next if(!defined $hIntsE->{$ensgHit1});

	    foreach my $ensgHit2 (sort keys %{$feats->{$pdbQuery2}}) {
#		print "4:$pdbQuery1\t$pdbQuery2\t$ensgHit1\t$ensgHit2\n";
		next if(!defined $hIntsE->{$ensgHit2});
		

		# we are not interested in homo-mers
		next if($ensgHit1 eq $ensgHit2);

		my $eH1 = $ensgHit1;
		my $eH2 = $ensgHit2;
 		my $pQ1 = $pdbQuery1;
		my $pQ2 = $pdbQuery2;
		if($eH1 ge $eH2) {
		    $eH1 = $ensgHit2;
		    $eH2 = $ensgHit1;
		    $pQ1 = $pdbQuery2;
		    $pQ2 = $pdbQuery1;
		}
		next if(defined $fusedFile and (defined $fused->{$pQ1} or defined $fused->{$pQ2}));

		if(!defined $isProbis){
		    # 1st return value is used to select best template for interaction!
		    my ($min1, $min2, $max2, $min3, $avg1, $avg2, $avg3);
		    if(defined $isNeedle) {
			($min1, $min2, $avg1, $avg2)= &getNeedleStats($identP->{$pQ1}->{$eH1},
								      $identP->{$pQ2}->{$eH2},
								      $scoreP->{$pQ1}->{$eH1},
								      $scoreP->{$pQ2}->{$eH2});
		    } elsif(defined $isBlast and defined $probisFile) {
			($min1, $max2, $min3, $avg1, $avg2, $avg3) = &getBlastProbisStats($eValueP->{$pQ1}->{$eH1},
											  $eValueP->{$pQ2}->{$eH2},
											  $bitscoreP->{$pQ1}->{$eH1},
											  $bitscoreP->{$pQ2}->{$eH2},
											  $zScoreP->{$pQ1}->{$eH1},
											  $zScoreP->{$pQ2}->{$eH2});
		    } elsif(defined $isBlast) {
			($min1, $max2, $avg1, $avg2) = &getBlastStats($eValueP->{$pQ1}->{$eH1},
								      $eValueP->{$pQ2}->{$eH2},
								      $bitscoreP->{$pQ1}->{$eH1},
								      $bitscoreP->{$pQ2}->{$eH2});
		    } elsif(defined $isHHsearch) {
			($min1, $min2, $min3, $avg1, $avg2, $avg3) = &getHHsearchStats($probP->{$pQ1}->{$eH1},
										       $probP->{$pQ2}->{$eH2},
										       $identP->{$pQ1}->{$eH1},
										       $identP->{$pQ2}->{$eH2},
										       $bitscoreP->{$pQ1}->{$eH1},
										       $bitscoreP->{$pQ2}->{$eH2}); 
		    }
		    
		    if(defined $result{$eH1}->{$eH2}) {
			# replace interaction with new template, if new template has better values. Each interaction should
			# should be described by the best template possible.
			if($min1 > $best{$eH1}->{$eH2}){
			    $best{$eH1}->{$eH2} = $min1;
			    $result{$eH1}->{$eH2} = "$eH1\t$eH2\t".
				(defined $hIntsG->{$eH1} ? $hIntsG->{$eH1}."\t" : "-\t").
				(defined $hIntsG->{$eH2} ? $hIntsG->{$eH2}."\t" : "-\t").
				"$pQ1\t$pQ2\t".
				"$min1\t".
				(defined $max2 ? "$max2\t" : "$min2\t").
				(defined $min3 ? "$min3\t" : "").
				"$avg1\t".
				"$avg2\t".
				(defined $avg3 ? "$avg3\t" : "").
				(defined $hIntsE->{$eH1}->{$eH2} ? "1\n" : "0\n");
			}
		    } else {
			$best{$eH1}->{$eH2} = $min1;
			$result{$eH1}->{$eH2} = "$eH1\t$eH2\t".
			    (defined $hIntsG->{$eH1} ? $hIntsG->{$eH1}."\t" : "-\t").
			    (defined $hIntsG->{$eH2} ? $hIntsG->{$eH2}."\t" : "-\t").
			    "$pQ1\t$pQ2\t".
			    "$min1\t".
			    (defined $max2 ? "$max2\t" : "$min2\t").
			    (defined $min3 ? "$min3\t" : "").
			    "$avg1\t".
			    "$avg2\t".
			    (defined $avg3 ? "$avg3\t" : "").
			    (defined $hIntsE->{$eH1}->{$eH2} ? "1\n" : "0\n");

		    }
		} elsif(defined $isProbis) {
		    for(my $i=0; $i < @{$pdbHit->{$pQ1}->{$eH1}}; $i++) {
			my $pH1 = $pdbHit->{$pQ1}->{$eH1}[$i];
			for(my $j=0; $j < @{$pdbHit->{$pQ2}->{$eH2}}; $j++) {
			    my $pH2 = $pdbHit->{$pQ2}->{$eH2}[$j];

# one should allow for matches to human proteins that are not in complex, to allow for conformational changes			
#			    next if(defined $hIntsE->{$eH1}->{$eH2} and !defined $hIntsP->{$pH1}->{$pH2});	
# to have a more fair comparison, skip 2nd check for fused proteins, as BLAST doesn't have this check, too.
#			    next if(defined $fusedFile and (defined $fused->{$pH1} or defined $fused->{$pH2}));

			    my ($minZscore, $maxRmsd, $minPvsN, $avgZscore, $avgRmsd, $avgPvsN) = &getProbisStats($zScoreP->{$pQ1}->{$eH1}[$i],
														  $zScoreP->{$pQ2}->{$eH2}[$j],
														  $rmsdP->{$pQ1}->{$eH1}[$i],
														  $rmsdP->{$pQ2}->{$eH2}[$j],
														  $ratioP->{$pQ1}->{$eH1}[$i],
														  $ratioP->{$pQ2}->{$eH2}[$j]); 
			    if(defined $result{$eH1}->{$eH2}) {
				# see comment above
				if($minZscore > $best{$eH1}->{$eH2}){
				    $best{$eH1}->{$eH2} = $minZscore;
				    $result{$eH1}->{$eH2} = "$eH1\t$eH2\t".
					(defined $hIntsG->{$eH1} ? $hIntsG->{$eH1}."\t" : "-\t").
					(defined $hIntsG->{$eH2} ? $hIntsG->{$eH2}."\t" : "-\t").
					"$pQ1\t$pQ2\t".
					"$pdbHit->{$pQ1}->{$eH1}[$i]\t$pdbHit->{$pQ2}->{$eH2}[$j]\t".
					"$minZscore\t".
					"$maxRmsd\t".
					"$minPvsN\t".
					"$avgZscore\t".
					"$avgRmsd\t".
					"$avgPvsN\t".
					(defined $hIntsE->{$eH1}->{$eH2} ? "1\n" : "0\n");
				}
			    } else {
				$best{$eH1}->{$eH2} = $minZscore;
				$result{$eH1}->{$eH2} = "$eH1\t$eH2\t".
				    (defined $hIntsG->{$eH1} ? $hIntsG->{$eH1}."\t" : "-\t").
				    (defined $hIntsG->{$eH2} ? $hIntsG->{$eH2}."\t" : "-\t").
				    "$pQ1\t$pQ2\t".
				    "$pdbHit->{$pQ1}->{$eH1}[$i]\t$pdbHit->{$pQ2}->{$eH2}[$j]\t".
				    "$minZscore\t".
				    "$maxRmsd\t".
				    "$minPvsN\t".
				    "$avgZscore\t".
				    "$avgRmsd\t".
				    "$avgPvsN\t".
				    (defined $hIntsE->{$eH1}->{$eH2} ? "1\n" : "0\n");
			    }
			}
		    }
		}
	    }
	}
    }
}

foreach my $eH1 (sort keys %result) {
    foreach my $eH2 (sort keys %{$result{$eH1}}) {
	print $result{$eH1}->{$eH2};
    }
}