#............ gene match stats from tblastn
# e.g. test daphnia genes x ixodes genome

set spt=ixod-dpulexprot
set spt=acyr1-dpulexgnomon
set spt=acyr1-nasvprot

# gtar -Ozxf $spt-tblastn.tgz | perl -n $bg/blast/blast9protstats.pl -p loc > $spt.blexons &
# -- or --
gunzip -c $spt.tblastn.gz | perl -n $bg/blast/blast9protstats.pl -p loc > $spt.blexons &


set ble=$spt.blexons
set blo=`echo $ble | sed -e's/blexons/blexonsort/'`
echo ====== $ble : $blo
cat $ble | perl -ne\
'($id)=m/tkey=(\S+)\-HSP:/; ($r)=m/^\w+:(\S+)/; ($tb,$te)=m/tloc=(\d+),(\d+)/;\
@v=split; $db=$v[1]; $sbin=10000*int($v[3]/10000);  \
print join("\t",$db,$r,$sbin,$v[5],$v[3],$v[4],$id,$tb,$te),"\n";'  \
| sort -k1,2 -k3,3n -k4,4nr \
> $blo

#.........

# produces .bltab7 (stats table), .blexonf (gff of all distinct gene matches)
gunzip -c $spt.tblastn.gz | perl -n $bg/blast/blast9protstats.pl -over 0.5 -best $spt.blexonsort -out $spt.bltab7

## fixup .bltab7 : col1,2 species,DB ; .blexonf annots
# perl -pi -e's/lcl\t/nasv_NCBI_GNO\t/; s/tkey=/Target=/; s/-HSP:\d+-\d+;tloc=/:/;  s/=hmm/=nasv_hmm/g;' $spt.blexonf

perl -pi -e's/\.\tHSP/daphnia_NCBI_GNO\tHSP/; s/tkey=/Target=/; s/-HSP:\d+-\d+;tloc=/:/; ' $spt.blexonf
perl -pi -e's/^SCAFFOLD\d+\t\t/acyr\tdaph_GNO\t/;  ' $spt.bltab7

# pull all _G1 best-match locs to 'genes' gff
# problem w/ multi-sameloc matches; keep only best of such genes?
# add HSP bitscores for best?

gunzip -c $spt.blexonf.gz | grep '_G1' | perl -ne\
'@v=split"\t"; ($r,$b,$e,$bt)=@v[0,3,4,5]; $v[8]=~s/;Target=.+$//;  \
 ($p)=m/Parent=(\w+)/; $p=~s/_G1//; if($p eq $lp) { $lr=$r; $le=$e; $bs+=$bt;} \
 else { if(@lv){ $lv[2]="mRNA"; $lv[4]=$le; $lv[5]=$bs; $lv[8]=~s/Parent=\w+/ID=$lp/; \
  print join("\t",@lv);}  @lv=@v; $lr=$r; $lb=$b; $lp=$p; $le=$e; $bs=$bt; }' \
 > $spt.genes


## fixme; this should allow for inside genes, not overlapping hsps
## redo for non-overlap exons then pull distinct gene ids?
cat $spt.genes | sort -k1,1 -k4,4n -k5,5nr -k6,6nr | perl -ne '@v=split"\t";\
($r,$b,$e,$s)=@v[0,3,4,5]; $p=1; if($r eq $lr and $b<$le) {$p=0;}\
print if($p); ($lr,$le,$ls)=($r,$e,$s);'\
> $spt.gene1

gunzip -c $spt.blexonf.gz | grep '_G1' | sort -k1,1 -k4,4n -k5,5nr -k6,6nr |\
perl -ne's/;Target=.+$//; @v=split; ($r,$b,$e,$bt)=@v[0,3,4,5]; \
 ($p)=m/Parent=(\w+)/; $p=~s/_G1//; if($r ne $le or $b > $le){ print "$p\n";}\
 ($lr,$lb,$le)=($r,$b,$e);' | sort | uniq > $spt.gene1.ids
 
# or this way; now includes HSP and mRNA; score sum of HSP scores
set spt=acyr1-dpulexgnomon
gunzip -c $spt.genes.gz $spt.blexonf.gz | egrep 'mRNA|_G1' | perl -pe's/_G1//;' |\
  $td/overbestgene2.perl -in stdin > $spt.gene2
#collect_gff=50501, ngene=13989
#done kept=7086, skipped=7044

set spt=acyr1-nasvprot
gunzip -c $spt.genes.gz $spt.blexonf.gz | egrep 'mRNA|_G1' | perl -pe's/_G1//;' | \
  $td/overbestgene2.perl -in stdin > $spt.gene2
#collect_gff=54428, ngene=15810
#done kept=8249, skipped=7740


# rough count of tandem dupl. genes ; can also run R stats on .bltab7
# ~/Desktop/dspp-work/eugenomestats/blast9protstats.R  : plots of mean gene stats / species

gunzip -c $spt.blexonf.gz | perl -pe ' s/_(G|S|o)\d+;/;/;' | \
 perl $td/tblastnear.perl -in stdin -over $spt.gene1 -bits 50 > & $spt.tblastnear.txt