#!/usr/bin/perl
# protmcl.perl
# wrapper for blastp > orthomcl > protdup.idtable (id, omclid, count)
# use:  in: org.aa, org.aa.blastp  ; out: org.mcl.idtable

$root="/bio/bio-grid";
$omcldir="$root/mb/orthomcl";
$bl2omcl="$root/blast/blast92orthomcl.pl";
$pval="1e-40";
$pmatch="20";

$orgid="orgdb";
$aadeflines="$orgid.deflines";

# getopt(...);

$aain= "$orgid.aa";
$blastin="$aain.blastp";
$orthodat="$aain_mcl";
$orthoout="";
$idtable="$aain_mcl.idtable";

# save in memory instead
my($ld,$naa);
$aainH={}; open($aainH,$aain);

open(AAD,">$aadeflines");
while(<$aainH>) { # in.aa
  if(/^>(\S+)/){ my $d=$1; print AAD ">$orgid:$ld size=$naa\n" if($ld); $naa=0; $ld=$d; }
  else{ chomp; $naa +=length($_); }
}
print AAD ">$orgid:$ld size=$naa\n" if($ld);
close(AAD); close($aainH);

$blastinH={}; open($blastinH,$blastin);
open(TOOMCL,"| perl $bl2omcl -org $orgid -defline $aadeflines -in stdin -out $orthodata");
while(<$blastinH>) {
  s/(IDPATT)/$orgid:$1/g;
  print TOOMCL $_;
}

system("perl -I$omcldir $omcldir/orthomcl.pl -mode 4"
 ." -pmatch_cutoff=$pmatch -pv_cutoff=$pval"
 ." -usr_bpo_file=$orthodat.bpo -usr_gg_file=$orthodat.gg");

open(OROUT,$orthoout);
open(IDTAB,">$idtable");
while(<OROUT>){
  my($oid) = m/^ORTHOMCL(\w+)/; 
  my($n)   = m/(\d+) genes/; 
  s/^.*:\s+//; 
  s/$orgid://g; 
  my @gn = split; 
  map{ $_.="\tOmcl$oid\t$n" } @gn; 
  print IDTAB join("\n",@gn),"\n";
}
close(IDTAB); close(OROUT);

> dpulex1_gnomon_mcl.idtable


__END__

melon.% more ../../orthomcl/info

gzcat $td/tanprot/dpulex1_gnomon.aa.gz | perl -ne\
'if(/^>(\S+)/){$d=$1; print ">Dappu1:$ld size=$n\n" if($ld); $n=0; $ld=$d; }\
else{ chomp; $n+=length($_); }' > dpulex1_gnomon_aa.deflines

 gzcat dpulex1_gnomon_aa.blastp.gz | perl -pi -e's/NCBI_G/Dappu1:NCBI_G/g;' \
 | perl $bg/blast/blast92orthomcl.pl -defli=dpulex1_gnomon_aa.deflines \
 -org=Dappu1 -in=stdin -orthomcl -out=dpulex1_gnomon_mcl


# add pval 1e-40 ?? default is 1e-5
perl -I$or $or/orthomcl.pl -mode 4 -usr_bpo_file=dpulex1_gnomon_mcl.bpo -usr_gg_file=dpulex1_gnomon_mcl.gg
 > & log.ormcl &

 pv_cutoff=<Float>         P-Value Cutoff used in BLAST search and/or
                           identification of inparalogs and orthologs, 1e-5
                           (DEFAULT)
   ^^ use =1e-30
 pi_cutoff=<Int>           Percent Identity Cutoff <0-100> used in identification
                           of inparalogs and orthologs, 0 (DEFAULT)
 pmatch_cutoff=<Int>       Percent Match Cutoff <0-100> used in identification of
                           inparalogs and orthologs, 0 (DEFAULT)
  ^^ use =50

cat Aug_7/all_orthomcl.out | perl -ne\
'($oid)=m/^ORTHOMCL(\w+)/; ($n)=m/(\d+) genes/; s/^.*:\s+//; s/Dappu1://g; @gn=split; \
map{$_.="\tOmcl$oid\t$n"}@gn; print join("\n",@gn),"\n";' > dpulex1_gnomon_mcl.idtable