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

Annotation of /FigKernelScripts/extract_metagene_contigs.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : arodri7 1.1 # -*- perl -*-
2 :     ########################################################################
3 :     # Copyright (c) 2003-2006 University of Chicago and Fellowship
4 :     # for Interpretations of Genomes. All Rights Reserved.
5 :     #
6 :     # This file is part of the SEED Toolkit.
7 :     #
8 :     # The SEED Toolkit is free software. You can redistribute
9 :     # it and/or modify it under the terms of the SEED Toolkit
10 :     # Public License.
11 :     #
12 :     # You should have received a copy of the SEED Toolkit Public License
13 :     # along with this program; if not write to the University of Chicago
14 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
15 :     # Genomes at veronika@thefig.info or download a copy from
16 :     # http://www.theseed.org/LICENSE.TXT.
17 :     ########################################################################
18 :    
19 :     use FIG;
20 :     use File::Basename;
21 :     use strict;
22 :     use Pod::Text;
23 :    
24 :     use strict;
25 :     use Pod::Text;
26 :     # use File::Basename; $this_tool_name = basename($0);
27 :     if ((@ARGV == 1) && ($ARGV[0] =~ m/-help/)) {
28 :     pod2text($0); exit(0);
29 :     }
30 :    
31 :     =pod
32 :    
33 :     =over 5
34 :    
35 : arodri7 1.4 =item Usage: extract_metagene_contigs [-help] [-output-directory=output directory] fasta_in fasta_out
36 : arodri7 1.1
37 :     =item Function: Given a contig fasta input file, metagene tool is run to find any contig which may have multiple gene calls.
38 :     If a contig with multiple gene calls is found, multiple new sequences are created in the output fasta file,
39 :     otherwise the original contig is written to the output fasta file.
40 :     Three output files will be generated:
41 :     - The raw metagene output
42 :     - The parsed metagene output of those sequences which have multiple genes, and their start and end locations.
43 :     - The new fasta file generated which include the subsequences generated from the metagene output.
44 :    
45 :     -output-directory Directory address where the output files should be written.
46 : arodri7 1.4 -peg Default is false. If true [peg=T], then the new ids assigned will be peg ids
47 :     -genome genome id for peg sequences
48 :     -AA Default is false. If true [AA=T], the output sequences will be given as amino acid sequences.
49 :     -exec Default is metagene. Metagene version to run. Options are -exec=mga or -exec=metagene
50 :     -featuremap File where you wish the new feature to contig map to be written to. Default is "metagene.map"
51 :     -startpeg when creating peg feature id, this would be the number to start at. Default is 1.
52 :     -called_by if present, it will write a file feature.called_by which states that metagene called the gene
53 : arodri7 1.1
54 :     =back
55 :    
56 :     =cut
57 :    
58 :     my $outD = "";
59 :     my $trouble = 0;
60 :     my $input_eol_marker = $/;
61 :     my $fasta_in = "";
62 :     my $fasta_out = "";
63 :     my $first_line = 1;
64 : arodri7 1.4 my $featuremap;
65 :     my $startpeg = 1;
66 :     my $featurecalled;
67 : arodri7 1.1
68 :     my $badchar = 0;
69 :     my $max_bad = 500;
70 :     my $width = 50;
71 : arodri7 1.4 my ($peg_flag,$genome_id,$amino,$tool);
72 :     my $exec = "$FIG_Config::metagene_bin/metagene";
73 :     my $tool="metagene";
74 : arodri7 1.1
75 :     while ($ARGV[0] =~ m/^-/) {
76 :     if ($ARGV[0] =~ m/-output-directory=(\S+)/) {
77 :     $outD = $1;
78 :     }
79 : arodri7 1.4 elsif ($ARGV[0] =~ m/-peg=(\S+)/) {
80 :     $peg_flag=1 if ($1 eq "T");
81 :     }
82 :     elsif ($ARGV[0] =~ m/-genome=(\S+)/) {
83 :     $genome_id=$1;
84 :     }
85 :     elsif ($ARGV[0] =~ m/-AA=(\S+)/) {
86 :     $amino=1 if ($1 eq "T");
87 :     }
88 :     elsif ($ARGV[0] =~ m/-featuremap=(\S+)/) {
89 :     $featuremap=$1;
90 :     }
91 :     elsif ($ARGV[0] =~ m/-called_by=(\S+)/) {
92 :     $featurecalled=$1;
93 :     }
94 :     elsif ($ARGV[0] =~ m/-startpeg=(\S+)/) {
95 :     $startpeg=$1;
96 :     }
97 :     elsif ($ARGV[0] =~ m/-exec=(\S+)/) {
98 :     if ($1 eq "mga") {
99 :     $tool="mga";
100 : gdpusch 1.12 if (-x "$FIG_Config::mga_bin/mga_exec") {
101 :     $exec="$FIG_Config::mga_bin/mga_exec";
102 :     }
103 :     elsif (-x "$FIG_Config::mga_bin/mga") {
104 :     $exec="$FIG_Config::mga_bin/mga";
105 :     }
106 :     else {
107 :     $trouble = 1;
108 :     }
109 : arodri7 1.4 }
110 :     elsif ($1 eq "metagene") {
111 : gdpusch 1.12 $tool="metagene";
112 : arodri7 1.4 $exec = "$FIG_Config::metagene_bin/$1";
113 :     }
114 :     else {
115 :     $trouble = 1;
116 : gdpusch 1.12 warn "Can't find MGA";
117 : arodri7 1.4 }
118 :     }
119 : arodri7 1.1 else {
120 :     $trouble = 1;
121 :     print STDERR "bad arg $ARGV[0]\n\n";
122 :     }
123 :    
124 :     shift;
125 :     }
126 :    
127 :     if (@ARGV == 2) {
128 :     $fasta_in = shift;
129 :     $fasta_out = shift;
130 :    
131 :     if (!-s $fasta_in) {
132 :     $trouble = 1;
133 :     warn "Input file $fasta_in does not exist\n";
134 :     }
135 :     }
136 :    
137 :     if ($trouble || @ARGV) {
138 :     warn qq(There were invalid arguments: ), join(" ", @ARGV), qq(\n\n);
139 :     pod2text($0);
140 :     die "aborted";
141 :     }
142 :    
143 :     my $fig = new FIG;
144 : arodri7 1.4 if (!$outD){
145 :     pod2text($0); exit(0);
146 :     }
147 : arodri7 1.1
148 :     # run metagene on the raw input file
149 : arodri7 1.4 my @metagene_args;
150 :     if ($tool eq "mga") {
151 :     push (@metagene_args, "-s");
152 :     }
153 :     push (@metagene_args, ($fasta_in, "> $outD/metagene.out"));
154 :    
155 :     $fig->run(join(" ", $exec, @metagene_args));
156 :    
157 : arodri7 1.1
158 : arodri7 1.4 if (!$featuremap)
159 :     {
160 :     $featuremap = "$outD/metagene_seqs.map";
161 :     }
162 :    
163 :     if ($featurecalled)
164 :     {
165 :     open (CALLED_BY, ">>$featurecalled");
166 :     }
167 :    
168 :     open (FH_MG, ">>$featuremap"); # write the sequences that need to be divided into subsequences in this file
169 : arodri7 1.1
170 :     # parse the metagene output file
171 :     my $metagene_hits = {};
172 : arodri7 1.4 if ($tool eq "metagene")
173 :     {
174 :     $/ = "\n\n";
175 :     }
176 : arodri7 1.1 open (FH, "$outD/metagene.out");
177 : arodri7 1.4 my $count=$startpeg;
178 : arodri7 1.3
179 : gdpusch 1.8 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
180 :     # ... Parse Metagene / MetageneAnnotator ouput ...
181 :     #-----------------------------------------------------------------------
182 :     my @mg_hits = ();
183 : arodri7 1.4 if ($tool eq "metagene")
184 :     {
185 :     while (my $block = <FH>){
186 : gdpusch 1.8 my ($contig_id) = ($block) =~ m/^\# (.*?)\n/;
187 : arodri7 1.4 my @lines = split (/\n/, $block);
188 :     if (@lines >= 4){ # write into a file the new sequences
189 :     foreach my $line (@lines){
190 :     next if ($line =~ /^\#/);
191 :    
192 : gdpusch 1.8 chomp $line;
193 :     my ($start, $end, $strand, undef, $trunc) = split (/\t/, $line);
194 :     push @mg_hits, [$contig_id, $start, $end, $strand, $trunc, $line];
195 :    
196 :     # $metagene_hits->{$contig_id}->{$new_fid}->{start} = $map_start;
197 :     # $metagene_hits->{$contig_id}->{$new_fid}->{stop} = $map_end;
198 :     # $metagene_hits->{$contig_id}->{$new_fid}->{strand} = $strand;
199 : arodri7 1.4 }
200 :     }
201 :     }
202 :     }
203 :     elsif ($tool eq "mga")
204 :     {
205 :     my $line = <FH>;
206 :     while ($line && $line =~ m/^\#\s+(\S+)/)
207 :     {
208 : gdpusch 1.8 my $contig_id = $1;
209 : arodri7 1.4 $line= <FH>; $line= <FH>;
210 :     $line = <FH>;
211 :     while ($line && $line !~ /^\#/)
212 :     {
213 : gdpusch 1.8 chomp $line;
214 : arodri7 1.10 my (undef, $start, $end, $strand, $frame, $trunc) = split (/\t/, $line);
215 :     push @mg_hits, [$contig_id, $start, $end, $strand, $frame, $trunc, $line];
216 : gdpusch 1.8
217 :     # $metagene_hits->{$contig_id}->{$new_fid}->{start} = $map_start;
218 :     # $metagene_hits->{$contig_id}->{$new_fid}->{stop} = $map_end;
219 :     # $metagene_hits->{$contig_id}->{$new_fid}->{strand} = $strand;
220 :    
221 : arodri7 1.4 $line = <FH>;
222 : arodri7 1.1 }
223 :     }
224 :     }
225 :     close (FH);
226 : arodri7 1.4 $/ = "\n" if ($tool eq "metagene");
227 : gdpusch 1.8
228 :    
229 :     foreach my $call (@mg_hits) {
230 : arodri7 1.10 my ($contig_id, $start, $end, $strand, $frame, $trunc, $line) = @$call;
231 : gdpusch 1.8
232 :     my ($map_start, $map_end, $sign);
233 :    
234 :     if ($strand eq "+") {
235 :     $map_start = $start;
236 :     $map_end = $end;
237 :     }
238 :     else {
239 :     $map_start = $end;
240 :     $map_end = $start;
241 :     }
242 :    
243 :     if ($tool eq qq(metagene)) {
244 :     if ($trunc ne qq(complete)) {
245 :     warn qq(Skipping METAGENE partial: $line\n);
246 :     next;
247 :     }
248 :     }
249 :     elsif ($tool eq qq(mga)) {
250 :     $sign = ($map_end <=> $map_start);
251 :     if ($trunc eq qq(01)) {
252 :     $map_start = $map_end - $sign * ( 3 * int( (1 + abs($map_end-$map_start)) / 3) - 1);
253 :     warn qq(Fixed truncated START:\t($start, $end) ==> ($map_start, $map_end):\t$line\n);
254 :     }
255 :     elsif ($trunc eq qq(10)) {
256 :     $map_end = $map_start + $sign * ( 3 * int( (1 + abs($map_end-$map_start)) / 3) - 1);
257 :     warn qq(Fixed truncated STOP:\t($start, $end) ==> ($map_start, $map_end):\t$line\n);
258 :     }
259 :     elsif ($trunc eq qq(00)) {
260 : arodri7 1.10 if ($strand eq "+")
261 :     {
262 :     $map_start = $map_start + $frame;
263 :     }
264 :     elsif ($strand eq "-")
265 :     {
266 :     $map_start = $map_start - $frame;
267 :     }
268 :    
269 :     $map_end = $map_start + $sign * ( 3 * int( (1 + abs($map_end-$map_start)) / 3) - 1);
270 :     #warn qq(Skipping MGA double-truncated partial: $line\n);
271 :     warn qq(Fixed MGA double-truncated partial: $line\n);
272 : arodri7 1.11 #next;
273 : gdpusch 1.8 }
274 :     }
275 :     else {
276 :     die qq(Unrecognized tool: $tool);
277 :     }
278 :    
279 :     my $new_fid;
280 :     if ($peg_flag) {
281 :     $new_fid = "fig|".$genome_id.".peg.".$count;
282 :     $count++;
283 :     }
284 :     else {
285 :     $new_fid = $contig_id . "_mg_" . $map_start . "_" . $map_end;
286 :     }
287 :    
288 : arodri7 1.9 $metagene_hits->{$contig_id}->{$new_fid}->{start} = $map_start;
289 :     $metagene_hits->{$contig_id}->{$new_fid}->{stop} = $map_end;
290 :     $metagene_hits->{$contig_id}->{$new_fid}->{strand} = $strand;
291 :    
292 : gdpusch 1.8 print FH_MG "$new_fid\t$contig_id"."_".$map_start."_". $map_end."\n";
293 :     print CALLED_BY "$new_fid\t$tool\n" if ($featurecalled);
294 :     }
295 : arodri7 1.1 close (FH_MG);
296 : arodri7 1.4 close (CALLED_BY) if ($featurecalled);
297 : arodri7 1.1
298 :    
299 : arodri7 1.5 if ($amino)
300 :     {
301 :     $fig->run("get_fasta_for_tbl_entries $fasta_in < $featuremap > $fasta_out");
302 :     }
303 :     else
304 :     {
305 : arodri7 1.9 # die qq(Disabled option --- aborting);
306 : gdpusch 1.8
307 : arodri7 1.5 my $fh_in = \*STDIN;
308 :     my $fh_out = \*STDOUT;
309 :    
310 : gdpusch 1.8 open($fh_in, "<$fasta_in") || die "Could not read-open $fasta_in";
311 :     open($fh_out, ">>$fasta_out") || die "Could not append-open $fasta_out";
312 : arodri7 1.5
313 :     while (my($head, $seqP) = &get_fasta_record($fh_in)) {
314 : gdpusch 1.8 # my $contig_len = length($$seqP);
315 :    
316 : arodri7 1.5 if (defined ($metagene_hits->{$head})){
317 : gdpusch 1.8 my @new_fids;
318 : arodri7 1.5 if ($peg_flag) {
319 : gdpusch 1.8 @new_fids = sort { FIG::by_fig_id($a,$b) } keys %{$metagene_hits->{$head}};
320 : arodri7 1.4 }
321 : arodri7 1.5 else {
322 : gdpusch 1.8 @new_fids = keys %{$metagene_hits->{$head}};
323 : arodri7 1.5 }
324 :    
325 : gdpusch 1.8 foreach my $new_fid ( @new_fids){
326 : arodri7 1.9 #if ($new_fid eq "fig|441534.3.peg.1805")
327 :     #{
328 :     # print "debug_time\n";
329 :     #}
330 : gdpusch 1.8 my $start = $metagene_hits->{$head}->{$new_fid}->{start};
331 :     my $stop = $metagene_hits->{$head}->{$new_fid}->{stop};
332 :     my $strand = $metagene_hits->{$head}->{$new_fid}->{strand};
333 : arodri7 1.5
334 : arodri7 1.9 my ($mg_contig);
335 :     if ($start > $stop)
336 :     {
337 :     $mg_contig = substr($$seqP, $stop-1, $start-($stop-1));
338 :     }
339 :     else
340 :     {
341 :     $mg_contig = substr($$seqP, $start-1, $stop-($start-1));
342 :     }
343 : gdpusch 1.8 if ($strand eq "-") {
344 :     $mg_contig = &FIG::rev_comp( \$mg_contig );
345 : arodri7 1.9 $mg_contig = $$mg_contig;
346 : gdpusch 1.8 }
347 :     &display_id_and_seq( $fig, $new_fid, \$mg_contig, $fh_out, $amino );
348 : arodri7 1.5 }
349 :     }
350 :     elsif (!$peg_flag)
351 :     {
352 :     &display_id_and_seq( $fig, $head, $seqP, $fh_out, $amino );
353 : arodri7 1.1 }
354 :     }
355 : arodri7 1.5 close ($fh_in);
356 :     close ($fh_out);
357 : arodri7 1.1 }
358 :    
359 :     sub get_fasta_record {
360 :     my ( $fh ) = @_;
361 :     my ( $old_eol, $entry, @lines, $head, $seq);
362 :    
363 :     if (not defined($fh)) { $fh = \*STDIN; }
364 :     $old_eol = $/;
365 :     $/ = "$input_eol_marker>";
366 :    
367 :     my @record = ();
368 :     if (defined($entry = <$fh>)) {
369 :     chomp $entry;
370 :     @lines = split( /$input_eol_marker/, $entry );
371 :     while (@lines and not $lines[0]) { shift @lines; }
372 :    
373 :     $head = shift @lines;
374 :     if ($first_line) {
375 :     $first_line = 0;
376 :    
377 :     if (not ($head =~ s/^\s*>//)) {
378 :     $trouble = 1;
379 :     warn $head;
380 :     die "ERROR: File does not appear to be in FASTA format\n";
381 :     }
382 :     }
383 :     else {
384 :     if ($head =~ s/^\s*>//) {
385 :     $trouble = 1;
386 :     warn $head;
387 :     die "Spurious beginning-of record mark found in record $.\n";
388 :     }
389 :     }
390 :    
391 :     foreach my $ln (@lines) {
392 :     $_ = $ln;
393 :     $ln =~ s/\s//g;
394 :    
395 :     print STDERR "$head: contains X's\n" if ($ln =~ s/x/n/ig);
396 :     print STDERR "$head: contains colons\n" if ($ln =~ s/://g);
397 :    
398 :     while ($ln =~ s/([^ACGTUMRWSYKBDHVN]+)/n/i) {
399 :     $trouble = 1;
400 :     $badchar++;
401 :     print STDERR ">$head:\tbad char $1 at ", pos($ln), " at line $.\n";
402 :     }
403 :     }
404 :    
405 :     $seq = join( "", @lines );
406 :     $seq =~ s/\cM//g;
407 :     $seq =~ tr/a-z/A-Z/;
408 :     @record = ($head, \$seq);
409 :     }
410 :    
411 :     $/ = $old_eol;
412 :     return @record;
413 :     }
414 :    
415 :    
416 :     sub display_id_and_seq {
417 : arodri7 1.4 my( $fig, $id, $seq, $fh, $amino ) = @_;
418 : arodri7 1.1 my ( $i, $n, $ln );
419 :    
420 :     if (! defined($fh) ) { $fh = \*STDOUT; }
421 :    
422 :     print $fh ">$id\n";
423 : arodri7 1.4 my $dna_seq = $$seq;
424 :     my $sequence;
425 :    
426 :     if ($amino) {
427 :     my $code_number = 11;
428 :     my $change_start_to_M=1;
429 :     #my $genetic_code_hash = $fig->genetic_code($code_number);
430 :     my $genetic_code_hash = $fig->standard_genetic_code();
431 :     $sequence = $fig->translate($dna_seq, $genetic_code_hash, $change_start_to_M);
432 :     }
433 :     else {
434 :     $sequence = $dna_seq;
435 :     }
436 : arodri7 1.1
437 : arodri7 1.4 $n = length($sequence);
438 : arodri7 1.1 # confess "zero-length sequence ???" if ( (! defined($n)) || ($n == 0) );
439 :     for ($i=0; ($i < $n); $i += $width) {
440 :     if (($i + $width) <= $n) {
441 : arodri7 1.4 $ln = substr($sequence,$i,$width);
442 : arodri7 1.1 }
443 :     else {
444 : arodri7 1.4 $ln = substr($sequence,$i,($n-$i));
445 : arodri7 1.1 }
446 :    
447 :     print $fh "$ln\n";
448 :     }
449 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3