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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3