Sequence alignment of extended IUPAC sequences

Question

I want to align a short sequence to a reference, and both are encoded in the extended IUPAC notation (i.e. 16 letters).

What alignment software can I use for that?

For example, in Bowtie "all non-A/C/G/T characters in reads become Ns", so it ignores which extended letters are matching.

Answer 1

You can try using Biostrings, depends on how large is your reference, if it is small, you can just pairwiseAlignment with the nucleotideSubstitutionMatrix which supposedly :

library(Biostrings)
data(phiX174Phage)

x = DNAString("TCCAWRSCGCAGAAGT")

lapply(phiX174Phage,function(y)pairwiseAlignment(x,y,
substitutionMatrix = nucleotideSubstitutionMatrix()))

$Genbank
Global PairwiseAlignmentsSingleSubject (1 of 1)
pattern: -------------TCCAWRSCGCAGAAGT...----------------------------
subject: GAGTTTTATCGCTTCCATGACGCAGAAGT...TAAAAATGATTGGCGTATCCAACCTGCA
score: -21491.92 

$RF70s
Global PairwiseAlignmentsSingleSubject (1 of 1)
pattern: -------------TCCAWRSCGCAGAAGT...----------------------------
subject: GAGTTTTATCGCTTCCATGACGCAGAAGT...TAAAAATGATTGGCGTATCCAACCTGCA
score: -21491.92 

$SS78
Global PairwiseAlignmentsSingleSubject (1 of 1)
pattern: -------------TCCAWRSCGCAGAAGT...----------------------------
subject: GAGTTTTATCGCTTCCATGACGCAGAAGT...TAAAAATGATTGGCGTATCCAACCTGCA
score: -21491.92 

$Bull
Global PairwiseAlignmentsSingleSubject (1 of 1)
pattern: -------------TCCAWRSCGCAGAAGT...----------------------------
subject: GAGTTTTATCGCTTCCATGACGCAGAAGT...TAAAAATGATTGGCGTATCCAACCTGCA
score: -21491.92 

$G97
Global PairwiseAlignmentsSingleSubject (1 of 1)
pattern: -------------TCCAWRSCGCAGAAGT...----------------------------
subject: GAGTTTTATCGCTTCCATGACGCAGAAGT...TAAAAATGATTGGCGTATCCAACCTGCA
score: -21491.92 

$NEB03
Global PairwiseAlignmentsSingleSubject (1 of 1)
pattern: -------------TCCAWRSCGCAGAAGT...----------------------------
subject: GAGTTTTATCGCTTCCATGACGCAGAAGT...TAAAAATGATTGGCGTATCCAACCTGCA
score: -21491.92 

If it is much larger, you can convert the sequence of interest to a pwm and find it in the reference fasta, like shown in this example