check_BDBHs.pl

#!/usr/bin/env perl

# 2013-20 Bruno Contreras-Moreira (1) and Pablo Vinuesa (2):
# 1: http://www.eead.csic.es/compbio (Laboratory of Computational Biology, EEAD/CSIC, Spain)
# 2: http://www.ccg.unam.mx/~vinuesa (Center for Genomic Sciences, UNAM, Mexico)

# This script can be used to retrieve BLAST (bidirectional) best hits produced by get_homologues.pl and get_homologues-est.pl

$|=1;

use strict;
use warnings;
use Getopt::Std;

#use File::Basename;
#use Cwd;
use FindBin '$Bin';
use lib "$Bin/lib";
use lib "$Bin/lib/bioperl-1.5.2_102/";
use marfil_homology;

my $FULLENGTHFLAG = 'flcdna'; # input sequences in files with names containing this flag are considered full length

my (%opts,$INP_dir,$INP_id,$INP_guess,$INP_maxhits,$INP_pfam,$INP_self,$INP_taxon,$INP_est);

getopts('hegsDd:i:n:t:', \%opts);

if(($opts{'h'})||(scalar(keys(%opts))==0))
{
  print   "\nusage: $0 [options]\n\n";
  print   "-h this message\n";
  print   "-d directory with get_homologues output    (example: /path/to/genomes_homologues)\n";
  print   "-i sequence identifier\n";
  print   "-s show self hits                          (optional)\n";
  print   "-n number of hits to show per species      (optional, default: -n 1)\n";
  print   "-D show Pfam protein domains if available  (optional)\n";
  print   "-g guess sequence identifier for -i label  (optional, find out ids for sequence names, such as 'rpoB')\n";
  print   "-t show only hits from this taxon          (optional, example: -t genome1.faa)\n";
  print   "-e sequences are (possibly redundant) ESTs (optional, get coverage of short isoforms if not labelled as full-length)\n\n";
  exit(0);
}

if(defined($opts{'d'})){  $INP_dir = $opts{'d'}; }
else{ die "# EXIT : need a -d directory\n"; }

if(defined($opts{'i'})){  $INP_id = $opts{'i'}; }
else{ die "# EXIT : need parameter -i\n"; }

if(defined($opts{'n'})){  $INP_maxhits = $opts{'n'}; }
else{ $INP_maxhits = 1 }

if(defined($opts{'D'})){  $INP_pfam = 1; }
else{ $INP_pfam = 0 }

if(defined($opts{'s'})){  $INP_self = 1; }
else{ $INP_self = 0 }

if(defined($opts{'t'})){  $INP_taxon = $opts{'t'}; }
else{ $INP_taxon = '' }

if(defined($opts{'g'})){  $INP_guess = 1; }

if(defined($opts{'e'})){  $INP_est = 1; }

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

my ($pQid,$querylen,$pSid,$subjectlen,$pEvalue,$ppercID,$simspan,$pbits);
my ($hit,$rest,$taxon,$idxtaxon,$fullname,@hits,%taxaids,%taxonhits,%names,%revhits);
my ($rhit,$rtaxon,$querytaxon,$cover,$use_short_sequence,$red,$red_parseOK);
my ($qcov,$scov,$input_order_file,%input_order,@sorted_taxa);

constructDirectory($INP_dir); # this sets $bpo_file and $p2ofilename

open(P2O,$p2ofilename) || die "# ERROR : cannot find $p2ofilename (generated by get_homologues), cannot proceed\n";
while(<P2O>)
{

  #8,Agos1.fasta,AgosADL393W
  ($pQid,$taxon,$fullname) = split(/,/,$_);
  if(!$taxaids{$taxon}){ $taxaids{$taxon}[0] = $pQid }
  else{ $taxaids{$taxon}[1] = $pQid }
  $names{$pQid} = $fullname;
}
close(P2O);

if($INP_guess)
{
  print "# Sequences containing label $INP_id:\n\n";
  open(P2O,$p2ofilename);
  while(<P2O>)
  {
    print if(/$INP_id/i);
  }
  close(P2O);
  exit;
}

# read input order
$input_order_file = $INP_dir."/input_order.txt";
if(-s $input_order_file)
{
  open(ORDER,$input_order_file) || warn "# ERROR : cannot read $input_order_file\n";  
  while(<ORDER>)
  {
    #0	E_coli_ST131_plasmid_pKC394.gb
    #1	E_coli_plasmid_pMUR050.gb
    #...
    if(/^(\d+)\t(\S+)/){ $input_order{$2} = $1 }
  }
  close(ORDER);
}  

# find hits of selected query >
$querytaxon = find_taxon($INP_id,\%taxaids);
open(BPO,$bpo_file) || die "# ERROR : cannot find $bpo_file (generated by get_homologues), cannot proceed\n";
flock(BPO,1);
if($taxa_bpo_index{$querytaxon}[0]){ seek(BPO,$taxa_bpo_index{$querytaxon}[0],0) }
while(<BPO>)
{
  #$pQid\t$pSid\t$pEvalue\t$ppercID\t$Qcov\t$Scov\t$pQlength\t$pSlength\t$simspan\t$pbits
  if(/^(\d+)\t(.*?)\n/)
  {
    ($pQid,$rest) = ($1,$2);

    if($pQid > $INP_id){ last }
    elsif($pQid < $INP_id){ next }
    
    ($pSid,$pEvalue,$ppercID,$qcov,$scov,$querylen,$subjectlen,$simspan,$pbits) = split("\t",$rest);
    
    next if(!$INP_self && $pSid == $pQid); # avoid selfmatches
    $taxon = find_taxon($pSid,\%taxaids);

    if($INP_est && ($querytaxon !~ /$FULLENGTHFLAG/ ||
        $taxon !~ /$FULLENGTHFLAG/)){ $use_short_sequence = 1 }
    else{ $use_short_sequence = 0 }

    next if($taxonhits{$taxon} && scalar(@{$taxonhits{$taxon}}) == $INP_maxhits);

	 if($use_short_sequence)
    {
		if($querylen < $subjectlen){ $cover = $qcov }
		else{ $cover = $scov }
    }
    else
    {
      if($querylen < $subjectlen){ $cover = $scov }
      else{ $cover = $qcov }
    }

    push(@{$taxonhits{$taxon}},$pSid);
    push(@hits,[$taxon,$pQid,$querylen,$pSid,$subjectlen,$pEvalue,$ppercID,$simspan,$pbits,$cover]);
  }
}
close(BPO);

if(!@hits)
{
  die "# ERROR: cannot find blast results for a sequence with identifier: $INP_id\n";
}
else{ construct_taxa_indexes($bpo_file); }

if($INP_pfam)
{
  if(!-s $pfam_file) # global, imported from marfil_homology
  {
    print "# WARNING: cannot retrieve Pfam results, should have been calculated previously\n";
  }
  else
  {
    my $n_of_pfam_parsed_lines = 0;
    open(PFAM,$pfam_file) || die "# cannot open $pfam_file\n";
    $n_of_pfam_parsed_lines = tr/\n/\n/ while sysread(PFAM,$_,2 ** 16);
    close(PFAM);
    printf("\n# creating Pfam indexes, this might take some time (lines=%1.2e) ...\n\n",$n_of_pfam_parsed_lines);
    construct_Pfam_hash($pfam_file,1); # %pfam_hash is now populated
  }
}

if($INP_est)
{
  if(-s $redundant_file) # global, imported from marfil_homology
  {
    printf("\n# reading redundant isoforms of last run of get_homologues-est ...\n\n");
    construct_redundant_hash($redundant_file,$INP_id); # fills global %redundant
    $red_parseOK = 1;
  }
  else
  {
    $red_parseOK = 0;
    print "# WARNING: cannot retrieve redundant isoforms, should have been calculated previously with get_homologues-est\n";
  }
}

print "\n# query = $INP_id\n";
$fullname = $names{$INP_id};
chomp($fullname);
print "# query fullname = $fullname ";

if($INP_est && $red_parseOK)
{
  $red = $redundant{$INP_id} || '';
  if($red){ print " (redundant with $red)\n"; }
  else{ print " (nr)\n"; }
}
else{ print "\n" }

print "\n# list of bidirectional best-hits:\n";
print "dir\tquery\tsbjct\tbits\tEval\t\%ident\tcover\tPfam\tannotation\n";

if(keys(%input_order)){ @sorted_taxa = sort {$input_order{$a}<=>$input_order{$b}} keys(%taxonhits) }
else{ @sorted_taxa = sort(keys(%taxonhits)) }

foreach $taxon (@sorted_taxa)
{
  print ":\t\t\t\t\t\t[$taxon]\n";
  
  next if($INP_taxon && $taxon !~ $INP_taxon);

  # print query hits >
  foreach $hit (@hits)
  {
    next if($hit->[0] ne $taxon);

    $fullname = $names{$hit->[3]};
    chomp($fullname);

    printf(">\t%s\t%s\t%d\t%g\t%1.1f\t%1.1f\t%s\t%s\n\n",
      $hit->[1],$hit->[3],$hit->[8],$hit->[5],$hit->[6],$hit->[9],
      $pfam_hash{$hit->[1]} || 'NA', $fullname);
  }

  # print reverse direction hits <
  my @sortedhits = sort {$a<=>$b} (@{$taxonhits{$taxon}});
  open(BPO,$bpo_file);
  flock(BPO,1);

  # find BPO block for this taxon
  $idxtaxon = $taxon;
  if(!defined($taxa_bpo_index{$idxtaxon})){ $idxtaxon .= '.nucl' } # EST jobs
  if(!defined($taxa_bpo_index{$idxtaxon}))
  {
    die "# ERROR: cannot find BPO index for taxon $idxtaxon\n";
  }
  seek(BPO,$taxa_bpo_index{$idxtaxon}[0],0);

  while($hit=<BPO>)
  {
    if($hit =~ /(\d+)\t(.*?)\n/)
    {
      ($pQid,$rest) = ($1,$2);
      if($pQid != $sortedhits[0])
      {
        next if(!$revhits{$sortedhits[0]});
        shift(@sortedhits);
        if(!@sortedhits){ last }
        else{ next }
      }

		($pSid,$pEvalue,$ppercID,$qcov,$scov,$querylen,$subjectlen,$simspan,$pbits) = split("\t",$rest);

      next if(!$INP_self && $pSid == $pQid); # avoid selfmatches
      $rtaxon = find_taxon($pSid,\%taxaids);
      next if($rtaxon ne $querytaxon);

      if($INP_est && ($querytaxon !~ /$FULLENGTHFLAG/ || 
        $taxon !~ /$FULLENGTHFLAG/)){ $use_short_sequence = 1 }
      else{ $use_short_sequence = 0 }

      next if($revhits{$sortedhits[0]} && scalar(@{$revhits{$sortedhits[0]}}) == $INP_maxhits);
      $cover = simspan_hsps($querylen,$subjectlen,$simspan,$use_short_sequence);
      push(@{$revhits{$sortedhits[0]}},
        [$pQid,$querylen,$pSid,$subjectlen,$pEvalue,$ppercID,$simspan,$pbits,$cover]);
    }
  }
  close(BPO);

  foreach $hit (@{$taxonhits{$taxon}})
  {
    foreach $rhit (@{$revhits{$hit}})
    {
      if($rhit->[2] eq $INP_id){ $fullname = '' }
      else{ $fullname = $names{$rhit->[2]}; chomp($fullname) }

      printf("<\t%s\t%s\t%d\t%g\t%1.1f\t%1.1f\t%s\n\n",
        $hit,$rhit->[2],$rhit->[7],$rhit->[4],$rhit->[5],$rhit->[8],
        $pfam_hash{$hit} || 'NA');
    }
  }
}

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

sub find_taxon
{
  my ($id,$ref_taxaids) = @_;
  foreach my $t (keys(%$ref_taxaids))
  {
    if($ref_taxaids->{$t}[0] <= $id && $ref_taxaids->{$t}[1] >= $id)
    {
      return $t;
    }
  }
  return 'unknown';
}