| 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 ************** |
