[Bio] / StandaloneTools / verify_genome_directory_standalone Repository:
ViewVC logotype

Annotation of /StandaloneTools/verify_genome_directory_standalone

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (view) (download)

1 : olson 1.1 #!/usr/bin/perl -w
2 :    
3 :     # usage: verify_genome_directory Dir
4 :    
5 :     ($dir = shift @ARGV)
6 :     || die "\n usage: verify_genome_directory Dir\n\n";
7 :    
8 :     if (!-d $dir)
9 :     {
10 :     die "Organism directory $dir does not exist";
11 :     }
12 :    
13 :     if (! (-s "$dir/GENOME")) { print "WARNING: missing GENOME file\n" } else { $genome = `cat $dir/GENOME`; chomp $genome; }
14 :     if (! (-s "$dir/PROJECT")) { print "WARNING: missing PROJECT file\n" }
15 :     if (! (-s "$dir/TAXONOMY")) { print "WARNING: missing TAXONOMY file\n" } else { $taxonomy = `cat $dir/TAXONOMY`; chomp $taxonomy; }
16 :    
17 :     $taxonomy =~ m/^([^;]+)/;
18 :     $domain = $1;
19 :    
20 :     if (!-d "$dir/Features")
21 :     {
22 :     die "ERROR: Features directory $dir/Features does not exist";
23 :     }
24 :     else
25 :     {
26 :     if (! opendir(TMP,$dir))
27 :     {
28 :     ++$bad{no_org_dir};
29 :     print "ERROR: could not open $dir";
30 :     }
31 :     else
32 :     {
33 :     @contigs = grep { $_ =~ /contig/ } readdir(TMP);
34 :     closedir(TMP) || warn "Could not close directory $dir";
35 :    
36 :     if (@contigs == 0)
37 :     {
38 :     print "missing contigs\n";
39 :     }
40 :     else
41 :     {
42 :     $size = 0;
43 :     foreach $file (@contigs)
44 :     {
45 :     if (!open(TMP, "<$dir/$file"))
46 :     {
47 :     ++$bad{unopenable_contigs};
48 :     print "ERROR: Could not read-open $dir/$file";
49 :     }
50 :     else
51 :     {
52 :     while (($id, $seqP) = &read_fasta_record(\*TMP))
53 :     {
54 :     $size += $len_of_contig{$id} = length($$seqP);
55 :     }
56 :     }
57 :     }
58 :    
59 :     if ((-s "$dir/COMPLETE") && ($size < 300000))
60 :     {
61 :     print "WARNING: contigs for $dir contain only $size bp of data\n";
62 :     }
63 :    
64 :     if (!-d "$dir/Features/peg")
65 :     {
66 :     ++$bad{peg_no_directory};
67 :     print "ERROR: No peg directory $dir/Features/peg\n";
68 :     }
69 :     else
70 :     {
71 :     if (!-e "$dir/Features/peg/tbl")
72 :     {
73 :     ++$bad{peg_no_tbl};
74 :     print "ERROR: No peg tbl file $dir/Features/peg/tbl\n";
75 :     }
76 :     else
77 :     {
78 :     if (!-s "$dir/Features/peg/tbl")
79 :     {
80 :     ++$bad{peg_zero_sz_tbl};
81 :     print "ERROR: Zero-size peg tbl file $dir/Features/peg/tbl\n";
82 :     }
83 :     else
84 :     {
85 :     if (!open(TMP, "<$dir/Features/peg/tbl"))
86 :     {
87 :     ++$bad{peg_unopenable_tbl};
88 :     print "ERROR: Could not read-open $dir/Features/peg/tbl";
89 :     }
90 :     else
91 :     {
92 :     while (defined($entry = <TMP>))
93 :     {
94 :     chomp $entry;
95 :     ($id, $loc) = split(/\t/, $entry);
96 :     $loc_of_peg{$id} = $loc;
97 :     foreach $exon (split(/,/, $loc))
98 :     {
99 :     $exon =~ m/^(\S+)_(\d+)_(\d+)$/;
100 :     ($contig, $beg, $end) = ($1, $2, $3);
101 :    
102 :     if ($taxonomy !~ m/Eukaryota/)
103 :     {
104 :     if (not defined($len_of_contig{$contig}))
105 :     {
106 :     ++$bad{peg_undef_contig_ref};
107 :     print "ERROR: for $id, $loc refers to undefined contig $contig\n";
108 :     }
109 :     else
110 :     {
111 :     if (&max($beg,$end) > $len_of_contig{$contig})
112 :     {
113 :     ++$bad{peg_out_of_bounds};
114 :     print "ERROR: for $id, $exon in $loc exceeds contig length $len_of_contig{$contig}\n";
115 :     }
116 :     }
117 :     }
118 :     }
119 :     }
120 :     close(TMP) || warn "could not close $dir/Features/peg/tbl";
121 :     }
122 :     }
123 :    
124 :     if (!open(TMP, "<$dir/Features/peg/fasta"))
125 :     {
126 :     ++$bad{peg_unopenable_fasta};
127 :     print "ERROR: Could not read-open $dir/Features/peg/fasta";
128 :     }
129 :     else
130 :     {
131 :     while (($id, $seqP) = &read_fasta_record(\*TMP))
132 :     {
133 :     $len_of_peg{$id} = length($$seqP);
134 :    
135 :     if (!defined($loc_of_peg{$id}))
136 :     {
137 :     ++$bad{peg_fasta_no_tbl};
138 :     print "ERROR: $id is in peg fasta, but not in tbl\n";
139 :     }
140 :     }
141 :     close(TMP) || warn "could not close $dir/Features/peg/fasta";
142 :     print STDERR "Read ", (scalar keys %len_of_peg), " pegs from $dir/Features/peg/fasta\n"
143 :     if $ENV{SEED_VERBOSE};
144 :     }
145 :    
146 :     foreach $id (sort keys %loc_of_peg)
147 :     {
148 :     if (!defined($len_of_peg{$id}))
149 :     {
150 :     ++$bad{peg_tbl_no_fasta};
151 :     print "ERROR: $id is in peg tbl, but not in fasta\n";
152 :     }
153 :     }
154 :     }
155 :    
156 :     if (!-d "$dir/Features/rna")
157 :     {
158 :     print "WARNING: no rna directory for $dir\n";
159 :     }
160 :     else
161 :     {
162 :     if (!-e "$dir/Features/rna/tbl")
163 :     {
164 :     ++$bad{rna_no_tbl};
165 :     print "ERROR: No rna tbl file $dir/Features/rna/tbl\n";
166 :     }
167 :     else
168 :     {
169 :     if (!-s "$dir/Features/rna/tbl")
170 :     {
171 :     ++$bad{rna_zero_sz_tbl};
172 :     print "ERROR: Zero-size rna tbl file $dir/Features/rna/tbl\n";
173 :     }
174 :     else
175 :     {
176 :     if (!open(TMP, "<$dir/Features/rna/tbl"))
177 :     {
178 :     ++$bad{rna_unopenable_tbl};
179 :     warn "Could not read-open $dir/Features/rna/tbl";
180 :     }
181 :     else
182 :     {
183 :     while (defined($entry = <TMP>))
184 :     {
185 :     chomp $entry;
186 :     ($id, $loc) = split(/\t/, $entry);
187 :     $loc_of_rna{$id} = $loc;
188 :     foreach $exon (split(/,/, $loc))
189 :     {
190 :     $exon =~ m/^(\S+)_(\d+)_(\d+)$/;
191 :     ($contig, $beg, $end) = ($1, $2, $3);
192 :     if (!defined($len_of_contig{$contig}))
193 :     {
194 :     ++$bad{rna_undef_contig_rna};
195 :     print "ERROR: for $id, $loc refers to undefined contig $contig\n";
196 :     }
197 :     else
198 :     {
199 :     if (&max($beg,$end) > $len_of_contig{$contig})
200 :     {
201 :     ++$bad{rna_out_of_bounds};
202 :     print "ERROR: for $id, $exon in $loc exceeds contig length $len_of_contig{$contig}\n";
203 :     }
204 :     }
205 :     }
206 :     }
207 :     close(TMP) || warn "could not close $dir/Features/rna/tbl";
208 :     }
209 :    
210 :     open(TMP, "<$dir/Features/rna/fasta") || warn "Could not read-open $dir/Features/rna/fasta";
211 :     while (($id, $seqP) = &read_fasta_record(\*TMP))
212 :     {
213 :     $len_of_rna{$id} = length($$seqP);
214 :    
215 :     if (!defined($loc_of_rna{$id}))
216 :     {
217 :     ++$bad{rna_fasta_no_tbl};
218 :     print "ERROR: $id is in rna fasta, but not in tbl\n";
219 :     }
220 :     }
221 :     close(TMP) || warn "could not close $dir/Features/rna/fasta";
222 :     print STDERR "Read ", (scalar keys %len_of_rna), " rnas from $dir/Features/rna/fasta\n"
223 :     if $ENV{SEED_VERBOSE};
224 :    
225 :     foreach $id (sort keys %loc_of_rna)
226 :     {
227 :     if (!defined($len_of_rna{$id}))
228 :     {
229 :     ++$bad{rna_tbl_no_fasta};
230 :     print "ERROR: $id is in rna tbl, but not in fasta\n";
231 :     }
232 :     }
233 :     }
234 :     }
235 :     }
236 :     }
237 :     }
238 :     }
239 :     }
240 :    
241 :     if (keys %bad)
242 :     {
243 :     print STDERR "$dir is corrupt (", join(", ", map { "$_=$bad{$_}" } sort keys %bad), ")\t$domain\t$genome\n";
244 :     }
245 :     else
246 :     {
247 :     print STDERR "$dir is O.K.\n";
248 :     }
249 :     # print STDERR "$0 has completed\n";
250 :    
251 :     exit(scalar keys %bad);
252 :    
253 :    
254 :     sub min {
255 :     my(@x) = @_;
256 :     my($min,$i);
257 :    
258 :     (@x > 0) || return undef;
259 :     $min = $x[0];
260 :     for ($i=1; ($i < @x); $i++)
261 :     {
262 :     $min = ($min > $x[$i]) ? $x[$i] : $min;
263 :     }
264 :     return $min;
265 :     }
266 :    
267 :     sub max {
268 :     my(@x) = @_;
269 :     my($max,$i);
270 :    
271 :     (@x > 0) || return undef;
272 :     $max = $x[0];
273 :     for ($i=1; ($i < @x); $i++)
274 :     {
275 :     $max = ($max < $x[$i]) ? $x[$i] : $max;
276 :     }
277 :     return $max;
278 :     }
279 :    
280 :    
281 :     sub read_fasta_record
282 :     {
283 :     my ($file_handle) = @_;
284 :     my ( $old_end_of_record, $fasta_record, @lines, $head, $sequence, $seq_id, $comment, @parsed_fasta_record );
285 :    
286 :     if (not defined($file_handle)) { $file_handle = \*STDIN; }
287 :    
288 :     $old_end_of_record = $/;
289 :     $/ = "\n>";
290 :    
291 :     if (defined($fasta_record = <$file_handle>))
292 :     {
293 :     chomp $fasta_record;
294 :     @lines = split( /\n/, $fasta_record );
295 :     $head = shift @lines;
296 :     $head =~ s/^>?//;
297 :     $head =~ m/^(\S+)/;
298 :     $seq_id = $1;
299 :    
300 :     if ($head =~ m/^\S+\s+(.*)$/) { $comment = $1; } else { $comment = ""; }
301 :    
302 :     $sequence = join( "", @lines );
303 :    
304 :     @parsed_fasta_record = ( $seq_id, \$sequence, $comment );
305 :     }
306 :     else
307 :     {
308 :     @parsed_fasta_record = ();
309 :     }
310 :    
311 :     $/ = $old_end_of_record;
312 :    
313 :     return @parsed_fasta_record;
314 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3