<p>The technique is phylogenetics is called likelihood ratio and for a regular amino acid matrix is 2Δ<em>ln</em> which is assessed on Chi<sup>2</sup> tables where Δ<em>ln</em> is difference in likelihood between model 1 and model 2. The degrees of freedom hmmmm... is the number of parameters that are constrained, this involves the number of taxa minus the number of parameter differences between the models. This makes performing 50 000 taxa ... tricky because the degrees of freedom will be outside the Chi<sup>2</sup> table.</p>
<p>iqtree will perform a tree and provide a likelihood for 50+K taxa. There are two limitations,</p>
<ul>
<li>No PHAT model</li>
<li>degrees of freedom (as mentioned)</li>
</ul>
<p>Overcoming df is simply to say whatever is the lowest likelihood is the best model and simply state Δ<em>ln</em> without bothering with significance testing. If the difference between models is large you'll get away with that (needs to be big) and for 50+K taxa it will be need to be huge.</p>
<p>Getting PHAT (the model below the ---- might be PHAT) into iqtree is certainly doable because of <a href="https://academic.oup.com/sysbio/article/70/5/1046/6146362" rel="nofollow noreferrer">QMaker</a>. However, how this is done I don't know, because I never done this calculation with iqtree, you'll need to read the manual. For traditional packages such as PAML it's easy ... you make a tree then specific the custom amino acid model and it spits out the likelihood. For the number of taxa involved, PAML is not possible, it's impossibly big for PAML, even 1000 is a lot of PAML.</p>
<p>iqtree does have LG, WAG and JTT amino acid models and will almost certainly have all the BLOSUM matrices and you just specify those at the command line. Again a comparison with PHAT needs to specify a custom matrix and I don't know how to do that in iqtree. The only thing I know is it must be doable.</p>
<p>iqtree is the most efficient of the maximum likelihood phylogenetic methods by a very long way. Hence for dealing with 50+K is why I'm recommending it here.</p>
<p><strong>Note</strong> The theory about whats going on and why the whole thing needs phylogenetics is very complicated. Basically its an entire lecture on phylogenetic theory.</p>
<hr />
<p>This I believe is PHAT (needs checking)</p>
<p><a href="https://github.com/biopython/biopython/pull/283" rel="nofollow noreferrer">https://github.com/biopython/biopython/pull/283</a></p>
<p>Note the code is copied and pasted from the above link, therefore it should be strictly regarded as a quote for copyright purposes.</p>
<pre><code>#!/usr/bin/env python
matrix = {
('W', 'F') :  0.0, ('L', 'R') : -6.0, ('I', 'I') :  5.0, ('Q', 'Q') :  9.0,
('W', 'N') : -5.0, ('V', 'I') :  3.0, ('H', 'T') : -4.0, ('H', 'P') : -6.0,
('W', 'V') : -4.0, ('Q', 'E') :  1.0, ('W', 'R') : -7.0, ('Q', 'A') : -3.0,
('H', 'H') : 11.0, ('H', 'D') : -1.0, ('L', 'N') : -3.0, ('Y', 'M') : -2.0,
('Y', 'I') : -3.0, ('Y', 'E') : -2.0, ('E', 'S') : -3.0, ('Y', 'A') : -3.0,
('Y', 'Y') : 11.0, ('T', 'C') : -1.0, ('E', 'C') : -7.0, ('Y', 'Q') :  0.0,
('E', 'G') : -3.0, ('V', 'A') :  1.0, ('C', 'C') :  7.0, ('M', 'R') : -6.0,
('P', 'T') : -4.0, ('V', 'E') : -5.0, ('P', 'P') : 13.0, ('I', 'T') : -1.0,
('K', 'S') : -5.0, ('R', 'G') : -5.0, ('I', 'P') : -4.0, ('R', 'C') : -8.0,
('A', 'T') :  0.0, ('K', 'K') :  5.0, ('A', 'P') :  3.0, ('V', 'M') :  1.0,
('I', 'D') : -5.0, ('K', 'C') :-10.0, ('K', 'G') : -5.0, ('R', 'S') : -6.0,
('F', 'Q') : -2.0, ('F', 'A') : -1.0, ('V', 'V') :  4.0, ('M', 'N') : -2.0,
('F', 'E') : -5.0, ('D', 'N') :  2.0, ('F', 'I') :  0.0, ('F', 'M') :  0.0,
('M', 'S') : -2.0, ('S', 'S') :  6.0, ('L', 'Q') : -3.0, ('W', 'E') : -7.0,
('W', 'A') : -4.0, ('W', 'M') : -4.0, ('H', 'S') : -2.0, ('W', 'I') : -4.0,
('S', 'C') :  1.0, ('L', 'A') : -1.0, ('L', 'E') : -5.0, ('W', 'Q') :  1.0,
('H', 'G') : -4.0, ('Q', 'N') :  2.0, ('H', 'C') : -7.0, ('L', 'M') :  2.0,
('W', 'Y') :  1.0, ('Y', 'N') :  2.0, ('E', 'P') : -5.0, ('Y', 'F') :  4.0,
('E', 'T') : -5.0, ('A', 'A') :  5.0, ('I', 'N') : -3.0, ('G', 'A') :  1.0,
('Y', 'V') : -3.0, ('E', 'D') :  6.0, ('W', 'H') : -3.0, ('Y', 'R') : -6.0,
('M', 'Q') : -1.0, ('P', 'S') : -3.0, ('R', 'H') : -4.0, ('A', 'C') :  1.0,
('R', 'D') : -7.0, ('K', 'P') : -4.0, ('L', 'D') : -5.0, ('K', 'T') : -6.0,
('V', 'N') : -3.0, ('M', 'A') : -1.0, ('K', 'H') : -5.0, ('V', 'R') : -7.0,
('P', 'C') : -8.0, ('M', 'E') : -5.0, ('A', 'S') :  2.0, ('T', 'T') :  3.0,
('R', 'T') : -6.0, ('I', 'G') : -2.0, ('R', 'P') : -7.0, ('K', 'D') : -5.0,
('I', 'C') : -3.0, ('F', 'R') : -7.0, ('F', 'V') : -1.0, ('L', 'C') : -2.0,
('F', 'F') :  6.0, ('D', 'A') : -5.0, ('F', 'N') : -1.0, ('W', 'D') : -7.0,
('L', 'P') : -5.0, ('Q', 'S') : -1.0, ('N', 'C') : -2.0, ('N', 'G') : -1.0,
('H', 'N') :  4.0, ('W', 'T') : -7.0, ('Q', 'G') : -2.0, ('W', 'P') : -6.0,
('Q', 'C') : -5.0, ('N', 'S') :  1.0, ('L', 'H') : -4.0, ('L', 'L') :  4.0,
('G', 'T') : -1.0, ('M', 'M') :  6.0, ('G', 'P') : -3.0, ('Y', 'K') : -4.0,
('Y', 'G') : -3.0, ('Y', 'C') : -1.0, ('E', 'A') : -5.0, ('E', 'E') : 12.0,
('Y', 'S') : -2.0, ('M', 'P') : -5.0, ('V', 'C') : -2.0, ('M', 'T') :  0.0,
('V', 'G') : -2.0, ('R', 'E') : -6.0, ('V', 'K') : -8.0, ('K', 'Q') : -1.0,
('R', 'A') : -6.0, ('I', 'R') : -6.0, ('N', 'A') : -2.0, ('V', 'S') : -2.0,
('M', 'D') : -5.0, ('M', 'H') : -4.0, ('K', 'A') : -7.0, ('R', 'Q') : -2.0,
('K', 'E') : -4.0, ('F', 'S') : -2.0, ('I', 'K') : -7.0, ('D', 'P') : -5.0,
('D', 'T') : -5.0, ('I', 'M') :  3.0, ('F', 'C') :  0.0, ('W', 'L') : -3.0,
('F', 'G') : -2.0, ('F', 'K') : -7.0, ('F', 'T') : -2.0, ('D', 'D') : 12.0,
('Q', 'T') : -3.0, ('W', 'G') : -5.0, ('Q', 'P') : -3.0, ('W', 'C') : -4.0,
('W', 'K') : -8.0, ('H', 'Q') :  2.0, ('Q', 'D') :  0.0, ('W', 'W') : 11.0,
('V', 'L') :  1.0, ('L', 'G') : -2.0, ('W', 'S') : -5.0, ('L', 'K') : -7.0,
('N', 'P') : -4.0, ('H', 'E') : -1.0, ('N', 'T') : -1.0, ('H', 'A') : -3.0,
('Y', 'L') : -2.0, ('Y', 'H') :  3.0, ('G', 'S') :  1.0, ('Y', 'D') : -4.0,
('V', 'Q') : -3.0, ('L', 'T') : -1.0, ('G', 'G') :  9.0, ('G', 'C') : -2.0,
('E', 'N') :  0.0, ('Y', 'T') : -3.0, ('Y', 'P') : -5.0, ('R', 'N') : -3.0,
('V', 'D') : -5.0, ('K', 'R') : -1.0, ('V', 'H') : -5.0, ('I', 'Q') : -3.0,
('V', 'P') : -4.0, ('M', 'C') : -2.0, ('K', 'N') : -2.0, ('V', 'T') :  0.0,
('M', 'G') : -1.0, ('T', 'S') :  1.0, ('I', 'E') : -5.0, ('M', 'K') : -6.0,
('I', 'A') :  0.0, ('N', 'N') : 11.0, ('R', 'R') :  9.0, ('F', 'P') : -5.0,
('L', 'I') :  2.0, ('I', 'S') : -2.0, ('D', 'S') : -4.0, ('L', 'S') : -2.0,
('I', 'H') : -5.0, ('F', 'D') : -5.0, ('D', 'C') : -7.0, ('F', 'H') : -2.0,
('D', 'G') : -2.0, ('F', 'L') :  1.0
}

chars = [aa[0] for aa in matrix] + [aa[1] for aa in matrix]
chars = sorted(list(set(chars)))

output = &quot; ,%s\n&quot; % &quot;,&quot;.join(chars)

for i, aa1 in enumerate(chars):
    output += aa1
    for j, aa2 in enumerate(chars):
        if (aa1, aa2) in matrix:
            output += &quot;,&quot; + str(matrix[(aa1, aa2)])
        else:
            output += &quot;,&quot; + str(matrix[(aa2, aa1)])
    output += &quot;\n&quot;
print(output)
</code></pre>

CC BY-SA 4.0
2023-10-09T14:20:10.240Z
2023-10-09T14:26:00.857Z
2023-10-09T14:26:00.857Z
M__
<p>I am looking for methods to compare substitution matrices to determine which one is better for a given set of data / type of protein.<br />
Do you know such methods / tools?</p>
<p>e.g.  BLOSUM62 vs PHAT for 50k - 100k+ sequences</p>

<phylogenetics><phylogeny><matrix><substitution><amino-acids>
Compare substitution matrices
# Compare substitution matrices

(Asked by: username_0 on 2023-10-06)

I am looking for methods to compare substitution matrices to determine which one is better for a given set of data / type of protein.  
Do you know such methods / tools?

e.g. BLOSUM62 vs PHAT for 50k - 100k+ sequences

**Tags:** phylogenetics, phylogeny, matrix, substitution, amino-acids

---

## 1 answers

### Answer by username_1 on 2023-10-09. Score: 2

The technique is phylogenetics is called likelihood ratio and for a regular amino acid matrix is 2Δln which is assessed on Chi2 tables where Δln is difference in likelihood between model 1 and model 2. The degrees of freedom hmmmm... is the number of parameters that are constrained, this involves the number of taxa minus the number of parameter differences between the models. This makes performing 50 000 taxa ... tricky because the degrees of freedom will be outside the Chi2 table.
iqtree will perform a tree and provide a likelihood for 50+K taxa. There are two limitations,

No PHAT model
degrees of freedom (as mentioned)

Overcoming df is simply to say whatever is the lowest likelihood is the best model and simply state Δln without bothering with significance testing. If the difference between models is large you'll get away with that (needs to be big) and for 50+K taxa it will be need to be huge.
Getting PHAT (the model below the ---- might be PHAT) into iqtree is certainly doable because of QMaker. However, how this is done I don't know, because I never done this calculation with iqtree, you'll need to read the manual. For traditional packages such as PAML it's easy ... you make a tree then specific the custom amino acid model and it spits out the likelihood. For the number of taxa involved, PAML is not possible, it's impossibly big for PAML, even 1000 is a lot of PAML.
iqtree does have LG, WAG and JTT amino acid models and will almost certainly have all the BLOSUM matrices and you just specify those at the command line. Again a comparison with PHAT needs to specify a custom matrix and I don't know how to do that in iqtree. The only thing I know is it must be doable.
iqtree is the most efficient of the maximum likelihood phylogenetic methods by a very long way. Hence for dealing with 50+K is why I'm recommending it here.
Note The theory about whats going on and why the whole thing needs phylogenetics is very complicated. Basically its an entire lecture on phylogenetic theory.

This I believe is PHAT (needs checking)
https://github.com/biopython/biopython/pull/283
Note the code is copied and pasted from the above link, therefore it should be strictly regarded as a quote for copyright purposes.
#!/usr/bin/env python
matrix = {
('W', 'F') :  0.0, ('L', 'R') : -6.0, ('I', 'I') :  5.0, ('Q', 'Q') :  9.0,
('W', 'N') : -5.0, ('V', 'I') :  3.0, ('H', 'T') : -4.0, ('H', 'P') : -6.0,
('W', 'V') : -4.0, ('Q', 'E') :  1.0, ('W', 'R') : -7.0, ('Q', 'A') : -3.0,
('H', 'H') : 11.0, ('H', 'D') : -1.0, ('L', 'N') : -3.0, ('Y', 'M') : -2.0,
('Y', 'I') : -3.0, ('Y', 'E') : -2.0, ('E', 'S') : -3.0, ('Y', 'A') : -3.0,
('Y', 'Y') : 11.0, ('T', 'C') : -1.0, ('E', 'C') : -7.0, ('Y', 'Q') :  0.0,
('E', 'G') : -3.0, ('V', 'A') :  1.0, ('C', 'C') :  7.0, ('M', 'R') : -6.0,
('P', 'T') : -4.0, ('V', 'E') : -5.0, ('P', 'P') : 13.0, ('I', 'T') : -1.0,
('K', 'S') : -5.0, ('R', 'G') : -5.0, ('I', 'P') : -4.0, ('R', 'C') : -8.0,
('A', 'T') :  0.0, ('K', 'K') :  5.0, ('A', 'P') :  3.0, ('V', 'M') :  1.0,
('I', 'D') : -5.0, ('K', 'C') :-10.0, ('K', 'G') : -5.0, ('R', 'S') : -6.0,
('F', 'Q') : -2.0, ('F', 'A') : -1.0, ('V', 'V') :  4.0, ('M', 'N') : -2.0,
('F', 'E') : -5.0, ('D', 'N') :  2.0, ('F', 'I') :  0.0, ('F', 'M') :  0.0,
('M', 'S') : -2.0, ('S', 'S') :  6.0, ('L', 'Q') : -3.0, ('W', 'E') : -7.0,
('W', 'A') : -4.0, ('W', 'M') : -4.0, ('H', 'S') : -2.0, ('W', 'I') : -4.0,
('S', 'C') :  1.0, ('L', 'A') : -1.0, ('L', 'E') : -5.0, ('W', 'Q') :  1.0,
('H', 'G') : -4.0, ('Q', 'N') :  2.0, ('H', 'C') : -7.0, ('L', 'M') :  2.0,
('W', 'Y') :  1.0, ('Y', 'N') :  2.0, ('E', 'P') : -5.0, ('Y', 'F') :  4.0,
('E', 'T') : -5.0, ('A', 'A') :  5.0, ('I', 'N') : -3.0, ('G', 'A') :  1.0,
('Y', 'V') : -3.0, ('E', 'D') :  6.0, ('W', 'H') : -3.0, ('Y', 'R') : -6.0,
('M', 'Q') : -1.0, ('P', 'S') : -3.0, ('R', 'H') : -4.0, ('A', 'C') :  1.0,
('R', 'D') : -7.0, ('K', 'P') : -4.0, ('L', 'D') : -5.0, ('K', 'T') : -6.0,
('V', 'N') : -3.0, ('M', 'A') : -1.0, ('K', 'H') : -5.0, ('V', 'R') : -7.0,
('P', 'C') : -8.0, ('M', 'E') : -5.0, ('A', 'S') :  2.0, ('T', 'T') :  3.0,
('R', 'T') : -6.0, ('I', 'G') : -2.0, ('R', 'P') : -7.0, ('K', 'D') : -5.0,
('I', 'C') : -3.0, ('F', 'R') : -7.0, ('F', 'V') : -1.0, ('L', 'C') : -2.0,
('F', 'F') :  6.0, ('D', 'A') : -5.0, ('F', 'N') : -1.0, ('W', 'D') : -7.0,
('L', 'P') : -5.0, ('Q', 'S') : -1.0, ('N', 'C') : -2.0, ('N', 'G') : -1.0,
('H', 'N') :  4.0, ('W', 'T') : -7.0, ('Q', 'G') : -2.0, ('W', 'P') : -6.0,
('Q', 'C') : -5.0, ('N', 'S') :  1.0, ('L', 'H') : -4.0, ('L', 'L') :  4.0,
('G', 'T') : -1.0, ('M', 'M') :  6.0, ('G', 'P') : -3.0, ('Y', 'K') : -4.0,
('Y', 'G') : -3.0, ('Y', 'C') : -1.0, ('E', 'A') : -5.0, ('E', 'E') : 12.0,
('Y', 'S') : -2.0, ('M', 'P') : -5.0, ('V', 'C') : -2.0, ('M', 'T') :  0.0,
('V', 'G') : -2.0, ('R', 'E') : -6.0, ('V', 'K') : -8.0, ('K', 'Q') : -1.0,
('R', 'A') : -6.0, ('I', 'R') : -6.0, ('N', 'A') : -2.0, ('V', 'S') : -2.0,
('M', 'D') : -5.0, ('M', 'H') : -4.0, ('K', 'A') : -7.0, ('R', 'Q') : -2.0,
('K', 'E') : -4.0, ('F', 'S') : -2.0, ('I', 'K') : -7.0, ('D', 'P') : -5.0,
('D', 'T') : -5.0, ('I', 'M') :  3.0, ('F', 'C') :  0.0, ('W', 'L') : -3.0,
('F', 'G') : -2.0, ('F', 'K') : -7.0, ('F', 'T') : -2.0, ('D', 'D') : 12.0,
('Q', 'T') : -3.0, ('W', 'G') : -5.0, ('Q', 'P') : -3.0, ('W', 'C') : -4.0,
('W', 'K') : -8.0, ('H', 'Q') :  2.0, ('Q', 'D') :  0.0, ('W', 'W') : 11.0,
('V', 'L') :  1.0, ('L', 'G') : -2.0, ('W', 'S') : -5.0, ('L', 'K') : -7.0,
('N', 'P') : -4.0, ('H', 'E') : -1.0, ('N', 'T') : -1.0, ('H', 'A') : -3.0,
('Y', 'L') : -2.0, ('Y', 'H') :  3.0, ('G', 'S') :  1.0, ('Y', 'D') : -4.0,
('V', 'Q') : -3.0, ('L', 'T') : -1.0, ('G', 'G') :  9.0, ('G', 'C') : -2.0,
('E', 'N') :  0.0, ('Y', 'T') : -3.0, ('Y', 'P') : -5.0, ('R', 'N') : -3.0,
('V', 'D') : -5.0, ('K', 'R') : -1.0, ('V', 'H') : -5.0, ('I', 'Q') : -3.0,
('V', 'P') : -4.0, ('M', 'C') : -2.0, ('K', 'N') : -2.0, ('V', 'T') :  0.0,
('M', 'G') : -1.0, ('T', 'S') :  1.0, ('I', 'E') : -5.0, ('M', 'K') : -6.0,
('I', 'A') :  0.0, ('N', 'N') : 11.0, ('R', 'R') :  9.0, ('F', 'P') : -5.0,
('L', 'I') :  2.0, ('I', 'S') : -2.0, ('D', 'S') : -4.0, ('L', 'S') : -2.0,
('I', 'H') : -5.0, ('F', 'D') : -5.0, ('D', 'C') : -7.0, ('F', 'H') : -2.0,
('D', 'G') : -2.0, ('F', 'L') :  1.0
}

chars = [aa[0] for aa in matrix] + [aa[1] for aa in matrix]
chars = sorted(list(set(chars)))

output = " ,%s\n" % ",".join(chars)

for i, aa1 in enumerate(chars):
    output += aa1
    for j, aa2 in enumerate(chars):
        if (aa1, aa2) in matrix:
            output += "," + str(matrix[(aa1, aa2)])
        else:
            output += "," + str(matrix[(aa2, aa1)])
    output += "\n"
print(output)


---

