[Bio] / FigKernelScripts / verify_pegs_in_tbl.pl Repository:
ViewVC logotype

Annotation of /FigKernelScripts/verify_pegs_in_tbl.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (view) (download) (as text)

1 : efrank 1.1 use FIG;
2 :    
3 :    
4 :     # usage: verify_pegs_in_tbl Contigs < tbl > report
5 :    
6 :     ($contigs = shift @ARGV)
7 :     || die "usage: verify_pegs_in_tbl Contigs < tbl > report";
8 :    
9 :     open(CONTIGS,$contigs) || die "could not open $contigs";
10 :     $/ = "\n>";
11 :     while (defined($_ = <CONTIGS>))
12 :     {
13 :     chomp;
14 :     if ($_ =~ /^>?(\S+)[^\n]*\n(.*)/s)
15 :     {
16 :     $id = $1;
17 :     $seq = $2;
18 :     $seq =~ s/\s//g;
19 :     $contigs{$id} = uc $seq;
20 :     }
21 :     }
22 :     close(CONTIGS);
23 :     $/ = "\n";
24 :    
25 :     while (defined($_ = <STDIN>))
26 :     {
27 :     if (($_ =~ /^(fig\|\d+\.\d+\.peg\.\d+)\t(\S+)/) &&
28 :     (($rid,$loc) = ($1,$2)) &&
29 :     ($loc =~ /^([^,]+)_(\d+)_\d+(,|$)/) &&
30 :     (($contig,$beg) = ($1,$2)) &&
31 :     ($loc =~ /$contig\_\d+_(\d+)$/) &&
32 :     ($end = $1))
33 :     {
34 :     $seq = &FIG::extract_seq(\%contigs,$loc);
35 :     if ($contigs{$contig} && ($beg < $end))
36 :     {
37 :     $start = substr($contigs{$contig},$beg-1,3);
38 :     $stop = substr($contigs{$contig},$end,3);
39 :     for ($i=0; ($i < (length($seq) - 3)) && (substr($seq,$i,3) !~ /(TAG|TGA|TAA)/); $i += 3) {}
40 :     if ($i < (length($seq)-3))
41 :     {
42 :     print "$rid has an embedded stop codon\n";
43 :     }
44 :     }
45 :     elsif ($contigs{$contig})
46 :     {
47 :     $start = &reverse_comp(substr($contigs{$contig},$beg-3,3));
48 :     $stop = &reverse_comp(substr($contigs{$contig},$end-4,3));
49 :     for ($i=0; ($i < (length($seq) - 3)) && (substr($seq,$i,3) !~ /(TAG|TGA|TAA)/); $i += 3) {}
50 :     if ($i < (length($seq)-3))
51 :     {
52 :     print "$rid has an embedded stop codon\n";
53 :     }
54 :     }
55 :     else
56 :     {
57 :     print STDERR "invalid: $_";
58 :     }
59 :    
60 :     if (($start !~ /^[ATG]TG/) && ($beg > 40) && ($beg < (length($contigs{$contig})-40)))
61 :     {
62 :     print "bad start for $rid: $start\n";
63 :     }
64 :     if (($stop !~ /(TAG|TGA|TAA)/) && ($end > 40) && ($end < (length($contigs{$contig}) - 40)))
65 :     {
66 :     print "bad stop for $rid: $stop\n";
67 :     }
68 :     $strand = ($beg < $end) ? "+" : "-";
69 :    
70 :     push(@{ $on{$contig} }, [&FIG::min($beg,$end),&FIG::max($beg,$end),$strand,$rid]);
71 :     }
72 :     }
73 :    
74 :     foreach $contig (sort keys(%on))
75 :     {
76 :     $x = $on{$contig};
77 :     @on = sort { ($a->[0] <=> $b->[0]) or ($a->[1] <=> $b->[1]) or ($a->[2] cmp $b->[2]) } @$x;
78 :    
79 :     for ($i=0; ($i < $#on); $i++)
80 :     {
81 :     for ($j=$i+1; ($j < @on) && ($on[$i]->[1] > ($on[$j]->[0]) + 10); $j++)
82 :     {
83 :     $overlap = ($on[$j]->[1] > $on[$i]->[1]) ?
84 :     ($on[$i]->[1] - $on[$j]->[0]) + 1 :
85 :     ($on[$j]->[1] - $on[$j]->[0]) + 1;
86 :     push(@overlaps,[$overlap,"$overlap\t$on[$i]->[2],$on[$j]->[2]\t$on[$i]->[3] overlaps $on[$j]->[3]\n"]);
87 :     }
88 :     }
89 :     }
90 :    
91 :     foreach $_ (sort { $b->[0] <=> $a->[0] } @overlaps)
92 :     {
93 :     print $_->[1];
94 :     }
95 :    
96 :     sub reverse_comp {
97 :     my( $seq ) = @_;
98 :     my( $rev );
99 :    
100 :     $rev = reverse( $seq );
101 :     $rev =~ tr/a-z/A-Z/;
102 :     $rev =~ tr/ACGTUMRWSYKBDHV/TGCAAKYWSRMVHDB/;
103 :     return $rev;
104 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3