[Bioperl-l] need help urgently - needle output parsing

[neeti somaiya]

Hi,

I am using needle alignment tool (standalone, on a linux machine), and then
I am using Bioperl to parse the output.
All data - sequence files and alignment outputs are attached with this mail.

I have 2 small sequences :- 693.seq and revcomp693.seq
I have 2 big sequences :- 80768-4291-5639.84809_84810_84809_1.scf.seq and
80768-4291-5639.84809_84810_84810_1.scf.seq
All these are in fasta format

Now I am doing the following :-
1) Aligning 80768-4291-5639.84809_84810_84809_1.scf.seq and 693.seq - output
file is 80768-4291-5639.84809_84810_84809_1.scf.out
parsing the output gives me the alignment start in 'traceseq' as 97
2) Aligning 80768-4291-5639.84809_84810_84809_1.scf.seq and revcomp693.seq -
output file is 80768-4291-5639.84809_84810_84809_1.scf.comp.out
parsing the output gives me the alignment start in 'traceseq' as 91

All this is correct.

Now I am doing the following :-
1) Aligning 80768-4291-5639.84809_84810_84810_1.scf.seq and 693.seq - output
file is 80768-4291-5639.84809_84810_84810_1.scf.out
parsing the output gives me the alignment start in 'traceseq' as 341 (this
is correct)
2) Aligning 80768-4291-5639.84809_84810_84810_1.scf.seq and revcomp693.seq -
output file is 80768-4291-5639.84809_84810_84810_1.scf.comp.out
parsing the output gives me the alignment start in 'traceseq' as 341 (this
is incorrect, correct position is 330)


Part of my code is as follows :-
---------------------------------------------

# running needle
`$needle_path./needle $trace.seq $snp_position_on_con.seq -gapopen
10.0-gapextend
0.5 $output`;

# parsing needle output
my $str = Bio::AlignIO->new(-format => 'emboss',-file => $output);
my $aln = $str->next_aln();
my $pos = $aln->column_from_residue_number('original',1);

$logger->info("Alignment pos is $pos");

####################################

 # running needle
`$needle_path./needle $trace.seq revcomp$snp_position_on_con.seq -gapopen
10.0 -gapextend 0.5 $comp_output`;

# parsing needle output
my $comp_str = Bio::AlignIO->new(-format => 'emboss',-file => $comp_output);
my $comp_aln = $comp_str->next_aln();
my $comp_pos = $comp_aln->column_from_residue_number('revcomp',1);

$logger->info("Alignment pos is $comp_pos");


Can someone please tell me what is going wrong here?


-- 
-Neeti
Even my blood says, B positive
-------------- next part --------------
A non-text attachment was scrubbed...
Name: data.zip
Type: application/zip
Size: 4456 bytes
Desc: not available
URL: <http://lists.open-bio.org/pipermail/bioperl-l/attachments/20070223/21658b7d/attachment-0004.zip>