[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