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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3