1 |
gpertea |
23 |
#!/usr/bin/perl |
2 |
|
|
use strict; |
3 |
|
|
use Getopt::Std; |
4 |
|
|
use FindBin;use lib $FindBin::Bin; |
5 |
|
|
|
6 |
|
|
my $usage = q/Usage: |
7 |
|
|
utrpaste.pl [-a 'newID'] [-o <output.gff3>] <jigsaw_out.jgff> <exon_mappings.gff3> |
8 |
|
|
|
9 |
|
|
Augment the gene models found in <jigsaw_out.jgff> with |
10 |
|
|
exon mappings found in <mappings.gff3> (only initial\/terminal exons, |
11 |
|
|
if they overlap the CDS) |
12 |
|
|
|
13 |
|
|
/; |
14 |
|
|
umask 0002; |
15 |
|
|
getopts('a:o:') || die($usage."\n"); |
16 |
|
|
my $outfile=$Getopt::Std::opt_o; |
17 |
|
|
my $newacc=$Getopt::Std::opt_a; |
18 |
|
|
my ($jgff,$mapgff)=@ARGV; |
19 |
|
|
foreach (@ARGV) { |
20 |
|
|
die ("Error locating file $_\n") unless -f $_; |
21 |
|
|
} |
22 |
|
|
open(MGFF, $mapgff) || die ("Error opening mappings file $mapgff\n"); |
23 |
|
|
# 0 1 2 3 4 |
24 |
|
|
my @mapdata; #list of [ ID, strand, start, end, [ @exons] ] |
25 |
|
|
# where each exon element is [exon_start, exon_end] |
26 |
|
|
my $mi=-1; #current index in @mapdata |
27 |
|
|
my $curid; |
28 |
|
|
while (<MGFF>) { |
29 |
|
|
my ($chr, $v, $f, $fstart, $fend, $fscore, $strand, $frame, $info)=split(/\t/); |
30 |
|
|
next unless ($f eq 'exon'); #! ignore anything else but exons |
31 |
|
|
($fstart, $fend)=($fend, $fstart) if ($fstart>$fend); |
32 |
|
|
my ($id)=($info=~m/Parent=([^;]+)/); |
33 |
|
|
die("Error getting mapping ID at:\n$_\n") unless $id; |
34 |
|
|
if ($id ne $curid) { |
35 |
|
|
#start a new record; |
36 |
|
|
push(@mapdata, [ $id, $strand, 0,0, []]); |
37 |
|
|
$mi=$#mapdata; |
38 |
|
|
$curid=$id; |
39 |
|
|
} |
40 |
|
|
push( @{$mapdata[$mi]->[4]} , [$fstart, $fend] ); |
41 |
|
|
} |
42 |
|
|
close(MGFF); |
43 |
|
|
|
44 |
|
|
foreach my $md (@mapdata) { |
45 |
|
|
my @exons=sort { $main::a->[0] <=> $main::b->[0] } @{$md->[4]}; |
46 |
|
|
$md->[2]=$exons[0]->[0]; |
47 |
|
|
$md->[3]=$exons[-1]->[1]; |
48 |
|
|
$md->[4]=[@exons]; |
49 |
|
|
} |
50 |
|
|
|
51 |
|
|
# sort mappings by start coordinate |
52 |
|
|
@mapdata=sort { $main::a->[2] <=> $main::b->[2] } @mapdata; |
53 |
|
|
|
54 |
|
|
if ($outfile) { |
55 |
|
|
open(FOUT, '>'.$outfile) || die("Error creating file $outfile\n"); |
56 |
|
|
select(FOUT); |
57 |
|
|
} |
58 |
|
|
#now parse and augment the jgff |
59 |
|
|
open(JGFF, $jgff) || die ("Error opening mappings file $jgff\n"); |
60 |
|
|
$curid=''; |
61 |
|
|
my @cds=0; #current CDS chain: list of [$start, $end, $frame] |
62 |
|
|
my @curdata; |
63 |
|
|
while (<JGFF>) { |
64 |
|
|
next if length($_)<4 || m/^# /; |
65 |
|
|
my ($chr, $track, $f, $fstart, $fend, $fscore, $strand, $phase, $jdata)=split(/\t/); |
66 |
|
|
next unless (($f eq 'CDS') || ($track=~m/jigsaw/ && $f=~m/\-exon$/)); |
67 |
|
|
my ($id)=($jdata=~m/^([^;]+); gene_score=/); |
68 |
|
|
unless ($id) { |
69 |
|
|
($id)=($jdata=~m/Parent=([^;]+)/); |
70 |
|
|
unless ($id) { #try GTF too |
71 |
|
|
($id)=($jdata=~m/transcript_id "([^"]+)/); |
72 |
|
|
} |
73 |
|
|
die("Error getting model ID at:\n$_\n") unless $id; |
74 |
|
|
} |
75 |
|
|
if ($id ne $curid) { |
76 |
|
|
extendCDS(@curdata, \@cds) if $curid; #checks @mapdata |
77 |
|
|
@cds=(); |
78 |
|
|
$curid=$id; |
79 |
|
|
my $wid = $newacc || $curid; |
80 |
|
|
my ($jnum)=($curid=~m/\.(\d+)$/); |
81 |
|
|
$jnum=1 unless $jnum; |
82 |
|
|
@curdata=($wid, $jnum, $chr, $track, $strand); |
83 |
|
|
} |
84 |
|
|
push(@cds, [$fstart, $fend, $phase]); |
85 |
|
|
} |
86 |
|
|
# last record: |
87 |
|
|
extendCDS(@curdata, \@cds) if @cds>0; #checks @mapdata |
88 |
|
|
close(JGFF); |
89 |
|
|
|
90 |
|
|
if ($outfile) { close(FOUT); } |
91 |
|
|
|
92 |
|
|
#================ SUBROUTINES ============ |
93 |
|
|
|
94 |
|
|
sub extendCDS { |
95 |
|
|
my ($id, $jnum, $chr, $track, $strand, $cdsl)=@_; |
96 |
|
|
#find mapdata overlapping first and last exons of this model |
97 |
|
|
my @cds = sort { $main::a->[0] <=> $main::b->[0] } @$cdsl; |
98 |
|
|
my $cstart=$cds[0]->[0]; |
99 |
|
|
my $cend=$cds[-1]->[1]; |
100 |
|
|
my $mrnaid=$id.'.j'.$jnum; |
101 |
|
|
|
102 |
|
|
my $found=0; |
103 |
|
|
my @exons; |
104 |
|
|
|
105 |
|
|
foreach my $m (@mapdata) { |
106 |
|
|
my ($mid, $mstrand, $mstart, $mend, $mexons)=@$m; |
107 |
|
|
#mapping should overlap on the same strand |
108 |
|
|
#print STDERR "===> comparing $mstart-$mend ($mstrand) w $cstart-$cend ($strand)\n"; |
109 |
|
|
next if ($mstrand ne $strand || $mstart>$cend || $cstart>$mend); |
110 |
|
|
#mapping ends should cover a larger interval than CDS |
111 |
|
|
next if (($mstart>$cstart) || ($mend<$cend)); |
112 |
|
|
# at least one of the ends should stretch farther |
113 |
|
|
next unless (($mstart<$cstart) || ($mend>$cend)); |
114 |
|
|
# -- found it: |
115 |
|
|
$found=1; |
116 |
|
|
# augment here, and print the mRNA line too, then exons |
117 |
|
|
# left end : |
118 |
|
|
my @lexons; |
119 |
|
|
my $lovl; |
120 |
|
|
for (my $i=0;$i<@$mexons;$i++) { |
121 |
|
|
my ($estart, $eend)=@{$$mexons[$i]}; |
122 |
|
|
#print STDERR "left end cmp: $estart-$eend vs CDS $cstart-$cds[0]->[1]\n"; |
123 |
|
|
if ($estart<=$cstart && $eend>$cstart) { |
124 |
|
|
# if they really overlapped properly, $estart is <= cds start |
125 |
|
|
#print STDERR " YES, adding exon $estart-$cds[0]->[1]^^!\n"; |
126 |
|
|
push(@lexons, [$estart, $cds[0]->[1]]); |
127 |
|
|
$lovl=1; |
128 |
|
|
last; |
129 |
|
|
} |
130 |
|
|
push(@lexons,[$estart, $eend]); |
131 |
|
|
} |
132 |
|
|
#right end: |
133 |
|
|
my @rexons; |
134 |
|
|
my $rovl; |
135 |
|
|
for (my $i=@$mexons-1;$i>=0;$i--) { |
136 |
|
|
my ($estart, $eend)=@{$$mexons[$i]}; |
137 |
|
|
#print STDERR "right end cmp: $estart-$eend vs CDS $cds[-1]->[0]-$cend\n"; |
138 |
|
|
if ($eend>=$cend && $estart<$cend) { |
139 |
|
|
# if they really overlapped properly, $eend is >= cds end |
140 |
|
|
#print STDERR " YES, adding exon $cds[-1]->[0]-$eend!\n"; |
141 |
|
|
unshift(@rexons, [$cds[-1]->[0], $eend]); |
142 |
|
|
$rovl=1; |
143 |
|
|
last; |
144 |
|
|
} |
145 |
|
|
unshift(@rexons,[$estart, $eend]); |
146 |
|
|
} |
147 |
|
|
|
148 |
|
|
if ($rovl && $lovl) { |
149 |
|
|
#special single exon case: |
150 |
|
|
if ( @cds==1 ) { |
151 |
|
|
#we have to merge rexon and lexons |
152 |
|
|
$lexons[-1]->[1]=$rexons[0]->[1]; |
153 |
|
|
shift(@rexons); #remove the first rexon |
154 |
|
|
push(@exons, @lexons, @rexons); |
155 |
|
|
} |
156 |
|
|
else { #normal multi-exon case |
157 |
|
|
push(@exons, @lexons); |
158 |
|
|
foreach my $cd (@cds) { |
159 |
|
|
push(@exons, [$$cd[0], $$cd[1]]) |
160 |
|
|
if ($cd->[0]>$lexons[-1]->[1] && $cd->[1]<$rexons[0]->[0]) |
161 |
|
|
} |
162 |
|
|
#for (my $i=1;$i<@cds;$i++) { |
163 |
|
|
# push(@exons, [$cds[$i]->[0], $cds[$i]->[1]]); |
164 |
|
|
# } |
165 |
|
|
push(@exons, @rexons); |
166 |
|
|
} |
167 |
|
|
# print mRNA and exons |
168 |
|
|
print join("\t", $chr, $track, 'mRNA', $exons[0]->[0], $exons[-1]->[1], '.', $strand, '.', |
169 |
|
|
'ID='.$mrnaid)."\n"; |
170 |
|
|
foreach my $e (@exons) { |
171 |
|
|
print join("\t", $chr, $track, 'exon', $$e[0], $$e[1], '.', $strand, '.', |
172 |
|
|
'Parent='.$mrnaid)."\n"; |
173 |
|
|
} |
174 |
|
|
last; |
175 |
|
|
} |
176 |
|
|
else { $found = 0; } |
177 |
|
|
} |
178 |
|
|
if (!$found) { |
179 |
|
|
#print mRNA line and the old CDS as is |
180 |
|
|
print join("\t", $chr, $track, 'mRNA', $cstart, $cend, '.', $strand, '.', |
181 |
|
|
'ID='.$mrnaid)."\n"; |
182 |
|
|
} |
183 |
|
|
#print CDS segments: |
184 |
|
|
foreach my $c (@cds) { |
185 |
|
|
print join("\t", $chr, $track, 'CDS', $$c[0], $$c[1], '.', $strand, $$c[2], |
186 |
|
|
'Parent='.$mrnaid)."\n"; |
187 |
|
|
|
188 |
|
|
} |
189 |
|
|
} |