[Bioperl-l] Writing .scf files
Adam Sjøgren
asjo at koldfront.dk
Sat Jul 18 20:00:25 UTC 2009
Hi.
I'm trying to write an .scf file from just a sequence (to create
test-input). I can't seem to figure out how to do it correctly. This is
what I do:
#!/usr/bin/perl
use strict;
use warnings;
use Bio::Seq::Quality;
use Bio::SeqIO;
print "BioPerl version " . $Bio::Root::Version::VERSION . "\n";
print "Creating Bio::Seq::Quality object\n";
my $seq=Bio::Seq::Quality->new(
-qual=>'65 66 67 68 69 70 71 72 73 74 75 76',
-seq =>'atcgatcgatcg',
);
print "Writing .scf\n";
my $out=Bio::SeqIO->new(-file=>'>write.scf', -format=>'scf');
$out->write_seq(-target=>$seq);
print "Reading .scf\n";
my $in=Bio::SeqIO->new(-file=>'write.scf', -format=>'scf');
my $in_seq=$in->next_seq;
print "Done\n";
Basically following the pod of Bio::Seq::Quality and Bio::SeqIO::scf
tells me to.
Reading the scf gives me an exception; this is the output of the script:
Creating Bio::Seq::Quality object
Writing .scf
Reading .scf
--------------------- WARNING ---------------------
MSG: seq doesn't validate with [0-9A-Za-z\*\-\.=~\\/\?], mismatch is ,,
---------------------------------------------------
------------- EXCEPTION: Bio::Root::Exception -------------
MSG: Attempting to set the sequence to [AEI] which does not look healthy
STACK: Error::throw
STACK: Bio::Root::Root::throw /home/adsj/work/bioperl/bioperl-live/Bio/Root/Root.pm:368
STACK: Bio::PrimarySeq::seq /home/adsj/work/bioperl/bioperl-live/Bio/PrimarySeq.pm:283
STACK: Bio::PrimarySeq::new /home/adsj/work/bioperl/bioperl-live/Bio/PrimarySeq.pm:234
STACK: Bio::LocatableSeq::new /home/adsj/work/bioperl/bioperl-live/Bio/LocatableSeq.pm:122
STACK: Bio::Seq::Meta::Array::new /home/adsj/work/bioperl/bioperl-live/Bio/Seq/Meta/Array.pm:180
STACK: Bio::Seq::Quality::new /home/adsj/work/bioperl/bioperl-live/Bio/Seq/Quality.pm:206
STACK: Bio::SeqIO::scf::next_seq /home/adsj/work/bioperl/bioperl-live/Bio/SeqIO/scf.pm:245
STACK: ./t/nzdb/write_scf:23
-----------------------------------------------------------
If I obtain an .scf file from an ab1-file (using stadens convert_trace
to convert it), I can read the resulting .scf fine and I can write it
anew and read the file I wrote too.
So I guess I'm doing something wrong when creating my test-file, but I
must have stared me blind on the code, I can't seem to see it.
I am using a fresh BioPerl from svn (r15859.)
Any ideas?
Best regards,
Adam
--
"I said to Hank Williams: How lonely does it get? Adam Sjøgren
Hank Williams hasn't answered yet" asjo at koldfront.dk
More information about the Bioperl-l
mailing list