[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.2 - (view) (download) (as text)

1 : olson 1.2 #
2 :     # Copyright (c) 2003-2006 University of Chicago and Fellowship
3 :     # for Interpretations of Genomes. All Rights Reserved.
4 :     #
5 :     # This file is part of the SEED Toolkit.
6 :     #
7 :     # The SEED Toolkit is free software. You can redistribute
8 :     # it and/or modify it under the terms of the SEED Toolkit
9 :     # Public License.
10 :     #
11 :     # You should have received a copy of the SEED Toolkit Public License
12 :     # along with this program; if not write to the University of Chicago
13 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
14 :     # Genomes at veronika@thefig.info or download a copy from
15 :     # http://www.theseed.org/LICENSE.TXT.
16 :     #
17 :    
18 : efrank 1.1 use FIG;
19 :    
20 :    
21 :     # usage: verify_pegs_in_tbl Contigs < tbl > report
22 :    
23 :     ($contigs = shift @ARGV)
24 :     || die "usage: verify_pegs_in_tbl Contigs < tbl > report";
25 :    
26 :     open(CONTIGS,$contigs) || die "could not open $contigs";
27 :     $/ = "\n>";
28 :     while (defined($_ = <CONTIGS>))
29 :     {
30 :     chomp;
31 :     if ($_ =~ /^>?(\S+)[^\n]*\n(.*)/s)
32 :     {
33 :     $id = $1;
34 :     $seq = $2;
35 :     $seq =~ s/\s//g;
36 :     $contigs{$id} = uc $seq;
37 :     }
38 :     }
39 :     close(CONTIGS);
40 :     $/ = "\n";
41 :    
42 :     while (defined($_ = <STDIN>))
43 :     {
44 :     if (($_ =~ /^(fig\|\d+\.\d+\.peg\.\d+)\t(\S+)/) &&
45 :     (($rid,$loc) = ($1,$2)) &&
46 :     ($loc =~ /^([^,]+)_(\d+)_\d+(,|$)/) &&
47 :     (($contig,$beg) = ($1,$2)) &&
48 :     ($loc =~ /$contig\_\d+_(\d+)$/) &&
49 :     ($end = $1))
50 :     {
51 :     $seq = &FIG::extract_seq(\%contigs,$loc);
52 :     if ($contigs{$contig} && ($beg < $end))
53 :     {
54 :     $start = substr($contigs{$contig},$beg-1,3);
55 :     $stop = substr($contigs{$contig},$end,3);
56 :     for ($i=0; ($i < (length($seq) - 3)) && (substr($seq,$i,3) !~ /(TAG|TGA|TAA)/); $i += 3) {}
57 :     if ($i < (length($seq)-3))
58 :     {
59 :     print "$rid has an embedded stop codon\n";
60 :     }
61 :     }
62 :     elsif ($contigs{$contig})
63 :     {
64 :     $start = &reverse_comp(substr($contigs{$contig},$beg-3,3));
65 :     $stop = &reverse_comp(substr($contigs{$contig},$end-4,3));
66 :     for ($i=0; ($i < (length($seq) - 3)) && (substr($seq,$i,3) !~ /(TAG|TGA|TAA)/); $i += 3) {}
67 :     if ($i < (length($seq)-3))
68 :     {
69 :     print "$rid has an embedded stop codon\n";
70 :     }
71 :     }
72 :     else
73 :     {
74 :     print STDERR "invalid: $_";
75 :     }
76 :    
77 :     if (($start !~ /^[ATG]TG/) && ($beg > 40) && ($beg < (length($contigs{$contig})-40)))
78 :     {
79 :     print "bad start for $rid: $start\n";
80 :     }
81 :     if (($stop !~ /(TAG|TGA|TAA)/) && ($end > 40) && ($end < (length($contigs{$contig}) - 40)))
82 :     {
83 :     print "bad stop for $rid: $stop\n";
84 :     }
85 :     $strand = ($beg < $end) ? "+" : "-";
86 :    
87 :     push(@{ $on{$contig} }, [&FIG::min($beg,$end),&FIG::max($beg,$end),$strand,$rid]);
88 :     }
89 :     }
90 :    
91 :     foreach $contig (sort keys(%on))
92 :     {
93 :     $x = $on{$contig};
94 :     @on = sort { ($a->[0] <=> $b->[0]) or ($a->[1] <=> $b->[1]) or ($a->[2] cmp $b->[2]) } @$x;
95 :    
96 :     for ($i=0; ($i < $#on); $i++)
97 :     {
98 :     for ($j=$i+1; ($j < @on) && ($on[$i]->[1] > ($on[$j]->[0]) + 10); $j++)
99 :     {
100 :     $overlap = ($on[$j]->[1] > $on[$i]->[1]) ?
101 :     ($on[$i]->[1] - $on[$j]->[0]) + 1 :
102 :     ($on[$j]->[1] - $on[$j]->[0]) + 1;
103 :     push(@overlaps,[$overlap,"$overlap\t$on[$i]->[2],$on[$j]->[2]\t$on[$i]->[3] overlaps $on[$j]->[3]\n"]);
104 :     }
105 :     }
106 :     }
107 :    
108 :     foreach $_ (sort { $b->[0] <=> $a->[0] } @overlaps)
109 :     {
110 :     print $_->[1];
111 :     }
112 :    
113 :     sub reverse_comp {
114 :     my( $seq ) = @_;
115 :     my( $rev );
116 :    
117 :     $rev = reverse( $seq );
118 :     $rev =~ tr/a-z/A-Z/;
119 :     $rev =~ tr/ACGTUMRWSYKBDHV/TGCAAKYWSRMVHDB/;
120 :     return $rev;
121 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3