<h2>Background</h2>
<p>I've dug into the source code for the <code>nucleotideSubstitutionMatrix()</code> function from the <code>Biostrings</code> package. Here is the code for the full function <a href="https://github.com/Bioconductor/Biostrings/blob/4f9031a1bc0e102148060cd39aa441d812d75ce7/R/pairwiseAlignment.R#L16" rel="nofollow noreferrer">from GitHub</a>:</p>
<pre><code>nucleotideSubstitutionMatrix &lt;- function(match = 1, mismatch = 0, baseOnly = FALSE, type = &quot;DNA&quot;)
{
  &quot;%safemult%&quot; &lt;- function(x, y) ifelse(is.infinite(x) &amp; y == 0, 0, x * y)
  type &lt;- match.arg(type, c(&quot;DNA&quot;, &quot;RNA&quot;))
  if (!isSingleNumber(match) || !isSingleNumber(mismatch))
    stop(&quot;'match' and 'mismatch' must be non-missing numbers&quot;)
  if (baseOnly)
    letters &lt;- IUPAC_CODE_MAP[DNA_BASES]
  else
    letters &lt;- IUPAC_CODE_MAP
  if (type == &quot;RNA&quot;)
    names(letters) &lt;- chartr(&quot;T&quot;, &quot;U&quot;, names(letters))
  nLetters &lt;- length(letters)
  splitLetters &lt;- strsplit(letters,split=&quot;&quot;)
  submat &lt;- 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] &lt;- submat[j,i] &lt;- mean(outer(splitLetters[[i]], splitLetters[[j]], &quot;==&quot;))
  abs(match) * submat - abs(mismatch) %safemult% (1 - submat)
}
</code></pre>
<p>Dissecting this, we can see that the matrix is constructed in four main steps.</p>
<p><strong>1</strong>. The <code>IUPAC_CODE_MAP</code> 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.</p>
<pre><code>&gt; IUPAC_CODE_MAP
     A      C      G      T      M      R      W      S      Y      K      V      H      D      B      N 
   &quot;A&quot;    &quot;C&quot;    &quot;G&quot;    &quot;T&quot;   &quot;AC&quot;   &quot;AG&quot;   &quot;AT&quot;   &quot;CG&quot;   &quot;CT&quot;   &quot;GT&quot;  &quot;ACG&quot;  &quot;ACT&quot;  &quot;AGT&quot;  &quot;CGT&quot; &quot;ACGT&quot; 
&gt; splitLetters$A
[1] &quot;A&quot;
&gt; splitLetters$V
[1] &quot;A&quot; &quot;C&quot; &quot;G&quot;
&gt; splitLetters$N
[1] &quot;A&quot; &quot;C&quot; &quot;G&quot; &quot;T&quot;
</code></pre>
<p><strong>2</strong>. <code>submat</code> is initiated as an <span class="math-container">$n × n$</span> matrix of zeros, where <span class="math-container">$n$</span> is <code>length(names(IUPAC_CODE_MAP))</code>.</p>
<pre><code>&gt; 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
</code></pre>
<p><strong>3</strong>. 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.</p>
<pre><code>&gt; for(i in 1:nLetters)
&gt;     for(j in i:nLetters)
&gt;         submat[i,j] &lt;- submat[j,i] &lt;- mean(outer(splitLetters[[i]], splitLetters[[j]], &quot;==&quot;))
&gt; 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
</code></pre>
<p><strong>4</strong>. Finally, the base substitution matrix is converted to a scoring matrix given the user-defined values for <code>match</code> and <code>mismatch</code>.</p>
<hr />
<h2>Answer</h2>
<p>Concerning step <strong>3</strong>, the <code>Biostrings</code> approach asks whether each element of each vector matches each element of the other vector, and then records the proportion of <code>TRUE</code> values in the  matrix of results. Note that the code <code>submat[i,j] &lt;- submat[j,i]</code> creates a symmetric matrix <span class="math-container">$A$</span>. However, even if <span class="math-container">$[A]_{ij}$</span> and <span class="math-container">$[A]_{ji}$</span> were generated by separate calls, the commutativity of the <code>mean(outer(i, j, &quot;==&quot;))</code> computation means that <span class="math-container">$A$</span> will always be symmetric:</p>
<pre><code>&gt; outer(splitLetters<span class="math-container">$K, splitLetters$</span>V, &quot;==&quot;)
      [,1]  [,2]  [,3]
[1,] FALSE FALSE  TRUE
[2,] FALSE FALSE FALSE
&gt; mean(outer(splitLetters<span class="math-container">$K, splitLetters$</span>V, &quot;==&quot;))
[1] 0.1666667

&gt; outer(splitLetters<span class="math-container">$V, splitLetters$</span>K, &quot;==&quot;)
      [,1]  [,2]
[1,] FALSE FALSE
[2,] FALSE FALSE
[3,]  TRUE FALSE
&gt; mean(outer(splitLetters<span class="math-container">$V, splitLetters$</span>K, &quot;==&quot;))
[1] 0.1666667
</code></pre>
<p>Instead of checking for pairwise equality of elements (<code>==</code>), we can check if a given element from the query is contained in the set of elements in the consensus / subject (<code>%in%</code>).</p>
<pre><code>&gt; outer(splitLetters<span class="math-container">$K, splitLetters$</span>V, &quot;%in%&quot;)
      [,1]  [,2]  [,3]
[1,]  TRUE  TRUE  TRUE
[2,] FALSE FALSE FALSE
&gt; mean(outer(splitLetters<span class="math-container">$K, splitLetters$</span>V, &quot;%in%&quot;))
[1] 0.5

&gt; outer(splitLetters<span class="math-container">$V, splitLetters$</span>K, &quot;%in%&quot;)
      [,1]  [,2]
[1,] FALSE FALSE
[2,] FALSE FALSE
[3,]  TRUE  TRUE
&gt; mean(outer(splitLetters<span class="math-container">$V, splitLetters$</span>K, &quot;%in%&quot;))
[1] 0.3333333
</code></pre>
<p>This is not a commutative operation, so the resulting substitution matrix is asymmetric.</p>
<pre><code>asym_NSM &lt;- function(match = 1, mismatch = 0, baseOnly = FALSE, type = &quot;DNA&quot;)
{
  &quot;%safemult%&quot; &lt;- function(x, y) ifelse(is.infinite(x) &amp; y == 0, 0, x * y)
  type &lt;- match.arg(type, c(&quot;DNA&quot;, &quot;RNA&quot;))
  if (!isSingleNumber(match) || !isSingleNumber(mismatch))
    stop(&quot;'match' and 'mismatch' must be non-missing numbers&quot;)
  if (baseOnly)
    letters &lt;- IUPAC_CODE_MAP[DNA_BASES]
  else
    letters &lt;- IUPAC_CODE_MAP
  if (type == &quot;RNA&quot;)
    names(letters) &lt;- chartr(&quot;T&quot;, &quot;U&quot;, names(letters))
  nLetters &lt;- length(letters)
  splitLetters &lt;- strsplit(letters,split=&quot;&quot;)
  submat &lt;- 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] &lt;- mean(outer(splitLetters[[i]], splitLetters[[j]], &quot;%in%&quot;))
      submat[j,i] &lt;- mean(outer(splitLetters[[j]], splitLetters[[i]], &quot;%in%&quot;))
    }
  }
  
  abs(match) * submat - abs(mismatch) %safemult% (1 - submat)

}
</code></pre>
<pre><code>&gt; 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
</code></pre>

CC BY-SA 4.0
2022-08-26T14:40:42.830Z
2022-08-26T14:40:42.830Z
acvill
<p>Additional answer on asymetric IUPAC score matrix.</p>
<p>To get the scores with an asymentric score matrix. This is also implemented in the Bioconductor package MSA2dist (<a href="https://www.bioconductor.org/packages/release/bioc/html/MSA2dist.html" rel="nofollow noreferrer">https://www.bioconductor.org/packages/release/bioc/html/MSA2dist.html</a>).</p>
<pre><code>library(Biostrings)
library(MSA2dist)

# consensus sequence and subject sequences
seqs &lt;- setNames(DNAStringSet(c(&quot;GTTGTGATTTGCTTTCRAATTAGTATCTTTGAACCATTGAAAACAAC&quot;,
&quot;GTTGTGATTTGCTTTCAAATTAGTATCTTTGAACCATTGAAAACAAC&quot;)), c(&quot;con&quot;, &quot;seq&quot;))

# define score matrix
scoreMatrix &lt;- MSA2dist::iupacMatrix()
scoreMatrix[&quot;A&quot;, &quot;R&quot;]
# [1] 0.5
scoreMatrix[&quot;R&quot;, &quot;A&quot;] &lt;- 1

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

# or here using the asymetric score matrix from Biostrings

ansm &lt;- 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)
</code></pre>
<p>Since <code>MSA2dist::dnastring2dist</code> 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.</p>

CC BY-SA 4.0
2023-05-21T21:18:28.477Z
2023-05-21T21:19:57.427Z
2023-05-21T21:19:57.427Z
kullrich
<p><strong>Update:</strong> I submitted a <a href="https://github.com/Bioconductor/Biostrings/pull/77" rel="nofollow noreferrer">pull request</a> to the Biostrings repo. The functionality I describe in my question and answer can now be implemented with <code>nucleotideSubstitutionMatrix(symmetric = FALSE)</code>.</p>
<hr />
<p>I am using the <code>pairwiseAlignment</code> function from the <code>Biostrings</code> package to calculate the distance between a consensus sequence and a set of other sequences. Some example code:</p>
<pre><code># R v4.2.0
library(Biostrings)   # v2.63.3
library(BiocGenerics) # v0.41.2
library(tidyverse)    # v1.3.1

# consensus sequence
con &lt;- DNAStringSet(c(&quot;GTTGTGATTTGCTTTCRAATTAGTATCTTTGAACCATTGAAAACAAC&quot;))

# subject sequences
seq &lt;- DNAStringSet(c(&quot;ATTGTGATTTGCTTTCAAATTAGTATCTTTAAACCATTGAAAACAGC&quot;,
                      &quot;GTTGTGATTTGCTTTCAAATTAGTATCTTTGAACCATTGAAAACAAC&quot;,
                      &quot;GTTGTGATTTGCTTTCAAATTAGTATCTTTGAACCATTGAAAACAAC&quot;,
                      &quot;GTTGTGATTTGCTTTCAAATTAGTATCTTTGAACCATTGAAAACAAC&quot;,
                      &quot;GTTGTGATTTGCTTTCAAATTAGTATCTTTGAACCATTGAAAACAAC&quot;))

# substitution matrix
sub &lt;- nucleotideSubstitutionMatrix(match = 0, mismatch = 1)

# get alignments
sapply(X = 1:length(seq),
       FUN = function(i) {
         pairwiseAlignment(pattern = con,
                           subject = seq[i],
                           substitutionMatrix = sub) %&gt;%
         score
      }
    ) %&gt;% as_tibble %&gt;% mutate(sequence = as.character(seq)) 
</code></pre>
<p>The resulting tibble looks like this:</p>
<pre><code># A tibble: 5 × 2
  value sequence                                       
  &lt;dbl&gt; &lt;chr&gt;                                          
1  -3.5 ATTGTGATTTGCTTTCAAATTAGTATCTTTAAACCATTGAAAACAGC
2  -0.5 GTTGTGATTTGCTTTCAAATTAGTATCTTTGAACCATTGAAAACAAC
3  -0.5 GTTGTGATTTGCTTTCAAATTAGTATCTTTGAACCATTGAAAACAAC
4  -0.5 GTTGTGATTTGCTTTCAAATTAGTATCTTTGAACCATTGAAAACAAC
5  -0.5 GTTGTGATTTGCTTTCAAATTAGTATCTTTGAACCATTGAAAACAAC
</code></pre>
<p>Something isn't quite right. Take a look at the alignment of <code>con</code> against <code>seq[2]</code>.</p>
<pre><code>con:    GTTGTGATTTGCTTTCRAATTAGTATCTTTGAACCATTGAAAACAAC
        |||||||||||||||| ||||||||||||||||||||||||||||||
seq[2]: GTTGTGATTTGCTTTCAAATTAGTATCTTTGAACCATTGAAAACAAC 
</code></pre>
<p>Note that the only difference is at the position in the consensus sequence where there is <a href="https://en.wikipedia.org/wiki/Nucleic_acid_notation#IUPAC_notation" rel="nofollow noreferrer">ambiguous IUPAC encoding</a>. Because <strong>R</strong> can represent <strong>A</strong> or <strong>G</strong>, the score for <code>seq[2]</code> should be 0 (a perfect match). If the ambiguous base was instead in <code>seq[2]</code> and <code>con</code> was unambiguous, then the score penalty would be appropriate. The discrepancy here is due to the symmetry of the substitution matrix <code>sub</code>:</p>
<pre><code>           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
</code></pre>
<p><a href="https://github.com/Bioconductor/Biostrings/issues/3#issuecomment-1095143260" rel="nofollow noreferrer">As part of a discussion on Github</a>, I provided the code to manually reassign elements of <code>sub</code> to create an asymmetric substitution matrix:</p>
<pre><code>sub[&quot;N&quot;, &quot;A&quot;] &lt;- 0
sub[&quot;N&quot;, &quot;T&quot;] &lt;- 0
sub[&quot;N&quot;, &quot;G&quot;] &lt;- 0
sub[&quot;N&quot;, &quot;C&quot;] &lt;- 0
sub[&quot;M&quot;, &quot;A&quot;] &lt;- 0
sub[&quot;M&quot;, &quot;C&quot;] &lt;- 0
sub[&quot;R&quot;, &quot;A&quot;] &lt;- 0
sub[&quot;R&quot;, &quot;G&quot;] &lt;- 0
sub[&quot;W&quot;, &quot;A&quot;] &lt;- 0
sub[&quot;W&quot;, &quot;T&quot;] &lt;- 0
sub[&quot;S&quot;, &quot;C&quot;] &lt;- 0
sub[&quot;S&quot;, &quot;G&quot;] &lt;- 0
sub[&quot;Y&quot;, &quot;C&quot;] &lt;- 0
sub[&quot;Y&quot;, &quot;T&quot;] &lt;- 0
sub[&quot;K&quot;, &quot;G&quot;] &lt;- 0
sub[&quot;K&quot;, &quot;T&quot;] &lt;- 0
sub[&quot;V&quot;, &quot;A&quot;] &lt;- 0
sub[&quot;V&quot;, &quot;C&quot;] &lt;- 0
sub[&quot;V&quot;, &quot;G&quot;] &lt;- 0
sub[&quot;H&quot;, &quot;A&quot;] &lt;- 0
sub[&quot;H&quot;, &quot;C&quot;] &lt;- 0
sub[&quot;H&quot;, &quot;T&quot;] &lt;- 0
sub[&quot;D&quot;, &quot;A&quot;] &lt;- 0
sub[&quot;D&quot;, &quot;G&quot;] &lt;- 0
sub[&quot;D&quot;, &quot;T&quot;] &lt;- 0
sub[&quot;B&quot;, &quot;C&quot;] &lt;- 0
sub[&quot;B&quot;, &quot;G&quot;] &lt;- 0
sub[&quot;B&quot;, &quot;T&quot;] &lt;- 0
</code></pre>
<p>This covers all cases where the pattern (<code>con</code>) contains ambiguous bases and the subjects (<code>seq</code>) 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; <em>e.g.</em> <code>[&quot;M&quot;,&quot;V&quot;]</code> and <code>[&quot;V&quot;,&quot;M&quot;]</code> have the same value associated with them, even though <strong>V</strong> (<code>A || C || G</code>) is a superset of <strong>M</strong> (<code>A || C</code>). 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 <strong>H</strong> (<code>A || C || T</code>) and <strong>D</strong> (<code>A || G || T</code>). 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. <strong>Can someone supply, as a standalone table or as part of a codebase, an asymmetric nucleotide substitution matrix that satisfies the criteria outlined herein?</strong></p>

<r><phylogenetics><bioconductor>
Where to find asymmetric nucleotide substitution matrix with IUPAC encodings?
# Where to find asymmetric nucleotide substitution matrix with IUPAC encodings?

(Asked by: username_0 on 2022-04-11)

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                                       
                                             
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?

**Tags:** r, phylogenetics, bioconductor

---

## 2 answers

### Answer by username_0 on 2022-08-26. Score: 1 | Accepted answer

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 by username_1 on 2023-05-21. Score: 1

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.


---

