[BioC] unexpected behaviour of pairwiseAlignment() in Biostrings

[Patrick Aboyoun]

Avril,
This is a printing issue that I haven't sorted out yet to everyone's 
satisfaction. I had some global alignments that produced long end gaps 
so I chose not to print end gaps out using the standard show method. 
Instead, I used the numbers i,  j

pattern: [i]
pattern: [j]

to denote the location of the first match/mismatch in the pattern and 
subject string. I could revisit this decision if this problem 
represented the exception rather than the rule. BTW it appears that 
needle isn't scoring the end gaps. The score from pairwiseAlignment (37) 
penalized the end gaps while the score from needle (45) did not.


Patrick


Coghlan, Avril wrote:
> Dear all,
>
> I have been using the pairwiseAlignment() function in the Biostrings
> library for creating global alignments.
>
> However, I noticed that the way it behaves is not quite like I expected.
>
>
> For example, I wanted to create a global alignment of sequences
> "HEIAKGKAL" and "HEIAKGKALIIIIEALKCLA" so I typed:
>   
>> x1 <- "HEIAKGKAL"
>> y1 <- "HEIAKGKALIIIIEALKCLA"
>> pairwiseAlignment(x1,y1, substitutionMatrix = BLOSUM62, gapOpening =
>>     
> -1,gapExtension=-1)
> Global PairwiseAlignedFixedSubject (1 of 1)
> pattern: [1] HEIAKGKAL 
> subject: [1] HEIAKGKAL 
> score: 37
>
> The output that I got seems to be a local alignment, as it doesn't
> contain the whole of the two input sequences.
> Therefore, I am wondering if this function is using the Needleman-Wunsch
> algorithm to make the global alignment, or not?
> It looks to me like this actually the result of the Smith-Waterman
> algorithm for local alignment. 
>
> By the way, I still get the same result when I explicitly use the
> "type=global" option:
>   
>> pairwiseAlignment(x1,y1, substitutionMatrix = BLOSUM62, gapOpening =
>>     
> -1 ,gapExtension=-1,type='global')
> Global PairwiseAlignedFixedSubject (1 of 1)
> pattern: [1] HEIAKGKAL 
> subject: [1] HEIAKGKAL 
> score: 37
>
> Just to check that I am not going crazy, I tried aligning the same
> sequences using needle from the EMBOSS package (which does
> Needleman-Wunsch global alignments) at
> http://mobyle.pasteur.fr/cgi-bin/portal.py?form=needle with gapOpening
> and gapExtension penalties of -1, and the BLOSUM62 matrix. It gave the
> best global alignment as:
> x1                 1 HEIAKGKAL-------      9
>                      |||||||||       
> y1                 1 HEIAKGKALEALKCLA     16
> Score: 45.0
> This is something like I would have expected pairwiseAlignment() to give
> me, ie. a global alignment spanning the full lengths of the two
> sequences.
>
> Please could you explain why pairwiseAlignment is not giving me a global
> alignment? I am very confused, and am wondering if I am being very
> stupid and have misunderstood something (probably!)..
>
> Regards,
> Avril
>
> Avril Coghlan
> Dept. Microbiology
> University College Cork
> Ireland
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor
>