3
|
1 #!/usr/bin/perl
|
|
2
|
|
3 #######
|
|
4 # POD #
|
|
5 #######
|
|
6
|
|
7 =pod
|
|
8
|
|
9 =head1 NAME
|
|
10
|
|
11 C<cds_extractor.pl> - extract protein or DNA sequences from CDS features
|
|
12
|
|
13 =head1 SYNOPSIS
|
|
14
|
|
15 C<perl cds_extractor.pl -i seq_file.[embl|gbk] -p>
|
|
16
|
|
17 =head1 DESCRIPTION
|
|
18
|
|
19 Extracts protein or DNA sequences of CDS features from a
|
|
20 (multi)-RichSeq file (e.g. EMBL or GENBANK format) and writes them to a
|
|
21 multi-FASTA file. The FASTA headers for each CDS include either the
|
|
22 locus tag, if that's not available, protein ID, gene, or an internal
|
|
23 CDS counter as identifier (in this order). The organism info
|
|
24 includes also possible plasmid names. Pseudogenes (tagged by
|
|
25 B</pseudo>) are not included (except in the CDS counter).
|
|
26
|
|
27 In addition to the identifier, FASTA headers include gene (B<g=>),
|
|
28 product (B<p=>), organism (B<o=>), and EC numbers (B<ec=>), if these
|
|
29 are present for a CDS. Individual EC numbers are separated by
|
|
30 B<semicolons>. The location/position (B<l=>start..stop) of a CDS will
|
|
31 always be included. If gene is used as FASTA header ID
|
|
32 'B<g=>gene' will only be included with option B<-f>.
|
|
33
|
|
34 Fuzzy locations in the feature table of a sequence file are not
|
|
35 taken into consideration for B<l=>. If you set options B<-u>
|
|
36 and/or B<-d> and the feature location overlaps a B<circular>
|
|
37 replicon boundary, positions are marked with '<' or '>' in the
|
|
38 direction of the exceeded boundary. Features with overlapping
|
|
39 locations in B<linear> sequences (e.g. contigs) will be skipped and
|
|
40 are B<not> included in the output! A CDS feature is on the lagging
|
|
41 strand if start > stop in the location. In the special case of
|
|
42 overlapping circular sequence boundaries this is reversed.
|
|
43
|
|
44 Of course, the B<l=> positions are separate for each sequence in a
|
|
45 multi-sequence file. Thus, if you want continuous positions for the
|
|
46 CDSs run these files first through L<C<cat_seq.pl>|/cat_seq>.
|
|
47
|
|
48 Optionally, a file with locus tags can be given to extract only
|
|
49 these CDS features with option B<-l> (each locus tag in a new line).
|
|
50
|
|
51 =head1 OPTIONS
|
|
52
|
|
53 =head2 Mandatory options
|
|
54
|
|
55 =over 23
|
|
56
|
|
57 =item B<-i>=I<str>, B<-input>=I<str>
|
|
58
|
|
59 Input RichSeq sequence file including CDS annotation (e.g. EMBL or
|
|
60 GENBANK format)
|
|
61
|
|
62 =item B<-p>, B<-protein>
|
|
63
|
|
64 Extract B<protein> sequence for each CDS feature, excludes option B<-n>
|
|
65
|
|
66 B<or>
|
|
67
|
|
68 =item B<-n>, B<-nucleotide>
|
|
69
|
|
70 Extract B<nucleotide> sequence for each CDS feature, excludes option
|
|
71 B<-p>
|
|
72
|
|
73 =back
|
|
74
|
|
75 =head2 Optional options
|
|
76
|
|
77 =over 28
|
|
78
|
|
79 =item B<-h>, B<-help>
|
|
80
|
|
81 Help (perldoc POD)
|
|
82
|
|
83 =item B<-u>=I<int>, B<-upstream>=I<int>
|
|
84
|
|
85 Include given number of flanking nucleotides upstream of each CDS
|
|
86 feature, forces option B<-n>
|
|
87
|
|
88 =item B<-d>=I<int>, B<-downstream>=I<int>
|
|
89
|
|
90 Include given number of flanking nucleotides downstream of each CDS
|
|
91 feature, forces option B<-n>
|
|
92
|
|
93 =item B<-c>=I<str>, B<-cds_prefix>=I<str>
|
|
94
|
|
95 Prefix for the internal CDS counter [default = 'CDS']
|
|
96
|
|
97 =item B<-l>=I<str>, B<-locustag_list>=I<str>
|
|
98
|
|
99 List of locus tags to extract only those (each locus tag on a new line)
|
|
100
|
|
101 =item B<-f>, B<-full_header>
|
|
102
|
|
103 If gene is used as ID include additionally 'B<g=>gene' in FASTA
|
|
104 headers, so downstream analyses can recognize the gene tag (e.g.
|
|
105 L<C<prot_finder.pl>|/prot_finder>).
|
|
106
|
|
107 =item B<-v>, B<-version>
|
|
108
|
|
109 Print version number to C<STDERR>
|
|
110
|
|
111 =back
|
|
112
|
|
113 =head1 OUTPUT
|
|
114
|
|
115 =over 23
|
|
116
|
|
117 =item F<*.faa>
|
|
118
|
|
119 Multi-FASTA file of CDS protein sequences
|
|
120
|
|
121 B<or>
|
|
122
|
|
123 =item F<*.ffn>
|
|
124
|
|
125 Multi-FASTA file of CDS DNA sequences
|
|
126
|
|
127 =item (F<no_annotation_err.txt>)
|
|
128
|
|
129 Lists input files missing CDS annotation, script exited with B<fatal error> i.e. no FASTA output file
|
|
130
|
|
131 =item (F<double_id_err.txt>)
|
|
132
|
|
133 Lists input files with ambiguous FASTA IDs, script exited with B<fatal error> i.e. no FASTA output file
|
|
134
|
|
135 =item (F<locus_tag_missing_err.txt>)
|
|
136
|
|
137 Lists CDS features without locus tags
|
|
138
|
|
139 =item (F<linear_seq_cds_overlap_err.txt>)
|
|
140
|
|
141 Lists CDS features overlapping sequence border of a B<linear>
|
|
142 molecule, which are not included in the result multi-FASTA file
|
|
143
|
|
144 =back
|
|
145
|
|
146 =head1 EXAMPLES
|
|
147
|
|
148 =over
|
|
149
|
|
150 =item C<perl cds_extractor.pl -i seq_file.gbk -p -l locus_tags.txt>
|
|
151
|
|
152 =item C<perl cds_extractor.pl -i seq_file.embl -n -l locus_tags.txt -u 100 -d 20>
|
|
153
|
|
154 =item C<perl cds_extractor.pl -i Ecoli_MG1655.gbk -p -f -c MG1655>
|
|
155
|
|
156 =back
|
|
157
|
|
158 =head1 DEPENDENCIES
|
|
159
|
|
160 =over
|
|
161
|
|
162 =item B<BioPerl (L<http://www.bioperl.org>)>
|
|
163
|
|
164 Tested with BioPerl version 1.006923
|
|
165
|
|
166 =back
|
|
167
|
|
168 =head1 VERSION
|
|
169
|
|
170 0.7.1 update: 26-10-2015
|
|
171 0.1 24-05-2012
|
|
172
|
|
173 =head1 AUTHOR
|
|
174
|
|
175 Andreas Leimbach aleimba[at]gmx[dot]de
|
|
176
|
|
177 =head1 LICENSE
|
|
178
|
|
179 This program is free software: you can redistribute it and/or modify
|
|
180 it under the terms of the GNU General Public License as published by
|
|
181 the Free Software Foundation; either version 3 (GPLv3) of the
|
|
182 License, or (at your option) any later version.
|
|
183
|
|
184 This program is distributed in the hope that it will be useful, but
|
|
185 WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
186 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
187 General Public License for more details.
|
|
188
|
|
189 You should have received a copy of the GNU General Public License
|
|
190 along with this program. If not, see L<http://www.gnu.org/licenses/>.
|
|
191
|
|
192 =cut
|
|
193
|
|
194
|
|
195 ########
|
|
196 # MAIN #
|
|
197 ########
|
|
198
|
|
199 use warnings;
|
|
200 use strict;
|
|
201 use autodie;
|
|
202 use Getopt::Long;
|
|
203 use Pod::Usage;
|
|
204 use Bio::SeqIO; # bioperl module to handle sequence input/output
|
|
205 # use Bio::Seq; # bioperl module to play with the sequence and its features ### apparently not needed, methods inherited
|
|
206 # use Bio::SeqFeatureI; # bioperl module to handle features in a sequence ### apparently not needed, methods inherited
|
|
207
|
|
208
|
|
209
|
|
210 ### Get options with Getopt::Long, works also abbreviated and with two "--": -i, --i, -input ...
|
|
211 my $Seq_File; # RichSeq sequence file including feature annotation
|
|
212 my $Opt_Protein; # extract protein sequences for each CDS feature; excludes option '-n'
|
|
213 my $Opt_Nucleotide; # extract nucleotide sequences for each CDS feature; excludes option '-p'
|
|
214 my $Upstream = 0; # include given number of flanking nucleotides upstream of each CDS feature; forces option '-n'
|
|
215 my $Downstream = 0; # include given number of flanking nucleotides downstream of each CDS feature; forces option '-n'
|
|
216 my $Prefix; # prefix for the internal CDS counter
|
|
217 my $Locustag_List; # list of locus_tags to extract only those
|
|
218 my $Opt_Full_Header; # include a full FASTA header for downstream 'prot_finder.pl' analysis (needed to identify /gene feature tag if used as identifier)
|
|
219 my $VERSION = '0.7.1';
|
|
220 my ($Opt_Version, $Opt_Help);
|
|
221 GetOptions ('input=s' => \$Seq_File,
|
|
222 'protein' => \$Opt_Protein,
|
|
223 'nucleotide' => \$Opt_Nucleotide,
|
|
224 'upstream:i' => \$Upstream,
|
|
225 'downstream:i' => \$Downstream,
|
|
226 'cds_prefix:s' => \$Prefix,
|
|
227 'locustag_list:s' => \$Locustag_List,
|
|
228 'full_header' => \$Opt_Full_Header,
|
|
229 'version' => \$Opt_Version,
|
|
230 'help|?' => \$Opt_Help);
|
|
231
|
|
232
|
|
233
|
|
234 ### Run perldoc on POD
|
|
235 pod2usage(-verbose => 2) if ($Opt_Help);
|
|
236 die "$0 $VERSION\n" if ($Opt_Version);
|
|
237 if (!$Seq_File) {
|
|
238 my $warning = "\n### Fatal error: Option '-i' or its argument is missing!\n";
|
|
239 pod2usage(-verbose => 1, -message => $warning, -exitval => 2);
|
|
240 }
|
|
241
|
|
242
|
|
243
|
|
244 ### Enforce mandatory or optional options
|
|
245 if (($Upstream || $Downstream) && !$Opt_Nucleotide) {
|
|
246 warn "Option '-d' and/or '-u' set, but not '-n'. Forcing option '-n'!\n";
|
|
247 $Opt_Nucleotide = 1;
|
|
248 undef $Opt_Protein;
|
|
249 }
|
|
250 if (!$Opt_Protein && !$Opt_Nucleotide) {
|
|
251 die "\n### Fatal error: None of the mandatory options ('-p' or '-n') given! Please choose one of the extraction methods!\n";
|
|
252 } elsif ($Opt_Protein && $Opt_Nucleotide) {
|
|
253 die "\n### Fatal error: Both mandatory options ('-p' and '-n') given! Choose only one of the extraction methods!\n";
|
|
254 }
|
|
255 $Prefix = 'CDS' if (!$Prefix); # default for internal CDS
|
|
256
|
|
257
|
|
258
|
|
259 ### Create a bioperl Bio::SeqIO object with $Seq_File
|
|
260 my $seqio_object = Bio::SeqIO->new(-file => "<$Seq_File"); # no '-format' to leave to bioperl guessing
|
|
261
|
|
262
|
|
263
|
|
264 ### Print the input and output file which are processed and written, respectively
|
|
265 print "\nInput: $Seq_File\t";
|
|
266 my $Ori_Filename = $Seq_File; # store original filename to give error statements
|
|
267 $Seq_File =~ s/(.+)\.\w+$/$1/; # strip filename extension
|
|
268 my $Out_Fh; # filehandle for output file
|
|
269 if ($Opt_Protein) {
|
|
270 $Seq_File = $Seq_File.'.faa';
|
|
271 print "Output: $Seq_File\n";
|
|
272 open ($Out_Fh, ">", "$Seq_File");
|
|
273 } elsif ($Opt_Nucleotide) {
|
|
274 $Seq_File = $Seq_File.'.ffn';
|
|
275 print "Output: $Seq_File\n";
|
|
276 open ($Out_Fh, ">", "$Seq_File");
|
|
277 }
|
|
278
|
|
279
|
|
280
|
|
281 ### Get the list of locus_tags if $Locustag_List defined
|
|
282 my %Locus_Tags; # store locus_tags and indicate if found in $Seq_File
|
|
283 if ($Locustag_List) {
|
|
284 open (my $locus_fh, "<", "$Locustag_List");
|
|
285 while (<$locus_fh>) {
|
|
286 chomp;
|
|
287 next if (/^$/); # skip empty lines
|
|
288 $Locus_Tags{$_} = 0; # changes to 1 if the locus tag was found in $Seq_File, see below
|
|
289 }
|
|
290 close $locus_fh;
|
|
291 }
|
|
292
|
|
293
|
|
294
|
|
295 ### Write output FASTA with respective header lines and either protein seq from /translation feature tag or nucleic acid subseq
|
|
296 my $Anno_Present = 0; # switches to true if CDS primary features present
|
|
297 my $Organism; # store /organism tag value to include in each header, include plasmid name
|
|
298 my $Cds_Count = sprintf("%04d", 0); # internal CDS counter, counts also 'pseudo' CDSs
|
|
299 my %Double_Id; # control if a locus_tag or protein_id is ambiguous in $Seq_File and die (they should be unique); see subroutine 'control_double'
|
|
300
|
|
301 my $No_Locus_Tag = 0; # count how many CDS features don't have a locus_tag
|
|
302 my $Locus_Tag_Miss_Err = 'locus_tag_missing_err.txt'; # error file that contains all alternative IDs (see subroutine 'locus_tag_missing')
|
|
303 my $Locus_Tag_Miss_Err_Fh; # filehandle
|
|
304
|
|
305 my $Linear_Seq_Overlap = 0; # count CDS features which overlap linear sequence borders (see subroutine 'linear_overlap')
|
|
306 my $Linear_Overlap_Err = 'linear_seq_cds_overlap_err.txt'; # error file that contains all affected IDs
|
|
307 my $Linear_Seq_Overlap_Fh; # filehandle
|
|
308
|
|
309 while (my $seq_object = $seqio_object->next_seq) { # a Bio::Seq object
|
|
310 my @feat_objects = $seq_object->get_SeqFeatures; # slurp all features from the seq to check CDS primary feature annotation
|
|
311
|
|
312 if (grep ($_->primary_tag eq 'CDS', @feat_objects) == 0) { # scalar context
|
|
313 next; # if no CDS features present skip to next $seq_object in (multi-)seq file
|
|
314 }
|
|
315 $Anno_Present = 1; # if not skipped CDS annotation present
|
|
316
|
|
317 $Organism = ''; # empty $Organism for each $seq_object
|
|
318 foreach my $feat_object (@feat_objects) { # a Bio::SeqFeatureI object
|
|
319 if ($Locustag_List) { # skip residual FEATURE objects if all locus_tags in $Locustag_List found
|
|
320 if (grep ($Locus_Tags{$_} == 1, keys %Locus_Tags) == keys %Locus_Tags) { # scalar context
|
|
321 last;
|
|
322 }
|
|
323 }
|
|
324
|
|
325 if ($feat_object->primary_tag eq 'source') {
|
|
326 $Organism = feature_value_eval($feat_object, 'organism'); # subroutine to evaluate existence of the tag and return value
|
|
327 if ($feat_object->has_tag('plasmid')) { # subroutine 'feature_value_eval' also possible, but method 'has_tag' more elegantly
|
|
328 my ($plasmid) = $feat_object->get_tag_values('plasmid'); # values always returned as ARRAYS
|
|
329 $plasmid =~ s/\s/_/g;
|
|
330 $Organism = $Organism.'-plasmid_'.$plasmid;
|
|
331 }
|
|
332 }
|
|
333
|
|
334 if ($feat_object->primary_tag eq 'CDS') {
|
|
335 $Cds_Count++;
|
|
336 if ($feat_object->has_tag('pseudo')) { # skip pseudogenes, they don't include '/translation'
|
|
337 next;
|
|
338 }
|
|
339 my ($start, $stop, $strand, $seq_stop, $location_stop) = get_location($feat_object, $seq_object); # subroutine to get start, stop and strand of feature; $seq_stop needed for sub 'print_seq' and with $location_stop needed for sub 'print_location'
|
|
340
|
|
341 if ($feat_object->has_tag('locus_tag')) { # LOCUS_TAG-ELSIF
|
|
342 my ($locus_tag) = $feat_object->get_tag_values('locus_tag');
|
|
343 control_double($locus_tag, 'locus_tag'); # subroutine to control uniqueness of identifier (here value of /locus_tag)
|
|
344
|
|
345 if ($Locustag_List) { # if '-locustag_list' is set only get those CDSs
|
|
346 my ($locus) = grep (/$locus_tag/i, keys %Locus_Tags); # case-insensitive for typing mistakes
|
|
347 if ($locus) {
|
|
348 $Locus_Tags{$locus} = 1;
|
|
349 if ($Opt_Nucleotide && ($start eq 'not_circular' || ($location_stop && $location_stop eq 'not_circular'))) { # CDS feature overlaps seq border of linear molecule
|
|
350 linear_overlap($locus_tag); # subroutine to write overlap CDSs to '$Linear_Overlap_Err'
|
|
351 next; # jump to the next feature object
|
|
352 }
|
|
353 print_fasta_ID($feat_object, $locus_tag); # subroutine to print the identifier (here /locus_tag), /gene (g=) and /product (p=) (both if present) to the FASTA header
|
|
354 print_location($start, $stop, $strand, $seq_stop, $location_stop); # subroutine to print feature location/position (l=) to FASTA header
|
|
355 print_org_ec($feat_object); # subroutine to print /organism (o=) and /EC_numbers (ec=) (both if present) to FASTA header
|
|
356 print_seq($feat_object, $start, $stop, $strand, $seq_stop, $seq_object); # subroutine to print the protein or nucleic sequence
|
|
357 }
|
|
358 next; # jump to the next feature object
|
|
359 }
|
|
360
|
|
361 ### '$Locustag_List' not defined
|
|
362 if ($Opt_Nucleotide && ($start eq 'not_circular' || ($location_stop && $location_stop eq 'not_circular'))) {
|
|
363 linear_overlap($locus_tag); # subroutine
|
|
364 next; # jump to the next feature object
|
|
365 }
|
|
366 print_fasta_ID($feat_object, $locus_tag); # subroutine
|
|
367
|
|
368 } elsif ($Locustag_List) { # in case a list of locus_tags is given, no need to look at CDSs that don't have a /locus_tag
|
|
369 next; # jump to the next feature object
|
|
370
|
|
371 } elsif ($feat_object->has_tag('protein_id')) { # PROTEIN_ID-ELSIF
|
|
372 my ($protein_id) = $feat_object->get_tag_values('protein_id');
|
|
373 control_double($protein_id, 'protein_id'); # subroutine
|
|
374 locus_tag_missing($protein_id, 'protein_id'); # subroutine to inform no /locus_tag is present for this CDS
|
|
375 $No_Locus_Tag++;
|
|
376 if ($Opt_Nucleotide && ($start eq 'not_circular' || ($location_stop && $location_stop eq 'not_circular'))) {
|
|
377 linear_overlap($protein_id); # subroutine
|
|
378 next;
|
|
379 }
|
|
380 print_fasta_ID($feat_object, $protein_id); # subroutine
|
|
381
|
|
382 } elsif ($feat_object->has_tag('gene')) { # GENE-ELSIF
|
|
383 my ($gene) = $feat_object->get_tag_values('gene');
|
|
384 control_double($gene, 'gene'); # subroutine
|
|
385 locus_tag_missing($gene, 'gene'); # subroutine
|
|
386 $No_Locus_Tag++;
|
|
387 if ($Opt_Nucleotide && ($start eq 'not_circular' || ($location_stop && $location_stop eq 'not_circular'))) {
|
|
388 linear_overlap($gene); # subroutine
|
|
389 next;
|
|
390 }
|
|
391 print_fasta_ID($feat_object, ''); # subroutine; give empty string for $tag to sub, to be able to determine that the sub call is from GENE-ELSIF and 'g=' should only be included if option '-f' is set
|
|
392
|
|
393 } else { # if none of the above tags are existent use the internal CDS counter; CDS_COUNT-ELSE
|
|
394 my $cds_prefix_count = $Prefix.'_'.$Cds_Count;
|
|
395 locus_tag_missing($cds_prefix_count, 'cds-counter'); # subroutine
|
|
396 $No_Locus_Tag++;
|
|
397 if ($Opt_Nucleotide && ($start eq 'not_circular' || ($location_stop && $location_stop eq 'not_circular'))) {
|
|
398 linear_overlap($cds_prefix_count); # subroutine
|
|
399 next;
|
|
400 }
|
|
401 print_fasta_ID($feat_object, $cds_prefix_count); # subroutine
|
|
402 }
|
|
403
|
|
404 ### Print location, organism, EC_number(s), and sequence for all $feat_object if '$Locustag_List' not defined
|
|
405 print_location($start, $stop, $strand, $seq_stop, $location_stop); # subroutine
|
|
406 print_org_ec($feat_object); # subroutine
|
|
407 print_seq($feat_object, $start, $stop, $strand, $seq_stop, $seq_object); # subroutine
|
|
408 }
|
|
409 }
|
|
410
|
|
411 if ($Locustag_List) { # skip ALSO residual SEQUENCE objects if all locus_tags are found
|
|
412 if (grep ($Locus_Tags{$_} == 1, keys %Locus_Tags) == keys %Locus_Tags) {
|
|
413 last;
|
|
414 }
|
|
415 }
|
|
416 }
|
|
417
|
|
418 close $Out_Fh;
|
|
419 if ($No_Locus_Tag > 0) {
|
|
420 warn "### $No_Locus_Tag CDS feature(s) don't have a locus tag in file '$Ori_Filename'. The respective ID(s) are written to the error file '$Locus_Tag_Miss_Err'!\n";
|
|
421 close $Locus_Tag_Miss_Err_Fh;
|
|
422 }
|
|
423 if ($Linear_Seq_Overlap > 0) {
|
|
424 warn "### $Linear_Seq_Overlap CDS feature(s) overlap a sequence border of a linear molecule in '$Ori_Filename' because of your settings for '-u' and/or '-d'. These CDS feature(s) are not included in the output file '$Seq_File'! The respective ID(s) are written to the error file '$Linear_Overlap_Err'!\n";
|
|
425 close $Linear_Seq_Overlap_Fh;
|
|
426 }
|
|
427
|
|
428
|
|
429
|
|
430 ### Exit with error and remove empty result file if no CDS features/annotation found in (multi-)seq file
|
|
431 if (!$Anno_Present) {
|
|
432 my $anno_err = 'no_annotation_err.txt';
|
|
433 err_file_exist($anno_err); # subroutine to test for file existence and give warning to STDERR
|
|
434 open (my $anno_err_fh, ">>", "$anno_err");
|
|
435 print $anno_err_fh "$Ori_Filename\n";
|
|
436 close $anno_err_fh;
|
|
437 unlink $Seq_File;
|
|
438 die "\n###Fatal error\nNo CDS annotation in '$Ori_Filename', the respective filename was written to '$anno_err'!\nExiting program!\n\n";
|
|
439 }
|
|
440
|
|
441
|
|
442
|
|
443 ### Print locus tags that were not found in seq_file
|
|
444 if ($Locustag_List) {
|
|
445 my @missed = grep ($Locus_Tags{$_} == 0, keys %Locus_Tags);
|
|
446 if (@missed) {
|
|
447 print "### The following locus tags were not found in '$Ori_Filename':\n";
|
|
448 foreach (sort @missed) {
|
|
449 print "$_\t";
|
|
450 }
|
|
451 print "\n";
|
|
452 }
|
|
453 }
|
|
454
|
|
455
|
|
456 exit;
|
|
457
|
|
458
|
|
459 ###############
|
|
460 # Subroutines #
|
|
461 ###############
|
|
462
|
|
463 ### Control if an identifier for the FASTA header is ambiguous in the $Seq_File and die (they should be unique/unambiguous)
|
|
464 sub control_double {
|
|
465 my ($tag, $type) = @_;
|
|
466 if ($Double_Id{$tag}) {
|
|
467 my $double_id_err = 'double_id_err.txt';
|
|
468 my $file_exist = err_file_exist($double_id_err); # subroutine
|
|
469 open (my $double_id_err_fh, ">>", "$double_id_err");
|
|
470 print $double_id_err_fh "Original_filename\tID-type\tDouble_ID\n" if (!$file_exist); # print header in file only if file doesn't exist
|
|
471 print $double_id_err_fh "$Ori_Filename\t$type\t$tag\n";
|
|
472 close $double_id_err_fh;
|
|
473 unlink $Seq_File;
|
|
474 die "\n###Fatal error!\n'$tag' of type '$type' exists at least two times in organism '$Organism' of file '$Ori_Filename', but should be unambiguous for downstream analyses! The error is written to the error file '$double_id_err'. Please modify all repetitive occurences.\nExiting program!\n\n";
|
|
475 } else {
|
|
476 $Double_Id{$tag} = 1;
|
|
477 }
|
|
478 return 1;
|
|
479 }
|
|
480
|
|
481
|
|
482
|
|
483 ### Test for error file existence and give warning to STDERR
|
|
484 sub err_file_exist {
|
|
485 my $file = shift;
|
|
486 if (-e $file) {
|
|
487 warn "\nThe error file '$file' exists already, the current errors will be appended to the existing file!\n";
|
|
488 return 1;
|
|
489 }
|
|
490 return 0;
|
|
491 }
|
|
492
|
|
493
|
|
494
|
|
495 ### Return value of a tag (replace whitespaces with '_') or catch error if non-existent to return empty value
|
|
496 sub feature_value_eval {
|
|
497 my ($feat_object, $tag) = @_;
|
|
498 my $value = '';
|
|
499 eval {$value = join(';', $feat_object->get_tag_values($tag));}; # values always returned as ARRAYS; catch error if tag doesn't exist. To seperate /EC_numbers ';' is used, the only feature tag needed with more than one occurrence in a single CDS
|
|
500 $value =~ s/\s/_/g; # replace spaces with '_' (guess only needed for /organism and /product values)
|
|
501 return $value;
|
|
502 }
|
|
503
|
|
504
|
|
505
|
|
506 ### Get the position/location of a feature
|
|
507 sub get_location {
|
|
508 my ($feat_object, $seq_object) = @_;
|
|
509 my ($start, $stop, $strand, $seq_stop, $location_stop); # $seq_stop needed for sub 'print_seq' and with $location_stop needed for sub 'print_location'
|
|
510
|
|
511 $strand = $feat_object->strand; # 1 = feature on leading strand, -1 = feature on lagging strand
|
|
512 if ($strand == 1) {
|
|
513 $start = $feat_object->start - $Upstream;
|
|
514 $stop = $feat_object->end + $Downstream;
|
|
515 } elsif ($strand == -1) {
|
|
516 $start = $feat_object->start - $Downstream;
|
|
517 $stop = $feat_object->end + $Upstream;
|
|
518 }
|
|
519
|
|
520 ### Adjust positions if they overlap replicon sequence borders
|
|
521 if ($start < 0) { # '$seq_obj->subseq|->trunc' (in subroutine 'print_seq') don't work with negative numbers, but with positions > '$seq_obj->length'
|
|
522 if ($seq_object->is_circular) { # true if molecule is circular
|
|
523 $start = $seq_object->length - abs($start);
|
|
524 } else { # if sequence linear wrong to print overlapping sequences (see sub 'linear_overlap')
|
|
525 $start = 'not_circular';
|
|
526 }
|
|
527 $seq_stop = $seq_object->length + $stop;
|
|
528 }
|
|
529 if ($stop > $seq_object->length) { # for sub 'print_location' a stop > '$seq_obj->length' is no use and should start again from the seq start
|
|
530 if ($seq_object->is_circular) {
|
|
531 $location_stop = $stop - $seq_object->length;
|
|
532 } else {
|
|
533 $location_stop = 'not_circular';
|
|
534 }
|
|
535 }
|
|
536
|
|
537 return ($start, $stop, $strand, $seq_stop, $location_stop);
|
|
538 }
|
|
539
|
|
540
|
|
541
|
|
542 ### Inform CDS locations with options '-u' and/or '-d' are overlapping a sequence border of a linear molecule
|
|
543 sub linear_overlap {
|
|
544 my $tag = shift;
|
|
545 if ($Linear_Seq_Overlap == 0) { # open error file to append errors only for first overlap error
|
|
546 my $file_exist = err_file_exist($Linear_Overlap_Err); # subroutine
|
|
547 open ($Linear_Seq_Overlap_Fh, ">>", "$Linear_Overlap_Err");
|
|
548 print $Linear_Seq_Overlap_Fh "Original_filename\tOrganism\tID\n" if (!$file_exist); # header of error file
|
|
549 }
|
|
550 $Linear_Seq_Overlap++;
|
|
551 print $Linear_Seq_Overlap_Fh "$Ori_Filename\t$Organism\t$tag\n";
|
|
552 return 1;
|
|
553 }
|
|
554
|
|
555
|
|
556
|
|
557 ### Inform if locus_tags are missing
|
|
558 sub locus_tag_missing {
|
|
559 my ($tag, $type) = @_;
|
|
560 if ($No_Locus_Tag == 0) { # open error file to append errors only for first missing locus_tag occurrence
|
|
561 my $file_exist = err_file_exist($Locus_Tag_Miss_Err); # subroutine
|
|
562 open ($Locus_Tag_Miss_Err_Fh, ">>", "$Locus_Tag_Miss_Err");
|
|
563 print $Locus_Tag_Miss_Err_Fh "Original_filename\tOrganism\tAlternative_ID_type\tAlternative_ID\n" if (!$file_exist); # header of error file
|
|
564 }
|
|
565 print $Locus_Tag_Miss_Err_Fh "$Ori_Filename\t$Organism\t$type\t$tag\n";
|
|
566 return 1;
|
|
567 }
|
|
568
|
|
569
|
|
570
|
|
571 ### Print identifier, /gene (g=) and /product (p=) (both if existent) of a CDS to FASTA header
|
|
572 sub print_fasta_ID {
|
|
573 my ($feat_object, $tag) = @_;
|
|
574 my $gene = feature_value_eval($feat_object, 'gene'); # subroutine
|
|
575 my $product = feature_value_eval($feat_object, 'product'); # subroutine
|
|
576 if ($tag) { # $tag is defined for LOCUS_TAG-ELSIF, PROTEIN_ID-ELSIF, and CDS_COUNT-ELSE
|
|
577 print $Out_Fh ">$tag";
|
|
578 } elsif (!$tag) { # subroutine call from GENE-ELSIF with empty string for $tag (see below 'g=')
|
|
579 print $Out_Fh ">$gene";
|
|
580 }
|
|
581 print $Out_Fh " g=$gene" if (($gene && $tag) || ($Opt_Full_Header && !$tag)); # print 'g='-gene only if /gene tag present; additionally, only print 'g=' for GENE-ELSIF if '-f' is set, otherwise can't determine if the identifier is a gene in downstream tools (e.g. 'prot_finder.pl')
|
|
582 print $Out_Fh " p=$product" if ($product);
|
|
583 return 1;
|
|
584 }
|
|
585
|
|
586
|
|
587
|
|
588 ### Print position/location (l=) of a CDS to FASTA header
|
|
589 sub print_location {
|
|
590 my ($start, $stop, $strand, $seq_stop, $location_stop) = @_;
|
|
591 $stop = $location_stop if ($location_stop); # see sub 'get_location'
|
|
592 print $Out_Fh " l=";
|
|
593 if ($strand == 1) {
|
|
594 print $Out_Fh ">" if ($seq_stop); # defined if feature start position overlaps replicon's seq start (with '-u' or '-d'), sub 'get_location'
|
|
595 print $Out_Fh "$start..$stop";
|
|
596 print $Out_Fh "<" if ($location_stop); # defined if feature stop position overlaps replicon's seq end (with '-u' or '-d')
|
|
597 return 1;
|
|
598 } elsif ($strand == -1) {
|
|
599 print $Out_Fh "<" if ($location_stop); # switched on lagging strand, because $stop and $start switched
|
|
600 print $Out_Fh "$stop..$start"; # switch $stop and $start to indicate position on lagging strand
|
|
601 print $Out_Fh ">" if ($seq_stop);
|
|
602 return -1;
|
|
603 }
|
|
604 return 0;
|
|
605 }
|
|
606
|
|
607
|
|
608
|
|
609 ### Print organism (o=) and EC numbers (ec=) (both if existent) of a CDS to FASTA header
|
|
610 sub print_org_ec {
|
|
611 my $feat_object = shift;
|
|
612 print $Out_Fh " o=$Organism" if ($Organism);
|
|
613 my $ec_numbers = feature_value_eval($feat_object, 'EC_number'); # subroutine
|
|
614 print $Out_Fh " ec=$ec_numbers" if ($ec_numbers);
|
|
615 print $Out_Fh "\n";
|
|
616 return 1;
|
|
617 }
|
|
618
|
|
619
|
|
620
|
|
621 ### Print the protein or nucleic sequence to the result file
|
|
622 sub print_seq {
|
|
623 my ($feat_object, $start, $stop, $strand, $seq_stop, $seq_object) = @_;
|
|
624 if ($Opt_Protein) {
|
|
625 print $Out_Fh $feat_object->get_tag_values('translation'), "\n";
|
|
626
|
|
627 } elsif ($Opt_Nucleotide) {
|
|
628 $stop = $seq_stop if ($seq_stop); # see sub 'get_location'
|
|
629 if ($strand == 1) {
|
|
630 my $subseq = $seq_object->subseq($start, $stop);
|
|
631 print $Out_Fh substr($subseq, 0, $Upstream); # upstream bases, not possible to include in '$feat_obj->spliced_seq'
|
|
632 print $Out_Fh $feat_object->spliced_seq->seq; # to include also features with "join"/fuzzy locations; '$feat_object->spliced_seq' returns a Bio::Seq object thus call funtion '->seq' for the seq string
|
|
633 print $Out_Fh substr($subseq, length($subseq) - $Downstream, $Downstream), "\n"; # downstream bases
|
|
634
|
|
635 } elsif ($strand == -1) {
|
|
636 my $trunc_obj = $seq_object->trunc($start, $stop); # a Bio::Seq object, needed for revcom
|
|
637 my $rev_obj = $trunc_obj->revcom; # a Bio::Seq object
|
|
638 print $Out_Fh substr($rev_obj->seq, 0, $Upstream);
|
|
639 print $Out_Fh $feat_object->spliced_seq->seq;
|
|
640 print $Out_Fh substr($rev_obj->seq, $rev_obj->length - $Downstream, $Downstream), "\n";
|
|
641 }
|
|
642 }
|
|
643 return 1;
|
|
644 }
|