[Bioperl-l] Understanding LocatableSeq

[Jason Stajich]

On Tue, 4 Nov 2003, Wes Barris wrote:

> Ewan Birney wrote:
>
> >
> > On Tue, 4 Nov 2003, Wes Barris wrote:
> >
> >
> >>Hi,
> >>
> >>I am trying to create an msf alignment of several LocatableSeq objects.
> >>I have tried setting the "start" and "end" attributes of each LocatableSeq
> >>object before adding it to the alignment but the resulting sequences are
> >>still not aligned in the SimpleAlign object.  What am I doing wrong?
> >>
> >
> >
> > SimpleAlign does not automagically do the alignment - you need to call out
> > to Bio::Tools::Clustalw or TCoffee or something else. SimpleAlign
> > *represents* alignments, doesn't make them.
>
> I know that SimpleAlign does not automagically do the alignment.  That
> is why I am setting the "start" and "end" attributes.  I have an ACE
> file from clustal.  I am trying to write an msf file.  I have the
> sequences and the alignment information.  Now I want to write this
> into a SimpleAlign object.
>

You still need to prefix/postfix with the requisite number of gaps.

The start/end describe where the sequence participating the alignment
COMES FROM not where they are in the alignment, so you have to
explicitly code their alignment by placing the right number of gaps.


> >
> >
> > (it would be cute if simple align did this magically but probably would
> > lead to some fascinating bug-hunts as simple align booted up a full
> > progressive alignment engine when a sequence was added. Hmmm. Wistful
> > thinking...)
> >
> >
> >
> >
> >
> >>#!/usr/local/bin/perl -w
> >>#
> >>#
> >>use strict;
> >>use Bio::AlignIO;
> >>#
> >>my $aln = new Bio::SimpleAlign();
> >>my $lseq;
> >>$lseq = new Bio::LocatableSeq();
> >>$lseq->seq('GATCGATC');
> >>$lseq->id('this');
> >>$lseq->start(1);
> >>$lseq->end(8);
> >>$aln->add_seq($lseq);
> >>
> >>$lseq = new Bio::LocatableSeq();
> >>$lseq->seq('ATCGAT');
> >>$lseq->id('that');
> >>$lseq->start(2);
> >>$lseq->end(7);
> >>$aln->add_seq($lseq);
> >>
> >>my $outstream = new Bio::AlignIO(-format=>'msf', -file=>">junk.msf");
> >>$outstream->write_aln($aln);
> >>undef $outstream;
> >>
> >>
> >>The output looks like this:
> >>
> >>
> >>NoName   MSF: 2  Type: N  Tue Nov  4 09:01:58 2003  Check: 00 ..
> >>
> >>  Name: this/1-8  Len:    8  Check:  2590  Weight:  1.00
> >>  Name: that/2-7  Len:    6  Check:  1547  Weight:  1.00
> >>
> >>//
> >>
> >>this/1-8              GATCGATC
> >>that/2-7              ATCGAT	<--- Should't this be shifted one position to the right?
> >>
> >>
> >>I am using bioperl-1.2.3.
> >>--
> >>Wes Barris
> >>E-Mail: Wes.Barris at csiro.au
> >>
> >>
> >>_______________________________________________
> >>Bioperl-l mailing list
> >>Bioperl-l at portal.open-bio.org
> >>http://portal.open-bio.org/mailman/listinfo/bioperl-l
> >>
>
>
>

--
Jason Stajich
Duke University
jason at cgt.mc.duke.edu