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

Annotation of /FigKernelScripts/run_glimmer3.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 # -*- perl -*-
2 :    
3 :     $ENV{PATH} = "$ENV{HOME}/Glimmer-3.01/bin:$ENV{PATH}";
4 :    
5 :     $0 =~ m/([\/]+)$/;
6 :     $usage = "$1 Org-ID contigs [-train=training_tbl[,training_contigs]] [-skip]";
7 :    
8 :     $train = $training_tbl = $training_contigs = $skip = "";
9 :    
10 :     $trouble = 0;
11 :     if (@ARGV >= 2)
12 :     {
13 :     if ($ARGV[0] && (!-f $ARGV[0])) {
14 :     $org = $ARGV[0];
15 :     print STDERR "Org-ID is $org\n";
16 :     } else {
17 :     $trouble = 1;
18 :     warn "Bad Org-ID $ARGV[0]";
19 :     }
20 :    
21 :     if ($ARGV[1] && (-s $ARGV[1])) {
22 :     $contigs = $ARGV[1];
23 :     print STDERR "Contigs file is $contigs\n";
24 :     } else {
25 :     $trouble = 1;
26 :     warn "Missing contigs file $ARGV[1]";
27 :     }
28 :     }
29 :     else
30 :     {
31 :     die "usage: $usage";
32 :     }
33 :    
34 :     $i = 2;
35 :     while ($i < @ARGV)
36 :     {
37 :     if ($ARGV[$i] =~ m/-train=(\S+)/)
38 :     {
39 :     ++$i;
40 :     $train = $1;
41 :     if ($train =~ m/^([^,]+)/)
42 :     {
43 :     $training_tbl = $1;
44 :     if (-s $training_tbl)
45 :     {
46 :     print STDERR "Training set will be taken from $training_tbl\n";
47 :     }
48 :     else
49 :     {
50 :     $trouble = 1;
51 :     warn "Training tbl $training_tbl does not exist";
52 :     }
53 :    
54 :     if ($train =~ m/,([^,]+)$/)
55 :     {
56 :     $training_contigs = $1;
57 :     if (-s $training_contigs)
58 :     {
59 :     print STDERR "Training set will be extracted from $training_contigs\n";
60 :     }
61 :     else
62 :     {
63 :     $trouble = 1;
64 :     warn "Training contigs $training_contigs do not exist";
65 :     }
66 :     }
67 :     }
68 :     else
69 :     {
70 :     $trouble = 1;
71 :     warn "Training argument $train is invalid";
72 :     }
73 :     }
74 :     elsif ($ARGV[$i] =~ m/-skip/)
75 :     {
76 :     $skip = $ARGV[$i];
77 :     print STDERR "Skip option is set\n";
78 :     ++$i;
79 :     }
80 :     else
81 :     {
82 :     die "Unrecognized argument $ARGV[$i]";
83 :     }
84 :     }
85 :     die "\n\tusage: $usage\n\n" if $trouble;
86 :    
87 :     if (not $training_contigs)
88 :     {
89 :     $training_contigs = $contigs;
90 :     }
91 :    
92 :     $tmp_contig = "$FIG_Config::temp/tmp$$.contig";
93 :     $tmp_coord = "$FIG_Config::temp/tmp$$.coord";
94 :     $tmp_train = "$FIG_Config::temp/tmp$$.train";
95 :     $tmp_model = "$FIG_Config::temp/tmp$$.model";
96 :     $tmp_predict = "$FIG_Config::temp/tmp$$.predict";
97 :    
98 :     $orf_num = 0;
99 :     if ($train)
100 :     {
101 :     @tbl = `cat $training_tbl`;
102 :     foreach $entry (@tbl)
103 :     {
104 :     ($fid, $locus) = split /\t/, $entry;
105 :     $fid =~ m/\.(\d+)$/;
106 :     $orf_num = ($1 > $orf_num) ? $1 : $orf_num;
107 :    
108 :     $locus =~ m/^(\S+)_(\d+)_(\d+)$/;
109 :     ($contig_id, $beg, $end) = ($1, $2, $3);
110 :    
111 :     if ($skip) { $skip{$contig_id} = 1; }
112 :    
113 :     if (not defined($tbl{$contig_id})) { $tbl{$contig_id} = []; }
114 :     $x = $tbl{$contig_id};
115 :     push @$x, [$beg, $end];
116 :     }
117 :    
118 :     $n = 0;
119 :     open(CONTIGS, "<$training_contigs")
120 :     || die "Could not read-open $training_contigs";
121 :     open(TRAIN, ">$tmp_train")
122 :     || die "Could not write-open $tmp_train";
123 :    
124 :     while (($contig_id, $seqP) = &read_fasta_record(\*CONTIGS))
125 :     {
126 :     $len_of{$contig_id} = length($$seqP);
127 :    
128 :     $x = $tbl{$contig_id};
129 :     foreach $y (@$x)
130 :     {
131 :     ($beg, $end) = @$y;
132 :    
133 :     if ($beg < $end)
134 :     {
135 :     $dna = substr($$seqP,$beg-1,($end+1-$beg));
136 :     }
137 :     else
138 :     {
139 :     $dna = substr($$seqP,$end-1,($beg+1-$end));
140 :     $dna = $ { &rev_comp(\$dna) };
141 :     }
142 :     $dna = lc($dna);
143 :    
144 :     ++$n;
145 :     $_ = "T$n";
146 :     print TRAIN $_ . (" " x (11-length($_))) . $dna . "\n";
147 :     }
148 :     }
149 :     close(CONTIGS) || die "Could not close $contigs";
150 :     }
151 :     else
152 :     {
153 :     print STDERR "\nFinding training ORFs using default procedure\n";
154 :     if (-s "$tmp_train") {
155 :     system("rm -f $tmp_train")
156 :     && die "Could not remove $tmp_train";
157 :     }
158 :    
159 :     open(CONTIGS, "<$training_contigs")
160 :     || die "could not read-open $training_contigs";
161 :    
162 :     while (($contig_id, $seqP) = &read_fasta_record(\*CONTIGS))
163 :     {
164 :     $len_of{$contig_id} = length($$seqP);
165 :    
166 :     open( TMP, ">$tmp_contig") || die "Could not write-open $tmp_contig";
167 :     &display_id_and_seq($contig_id, $seqP, \*TMP);
168 :     close(TMP) || die "Could not close $tmp_contig";
169 :    
170 :     print STDERR "\nScanning contig $contig_id\n";
171 :     system("long-orfs -l $tmp_contig - | perl -e \'while(<>) { if (/^Putative/) { \$x=1; next; } print if \$x; }\' > $tmp_coord")
172 :     && die "Could not extract training ORFs from contig $contig_id";
173 :     system("extract $tmp_contig $tmp_coord >> $tmp_train")
174 :     && die "Could not extract training sequences from contig $contig_id";
175 :     }
176 :     close(CONTIGS) || die "Could not close $contigs";
177 :    
178 :     if (-s "$tmp_coord") { unlink("$tmp_coord") || warn "Could not remove $tmp_coord"; }
179 :     if (-s "$tmp_contig") { unlink("$tmp_contig") || warn "Could not remove $tmp_contig"; }
180 :     }
181 :    
182 :     if (($_ = `grep -c "^>" $tmp_train`) && ($_ =~ m/^\s*(\d+)/))
183 :     {
184 :     print STDERR "\nExtracted $1 training sequences\n\n";
185 :     }
186 :     else
187 :     {
188 :     die "\nCould not extract any training sequences";
189 :     }
190 :    
191 :     print STDERR "Building interpolated context model in $tmp_model\n";
192 :    
193 :     if (-s "$tmp_model")
194 :     {
195 :     system("rm -f $tmp_model") && die "Could not remove $tmp_model";
196 :     }
197 :    
198 :     system("build-icm $tmp_model <$tmp_train")
199 :     && die "Failure at build-icm $tmp_model <$tmp_train";
200 :    
201 :     #open(TBL, ">tbl") || die "Could not write-open tbl";
202 :     open(CONTIGS, "<$contigs") || die "Could not read-open $contigs";
203 :     while (($contig_id, $seqP) = &read_fasta_record(\*CONTIGS))
204 :     {
205 :     if ($skip{$contig_id})
206 :     {
207 :     print STDERR "Skipping $contig_id\n" if $ENV{VERBOSE};
208 :     next;
209 :     }
210 :    
211 :     $len_of{$contig_id} = length($$seqP);
212 :    
213 :     open( TMP, ">$tmp_contig") || die "Could not write-open $tmp_contig";
214 :     &display_id_and_seq($contig_id, $seqP, \*TMP);
215 :     close(TMP) || die "Could not close $tmp_contig";
216 :    
217 :     print STDERR "\nFinding ORFs on contig $contig_id, len=", length($$seqP), "\n";
218 :     system("glimmer3 -l $tmp_contig $tmp_model $FIG_Config::temp/tmp$$")
219 :     && die "FAILED: glimmer3 -l $tmp_contig $tmp_model $FIG_Config::temp/tmp$$";
220 :    
221 :     $found = 0;
222 :     open(CALLS, "<$FIG_Config::temp/tmp$$.predict")
223 :     || die "Could not read-open $FIG_Config::temp/tmp$$.predict";
224 :     while (defined($call = <CALLS>))
225 :     {
226 :     next if ($call =~ m/^>/);
227 :    
228 :     if ($call =~ m/^\S+\s+(\d+)\s+(\d+)/)
229 :     {
230 :     ++$found;
231 :     ++$orf_num;
232 :     print STDOUT "fig|$org.peg.$orf_num\t$contig_id\_$1\_$2\t\n";
233 :     }
234 :     else
235 :     {
236 :     warn "Invalid call format: $call\n";
237 :     }
238 :     }
239 :    
240 :     print STDERR "glimmer3 found $found ORFs on contig $contig_id\n";
241 :     }
242 :     close(CONTIGS) || die "Could not close $contigs";
243 :    
244 :     if (-s "$tmp_contig") { unlink("$tmp_contig") || warn "Could not remove $tmp_contig"; }
245 :     if (-s "$tmp_train") { unlink("$tmp_train") || warn "Could not remove $tmp_train"; }
246 :     if (-s "$tmp_model") { unlink("$tmp_model") || warn "Could not remove $tmp_model"; }
247 :    
248 :     exit 0;
249 :    
250 :    
251 :    
252 :     sub read_fasta_record
253 :     {
254 :     my ($file_handle) = @_;
255 :     my ( $old_end_of_record, $fasta_record, @lines, $head, $sequence, $seq_id, $comment, @parsed_fasta_record );
256 :    
257 :     if (not defined($file_handle)) { $file_handle = \*STDIN; }
258 :    
259 :     $old_end_of_record = $/;
260 :     $/ = "\n>";
261 :    
262 :     if (defined($fasta_record = <$file_handle>))
263 :     {
264 :     chomp $fasta_record;
265 :     @lines = split( /\n/, $fasta_record );
266 :     $head = shift @lines;
267 :     $head =~ s/^>?//;
268 :     $head =~ m/^(\S+)/;
269 :     $seq_id = $1;
270 :    
271 :     if ($head =~ m/^\S+\s+(.*)$/) { $comment = $1; } else { $comment = ""; }
272 :    
273 :     $sequence = join( "", @lines );
274 :    
275 :     @parsed_fasta_record = ( $seq_id, \$sequence, $comment );
276 :     }
277 :     else
278 :     {
279 :     @parsed_fasta_record = ();
280 :     }
281 :    
282 :     $/ = $old_end_of_record;
283 :    
284 :     return @parsed_fasta_record;
285 :     }
286 :    
287 :     sub display_id_and_seq
288 :     {
289 :     my( $id, $seq, $fh ) = @_;
290 :    
291 :     if (! defined($fh) ) { $fh = \*STDOUT; }
292 :    
293 :     print $fh ">$id\n";
294 :     &display_seq($seq, $fh);
295 :     }
296 :    
297 :     sub display_seq
298 :     {
299 :     my ( $seq, $fh ) = @_;
300 :     my ( $i, $n, $ln );
301 :    
302 :     if (! defined($fh) ) { $fh = \*STDOUT; }
303 :    
304 :     $n = length($$seq);
305 :     # confess "zero-length sequence ???" if ( (! defined($n)) || ($n == 0) );
306 :     for ($i=0; ($i < $n); $i += 60)
307 :     {
308 :     if (($i + 60) <= $n)
309 :     {
310 :     $ln = substr($$seq,$i,60);
311 :     }
312 :     else
313 :     {
314 :     $ln = substr($$seq,$i,($n-$i));
315 :     }
316 :     print $fh "$ln\n";
317 :     }
318 :     }
319 :    
320 :     sub rev_comp
321 :     {
322 :     my( $seqP ) = @_;
323 :     my( $rev );
324 :    
325 :     $rev = reverse( $$seqP );
326 :     $rev =~ tr/a-z/A-Z/;
327 :     $rev =~ tr/ACGTUMRWSYKBDHV/TGCAAKYWSRMVHDB/;
328 :     return \$rev;
329 :     }
330 :    
331 :     sub min { my ($x, $y) = @_; return (($x < $y) ? $x : $y); }
332 :    
333 :     sub max { my ($x, $y) = @_; return (($x > $y) ? $x : $y); }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3