Repository: cogent3/cogent3
Pull request #1447 (2023-05-19 07:18:44)

ENH: more alignment quality measures


Patch:
diff --git a/changelog.d/20230519_171446_Gavin.Huttley.md b/changelog.d/20230519_171446_Gavin.Huttley.md
new file mode 100644
index 000000000..5cb288a15
--- /dev/null
+++ b/changelog.d/20230519_171446_Gavin.Huttley.md
@@ -0,0 +1,45 @@
+<!--
+A new scriv changelog fragment.
+
+Uncomment the section that is right (remove the HTML comment wrapper).
+-->
+
+<!--
+### Contributers
+
+- A bullet item for the Contributers category.
+
+-->
+### ENH
+
+- Added new alignment quality measures as apps, and the ability to invoke them
+  from the Alignment.alignment_quality() method. The new apps are the
+  Information content measure of Hertz and Stormo (denoted `ic_score`), a
+  variant on the the sum of pairs measure of Carillo and Lipman
+  (denoted `sp_score`), and the log-liklelihood produced by the cogent3
+  progressive-HMM aligner (denoted `cogent3_score`)
+
+<!--
+### BUG
+
+- A bullet item for the BUG category.
+
+-->
+<!--
+### DOC
+
+- A bullet item for the DOC category.
+
+-->
+<!--
+### Deprecations
+
+- A bullet item for the Deprecations category.
+
+-->
+<!--
+### Discontinued
+
+- A bullet item for the Discontinued category.
+
+-->
diff --git a/src/cogent3/app/align.py b/src/cogent3/app/align.py
index 3b4783cda..309992a20 100644
--- a/src/cogent3/app/align.py
+++ b/src/cogent3/app/align.py
@@ -1,8 +1,11 @@
 import warnings
 
 from bisect import bisect_left
+from itertools import combinations
 from typing import Union
 
+from numpy import array
+
 from cogent3 import make_tree
 from cogent3.align import (
     global_pairwise,
@@ -14,7 +17,9 @@ from cogent3.app import dist
 from cogent3.core.alignment import Aligned, Alignment
 from cogent3.core.location import gap_coords_to_map
 from cogent3.core.moltype import get_moltype
+from cogent3.evolve.fast_distance import get_distance_calculator
 from cogent3.evolve.models import get_model
+from cogent3.maths.util import safe_log
 
 from .composable import NotCompleted, define_app
 from .tree import quick_tree, scale_branches
@@ -528,3 +533,157 @@ class progressive_align:
                 result = NotCompleted("ERROR", self, err.args[0], source=seqs)
                 return result
         return result
+
+
+@define_app
+class ic_score:
+    """compute the Information Content alignment quality score
+
+    Returns
+    -------
+    The score or 0.0 if it cannot be computed.
+
+    Notes
+    -----
+    Based on eq. (2) in noted reference.
+
+    Hertz, G. Z. & Stormo, G. D. Identifying DNA and protein patterns with
+    statistically significant alignments of multiple sequences.
+    Bioinformatics vol. 15 563–577 (Oxford University Press, 1999)
+    """
+
+    def __init__(self, equifreq_mprobs=True):
+        """
+        Parameters
+        ----------
+        equifreq_mprobs : bool
+            If true, specifies equally frequent motif probabilities.
+        """
+        self._equi_frequent = equifreq_mprobs
+
+    def main(self, aln: AlignedSeqsType) -> float:
+        counts = aln.counts_per_pos(include_ambiguity=False, allow_gap=False)
+        if counts.array.max() == 0 or aln.num_seqs == 1:
+            return 0.0
+
+        motif_probs = aln.get_motif_probs(include_ambiguity=False, allow_gap=False)
+
+        if self._equi_frequent:
+            # we reduce motif_probs to observed states
+            motif_probs = {m: v for m, v in motif_probs.items() if v > 0}
+            num_motifs = len(motif_probs)
+            motif_probs = {m: 1 / num_motifs for m in motif_probs}
+
+        p = array([motif_probs.get(b, 0.0) for b in counts.motifs])
+
+        cols = p != 0
+        p = p[cols]
+        counts = counts.array[:, cols]
+        frequency = counts / aln.num_seqs
+        log_f = safe_log(frequency / p)
+        I_seq = log_f * frequency
+        return I_seq.sum()
+
+
+@define_app
+def cogent3_score(aln: AlignedSeqsType) -> float:
+    """returns the cogent3 log-likelihood as an alignment quality score
+
+    Parameters
+    ----------
+    aln
+        An alignment instance.
+
+    Returns
+    -------
+    returns the log-likelihood, or 0.0 if the alignment does not have the
+    score
+
+    Notes
+    -----
+    The cogent3 tree_align algorithm is a progressive-HMM. It records
+    the log-likelihood of the alignment (in addition to all other
+    the parameter values, including the guide tree) in
+    ``alignment.info['align_params']``. This can be used as a measure of
+    alignment quality.
+
+    The instance must have been aligned using cogent3 tree_align. In addition,
+    if the alignment has been saved, it has must have been serialised
+    using a format that preserves the score.
+    """
+    if aln.num_seqs == 1 or len(aln) == 0:
+        return 0.0
+
+    align_params = aln.info.get("align_params", {})
+    return align_params.get("lnL", 0.0)
+
+
+@define_app
+class sp_score:
+    """Compute a variant of Sum of Pairs alignment quality score
+
+    Returns
+    -------
+    The score or 0.0 if it cannot be computed.
+
+    Notes
+    -----
+    We first compute pairwise genetic distances by the user specified
+    method. For each pair, this genetic distance is converted into an
+    estimate of the number of unchanged positions (with the minimum being
+    zero). The result is a sequence similarity matrix.
+
+    We compute pairwise affine gap distance using the user specified
+    values. The first matrix minus the second gives the statistic for
+    individual pairs and their sum is the Sum of Pairs score.
+
+    Carrillo, H. & Lipman, D. The Multiple Sequence Alignment Problem in
+    Biology. Siam J Appl Math 48, 1073–1082 (1988).
+    """
+
+    def __init__(
+        self, calc: str = "JC69", gap_insert: float = 12.0, gap_extend: float = 1.0
+    ):
+        """
+        Parameters
+        ----------
+        calc
+            name of the pairwise distance calculator
+        gap_insert
+            gap insertion penalty
+        gap_extend
+            gap extension penalty
+
+        Notes
+        -----
+        see available_distances() for the list of available methods.
+        """
+        self._insert = gap_insert
+        self._extend = gap_extend
+        self._calc = get_distance_calculator(calc)
+        self._gap_calc = dist.gap_dist(gap_insert=gap_insert, gap_extend=gap_extend)
+
+    def main(self, aln: AlignedSeqsType) -> float:
+        # this is a distance matrix, we want a similarity matrix
+        # because genetic distances can exceed 1, we need to adjust by
+        # multiplying by the length of the alignment to get the estimated
+        # number of changes and subtract that from the alignment length
+        if aln.num_seqs == 1 or len(aln) == 0:
+            return 0.0
+        self._calc(aln, show_progress=False)
+        dmat = self._calc.dists
+        lengths = self._calc.lengths
+        for i, j in combinations(range(aln.num_seqs), 2):
+            n1, n2 = dmat.names[i], dmat.names[j]
+            length = lengths[n1, n2]
+            unchanged = length - dmat[i, j] * length
+            # for very divergent cases, unchanged could be < 0,
+            # limit this to 0
+            dmat[i, j] = dmat[j, i] = max(unchanged, 0.0)
+        # compute the pairwise gap as number of gaps * gap_insert + length
+        # of gaps * gap_extend
+        gap_scores = self._gap_calc(aln)
+        # make sure the two dmats in same name order
+        gap_scores = gap_scores.take_dists(dmat.names)
+        dmat.array -= gap_scores
+        return dmat.array.sum() / 2
diff --git a/src/cogent3/app/dist.py b/src/cogent3/app/dist.py
index 7e88840d9..d91fc1f4e 100644
--- a/src/cogent3/app/dist.py
+++ b/src/cogent3/app/dist.py
@@ -1,11 +1,12 @@
 import itertools
 
 from copy import deepcopy
+from itertools import combinations
 from typing import Union
 
 import numpy
 
-from numpy import log, polyval, triu_indices, zeros_like
+from numpy import polyval, triu_indices
 
 from cogent3 import get_moltype, make_tree
 from cogent3.evolve.fast_distance import (
@@ -14,6 +15,7 @@ from cogent3.evolve.fast_distance import (
 )
 from cogent3.evolve.models import get_model
 from cogent3.maths.distance_transform import jaccard
+from cogent3.util.misc import get_true_spans
 
 from .composable import define_app
 from .typing import (
@@ -239,3 +241,42 @@ def approx_jc69(
     arr.T[upper_indices] = arr[upper_indices]
     result.array = arr
     return result
+
+
+@define_app
+class gap_dist:
+    """compute the pairwise difference in gaps using affine gap score"""
+
+    def __init__(self, gap_insert: float = 12.0, gap_extend: float = 1.0):
+        """
+        Parameters
+        ----------
+        gap_insert
+            gap insertion penalty
+        gap_extend
+            gap extension penalty
+        """
+        self._insert = gap_insert
+        self._extend = gap_extend
+
+    def main(self, aln: AlignedSeqsType) -> PairwiseDistanceType:
+        gap_diffs = {}
+        gap_data = dict(zip(aln.names, aln.get_gap_array()))
+        # convert each gap vector to gap spans
+        for k, gap_vector in gap_data.items():
+            gap_data[k] = set(tuple(pr) for pr in get_true_spans(gap_vector).tolist())
+
+        for i, j in combinations(range(aln.num_seqs), 2):
+            n1, n2 = aln.names[i], aln.names[j]
+
+            # unique gaps
+            unique = gap_data[n1] ^ gap_data[n2]
+
+            # compute the score as
+            # number of gaps * gap_insert + total length of gaps * gap_extend
+            score = self._insert * len(unique) + self._extend * sum(
+                l for p, l in unique
+            )
+            gap_diffs[(n1, n2)] = gap_diffs[(n2, n1)] = score
+
+        return DistanceMatrix(gap_diffs)
diff --git a/src/cogent3/core/alignment.py b/src/cogent3/core/alignment.py
index 2af6ce3e7..43ab9d4cd 100644
--- a/src/cogent3/core/alignment.py
+++ b/src/cogent3/core/alignment.py
@@ -2560,43 +2560,22 @@ class AlignmentI(object):
     default_gap = "-"  # default gap character for padding
     gap_chars = dict.fromkeys("-?")  # default gap chars for comparisons
 
-    def alignment_quality(self, equifreq_mprobs=True):
+    def alignment_quality(self, app_name: str = "ic_score", **kwargs):
         """
-        Computes the alignment quality for an alignment based on eq. (2) in noted reference.
+        Computes the alignment quality using the indicated app
 
         Parameters
         ----------
-        equifreq_mprobs : bool
-            If true, specifies equally frequent motif probabilities.
+        app_name
+            name of an alignment score calculating app, e.g. 'ic_score',
+            'cogent3_score', 'sp_score'
 
-        Notes
-        -----
-        G. Z. Hertz, G. D. Stormo - Published 1999, Bioinformatics, vol. 15 pg. 563-577.
-
-        The alignment quality statistic is a log-likelihood ratio (computed using log2)
-        of the observed alignment column freqs versus the expected.
+        kwargs
+            keyword arguments to be passed to the app. Use
+            ``cogent3.app_help(app_name)`` to see the available options.
         """
-        counts = self.counts_per_pos()
-        if counts.array.max() == 0 or len(self.seqs) == 1:
-            return None
-
-        motif_probs = self.get_motif_probs()
-
-        if equifreq_mprobs:
-            # we reduce motif_probs to observed states
-            motif_probs = {m: v for m, v in motif_probs.items() if v > 0}
-            num_motifs = len(motif_probs)
-            motif_probs = {m: 1 / num_motifs for m in motif_probs}
-
-        p = array([motif_probs.get(b, 0.0) for b in counts.motifs])
-
-        cols = p != 0
-        p = p[cols]
-        counts = counts.array[:, cols]
-        frequency = counts / self.num_seqs
-        log_f = safe_log(frequency / p)
-        I_seq = log_f * frequency
-        return I_seq.sum()
+        app = cogent3.get_app(app_name, **kwargs)
+        return app(self)
 
     def iter_positions(self, pos_order=None):
         """Iterates over positions in the alignment, in order.
diff --git a/src/cogent3/util/misc.py b/src/cogent3/util/misc.py
index af3a72e5d..14c7f6376 100644
--- a/src/cogent3/util/misc.py
+++ b/src/cogent3/util/misc.py
@@ -1,5 +1,7 @@
 """Generally useful utility classes and methods.
 """
+from __future__ import annotations
+
 import inspect
 import os
 import re
@@ -11,7 +13,7 @@ from warnings import warn
 
 import numpy
 
-from numpy import array, finfo, float64
+from numpy import array, finfo, float64, ndarray, zeros
 
 
 def _adjusted_gt_minprob_vector(probs, minprob):
@@ -1034,3 +1036,45 @@ def in_jupyter() -> bool:
         val = False
 
     return val
+
+
+def get_true_spans(arr: ndarray, absolute_pos: bool = True) -> ndarray:
+    """
+    returns array of [[start position, length],...] of True runs
+
+    Parameters
+    ----------
+    arr
+        bool array
+    absolute_pos
+        the absolute position in self, otherwise the position
+        represents the index in the underlying series of False values
+        between which the True run occurs.
+
+    Notes
+    -----
+    Designed for use with individual records from Alignment.get_gap_array(). This
+    method returns a bool array with the gap state as True.
+    """
+    result = zeros((arr.shape[0], 2), dtype=int)
+    true_run = False
+    num_runs = 0
+    cum_sum = 0
+    for i, v in enumerate(arr):
+        if v and true_run:
+            result[num_runs][1] += 1
+            cum_sum += 1
+        elif not v and true_run:
+            # we just ended a run of False values
+            num_runs += 1
+            true_run = False
+        elif v:
+            result[num_runs][:] = [i if absolute_pos else i - cum_sum, 1]
+            true_run = True
+            cum_sum += 1
+
+    if true_run:
+        # finished in a run
+        num_runs += 1
+
+    return result[:num_runs]


Test patch:
diff --git a/tests/test_app/test_align.py b/tests/test_app/test_align.py
index c47aead52..50c9c064a 100644
--- a/tests/test_app/test_align.py
+++ b/tests/test_app/test_align.py
@@ -1,7 +1,14 @@
-from unittest import TestCase, main
+from unittest import TestCase
+
+import numpy
+import pytest
+
+from numpy import log2
+from numpy.testing import assert_allclose
 
 from cogent3 import (
     DNA,
+    get_app,
     get_moltype,
     make_aligned_seqs,
     make_tree,
@@ -19,7 +26,7 @@ from cogent3.app.align import (
     pairwise_to_multiple,
 )
 from cogent3.app.composable import NotCompleted
-from cogent3.core.alignment import Aligned, ArrayAlignment
+from cogent3.core.alignment import Aligned, Alignment, ArrayAlignment
 from cogent3.core.location import gap_coords_to_map
 
 
@@ -466,5 +473,110 @@ class GapOffsetTests(TestCase):
             self.assertEqual(aln_pos - aln2seq[aln_pos], seq_pos)
 
 
-if __name__ == "__main__":
-    main()
+@pytest.mark.parametrize("cls", (Alignment, ArrayAlignment))
+def test_information_content_score(cls):
+    """Tests that the alignment_quality generates the right alignment quality
+    value based on the Hertz-Stormo metric. expected values are hand calculated
+    using the formula in the paper."""
+    app_equifreq = get_app("ic_score", equifreq_mprobs=True)
+    app_not_equifreq = get_app("ic_score", equifreq_mprobs=False)
+
+    aln = cls(["AATTGA", "AGGTCC", "AGGATG", "AGGCGT"], moltype="dna")
+    got = app_equifreq(aln)
+    expect = log2(4) + (3 / 2) * log2(3) + (1 / 2) * log2(2) + (1 / 2) * log2(2)
+    assert_allclose(got, expect)
+    # should be the same with the default moltype too
+    aln = cls(["AATTGA", "AGGTCC", "AGGATG", "AGGCGT"])
+    got = app_equifreq(aln)
+    assert_allclose(got, expect)
+
+    aln = cls(["AAAC", "ACGC", "AGCC", "A-TC"], moltype="dna")
+    got = app_not_equifreq(aln)
+    expect = (
+        2 * log2(1 / 0.4)
+        + log2(1 / (4 * 0.4))
+        + (1 / 2) * log2(1 / (8 / 15))
+        + (1 / 4) * log2(1 / (4 / 15))
+    )
+    assert_allclose(got, expect)
+
+    # 1. Alignment just gaps - alignment_quality returns 0.0
+    aln = cls(["----", "----"])
+    got = app_equifreq(aln)
+    assert_allclose(got, 0.0)
+
+    # 2 Just one sequence - alignment_quality returns 0.0
+    aln = cls(["AAAC"])
+    got = app_equifreq(aln)
+    assert_allclose(got, 0.0)
+
+    # 3.1 Two seqs, one all gaps. (equifreq_mprobs=True)
+    aln = cls(["----", "ACAT"])
+    got = app_equifreq(aln)
+    assert_allclose(got, 1.1699250014423124)
+
+    # 3.2 Two seqs, one all gaps. (equifreq_mprobs=False)
+    aln = cls(["----", "AAAA"])
+    got = app_not_equifreq(aln)
+    assert_allclose(got, -2)
+
+
+@pytest.fixture(scope="function")
+def aln():
+    aligner = align_app.progressive_align(model="TN93", distance="TN93")
+    seqs = make_unaligned_seqs(_seqs, moltype=DNA)
+    return aligner(seqs)
+
+
+def test_cogent3_score(aln):
+    get_score = get_app("cogent3_score")
+    score = get_score(aln)
+    assert score < -100
+
+
+@pytest.mark.parametrize("del_all_params", (True, False))
+def test_cogent3_score_missing(aln, del_all_params):
+    get_score = get_app("cogent3_score")
+    if del_all_params:
+        aln.info.pop("align_params")
+    else:
+        aln.info["align_params"].pop("lnL")
+    score = get_score(aln)
+    assert score == 0.0
+
+
+def test_sp_score_exclude_gap():
+    # no gap penalty
+    app = get_app("sp_score", calc="pdist", gap_extend=0, gap_insert=0)
+    data = {"s1": "AAGAA-A", "s2": "-ATAATG", "s3": "C-TGG-G"}
+    # prop unchanged s1-s2, s1-s3
+    expect = sum([6 * 3 / 6, 0, 5 * 2 / 5])
+    aln = make_aligned_seqs(data, moltype="dna")
+    got = app.main(aln)
+    assert_allclose(got, expect)
+
+
+def test_sp_score_additive_gap():
+    # additive gap score
+    app = get_app("sp_score", calc="pdist", gap_extend=1, gap_insert=0)
+    data = {"s1": "AAGAA-A", "s2": "-ATAATG", "s3": "C-TGG-G"}
+    # match score
+    mscore = numpy.array([6 * 3 / 6, 0, 5 * 2 / 5])
+    # gap score
+    gscore = numpy.array([2, 1, 3])
+    aln = make_aligned_seqs(data, moltype="dna")
+    got = app.main(aln)
+    assert_allclose(got, (mscore - gscore).sum())
+
+
+def test_sp_score_affine_gap():
+    # affine gap score
+    app = get_app("sp_score", calc="pdist", gap_extend=1, gap_insert=2)
+    data = {"a": "AAGAA-A", "b": "-ATAATG", "c": "C-TGG-G"}
+    # match score
+    mscore = numpy.array([6 * 3 / 6, 0, 5 * 2 / 5])
+    # gap score
+    gscore = numpy.array([2 + 4, 2 + 1, 3 + 6])
+    aln = make_aligned_seqs(data, moltype="dna")
+    got = app.main(aln)
+    assert_allclose(got, (mscore - gscore).sum())
diff --git a/tests/test_app/test_dist.py b/tests/test_app/test_dist.py
index 1ef7e2866..ec5e90922 100644
--- a/tests/test_app/test_dist.py
+++ b/tests/test_app/test_dist.py
@@ -461,3 +461,64 @@ def test_symmetry_of_dists():
     app = approx_jc69()
     got = app(dists)
     assert_allclose(got.array, got.array.T)
+
+
+def test_gap_dist():
+    app = get_app("gap_dist", gap_insert=10, gap_extend=1)
+    # two sequences share a gap
+    data = {
+        "a": "TG----AATATGT------GAAAGAG",
+        "b": "TTGAAGAATATGT------GAAAGAG",
+        "c": "CTGAAGAACCTGTGAAAGTGAAAGAG",
+    }
+    aln = make_aligned_seqs(data, moltype="dna", array_align=True)
+    expect = {
+        ("a", "b"): 14.0,  # one gap diff of size 4
+        ("a", "c"): 30.0,
+        ("b", "c"): 16.0,
+    }
+    expect = DistanceMatrix(expect)
+    dmat = app.main(aln)
+    assert dmat.to_dict() == expect.to_dict()
+
+    # shared gap actually not shared, 3 events
+    data = {
+        "a": "TG----AATATGTA-----GAAAGAG",
+        "b": "TTGAAGAATATGTA------AAAGAG",
+        "c": "CTGAAGAACCTGTGAAAGTGAAAGAG",
+    }
+    aln = make_aligned_seqs(data, moltype="dna", array_align=True)
+    expect = {
+        ("a", "b"): 45.0,  # 3 gaps diff of size 15
+        ("a", "c"): 29.0,
+        ("b", "c"): 16.0,
+    }
+    expect = DistanceMatrix(expect)
+    dmat = app.main(aln)
+    assert dmat.to_dict() == expect.to_dict()
+
+    # additional gaps on either side of shared gap is two events
+    data = {
+        "a": "G--AG----A",
+        "b": "TGGAGT--GA",
+        "c": "TGGAGTGTGA",
+    }
+    aln = make_aligned_seqs(data, moltype="dna", array_align=True)
+    expect = {
+        ("a", "b"): 38,  # 3 gaps diff of size 8
+        ("a", "c"): 26.0,
+        ("b", "c"): 12.0,
+    }
+    expect = DistanceMatrix(expect)
+    dmat = app.main(aln)
+    assert dmat.to_dict() == expect.to_dict()
+    data = {"a": "AAGAA-A", "b": "-ATAATG", "c": "C-TGG-G"}
+    aln = make_aligned_seqs(data, moltype="dna", array_align=True)
+    expect = {
+        ("a", "b"): 22.0,  # 2 gaps diff of size 2
+        ("a", "c"): 11.0,
+        ("b", "c"): 33.0,  # 3 gaps diff of size 2
+    }
+    expect = DistanceMatrix(expect)
+    dmat = app.main(aln)
+    assert dmat.to_dict() == expect.to_dict()
diff --git a/tests/test_core/test_alignment.py b/tests/test_core/test_alignment.py
index 92b62601c..aa2ca4725 100644
--- a/tests/test_core/test_alignment.py
+++ b/tests/test_core/test_alignment.py
@@ -15,6 +15,7 @@ from numpy import array, log2, nan, transpose
 from numpy.testing import assert_allclose, assert_equal
 
 from cogent3 import (
+    get_app,
     get_code,
     load_aligned_seqs,
     load_seq,
@@ -1409,18 +1410,7 @@ class AlignmentBaseTests(SequenceCollectionBaseTests):
     """
 
     def test_alignment_quality(self):
-        """Tests that the alignment_quality generates the right alignment quality
-        value based on the Hertz-Stormo metric. expected values are hand calculated
-        using the formula in the paper."""
-        aln = self.Class(["AATTGA", "AGGTCC", "AGGATG", "AGGCGT"], moltype="dna")
-        got = aln.alignment_quality(equifreq_mprobs=True)
-        expect = log2(4) + (3 / 2) * log2(3) + (1 / 2) * log2(2) + (1 / 2) * log2(2)
-        assert_allclose(got, expect)
-        # should be the same with the default moltype too
-        aln = self.Class(["AATTGA", "AGGTCC", "AGGATG", "AGGCGT"])
-        got = aln.alignment_quality(equifreq_mprobs=True)
-        assert_allclose(got, expect)
-
+        """check alignment method correctly invokes underlying app"""
         aln = self.Class(["AAAC", "ACGC", "AGCC", "A-TC"], moltype="dna")
         got = aln.alignment_quality(equifreq_mprobs=False)
         expect = (
@@ -1431,26 +1421,6 @@ class AlignmentBaseTests(SequenceCollectionBaseTests):
         )
         assert_allclose(got, expect)
 
-        # 1. Alignment just gaps - alignment_quality returns None
-        aln = self.Class(["----", "----"])
-        got = aln.alignment_quality(equifreq_mprobs=True)
-        self.assertIsNone(got)
-
-        # 2 Just one sequence - alignment_quality returns None
-        aln = self.Class(["AAAC"])
-        got = aln.alignment_quality(equifreq_mprobs=True)
-        self.assertIsNone(got)
-
-        # 3.1 Two seqs, one all gaps. (equifreq_mprobs=True)
-        aln = self.Class(["----", "ACAT"])
-        got = aln.alignment_quality(equifreq_mprobs=True)
-        assert_allclose(got, 1.1699250014423124)
-
-        # 3.2 Two seqs, one all gaps. (equifreq_mprobs=False)
-        aln = self.Class(["----", "AAAA"])
-        got = aln.alignment_quality(equifreq_mprobs=False)
-        assert_allclose(got, -2)
-
     def make_and_filter(self, raw, expected, motif_length, drop_remainder):
         # a simple filter func
         aln = self.Class(raw, info={"key": "value"})
@@ -3471,3 +3441,69 @@ def test_get_translation_error(cls, seqs):
     alignment = cls(data=seqs)
     with pytest.raises(TypeError):
         alignment.get_translation()
+
+
+@pytest.mark.parametrize("cls", (Alignment, ArrayAlignment))
+@pytest.mark.parametrize("method", ("ic_score", "cogent3_score", "sp_score"))
+def test_alignment_quality_methods(cls, method):
+    data = {
+        "DogFaced": "TG----AATATGT------GAAAGAG",
+        "FreeTaile": "TTGAAGAATATGT------GAAAGAG",
+        "LittleBro": "CTGAAGAACCTGTGAAAGTGAAAGAG",
+    }
+    expected_score = dict(
+        cogent3_score=-123.0, ic_score=32.93032499, sp_score=-25.25687109
+    )[method]
+    aln = cls(
+        data,
+        moltype="dna",
+        info=dict(align_params=dict(lnL=-123.0)),
+    )
+    app = get_app(method)
+    score = app(aln)
+    assert_allclose(score, expected_score)
+
+
+@pytest.mark.parametrize("method", ("ic_score", "cogent3_score", "sp_score"))
+def test_alignment_quality_methods_oneseq(method):
+    data = {
+        "DogFaced": "TG----AATATGT------GAAAGAG",
+    }
+    aln = make_aligned_seqs(
+        data,
+        moltype="dna",
+        info=dict(align_params=dict(lnL=-123.0)),
+    )
+    app = get_app(method)
+    score = app(aln)
+    assert_allclose(score, 0.0)
+
+
+@pytest.mark.parametrize("method", ("ic_score", "cogent3_score", "sp_score"))
+def test_alignment_quality_methods_zero_length(method):
+    data = {
+        "a": "",
+        "b": "",
+        "c": "",
+    }
+    aln = make_aligned_seqs(
+        data,
+        moltype="dna",
+        info=dict(align_params=dict(lnL=-123.0)),
+    )
+    app = get_app(method)
+    score = app(aln)
+    assert_allclose(score, 0.0)
+
+
+def test_get_gap_array_equivalence():
+    # make sure produced gap arrays are identical between the
+    # two Alignment classes
+    data = {
+        "DogFaced": "TG----AATATGT------GAAAGAG",
+        "FreeTaile": "TTGAAGAATATGT------GAAAGAG",
+        "LittleBro": "CTGAAGAACCTGTGAAAGTGAAAGAG",
+    }
+    array_aln = make_aligned_seqs(data, moltype="dna", array_align=True)
+    aln = make_aligned_seqs(data, moltype="dna", array_align=False)
+    assert_allclose(array_aln.get_gap_array(), aln.get_gap_array())
diff --git a/tests/test_util/test_misc.py b/tests/test_util/test_misc.py
index 4c98fabed..d56bf6256 100644
--- a/tests/test_util/test_misc.py
+++ b/tests/test_util/test_misc.py
@@ -7,6 +7,7 @@ from unittest import TestCase
 
 import pytest
 
+from numpy import array
 from numpy.testing import assert_allclose
 
 from cogent3.util.misc import (
@@ -33,6 +34,7 @@ from cogent3.util.misc import (
     get_object_provenance,
     get_run_start_indices,
     get_setting_from_environ,
+    get_true_spans,
     identity,
     is_char,
     is_char_or_noniterable,
@@ -1250,3 +1252,35 @@ _body_expect = ["Notes", "-----", "body"]
 def test_docstring_to_summary_rest(foo, sum_exp, body_exp):
     summary, body = docstring_to_summary_rest(foo.__doc__)
     assert summary == sum_exp and body.split() == body_exp
+
+
+def test_get_true_spans_absolute():
+    got = get_true_spans(array([0, 1, 1, 0, 1]))
+    assert_allclose(got, array([[1, 2], [4, 1]]))
+    got = get_true_spans(array([0, 0]))
+    assert not len(got)
+
+    got = get_true_spans(array([0, 0, 1, 1]))
+    assert_allclose(got, array([[2, 2]]))
+    got = get_true_spans(array([1, 0, 0, 1, 1]))
+    assert_allclose(got, array([[0, 1], [3, 2]]))
+    got = get_true_spans(array([1, 0, 0, 1, 1, 1, 1, 0, 1, 1]))
+    assert_allclose(got, array([[0, 1], [3, 4], [8, 2]]))
+    got = get_true_spans(
+        # abs  0  1  2  3  4  5  6  7  8  9
+        array([1, 0, 0, 1, 1, 1, 1, 0, 1, 1]),
+        absolute_pos=False
+        # u       0  1              2
+    )
+    assert_allclose(got, array([[0, 1], [2, 4], [3, 2]]))
+
+    got = get_true_spans(array([0, 0, 0, 1, 1, 1, 0, 0]))
+    assert_allclose(got, array([(3, 3)]))
+
+
+def test_get_true_spans_not_absolute():
+    got = get_true_spans(array([0, 1, 1, 0, 1]), absolute_pos=False)
+    assert_allclose(got, array([[1, 2], [2, 1]]))
+
+    got = get_true_spans(array([1, 0, 0, 1, 1, 1, 1, 0, 1, 1]), absolute_pos=False)
+    assert_allclose(got, array([[0, 1], [2, 4], [3, 2]]))