Where to find asymmetric nucleotide substitution matrix with IUPAC encodings?

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).

---

I am using the pairwiseAlignment function from the Biostrings package to calculate the distance between a consensus sequence and a set of other sequences. Some example code:

# 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)) 

The resulting tibble looks like this:

# A tibble: 5 × 2
  value sequence                                       
  <dbl> <chr>                                          
1  -3.5 ATTGTGATTTGCTTTCAAATTAGTATCTTTAAACCATTGAAAACAGC
2  -0.5 GTTGTGATTTGCTTTCAAATTAGTATCTTTGAACCATTGAAAACAAC
3  -0.5 GTTGTGATTTGCTTTCAAATTAGTATCTTTGAACCATTGAAAACAAC
4  -0.5 GTTGTGATTTGCTTTCAAATTAGTATCTTTGAACCATTGAAAACAAC
5  -0.5 GTTGTGATTTGCTTTCAAATTAGTATCTTTGAACCATTGAAAACAAC

Something isn't quite right. Take a look at the alignment of con against seq[2].

con:    GTTGTGATTTGCTTTCRAATTAGTATCTTTGAACCATTGAAAACAAC
        |||||||||||||||| ||||||||||||||||||||||||||||||
seq[2]: GTTGTGATTTGCTTTCAAATTAGTATCTTTGAACCATTGAAAACAAC 

Note that the only difference is at the position in the consensus sequence where there is ambiguous IUPAC encoding. Because R can represent A or G, the score for seq[2] should be 0 (a perfect match). If the ambiguous base was instead in seq[2] and con was unambiguous, then the score penalty would be appropriate. The discrepancy here is due to the symmetry of the substitution matrix sub:

           A          C          G          T          M          R          W          S          Y          K
A  0.0000000 -1.0000000 -1.0000000 -1.0000000 -0.5000000 -0.5000000 -0.5000000 -1.0000000 -1.0000000 -1.0000000
C -1.0000000  0.0000000 -1.0000000 -1.0000000 -0.5000000 -1.0000000 -1.0000000 -0.5000000 -0.5000000 -1.0000000
G -1.0000000 -1.0000000  0.0000000 -1.0000000 -1.0000000 -0.5000000 -1.0000000 -0.5000000 -1.0000000 -0.5000000
T -1.0000000 -1.0000000 -1.0000000  0.0000000 -1.0000000 -1.0000000 -0.5000000 -1.0000000 -0.5000000 -0.5000000
M -0.5000000 -0.5000000 -1.0000000 -1.0000000 -0.5000000 -0.7500000 -0.7500000 -0.7500000 -0.7500000 -1.0000000
R -0.5000000 -1.0000000 -0.5000000 -1.0000000 -0.7500000 -0.5000000 -0.7500000 -0.7500000 -1.0000000 -0.7500000
W -0.5000000 -1.0000000 -1.0000000 -0.5000000 -0.7500000 -0.7500000 -0.5000000 -1.0000000 -0.7500000 -0.7500000
S -1.0000000 -0.5000000 -0.5000000 -1.0000000 -0.7500000 -0.7500000 -1.0000000 -0.5000000 -0.7500000 -0.7500000
Y -1.0000000 -0.5000000 -1.0000000 -0.5000000 -0.7500000 -1.0000000 -0.7500000 -0.7500000 -0.5000000 -0.7500000
K -1.0000000 -1.0000000 -0.5000000 -0.5000000 -1.0000000 -0.7500000 -0.7500000 -0.7500000 -0.7500000 -0.5000000
V -0.6666667 -0.6666667 -0.6666667 -1.0000000 -0.6666667 -0.6666667 -0.8333333 -0.6666667 -0.8333333 -0.8333333
H -0.6666667 -0.6666667 -1.0000000 -0.6666667 -0.6666667 -0.8333333 -0.6666667 -0.8333333 -0.6666667 -0.8333333
D -0.6666667 -1.0000000 -0.6666667 -0.6666667 -0.8333333 -0.6666667 -0.6666667 -0.8333333 -0.8333333 -0.6666667
B -1.0000000 -0.6666667 -0.6666667 -0.6666667 -0.8333333 -0.8333333 -0.8333333 -0.6666667 -0.6666667 -0.6666667
N -0.7500000 -0.7500000 -0.7500000 -0.7500000 -0.7500000 -0.7500000 -0.7500000 -0.7500000 -0.7500000 -0.7500000
           V          H          D          B     N
A -0.6666667 -0.6666667 -0.6666667 -1.0000000 -0.75
C -0.6666667 -0.6666667 -1.0000000 -0.6666667 -0.75
G -0.6666667 -1.0000000 -0.6666667 -0.6666667 -0.75
T -1.0000000 -0.6666667 -0.6666667 -0.6666667 -0.75
M -0.6666667 -0.6666667 -0.8333333 -0.8333333 -0.75
R -0.6666667 -0.8333333 -0.6666667 -0.8333333 -0.75
W -0.8333333 -0.6666667 -0.6666667 -0.8333333 -0.75
S -0.6666667 -0.8333333 -0.8333333 -0.6666667 -0.75
Y -0.8333333 -0.6666667 -0.8333333 -0.6666667 -0.75
K -0.8333333 -0.8333333 -0.6666667 -0.6666667 -0.75
V -0.6666667 -0.7777778 -0.7777778 -0.7777778 -0.75
H -0.7777778 -0.6666667 -0.7777778 -0.7777778 -0.75
D -0.7777778 -0.7777778 -0.6666667 -0.7777778 -0.75
B -0.7777778 -0.7777778 -0.7777778 -0.6666667 -0.75
N -0.7500000 -0.7500000 -0.7500000 -0.7500000 -0.75

As part of a discussion on Github, I provided the code to manually reassign elements of sub to create an asymmetric substitution matrix:

sub["N", "A"] <- 0
sub["N", "T"] <- 0
sub["N", "G"] <- 0
sub["N", "C"] <- 0
sub["M", "A"] <- 0
sub["M", "C"] <- 0
sub["R", "A"] <- 0
sub["R", "G"] <- 0
sub["W", "A"] <- 0
sub["W", "T"] <- 0
sub["S", "C"] <- 0
sub["S", "G"] <- 0
sub["Y", "C"] <- 0
sub["Y", "T"] <- 0
sub["K", "G"] <- 0
sub["K", "T"] <- 0
sub["V", "A"] <- 0
sub["V", "C"] <- 0
sub["V", "G"] <- 0
sub["H", "A"] <- 0
sub["H", "C"] <- 0
sub["H", "T"] <- 0
sub["D", "A"] <- 0
sub["D", "G"] <- 0
sub["D", "T"] <- 0
sub["B", "C"] <- 0
sub["B", "G"] <- 0
sub["B", "T"] <- 0

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?

Answer 1

Background

I've dug into the source code for the nucleotideSubstitutionMatrix() function from the Biostrings package. Here is the code for the full function from GitHub:

nucleotideSubstitutionMatrix <- function(match = 1, mismatch = 0, baseOnly = FALSE, type = "DNA")
{
  "%safemult%" <- function(x, y) ifelse(is.infinite(x) & y == 0, 0, x * y)
  type <- match.arg(type, c("DNA", "RNA"))
  if (!isSingleNumber(match) || !isSingleNumber(mismatch))
    stop("'match' and 'mismatch' must be non-missing numbers")
  if (baseOnly)
    letters <- IUPAC_CODE_MAP[DNA_BASES]
  else
    letters <- IUPAC_CODE_MAP
  if (type == "RNA")
    names(letters) <- chartr("T", "U", names(letters))
  nLetters <- length(letters)
  splitLetters <- strsplit(letters,split="")
  submat <- matrix(0, nrow = nLetters, ncol = nLetters, dimnames = list(names(letters), names(letters)))
  for(i in 1:nLetters)
    for(j in i:nLetters)
      submat[i,j] <- submat[j,i] <- mean(outer(splitLetters[[i]], splitLetters[[j]], "=="))
  abs(match) * submat - abs(mismatch) %safemult% (1 - submat)
}

Dissecting this, we can see that the matrix is constructed in four main steps.

1. The IUPAC_CODE_MAP named character vector is split into a list of character vectors, where each named element in the list gives an IUPAC character and its constituent DNA bases.

> IUPAC_CODE_MAP
     A      C      G      T      M      R      W      S      Y      K      V      H      D      B      N 
   "A"    "C"    "G"    "T"   "AC"   "AG"   "AT"   "CG"   "CT"   "GT"  "ACG"  "ACT"  "AGT"  "CGT" "ACGT" 
> splitLetters$A
[1] "A"
> splitLetters$V
[1] "A" "C" "G"
> splitLetters$N
[1] "A" "C" "G" "T"

2. submat is initiated as an $n × n$ matrix of zeros, where $n$ is length(names(IUPAC_CODE_MAP)).

> matrix(0, nrow = nLetters, ncol = nLetters, dimnames = list(names(letters), names(letters)))
  A C G T M R W S Y K V H D B N
A 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
C 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
G 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
T 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
M 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
R 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
W 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
S 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Y 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
K 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
V 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
H 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
D 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
B 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
N 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

3. Each element of the matrix is calculated by counting the overlapping elements in the arrays represented by the IUPAC labels given as row and column names.

> for(i in 1:nLetters)
>     for(j in i:nLetters)
>         submat[i,j] <- submat[j,i] <- mean(outer(splitLetters[[i]], splitLetters[[j]], "=="))
> submat
          A         C         G         T         M         R         W         S         Y         K         V         H         D
A 1.0000000 0.0000000 0.0000000 0.0000000 0.5000000 0.5000000 0.5000000 0.0000000 0.0000000 0.0000000 0.3333333 0.3333333 0.3333333
C 0.0000000 1.0000000 0.0000000 0.0000000 0.5000000 0.0000000 0.0000000 0.5000000 0.5000000 0.0000000 0.3333333 0.3333333 0.0000000
G 0.0000000 0.0000000 1.0000000 0.0000000 0.0000000 0.5000000 0.0000000 0.5000000 0.0000000 0.5000000 0.3333333 0.0000000 0.3333333
T 0.0000000 0.0000000 0.0000000 1.0000000 0.0000000 0.0000000 0.5000000 0.0000000 0.5000000 0.5000000 0.0000000 0.3333333 0.3333333
M 0.5000000 0.5000000 0.0000000 0.0000000 0.5000000 0.2500000 0.2500000 0.2500000 0.2500000 0.0000000 0.3333333 0.3333333 0.1666667
R 0.5000000 0.0000000 0.5000000 0.0000000 0.2500000 0.5000000 0.2500000 0.2500000 0.0000000 0.2500000 0.3333333 0.1666667 0.3333333
W 0.5000000 0.0000000 0.0000000 0.5000000 0.2500000 0.2500000 0.5000000 0.0000000 0.2500000 0.2500000 0.1666667 0.3333333 0.3333333
S 0.0000000 0.5000000 0.5000000 0.0000000 0.2500000 0.2500000 0.0000000 0.5000000 0.2500000 0.2500000 0.3333333 0.1666667 0.1666667
Y 0.0000000 0.5000000 0.0000000 0.5000000 0.2500000 0.0000000 0.2500000 0.2500000 0.5000000 0.2500000 0.1666667 0.3333333 0.1666667
K 0.0000000 0.0000000 0.5000000 0.5000000 0.0000000 0.2500000 0.2500000 0.2500000 0.2500000 0.5000000 0.1666667 0.1666667 0.3333333
V 0.3333333 0.3333333 0.3333333 0.0000000 0.3333333 0.3333333 0.1666667 0.3333333 0.1666667 0.1666667 0.3333333 0.2222222 0.2222222
H 0.3333333 0.3333333 0.0000000 0.3333333 0.3333333 0.1666667 0.3333333 0.1666667 0.3333333 0.1666667 0.2222222 0.3333333 0.2222222
D 0.3333333 0.0000000 0.3333333 0.3333333 0.1666667 0.3333333 0.3333333 0.1666667 0.1666667 0.3333333 0.2222222 0.2222222 0.3333333
B 0.0000000 0.3333333 0.3333333 0.3333333 0.1666667 0.1666667 0.1666667 0.3333333 0.3333333 0.3333333 0.2222222 0.2222222 0.2222222
N 0.2500000 0.2500000 0.2500000 0.2500000 0.2500000 0.2500000 0.2500000 0.2500000 0.2500000 0.2500000 0.2500000 0.2500000 0.2500000
          B    N
A 0.0000000 0.25
C 0.3333333 0.25
G 0.3333333 0.25
T 0.3333333 0.25
M 0.1666667 0.25
R 0.1666667 0.25
W 0.1666667 0.25
S 0.3333333 0.25
Y 0.3333333 0.25
K 0.3333333 0.25
V 0.2222222 0.25
H 0.2222222 0.25
D 0.2222222 0.25
B 0.3333333 0.25
N 0.2500000 0.25

4. Finally, the base substitution matrix is converted to a scoring matrix given the user-defined values for match and mismatch.

---

Answer

Concerning step 3, the Biostrings approach asks whether each element of each vector matches each element of the other vector, and then records the proportion of TRUE values in the matrix of results. Note that the code submat[i,j] <- submat[j,i] creates a symmetric matrix $A$. However, even if $[A]_{ij}$ and $[A]_{ji}$ were generated by separate calls, the commutativity of the mean(outer(i, j, "==")) computation means that $A$ will always be symmetric:

> outer(splitLetters$K, splitLetters$V, "==")
      [,1]  [,2]  [,3]
[1,] FALSE FALSE  TRUE
[2,] FALSE FALSE FALSE
> mean(outer(splitLetters$K, splitLetters$V, "=="))
[1] 0.1666667

> outer(splitLetters$V, splitLetters$K, "==")
      [,1]  [,2]
[1,] FALSE FALSE
[2,] FALSE FALSE
[3,]  TRUE FALSE
> mean(outer(splitLetters$V, splitLetters$K, "=="))
[1] 0.1666667

Instead of checking for pairwise equality of elements (==), we can check if a given element from the query is contained in the set of elements in the consensus / subject (%in%).

> outer(splitLetters$K, splitLetters$V, "%in%")
      [,1]  [,2]  [,3]
[1,]  TRUE  TRUE  TRUE
[2,] FALSE FALSE FALSE
> mean(outer(splitLetters$K, splitLetters$V, "%in%"))
[1] 0.5

> outer(splitLetters$V, splitLetters$K, "%in%")
      [,1]  [,2]
[1,] FALSE FALSE
[2,] FALSE FALSE
[3,]  TRUE  TRUE
> mean(outer(splitLetters$V, splitLetters$K, "%in%"))
[1] 0.3333333

This is not a commutative operation, so the resulting substitution matrix is asymmetric.

asym_NSM <- function(match = 1, mismatch = 0, baseOnly = FALSE, type = "DNA")
{
  "%safemult%" <- function(x, y) ifelse(is.infinite(x) & y == 0, 0, x * y)
  type <- match.arg(type, c("DNA", "RNA"))
  if (!isSingleNumber(match) || !isSingleNumber(mismatch))
    stop("'match' and 'mismatch' must be non-missing numbers")
  if (baseOnly)
    letters <- IUPAC_CODE_MAP[DNA_BASES]
  else
    letters <- IUPAC_CODE_MAP
  if (type == "RNA")
    names(letters) <- chartr("T", "U", names(letters))
  nLetters <- length(letters)
  splitLetters <- strsplit(letters,split="")
  submat <- matrix(0, nrow = nLetters, ncol = nLetters, dimnames = list(names(letters), names(letters)))
  for(i in 1:nLetters) {
    for(j in i:nLetters) {
      submat[i,j] <- mean(outer(splitLetters[[i]], splitLetters[[j]], "%in%"))
      submat[j,i] <- mean(outer(splitLetters[[j]], splitLetters[[i]], "%in%"))
    }
  }
  
  abs(match) * submat - abs(mismatch) %safemult% (1 - submat)

}

> asym_NSM()
          A         C         G         T         M         R         W         S         Y         K         V         H         D
A 1.0000000 0.0000000 0.0000000 0.0000000 1.0000000 1.0000000 1.0000000 0.0000000 0.0000000 0.0000000 1.0000000 1.0000000 1.0000000
C 0.0000000 1.0000000 0.0000000 0.0000000 1.0000000 0.0000000 0.0000000 1.0000000 1.0000000 0.0000000 1.0000000 1.0000000 0.0000000
G 0.0000000 0.0000000 1.0000000 0.0000000 0.0000000 1.0000000 0.0000000 1.0000000 0.0000000 1.0000000 1.0000000 0.0000000 1.0000000
T 0.0000000 0.0000000 0.0000000 1.0000000 0.0000000 0.0000000 1.0000000 0.0000000 1.0000000 1.0000000 0.0000000 1.0000000 1.0000000
M 0.5000000 0.5000000 0.0000000 0.0000000 1.0000000 0.5000000 0.5000000 0.5000000 0.5000000 0.0000000 1.0000000 1.0000000 0.5000000
R 0.5000000 0.0000000 0.5000000 0.0000000 0.5000000 1.0000000 0.5000000 0.5000000 0.0000000 0.5000000 1.0000000 0.5000000 1.0000000
W 0.5000000 0.0000000 0.0000000 0.5000000 0.5000000 0.5000000 1.0000000 0.0000000 0.5000000 0.5000000 0.5000000 1.0000000 1.0000000
S 0.0000000 0.5000000 0.5000000 0.0000000 0.5000000 0.5000000 0.0000000 1.0000000 0.5000000 0.5000000 1.0000000 0.5000000 0.5000000
Y 0.0000000 0.5000000 0.0000000 0.5000000 0.5000000 0.0000000 0.5000000 0.5000000 1.0000000 0.5000000 0.5000000 1.0000000 0.5000000
K 0.0000000 0.0000000 0.5000000 0.5000000 0.0000000 0.5000000 0.5000000 0.5000000 0.5000000 1.0000000 0.5000000 0.5000000 1.0000000
V 0.3333333 0.3333333 0.3333333 0.0000000 0.6666667 0.6666667 0.3333333 0.6666667 0.3333333 0.3333333 1.0000000 0.6666667 0.6666667
H 0.3333333 0.3333333 0.0000000 0.3333333 0.6666667 0.3333333 0.6666667 0.3333333 0.6666667 0.3333333 0.6666667 1.0000000 0.6666667
D 0.3333333 0.0000000 0.3333333 0.3333333 0.3333333 0.6666667 0.6666667 0.3333333 0.3333333 0.6666667 0.6666667 0.6666667 1.0000000
B 0.0000000 0.3333333 0.3333333 0.3333333 0.3333333 0.3333333 0.3333333 0.6666667 0.6666667 0.6666667 0.6666667 0.6666667 0.6666667
N 0.2500000 0.2500000 0.2500000 0.2500000 0.5000000 0.5000000 0.5000000 0.5000000 0.5000000 0.5000000 0.7500000 0.7500000 0.7500000
          B N
A 0.0000000 1
C 1.0000000 1
G 1.0000000 1
T 1.0000000 1
M 0.5000000 1
R 0.5000000 1
W 0.5000000 1
S 1.0000000 1
Y 1.0000000 1
K 1.0000000 1
V 0.6666667 1
H 0.6666667 1
D 0.6666667 1
B 1.0000000 1
N 0.7500000 1

Answer 2

Additional answer on asymetric IUPAC score matrix.

To get the scores with an asymentric score matrix. This is also implemented in the Bioconductor package MSA2dist (https://www.bioconductor.org/packages/release/bioc/html/MSA2dist.html).

library(Biostrings)
library(MSA2dist)

# consensus sequence and subject sequences
seqs <- setNames(DNAStringSet(c("GTTGTGATTTGCTTTCRAATTAGTATCTTTGAACCATTGAAAACAAC",
"GTTGTGATTTGCTTTCAAATTAGTATCTTTGAACCATTGAAAACAAC")), c("con", "seq"))

# define score matrix
scoreMatrix <- MSA2dist::iupacMatrix()
scoreMatrix["A", "R"]
# [1] 0.5
scoreMatrix["R", "A"] <- 1

MSA2dist::dnastring2dist(seqs, score=scoreMatrix)
MSA2dist::dnastring2dist(rev(seqs), score=scoreMatrix)

# or here using the asymetric score matrix from Biostrings

ansm <- nucleotideSubstitutionMatrix(symmetric = FALSE)

#note: MSA2dist score is normalized by analysed sites
MSA2dist::dnastring2dist(seqs, score=ansm)
MSA2dist::dnastring2dist(rev(seqs), score=ansm)

pairwiseAlignment(seqs[[1]], seqs[[2]], substitutionMatrix = ansm)
pairwiseAlignment(seqs[[2]], seqs[[1]], substitutionMatrix = ansm)

Since MSA2dist::dnastring2dist was until now calculating all possible pairwise combinations from top to bottom (leaving out the already seen combinations), now I will include the option to really calculate all pairwise combinations to take into account asymetric scores and sequence position.