#!/bin/tcsh
# intronannot.sh

set td=/bio/bio-grid/dpulex/prots/

# ## NOTE: want to use JGI v11 genes; EVM daphe has JGI FM5 - Gno - Sno set
set preds=(JGI NCBI_GNO DGIL_SNO)
set dp=dpulex1

foreach pred ($preds)
  set predexon=exon
  if( $pred == "DGIL_SNO" ) set predexon=CDS
 
  echo "# pred:$pred species:$dp :: intron annotation"

  gzcat *_predict.gff.gz | grep "$pred.$predexon" | \
    $td/overlapfilter.perl -act mark -mark=cdnapasa -typeover overlap -pctover 50 \
    -in stdin  -over *_estmatch.gff.gz \
   | $td/overlapfilter.perl -act mark -mark=modprot -typeover overlap -pctover 50 \
    -in stdin  -over *prothsp.gff.gz \
   | $td/overlapfilter.perl -act mark -mark=termexon -typeover overlap -pctover 50 \
    -in stdin  -over *terminal_exons.gff.gz \
   | env exon=$predexon $td/exon2intron.pl > ${dp}_introns.$pred.gff
  
  #? want to keep exon annot as file set?
    
  egrep 'intergenic_region|#' ${dp}_introns.$pred.gff  > ${dp}_intergene.$pred.gff
  egrep -v 'intergenic_region' ${dp}_introns.$pred.gff  > ${dp}_intron1.$pred.gff
  /bin/mv  ${dp}_intron1.$pred.gff  ${dp}_introns.$pred.gff
  
  cat ${dp}_introns.$pred.gff | \
    $td/overlapfilter.perl -act mark -mark=inest -typeover overlap -pctover 50 \
    -in stdin  -over *_estmatch.gff.gz \
   | $td/overlapfilter.perl -act mark -mark=inprot -typeover overlap -pctover 50 \
    -in stdin  -over *prothsp.gff.gz \
      > ${dp}_introns2.$pred.gff
  /bin/mv  ${dp}_introns2.$pred.gff  ${dp}_introns.$pred.gff
  
  cat ${dp}_intergene.$pred.gff | \
    $td/overlapfilter.perl -act mark -mark=inest -typeover overlap -pctover 50 \
    -in stdin  -over *_estmatch.gff.gz \
   | $td/overlapfilter.perl -act mark -mark=inprot -typeover overlap -pctover 50 \
    -in stdin  -over *prothsp.gff.gz \
      > ${dp}_intergene2.$pred.gff
  /bin/mv  ${dp}_intergene2.$pred.gff  ${dp}_intergene.$pred.gff

end

# also find uniq subsets of predictors: JGI, GNO-JGI, SNO-JGI-GNO
cat ${dp}_introns.JGI.gff |\
$td/overlapfilter.perl -act drop -typeover sameloc  \
  -in ${dp}_introns.NCBI_GNO.gff  -over stdin \
    > ${dp}_introns3.NCBI_GNO.gff

cat ${dp}_introns.JGI.gff ${dp}_introns.NCBI_GNO.gff |\
$td/overlapfilter.perl -act drop -typeover sameloc  \
  -in ${dp}_introns.DGIL_SNO.gff  -over stdin \
    > ${dp}_introns3.DGIL_SNO.gff

cat ${dp}_intergene.JGI.gff |\
$td/overlapfilter.perl -act drop -typeover sameloc  \
  -in ${dp}_intergene.NCBI_GNO.gff  -over stdin \
    > ${dp}_intergene3.NCBI_GNO.gff

cat ${dp}_intergene.JGI.gff ${dp}_intergene.NCBI_GNO.gff |\
$td/overlapfilter.perl -act drop -typeover sameloc  \
  -in ${dp}_intergene.DGIL_SNO.gff  -over stdin \
    > ${dp}_intergene3.DGIL_SNO.gff

exit;

set info=<<EOT

# want
1) introns for both gnomon and JGIv1.1 gene predictions
  -- add  SNAP not in above
2) exons for the gene models that were validated by your PASA predictions
3) introns for the gene models that were validated by PASA
4) intergenic spacers for both gnomon and JGIv1.1 predictions

steps:
  annot exons from predictions with: pasa-asm, est? (e.g. pasa rejected), prot-hsp,
     with overlapfilter -pctover=50 (75?, 90?)
    -- use overlap flags not ids for annots? e.g. haspasa=1; hasest=1; ...
    -- distinguish CDS and UTR exons; split exon to cds/utr parts?
    -- keep gene ID w/ exons    

  pull introns, attach annots of paired exons 
  pull intergenes
  annotate intron/intergene gff with pasa-asm, est, prothsp, ? tilex, other?


daphe/scaffold_4/*.gff.gz

  33965 Jun 21 12:37 daphe/scaffold_4/dpulex1_estmatch.gff.gz
        == alignAssembly-pasa_daphc        cDNA_match 
 238764 Jun 21 12:36 daphe/scaffold_4/dpulex1_predict.gff.gz
 121447 Jun 21 12:36 daphe/scaffold_4/dpulex1_prothsp.gff.gz
   7307 Jun 21 12:37 daphe/scaffold_4/dpulex1_terminal_exons.gff.gz

EOT