Repository: acorg/dark-matter

Issue:
Add compareAAReads function (to match compareDNAReads)

Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 997bd69..dacc548 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,9 @@
+## 3.0.63 Jan 14, 2019
+
+Added `compareAaReads` and `matchToString` to `aa.py`. Wrote tests in 
+`test_aa.py` for both. Moved `countPrint` to utils, used by `matchToString`
+in `dna.py` and `aa.py`. Added `compare-aa-sequences.py` to the bin. 
+
 ## 3.0.62 Dec 30, 2018
 
 Added `matchToString` to `dna.py` to allow printing of a DNA match.
diff --git a/bin/compare-aa-sequences.py b/bin/compare-aa-sequences.py
new file mode 100755
index 0000000..d19dc25
--- /dev/null
+++ b/bin/compare-aa-sequences.py
@@ -0,0 +1,152 @@
+#!/usr/bin/env python
+
+from __future__ import print_function, division
+
+import sys
+import argparse
+from os.path import join
+from math import log10
+
+from dark.aa import compareAaReads, matchToString
+from dark.fasta import FastaReads
+from dark.reads import (Reads, addFASTACommandLineOptions,
+                        parseFASTACommandLineOptions)
+from dark.process import Executor
+from dark.utils import parseRangeString
+
+
+def needle(reads):
+    """
+    Run a Needleman-Wunsch alignment and return the two sequences.
+
+    @param reads: An iterable of two reads.
+    @return: A C{Reads} instance with the two aligned sequences.
+    """
+    from tempfile import mkdtemp
+    from shutil import rmtree
+
+    dir = mkdtemp()
+
+    file1 = join(dir, 'file1.fasta')
+    with open(file1, 'w') as fp:
+        print(reads[0].toString('fasta'), end='', file=fp)
+
+    file2 = join(dir, 'file2.fasta')
+    with open(file2, 'w') as fp:
+        print(reads[1].toString('fasta'), end='', file=fp)
+
+    out = join(dir, 'result.fasta')
+
+    Executor().execute("needle -asequence '%s' -bsequence '%s' -auto "
+                       "-outfile '%s' -aformat fasta" % (
+                           file1, file2, out))
+
+    # Use 'list' in the following to force reading the FASTA from disk.
+    result = Reads(list(FastaReads(out)))
+    rmtree(dir)
+
+    return result
+
+
+parser = argparse.ArgumentParser(
+    formatter_class=argparse.ArgumentDefaultsHelpFormatter,
+    description=('Compare two amino acid sequences.'))
+
+parser.add_argument(
+    '--index1', type=int, default=1,
+    help='The (1-based) index in the input of the first sequence.')
+
+parser.add_argument(
+    '--index2', type=int, default=2,
+    help='The (1-based) index in the input of the second sequence.')
+
+parser.add_argument(
+    '--align', default=False, action='store_true',
+    help='If given, use needle to do a global alignment of the two sequences.')
+
+parser.add_argument(
+    '--alignmentFile',
+    help='The file to save the alignment to (implies --align).')
+
+parser.add_argument(
+    '--sites',
+    help=('Specify (1-based) sequence sites to keep. All other sites will '
+          'be ignored. The sites must be given in the form e.g., '
+          '24,100-200,260.'))
+
+parser.add_argument(
+    '--showDiffs', default=False, action='store_true',
+    help='Print (1-based) sites where the sequence amino acids differ.')
+
+addFASTACommandLineOptions(parser)
+args = parser.parse_args()
+
+keepSequences = set([args.index1 - 1, args.index2 - 1])
+
+reads = list(parseFASTACommandLineOptions(args).filter(
+    keepSequences=keepSequences))
+
+if len(reads) == 1:
+    if len(keepSequences) == 1:
+        # This is ok, they want to compare a sequence with itself.
+        reads = Reads([reads[0], reads[0]])
+    else:
+        print('Could not find both requested sequence indices. Exiting.')
+        sys.exit(1)
+elif len(reads) != 2:
+    print('Could not find both requested sequence indices. Exiting.')
+    sys.exit(1)
+
+if args.alignmentFile:
+    args.align = True
+
+if args.align:
+    len1, len2 = map(len, reads)
+    if len1 == len2:
+        print('Pre-alignment, sequence lengths were identical: %s' % len1)
+    else:
+        print('Pre-alignment, sequence lengths: %d, %d (difference %d)' % (
+            len1, len2, abs(len1 - len2)))
+
+    # Align.
+    reads = needle(reads)
+
+    if args.alignmentFile:
+        assert reads.save(args.alignmentFile) == 2
+
+offsets = (parseRangeString(args.sites, convertToZeroBased=True)
+           if args.sites else None)
+
+read1, read2 = reads
+len1, len2 = map(len, reads)
+identicalLengths = len1 == len2
+
+# Sanity check.
+if args.align:
+    assert identicalLengths
+
+match = compareAaReads(read1, read2, offsets=offsets)
+
+x = 'Post-alignment, sequence' if args.align else 'Sequence'
+if identicalLengths:
+    print('%s lengths are identical: %s' % (x, len1))
+else:
+    print('%s lengths: %d, %d (difference %d)' % (x, len1, len2,
+                                                  abs(len1 - len2)))
+
+print(matchToString(match, read1, read2, offsets=offsets))
+
+if args.showDiffs:
+    # Print all sites where the sequences differ.
+    width = int(log10(max(len1, len2))) + 1
+    headerPrinted = False
+    for site, (a, b) in enumerate(zip(read1.sequence, read2.sequence),
+                                  start=1):
+        if a != b:
+            if not headerPrinted:
+                print('Differences (site, %s, %s):' % (read1.id, read2.id))
+                headerPrinted = True
+            print('  %*d %s %s' % (width, site, a, b))
+
+    if not headerPrinted:
+        print('No sequence differences found.')
diff --git a/dark/__init__.py b/dark/__init__.py
index d9b0c2f..754f0cd 100644
--- a/dark/__init__.py
+++ b/dark/__init__.py
@@ -8,4 +8,4 @@ if sys.version_info < (2, 7):
 #
 # Remember to update ../CHANGELOG.md describing what's new in each version.
 
-__version__ = '3.0.62'
+__version__ = '3.0.63'
diff --git a/dark/aa.py b/dark/aa.py
index 8a2a849..4bfd478 100644
--- a/dark/aa.py
+++ b/dark/aa.py
@@ -1,3 +1,13 @@
+from dark.utils import countPrint
+try:
+    from itertools import zip_longest
+except ImportError:
+    # zip_longest does not exist in Python 2.7 itertools. We should be able
+    # to get it via from six.moves import zip_longest according to
+    # https://pythonhosted.org/six/index.html?highlight=zip_longest but
+    # that doesn't work for me.
+    from itertools import izip_longest as zip_longest
+
 # From https://en.wikipedia.org/wiki/Amino_acid
 #
 # Alanine          Ala     A
@@ -1093,3 +1103,140 @@ def clustersForSequence(sequence, propertyNames, missingAAValue=0):
     """
     return _propertiesOrClustersForSequence(
         sequence, propertyNames, PROPERTY_CLUSTERS, missingAAValue)
+
+
+def matchToString(aaMatch, read1, read2, indent='',
+                  offsets=None):
+    """
+    Format amino acid sequence match as a string.
+
+    @param aaMatch: A C{dict} returned by C{compareAaReads}.
+    @param read1: A C{Read} instance or an instance of one of its subclasses.
+    @param read2: A C{Read} instance or an instance of one of its subclasses.
+    @param indent: A C{str} to indent all returned lines with.
+    @param offsets: If not C{None}, a C{set} of offsets of interest that were
+        only considered when making C{match}.
+    @return: A C{str} describing the match.
+    """
+    match = aaMatch['match']
+    matchCount = match['matchCount']
+    gapMismatchCount = match['gapMismatchCount']
+    gapGapMismatchCount = match['gapGapMismatchCount']
+    nonGapMismatchCount = match['nonGapMismatchCount']
+
+    if offsets:
+        len1 = len2 = len(offsets)
+    else:
+        len1, len2 = map(len, (read1, read2))
+
+    result = []
+    append = result.append
+
+    append(countPrint('%sMatches' % indent, matchCount, len1, len2))
+    mismatchCount = (gapMismatchCount + gapGapMismatchCount +
+                     nonGapMismatchCount)
+    append(countPrint('%sMismatches' % indent, mismatchCount, len1, len2))
+    append(countPrint('%s  Not involving gaps (i.e., conflicts)' % (indent),
+                      nonGapMismatchCount, len1, len2))
+    append(countPrint('%s  Involving a gap in one sequence' % indent,
+                      gapMismatchCount, len1, len2))
+    append(countPrint('%s  Involving a gap in both sequences' % indent,
+                      gapGapMismatchCount, len1, len2))
+
+    for read, key in zip((read1, read2), ('read1', 'read2')):
+        append('%s  Id: %s' % (indent, read.id))
+        length = len(read)
+        append('%s    Length: %d' % (indent, length))
+        gapCount = len(aaMatch[key]['gapOffsets'])
+        append(countPrint('%s    Gaps' % indent, gapCount, length))
+        if gapCount:
+            append(
+                '%s    Gap locations (1-based): %s' %
+                (indent,
+                 ', '.join(map(lambda offset: str(offset + 1),
+                               sorted(aaMatch[key]['gapOffsets'])))))
+        extraCount = aaMatch[key]['extraCount']
+        if extraCount:
+            append(countPrint('%s    Extra nucleotides at end' % indent,
+                              extraCount, length))
+
+    return '\n'.join(result)
+
+
+def compareAaReads(read1, read2, gapChars='-', offsets=None):
+    """
+    Compare two amino acid sequences.
+
+    @param read1: A C{Read} instance or an instance of one of its subclasses.
+    @param read2: A C{Read} instance or an instance of one of its subclasses.
+    @param gapChars: An object supporting __contains__ with characters that
+        should be considered to be gaps.
+    @param offsets: If not C{None}, a C{set} of offsets of interest. Offsets
+        not in the set will not be considered.
+    @return: A C{dict} with information about the match and the individual
+        sequences (see below).
+    """
+    matchCount = 0
+    gapMismatchCount = nonGapMismatchCount = gapGapMismatchCount = 0
+    read1ExtraCount = read2ExtraCount = 0
+    read1GapOffsets = []
+    read2GapOffsets = []
+
+    for offset, (a, b) in enumerate(zip_longest(read1.sequence.upper(),
+                                                read2.sequence.upper())):
+        # Use 'is not None' in the following to allow an empty offsets set
+        # to be passed.
+        if offsets is not None and offset not in offsets:
+            continue
+        if a is None:
+            # b has an extra character at its end (it cannot be None).
+            assert b is not None
+            read2ExtraCount += 1
+            if b in gapChars:
+                read2GapOffsets.append(offset)
+        elif b is None:
+            # a has an extra character at its end.
+            read1ExtraCount += 1
+            if a in gapChars:
+                read1GapOffsets.append(offset)
+        else:
+            # We have a character from both sequences (they could still be
+            # gap characters).
+            if a in gapChars:
+                read1GapOffsets.append(offset)
+                if b in gapChars:
+                    # Both are gaps. This can happen (though hopefully not
+                    # if the sequences were pairwise aligned).
+                    gapGapMismatchCount += 1
+                    read2GapOffsets.append(offset)
+                else:
+                    # a is a gap, b is not.
+                    gapMismatchCount += 1
+            else:
+                if b in gapChars:
+                    # b is a gap, a is not.
+                    gapMismatchCount += 1
+                    read2GapOffsets.append(offset)
+                else:
+                    # Neither is a gap character.
+                    if a == b:
+                        matchCount += 1
+                    else:
+                        nonGapMismatchCount += 1
+
+    return {
+        'match': {
+            'matchCount': matchCount,
+            'gapMismatchCount': gapMismatchCount,
+            'gapGapMismatchCount': gapGapMismatchCount,
+            'nonGapMismatchCount': nonGapMismatchCount,
+        },
+        'read1': {
+            'extraCount': read1ExtraCount,
+            'gapOffsets': read1GapOffsets,
+        },
+        'read2': {
+            'extraCount': read2ExtraCount,
+            'gapOffsets': read2GapOffsets,
+        },
+    }
diff --git a/dark/dna.py b/dark/dna.py
index f189294..327328e 100644
--- a/dark/dna.py
+++ b/dark/dna.py
@@ -1,3 +1,4 @@
+from dark.utils import countPrint
 try:
     from itertools import zip_longest
 except ImportError:
@@ -32,44 +33,6 @@ BASES_TO_AMBIGUOUS = dict(
     (''.join(sorted(bases)), symbol) for symbol, bases in AMBIGUOUS.items())
 
 
-def _pct(a, b):
-    """
-    What percent of a is b?
-
-    @param a: a numeric value.
-    @param b: a numeric value.
-    @return: the C{float} percentage.
-    """
-    return 100.0 * a / b if b else 0.0
-
-
-def _pp(mesg, count, len1, len2=None):
-    """
-    Format a message followed by an integer count and a percentage (or
-    two, if the sequence lengths are unequal).
-
-    @param mesg: a C{str} message.
-    @param count: a numeric value.
-    @param len1: the C{int} length of sequence 1.
-    @param len2: the C{int} length of sequence 2. If C{None}, will
-        default to C{len1}.
-    @return: A C{str} for printing.
-    """
-    if count == 0:
-        return '%s: %d' % (mesg, count)
-    else:
-        len2 = len2 or len1
-        if len1 == len2:
-            return '%s: %d/%d (%.2f%%)' % (
-                mesg, count, len1, _pct(count, len1))
-        else:
-            return ('%s: %d/%d (%.2f%%) of sequence 1, '
-                    '%d/%d (%.2f%%) of sequence 2)' % (
-                        mesg,
-                        count, len1, _pct(count, len1),
-                        count, len2, _pct(count, len2)))
-
-
 def matchToString(dnaMatch, read1, read2, matchAmbiguous=True, indent='',
                   offsets=None):
     """
@@ -102,30 +65,32 @@ def matchToString(dnaMatch, read1, read2, matchAmbiguous=True, indent='',
     result = []
     append = result.append
 
-    append(_pp('%sExact matches' % indent, identicalMatchCount, len1, len2))
-    append(_pp('%sAmbiguous matches' % indent, ambiguousMatchCount, len1,
-               len2))
+    append(countPrint('%sExact matches' % indent, identicalMatchCount,
+                      len1, len2))
+    append(countPrint('%sAmbiguous matches' % indent, ambiguousMatchCount,
+                      len1, len2))
     if ambiguousMatchCount and identicalMatchCount:
         anyMatchCount = identicalMatchCount + ambiguousMatchCount
-        append(_pp('%sExact or ambiguous matches' % indent, anyMatchCount,
-                   len1, len2))
+        append(countPrint('%sExact or ambiguous matches' % indent,
+                          anyMatchCount, len1, len2))
+
     mismatchCount = (gapMismatchCount + gapGapMismatchCount +
                      nonGapMismatchCount)
-    append(_pp('%sMismatches' % indent, mismatchCount, len1, len2))
+    append(countPrint('%sMismatches' % indent, mismatchCount, len1, len2))
     conflicts = 'conflicts' if matchAmbiguous else 'conflicts or ambiguities'
-    append(_pp('%s  Not involving gaps (i.e., %s)' % (indent, conflicts),
-               nonGapMismatchCount, len1, len2))
-    append(_pp('%s  Involving a gap in one sequence' % indent,
-               gapMismatchCount, len1, len2))
-    append(_pp('%s  Involving a gap in both sequences' % indent,
-               gapGapMismatchCount, len1, len2))
+    append(countPrint('%s  Not involving gaps (i.e., %s)' % (indent,
+                      conflicts), nonGapMismatchCount, len1, len2))
+    append(countPrint('%s  Involving a gap in one sequence' % indent,
+                      gapMismatchCount, len1, len2))
+    append(countPrint('%s  Involving a gap in both sequences' % indent,
+                      gapGapMismatchCount, len1, len2))
 
     for read, key in zip((read1, read2), ('read1', 'read2')):
         append('%s  Id: %s' % (indent, read.id))
         length = len(read)
         append('%s    Length: %d' % (indent, length))
         gapCount = len(dnaMatch[key]['gapOffsets'])
-        append(_pp('%s    Gaps' % indent, gapCount, length))
+        append(countPrint('%s    Gaps' % indent, gapCount, length))
         if gapCount:
             append(
                 '%s    Gap locations (1-based): %s' %
@@ -133,11 +98,11 @@ def matchToString(dnaMatch, read1, read2, matchAmbiguous=True, indent='',
                  ', '.join(map(lambda offset: str(offset + 1),
                                sorted(dnaMatch[key]['gapOffsets'])))))
         ambiguousCount = len(dnaMatch[key]['ambiguousOffsets'])
-        append(_pp('%s    Ambiguous' % indent, ambiguousCount, length))
+        append(countPrint('%s    Ambiguous' % indent, ambiguousCount, length))
         extraCount = dnaMatch[key]['extraCount']
         if extraCount:
-            append(_pp('%s    Extra nucleotides at end' % indent, extraCount,
-                       length))
+            append(countPrint('%s    Extra nucleotides at end' % indent,
+                              extraCount, length))
 
     return '\n'.join(result)
 
diff --git a/dark/utils.py b/dark/utils.py
index 4862cba..810088f 100644
--- a/dark/utils.py
+++ b/dark/utils.py
@@ -214,3 +214,40 @@ def nucleotidesToStr(nucleotides, prefix=''):
             '%s%d: %s' % (prefix, offset,
                           baseCountsToStr(nucleotides[offset])))
     return '\n'.join(result)
+
+
+def countPrint(mesg, count, len1, len2=None):
+    """
+    Format a message followed by an integer count and a percentage (or
+    two, if the sequence lengths are unequal).
+
+    @param mesg: a C{str} message.
+    @param count: a numeric value.
+    @param len1: the C{int} length of sequence 1.
+    @param len2: the C{int} length of sequence 2. If C{None}, will
+        default to C{len1}.
+    @return: A C{str} for printing.
+    """
+    def percentage(a, b):
+        """
+        What percent of a is b?
+
+        @param a: a numeric value.
+        @param b: a numeric value.
+        @return: the C{float} percentage.
+        """
+        return 100.0 * a / b if b else 0.0
+
+    if count == 0:
+        return '%s: %d' % (mesg, count)
+    else:
+        len2 = len2 or len1
+        if len1 == len2:
+            return '%s: %d/%d (%.2f%%)' % (
+                mesg, count, len1, percentage(count, len1))
+        else:
+            return ('%s: %d/%d (%.2f%%) of sequence 1, '
+                    '%d/%d (%.2f%%) of sequence 2' % (
+                        mesg,
+                        count, len1, percentage(count, len1),
+                        count, len2, percentage(count, len2)))


Test patch:
diff --git a/test/test_aa.py b/test/test_aa.py
index fed28c2..43812f5 100644
--- a/test/test_aa.py
+++ b/test/test_aa.py
@@ -8,7 +8,8 @@ from dark.aa import (
     HYDROXYLIC, NEGATIVE, NONE, POLAR, SMALL, SULPHUR, TINY, NAMES,
     NAMES_TO_ABBREV1, ABBREV3, ABBREV3_TO_ABBREV1, CODONS, AA_LETTERS,
     find, AminoAcid, clustersForSequence, propertiesForSequence,
-    PROPERTY_CLUSTERS, PROPERTY_DETAILS_RAW, START_CODON, STOP_CODONS)
+    PROPERTY_CLUSTERS, PROPERTY_DETAILS_RAW, START_CODON, STOP_CODONS,
+    compareAaReads, matchToString)
 from dark.reads import AARead
 
 
@@ -822,3 +823,333 @@ class TestClustersForSequence(TestCase):
             },
             clustersForSequence(read, ['composition', 'hydropathy'],
                                 missingAAValue=10))
+
+
+class TestCompareAaReads(TestCase):
+    """
+    Test the compareAaReads function.
+    """
+    def testEmptySequences(self):
+        """
+        Two empty sequences must compare as expected.
+        """
+        self.assertEqual(
+            {
+                'match': {
+                    'matchCount': 0,
+                    'gapMismatchCount': 0,
+                    'gapGapMismatchCount': 0,
+                    'nonGapMismatchCount': 0,
+                },
+                'read1': {
+                    'extraCount': 0,
+                    'gapOffsets': [],
+                },
+                'read2': {
+                    'extraCount': 0,
+                    'gapOffsets': [],
+                },
+            },
+            compareAaReads(AARead('id1', ''),
+                           AARead('id2', '')))
+
+    def testExactMatch(self):
+        """
+        Two sequences that match exactly must compare as expected.
+        """
+        self.assertEqual(
+            {
+                'match': {
+                    'matchCount': 5,
+                    'gapMismatchCount': 0,
+                    'gapGapMismatchCount': 0,
+                    'nonGapMismatchCount': 0,
+                },
+                'read1': {
+                    'extraCount': 0,
+                    'gapOffsets': [],
+                },
+                'read2': {
+                    'extraCount': 0,
+                    'gapOffsets': [],
+                },
+            },
+            compareAaReads(AARead('id1', 'GALHN'),
+                           AARead('id2', 'GALHN')))
+
+    def testMismatch(self):
+        """
+        If the sequences have mismatched bases, their count
+        must be given correctly in the nonGapMismatchCount.
+        """
+        self.assertEqual(
+            {
+                'match': {
+                    'matchCount': 3,
+                    'gapMismatchCount': 0,
+                    'gapGapMismatchCount': 0,
+                    'nonGapMismatchCount': 2,
+                },
+                'read1': {
+                    'extraCount': 0,
+                    'gapOffsets': [],
+                },
+                'read2': {
+                    'extraCount': 0,
+                    'gapOffsets': [],
+                },
+            },
+            compareAaReads(AARead('id1', 'GALYY'),
+                           AARead('id2', 'GALHN')))
+
+    def testNoOffsets(self):
+        """
+        If an empty set of wanted offsets is passed, the result must be empty.
+        """
+        self.assertEqual(
+            {
+                'match': {
+                    'matchCount': 0,
+                    'gapMismatchCount': 0,
+                    'gapGapMismatchCount': 0,
+                    'nonGapMismatchCount': 0,
+                },
+                'read1': {
+                    'extraCount': 0,
+                    'gapOffsets': [],
+                },
+                'read2': {
+                    'extraCount': 0,
+                    'gapOffsets': [],
+                },
+            },
+            compareAaReads(AARead('id1', 'GAL-N'),
+                           AARead('id2', 'G-LHN'), offsets=set()))
+
+    def testOffsets(self):
+        """
+        If a set of wanted offsets is passed, the result must be restricted to
+        just those offsets.
+        """
+        self.assertEqual(
+            {
+                'match': {
+                    'matchCount': 1,
+                    'gapMismatchCount': 0,
+                    'gapGapMismatchCount': 0,
+                    'nonGapMismatchCount': 1,
+                },
+                'read1': {
+                    'extraCount': 0,
+                    'gapOffsets': [],
+                },
+                'read2': {
+                    'extraCount': 0,
+                    'gapOffsets': [],
+                },
+            },
+            compareAaReads(AARead('id1', 'GAL-L'),
+                           AARead('id2', 'G-LHN'), offsets=set([0, 4])))
+
+    def testGapInFirst(self):
+        """
+        A gap in the first sequence must be dealt with correctly.
+        """
+        self.assertEqual(
+            {
+                'match': {
+                    'matchCount': 4,
+                    'gapMismatchCount': 1,
+                    'gapGapMismatchCount': 0,
+                    'nonGapMismatchCount': 0,
+                },
+                'read1': {
+                    'extraCount': 0,
+                    'gapOffsets': [3],
+                },
+                'read2': {
+                    'extraCount': 0,
+                    'gapOffsets': [],
+                },
+            },
+            compareAaReads(AARead('id1', 'GAL-N'),
+                           AARead('id2', 'GALHN')))
+
+    def testGapInSecond(self):
+        """
+        A gap in the second sequence must be dealt with correctly.
+        """
+        self.assertEqual(
+            {
+                'match': {
+                    'matchCount': 3,
+                    'gapMismatchCount': 2,
+                    'gapGapMismatchCount': 0,
+                    'nonGapMismatchCount': 0,
+                },
+                'read1': {
+                    'extraCount': 0,
+                    'gapOffsets': [],
+                },
+                'read2': {
+                    'extraCount': 0,
+                    'gapOffsets': [1, 2],
+                },
+            },
+            compareAaReads(AARead('id1', 'GALHN'),
+                           AARead('id2', 'G--HN')))
+
+    def testNonDefaultGapChars(self):
+        """
+        We must be able to specify the gap characters.
+        """
+        for gap in '+$':
+            self.assertEqual(
+                {
+                    'match': {
+                        'matchCount': 3,
+                        'gapMismatchCount': 2,
+                        'gapGapMismatchCount': 0,
+                        'nonGapMismatchCount': 0,
+                    },
+                    'read1': {
+                        'extraCount': 0,
+                        'gapOffsets': [2],
+                    },
+                    'read2': {
+                        'extraCount': 0,
+                        'gapOffsets': [0],
+                    },
+                },
+                compareAaReads(AARead('id1', 'GA%sHN' % gap),
+                               AARead('id2', '%sALHN' % gap), gapChars='+$'))
+
+    def testGapGap(self):
+        """
+        Coinciding gaps in the sequences must be dealt with correctly
+        """
+        self.assertEqual(
+            {
+                'match': {
+                    'matchCount': 2,
+                    'gapMismatchCount': 2,
+                    'gapGapMismatchCount': 1,
+                    'nonGapMismatchCount': 0,
+                },
+                'read1': {
+                    'extraCount': 0,
+                    'gapOffsets': [2, 3],
+                },
+                'read2': {
+                    'extraCount': 0,
+                    'gapOffsets': [1, 2],
+                },
+            },
+            compareAaReads(AARead('id1', 'GA--N'),
+                           AARead('id2', 'G--HN')))
+
+    def testExtraInFirst(self):
+        """
+        If the first sequence has extra bases, they must be indicated in the
+        extraCount.
+        """
+        self.assertEqual(
+            {
+                'match': {
+                    'matchCount': 5,
+                    'gapMismatchCount': 0,
+                    'gapGapMismatchCount': 0,
+                    'nonGapMismatchCount': 0,
+                },
+                'read1': {
+                    'extraCount': 2,
+                    'gapOffsets': [],
+                },
+                'read2': {
+                    'extraCount': 0,
+                    'gapOffsets': [],
+                },
+            },
+            compareAaReads(AARead('id1', 'GALHNHN'),
+                           AARead('id2', 'GALHN')))
+
+    def testExtraInSecond(self):
+        """
+        If the second sequence has extra bases, they must be indicated in the
+        extraCount.
+        """
+        self.assertEqual(
+            {
+                'match': {
+                    'matchCount': 5,
+                    'gapMismatchCount': 0,
+                    'gapGapMismatchCount': 0,
+                    'nonGapMismatchCount': 0,
+                },
+                'read1': {
+                    'extraCount': 0,
+                    'gapOffsets': [],
+                },
+                'read2': {
+                    'extraCount': 2,
+                    'gapOffsets': [],
+                },
+            },
+            compareAaReads(AARead('id1', 'GALHN'),
+                           AARead('id2', 'GALHNHN')))
+
+
+class TestMatchToString(TestCase):
+    """
+    Test the matchToString function.
+    """
+    def testMismatchAndMatch(self):
+        """
+        Two sequences containing matches and mismatches must compare as
+        expected.
+        """
+        read1 = AARead('id1', 'GALHNG')
+        read2 = AARead('id2', 'GALHNA')
+        match = compareAaReads(read1, read2)
+
+        self.assertEqual(
+            '''\
+Matches: 5/6 (83.33%)
+Mismatches: 1/6 (16.67%)
+  Not involving gaps (i.e., conflicts): 1/6 (16.67%)
+  Involving a gap in one sequence: 0
+  Involving a gap in both sequences: 0
+  Id: id1
+    Length: 6
+    Gaps: 0
+  Id: id2
+    Length: 6
+    Gaps: 0''',
+            matchToString(match, read1, read2)
+        )
+
+    def testGapAndMatch(self):
+        """
+        Two sequences containing matches and gaps must compare as
+        expected.
+        """
+        read1 = AARead('id1', 'GALHN-')
+        read2 = AARead('id2', 'GALHNA')
+        match = compareAaReads(read1, read2)
+
+        self.assertEqual(
+            '''\
+Matches: 5/6 (83.33%)
+Mismatches: 1/6 (16.67%)
+  Not involving gaps (i.e., conflicts): 0
+  Involving a gap in one sequence: 1/6 (16.67%)
+  Involving a gap in both sequences: 0
+  Id: id1
+    Length: 6
+    Gaps: 1/6 (16.67%)
+    Gap locations (1-based): 6
+  Id: id2
+    Length: 6
+    Gaps: 0''',
+            matchToString(match, read1, read2)
+        )
diff --git a/test/test_utils.py b/test/test_utils.py
index f1db1b1..38be676 100644
--- a/test/test_utils.py
+++ b/test/test_utils.py
@@ -15,7 +15,7 @@ from .mocking import mockOpen, File
 
 from dark.utils import (
     numericallySortFilenames, median, asHandle, parseRangeString, StringIO,
-    baseCountsToStr, nucleotidesToStr)
+    baseCountsToStr, nucleotidesToStr, countPrint)
 
 
 class TestNumericallySortFilenames(TestCase):
@@ -355,3 +355,30 @@ class TestNucleotidesToStr(TestCase):
                 }
             )
         )
+
+
+class TestCountPrint(TestCase):
+    """
+    Test the countPrint function and the contained percentage function.
+    """
+    def testSimple(self):
+        """
+        A simple example must work as expected.
+        """
+        count = 2
+        len1 = 10
+        self.assertEqual('Count is: 2/10 (20.00%)',
+                         countPrint('Count is', count, len1))
+
+    def testTwoSequences(self):
+        """
+        An example involving two different lengths must work as expected.
+        """
+        count = 2
+        len1 = 10
+        len2 = 8
+        self.assertEqual(
+            'Count is: 2/10 (20.00%) of sequence 1,'
+            ' 2/8 (25.00%) of sequence 2',
+            countPrint('Count is', count, len1, len2)
+        )