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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3