| 1 |
gpertea |
23 |
#!/bin/tcsh -f |
| 2 |
|
|
|
| 3 |
|
|
#set gmap = 'refseq_exon.gff3' |
| 4 |
|
|
set gmap = 'sim4.gff3' |
| 5 |
|
|
#set pmap = 'unipr.blat.gff3' |
| 6 |
|
|
set pmap = 'pexo_exon.gff3' |
| 7 |
|
|
|
| 8 |
|
|
set mtrack='sim4cc' |
| 9 |
|
|
set ptrack='p2g_exo' |
| 10 |
|
|
|
| 11 |
|
|
if ($2'x' == 'x') then |
| 12 |
|
|
echo "Usage: tri_jigsaw.csh <dirlist_file> <jigsaw_output_suffix>" |
| 13 |
|
|
echo " <dirlist_file> is the one corresponding to the 'both strands' run" |
| 14 |
|
|
echo " <jigsaw_output_suffix> is the one that was used when jigsaw was run" |
| 15 |
|
|
echo " (e.g. refseq.jgff) " |
| 16 |
|
|
echo ' WARNING: make sure variables $pmap and $gmap are correct!' |
| 17 |
|
|
echo " (currently they are: $pmap and $gmap)" |
| 18 |
|
|
exit 1 |
| 19 |
|
|
endif |
| 20 |
|
|
|
| 21 |
|
|
set dirlist=$1 |
| 22 |
|
|
set jsuf=$2 |
| 23 |
|
|
echo "..running:\ntri_jigsaw.csh $*" |
| 24 |
|
|
set protdb = '~/ann/protdb/unipr_mammals.fa.cidx' |
| 25 |
|
|
set predcount=0 |
| 26 |
|
|
foreach d ( `cat $dirlist` ) |
| 27 |
|
|
set bname=$d:t |
| 28 |
|
|
set dbase=$d:h |
| 29 |
|
|
#echo "dirbase=$dbase (bname=$bname)" |
| 30 |
|
|
set fname=$bname'f' |
| 31 |
|
|
set rname=$bname'r' |
| 32 |
|
|
cd $dbase |
| 33 |
|
|
|
| 34 |
|
|
set j_gff3 = $bname.recon_$jsuf.gff3 |
| 35 |
|
|
set j_anngff3 = $bname.ann.recon_$jsuf.gff3 |
| 36 |
|
|
echo "recon_jigsaw.pl -P $bname/$bname.fa -o $j_gff3 $bname/$bname.$jsuf $fname/$fname.$jsuf $rname/$rname.$jsuf" |
| 37 |
|
|
recon_jigsaw.pl -P $bname/$bname.fa -o $j_gff3 $bname/$bname.$jsuf $fname/$fname.$jsuf $rname/$rname.$jsuf |
| 38 |
|
|
#ls -al $bname/$bname.pmap.fltOK.gtf |
| 39 |
|
|
set pcount=`cat $bname/$bname.$jsuf | grep -v '^#' | cut -f9 | cut -f1 -d ';' | sort -u | wc -l` |
| 40 |
|
|
@ predcount = $predcount + $pcount |
| 41 |
|
|
set d=$bname |
| 42 |
|
|
if (-s $j_gff3) then |
| 43 |
|
|
#--validate the resulting gff |
| 44 |
|
|
# set badcds=$bname.recon_$jsuf.badCDS.gtf |
| 45 |
|
|
# gffilter -g $bname/$bname.fa -b $badcds $j_gff3 |
| 46 |
|
|
# if (-s $badcds) then |
| 47 |
|
|
# echo "**** ERROR found at CDS validation for $j_gff3! (see $badcds)" |
| 48 |
|
|
# else |
| 49 |
|
|
# /bin/rm -f $badcds |
| 50 |
|
|
# endif |
| 51 |
|
|
/bin/rm -f $d.[pg]map.i{it,fa} |
| 52 |
|
|
/bin/rm -f $d.pmap.iit $d.pmap.ifa |
| 53 |
|
|
#gff2iit -o $d.pmap -t $ptrack $d/$d.$pmap |
| 54 |
|
|
gtf2gff -t $ptrack $d/$d.$pmap > $d.iit_pmap.gff3 |
| 55 |
|
|
gtf2gff -t 'pmap' $d/$d.pmap_exon.gff3 >> $d.iit_pmap.gff3 |
| 56 |
|
|
gff2iit -o $d.pmap $d.iit_pmap.gff3 |
| 57 |
|
|
/bin/rm -f $d.gmap.iit $d.gmap.ifa |
| 58 |
|
|
#gff2iit -o $d.gmap -t $mtrack $d/$d.$gmap |
| 59 |
|
|
gtf2gff -t $mtrack $d/$d.$gmap | perl -pe 's/\.mrna(\d+)/\.m$1/' > $d.iit_gmap.gff3 |
| 60 |
|
|
gtf2gff -t 'gmap' $d/$d.gmap_exon.gff3 >> $d.iit_gmap.gff3 |
| 61 |
|
|
gff2iit -o $d.gmap -t $mtrack $d.iit_gmap.gff3 |
| 62 |
|
|
|
| 63 |
|
|
#echo "gffann.pl -o $j_anngff3 -t $d -m $d.gmap.iit -p $d.pmap.iit -P $protdb $j_gff3" |
| 64 |
|
|
gffann.pl -o $j_anngff3 -t $d -m $d.gmap.iit -p $d.pmap.iit -P $protdb $j_gff3 |
| 65 |
|
|
# -- temporary: throw in all the mapping evidence |
| 66 |
|
|
gtf2gff -t $mtrack $d/$d.$gmap > $d.top5gmap.gff3 |
| 67 |
|
|
gtf2gff -t $ptrack $d/$d.$pmap > $d.top5pmap.gff3 |
| 68 |
|
|
gtf2gff -t 'gmap' $d/$d.gmap_exon.gff3 >> $d.top5gmap.gff3 |
| 69 |
|
|
|
| 70 |
|
|
gtf2gff -t 'pmap' $d/$d.pmap_exon.gff3 >> $d.top5pmap.gff3 |
| 71 |
|
|
gtf2gff -t 'blat' $d/$d.pblat.gff3 >> $d.top5pmap.gff3 |
| 72 |
|
|
|
| 73 |
|
|
else #no predictions.. just show whatever evidence is there |
| 74 |
|
|
cp $j_gff3 $j_anngff3 |
| 75 |
|
|
gtf2gff -t $mtrack $d/$d.$gmap | perl -pe 's/\.mrna(\d+)/\.m$1/' > $d.top5gmap.gff3 |
| 76 |
|
|
gtf2gff -t $ptrack $d/$d.$pmap > $d.top5pmap.gff3 |
| 77 |
|
|
gtf2gff -t 'gmap' $d/$d.gmap_exon.gff3 >> $d.top5gmap.gff3 |
| 78 |
|
|
gtf2gff -t 'pmap' $d/$d.pmap_exon.gff3 >> $d.top5pmap.gff3 |
| 79 |
|
|
gtf2gff -t 'blat' $d/$d.pblat.gff3 >> $d.top5pmap.gff3 |
| 80 |
|
|
endif |
| 81 |
|
|
|
| 82 |
|
|
cd .. |
| 83 |
|
|
end |
| 84 |
|
|
|
| 85 |
|
|
echo "Total: $predcount gene models predicted." |
