Problem with pairwiseAlignment in R

Question

[l.roca]

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

[Devon Ryan]

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

[l.roca]

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

Match score = +1.

Mismatch score = -1.

Gap score = -2.

Answer 1

[Devon Ryan]

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.

[l.roca]

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

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

GGAATGG
--A-TG-

[poisonAlien]

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

[Devon Ryan]

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