[BioC] Putative bug in Biostrings

[Hervé Pagès]

Hi Watal,

Hope you don't mind that I'm cc'ing the Bioconductor mailing list.

On 07/10/2014 10:50 PM, Watal M. Iwasaki wrote:
> Dear Mr. Pages,
>
> I encountered a strange behavior in Biostrings::mismatchTable(). It
> seems to throw an error when the PairwiseAlignments argument contains
> the sequences that differ in the last two consecutive nucleotides. For
> example,
>
>> mismatchTable(PairwiseAlignments('ATGC', 'ATAT'))
> Error in fromXStringViewsToStringSet(x, out.of.limits = out.of.limits,  :
>    'x' has "out of limits" views

You are calling the PairwiseAlignments() constructor function but I
suspect that what you really want is the pairwiseAlignment() function
to perform the alignment. My understanding is that the
PairwiseAlignments() constructor function is a low-level utility
that is almost never intended to be used directly. Always use
pairwiseAlignment():

   > pa <- pairwiseAlignment('ATGC', 'ATAT')

   > writePairwiseAlignments(pa)
   ########################################
   # Program: Biostrings (version 2.33.13), a Bioconductor package
   # Rundate: Tue Jul 15 19:09:10 2014
   ########################################
   #=======================================
   #
   # Aligned_sequences: 2
   # 1: P1
   # 2: S1
   # Matrix: NA
   # Gap_penalty: 14.0
   # Extend_penalty: 4.0
   #
   # Length: 4
   # Identity:       2/4 (50.0%)
   # Similarity:    NA/4 (NA%)
   # Gaps:           0/4 (0.0%)
   # Score: -7.835059
   #
   #
   #=======================================

   P1                 1 ATGC      4
                        ||
   S1                 1 ATAT      4


   #---------------------------------------
   #---------------------------------------

   > mismatchTable(pa)
     PatternId PatternStart PatternEnd PatternSubstring PatternQuality
   1         1            3          3                G              7
   2         1            4          4                C              7
     SubjectStart SubjectEnd SubjectSubstring SubjectQuality
   1            3          3                A              7
   2            4          4                T              7

The documentation could certainly be a little bit clearer about this.

One could fairly argue that the PairwiseAlignments() constructor
should not return invalid PairwiseAlignments objects (which is the
reason why mismatchTable() then choke on it) but this is only one
of the many things that would need some refactoring in the
pairwiseAlignment/PairwiseAlignments code (this code is big and
complex) and unfortunately I can't dedicate time to this at the
moment.

Hope this helps though.

Cheers,
H.

>
> Thanks,
> Watal
>

-- 
Hervé Pagès

Program in Computational Biology
Division of Public Health Sciences
Fred Hutchinson Cancer Research Center
1100 Fairview Ave. N, M1-B514
P.O. Box 19024
Seattle, WA 98109-1024

E-mail: hpages at fhcrc.org
Phone:  (206) 667-5791
Fax:    (206) 667-1319