Return to Question

Update: I submitted a pull request to the Biostrings repo. The functionality I describe in my question and answer can now be implemented with nucleotideSubstitutionMatrix(symmetric = FALSE).

# R v4.2.0
library(Biostrings)   # v2.63.3
library(BiocGenerics) # v0.41.2
library(tidyverse)    # v1.3.1

# consensus sequence
con <- DNAStringSet(c("GTTGTGATTTGCTTTCRAATTAGTATCTTTGAACCATTGAAAACAAC"))

# subject sequences
seq <- DNAStringSet(c("ATTGTGATTTGCTTTCAAATTAGTATCTTTAAACCATTGAAAACAGC",
                      "GTTGTGATTTGCTTTCAAATTAGTATCTTTGAACCATTGAAAACAAC",
                      "GTTGTGATTTGCTTTCAAATTAGTATCTTTGAACCATTGAAAACAAC",
                      "GTTGTGATTTGCTTTCAAATTAGTATCTTTGAACCATTGAAAACAAC",
                      "GTTGTGATTTGCTTTCAAATTAGTATCTTTGAACCATTGAAAACAAC"))

# substitution matrix
sub <- nucleotideSubstitutionMatrix(match = 0, mismatch = 1)

# get alignments
sapply(X = 1:length(seq),
       FUN = function(i) {
         pairwiseAlignment(pattern = con,
                           subject = seq[i],
                           substitutionMatrix = sub) %>%
         score
      }
    ) %>% as_tibble %>% mutate(seq = as.character(seq)) 

# R v4.2.0
library(Biostrings)   # v2.63.3
library(BiocGenerics) # v0.41.2
library(tidyverse)    # v1.3.1

# consensus sequence
con <- DNAStringSet(c("GTTGTGATTTGCTTTCRAATTAGTATCTTTGAACCATTGAAAACAAC"))

# subject sequences
seq <- DNAStringSet(c("ATTGTGATTTGCTTTCAAATTAGTATCTTTAAACCATTGAAAACAGC",
                      "GTTGTGATTTGCTTTCAAATTAGTATCTTTGAACCATTGAAAACAAC",
                      "GTTGTGATTTGCTTTCAAATTAGTATCTTTGAACCATTGAAAACAAC",
                      "GTTGTGATTTGCTTTCAAATTAGTATCTTTGAACCATTGAAAACAAC",
                      "GTTGTGATTTGCTTTCAAATTAGTATCTTTGAACCATTGAAAACAAC"))

# substitution matrix
sub <- nucleotideSubstitutionMatrix(match = 0, mismatch = 1)

# get alignments
sapply(X = 1:length(seq),
       FUN = function(i) {
         pairwiseAlignment(pattern = con,
                           subject = seq[i],
                           substitutionMatrix = sub) %>%
         score
      }
    ) %>% as_tibble %>% mutate(sequence = as.character(seq)) 

This covers all cases where the pattern (con) contains ambiguous bases and the subjects (seq) are unambiguous, changing the matrix so the subjects are not penalized. However, these edits do not account for cases where both the pattern and the subject have ambiguous encodings at the same position; e.g. ["M","V"] and ["V","M"] have the same value associated with them, even though V is a superset of M. There are also more nuanced cases where the set of bases represented by a pair of ambiguous encodings overlap incompletely in both directions, as with ["H","D"] and ["D","H"]. Given enough time, I could probably figure out the appropriate scores for each case to complete the asymmetric substitution matrix, but I'm surely reinventing the wheel. Can someone supply, as a standalone table or as part of a codebase, an asymmetric nucleotide substitution matrix that satisfies the criteria outlined herein?

This covers all cases where the pattern (con) contains ambiguous bases and the subjects (seq) are unambiguous, changing the matrix so the subjects are not penalized. However, these edits do not account for cases where both the pattern and the subject have ambiguous encodings at the same position; e.g. ["M","V"] and ["V","M"] have the same value associated with them, even though V (A || C || G) is a superset of M (A || C). There are also more nuanced cases where the set of bases represented by a pair of ambiguous encodings overlap incompletely in both directions, as with H (A || C || T) and D (A || G || T). Given enough time, I could probably figure out the appropriate scores for each case to complete the asymmetric substitution matrix, but I'm surely reinventing the wheel. Can someone supply, as a standalone table or as part of a codebase, an asymmetric nucleotide substitution matrix that satisfies the criteria outlined herein?