[Bioperl-l] create scf file from SimpleAlign object

Eckhard Lehmann ecky@e-lehmann.de
Mon, 8 Apr 2002 21:22:06 +0200


Hi,

I have got a little problem:

There are two fasta sequences that are very similar to each other except of 
two or three nucleotide variants.
These two sequences should be aligned and at the end there should be a .scf 
file. When I open this .scf file with trev or another .scf viewer, there 
sould be two peaks at the positions where a nucleotide variant occurs.
The gap4 program from the staden package performs such files and I have to 
automatise the generation of such .scf files for many sequences...

With the following code snippets I get the two sequences aligned (with 
ClustalW)
-----------------
$fas1 = Bio::SeqIO->new('-file'=>"seq1.fasta");
$fas2 = Bio::SeqIO->new('-file'=>"seq2.fasta");

@seqs = ($fas1->next_seq(), $fas2->next_seq());
$seqs_ref = \@seqs;

@params = ('ktuple'=>2, 'matrix'=>'BLOSUM');

$factory = Bio::Tools::Run::Alignment::Clustalw->new(@params);
$aln = $factory->align($seqs_ref);
-----------------

Then I tried to write a .scf file this way:
-----------------
$q = "56 " x $aln->length();
$s = $aln->consensus_iupac();
$i = "aligned_sequences";

$qualSeq = Bio::Seq::SeqWithQuality->new(-qual=>$q, -seq=>$s, -id=>$i);

$out = Bio::SeqIO->new('-file'=>">output.scf", '-format'=>'scf');
$out->write_seq(-SeqWithQuality=>$qualSeq);
------------------

There comes out a .scf file but at the positions where the two peaks for the 
variants should appear, there is nothing (not even one peak). Probably the 
SeqIO::scf object does not recognize the M's and Y's and so on in the 
consensus_iupac - sequence from the SimpleAlign object.

Is there another posibility to get out an .scf file from two aligned 
sequences with bioperl, which contains exactly the same information as the 
.scf files generated from gap4/staden package?

Thank you in advance

Eckhard Lehmann