Putative bug in Biostrings
1
0
Entering edit mode
@herve-pages-1542
Last seen 2 hours ago
Seattle, WA, United States
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
Alignment Cancer Biostrings Alignment Cancer Biostrings • 1.2k views
ADD COMMENT
0
Entering edit mode
@watal-m-iwasaki-6651
Last seen 9.0 years ago
Japan
Thanks Herv?, Using pairwiseAlignment() instead of PairwiseAlignments() solved the problem. Sorry that I didn't read the documentation thoroughly and bothered you. Aside from refactoring in those functions and the documentation, it could be reasonable to exclude such low-level (and confusable) utility from NAMESPACE. Anyway, thank you again for your help. Best, Watal On Wed, Jul 16, 2014 at 11:17 AM, Hervé Pagès <hpages at="" fhcrc.org=""> wrote: > 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 -- Watal M. Iwasaki / ?? ? Graduate University for Advanced Studies, Hayama, Kanagawa 240-0193, Japan +81-46-858-1576 heavy.watal at gmail.com http://meme.biology.tohoku.ac.jp/students/iwasaki/
ADD COMMENT

Login before adding your answer.

Traffic: 558 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6