This is a test version of Biostars. For the public version, visit https://www.biostars.org.
R/Biostring: Print Sequence Alignment To File

I do a simple pairwise DNA sequence alignment with pairwiseAlignment from the Biostrings package in Bioconductor:

library('Biostrings')
seq1 = 'ATGCTA'
seq2 = 'ATGTA'
pairwiseAlignment(pattern = seq1, subject = seq2)

The output looks as follows:

Global PairwiseAlignedFixedSubject (1 of 1)
pattern: [1] ATGCTA 
subject: [1] ATG-TA 
score: -4.091219

For very long sequences, the output is truncated and only one line is shown:

Global PairwiseAlignedFixedSubject (1 of 1)
pattern: [1] AT-----------------------------...----------------TGTCTTCCAKATCTGGCGCGCCTGGGTTGATATC 
subject: [1] ATTGGCGGCCGCGCCACCATGCCAGAGCCAG...GAAGGCTGTATGCTGTTGTCTTCAAGATCTGGTACCGCTGGGTTGATATC 
score: -29418.8

How can I output the complete alignment to a text file?

bioconductor sequence alignment
• 6,578 views • link
updated by Steve Lianoglou Steve Lianoglou • written by Emempe • Emempe

Looking briefly to your long alignment shows an AT base pair much closer to the rest of the sequence. In fact alignment start at 1 for each sequence. Question to anyone: how to avoid such a suboptimal alignment ?

• 0 views • link
written by Mycroft34 • Mycroft34

1 answer

You can get to the "text representation" of the pairwise alignment by using the pattern and subject accessors, then convert these to normal character object that you can write "as usual" eg:

library('Biostrings')
seq1 <- 'ATGCTA'
seq2 <- 'ATGTA'
pa <- pairwiseAlignment(pattern = seq1, subject = seq2)

as.character(pattern(pa))
[1] "ATGCTA"

as.character(subject(pa))
[1] "ATG-TA"

Not sure how you want to write it out, but there's (at least) a start.

• 0 views • link
written by Steve Lianoglou Steve Lianoglou

Thanks! I think I have to write an output function myself ... one would have to put it in blocks of 60, 80 or whatever nucleotides.

• 0 views • link
written by Emempe • Emempe

Log in to answer this question.

More posts like this