Printing sequences of a FASTA file using BioPerl

Question

I have this code

use warnings;
use Getopt::Long;
use Bio::SeqIO;

GetOptions("in=s" => \$file) or die("Error in command line arguments\n");
open $new3, ">", "sequences_tmp.tab";
$seqin = Bio::SeqIO->new(-file => $file, -format => "Fasta");
$seqout = Bio::SeqIO->new(-file => ">$new3", -format => "tab");

while ($seq = $seqin->next_seq()) {
    $seqout->width($seq->length);
    $obj = $seq->id ."\t".$seq->seq()."\n";
    $seqout->write_seq($obj);
}

close $new3;

expecting to print the sequences this way seq_id TAB sequence. However, this code prints an empty file. Do you know what's going on?

The docs suggest the -file flag needs a filename rather than a filehandle — RobEarl, Jun 12 '14 at 13:40
You should *always* `use strict` as well as `use warnings` at the top of your Perl programs. It is only good manners when you are asking for help with your code to have availed yourself of this obvious aid to debugging built into the language itself. You should also test whether every `open` call succeeded, and use a `die` string that incloudes the `$!` built-in variable so that it is clear *why* the operation failed. In many cases it is simplest to simply `use autodie`, when the necessary checks will be included implicitly — Borodin, Jun 12 '14 at 14:36

BarneySchmale · Answer 1 · 2014-06-12T16:29:16.807

The $obj variable looks useless to me. It is a string, not a sequence object. As you only want to reformat the sequence, you can simple pass $seq to the write_seq()method.

So I wondering wether you are executing the loop body at all. You can print debug output for verifying that. If the loop body is not executed, than make sure that your input file really contains a sequence in FASTA format.

Also please declare use strict;on top of your script. It will help you to avoid many pitfalls.

score 0 · Answer 2 · answered Jun 12 '14 at 17:58

You open a filehandle $new3, but then stringify it and use it as a filename in your -file parameter. This is a bug.

open $new3, ">", "sequences_tmp.tab";
$seqout = Bio::SeqIO->new(
    -file    => ">$new3",      # <--- Not what you want
    -format  => "tab",
);

Bio::SeqIO->new can accept either a filehandle -fh or a filename -file as an initialization parameter. Therefore either of the following would likely work for you:

my $seqout = Bio::SeqIO->new(
    -fh      => $new3,
    -format  => "tab",
);

      #or#

my $seqout = Bio::SeqIO->new(
    -file    => '>sequences_tmp.tab',
    -format  => "tab",
);

Your code could use further cleaning up as well:

Include use strict; and use warnings; in EVERY perl script.
Include use autodie; anytime you're doing file processing.

Applying these and removing the likely debugging artifacts from your code reduces it to:

use strict;
use warnings;
use autodie;

use Getopt::Long;
use Bio::SeqIO;

GetOptions(
    "in=s"    => \my $infile,
) or die "Error in command line arguments\n";

my $outfile = "sequences_tmp.tab";

my $seqin  = Bio::SeqIO->new(-file => $infile,     -format => "Fasta");
my $seqout = Bio::SeqIO->new(-file => ">$outfile", -format => "tab");

while (my $seq = $seqin->next_seq()) {
    $seqout->write_seq($seq);
}

Printing sequences of a FASTA file using BioPerl

2 Answers2