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

Annotation of /FigKernelScripts/bielefeld_tarfile_to_seed.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 # -*- perl -*-
2 : olson 1.4 #
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 : overbeek 1.1
20 :     use FIG;
21 :     $fig = new FIG;
22 :    
23 :     use File::Path;
24 :    
25 : overbeek 1.2 $usage = "bielefeld_gff_to_seed.pl [-update] [-discard] input_tarfile output_dir";
26 : overbeek 1.1
27 : overbeek 1.2 $trouble = 0;
28 :     while ($ARGV[0] =~ m/^-/)
29 :     {
30 :     if ($ARGV[0] =~ m/-update/)
31 :     {
32 :     $update = shift @ARGV;
33 :     }
34 :     elsif ($ARGV[0] =~ m/-discard/)
35 :     {
36 :     $discard = shift @ARGV;
37 :     }
38 :     else
39 :     {
40 :     warn "\tInvalid argument $ARGV[0]\n";
41 :     $trouble = shift @ARGV;
42 :     }
43 :     }
44 :    
45 :     if ($trouble) { die "\nThere were invalid arguments\n\n\tusage = $usage\n\n"; }
46 : overbeek 1.1
47 :     (($tarfile = shift) && (-s $tarfile))
48 :     || die "Input tarfile $tarfile does not exist.\n\n\tusage: $usage\n\n";
49 : overbeek 1.3 print STDERR "Reading data from $tarfile\n";
50 : overbeek 1.1
51 :     ($output_dir = shift @ARGV) || die "no Output-Dir given\n\usage:$usage\n\n";
52 :    
53 :     if ($output_dir =~ m/(\d+\.\d+)$/)
54 :     {
55 :     $tax_id = $1;
56 :     $org = quotemeta($tax_id);
57 :     }
58 :     else
59 :     {
60 :     die "Output-dir $output_dir must end in a valid format Taxon-ID, e.g. \"83333.1\"";
61 :     }
62 :    
63 :     ($base_dir = `pwd`) || die "Could not get current working directory";
64 :     chomp $base_dir;
65 :     print STDERR "\nSetting base_dir = $base_dir\n";
66 :    
67 :     $dir = "$base_dir/$output_dir/RAW_DATA";
68 :     (-d $dir) || mkpath($dir, 1, 0777) || die "could not create $dir";
69 : overbeek 1.2 # chdir($dir) || die "Could not change working-directory to $dir";
70 : overbeek 1.1
71 : overbeek 1.3 if ($tarfile !~ m{^\/})
72 : overbeek 1.2 {
73 : overbeek 1.3 if (-s "$base_dir/$tarfile")
74 :     {
75 :     $tarfile = "$base_dir/$tarfile";
76 :     }
77 :     else
78 :     {
79 :     die "Could not locate $base_dir/$tarfile";
80 :     }
81 : overbeek 1.2 }
82 :     else
83 :     {
84 : overbeek 1.3 print STDERR "Reading data from $tarfile\n";
85 : overbeek 1.2 }
86 :    
87 :     &FIG::run("tar xpzvf $tarfile -C $dir");
88 : overbeek 1.1
89 :     chdir("$base_dir/$output_dir") || "Could not change working-directory to $output_dir";
90 : overbeek 1.3 &FIG::run("cat RAW_DATA/*\.fas | reformat_contigs > contigs");
91 : overbeek 1.1 open(TMP, "<contigs") || die "could not read-open contigs";
92 :     while (($id, $seqP) = &read_fasta_record(\*TMP)) { $len{$id} = length($$seqP); }
93 :     close(TMP) || die "could not close contigs";
94 :    
95 : overbeek 1.3 if ($num_contigs = (scalar keys %len))
96 :     {
97 :     print STDERR "Read $num_contigs contigs\n";
98 :     }
99 :     else
100 :     {
101 :     die "ERROR: no contigs found in contigs file\n";
102 :     }
103 :    
104 :    
105 : overbeek 1.1 $dir = "Features/peg";
106 :     (-d $dir) || mkpath($dir, 1, 0777) || die "could not create $dir";
107 :     open(CDS_TBL, ">$dir/tbl") || die "could not open $dir/tbl";
108 :    
109 :     $dir = "Features/rna";
110 :     (-d $dir) || mkpath($dir, 1, 0777) || die "could not create $dir";
111 :     open(RNA_TBL, ">$dir/tbl") || die "could not open $dir/tbl";
112 :    
113 :     $peg_num = 0;
114 :     if ($update && (-d ($dir = "$FIG_Config::organisms/$tax_id/Features/peg")))
115 :     {
116 :     open(TMP, "<$dir/tbl") || die "Could not open $dir/tbl";
117 :     while (defined($line = <TMP>))
118 :     {
119 :     if ($line =~ m/^fig\|$org\.peg\.(\d+)/) { $peg_num = &FIG::max($peg_num, $1); } else { die "PEG match failed for $line"; }
120 :     }
121 :     close(TMP) || die "Could not close $FIG_Config::organisms/$tax_id/Features/peg/tbl";
122 :     }
123 :    
124 :     $rna_num = 0;
125 :     if ($update && (-d ($dir = "$FIG_Config::organisms/$tax_id/Features/rna")))
126 :     {
127 :     open(TMP, "<$dir/tbl") || die "Could not open $dir/tbl";
128 :     while (defined($line = <TMP>))
129 :     {
130 :     if ($line =~ m/^fig\|\d+\.\d+\.rna\.(\d+)/) { $rna_num = &FIG::max($rna_num, $1); } else { die "RNA match failed for $line"; }
131 :     }
132 :     close(TMP) || die "Could not close $FIG_Config::organisms/$tax_id/Features/rna/tbl";
133 :     }
134 :    
135 :     opendir(RAW_DATA, "RAW_DATA") || die "Could not opendir $dir/RAW_DATA";
136 :     (@gff_files = grep /_genes.gff$/, readdir(RAW_DATA)) || die "Could not readdir GFF files";
137 :     foreach $file (@gff_files)
138 :     {
139 :     print STDERR "\nParsing RAW_DATA/$file ...\n";
140 :     open(GFF, "<RAW_DATA/$file") || die "Could not read-open RAW_DATA/$file";
141 :     while (defined($line = <GFF>))
142 :     {
143 :     chomp $line;
144 :    
145 :     next if ($line =~ m/^\#/);
146 :     next if ($line !~ m/\S/);
147 :    
148 :     ($contig, undef, $type, $left, $right, undef, $strand, undef, $comment) = split /\t/, $line;
149 :    
150 : overbeek 1.2 next if ($discard && ($comment =~ m/^discarded/));
151 :    
152 : overbeek 1.1 if ($left <= 0)
153 :     {
154 :     ++$bad;
155 :     while ($left <= 0) { $left += 3; } #...Hack
156 :     }
157 :    
158 :     if ($right > $len{$contig})
159 :     {
160 :     ++$bad;
161 :     while ($right > $len{$contig}) { $right -= 3; } #...Hack
162 :     }
163 :    
164 :     if ($strand eq '+')
165 :     {
166 :     ($beg, $end) = ($left, $right);
167 :     }
168 :     elsif ($strand eq '-')
169 :     {
170 :     ($beg, $end) = ($right, $left);
171 :     }
172 :     else
173 :     {
174 :     die "Invalid strand \'$strand\' in line $.: $line\n";
175 :     }
176 :    
177 :     if ($type =~ m/CDS/)
178 :     {
179 :     ++$peg_num;
180 :     print CDS_TBL "fig|$tax_id.peg.$peg_num\t$contig\_$beg\_$end\t\n";
181 :     }
182 :     elsif ($type =~ m/RNA/)
183 :     {
184 :     ++$rna_num;
185 :     $comment =~ m/^([^\;]+)/;
186 :     print RNA_TBL "fig|$tax_id.rna.$rna_num\t$contig\_$beg\_$end\t$1\n";
187 :     # die "$line\n";
188 :     }
189 :     else
190 :     {
191 :     print STDERR "Could not handle $line\n\n";
192 :     }
193 :     }
194 :     }
195 :     print STDERR "bad=$bad\n" if $bad;
196 :    
197 : overbeek 1.2 close(CDS_TBL) || die "Could not close peg tbl";
198 :     close(RNA_TBL) || die "Could not close rna tbl";
199 :    
200 : overbeek 1.1 $dir = "$base_dir/$output_dir";
201 :     if (-s "$dir/Features/rna/tbl")
202 :     {
203 :     print STDERR "\nGetting RNA FASTA sequences ...\n";
204 :     &FIG::run("get_dna $dir/contigs $dir/Features/rna/tbl > $dir/Features/rna/fasta");
205 :     }
206 :     else
207 :     {
208 :     warn "\'$output_dir/Features/rna/tbl\' is empty --- deleting Features/rna/\n";
209 :     # system("rm -R $output_dir/Features/rna") && die("could not delete empty Features/rna");
210 :     }
211 :    
212 :     if (-s "$dir/Features/peg/tbl")
213 :     {
214 :     print STDERR "\nGetting PEG FASTA sequences ...\n";
215 :     &FIG::run("get_fasta_for_tbl_entries $dir/contigs < $dir/Features/peg/tbl > $dir/Features/peg/fasta");
216 :     }
217 :     else
218 :     {
219 :     die "Features/peg/tbl is empty --- something is wrong";
220 :     }
221 :     print STDERR "\nSuccessfully completed processing $tarfile\n\n";
222 :     exit(0);
223 :    
224 :    
225 :    
226 :     sub read_fasta_record
227 :     {
228 :     my ($file_handle) = @_;
229 :     my ( $old_end_of_record, $fasta_record, @lines, $head, $sequence, $seq_id, $comment, @parsed_fasta_record );
230 :    
231 :     if (not defined($file_handle)) { $file_handle = \*STDIN; }
232 :    
233 :     $old_end_of_record = $/;
234 :     $/ = "\n>";
235 :    
236 :     if (defined($fasta_record = <$file_handle>))
237 :     {
238 :     chomp $fasta_record;
239 :     @lines = split( /\n/, $fasta_record );
240 :     $head = shift @lines;
241 :     $head =~ s/^>?//;
242 :     $head =~ m/^(\S+)/;
243 :     $seq_id = $1;
244 :    
245 :     if ($head =~ m/^\S+\s+(.*)$/) { $comment = $1; } else { $comment = ""; }
246 :    
247 :     $sequence = join( "", @lines );
248 :    
249 :     @parsed_fasta_record = ( $seq_id, \$sequence, $comment );
250 :     }
251 :     else
252 :     {
253 :     @parsed_fasta_record = ();
254 :     }
255 :    
256 :     $/ = $old_end_of_record;
257 :    
258 :     return @parsed_fasta_record;
259 :     }
260 :    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3