#
## source this script with following params set:
# set dp=dana
# set species=ananassae
# set dpid=${dp}_caf051209
# set scd=$sc/${dp}3
# (source dsppblout2gff.script )>& log.bl2gff.$dpid &
#----------------
echo dspp blastout 2 gff: $dpid

  
echo "# 1. create blast index for server: $dp-chromosome"
  $bg/blast/genoblastformat.pl -dnablast -in=$scd/${dpid}.fa.gz -blid=${dpid} -out=$dp-chromosome

echo "# 2. convert blast output to gff for species x dmel genome dna"

echo '#' 2a.  ${dp3}-dmel-dna.gff 
  gtar -Ozxf $soc/${dp3}dmelc.blout.tgz |\
  $bg/blast/blast92gff.pl -swap_q2s -in=stdin -gff  -deshredquery -add_match -uniq -nodebug \
    -chrin=$sc/dmel/dmel-chromosomes.gff -qdb=${dpid} -out=${dp3}-dmel-dna.gff 

echo '#' 2a.rev  dmel-${dp3}-dna.gff 
  gtar -Ozxf $soc/${dp3}dmelc.blout.tgz |\
   $bg/blast/blast92gff.pl -noswap -in=stdin -gff -deshredquery -add_match -uniq -nodebug \
    -chrin=$sc/$dp3/${dp}-scaffolds.gff  -qdb=${dpid} -out=dmel-${dp3}-dna.gff 
 
echo '#' 2b. ${dp3}-dmel-algn.gff , ${dp3}-dmel-algn.stats 
  $bg/blast/chrgff2align.pl -stats -in=${dp3}-dmel-dna.gff -out=${dp3}-dmel-algn.gff > & ${dp3}-dmel-algn.stats 
  $bg/blast/chrgff2align.pl -stats -in=dmel-${dp3}-dna.gff -out=dmel-${dp3}-algn.gff > & dmel-${dp3}-algn.stats 
    
echo "# 3. convert blast output to gff for species genome x MOD proteomes"
 
echo '#' 3a. ${dp3}-prot9-weaker match filter
  gtar -Ozxf $soc/${dp3}prot9.blout.tgz |\
  $bg/blast/blast92gff.pl -nostringent -gff -in=stdin -out=${dp3}-prot9-weaker.gff -nodebug
  
echo '#' 3b. ${dp3}-prot9-strong2.gff stringent match filter 
  gtar -Ozxf $soc/${dp3}prot9.blout.tgz |\
  $bg/blast/blast92gff.pl -evalmin=1e-10 -stringent -gff -in=stdin \
  -out=${dp3}-prot9-strong2.gff -nodebug 

echo "# 4. make GFF for Markers, Microsatts, etc."

echo '#' 4a. dmel ${dp3}-markers.gff
  egrep 'modDM|#' ${dp3}-prot9-strong2.gff | grep -v match_part | perl -ne\
   'BEGIN{open(F,"'$sc/dmel/dmel-best-gene3.ids'") or die "dmel-best-gene3.ids";\
   %best=map{chomp;($d,$n)=split;$d,$n;} <F>;}\
   ($id)=m/ID=(\w+)/; if(/^#/){print} elsif($best{$id}){ s/tblastn:modDM/marker:modDM/; \
   s/ID=[^;]+;/ID=$best{$id};/; s/;loc=[^;\n]+//; s/[\n]/;id2=$id\n/ unless(/$id/); print;}' \
   > ${dp3}-markers.gff
   
echo '#' 4b. dros-${dp3}-micsat.gff
  gtar -Ozxf $soc/dmicsat${dp3}.blout.tgz | \
  $bg/blast/blast92gff.pl -in=stdin -gff -out=stdout -nodebug -evalmin=1e-10 | \
  perl -pe's,blastn:\w+,blastn:drosmicsat,;' > dros-${dp3}-micsat.gff

gzip *.gff