Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-08-26T14:40:42.830Z
Answers.LastActivityDate: 2022-08-26T14:40:42.830Z
Answers.UserDisplayName: acvill
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-05-21T21:18:28.477Z
Answers.LastActivityDate: 2023-05-21T21:19:57.427Z
Answers.LastEditDate: 2023-05-21T21:19:57.427Z
Answers.UserDisplayName: kullrich
Body: <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>
  
Tags: <r><phylogenetics><bioconductor>
Title: Where to find asymmetric nucleotide substitution matrix with IUPAC encodings?
ThreadText: # Where to find asymmetric nucleotide substitution matrix with IUPAC encodings?
  
  (Asked by: username_0 on 2022-04-11)
  
  **Update:** I submitted a [pull request](https://github.com/Bioconductor/Biostrings/pull/77) 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](https://en.wikipedia.org/wiki/Nucleic_acid_notation#IUPAC_notation). 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](https://github.com/Bioconductor/Biostrings/issues/3#issuecomment-1095143260), 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](https://github.com/Bioconductor/Biostrings/blob/4f9031a1bc0e102148060cd39aa441d812d75ce7/R/pairwiseAlignment.R#L16):
  
  ```
  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.
  
  
  ---
  
  