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

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3