[Bioperl-l] Bio::Seq::Quality description line problem

T.D. Houfek thoufek at pngg.org
Thu May 4 16:50:44 UTC 2006


Using Bioperl 1.5, having trouble with writing FASTA-style quality files 
using Bio::Seq::Quality.

I create the Bio::Seq::Quality object, giving its constructor an ID, a 
description, a nucleotide sequence, and a quality sequence. I then write 
the sequence FASTA and the quality FASTA. The description string will 
appear in the header line of the sequence FASTA, but not in the header 
line of the quality FASTA.

Can anybody help me figure out how to fix this? I've attached a sample 
script and output.

-T.D.

------------------- sample script follows 
---------------------------------------

#!/usr/bin/perl
use strict;
use Bio::Seq::Quality;
use Bio::SeqIO;

my $id = "bogus_id";
my $desc = "bogus description";
my $seq = "ATTATTATTATTATT";
my $qual = "10 20 30 10 20 30 10 20 30 10 20 30 10 20 30";

my $sequal_obj = Bio::Seq::Quality->new(
-display_id => $id,
-desc => $desc,
-seq => $seq,
-qual => $qual
);

my $qualout = Bio::SeqIO->new(
-file => ">myfile.qual",
-format => 'qual'
);
my $seqout = Bio::SeqIO->new(
-file => ">myfile.seq",
-format => 'Fasta'
);

$seqout->write_seq($sequal_obj);
$qualout->write_seq($sequal_obj);


------------------ sample output follows 
---------------------------------------

tdhoufek at aether:~$ cat myfile.seq
 >bogus_id bogus description
ATTATTATTATTATT
tdhoufek at aether:~$ cat myfile.qual
 >bogus_id
10 20 30 10 20 30 10 20 30 10 20 30 10 20 30

--------------------------------------------------------------------------------------------------




-- 
T.D. Houfek
senior bioinformatics developer
plant nematode genetics group
north carolina state university
Email: thoufek at pngg.org
----------------------------------------------------------
use Bio::Seq; @a =qw/NNN CCT GAG CAT GCG TGT AAG AAC TAG/;
$u=seq;$r=Bio::Seq;sub c{$c=$r->new(-$u=>"@_[0]")->revcom;
$t=$c->$u;}map{m/\d/?$g=c($a[$_]):tr/a-i/1-9/&&($g=$a[$_])
;$x[$i++]=$g;} split //,"dgh5cb40ab120cdefb4";$z=$r->new(-
$u=>(join"", at x))->translate()->$u;$z =~s/X/ /g;print"$z\n"





More information about the Bioperl-l mailing list