Question: Hello!

I am parsing a Genebank flat file from NCBI (the .gbff.gz found in [here](ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA_000020585.3_ASM2058v3/)). To do this I'm using Bio::SeqIO and I am finding a problem in the very first CDS which looks like this

source 1..4480937
/organism="Alteromonas macleodii str. 'Deep ecotype'"
/mol_type="genomic DNA"
/strain="Deep ecotype"
/db_xref="taxon:314275"
gene join(4480899..4480937,1..1626)
/locus_tag="MADE_000001020675"
CDS join(4480899..4480937,1..1626)
/locus_tag="MADE_000001020675"
/inference="EXISTENCE: similar to AA
sequence:RefSeq:YP_006797084.1"
/note="Derived by automated computational analysis using
gene prediction method: Protein Homology."
/codon_start=1
/transl_table=11
/product="chromosomal replication initiator protein DnaA"
/protein_id="AGV54147.1"
/db_xref="GI:543176245"
/translation="MSLWN...TLSS"

My concern is with this line:

CDS join(4480899..4480937,1..1626)

However, when I look at the nucleotide sequence that is outputed by Bio::SeqIO, it is outputting ~4.5 million nucleotides whereas if you look at the highlighted area it clearly shows that the coding sequence has ~1600ish nucleotides.What am I doing wrong?

Below you can find a snippet of the code that I am using. The reason why I'm using

next unless (($feat->primary_tag eq "CDS") and $feat->has_tag("translation"));

is that I am only interested in getting protein sequences and their corresponding gene sequences.

Thank you for all help!

#!/usr/local/bin/perl
use strict;
use warnings;
use autodie;
use Data::Dump;
use File::Basename;
use Bio::SeqIO;
use Data::Dumper;

[$gbffpath and other non-relevant stuff hidden]

opendir(DIR, "$gbffpath");
my @files= grep {/\.gbff$/ } readdir(DIR);
closedir(DIR);

my $k=1;
my ($i,$j);
foreach my $f (@files) {
$i=0;$j=0;
my $in = Bio::SeqIO -> new( -file => "<$gbffpath\\$f", -format => 'genbank');

while (my $seq = $in->next_seq) {
for my $feat ($seq -> get_SeqFeatures) {
next unless (($feat->primary_tag eq "CDS") and $feat->has_tag("translation"));
if ($feat->primary_tag eq "CDS") {#use to print DNA

print $geneout "\>";
# construct fasta header based on information present for each CDS:
print $geneout $feat->get_tag_values('locus_tag'),"\|" if ( $feat->has_tag('locus_tag') );
print $geneout $feat->get_tag_values('db_xref') ,"\|" if ( $feat->has_tag('db_xref') );
print $geneout $feat->get_tag_values('protein_id'),"\|" if ( $feat->has_tag('protein_id') );
print $geneout $feat->get_tag_values('gene'),"\|" if ( $feat->has_tag('gene') );
print $geneout $feat->get_tag_values('product'),"\|" if ( $feat->has_tag('product') );
print $protout $feat->get_tag_values('note'), if ($feat->has_tag('note') );

# print sequence:
print $geneout "\n",$feat->seq->seq,"\n";
$i++;
}
}
}
}

if (($i ne $j) or ($i==0) or ($j==0)){
print "\t\t$f\t$i CDS\t$j prots\t\n" ;
print "\t\t\>\>\>$f\t$i CDS\t$j prots\t\n";
}
else{
print "Extracted $i\/$j CDS\/proteins from file $f ($k/".($#files+1).")\n";
}
$k++;
}

Answer: Given how the `seq()` function is supposed to start from the beginning of the first range and go to the end of the last, you have an edge case here where your sequence fragments are at the ends. The method sees your `1..1626` and `4480899..4480937` and goes `1..4480937`

I think you'd be better off using `spliced_seq()` and joining the result from that call.

Source: http://www.bioperl.org/wiki/HOWTO:Feature-Annotation