| 1 |
gpertea |
23 |
#!/usr/bin/perl |
| 2 |
|
|
use strict; |
| 3 |
|
|
use Getopt::Std; |
| 4 |
|
|
use FindBin;use lib $FindBin::Bin; |
| 5 |
|
|
umask 0002; |
| 6 |
|
|
my $usage=q/ |
| 7 |
|
|
Usage: |
| 8 |
|
|
twinscan2gff.pl [-p <transcript_suffix>] < twinscan_output.gtf > gff3 |
| 9 |
|
|
/; |
| 10 |
|
|
getopts('p:') || die($usage."\n"); |
| 11 |
|
|
my $prefix=$Getopt::Std::opt_p || "winscan"; |
| 12 |
|
|
my $curgene; #current gene name |
| 13 |
|
|
my $curtran; #current transcript |
| 14 |
|
|
my ($curstrand, $curchr); #chromosome and strand for current model |
| 15 |
|
|
my @exd; #exons for current model |
| 16 |
|
|
my %stops; #transcripts with stop_codon info transcript_id=>stop_start |
| 17 |
|
|
|
| 18 |
|
|
while (<>) { |
| 19 |
|
|
next if m/^\s*#/; |
| 20 |
|
|
chomp; |
| 21 |
|
|
next unless $_; |
| 22 |
|
|
my ($chr, $prog, $ftype, $fstart, $fend, $fscore, $strand, $frame, $info)= |
| 23 |
|
|
split(/\t/); |
| 24 |
|
|
my ($gid)=($info=~m/gene_id\s+\"([\w\.\:\;\'\|\#\-]+)\"/); |
| 25 |
|
|
my ($tid)=($info=~m/transcript_id\s+\"([\w\.\:\;\'\|\#\-]+)\"/); |
| 26 |
|
|
die ("Error parsing gene/transcript id from $_\n") unless $gid && $tid; |
| 27 |
|
|
if ($ftype eq 'stop_codon') { |
| 28 |
|
|
$stops{$tid}=($strand eq '-')? $fend : $fstart; |
| 29 |
|
|
next; |
| 30 |
|
|
} |
| 31 |
|
|
next unless $ftype eq 'CDS'; |
| 32 |
|
|
|
| 33 |
|
|
if ($tid ne $curtran) { |
| 34 |
|
|
&writeTranscript() if $curtran; |
| 35 |
|
|
$curtran=$tid; |
| 36 |
|
|
$curchr=$chr; |
| 37 |
|
|
$curstrand=$strand; |
| 38 |
|
|
@exd=(); |
| 39 |
|
|
} |
| 40 |
|
|
push(@exd, [$fstart, $fend, $frame]); |
| 41 |
|
|
} |
| 42 |
|
|
|
| 43 |
|
|
&writeTranscript() if $curtran; |
| 44 |
|
|
|
| 45 |
|
|
|
| 46 |
|
|
sub writeTranscript { |
| 47 |
|
|
my @ex= sort { $main::a->[0] <=> $main::b->[0] } @exd; |
| 48 |
|
|
my ($mstart, $mend)=($ex[0]->[0], $ex[-1]->[1]); |
| 49 |
|
|
if (my $stpos=$stops{$curtran}) { |
| 50 |
|
|
if ($curstrand eq '-') { |
| 51 |
|
|
die ("Invalid stop codon position for $curtran! ($stpos vs $mstart)\n") |
| 52 |
|
|
unless $stpos==$mstart-1; |
| 53 |
|
|
$mstart-=3; #include stop codon (twinscan doesn't..) |
| 54 |
|
|
$ex[0]->[0]-=3; |
| 55 |
|
|
} |
| 56 |
|
|
else { |
| 57 |
|
|
die ("Invalid stop codon position for $curtran! ($stpos vs $mend)\n") |
| 58 |
|
|
unless $stpos==$mend+1; |
| 59 |
|
|
$mend+=3; #include stop codon! |
| 60 |
|
|
$ex[-1]->[1]+=3; |
| 61 |
|
|
} |
| 62 |
|
|
} |
| 63 |
|
|
my ($tnum, $anum)=($curtran=~m/\.(\d+)\.(\d+)$/); |
| 64 |
|
|
my $geneid='t'.$prefix.$tnum; |
| 65 |
|
|
$curtran=$geneid."tws$anum"; |
| 66 |
|
|
print join("\t",$curchr, 'twinscan', 'mRNA', $mstart, $mend, '.', |
| 67 |
|
|
$curstrand, '.', "ID=$curtran;Name=$geneid")."\n"; |
| 68 |
|
|
my $i=1; |
| 69 |
|
|
foreach my $exon (@ex) { |
| 70 |
|
|
my ($estart, $eend, $eframe)=@$exon; |
| 71 |
|
|
print join("\t",$curchr, 'twinscan', 'CDS', $estart, $eend, '.', |
| 72 |
|
|
$curstrand, $eframe, |
| 73 |
|
|
"Parent=$curtran")."\n"; |
| 74 |
|
|
$i++; |
| 75 |
|
|
} |
| 76 |
|
|
} |
