#!/usr/local/bin/perl
# anoblast.pl
# jun02 q'n'd pep blast of anopheles (Ensembl peptides)

=head1 NOTES
# use this sort on hgtable to get top mosq hits for each FBgn
  
grep '^#' hgtable > ! shgtab ; grep -v '^#' hgtable | sort -k2,2 -k5,5nr | uniq >> shgtab

=cut

use Getopt::Long;    
use flybase::sym2id;

my $workpath= '/c7/eugenes/genomew/';  
$makequery= 0;
$doblast= 0;
$dohgtable= 0; # see longer blastout.pl parsing non-table blast output

@sys= `uname -X`;
%sys= map { chomp; split(' = '); } @sys;
if ($sys{Node} eq 'oat') { $TMPD= "/c7/eugenes/tmp/" ; }
# elsif ($sys{Node} eq 'kalo')  { $TMPD= "/c5/tmp"; } #??

$refsq='/c7/eugenes/refseq1/';
$droaadb= $refsq.'aagadfly2.fa'; # 9191669b, bdgpseqs/aa_gadfly.dros.RELEASE2
$anoaadb= $refsq.'Anopheles_gambiae.pep.fa'; # 8644017b
$blastout= $refsq.'anogad.blastout';
$flyfeats= "/c7/eugenes/genomes/fly/features*.tsv";
$flyidmap= "/c7/eugenes/genomes/fly/idmap.tsv";
$flyacode= "/bio/meow-pub/eugenes/server/fly/acode";
$savefbgn= $refsq."aagadfly2.missingfbgn";
$mosqidmap= $workpath."mosquito/idmap.tsv";

# > entries
# Anopheles_gambiae.pep.fa:15101
# aagadfly2.fa:14335

$blast= "/bio/mb/blast2/blastall";
$formatano="/bio/mb/blast2/formatdb  -p T -t 'Anopheles Emsembl peptides' -i $anoaadb ";
$formatgad="/bio/mb/blast2/formatdb  -p T -t 'Drosophila GadFly2 peptides' -i $droaadb ";

  # -a #cpus  -v 100 limit #matches below 500 default?
$blastopts = "-p blastp -v 5 -b 5 -F F  -e 1e-30 -m 9 -a $sys{NumCPU}";
$blastvers= 'BLASTP 2.2.2 [Jan-08-2002]';  # get from blast output !

$qorg= 'fly';  # query
$sorg= 'mosquito'; # subject

@blastcmd= ("$blast $opts -d $anoaadb -i $droaadb -o $blastout");

GetOptions( 
'workpath=s' => \$workpath,
'run!' => \$dorun,
'blast!' => \$doblast,
'hgtable!' => \$dohgtable,
'cleantable!' => \$cleantable,
'syntable!' => \$dosyntab,
'debug!' => \$debug,
);
#'dnafile=s' => \$dnafile,
#'subrange=s' => \$subrange,
#'org=s' => \$org,
#'csomes=s' => \@csome,

$usage= <<USAGE;
$0 -- blast anopheles against fly seqs (proteins)
options:
'workpath=s' => $workpath,
'blast!' => $doblast,
'hgtable!' => $dohgtable,
'cleantable!' => $cleantable,
'syntable!' => $dosyntab,
'debug!' => $debug,
 blastcmd => @blastcmd,
USAGE


die $usage unless($dosyntab || $doblast || $dohgtable || $cleantable);

@blastcmd= ("$blast $opts -d $anoaadb -i $droaadb -o $blastout");

$cleantable= $cleantable || $dohgtable;

if ($doblast) { system(@blastcmd); }

if ($dohgtable) { makeHgTable(); }

if ($cleantable) {	cleanHtTable(); }

if ($dosyntab) { syntenyTab(); }


exit;
#------


sub makeHgTable {
	my $ofile= "$workpath/$sorg/hgtable";
	# save old hgtable
	if (-f $ofile) { rename($ofile,"$ofile.old"); }
	open( O,">>$ofile") || die "cannot write $ofile";
	print O &hgtabheader("$sorg/hgtable",$org,$fdate, $blstats) unless($needstats);
	
	initSymbols('fly',$flyacode);
	getMosP_G_Ids($anoaadb); # $pg{$p}
	getMissingFBgn($droaadb); # cgfb
	
	$sym= '--';
	open(IN,$blastout) || die $blastout;
	while(<IN>) {
		chomp;
		if (/^#/) {
			# parse some of this
			}
		else {
		  @v= split "\t";
			my ($query, $sid, $pctident, $alignment_length, $mismatches, $gap_openings, 
				 $q_start, $q_end, $s_start, $s_end, $prob, $bit_score ) = @v;
			
			$query =~ m/(FBgn\d+)/; $qid= $1;
			$query =~ m/(FBan\d+)/; $fban= $1;
			$refsq=  ($fban) ? "GadFly:$1" :  '';
			unless($qid) {
				$qid= $cgfb{$fban};
				}
				
			# $dbid= "FlyBase:$qid"; ##??
			# $sym= id2sym($qid);
			$pctident= int($pctident);
			$bit_score= int($bit_score);
			$sid= $pg{$sid} || $sid;

			$sym= flybase::sym2id::id2symbol($qid) || '--';
			
#			$oline= join("\t", $sid, $qid, $qorg, $sym, 
#				$bit_score, $prob, $pctident, $refsq, $dbid );
			$oline= "$sid\t$qid\t$qorg\t$sym\t$bit_score\t$prob\t$pctident\t$refsq\t$dbid\n";
			print O $oline; 
			}
		}
	close(O); close(IN);

}

sub cleanHtTable {
 #? add fly gene location to hgclean ?
 my $hgtable= "$workpath/$sorg/hgtable";
 my $shgtab= "$workpath/$sorg/hgclean";
 unlink($shgtab);
 my $cmd= "grep '^#' $hgtable > $shgtab ; grep -v '^#' $hgtable | sort -k2,2 -k5,5nr | uniq >> $shgtab";
 warn $cmd if $debug;
 system($cmd);
 
 my $h= readIdMap($flyidmap);
 %flyidmap= %$h;
 $h= readIdMap($mosqidmap);
 %mosqidmap= %$h;
 
 open(F,$shgtab) || die $shgtab;
 open(O,">$shgtab.1") || die "$shgtab.1";
 while(<F>){
 	if (/^\w/) {
 		my @v= split "\t";
 		my $mid= $v[0];
 		my $fid= $v[1];
 		my $bits= $v[4];
 		next if ($fid eq $lastid && $bits < $minbits);
 		chomp; print O $_;
 		
 		##? add these locations?
 		my $floc;
 		$floc = $mosqidmap{$mid} || [ '--', '--' ];
 		if (ref $floc) { print O join("\t",@ $floc),"\t"; }
 		$floc= $flyidmap{$fid} || [ '--', '--' ];
 		if (ref $floc) { print O join("\t",@ $floc); }
 		
 		print O "\n";
 		if ($fid ne $lastid) {
 			$lastid= $fid; 
 			$minbits= $bits/2;
 			}
 		}
 	else { print O $_; }
 	}
 	close(O); close(F);
 	rename("$shgtab.1","$shgtab");
}
# also grep -v '^#' hgclean | sort -k9,9 -k10,10n -k2,2 -k5,5nr > hgloc



sub syntenyTab {
	my $shgtab= "$workpath/$sorg/hgloc"; # loc sorted hgclean
	my $syntab= "$workpath/$sorg/synt-fly.tab";
	open(F,$shgtab) || die $shgtab;
	open(O,">$syntab") || die "$syntab";
 	print O join("\t",
 		qw(dID	dChr	dBstart	dBend	dPct	wID	wChr	wBstart	wBend	wPct)),"\n";
	while(<F>){
	 	if (/^\w/) {
	 		chomp;
 			my @v= split "\t";
			my($did, $wid, $qorg, $sym, 
 				$bit_score, $prob, $pct, $refsq,  
 				$dChr, $dLoc, $wChr, $wLoc )= @v;
 			next unless($dChr ne '--' && $wChr ne '--');
 			$dLoc =~ m/(\d+)\.\.(\d+)/; my($dStart,$dEnd)= ($1,$2);
 			$wLoc =~ m/(\d+)\.\.(\d+)/; my($wStart,$wEnd)= ($1,$2);
			print O join("\t",
				$did,$dChr,$dStart,$dEnd,$pct,
				$wid,$wChr,$wStart,$wEnd,$pct),"\n";
			}
		}
	close(O); close(F);
}


sub bestLocmatch {
 my $shgtab= "$workpath/$sorg/hgclean";
 
# my $h= readIdMap($flyidmap);
# %flyidmap= %$h;
# $h= readIdMap($mosqidmap);
# %mosqidmap= %$h;
 
 open(F,$shgtab) || die $shgtab;
 open(O,">$shgtab.best") || die "$shgtab.best";
 while(<F>){
 	if (/^\w/) {
 		my @v= split "\t";
 		my $mid= $v[0];
 		my $fid= $v[1];
 		my $bits= $v[4];
 		my $pct= $v[4];

 	#	next if ($fid eq $lastid && $bits < $minbits);
 	#	chomp; print O $_;
# 		print O "\n";
# 		if ($fid ne $lastid) {
# 			$lastid= $fid; 
# 			$minbits= $bits/2;
# 			}
 		}
# 	else { print O $_; }
 	}
 	close(O); 
 	close(F);
}

=head2 bestloc note

test region:

mosq. 2L 3348009..3829138
fly   3L 3376287   8653140..9622936   5544194..6203106
find runs of >pctident, m-chr,loc & f-chr,loc same  -- sorted by m-chr,loc
m-id f-id org f-sym bits eval pctident protid m-chr m-loc f-chr f-loc

AGgn0008188     FBgn0025630     fly     EG:22E5.3       194     2.9e-50 37      GadF
ly:FBan0004061  2L      3348009..3349170        X       1820937..1822277
AGgn0002690     FBgn0014388     fly     sty     293     1.1e-79 35      GadFly:FBan0
001921  2L      3392466..3394601        3L      3376287..3378408
AGgn0008279     FBgn0026259     fly     cIF2    643     0.0     58      GadFly:FBan0
010840  2L      3430391..3457978        3L      3430557..3434058
AGgn0008367     FBgn0027500     fly     BcDNA:LD24702   263     2.5e-70 33      GadF
ly:FBan0017286  2L      3471096..3473154        3L      16489812..16493649
AGgn0000975     FBgn0016081     fly     fry     1780    0.0     67      GadFly:FBan0
006774  2L      3540816..3578673        3L      9551330..9557158
AGgn0000975     FBgn0016081     fly     fry     1748    0.0     63      GadFly:FBan0
006780  2L      3540816..3578673        3L      9551330..9557158
AGgn0010231     FBgn0036028     fly     CG16717 490     3e-139  73      GadFly:FBan0
016717  2L      3593330..3594307        3L      9599249..9600151
AGgn0010668     FBgn0036030     fly     CG6767  621     1e-178  93      GadFly:FBan0
006767  2L      3596581..3605307        3L      9603707..9616005
AGgn0010668     FBgn0036030     fly     CG6767  603     2e-173  91      GadFly:FBan0
006767  2L      3596581..3605307        3L      9603707..9616005
AGgn0010577     FBgn0015321     fly     UbcD4   263     3.6e-71 67      GadFly:FBan0
008284  2L      3604886..3606647        3L      9617683..9619243
AGgn0010732     FBgn0036031     fly     CG6761  591     3e-169  51      GadFly:FBan0
006761  2L      3607503..3609206        3L      9620042..9622936
AGgn0010716     FBgn0005533     fly     RpS17   206     2.2e-54 82      GadFly:FBan0
003922  2L      3610849..3611768        3L      9347928..9348964
AGgn0010662     FBgn0027066     fly     Eb1     339     9.2e-94 61      GadFly:FBan0
003265  2L      3620637..3625516        2R      1792201..1799006
AGgn0010662     FBgn0034403     fly     CG18190 155     1.0e-38 60      GadFly:FBan0
018190  2L      3620637..3625516        2R      13802384..13803175
AGgn0010772     FBgn0033100     fly     CG3420  149     3.5e-37 63      GadFly:FBan0
003420  2L      3627654..3628349        2R      1799257..1799928
AGgn0000990     FBgn0033757     fly     CG8811  867     0.0     53      GadFly:FBan0
008811  2L      3628838..3633025        2R      7629934..7633271
AGgn0000991     FBgn0035622     fly     CG10590 916     0.0     80      GadFly:FBan0
010590  2L      3637553..3646472        3L      5544194..5546362
AGgn0010176     FBgn0035696     fly     Best2   568     2e-162  79      GadFly:FBan0
010173  2L      3650766..3655896        3L      6187069..6192573
AGgn0010176     FBgn0036491     fly     Best4   331     2.2e-91 49      GadFly:FBan0
007259  2L      3650766..3655896        3L      15197994..15199585
AGgn0010176     FBgn0036492     fly     Best3   397     5e-111  47      GadFly:FBan0
012327  2L      3650766..3655896        3L      15199649..15201433
AGgn0010176     FBgn0040238     fly     Best1   436     2e-122  51      GadFly:FBan0
006264  2L      3650766..3655896        3R      5957792..5966551
AGgn0010578     FBgn0001258     fly     ImpL3   510     4e-145  77      GadFly:FBan0
010160  2L      3673465..3676506        3L      6199921..6203106
AGgn0010578     FBgn0033856     fly     CG13334 248     2.1e-66 39      GadFly:FBan0
013334  2L      3673465..3676506        2R      8543853..8544908
AGgn0000994     FBgn0042105     fly     CG18748 132     2.0e-31 33      GadFly:FBan0
018748  2L      3701601..3702299        3R      3598536..3599735
AGgn0010674     FBgn0003149     fly     Prm     1219    0.0     71      GadFly:FBan0
005939  2L      3716649..3729920        3L      8653140..8665597
AGgn0000999     FBgn0022702     fly     Cht2    282     1.5e-76 41      GadFly:FBan0
002054  2L      3734884..3736011        3L      1730959..1734036
AGgn0000999     FBgn0034580     fly     CG9357  340     4.4e-94 45      GadFly:FBan0
009357  2L      3734884..3736011        2R      15984068..15985619
AGgn0000999     FBgn0034582     fly     CG10531 306     6.9e-84 42      GadFly:FBan0
010531  2L      3734884..3736011        2R      15990454..15991694
AGgn0001001     FBgn0035925     fly     CG5797  815     0.0     75      GadFly:FBan0
005797  2L      3781624..3810795        3L      8697258..8711721
AGgn0001001     FBgn0035925     fly     CG5797  623     1e-178  62      GadFly:FBan0
005797  2L      3781624..3810795        3L      8697258..8711721
AGgn0010581     FBgn0035928     fly     CG13310 524     2e-149  66      GadFly:FBan0
013310  2L      3819805..3821023        3L      8719139..8720767
AGgn0010581     FBgn0038695     fly     CG14280 149     2.3e-36 29      GadFly:FBan0
014280  2L      3819805..3821023        3R      15018468..15020644
AGgn0001009     FBgn0030992     fly     CG7633  300     4.5e-82 71      GadFly:FBan0
007633  2L      3825202..3829138        X       18923962..18925532
AGgn0001009     FBgn0030993     fly     CG7635  325     1.8e-89 68      GadFly:FBan0

=cut


sub readIdMap {
	my $idmap= shift;
	my %idmap= ();
	open(ID,$idmap) || die $idmap;
	while(<ID>) {
		chomp; my($id,$csome,$bloc)= split "\t";
		$idmap{$id}= [$csome,$bloc];
		}
	close(ID);
	return \%idmap;
}

sub getMosP_G_Ids {
	my $x= shift;
 	open(F,"grep '^>' $x|") || die $x;
 	while(<F>) {
 	 	 m/(ENSANGP0+\d+)\s+Gene:ENSANGG0+(\d+)/;
 	 	 $p=$1; $g=$2;
 	 	 $pg{$p}= 'AGgn'.sprintf("%07d",$g);
 		}	
 	close(F);
}

sub getMissingFBgn {
	my $x= shift;
	if (-r $savefbgn) {
		open(F,"$savefbgn");
		while(<F>){ chomp; ($c,$f)= split "\t"; $cgfb{$c}= $f; }
		close(F); return;
	}

  $cgfb{'FBan0017696'}= 'FBgn0010247'; ## missing from feats
  $cgfb{'FBan0017718'}= 'FBgn0010247'; ## missing from feats - has 3 fbans FBan0017685
 	open(F,"grep '^>' $x | grep -v FBgn|") || die $x;
 	while(<F>) {
		/(CG\d+)/; $c= $1;
		/(FBan\d+)/; $a= $1;
		$v= `fgrep $c $flyfeats`;
		if ($v =~ /(FBgn\d+)/) { $cgfb{$a}= $1; }
		}
	close(F);
	if ($savefbgn) {
		open(F,">$savefbgn");
		foreach $cg (keys %cgfb) { print F "$cg\t$cgfb{$cg}\n"; }
		close(F);
	  }
}			

sub initSymbols {
	my( $org, $orgacode)= @_;
	# $myorg= $org if ($org);
	flybase::sym2id::set4meow($isMeowId);
	# flybase::sym2id::setcaseless(1) if ($self->{org} =~ m/weed/);
	flybase::sym2id::readSym2Id($orgacode, $org) if (-f $orgacode);
}
	
sub hgtabheader {
  my ($fname,$org,$date,$blstats)= @_;
	my $o= <<TEOF;
\# $fname
\# Homologous genes for $org computed $date with $blastvers 
TEOF

	$o .= "\# BLAST command: ".join("\# ",@blastcmd) if (scalar(@blastcmd));   	
	$o .= "\n";
	$o .= $blstats if ($blstats);

	$o .= <<TEOF;
\# ID lookup at http://iubio.bio.indiana.edu/meow/.bin/moidq.html?id
\# Includes source-id == homolog-id to document source
\#
\# tab-separated-values
\# source  	homolog   	organ-	gene	blast	blast	pct 	homolog	 ref. db
\#  id     	  id      	 ism  	sym.	score	prob 	ident	refseq	    id

TEOF
	return $o;
}

=head1 Sample blast output

# BLASTP 2.2.2 [Jan-08-2002]
# Database: /c7/eugenes/refseq1/Anopheles_gambiae.pep.fa
# Query: puromycin-sensitive-aminopeptidase|FBgn0035226|pp-CT1016|FBan0001009 GO
:[] mol_weight=99328  located on: 3L 62A5-62A5; 
# Fields: Query id, Subject id, % identity, alignment length, mismatches, gap op
enings, q. start, q. end, s. start, s. end, e-value, bit score
puromycin-sensitive-aminopeptidase|FBgn0035226|pp-CT1016|FBan0001009    ENSANGP0
0000004374      65.33   845     255     5       9       851     1       809     
0.0     1097.8
puromycin-sensitive-aminopeptidase|FBgn0035226|pp-CT1016|FBan0001009    ENSANGP0
0000020286      34.97   895     526     18      9       864     108     985     
8e-138  487.3

FIXME:
	$blstats = "\# No. sequences=$nseqs; No. letters=$nletters; Expectation cutoff=$blastmaxprob \n";
	$blstats .= "\# Matrix=$matrix; Gap Existence=$gapo; Gap Extension=$gape \n";

check mosq. ensembl ENSANGG id range

grep '^gene' features-*.tsv |\
 perl -e '$b=99999;while(<>){ $d=$1 if(/ENSANGG0*(\d+)/);push(@d,$d); \
 $b=$d if($b>$d);$e=$d if($e<$d);} print "ENSANGG b=$b ; e=$e; n=".scalar(@d)."\n";'

2L only>> ENSANGG b=00000000994 ; e=00000019785; n=2821
all    >> ENSANGG b=00000000207 ; e=00000019857; n=12704
 -- safe to use ENSANGG as eugenes id#  AGgn
 
 -- check P vs G ids - not same.
 grep '^>' $e/refseq1/Anopheles_gambiae.pep.fa |\
 perl -e 'while(<>){ m/ENSANGP0+(\d+) Gene:ENSANGG0+(\d+)/;$p=$1;$g=$2; \
  if ($g == $p){$n++;}else{$x++;}} print "P=G id same=$n diff=$x\n";'
P=G id same=47 diff=15054
 
 -- missing FBgn s
 grep '^>' $e/refseq1/aagadfly2.fa | grep -v FBgn |wc
 == 99 - CG, FBan, no GO or FBgn
 
 $droaadb= $refsq.'aagadfly2.fa'; # 9191669b, bdgpseqs/aa_gadfly.dros.RELEASE2

grep '^#' hgtable > ! shgtab ; grep -v '^#' hgtable | sort -k2,2 -k5,5nr | uniq >> shgtab


== sort hgclean by fly loc:
 my $cmd= "grep '^#' $hgtable > $shgtab ; 
 
 # by mosq loc
 grep -v '^#' hgclean | sort -k9,9 -k10,10n -k2,2 -k5,5nr > hgloc
 # by fly loc
 grep -v '^#' hgclean | sort -k11,11 -k12,12n -k2,2 -k5,5nr > hgloc

AGgn0018733     FBgn0000008     fly     a       134     1.9e-31 73      GadFly:FBan0006741      2R      17087053..17089167
AGgn0005289     FBgn0000011     fly     ab      478     3e-135  47      GadFly:FBan0004807      2L      11133135..11158800
 


=cut