[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.3 - (view) (download) (as text)

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3