[Bioperl-l] Writing .scf files

[Adam Sjøgren]

  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