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

Annotation of /FigKernelScripts/parse_genbank.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.2 - (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 :     $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 :     &verify_dir("$dir/Features/peg");
14 :     &verify_dir("$dir/Features/rna");
15 :     open(CONTIGS,">$dir/contigs") "" die "could not open $dir/contigs";
16 :     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 :     open(ASSIGNMENTS,">$dir/assigned_functions") "" die "could not open $dir/assigned_functions";
19 :     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 : efrank 1.1 $/ = "\nLOCUS";
23 :     while ($contig = <STDIN>)
24 :     {
25 :     if ($contig =~ /\nACCESSION\s+(\S+)/s)
26 :     {
27 :     $id = $1;
28 :     if ($contig =~ /ORIGIN(.*)(\/\/|LOCUS)/s)
29 :     {
30 :     $seq = $1;
31 :     $seq =~ s/\s//gs;
32 :     $seq =~ s/\d//g;
33 :     undef $contigs;
34 :     $contigs->{$id} = $seq;
35 :     &FIG::display_id_and_seq($id,\$seq,\*CONTIGS);
36 :     }
37 :     else
38 :     {
39 :     print STDERR &Dumper($contig);
40 :     die "could not find sequence";
41 :     }
42 :    
43 :     while ($contig =~ /\n\s{4,6}CDS\s+([^\n]+(\n\s{20,22}\S[^\n]*)+)/gs)
44 :     {
45 :     $cds = $1;
46 : efrank 1.2 &process_cds($id,\$cds,$prefix,\$idN,$contigs,\*TBLPEG,\*FASTAPEG,\*ASSIGNMENTS);
47 :     }
48 :    
49 :     while ($contig =~ /\n\s{3,6}[tr]RNA\s+([^\n]+(\n\s{20,22}\S[^\n]*)+)/gs)
50 :     {
51 :     $rna = $1;
52 :     &process_rna($id,\$rna,$prefixR,\$idNr,$contigs,\*TBLRNA,\*FASTARNA);
53 : efrank 1.1 }
54 :     }
55 :     }
56 : efrank 1.2 close(CONTIGS);
57 :     close(TBLPEG);
58 :     close(FASTPEG);
59 :     close(ASSIGNMENTS);
60 :     close(TBLRNA);
61 :     close(FASTARNA);
62 :    
63 :     if (! -s "$dir/contigs") { unlink("$dir/contigs"); print STDERR "no contigs in $dir\n"; }
64 :     if (! -s "$dir/assigned_functions") { unlink("$dir/assigned_functions"); print STDERR "no assigned_functions in $dir }
65 :     if (! -s "$dir/Features/peg/tbl") { system "rm -rf $dir/Features/peg"; print STDERR "no PEGs in $dir }
66 :     if (! -s "$dir/Features/rna/tbl") { system "rm -rf $dir/Features/rna"; print STDERR "no RNAs in $dir }
67 :    
68 : efrank 1.1
69 :     sub process_cds {
70 :     my($contigID,$cdsP,$prefix,$idNP,$contigs,$fh_tbl,$fh_fasta,$fh_ass,$fh_dna) = @_;
71 :     my($loc,@aliases,$func,$trans,$id,$precise);
72 :    
73 :     ($loc,$precise) = &get_loc($contigID,$cdsP);
74 :     @aliases = &get_aliases($cdsP);
75 :     $func = &get_func($cdsP);
76 :     $trans = &get_trans($cdsP);
77 :    
78 :     if ($loc || $trans)
79 :     {
80 :     my $dna = &FIG::extract_seq($contigs,$loc);
81 :     if ($precise)
82 :     {
83 :     my $prot = &FIG::translate($dna,undef,1);
84 :     if ($trans && (uc $prot ne uc $trans))
85 :     {
86 :     print STDERR "BAD TRANSLATION: $contigID $$cdsP\n";
87 :     print STDERR &Dumper($trans,$prot,$dna),"\n";
88 :     }
89 :     }
90 :    
91 :     $id = $prefix . "$$idNP";
92 :     $$idNP++;
93 :     print $fh_tbl "$id\t$loc\t",join("\t",@aliases),"\n";
94 :     &FIG::display_id_and_seq($id,\$dna,$fh_dna);
95 :     if ($func)
96 :     {
97 :     print $fh_ass "$id\t$func\n";
98 :     }
99 :     if ($trans)
100 :     {
101 :     &FIG::display_id_and_seq($id,\$trans,$fh_fasta);
102 :     }
103 :     }
104 :     else
105 :     {
106 :     print &Dumper($cdsP); die "aborted";
107 :     }
108 :     }
109 :    
110 :     sub get_loc {
111 :     my($contigID,$cdsP) = @_;
112 :     my($beg,$end,$loc,$locS,$iter,$n,$args,@pieces);
113 :     my($precise);
114 :    
115 : efrank 1.2 if ($$cdsP =~ /^(([^\n]+)((\n\s{21,23}[^\/ ][^\n]+)+)?)/s)
116 : efrank 1.1 {
117 :     @loc = ();
118 :     $locS = $1;
119 :     $locS =~ s/\s//g;
120 :     $precise = ($locS !~ /[<>]/);
121 :    
122 :     @pieces = ();
123 :     $n = 0;
124 :     $iter = 500;
125 :     while (($locS !~ /^\[\d+\]$/) && $iter)
126 :     {
127 :     $iter--;
128 :     if ($locS =~ s/[<>]?(\d+)\.\.[<>]?(\d+)/\[$n\]/)
129 :     {
130 :     push(@pieces,["loc",$1,$2]);
131 :     $n++;
132 :     }
133 :    
134 :     if ($locS =~ s/([,\(])(\d+)([,\)])/$1\[$n\]$3/)
135 :     {
136 :     push(@pieces,["loc",$2,$2]);
137 :     $n++;
138 :     }
139 :    
140 :     if ($locS =~ s/join\((\[\d+\](,\[\d+\])*)\)/\[$n\]/)
141 :     {
142 :     $args = $1;
143 :     push(@pieces,["join",map { $_ =~ /^\[(\d+)\]$/; $1 } split(/,/,$args)]);
144 :     $n++;
145 :     }
146 :    
147 :     if ($locS =~ s/complement\((\[\d+\](,\[\d+\])*)\)/\[$n\]/g)
148 :     {
149 :     $args = $1;
150 :     push(@pieces,["complement",map { $_ =~ /^\[(\d+)\]$/; $1 } split(/,/,$args)]);
151 :     $n++;
152 :     }
153 :     }
154 :     if ($locS =~ /^\[(\d+)\]$/)
155 :     {
156 :     $loc = &conv(\@pieces,$1,$contigID);
157 :     }
158 :     else
159 :     {
160 :     print STDERR &Dumper(["could not parse $locS $iter",$cdsP]);
161 :     die "aborted";
162 :     }
163 :    
164 :     my @locs = split(/,/,$loc);
165 :     &trim_stop(\@locs);
166 :     $loc = join(",",@locs);
167 :     }
168 :     else
169 :     {
170 :     print STDERR &Dumper($cdsP); die "could not parse location";
171 :     die "aborted";
172 :     }
173 :     return ($loc,$precise);
174 :     }
175 :    
176 :     sub trim_stop {
177 :     my($locs) = @_;
178 :     my($to_trim,$n);
179 :    
180 :     $to_trim = 3;
181 :     while ($to_trim && (@$locs > 0))
182 :     {
183 :     $n = @$locs - 1;
184 :     if ($locs->[$n] =~ /^(\S+)_(\d+)_(\d+)$/)
185 :     {
186 :     if ($2 <= ($3-$to_trim))
187 :     {
188 :     $locs->[$n] = join("_",($1,$2,$3-$to_trim));
189 :     $to_trim = 0;
190 :     }
191 :     elsif ($2 >= ($3 + $to_trim))
192 :     {
193 :     $locs->[$n] = join("_",($1,$2,$3+$to_trim));
194 :     $to_trim = 0;
195 :     }
196 :     else
197 :     {
198 :     $to_trim -= abs($3-$2) + 1;
199 :     pop @$locs;
200 :     }
201 :     }
202 :     else
203 :     {
204 :     die "could not parse $locs->[$n]";
205 :     }
206 :     }
207 :     }
208 :    
209 :    
210 :     sub conv {
211 :     my($pieces,$n,$contigID) = @_;
212 :    
213 :     if ($pieces->[$n]->[0] eq "loc")
214 :     {
215 :     return join("_",$contigID,@{$pieces->[$n]}[1..2]);
216 :     }
217 :     elsif ($pieces->[$n]->[0] eq "join")
218 :     {
219 :     return join(",",map { &conv($pieces,$_,$contigID) } @{$pieces->[$n]}[1..$#{$pieces->[$n]}]);
220 :     }
221 :     elsif ($pieces->[$n]->[0] eq "complement")
222 :     {
223 :     return join(",",&complement(join(",",map { &conv($pieces,$_,$contigID) } @{$pieces->[$n]}[1..$#{$pieces->[$n]}])));;
224 :     }
225 :     }
226 :    
227 :     sub complement {
228 :     my($loc) = @_;
229 :     my(@locs);
230 :    
231 :     @locs = reverse split(/,/,$loc);
232 :     foreach $loc (@locs)
233 :     {
234 :     $loc =~ /^(\S+)_(\d+)_(\d+)$/;
235 :     $loc = join("_",($1,$3,$2));
236 :     }
237 :     return join(",",@locs);
238 :     }
239 :    
240 :     sub get_aliases {
241 :     my($cdsP) = @_;
242 :     my($db_ref);
243 :    
244 :     my @aliases = ();
245 :     if ($$cdsP =~ /\/(gene|locus_tag)=\"([^\"]+)\"/s)
246 :     {
247 :     push(@aliases,$2);
248 :     }
249 :    
250 :     if ($$cdsP =~ /\/db_xref=\"([^\"]+)\"/s)
251 :     {
252 :     $db_ref = $1;
253 :     $db_ref =~ s/^GI:/gi\|/;
254 :     push(@aliases,$db_ref);
255 :     }
256 :    
257 :     return @aliases;
258 :     }
259 :    
260 :     sub get_trans {
261 :     my($cdsP) = @_;
262 :     my $tran = "";
263 :    
264 :     if ($$cdsP =~ /\/translation=\"([^"]*)\"/s)
265 :     {
266 :     $tran = $1;
267 :     $tran =~ s/\s//gs;
268 :     }
269 :     return $tran;
270 :     }
271 :    
272 :     sub get_func {
273 :     my($cdsP) = @_;
274 :     my $func = "";
275 :    
276 :     if ($$cdsP =~ /\/product=\"([^"]*)\"/s)
277 : efrank 1.2 {
278 :     $func = $1;
279 :     $func =~ s/\s+/ /gs;
280 :     }
281 :     return $func;
282 :     }
283 :    
284 :     sub process_rna {
285 :     my($contigID,$cdsP,$prefix,$idNP,$contigs,$fh_tbl,$fh_dna) = @_;
286 :     my($loc,@aliases,$func,$trans,$id,$precise);
287 :    
288 :     ($loc,$precise) = &get_loc($contigID,$cdsP);
289 :     $func = &get_rna_func($cdsP);
290 :    
291 :     if ($loc)
292 :     {
293 :     my $dna = &FIG::extract_seq($contigs,$loc);
294 :     $id = $prefix . "$$idNP";
295 :     $$idNP++;
296 :     if (! $func )
297 :     {
298 :     warn "could not get func $$cdsP\n";
299 :     }
300 :     else
301 :     {
302 :     print $fh_tbl "$id\t$loc\t$func\n";
303 :     &FIG::display_id_and_seq($id,\$dna,$fh_dna);
304 :     }
305 :     }
306 :     else
307 :     {
308 :     print &Dumper($cdsP); die "aborted";
309 :     }
310 :     }
311 :    
312 :     sub get_rna_func {
313 :     my($cdsP) = @_;
314 :     my $func = "";
315 :    
316 :     if ($$cdsP =~ /\/product=\"([^"]*)\"/s)
317 :     {
318 :     $func = $1;
319 :     $func =~ s/\s+/ /gs;
320 :     }
321 :     elsif ($$cdsP =~ /\/gene=\"([^"]*)\"/s)
322 :     {
323 :     $func = $1;
324 :     $func =~ s/\s+/ /gs;
325 :     }
326 :     elsif ($$cdsP =~ /\/note=\"([^"]*)\"/s)
327 : efrank 1.1 {
328 :     $func = $1;
329 :     $func =~ s/\s+/ /gs;
330 :     }
331 :     return $func;
332 :     }
333 :    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3