#!/bin/tcsh
# megablast find predicted exons in genomes; run per genome dir

set nb=/bio/bio-grid/ncbi-20060507/ncbi/bin

# dmel case:
#set scs=( 4 2L 2R 3L 3R X )
# dpse case:
# set scs=( 2 3 4* 5* X* )
# dyak case:
# set scs=( chr* )
# cele case:
# set scs=( I II III IV V X )
# dmoj fix
set scs=(scaffold_6308)
# most cases:
#set scs=(scaffold_*)

#for filetest w/ globs
set nonomatch=1
foreach dir ($scs)

  set fex=$dir/*_exons.fa.gz
  echo "# split $fex by predictors"
  gunzip -c $fex | env fin=$fex perl -ne 'use FileHandle;\
    if(/^>(\D+)/){ $gr=$1; $gr=~s/_$//; $fh=undef;\
    $skip= ($gr=~m/DGIL_SNO|BATZ_CNA|RGUI_GID|^TRd\w+/); next if($skip);\
    $fn=$ENV{fin}; $fn=~s,/[^/]*$,/$gr-pexons.fa, or die;  \
    $fh=$fhs{$fn}; unless($fh){ $fhs{$fn}= $fh= new FileHandle "> $fn"; } }\
    print $fh $_ unless($skip);'

 # fixme: if( ! test $dir/*-pexons.fa ) continue
 set pfe=`filetest -s $dir/*-pexons.fa` ; if ( $pfe[1] != 1 ) continue

 # fixme _caf > other, but not all *.fa.gz > dsec_br051028.fa.gz
 # set faz="$dir/d*_caf060210.fa.gz"
 set faz=`/bin/ls -1 $dir/*.fa.gz | grep -v exons`
 set fna=`echo $faz | sed -e's/.gz//'`
 gzcat $faz | $nb/formatdb -pF -n $fna  -i stdin

 echo "# megablast start" $fna , `date`
 foreach predict ($dir/*-pexons.fa)
  echo "# megablast query=$predict"
  $nb/megablast -W11 -t16 -N0 -gT -e0.05 -D3 -i $predict -d $fna -o $predict.mbout
 end
 
 /bin/rm $fna.{nsq,nin,nhr}
 echo "# megablast done" $fna , `date`
end

exit
# --- formatdb ------ now in main loop
# --- exon split ---- now in main loop