1 |
#!/usr/bin/perl |
2 |
use strict; |
3 |
use Getopt::Std; |
4 |
use FindBin;use lib $FindBin::Bin; |
5 |
|
6 |
my $usage = q/Usage: |
7 |
orf_stats.pl <orf_check.tab> <ref_mrnas.iit> |
8 |
|
9 |
It will try to report various stats about the ORFs |
10 |
and junctions of transfrags given in <orf_check.tab> |
11 |
|
12 |
<orf_check.tab> is the output of orf_check.pl, with the columns: |
13 |
<transcript_id>,<xlocus>,<chromosome><strand>,<exons>,<orfs> |
14 |
|
15 |
ORF coordinates are local - i.e. relative to the transcript sequence |
16 |
/; |
17 |
umask 0002; |
18 |
getopts('o:') || die($usage."\n"); |
19 |
my $outfile=$Getopt::Std::opt_o; |
20 |
if ($outfile) { |
21 |
open(OUTF, '>'.$outfile) || die("Error creating output file $outfile\n"); |
22 |
select(OUTF); |
23 |
} |
24 |
# -- |
25 |
my ($forfs, $fiit)=@ARGV; |
26 |
die($usage."\n") unless -f $forfs; # && -f $fiit; |
27 |
my %ci; #cufflinks isoforms data |
28 |
# transcript_id => [$chr, $strand, $cstart, $cend, [ exondata..], [ orfdata.. ]] |
29 |
my $fiitqry=$forfs.'.iitqry'; |
30 |
|
31 |
open(QDATA, '>'.$fiitqry) || die("Error creating iit query file $fiitqry\n"); |
32 |
|
33 |
open(FORFS, $forfs) || die("$usage\n"); |
34 |
while(<FORFS>) { |
35 |
chomp; |
36 |
next unless $_; |
37 |
my ($t, $xloc, $chr, $exonstr, $orfstr)=split(/\t/); |
38 |
my $strand=chop($chr); |
39 |
$exonstr=~s/^exons=//; |
40 |
$orfstr=~s/^orfs=//; |
41 |
my @ex= map { [split(/\-/)] } split(/\,/,$exonstr); |
42 |
my @orfs; |
43 |
if ($orfstr) { |
44 |
@orfs=map { [split(/\:/)] } split(/\,/,$orfstr); |
45 |
} |
46 |
my @t2g; |
47 |
my @gorfs; |
48 |
my $x=0; #position on transcript sequence |
49 |
if ($strand eq '-') { |
50 |
foreach my $d (reverse @ex) { |
51 |
foreach my $gx (reverse($$d[0] .. $$d[1])) { $t2g[++$x]=$gx } |
52 |
} |
53 |
#@gorfs= map { [$t2g[$$_[0]+$$_[1]]-2, $t2g[$$_[0]]] } @orfs; |
54 |
# keep the order so we can more easily match GFF CDS end coordinate |
55 |
if ($orfstr) { |
56 |
@gorfs= map { [$t2g[$$_[0]], $t2g[$$_[0]+$$_[1]]-2] } @orfs; |
57 |
} |
58 |
} |
59 |
else { |
60 |
foreach my $d (@ex) { |
61 |
foreach my $gx ($$d[0] .. $$d[1]) { $t2g[++$x]=$gx } |
62 |
} |
63 |
if ($orfstr) { |
64 |
@gorfs= map { [$t2g[$$_[0]], $t2g[$$_[0]+$$_[1]]+2] } @orfs; |
65 |
} |
66 |
} |
67 |
# - debug print: |
68 |
#print join("\t", $t, $chr.$strand, join(',', (map { $$_[0].'-'.$$_[1] } @ex)), |
69 |
# join(',', (map { $$_[0].'|'.$$_[1] } @gorfs)))."\n"; |
70 |
print QDATA join(' ',$ex[0]->[0], $ex[-1]->[1], $chr.$strand, '<'.$t)."\n"; |
71 |
$ci{$t}=[$chr, $strand, $ex[0]->[0], $ex[-1]->[1], [@ex], [@gorfs]]; |
72 |
} |
73 |
close(QDATA); |
74 |
close(FORFS); |
75 |
# -- now query the iit_file with $fiitqry |
76 |
goto ENDP unless ($fiit && -f $fiit); |
77 |
|
78 |
my $cmd="iit_get $fiit < $fiitqry"; |
79 |
open(IITGET, "$cmd |") || die ("Error opening pipe $cmd |\n"); |
80 |
my @refs; #refs overlapping query transfrag -- forming this cluster |
81 |
my $tf; # current transfrag; |
82 |
my $strand; #strand for current cluster |
83 |
my $lref; #last ref id |
84 |
my $lchrstrand; #last ref chr and strand |
85 |
my %accs; #all acceptor positions for this cluster |
86 |
my %dons; #all donor positions for this cluster |
87 |
my %cdstops; #all CDS stop locations |
88 |
my $jrecomb; # counter for transfrags that have all junctions present in the cluster |
89 |
my $samestop; #counter for transfrags that have an orf stopping |
90 |
# at a known stop in the cluster (e.g. in %cdstops) |
91 |
my $haveORFs; |
92 |
while (<IITGET>) { |
93 |
chomp; |
94 |
next unless $_; |
95 |
#print ">> $_\n"; |
96 |
if (m/^# Query: /) { |
97 |
($strand, $tf)=(m/([\+\-]) <([^<]+)$/); |
98 |
die("Error parsing current transfrag in: $_\n") |
99 |
unless $strand && $tf && exists($ci{$tf}); |
100 |
$haveORFs=(@{$ci{$tf}->[5]}>0); |
101 |
next; |
102 |
} |
103 |
if (m/>(\S+) \d+ \d+ (\S+)/) { |
104 |
#header for a ref interval |
105 |
($lref, $lchrstrand)=($1, $2); |
106 |
} |
107 |
elsif (m/^\d+\-\d+/) { |
108 |
#exons line: |
109 |
if ($strand eq '-') { |
110 |
my @acc=(m/\-(\d+)\,/g); |
111 |
@accs{@acc}=(); |
112 |
my @don=(m/\,(\d+)\-/g); |
113 |
@dons{@don}=(); |
114 |
} |
115 |
else { |
116 |
my @acc=(m/\,(\d+)\-/g); |
117 |
@accs{@acc}=(); |
118 |
my @don=(m/\-(\d+)\,/g); |
119 |
@dons{@don}=(); |
120 |
} |
121 |
} |
122 |
elsif($haveORFs && m/^C:\d+\-/) { |
123 |
# CDS line |
124 |
my $cdstop; |
125 |
if ($strand eq '-') { |
126 |
($cdstop)=(m/^C:(\d+)/); |
127 |
#print STDERR ">>>> storing stop $cdstop\n"; |
128 |
} |
129 |
else { |
130 |
($cdstop)=(m/(\d+)$/); |
131 |
} |
132 |
$cdstops{$cdstop}=1; |
133 |
} |
134 |
if (m/^# End$/) { |
135 |
# collect stats here |
136 |
my ($tchr, $tstrand, $tmin, $tmax, $texs, $torfs)=@{$ci{$tf}}; |
137 |
die ("ChrStrand mismatch!\n") unless $tstrand eq $strand; |
138 |
my @acc; |
139 |
my @don; |
140 |
if ($tstrand eq '-') { |
141 |
@acc= map { $$_[1] } @$texs; |
142 |
pop(@acc); |
143 |
@don= map { $$_[0] } @$texs; |
144 |
shift(@don); |
145 |
} |
146 |
else { |
147 |
@don= map { $$_[1] } @$texs; |
148 |
pop(@don); |
149 |
@acc= map { $$_[0] } @$texs; |
150 |
shift(@acc); |
151 |
} |
152 |
# check if all acceptors and all donors are in the cluster |
153 |
my $alljmatch=0; |
154 |
my $print=0; |
155 |
my @pcol=('',''); |
156 |
foreach my $a (@acc) { |
157 |
goto JNOTFOUND if !exists($accs{$a}); |
158 |
} |
159 |
foreach my $d (@don) { |
160 |
goto JNOTFOUND if !exists($dons{$d}); |
161 |
} |
162 |
$alljmatch=1; |
163 |
JNOTFOUND: |
164 |
if ($alljmatch) { |
165 |
#print "$tr\tjrecomb\n"; |
166 |
$jrecomb++; |
167 |
$print=1; |
168 |
$pcol[0]='jrecomb'; |
169 |
} |
170 |
#now check if same stop codon |
171 |
foreach my $od (@$torfs) { |
172 |
#print STDERR ">>>> >> checking stop $$od[1]\n" if $strand eq '-'; |
173 |
if (exists($cdstops{$$od[1]})) { |
174 |
$print=1; |
175 |
$pcol[1]='samestop'; |
176 |
$samestop++; |
177 |
last; |
178 |
} |
179 |
} |
180 |
print join("\t",$tf,@pcol)."\n" if $print; |
181 |
# -- clear for next run |
182 |
($strand, $tf)=(); |
183 |
%accs=(); |
184 |
%dons=(); |
185 |
%cdstops=(); |
186 |
@refs=(); |
187 |
} |
188 |
} |
189 |
close(IITGET); |
190 |
print STDERR "$jrecomb transfrags are recombinations of known junctions.\n"; |
191 |
print STDERR "$samestop transfrags have an ORF that match a known stop codon.\n"; |
192 |
|
193 |
# -- |
194 |
ENDP: |
195 |
if ($outfile) { |
196 |
select(STDOUT); |
197 |
close(OUTF); |
198 |
} |
199 |
|
200 |
#************ Subroutines ************** |