This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Problem with pairwiseAlignment in R

Hi,

Based on two sequences from the book: Sequence Alignment Methods, Models, Concepts,and Strategies Edited by Michael S. Rosenberg.

I run the following script:

> library(Biostrings)
> s1 <- "GGAATGG"
> s2 <- "ATG"
> sigma <- nucleotideSubstitutionMatrix(match = 2, mismatch = -1, baseOnly = TRUE)
> globalAligns1s2 <- pairwiseAlignment(s1, s2, type="global", substitutionMatrix = sigma,  scoreOnly = FALSE)
> globalAligns1s2 # Print out the optimal alignment and its score

With the following output:

Global PairwiseAlignmentsSingleSubject (1 of 1)
pattern: [1] GGAATGG
subject: [1] AT----G
score: -26

Based on the book results I should get an alignment score of -5 and one of these alignment:

GGAATGG
---ATG-

Or

GGAATGG
---AT-G

Or

GGAATGG
--A-TG-

Or

GGAATGG
--A-T-G

Any help or idea what I did wrong?

Thanks

alignment r
• 4,882 views • link
updated by Ram Ram • written by l.roca • l.roca

Try tweaking the gapOpening and gapExtension parameters to be lower values, such as -2 and -1.

• 0 views • link
updated by Ram Ram • written by Devon Ryan Devon Ryan

I tried without success. The parameters in the book where:

Match score = +1.

Mismatch score = -1.

Gap score = -2.

• 0 views • link
updated by Ram Ram • written by l.roca • l.roca

1 answer 

library(Biostrings)
sm <- nucleotideSubstitutionMatrix(match = 1, mismatch = -1, baseOnly = TRUE)
s1 <- DNAString("GGAATGG")
s2 <- DNAString("ATG")
pairwiseAlignment(s1,s2, substitutionMatrix=sm, scoreOnly=F, gapOpening=0, gapExtension=-2)

I get one of the expected results and a score of -5. Note that pairwiseAlignment() takes a bit of playing with to see how its scoring really works.

• 0 views • link
updated by Ram Ram • written by Devon Ryan Devon Ryan

Thanks a lot, I am getting it (slowley :-)).

How can I output it so it will look like this:

GGAATGG
--A-TG-
• 0 views • link
written by l.roca • l.roca

No clue, I honestly don't use pairwiseAlignment() for anything real.

• 0 views • link
written by Devon Ryan Devon Ryan

use writePairwiseAlignments() to write the alignment object to a file. It would like you have shown above.

• 0 views • link
written by poisonAlien poisonAlien

Log in to answer this question.

More posts like this