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

Annotation of /FigKernelScripts/parse_genbank.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : efrank 1.1 use FIG;
2 :    
3 : efrank 1.2 $usage = "usage: parse_genbank Genome Dir < genbank.entry";
4 :     (
5 : efrank 1.3 ($genome = shift(@ARGV)) &&
6 :     ($dir = shift(@ARGV))
7 : efrank 1.1 )
8 :     || die $usage;
9 :    
10 : efrank 1.2 $prefixP = "fig|$genome.peg.";
11 :     $prefixR = "fig|$genome.rna.";
12 :    
13 : efrank 1.3 &FIG::verify_dir("$dir/Features/peg");
14 :     &FIG::verify_dir("$dir/Features/rna");
15 :     open(CONTIGS,">$dir/contigs") || die "could not open $dir/contigs";
16 : efrank 1.2 open(TBLPEG,">$dir/Features/peg/tbl") || die "could not open $dir/Features/peg/tbl";
17 :     open(FASTAPEG,">$dir/Features/peg/fasta") || die "could not open $dir/Features/peg/fasta";
18 : efrank 1.3 open(ASSIGNMENTS,">$dir/assigned_functions") || die "could not open $dir/assigned_functions";
19 : efrank 1.2 open(TBLRNA,">$dir/Features/rna/tbl") || die "could not open $dir/Features/rna/tbl";
20 :     open(FASTARNA,">$dir/Features/rna/fasta") || die "could not open $dir/Features/rna/fasta";
21 :    
22 : overbeek 1.4 open(TAX,">$dir/TAXONOMY") || die "could not open $dir/TAXONOMY";
23 :     open(GENOME,">$dir/GENOME") || die "could not open $dir/GENOME";
24 :     open(PROJECT,">$dir/PROJECT") || die "could not open $dir/PROJECT";
25 :     print PROJECT "Parsed NCBI Reference Sequences\n";
26 :     close(PROJECT);
27 :    
28 :     $idNp = 1;
29 :     $idNr = 1;
30 :    
31 : efrank 1.1 $/ = "\nLOCUS";
32 :     while ($contig = <STDIN>)
33 :     {
34 :     if ($contig =~ /\nACCESSION\s+(\S+)/s)
35 :     {
36 :     $id = $1;
37 : gdpusch 1.5 if ($contig =~ /\nORIGIN([acgtumrwsykbdhvnxACGTUMRWSYKBDHVNX0-9\s\n]*)(\n\/\/\n|LOCUS)/s)
38 : efrank 1.1 {
39 :     $seq = $1;
40 :     $seq =~ s/\s//gs;
41 :     $seq =~ s/\d//g;
42 :     undef $contigs;
43 :     $contigs->{$id} = $seq;
44 :     &FIG::display_id_and_seq($id,\$seq,\*CONTIGS);
45 :     }
46 :     else
47 :     {
48 : overbeek 1.4 warn "could not find sequence for $id\n";
49 :     next;
50 :     }
51 :    
52 :     if ($contig =~ /\nSOURCE\s{4,8}(\S[^\n]+\S)\n\s+ORGANISM\s+(\S[^\n]+\S)((\n\s{10,14}\S[^\n]+\S)+)/s)
53 :     {
54 :     $genome = $2;
55 :     $tax = $3;
56 :     $tax =~ s/\n\s+//g;
57 :     if (! $written_genome)
58 :     {
59 :     print GENOME "$genome\n";
60 :     close(GENOME);
61 :     print TAX "$tax\n";
62 :     close(TAX);
63 :    
64 :     $written_genome = "$tax\t$genome";
65 :     }
66 :     elsif ($written_genome ne "$tax\t$genome")
67 :     {
68 :     print STDERR "serious mismatch in GENOME/TAX for $id\n$written_genome\n$tax\t$genome\n\n";
69 :     }
70 : efrank 1.1 }
71 :    
72 : overbeek 1.6 while ($contig =~ /\n\s{4,6}CDS\s+([^\n]+(\n {20,}[^\n]*)+)/gs)
73 : efrank 1.1 {
74 :     $cds = $1;
75 : overbeek 1.4 &process_cds($id,\$cds,$prefixP,\$idNp,$contigs,\*TBLPEG,\*FASTAPEG,\*ASSIGNMENTS);
76 : efrank 1.2 }
77 :    
78 :     while ($contig =~ /\n\s{3,6}[tr]RNA\s+([^\n]+(\n\s{20,22}\S[^\n]*)+)/gs)
79 :     {
80 :     $rna = $1;
81 :     &process_rna($id,\$rna,$prefixR,\$idNr,$contigs,\*TBLRNA,\*FASTARNA);
82 : efrank 1.1 }
83 :     }
84 :     }
85 : efrank 1.2 close(CONTIGS);
86 :     close(TBLPEG);
87 : overbeek 1.4 close(FASTAPEG);
88 : efrank 1.2 close(ASSIGNMENTS);
89 :     close(TBLRNA);
90 :     close(FASTARNA);
91 :    
92 :     if (! -s "$dir/contigs") { unlink("$dir/contigs"); print STDERR "no contigs in $dir\n"; }
93 : efrank 1.3 if (! -s "$dir/assigned_functions") { unlink("$dir/assigned_functions"); print STDERR "no assigned_functions in $dir\n"; }
94 :     if (! -s "$dir/Features/peg/tbl") { system "rm -rf $dir/Features/peg"; print STDERR "no PEGs in $dir\n"; }
95 :     if (! -s "$dir/Features/rna/tbl") { system "rm -rf $dir/Features/rna"; print STDERR "no RNAs in $dir\n"; }
96 : efrank 1.2
97 : efrank 1.1
98 :     sub process_cds {
99 : overbeek 1.4 my($contigID,$cdsP,$prefix,$idNp,$contigs,$fh_tbl,$fh_fasta,$fh_ass) = @_;
100 : efrank 1.1 my($loc,@aliases,$func,$trans,$id,$precise);
101 :    
102 :     ($loc,$precise) = &get_loc($contigID,$cdsP);
103 :     @aliases = &get_aliases($cdsP);
104 :     $func = &get_func($cdsP);
105 :     $trans = &get_trans($cdsP);
106 :    
107 :     if ($loc || $trans)
108 :     {
109 :     my $dna = &FIG::extract_seq($contigs,$loc);
110 :     if ($precise)
111 :     {
112 :     my $prot = &FIG::translate($dna,undef,1);
113 :     if ($trans && (uc $prot ne uc $trans))
114 :     {
115 : overbeek 1.4 # print STDERR "BAD TRANSLATION: $contigID $$cdsP\n";
116 :     # print STDERR &Dumper($trans,$prot,$dna),"\n";
117 : efrank 1.1 }
118 :     }
119 :    
120 : overbeek 1.4 $id = $prefix . "$$idNp";
121 :     $$idNp++;
122 : efrank 1.1 print $fh_tbl "$id\t$loc\t",join("\t",@aliases),"\n";
123 : overbeek 1.4
124 : efrank 1.1 if ($func)
125 :     {
126 :     print $fh_ass "$id\t$func\n";
127 :     }
128 :     if ($trans)
129 :     {
130 :     &FIG::display_id_and_seq($id,\$trans,$fh_fasta);
131 :     }
132 :     }
133 :     else
134 :     {
135 :     print &Dumper($cdsP); die "aborted";
136 :     }
137 :     }
138 :    
139 :     sub get_loc {
140 :     my($contigID,$cdsP) = @_;
141 :     my($beg,$end,$loc,$locS,$iter,$n,$args,@pieces);
142 :     my($precise);
143 :    
144 : efrank 1.2 if ($$cdsP =~ /^(([^\n]+)((\n\s{21,23}[^\/ ][^\n]+)+)?)/s)
145 : efrank 1.1 {
146 :     $locS = $1;
147 :     $locS =~ s/\s//g;
148 :     $precise = ($locS !~ /[<>]/);
149 :    
150 :     @pieces = ();
151 :     $n = 0;
152 :     $iter = 500;
153 :     while (($locS !~ /^\[\d+\]$/) && $iter)
154 :     {
155 :     $iter--;
156 :     if ($locS =~ s/[<>]?(\d+)\.\.[<>]?(\d+)/\[$n\]/)
157 :     {
158 :     push(@pieces,["loc",$1,$2]);
159 :     $n++;
160 :     }
161 :    
162 :     if ($locS =~ s/([,\(])(\d+)([,\)])/$1\[$n\]$3/)
163 :     {
164 :     push(@pieces,["loc",$2,$2]);
165 :     $n++;
166 :     }
167 :    
168 :     if ($locS =~ s/join\((\[\d+\](,\[\d+\])*)\)/\[$n\]/)
169 :     {
170 :     $args = $1;
171 :     push(@pieces,["join",map { $_ =~ /^\[(\d+)\]$/; $1 } split(/,/,$args)]);
172 :     $n++;
173 :     }
174 :    
175 :     if ($locS =~ s/complement\((\[\d+\](,\[\d+\])*)\)/\[$n\]/g)
176 :     {
177 :     $args = $1;
178 :     push(@pieces,["complement",map { $_ =~ /^\[(\d+)\]$/; $1 } split(/,/,$args)]);
179 :     $n++;
180 :     }
181 :     }
182 :     if ($locS =~ /^\[(\d+)\]$/)
183 :     {
184 :     $loc = &conv(\@pieces,$1,$contigID);
185 :     }
186 :     else
187 :     {
188 :     print STDERR &Dumper(["could not parse $locS $iter",$cdsP]);
189 :     die "aborted";
190 :     }
191 :    
192 :     my @locs = split(/,/,$loc);
193 :     &trim_stop(\@locs);
194 :     $loc = join(",",@locs);
195 :     }
196 :     else
197 :     {
198 :     print STDERR &Dumper($cdsP); die "could not parse location";
199 :     die "aborted";
200 :     }
201 :     return ($loc,$precise);
202 :     }
203 :    
204 :     sub trim_stop {
205 :     my($locs) = @_;
206 :     my($to_trim,$n);
207 :    
208 :     $to_trim = 3;
209 :     while ($to_trim && (@$locs > 0))
210 :     {
211 :     $n = @$locs - 1;
212 :     if ($locs->[$n] =~ /^(\S+)_(\d+)_(\d+)$/)
213 :     {
214 :     if ($2 <= ($3-$to_trim))
215 :     {
216 :     $locs->[$n] = join("_",($1,$2,$3-$to_trim));
217 :     $to_trim = 0;
218 :     }
219 :     elsif ($2 >= ($3 + $to_trim))
220 :     {
221 :     $locs->[$n] = join("_",($1,$2,$3+$to_trim));
222 :     $to_trim = 0;
223 :     }
224 :     else
225 :     {
226 :     $to_trim -= abs($3-$2) + 1;
227 :     pop @$locs;
228 :     }
229 :     }
230 :     else
231 :     {
232 :     die "could not parse $locs->[$n]";
233 :     }
234 :     }
235 :     }
236 :    
237 :    
238 :     sub conv {
239 :     my($pieces,$n,$contigID) = @_;
240 :    
241 :     if ($pieces->[$n]->[0] eq "loc")
242 :     {
243 :     return join("_",$contigID,@{$pieces->[$n]}[1..2]);
244 :     }
245 :     elsif ($pieces->[$n]->[0] eq "join")
246 :     {
247 :     return join(",",map { &conv($pieces,$_,$contigID) } @{$pieces->[$n]}[1..$#{$pieces->[$n]}]);
248 :     }
249 :     elsif ($pieces->[$n]->[0] eq "complement")
250 :     {
251 :     return join(",",&complement(join(",",map { &conv($pieces,$_,$contigID) } @{$pieces->[$n]}[1..$#{$pieces->[$n]}])));;
252 :     }
253 :     }
254 :    
255 :     sub complement {
256 :     my($loc) = @_;
257 :     my(@locs);
258 :    
259 :     @locs = reverse split(/,/,$loc);
260 :     foreach $loc (@locs)
261 :     {
262 :     $loc =~ /^(\S+)_(\d+)_(\d+)$/;
263 :     $loc = join("_",($1,$3,$2));
264 :     }
265 :     return join(",",@locs);
266 :     }
267 :    
268 :     sub get_aliases {
269 :     my($cdsP) = @_;
270 :     my($db_ref);
271 :    
272 :     my @aliases = ();
273 : overbeek 1.4 while ($$cdsP =~ /\/(protein_id|gene|locus_tag)=\"([^\"]+)\"/sg)
274 : efrank 1.1 {
275 :     push(@aliases,$2);
276 :     }
277 :    
278 : overbeek 1.4 while ($$cdsP =~ /\/db_xref=\"([^\"]+)\"/sg)
279 : efrank 1.1 {
280 :     $db_ref = $1;
281 :     $db_ref =~ s/^GI:/gi\|/;
282 : overbeek 1.4 $db_ref =~ s/SWISS-PROT:/sp\|/;
283 : efrank 1.1 push(@aliases,$db_ref);
284 :     }
285 :    
286 :     return @aliases;
287 :     }
288 :    
289 :     sub get_trans {
290 :     my($cdsP) = @_;
291 :     my $tran = "";
292 :    
293 :     if ($$cdsP =~ /\/translation=\"([^"]*)\"/s)
294 :     {
295 :     $tran = $1;
296 :     $tran =~ s/\s//gs;
297 :     }
298 :     return $tran;
299 :     }
300 :    
301 :     sub get_func {
302 :     my($cdsP) = @_;
303 :     my $func = "";
304 :    
305 : overbeek 1.4 if (($$cdsP =~ /\/function=\"([^"]*)\"/s) && ($func = $1) && &ok_func($func))
306 :     {
307 :     $func =~ s/\s+/ /gs;
308 :     }
309 :     elsif ((($$cdsP =~ /\/product=\"([^"]*)\"/s) && ($func = $1)) && &ok_func($func) &&
310 :     (($func =~ / /) || ($$cdsP !~ /\/note/)))
311 :     {
312 :     $func =~ s/\s+/ /gs;
313 :     }
314 :     elsif (($$cdsP =~ /\/note=\"([^"]*)\"/s) && ($func = $1) && &ok_func($func))
315 : efrank 1.2 {
316 :     $func =~ s/\s+/ /gs;
317 :     }
318 : overbeek 1.4 $func = &fixup_func($func);
319 :     return $func;
320 :     }
321 :    
322 :     sub fixup_func {
323 :     my($func) = @_;
324 :    
325 :     $func =~ s/^COG\d+:\s+//i;
326 :     $func =~ s/^ORFID:\S+\s+//i;
327 : efrank 1.2 return $func;
328 :     }
329 :    
330 : overbeek 1.4 sub ok_func {
331 :     my($func) = @_;
332 :    
333 :     return (
334 :     ($func !~ /^[a-zA-Z]{1,3}\d+[a-zA-Z]{0,3}(\.\d+([a-zA-Z])?)?$/) &&
335 :     ($func !~ /^\d+$/) &&
336 :     ($func !~ /ensangp/i) &&
337 :     ($func !~ /agr_/i) &&
338 :     ($func !~ /^SC.*:/) &&
339 :     ($func !~ /^RIKEN/) &&
340 :     ($func !~ /\slen: \d+/) &&
341 :     ($func !~ /^similar to/i) &&
342 :     ($func !~ /^CPX/i) &&
343 :     ($func !~ /^\S{3,4}( or \S+ protein)?$/i) &&
344 :     ($func !~ /^putative$/) &&
345 :     ($func !~ /\scomes from/i) &&
346 :     ($func !~ /^unknown( function)?$/i) &&
347 :     ($func !~ /^hypothetical( (function|protein))?$/i) &&
348 :     ($func !~ /^orf /i) &&
349 :     ($func !~ /^ORF\s?\d+[lr]?$/i) &&
350 :     ($func !~ /^[a-z]{1,3}\s?\d+$/i)
351 :     )
352 :     }
353 :    
354 : efrank 1.2 sub process_rna {
355 : overbeek 1.4 my($contigID,$cdsP,$prefix,$idNr,$contigs,$fh_tbl,$fh_dna) = @_;
356 : efrank 1.2 my($loc,@aliases,$func,$trans,$id,$precise);
357 :    
358 :     ($loc,$precise) = &get_loc($contigID,$cdsP);
359 :     $func = &get_rna_func($cdsP);
360 :    
361 :     if ($loc)
362 :     {
363 :     my $dna = &FIG::extract_seq($contigs,$loc);
364 : overbeek 1.4 $id = $prefix . "$$idNr";
365 :     $$idNr++;
366 : efrank 1.2 if (! $func )
367 :     {
368 :     warn "could not get func $$cdsP\n";
369 :     }
370 :     else
371 :     {
372 :     print $fh_tbl "$id\t$loc\t$func\n";
373 :     &FIG::display_id_and_seq($id,\$dna,$fh_dna);
374 :     }
375 :     }
376 :     else
377 :     {
378 :     print &Dumper($cdsP); die "aborted";
379 :     }
380 :     }
381 :    
382 :     sub get_rna_func {
383 :     my($cdsP) = @_;
384 :     my $func = "";
385 :    
386 :     if ($$cdsP =~ /\/product=\"([^"]*)\"/s)
387 :     {
388 :     $func = $1;
389 :     $func =~ s/\s+/ /gs;
390 :     }
391 :     elsif ($$cdsP =~ /\/gene=\"([^"]*)\"/s)
392 :     {
393 :     $func = $1;
394 :     $func =~ s/\s+/ /gs;
395 :     }
396 :     elsif ($$cdsP =~ /\/note=\"([^"]*)\"/s)
397 : efrank 1.1 {
398 :     $func = $1;
399 :     $func =~ s/\s+/ /gs;
400 :     }
401 :     return $func;
402 :     }
403 :    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3