Repository: acorg/dark-matter

Issue:
Add a script to make a protein database

Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 67f29d6..9e16187 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,11 @@
+## 3.1.0 Aug 31, 2019
+
+Added `make-protein-database.py`, `download-genbank.sh`, and
+`parse-genbank-flat-file.py` scripts, as well as `doc/protein-database.md`
+with some instructions on how to make a protein database. Added CIV
+(Charite Institute of Virology) scripts `proteins-to-pathogens-civ.py` and
+`alignment-panel-civ.py`.
+
 ## 3.0.80 Aug 19, 2019
 
 Added `bin/create-newick-relabeling-output.py` and  `bin/relabel-newick-tree.py`.
@@ -60,6 +68,7 @@ Added `--title` and `--preamble` args to output from
 `proteins-to-pathogens.py`.  Fixed `ProteinGrouper` HTML NCBI protein link
 and added genome link. Added positive and negative filtering by regex to
 `TitlesAlignments` and tests. Improved NCBI link generation and tests.
+>>>>>>> master
 
 ## 3.0.68 Jun 9, 2019
 
@@ -77,6 +86,7 @@ Added `pysam` to `setup.py` `install_requires` list. Removed `cffi`. Fixed
 tests that were failing under Linux (apart from pyfaidx tests which are now
 skipped on Linux). Removed mocking `File` class and replaced it with
 `StringIO`.
+>>>>>>> master
 
 ## 3.0.65 Jan 31, 2019
 
diff --git a/bin/aa-to-dna.py b/bin/aa-to-dna.py
new file mode 100755
index 0000000..1e27cb6
--- /dev/null
+++ b/bin/aa-to-dna.py
@@ -0,0 +1,38 @@
+#!/usr/bin/env python
+
+from __future__ import print_function
+
+from random import choice
+from operator import itemgetter
+
+from dark.aa import CODONS
+from dark.reads import (
+    Read, addFASTACommandLineOptions, parseFASTACommandLineOptions)
+
+
+if __name__ == '__main__':
+    import argparse
+
+    parser = argparse.ArgumentParser(
+        description=('Given amino acid FASTA on stdin, write FASTA with DNA '
+                     'sequences to stdout, where the DNA sequences would '
+                     'translate to the given amino acid seqeunces. The DNA '
+                     'bases used for each amino acid is either fixed or is '
+                     'chosen at random from the corresponding codons.'))
+
+    parser.add_argument(
+        '--random', action='store_true',
+        help=('Use codons chosen at random for each amino acid. Otherwise a '
+              'fixed codon will be used for each amino acid so that the '
+              'output for a given amino acid sequence is always identical.'))
+
+    addFASTACommandLineOptions(parser)
+    args = parser.parse_args()
+    reads = parseFASTACommandLineOptions(args)
+
+    getCodon = choice if args.random else itemgetter(0)
+
+    for aaRead in reads:
+        dnaRead = Read(aaRead.id,
+                       ''.join(getCodon(CODONS[aa]) for aa in aaRead.sequence))
+        print(dnaRead.toString('fasta'), end='')
diff --git a/bin/accession-to-taxonomy.py b/bin/accession-to-taxonomy.py
new file mode 100755
index 0000000..5529616
--- /dev/null
+++ b/bin/accession-to-taxonomy.py
@@ -0,0 +1,91 @@
+#!/usr/bin/env python
+
+from __future__ import print_function
+
+import sys
+import argparse
+import re
+
+from dark.taxonomy import AccessionLineageFetcher, formatLineage
+
+VERSION_REGEX = re.compile(r'\.\d+$')
+
+
+def taxonomyInfo(accession, db, namesOnly, separator):
+    """
+    Convert the lineage information for C{accession} into a formatted string.
+
+    @param accession: The C{str} accession number to look up.
+    @param db: The sqlite3 taxonomy database to consult.
+    @param namesOnly: If C{True} only print taxonomic names.
+    @param separator: The C{str} separator to put between fields. If C{None}
+        return space-padded aligned columns.
+    @return: A formatted C{str} for printing or C{None} if C{accession} is
+        not found in the taxonomy database.
+    """
+    lineage = db.lineage(accession)
+
+    if lineage is None and VERSION_REGEX.search(accession) is None:
+        # Try adding a version number.
+        lineage = db.lineage(accession + '.1')
+
+    if lineage:
+        return formatLineage(lineage, namesOnly, separator)
+
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser(
+        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
+        description='Print taxonomy information for accession numbers')
+
+    parser.add_argument(
+        'accessions', nargs='*',
+        help=('The accession number(s) to print taxonomy information for. If '
+              'not given, accession numbers are read from standard input, one '
+              'per line.'))
+
+    parser.add_argument(
+        '--separator',
+        help=('The string to separate output columns with. Default is to '
+              'print aligned output padded with multiple spaces (or a TAB if '
+              '--namesOnly is used).'))
+
+    parser.add_argument(
+        '--database', required=True,
+        help=('The file holding the sqlite3 taxonomy database. See '
+              'https://github.com/acorg/ncbi-taxonomy-database for how to '
+              'build one.'))
+
+    parser.add_argument(
+        '--namesOnly', default=False, action='store_true',
+        help='If specified, only print the taxonomic names.')
+
+    parser.add_argument(
+        '--printAccession', default=False, action='store_true',
+        help='If specified, also print the accession number.')
+
+    args = parser.parse_args()
+
+    db = AccessionLineageFetcher(args.database)
+
+    if args.accessions:
+        accessions = args.accessions
+    else:
+        accessions = (line[:-1] for line in sys.stdin)
+
+    separator = args.separator
+
+    if args.namesOnly and separator is None:
+        separator = '\t'
+
+    printAccession = args.printAccession
+
+    for accession in accessions:
+        if printAccession:
+            print(accession + ':')
+        result = taxonomyInfo(accession, db, args.namesOnly, separator)
+        if result:
+            print(result)
+        else:
+            print('%r was not found in the taxonomy database' % accession,
+                  file=sys.stderr)
diff --git a/bin/alignment-panel-civ.py b/bin/alignment-panel-civ.py
new file mode 100755
index 0000000..03c66c3
--- /dev/null
+++ b/bin/alignment-panel-civ.py
@@ -0,0 +1,385 @@
+#!/usr/bin/env python
+
+"""
+Given a BLAST or DIAMOND JSON output files, the corresponding FASTA (or FASTQ)
+sequence files, and filtering criteria, produce a summary of matched titles
+and (optionally) an alignment panel.
+
+Run with --help for help.
+"""
+
+from __future__ import print_function
+
+import os
+import sys
+import argparse
+from json import dump
+from collections import defaultdict
+from itertools import chain
+
+# It's not clear that the PDF backend is the right choice here, but it
+# works (i.e., the generation of PNG images works fine).
+import matplotlib
+matplotlib.use('PDF')
+
+# These imports are here because dark.graphics imports matplotlib.pyplot
+# and we need to set the matplotlib backend (see above) before that import
+# happens. So please don't move these imports higher in this file.
+from dark.titles import TitlesAlignments
+from dark.fasta import FastaReads
+from dark.fastq import FastqReads
+from dark.civ.graphics import (
+    DEFAULT_LOG_LINEAR_X_AXIS_BASE, alignmentPanelHTML)
+from dark.civ.proteins import SqliteIndex
+from dark.utils import numericallySortFilenames
+
+
+def parseColors(colors, args):
+    """
+    Parse read id color specification.
+
+    @param colors: A C{list} of C{str}s. Each item is of the form, e.g.,
+        'green X Y Z...', where each of X, Y, Z, ... etc. is either a read
+        id or the name of a FASTA or FASTQ file containing reads whose ids
+        should be displayed with the corresponding color. Note that if read
+        ids contain spaces you will need to use the latter (i.e. FASTA/Q file
+        name) approach because C{args.colors} is split on whitespace.
+    @param args: The argparse C{Namespace} instance holding the other parsed
+        command line arguments.
+    @return: A C{dict} whose keys are colors and whose values are sets of
+        read ids.
+    """
+    result = defaultdict(set)
+    for colorInfo in colors:
+        readIds = colorInfo.split()
+        color = readIds.pop(0)
+        for readId in readIds:
+            if os.path.isfile(readId):
+                filename = readId
+                if args.fasta:
+                    reads = FastaReads(filename)
+                else:
+                    reads = FastqReads(filename)
+                for read in reads:
+                    result[color].add(read.id)
+            else:
+                result[color].add(readId)
+    return result
+
+
+if __name__ == '__main__':
+
+    # We do not use the addFASTACommandLineOptions and
+    # parseFASTACommandLineOptions utility functions below because we allow
+    # multiple FASTA or FASTQ files on the command line, which we specify
+    # by --fasta and --fastq. And those names clash with the option names
+    # used by those utility functions.
+
+    parser = argparse.ArgumentParser(
+        description='Produce an alignment panel in HTML',
+        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
+        epilog=('Given BLAST or DIAMOND JSON output files, the '
+                'corresponding FASTA (or FASTQ) sequence files, and '
+                'filtering criteria, produce an alignment panel.'))
+
+    parser.add_argument(
+        '--proteinGenomeDatabase', required=True,
+        help=('The filename of an Sqlite3 database holding protein and '
+              'genome information, as built by make-protein-database.py'))
+
+    parser.add_argument(
+        '--earlyExit', default=False, action='store_true',
+        help=('If True, just print the number of interesting matches, but do '
+              'not create the alignment panel.'))
+
+    parser.add_argument(
+        '--matcher', default='blast', choices=('blast', 'diamond'),
+        help='The matching algorithm that was used to produce the JSON.')
+
+    parser.add_argument(
+        '--json', metavar='JSON-file', nargs='+', action='append',
+        required=True, help='the JSON file(s) of BLAST or DIAMOND output.')
+
+    # A mutually exclusive group for either FASTA or FASTQ files.
+    group = parser.add_mutually_exclusive_group(required=True)
+
+    group.add_argument(
+        '--fasta', metavar='FASTA-file', nargs='+', action='append',
+        help=('the FASTA file(s) of read sequences that were given to BLAST '
+              'or DIAMOND.'))
+
+    group.add_argument(
+        '--fastq', metavar='FASTQ-file', nargs='+', action='append',
+        help=('the FASTQ file(s) of read sequences that were given to BLAST '
+              'or DIAMOND.'))
+
+    # Args for filtering on ReadsAlignments.
+    parser.add_argument(
+        '--minStart', type=int, default=None,
+        help='Reads that start before this subject offset should not be '
+        'shown.')
+
+    parser.add_argument(
+        '--maxStop', type=int, default=None,
+        help='Reads that end after this subject offset should not be shown.')
+
+    parser.add_argument(
+        '--oneAlignmentPerRead', default=False, action='store_true',
+        help='If True, only keep the best alignment for each read.')
+
+    parser.add_argument(
+        '--maxAlignmentsPerRead', type=int, default=None,
+        help=('Reads with more than this many alignments will be elided. Pass '
+              'zero to only keep reads with no matches (alignments).'))
+
+    parser.add_argument(
+        '--scoreCutoff', type=float, default=None,
+        help=('A float score. Matches with scores worse than this will be '
+              'ignored.'))
+
+    parser.add_argument(
+        '--maxHspsPerHit', type=int, default=None,
+        help='A numeric max number of HSPs to show for each hit on hitId.')
+
+    parser.add_argument(
+        '--whitelist', nargs='+', default=None, action='append',
+        help='sequence titles that should be whitelisted')
+
+    parser.add_argument(
+        '--whitelistFile', default=None,
+        help=('A file containing sequence titles (one per line) to be '
+              'whitelisted'))
+
+    parser.add_argument(
+        '--blacklist', nargs='+', default=None, action='append',
+        help='sequence titles that should be blacklisted')
+
+    parser.add_argument(
+        '--blacklistFile', default=None,
+        help=('A file containing sequence titles (one per line) to be '
+              'blacklisted'))
+
+    parser.add_argument(
+        '--titleRegex', default=None,
+        help='a regex that sequence titles must match.')
+
+    parser.add_argument(
+        '--negativeTitleRegex', default=None,
+        help='a regex that sequence titles must not match.')
+
+    parser.add_argument(
+        '--truncateTitlesAfter', default=None,
+        help=('a string that titles will be truncated beyond. If the '
+              'truncated version of a title has already been seen, '
+              'that title will be skipped.'))
+
+    parser.add_argument(
+        '--minSequenceLen', type=int, default=None,
+        help='sequences of lesser length will be elided.')
+
+    parser.add_argument(
+        '--maxSequenceLen', type=int, default=None,
+        help='sequences of greater length will be elided.')
+
+    parser.add_argument(
+        '--taxonomy', default=None,
+        help=('a string of the taxonomic group on which should be '
+              'filtered. eg "Vira" will filter on viruses.'))
+
+    # Args for filtering on TitlesAlignments.
+    parser.add_argument(
+        '--minMatchingReads', type=int, default=None,
+        help='sequences that are matched by fewer reads will be elided.')
+
+    parser.add_argument(
+        '--minMedianScore', type=float, default=None,
+        help=('sequences that are matched with a median score that is '
+              'worse will be elided.'))
+
+    parser.add_argument(
+        '--withScoreBetterThan', type=float, default=None,
+        help=('sequences that are matched without at least one score '
+              'at least this good will be elided.'))
+
+    parser.add_argument(
+        '--minNewReads', type=float, default=None,
+        help=('The fraction of its reads by which a new read set must differ '
+              'from all previously seen read sets in order to be considered '
+              'acceptably different.'))
+
+    parser.add_argument(
+        '--maxTitles', type=int, default=None,
+        help=('The maximum number of titles to keep. If more titles than '
+              'this result from the filtering, titles will be sorted '
+              '(according to the --sortOn value) and only the best will be '
+              'retained.'))
+
+    parser.add_argument(
+        '--minCoverage', type=float, default=None,
+        help=('The (0.0 to 1.0) minimum fraction of a subject sequence that '
+              'must be matched by at least one read.'))
+
+    # Args for the alignment panel
+    parser.add_argument(
+        '--sortOn', default='maxScore',
+        choices=('maxScore', 'medianScore', 'readCount', 'length', 'title'),
+        help='The attribute to sort subplots on.')
+
+    parser.add_argument(
+        '--rankValues', type=bool, default=False,
+        help=('If True, display reads with a Y axis coord that is the rank of '
+              'the score.'))
+
+    parser.add_argument(
+        '--outputDir', default=None, required=True,
+        help='Specifies a directory to write the HTML summary to. Required.')
+
+    parser.add_argument(
+        '--color', action='append',
+        help=('a string which has a color as the first element and read ids '
+              'and/or FASTA file names as the following element(s), separated '
+              'by spaces.'))
+
+    parser.add_argument(
+        '--equalizeXAxes', default=False, action='store_true',
+        help=('If True, all alignment graphs will have their X axes drawn '
+              'with the same range.'))
+
+    parser.add_argument(
+        '--xRange', default='subject',
+        choices=('reads', 'subject'),
+        help=('Set the X axis range to show either the subject or the extent '
+              'of the reads that hit the subject.'))
+
+    parser.add_argument(
+        '--logLinearXAxis', default=False, action='store_true',
+        help=('If True, convert read offsets so that empty regions in the '
+              'alignment panel plots will only be as wide as their logged '
+              'actual values'))
+
+    parser.add_argument(
+        '--logBase', type=float, default=DEFAULT_LOG_LINEAR_X_AXIS_BASE,
+        help='The base of the logarithm to use if logLinearXAxis is True')
+
+    parser.add_argument(
+        '--showFeatures', default=False, action='store_true',
+        help=('If specified, look up features for the individual images in '
+              'the alignment panel.'))
+
+    parser.add_argument(
+        '--sortFilenames', default=False, action='store_true',
+        help=('If specified, the JSON and FASTA/Q file names will be '
+              'processed in sorted order. The sorting is based on finding '
+              'a numerical prefix in the filename. This can be useful when '
+              'processing output files produced by systems like HTCondor, '
+              'which makes files with names like 1.out, 10.out, etc. that do '
+              'not sort properly and so cannot conveniently be given to this '
+              'program along with a single FASTA/Q file (because the order of '
+              'the results in the files from HTCondor does not match the '
+              'order of sequences in the FASTA/Q file.'))
+
+    parser.add_argument(
+        '--titlesJSONFile',
+        help=('Give a file name for JSON holding information about titles to '
+              'be written to'))
+
+    args = parser.parse_args()
+
+    # Flatten lists of lists that we get from using both nargs='+' and
+    # action='append'. We use both because it allows people to use (e.g.)
+    # --json on the command line either via "--json file1 --json file2" or
+    # "--json file1 file2", or a combination of these. That way it's not
+    # necessary to remember which way you're supposed to use it and you also
+    # can't be hit by the subtle problem encountered in
+    # https://github.com/acorg/dark-matter/issues/453
+    jsonFiles = list(chain.from_iterable(args.json))
+    whitelist = (
+        set(chain.from_iterable(args.whitelist)) if args.whitelist else None)
+    blacklist = (
+        set(chain.from_iterable(args.blacklist)) if args.blacklist else None)
+
+    # TODO: Add a --readClass command-line option in case we want to
+    # process FASTA containing AA sequences.
+    if args.fasta:
+        if args.sortFilenames:
+            files = numericallySortFilenames(chain.from_iterable(args.fasta))
+        else:
+            files = list(chain.from_iterable(args.fasta))
+        reads = FastaReads(files)
+    else:
+        if args.sortFilenames:
+            files = numericallySortFilenames(chain.from_iterable(args.fastq))
+        else:
+            files = list(chain.from_iterable(args.fastq))
+        reads = FastqReads(files)
+
+    if args.matcher == 'blast':
+        from dark.blast.alignments import BlastReadsAlignments
+        readsAlignments = BlastReadsAlignments(
+            reads, jsonFiles,
+            sqliteDatabaseFilename=args.sqliteDatabaseFilename,
+            sortBlastFilenames=args.sortFilenames)
+    else:
+        # Must be 'diamond' (due to parser.add_argument 'choices' argument).
+        assert args.matcher == 'diamond'
+        from dark.diamond.alignments import DiamondReadsAlignments
+        readsAlignments = DiamondReadsAlignments(
+            reads, jsonFiles, sortFilenames=args.sortFilenames)
+
+    readsAlignments.filter(
+        maxAlignmentsPerRead=args.maxAlignmentsPerRead,
+        minSequenceLen=args.minSequenceLen,
+        maxSequenceLen=args.maxSequenceLen,
+        minStart=args.minStart, maxStop=args.maxStop,
+        oneAlignmentPerRead=args.oneAlignmentPerRead,
+        maxHspsPerHit=args.maxHspsPerHit,
+        scoreCutoff=args.scoreCutoff,
+        whitelist=whitelist, blacklist=blacklist,
+        whitelistFile=args.whitelistFile, blacklistFile=args.blacklistFile,
+        titleRegex=args.titleRegex, negativeTitleRegex=args.negativeTitleRegex,
+        truncateTitlesAfter=args.truncateTitlesAfter, taxonomy=args.taxonomy)
+
+    titlesAlignments = TitlesAlignments(readsAlignments).filter(
+        minMatchingReads=args.minMatchingReads,
+        minMedianScore=args.minMedianScore,
+        withScoreBetterThan=args.withScoreBetterThan,
+        minNewReads=args.minNewReads, maxTitles=args.maxTitles,
+        sortOn=args.sortOn, minCoverage=args.minCoverage)
+
+    nTitles = len(titlesAlignments)
+    print('Found %d interesting title%s.' %
+          (nTitles, '' if nTitles == 1 else 's'), file=sys.stderr)
+
+    # If we've been asked to save the matched title information in JSON,
+    # write it out.
+    if args.titlesJSONFile:
+        with open(args.titlesJSONFile, 'w') as fp:
+            dump(titlesAlignments.toDict(), fp, sort_keys=True, indent=2,
+                 separators=(',', ': '))
+
+    if nTitles:
+        # Use flush=True on the print of the titles so the output is
+        # definitely written out. This is because in some cases the
+        # following aligment panel generation takes a long time and exceeds
+        # its CPU or memory limit on an HPCS system.  Also, it can be good
+        # to see the matching title details while waiting for the panel
+        # plots to be generated.
+        print(titlesAlignments.tabSeparatedSummary(sortOn=args.sortOn),
+              flush=True)
+
+    if args.earlyExit:
+        sys.exit(0)
+
+    if nTitles == 0:
+        print('No alignment panel generated due to no matching titles.',
+              file=sys.stderr)
+        sys.exit(0)
+
+    idList = parseColors(args.color, args) if args.color else None
+
+    alignmentPanelHTML(
+        titlesAlignments, SqliteIndex(args.proteinGenomeDatabase),
+        sortOn=args.sortOn, outputDir=args.outputDir,
+        idList=idList, equalizeXAxes=args.equalizeXAxes, xRange=args.xRange,
+        logLinearXAxis=args.logLinearXAxis, logBase=args.logBase,
+        showFeatures=args.showFeatures)
diff --git a/bin/convert-diamond-to-sam.py b/bin/convert-diamond-to-sam.py
index fbab91f..46bc770 100755
--- a/bin/convert-diamond-to-sam.py
+++ b/bin/convert-diamond-to-sam.py
@@ -7,14 +7,14 @@ from __future__ import print_function, division
 
 import sys
 import argparse
-from os.path import basename
-from tempfile import TemporaryFile
+from os.path import join
+from resource import getrlimit, RLIMIT_NOFILE
 from functools import partial
 
-from dark import __version__ as VERSION
-from dark.btop import btop2cigar
-from dark.diamond.conversion import diamondTabularFormatToDicts, FIELDS
-from dark.reads import DNARead
+from dark.diamond.conversion import FIELDS
+from dark.diamond.sam import (
+    SimpleDiamondSAMWriter, PerReferenceDiamondSAMWriter)
+from dark.proteins import SqliteIndex
 
 parser = argparse.ArgumentParser(
     formatter_class=argparse.ArgumentDefaultsHelpFormatter,
@@ -42,10 +42,38 @@ parser.add_argument(
 parser.add_argument(
     '--keepDescriptions', action='store_true', default=False,
     help=('Do not discard text after the first space in query or subject '
-          'sequence ids. Note that. Note that this violates the SAM '
-          'specification, but since SAM files are TAB-separated there '
-          'is probably only a small chance this will cause any problems '
-          'downstream.'))
+          'sequence ids. Note that this violates the SAM specification, but '
+          'since SAM files are TAB-separated it may not cause problems '
+          'downstream, depending on what tools you use on the SAM.'))
+
+parser.add_argument(
+    '--genomesProteinsDatabaseFilename', required=True,
+    help=('The name of the file containing an sqlite3 database with '
+          'information on genomes and proteins, as made by '
+          'make-protein-database.py'))
+
+parser.add_argument(
+    '--outputDirectory', default='.',
+    help=('The name of the directory to create per-reference SAM files in. '
+          'This is only used when --perReferenceOutput is given, otherwise '
+          'the output is written to stdout as a single SAM file.'))
+
+parser.add_argument(
+    '--perReferenceOutput', action='store_true', default=False,
+    help=('If given, a separate SAM file is written for matches against each '
+          'reference sequence. The files will have names corresponding to '
+          'the accession number of the reference, with a .sam suffix '
+          'and will be put into the directory specified by '
+          '--outputDirectory.'))
+
+parser.add_argument(
+    '--filePointerCacheMaxSize', type=int,
+    default=getrlimit(RLIMIT_NOFILE)[0] >> 1,
+    help=('The maximum size for the file pointer cache. SAM output may '
+          'be being written to an arbitrary number of output files so we '
+          'need to make sure the OS limit on open files is not exceeded. '
+          'The default value is size half the maximum number of open files '
+          'will be used.'))
 
 args = parser.parse_args()
 
@@ -53,111 +81,44 @@ if args.printFields:
     print(FIELDS)
     sys.exit(0)
 
-idOnly = not args.keepDescriptions
-mappingQuality = args.mappingQuality
-ram = args.ram
-
-if 0 > mappingQuality > 255:
+if 0 > args.mappingQuality > 255:
     raise ValueError('Mapping quality must be between 0 and 255 (inclusive)')
 
-referenceLengths = {}
+genomesProteins = SqliteIndex(args.genomesProteinsDatabaseFilename)
 
-if ram:
-    nonHeaderLines = []
-    emit = nonHeaderLines.append
-else:
-    tf = TemporaryFile(mode='w+t', encoding='utf-8')
-    emit = partial(print, file=tf)
-
-for match in diamondTabularFormatToDicts(sys.stdin):
-    qseqid = match['qseqid'].split()[0] if idOnly else match['qseqid']
-    stitle = match['stitle'].split()[0] if idOnly else match['stitle']
-
-    # The subject length is ALWAYS in amino acids in DIAMOND.
-    referenceLengths[stitle] = 3 * match['slen']
-
-    # If the query frame is less than zero, the match was with a reverse
-    # complemented translation of the query. Put the reverse complement
-    # into the SAM output, which seems to be standard / accepted practice
-    # based on my web searches. See e.g., https://www.biostars.org/p/131891/
-    # for what Bowtie2 does and for some comments on this issue for SAM/BAM
-    # files in general.
-    if match['qframe'] > 0:
-        flag = 0
-        qseq = match['qseq']
-        qqual = match['qqual'] or '*'
-    else:
-        flag = 16
-        qseq = DNARead('id', match['qseq']).reverseComplement().sequence
-        qqual = match['qqual'][::-1] if match['qqual'] else '*'
-
-    # Make a CIGAR string, including hard-clipped bases at the start and
-    # end of the query (DIAMOND outputs a hard-clipped query sequence).
-    startClipCount = match['qstart'] - 1
-    endClipCount = match['qlen'] - match['qend']
-
-    assert startClipCount >= 0
-    assert endClipCount >= 0, (
-        'Query sequence %s has length %d but the qend value is %d' %
-        (qseq, len(match['qseq']), match['qend']))
-
-    cigar = (
-        ('%dH' % startClipCount if startClipCount else '') +
-        btop2cigar(match['btop'], concise=False, aa=True) +
-        ('%dH' % endClipCount if endClipCount else ''))
-
-    emit('\t'.join(map(str, [
-        # 1. QNAME
-        qseqid,
-        # 2. FLAG
-        flag,
-        # 3. RNAME
-        stitle,
-        # 4. POS. This needs to be a 1-based offset into the
-        # nucleotide-equivalent of the DIAMOND subject sequence (which was
-        # a protein since that is how DIAMOND operates). Because DIAMOND
-        # gives back a 1-based protein location, we adjust to 0-based,
-        # multiply by 3 to get to nucleotides, then adjust to 1-based.
-        3 * (match['sstart'] - 1) + 1,
-        # 5. MAPQ
-        mappingQuality,
-        # 6. CIGAR
-        cigar,
-        # 7. RNEXT
-        '*',
-        # 8. PNEXT
-        0,
-        # 9. TLEN
-        0,
-        # 10. SEQ
-        qseq,
-        # 11. QUAL
-        qqual,
-        # 12. Alignment score
-        'AS:i:%d' % int(match['bitscore'])])))
-
-
-progName = basename(sys.argv[0])
-
-# Print SAM headers.
-print('\n'.join(
-    [
-        '@PG\tID:DIAMOND\tPN:DIAMOND',
-        '@PG\tID:%s\tPN:%s (version %s)\tCL:%s %s\tPP:DIAMOND' %
-        (progName, progName, VERSION, progName, ' '.join(sys.argv[1:])),
-        '@CO\t%s is from the dark-matter package '
-        '(https://github.com/acorg/dark-matter/)' % progName,
-    ] +
-    [
-        '@SQ\tSN:%s\tLN:%d' % (name, referenceLengths[name])
-        for name in sorted(referenceLengths)
-    ]))
-
-# Print non-header lines.
-if ram:
-    print('\n'.join(nonHeaderLines))
+if args.perReferenceOutput:
+
+    def referenceTitleToBasename(subjectTitle):
+        """
+        Make a basename for a SAM output file for a DIAMOND subject match.
+
+        @param subjectTitle: A C{str} database subject id. This must have a
+            form like "acc|GENBANK|YP_009137150.1|GENBANK|NC_001798.2|..."
+            where the accession number of the matched sequence is in the third
+            '|'-separated field (which must not contain any characters that
+            cannot be used in a filename).
+        @return: A C{str} pathname to store the SAM for the reference to.
+        """
+        return join(args.outputDirectory,
+                    genomesProteins.genomeAccession(subjectTitle))
+
+    writer = PerReferenceDiamondSAMWriter(
+        genomesProteins, mappingQuality=args.mappingQuality,
+        ram=args.ram, keepDescriptions=args.keepDescriptions,
+        baseFilenameFunc=referenceTitleToBasename,
+        fpcMaxsize=args.filePointerCacheMaxSize)
+
+    save = writer.save
 else:
-    tf.seek(0)
-    for line in tf:
-        print(line, end='')
-    tf.close()
+    writer = SimpleDiamondSAMWriter(
+        genomesProteins, mappingQuality=args.mappingQuality,
+        ram=args.ram, keepDescriptions=args.keepDescriptions)
+
+    save = partial(writer.save, filename=sys.stdout)
+
+addMatch = writer.addMatch
+
+for line in sys.stdin:
+    addMatch(line)
+
+save()
diff --git a/bin/download-genbank.sh b/bin/download-genbank.sh
new file mode 100755
index 0000000..81fe766
--- /dev/null
+++ b/bin/download-genbank.sh
@@ -0,0 +1,170 @@
+#!/bin/bash
+
+# This script downloads flat files (in GenBank format) of sequence info
+# from NCBI. It is one of a set of utilities that we use to make a protein
+# database.
+#
+# Usage is e.g.:
+#
+#     download-genbank.sh [-n] [id-dir] < accession-numbers-file
+#
+# where:
+#
+#     accession-numbers-file is a file containing GenBank accession numbers,
+#         one per line. These should include the version number.
+#
+#     id-dir is the name of a directory in which to store files of sequence
+#         ids, parameters to cURL, and GenBank files downloaded from
+#         NCBI. If not given, defaults to the current directory.
+#
+#     -n option can be used to make the script just print out what it would
+#         download, rather than actually doing the downloads. In this case,
+#         a temporary id directory will be created and used and its name
+#         will be printed on standard output.
+#
+# See ../doc/protein-database.md for information on the various tools
+# needed to make a protein database and ../misc/Makefile-protein-database
+# for a Makefile that can do the work of calling the scripts in order with
+# a reasonable set of output file and directory names.
+
+# Note that the number of [a-z] repeats in filenames in various places in
+# the code below matches the 5 in the call to split.
+
+set -Eeuo pipefail
+
+URL=https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi
+maxIdsPerRequest=200
+
+# These are accession numbers that correspond to dead (suppressed or
+# withdrawn) sequence records (see ../doc/protein-database.md) for more
+# info. These were all in the RVDB clustered database in early 2019.
+DEAD_REGEX='JX873962|KX912841|NC_021196|NC_036584'
+
+dryRun=0
+
+while [ $# -gt 0 ]
+do
+    case "$1" in
+        -n)
+            dryRun=1
+            shift
+            ;;
+        *)
+            break
+            ;;
+    esac
+done
+
+case $# in
+    0) idDir=.;;
+    1) idDir=$1;;
+    *) echo "Usage: $(basename $0) [-n] [id-dir]" >&2; exit 1;;
+esac
+
+if [ $dryRun -eq 0 ]
+then
+    test -d $idDir || mkdir "$idDir"
+else
+    # Use a temporary id directory for a dry run.
+    idDir=$(mktemp -d)
+fi
+
+# Read stdin for accession numbers and split them into files we'll use in
+# entrez requests to NCBI.
+egrep -v "$DEAD_REGEX" | split -a 5 -l $maxIdsPerRequest - $idDir/ids_
+count=$(cat $idDir/ids_* | wc -l)
+
+# If the user has an NCBI API key, we can send requests more frequently.
+# See https://www.ncbi.nlm.nih.gov/books/NBK25497/ and man parallel.
+if [ -z "$NCBI_API_KEY" ]
+then
+    delayArg='--delay 0.35s'
+    apiParam=
+else
+    delayArg='--delay 0.1s'
+    apiParam="api_key=${NCBI_API_KEY}&"
+fi
+
+failureCount=0
+attempt=0
+attemptLimit=5
+
+# Note 5 repetitions of [a-z] here (see above comment).
+lastFile=$(ls $idDir/ids_[a-z][a-z][a-z][a-z][a-z] | tail -n 1)
+
+commandFile=$(mktemp)
+
+if [ $dryRun -eq 0 ]
+then
+    # Clean up on exit.
+    # Note 5 repetitions of [a-z] here (see above comment).
+    trap 'rm -f $commandFile $idDir/ids_[a-z][a-z][a-z][a-z][a-z] $idDir/ids_[a-z][a-z][a-z][a-z][a-z].params' 0 1 2 3 15
+fi
+
+while [ $attempt -eq 0 -o $failureCount -gt 0 ]
+do
+    attempt=$((attempt + 1))
+    
+    if [ $attempt -gt $attemptLimit ]
+    then
+        echo "Attempt limit ($attemptLimit) reached. $failureCount GenBank files still not downloaded. Exiting." 2>&1
+        break
+    elif [ $attempt -gt 1 ]
+    then
+        echo "Re-running (attempt $attempt of $attemptLimit) to retry $failureCount incompletely downloaded files." 2>&1
+    fi
+
+    failureCount=0
+
+    # Note 5 repetitions of [a-z] here (see above comment).
+    for i in $idDir/ids_[a-z][a-z][a-z][a-z][a-z]
+    do
+        out=$i.gb
+
+        if [ -f $out ]
+        then
+            # We know how many GenBank records there should be in each file.
+            if [ $i = $lastFile ]
+            then
+                expectedCountArg=$(wc -l < $i | awk '{print $1}')
+            else
+                expectedCountArg=$maxIdsPerRequest
+            fi
+
+            if parse-genbank-flat-file.py --quiet --expectedCount $expectedCountArg $out
+            then
+                # $out already exists, can be parsed, and has the right number of records. Skip.
+                continue
+            else
+                echo "$out already exists but cannot be parsed. Removing and re-downloading." >&2
+                rm -f $out
+                failureCount=$((failureCount + 1))
+            fi
+        fi
+
+        echo -n "${apiParam}db=nucleotide&rettype=gb&id=" > $i.params
+        cat $i | tr '\n' , | sed -e 's/,$//' >> $i.params
+        echo "curl --silent -d @$i.params $URL > $out"
+    done > $commandFile
+
+    if [ $dryRun -eq 0 ]
+    then
+        parallel --bar $delayArg < $commandFile
+    fi
+done
+
+if [ $dryRun -eq 1 ]
+then
+    echo "$count Genbank records would be downloaded"
+    echo "The dry run id directory was $idDir"
+    echo "Download commands are in $commandFile"
+    # Just show what we'd do. We could instead cat the command file,
+    # but that can be very verbose. Instead just let the user see the
+    # file name and they can look in it if they care.
+    echo "These would be executed using: parallel --bar $delayArg < $commandFile"
+    echo "To clean up, run: rm -r '$idDir' '$commandFile'"
+else
+    echo "$count Genbank records downloaded"
+fi
+
+exit $failureCount
diff --git a/bin/make-fasta-database.py b/bin/make-fasta-database.py
index b4ffcdc..d3b7920 100755
--- a/bin/make-fasta-database.py
+++ b/bin/make-fasta-database.py
@@ -18,7 +18,7 @@ if __name__ == '__main__':
 
     parser.add_argument(
         '--out', required=True,
-        help=('the output file. This file must not exist (use --force to '
+        help=('The output file. This file must not exist (use --force to '
               'overwrite).'))
 
     parser.add_argument(
diff --git a/bin/make-protein-database.py b/bin/make-protein-database.py
new file mode 100755
index 0000000..1dc05b5
--- /dev/null
+++ b/bin/make-protein-database.py
@@ -0,0 +1,187 @@
+#!/usr/bin/env python
+
+from __future__ import print_function, division
+
+import sys
+import warnings
+from time import time
+from itertools import chain
+import argparse
+
+from dark.civ.proteins import SqliteIndexWriter
+from dark.taxonomy import AccessionLineageFetcher
+
+
+parser = argparse.ArgumentParser(
+    formatter_class=argparse.ArgumentDefaultsHelpFormatter,
+    description=('Create a genome/protein sqlite3 database from GenBank '
+                 'and JSON files. The protein sequences for the database are '
+                 'printed to standard output for indexing by other tools.'))
+
+parser.add_argument(
+    '--databaseFile', required=True,
+    help=('The output file. This file must not exist (use --force to '
+          'overwrite).'))
+
+parser.add_argument(
+    '--databaseName',
+    help=('The database that the records in the (--gb) GenBank files came '
+          'from (e.g., "refseq" or "RVDB").'))
+
+parser.add_argument(
+    '--duplicationPolicy', choices=('error', 'ignore'), default='ignore',
+    help=('What to do if a to-be-inserted accession number is already '
+          'present in the database. "error" results in a ValueError being '
+          'raised, "ignore" will ignore the duplicate. It should also be '
+          'possible to update (i.e., replace) but that is not supported yet.'))
+
+parser.add_argument(
+    '--progress', default=False, action='store_true',
+    help='Print indexing progress.')
+
+parser.add_argument(
+    '--logFile', type=argparse.FileType('w'),
+    help='Write indexing details to a log file.')
+
+parser.add_argument(
+    '--noWarnings', default=False, action='store_true',
+    help='Do not print warnings about unparseable GenBank or JSON records.')
+
+parser.add_argument(
+    '--rnaOnly', default=False, action='store_true',
+    help='If given, only include RNA viruses.')
+
+parser.add_argument(
+    '--gb', metavar='GenBank-file', nargs='+', action='append',
+    help=('The GenBank file(s) to make the database from. These may be '
+          'uncompressed, or compressed with bgzip (from samtools), with '
+          'a .gz suffix.'))
+
+parser.add_argument(
+    '--json', metavar='JSON-file', nargs='+', action='append',
+    help=('The JSON file(s) to make the database from. These contain '
+          'genome and protein information for cases where sequences that '
+          'are not in GenBank should be added.'))
+
+parser.add_argument(
+    '--taxonomyDatabase',
+    help=('The file holding the sqlite3 taxonomy database. See '
+          'https://github.com/acorg/ncbi-taxonomy-database for how to '
+          'build one.'))
+
+parser.add_argument(
+    '--proteinSource', default='GENBANK',
+    help=('The source of the accession numbers for the proteins found in the '
+          'input files. This becomes part of the sequence id printed in the '
+          'protein FASTA output.'))
+
+parser.add_argument(
+    '--genomeSource', default='GENBANK',
+    help=('The source of the accession numbers for the genomes in the input '
+          'files. This becomes part of the sequence id printed in the '
+          'protein FASTA output.'))
+
+
+def filenamesAndAdders(args, db):
+    """
+    Get the filenames and database adding functions.
+
+    @param args: A C{Namespace} instance as returned by argparse.
+    @param db: An C{SqliteIndexWriter} instance.
+    @return: A generator yielding 2-tuples containing a filename and a
+        database adder function that can read the file and incorporate it.
+    """
+    # Use chain to flatten the lists of lists that we get from using both
+    # nargs='+' and action='append'. We use both because it allows people
+    # to use (e.g.) --gb on the command line either via "--gb file1 --gb
+    # file2" or "--gb file1 file2", or a combination of these. That way
+    # it's not necessary to remember which way you're supposed to use it
+    # and you also can't be hit by the subtle problem encountered in
+    # https://github.com/acorg/dark-matter/issues/453
+
+    if not (args.gb or args.json):
+        print('At least one of --gb or --json must be used to indicate input '
+              'files to process.', file=sys.stderr)
+        sys.exit(1)
+
+    if args.gb:
+        for filename in chain.from_iterable(args.gb):
+            yield filename, db.addGenBankFile
+
+    if args.json:
+        for filename in chain.from_iterable(args.json):
+            yield filename, db.addJSONFile
+
+
+def main(args):
+    """
+    Build the protein database.
+
+    @param args: The namespace of command-line arguments returned by
+        argparse.parse_args()
+    """
+    if args.rnaOnly:
+        if args.taxonomyDatabase:
+            lineageFetcher = AccessionLineageFetcher(
+                args.taxonomyDatabase).lineage
+        else:
+            print('If you specify --rnaOnly, you must also give a taxonomy '
+                  'database file with --taxonomyDatabase', file=sys.stderr)
+            sys.exit(1)
+    else:
+        lineageFetcher = None
+
+    progress = args.progress
+
+    if progress:
+        overallStart = time()
+        totalGenomeCount = totalProteinCount = 0
+
+    with SqliteIndexWriter(args.databaseFile) as db:
+        for fileCount, (filename, addFunc) in enumerate(
+                filenamesAndAdders(args, db), start=1):
+            if progress:
+                print("Indexing '%s' ... " % filename, end='', file=sys.stderr)
+                start = time()
+
+            genomeCount, proteinCount = addFunc(
+                filename, rnaOnly=args.rnaOnly, databaseName=args.databaseName,
+                lineageFetcher=lineageFetcher,
+                proteinSource=args.proteinSource,
+                genomeSource=args.genomeSource,
+                duplicationPolicy=args.duplicationPolicy, logfp=args.logFile)
+
+            if not genomeCount:
+                if addFunc is db.addGenBankFile:
+                    print('WARNING: no genomes added from %r. Did the GenBank '
+                          'download fail on that file?' % filename,
+                          file=sys.stderr)
+                else:
+                    print('WARNING: no genomes added from JSON file %r.' %
+                          filename, file=sys.stderr)
+
+            if progress:
+                totalGenomeCount += genomeCount
+                totalProteinCount += proteinCount
+                elapsed = time() - start
+                print('indexed %3d genome%s (%5d protein%s) in %.2f seconds.' %
+                      (genomeCount, '' if genomeCount == 1 else 's',
+                       proteinCount, '' if proteinCount == 1 else 's',
+                       elapsed), file=sys.stderr)
+
+    if progress:
+        elapsed = time() - overallStart
+        print('%d files (containing %d genomes and %d proteins) '
+              'processed in %.2f seconds (%.2f mins).' %
+              (fileCount, totalGenomeCount, totalProteinCount, elapsed,
+               elapsed / 60), file=sys.stderr)
+
+
+args = parser.parse_args()
+
+if args.noWarnings:
+    with warnings.catch_warnings():
+        warnings.simplefilter('ignore')
+        main(args)
+else:
+    main(args)
diff --git a/bin/noninteractive-alignment-panel.py b/bin/noninteractive-alignment-panel.py
index ccd7955..44152fc 100755
--- a/bin/noninteractive-alignment-panel.py
+++ b/bin/noninteractive-alignment-panel.py
@@ -78,7 +78,7 @@ if __name__ == '__main__':
         formatter_class=argparse.ArgumentDefaultsHelpFormatter,
         epilog=('Given BLAST or DIAMOND JSON output files, the '
                 'corresponding FASTA (or FASTQ) sequence files, and '
-                'filtering  criteria, produce an alignment panel.'))
+                'filtering criteria, produce an alignment panel.'))
 
     parser.add_argument(
         '--earlyExit', default=False, action='store_true',
diff --git a/bin/parse-genbank-flat-file.py b/bin/parse-genbank-flat-file.py
new file mode 100755
index 0000000..e2c5fea
--- /dev/null
+++ b/bin/parse-genbank-flat-file.py
@@ -0,0 +1,52 @@
+#!/usr/bin/env python
+
+from __future__ import print_function
+
+import sys
+import argparse
+from Bio import SeqIO
+
+
+parser = argparse.ArgumentParser(
+    formatter_class=argparse.ArgumentDefaultsHelpFormatter,
+    description='Parse GenBank flat files. Exit non-zero if parsing fails.')
+
+parser.add_argument(
+    '--expectedCount', type=int,
+    help=('The expected number of records in each GenBank file. If given '
+          'and this number is not present in any file, exit non-zero.'))
+
+parser.add_argument(
+    'files', nargs='+',
+    help='The files to parse.')
+
+parser.add_argument(
+    '--quiet', default=False, action='store_true',
+    help='If given, write no output.')
+
+args = parser.parse_args()
+
+expectedCount = args.expectedCount
+totalCount = 0
+
+for i in args.files:
+    count = 0
+    try:
+        records = SeqIO.parse(open(i), 'gb')
+        for record in records:
+            count += 1
+    except ValueError:
+        if not args.quiet:
+            print('Could not parse %s' % i)
+        sys.exit(1)
+    else:
+        if not args.quiet:
+            print('Read %d records from %s' % (count, i))
+        if expectedCount is not None and count != expectedCount:
+            if not args.quiet:
+                print('Expected %d records. Exiting.' % expectedCount)
+            sys.exit(1)
+        totalCount += count
+
+if not args.quiet:
+    print('Total records read: %d' % totalCount)
diff --git a/bin/proteins-to-pathogens-civ.py b/bin/proteins-to-pathogens-civ.py
new file mode 100755
index 0000000..d433672
--- /dev/null
+++ b/bin/proteins-to-pathogens-civ.py
@@ -0,0 +1,217 @@
+#!/usr/bin/env python
+
+"""
+Read protein match output produced by alignment-panel-civ.py and group it by
+pathogen (either virus or bacteria).
+
+This is currently only useful when you are matching against a subject protein
+database whose titles have a pathogen accession number in the 5th |-separated
+field, like this:
+
+civ|GENBANK|YP_008686600.1|GENBANK|NC_022580.1|glycoprotein [pathogen name]
+civ|GENBANK|YP_008686601.1|GENBANK|NC_022580.1|L protein [pathogen name]
+civ|GENBANK|YP_008686602.1|GENBANK|NC_022581.1|nucleoprotein [pathogen name]
+
+In this case, the first two matched subjects are from the same pathogen
+(accession NC_022580.1). This script will gather those matches under their
+common accession number and provides methods to print them.
+
+Files with names in this format are generated by the make-protein-database.py
+script in this directory and are used in databases such as RVDB.
+
+Note that the *only* part of the title that is used in this script is the
+pathogen accession number. The product (e.g., nucleoprotein) and pathogen name
+are just there to ease human reading of the matches. The accession number is
+looked up in a passed database (created with make-protein-database.py)
+
+The script reads *file names* from standard input, and writes to standard
+output.  Alternately, you can also provide file names on the command line.
+
+Typical usage:
+
+  $ find . -name summary-proteins |
+        proteins-to-pathogens-civ.py --html [options] > index.html
+
+Input files must contain lines in the following format:
+
+0.01 58.9 58.9 1 1 1772 title
+0.01 58.9 58.9 1 1 1772 title
+0.09 51.8 57.4 2 2 3481 title
+
+Fields must be whitespace separated. The seven fields are:
+
+    Coverage
+    Median bit score
+    Best bit score
+    Read count
+    HSP count
+    Protein length (in amino acids)
+    Title (in the above-mentioned "protein name [pathogen name]" format)
+
+"""
+
+from __future__ import print_function
+import argparse
+import sys
+
+# It's not clear that the PDF backend is the right choice here, but it
+# works (i.e., the generation of PNG images works fine).
+import matplotlib
+matplotlib.use('PDF')
+
+# These imports are here because dark.civ.proteins imports
+# matplotlib.pyplot and we need to set the matplotlib backend before the
+# import. So please don't move this import higher in this file.
+
+from dark.civ.proteins import ProteinGrouper, SqliteIndex
+
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser(
+        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
+        description="Group proteins by the pathogen they're from.")
+
+    parser.add_argument(
+        'filenames', nargs='*',
+        help=('File names to read input from. The input will typically be '
+              'generated by alignment-panel-civ.py'))
+
+    parser.add_argument(
+        '--proteinGenomeDatabase', required=True,
+        help=('The filename of an Sqlite3 database holding protein and '
+              'genome information, as built by make-protein-database.py'))
+
+    # A mutually exclusive group for either --sampleName or --sampleNameRegex
+    group = parser.add_mutually_exclusive_group()
+
+    group.add_argument(
+        '--sampleName',
+        help=('An (optional) sample name. Use when all input files are for a '
+              'single sample. Cannot be used with --sampleNameRegex.'))
+
+    group.add_argument(
+        '--sampleNameRegex',
+        help=('An (optional) regular expression that can be used to extract a '
+              'short sample name from full sample file name.  The regular '
+              'expression must have a matching group (delimited by '
+              'parentheses) that captures the part of the file name that '
+              'should be used as the sample name.'))
+
+    parser.add_argument(
+        '--pathogenPanelFilename',
+        help=('An (optional) filename to write a pathogen-sample panel PNG '
+              'image to.'))
+
+    parser.add_argument(
+        '--sampleIndexFilename',
+        help=('An (optional) filename to write a sample index file to. '
+              'Lines in the file will have an integer index, a space, and '
+              'then the sample name. Only produced if --html is used '
+              '(because the pathogen-NNN-sample-MMM.fastq are only written '
+              'in that case).'))
+
+    parser.add_argument(
+        '--html', default=False, action='store_true',
+        help='If specified, output HTML instead of plain text.')
+
+    parser.add_argument(
+        '--format', default='fasta', choices=('fasta', 'fastq'),
+        help=('Give the format of the sequence files written by '
+              'alignment-panel-civ.py.'))
+
+    parser.add_argument(
+        '--minProteinFraction', type=float, default=0.0,
+        help=('The minimum fraction of proteins in a pathogen that must be '
+              'matched by a particular sample in order for that pathogen to '
+              'be displayed for that sample.'))
+
+    parser.add_argument(
+        '--pathogenType', default='viral', choices=('bacterial', 'viral'),
+        help=('Specify the pathogen type. This option only affects the '
+              'language used in HTML output.'))
+
+    parser.add_argument(
+        '--showReadLengths', default=False, action='store_true',
+        help=('If specified, the HTML output (use --html to get this) will '
+              'contain the lengths of all reads that match proteins for a '
+              'pathogen.'))
+
+    parser.add_argument(
+        '--assetDir', default='out',
+        help=('The output directory where noninteractive-alignment-panel.py '
+              'put its plots and FASTA/FASTQ files.'))
+
+    parser.add_argument(
+        '--pathogenDataDir', default='pathogen-data',
+        help=('The directory where per-pathogen information (e.g., collected '
+              'reads across all samples) should be written.'))
+
+    parser.add_argument(
+        '--title', default='Summary of pathogens',
+        help='The title to show at the top of the output.')
+
+    parser.add_argument(
+        '--preamble',
+        help='Optional preamble text to show after the title.')
+
+    parser.add_argument(
+        '--titleRegex', default=None,
+        help='A regex that pathogen names must match.')
+
+    parser.add_argument(
+        '--negativeTitleRegex', default=None,
+        help='a regex that pathogen names must not match.')
+
+    parser.add_argument(
+        '--omitVirusLinks', default=False, action='store_true',
+        help=('If specified, the HTML output (use --html to get this) for '
+              'viruses will not contain links to ICTV and ViralZone. '
+              'This should be used when working with viruses that do not yet '
+              'have names that can be looked up.'))
+
+    args = parser.parse_args()
+
+    if not args.html:
+        if args.sampleIndexFilename:
+            print('It does not make sense to use --sampleIndexFilename '
+                  'without also using --html', file=sys.stderr)
+            sys.exit(1)
+        if args.omitVirusLinks:
+            print('It does not make sense to use --omitVirusLinks '
+                  'without also using --html', file=sys.stderr)
+            sys.exit(1)
+
+    if args.omitVirusLinks and args.pathogenType != 'viral':
+        print('The --omitVirusLinks option only makes sense with '
+              '--pathogenType viral', file=sys.stderr)
+        sys.exit(1)
+
+    grouper = ProteinGrouper(SqliteIndex(args.proteinGenomeDatabase),
+                             assetDir=args.assetDir,
+                             sampleName=args.sampleName,
+                             sampleNameRegex=args.sampleNameRegex,
+                             format_=args.format,
+                             saveReadLengths=args.showReadLengths,
+                             titleRegex=args.titleRegex,
+                             negativeTitleRegex=args.negativeTitleRegex,
+                             pathogenDataDir=args.pathogenDataDir)
+
+    if args.filenames:
+        filenames = args.filenames
+    else:
+        filenames = (line[:-1] for line in sys.stdin)
+
+    for filename in filenames:
+        with open(filename) as fp:
+            grouper.addFile(filename, fp)
+
+    if args.html:
+        print(grouper.toHTML(
+            args.pathogenPanelFilename,
+            minProteinFraction=args.minProteinFraction,
+            pathogenType=args.pathogenType,
+            title=args.title, preamble=args.preamble,
+            sampleIndexFilename=args.sampleIndexFilename,
+            omitVirusLinks=args.omitVirusLinks))
+    else:
+        print(grouper.toStr())
diff --git a/bin/proteins-to-pathogens.py b/bin/proteins-to-pathogens.py
index 3daa58e..664ac4e 100755
--- a/bin/proteins-to-pathogens.py
+++ b/bin/proteins-to-pathogens.py
@@ -74,7 +74,7 @@ if __name__ == '__main__':
         '--sampleName',
         help=('An (optional) sample name. This is only used in producing '
               'HTML output. Should be used when all input files are for a '
-              'single sample.  Cannot be used with --sampleNameRegex.'))
+              'single sample. Cannot be used with --sampleNameRegex.'))
 
     parser.add_argument(
         '--sampleNameRegex',
diff --git a/bin/subset-protein-database.py b/bin/subset-protein-database.py
new file mode 100755
index 0000000..50e9c48
--- /dev/null
+++ b/bin/subset-protein-database.py
@@ -0,0 +1,82 @@
+#!/usr/bin/env python
+
+from __future__ import print_function, division
+
+import argparse
+
+from dark.proteins import SqliteIndex
+
+
+parser = argparse.ArgumentParser(
+    formatter_class=argparse.ArgumentDefaultsHelpFormatter,
+    description=('Query an sqlite3 protein/genome database and write '
+                 'the matching sequences to standard output as FASTA.'))
+
+parser.add_argument(
+    '--databaseFile', required=True,
+    help=('The protein/genome database file. This must be an sqlite3 database '
+          'file created by make-protein-database.py'))
+
+parser.add_argument(
+    '--query',
+    help=('An SQL database query. '
+          'The query must result in rows with two columns, the sequence id '
+          'and the sequence. These can be protein or nucleotide sequences, '
+          'but it is probably better not to mix them!'))
+
+parser.add_argument(
+    '--queryFile',
+    help=('The name of a file containing an SQL database query. '
+          'The query must result in rows with two columns, the sequence id '
+          'and the sequence. These can be protein or nucleotide sequences, '
+          'but it is probably better not to mix them!'))
+
+parser.add_argument(
+    '--keepDuplicateIds', default=False, action='store_true',
+    help='If specified, do not remove duplicate ids.')
+
+parser.add_argument(
+    '--idsOnly', default=False, action='store_true',
+    help='If specified, only print matching sequence ids.')
+
+args = parser.parse_args()
+
+db = SqliteIndex(args.databaseFile)
+
+if args.queryFile:
+    query = open(args.queryFile).read()
+elif args.query:
+    query = args.query
+else:
+    raise RuntimeError('You must use either --query or --queryFile')
+
+cur = db.execute(query)
+
+removeDuplicates = not args.keepDuplicateIds
+idsOnly = args.idsOnly
+
+if removeDuplicates:
+    seen = set()
+
+while True:
+    rows = cur.fetchmany()
+    if rows:
+        for row in rows:
+            try:
+                id_, sequence, genomeName = row
+            except ValueError:
+                id_, sequence = row
+            else:
+                id_ = '%s [%s]' % (id_, genomeName)
+
+            if removeDuplicates:
+                if id_ in seen:
+                    continue
+                else:
+                    seen.add(id_)
+            if idsOnly:
+                print(id_)
+            else:
+                print('>%s\n%s' % (id_, sequence))
+    else:
+        break
diff --git a/dark/__init__.py b/dark/__init__.py
index 6a38d9a..5533064 100644
--- a/dark/__init__.py
+++ b/dark/__init__.py
@@ -7,5 +7,4 @@ if sys.version_info < (2, 7):
 # will not be found by the version() function in ../setup.py
 #
 # Remember to update ../CHANGELOG.md describing what's new in each version.
-
-__version__ = '3.0.80'
+__version__ = '3.1.0'
diff --git a/dark/btop.py b/dark/btop.py
index aa338da..3ee24ed 100644
--- a/dark/btop.py
+++ b/dark/btop.py
@@ -1,5 +1,4 @@
-# From https://samtools.github.io/hts-specs/SAMv1.pdf
-CINS, CDEL, CMATCH, CEQUAL, CDIFF = 'IDM=X'
+from dark.cigar import CINS, CDEL, CMATCH, CEQUAL, CDIFF
 
 
 def parseBtop(btopString):
diff --git a/dark/cigar.py b/dark/cigar.py
new file mode 100644
index 0000000..9606af2
--- /dev/null
+++ b/dark/cigar.py
@@ -0,0 +1,54 @@
+# From https://samtools.github.io/hts-specs/SAMv1.pdf
+CINS, CDEL, CMATCH, CEQUAL, CDIFF = 'IDM=X'
+
+
+def dna2cigar(s1, s2, concise=False):
+    """
+    Form a CIGAR string from two equal-length DNA sequences. Case is not
+    ignored, so if you want to ignore case you must convert the sequences
+    to a consistent case (upper or lower, doesn't matter) before calling.
+
+    @param s1: A C{str} of DNA 'A', 'G', 'C', 'T' characters.
+    @param s2: A C{str} of DNA 'A', 'G', 'C', 'T' characters.
+    @param concise: If C{True}, use 'M' for matches and mismatches instead
+        of the more specific 'X' and '='.
+    @raise ValueError: If the two sequences are not of equal length or if the
+        sequences are of length zero.
+    @return: A C{str} CIGAR string.
+    """
+    len1 = len(s1)
+    if len1 != len(s2):
+        raise ValueError('Sequences %r and %r of unequal length (%d != %d).' %
+                         (s1, s2, len1, len(s2)))
+
+    if len1 == 0:
+        raise ValueError('Two sequences of zero length were passed.')
+
+    if concise:
+        return '%d%s' % (len1, CMATCH)
+
+    result = []
+    length, operation = 0, None
+
+    for n1, n2 in zip(s1, s2):
+        if n1 == n2:
+            if operation == CEQUAL:
+                length += 1
+            else:
+                if length:
+                    assert operation == CDIFF
+                    result.append('%d%s' % (length, CDIFF))
+                length, operation = 1, CEQUAL
+        else:
+            if operation == CDIFF:
+                length += 1
+            else:
+                if length:
+                    assert operation == CEQUAL
+                    result.append('%d%s' % (length, CEQUAL))
+                length, operation = 1, CDIFF
+
+    # Append the final operation.
+    result.append('%d%s' % (length, operation))
+
+    return ''.join(result)
diff --git a/dark/civ/__init__.py b/dark/civ/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/dark/civ/graphics.py b/dark/civ/graphics.py
new file mode 100644
index 0000000..bc7c25a
--- /dev/null
+++ b/dark/civ/graphics.py
@@ -0,0 +1,537 @@
+import os
+from copy import deepcopy
+from stat import S_ISDIR
+from math import ceil
+from time import ctime, time
+from textwrap import fill
+from os.path import join
+
+try:
+    import matplotlib
+    if not os.environ.get('DISPLAY'):
+        # Use non-interactive Agg backend
+        matplotlib.use('Agg')
+    import matplotlib.pyplot as plt
+except ImportError:
+    import platform
+    if platform.python_implementation() == 'PyPy':
+        # PyPy doesn't have a version of matplotlib. Make fake classes and
+        # a Line2D function that raises if used. This allows us to use
+        # other 'dark' code that happens to import dark.civ.graphics but
+        # which does not use the functions that rely on matplotlib.
+        class plt(object):
+            def __getattr__(self, _):
+                raise NotImplementedError(
+                    'matplotlib is not supported under pypy')
+
+        gridspec = patches = plt
+
+        def Line2D(*args, **kwargs):
+            raise NotImplementedError('matplotlib is not supported under pypy')
+    else:
+        raise
+else:
+    from matplotlib.lines import Line2D
+    from matplotlib import gridspec
+
+from dark.dimension import dimensionalIterator
+from dark.html import NCBISequenceLinkURL
+from dark.civ.html import AlignmentPanelHTMLWriter
+from dark.intervals import ReadIntervals
+from dark.features import ProteinFeatureAdder, NucleotideFeatureAdder
+from dark.intervals import OffsetAdjuster
+from dark.score import HigherIsBetterScore
+
+
+# If we're making a plot that has a log-linear X axis, don't show
+# background light grey rectangles for any gap whose (logged) width is less
+# than SMALLEST_LOGGED_GAP_TO_DISPLAY.
+SMALLEST_LOGGED_GAP_TO_DISPLAY = 20
+
+# Y (Score) axis extra spacing above the best read match. The score of the
+# best read HSP will be multiplied by this value, the result converted to
+# an int, and then used as the upper bound on the Y axis. Adding 1% seems
+# to be a good heuristic.
+Y_AXIS_UPPER_PADDING = 1.01
+
+# The default base of the logarithm to use when logLinearXAxis is used to
+# produce an alignment graph.
+DEFAULT_LOG_LINEAR_X_AXIS_BASE = 1.1
+
+
+def report(msg):
+    print('%s: %s' % (ctime(time()), msg))
+
+
+def alignmentGraph(titlesAlignments, title, accession, addQueryLines=True,
+                   showFeatures=True, logLinearXAxis=False,
+                   logBase=DEFAULT_LOG_LINEAR_X_AXIS_BASE, rankScores=False,
+                   createFigure=True, showFigure=True,
+                   readsAx=None, imageFile=None, quiet=False, idList=False,
+                   xRange='subject'):
+    """
+    Align a set of matching reads against a BLAST or DIAMOND hit.
+
+    @param titlesAlignments: A L{dark.titles.TitlesAlignments} instance.
+    @param title: A C{str} sequence title that was matched. We plot the
+        reads that hit this title.
+    @param accession: The C{str} accession number of the matched title.
+    @param addQueryLines: if C{True}, draw query lines in full (these will then
+        be partly overdrawn by the HSP match against the subject). These are
+        the 'whiskers' that potentially protrude from each side of a query.
+    @param showFeatures: if C{True}, look online for features of the subject
+        sequence (given by hitId).
+    @param logLinearXAxis: if C{True}, convert read offsets so that empty
+        regions in the plot we're preparing will only be as wide as their
+        logged actual values.
+    @param logBase: The base of the logarithm to use if logLinearXAxis is
+        C{True}.
+    @param: rankScores: If C{True}, change the e-values and bit scores for the
+        reads for each title to be their rank (worst to best).
+    @param createFigure: If C{True}, create a figure and give it a title.
+    @param showFigure: If C{True}, show the created figure. Set this to
+        C{False} if you're creating a panel of figures or just want to save an
+        image (with C{imageFile}).
+    @param readsAx: If not None, use this as the subplot for displaying reads.
+    @param imageFile: If not None, specifies a filename to write the image to.
+    @param quiet: If C{True}, don't print progress / timing output.
+    @param idList: a dictionary. The keys is a color and the values is a list
+        of read identifiers that should be colored in the respective color.
+    @param xRange: set to either 'subject' or 'reads' to indicate the range of
+        the X axis.
+    """
+
+    startTime = time()
+
+    assert xRange in ('subject', 'reads'), (
+        'xRange must be either "subject" or "reads".')
+
+    if createFigure:
+        width = 20
+        figure = plt.figure(figsize=(width, 20))
+
+    createdReadsAx = readsAx is None
+
+    if showFeatures:
+        gs = gridspec.GridSpec(2, 1, height_ratios=[1, 1])
+        featureAx = plt.subplot(gs[0, 0])
+        readsAx = readsAx or plt.subplot(gs[1, 0])
+    else:
+        readsAx = readsAx or plt.subplot(111)
+
+    # Make a deep copy of the title alignments. We're potentially going to
+    # change the HSP scores, the X axis offsets, etc., and we don't want to
+    # interfere with the data we were passed.
+    titleAlignments = deepcopy(titlesAlignments[title])
+
+    readsAlignments = titlesAlignments.readsAlignments
+    subjectIsNucleotides = readsAlignments.params.subjectIsNucleotides
+
+    # Allow the class of titlesAlignments to adjust HSPs for plotting,
+    # if it has a method for doing so.
+    try:
+        adjuster = readsAlignments.adjustHspsForPlotting
+    except AttributeError:
+        pass
+    else:
+        adjuster(titleAlignments)
+
+    if rankScores:
+        reverse = titlesAlignments.scoreClass is not HigherIsBetterScore
+        for rank, hsp in enumerate(sorted(titleAlignments.hsps(),
+                                   reverse=reverse), start=1):
+            hsp.score.score = rank
+
+    if logLinearXAxis:
+        readIntervals = ReadIntervals(titleAlignments.subjectLength)
+        # Examine all HSPs so we can build an offset adjuster.
+        for hsp in titleAlignments.hsps():
+            readIntervals.add(hsp.readStartInSubject, hsp.readEndInSubject)
+        # Now adjust offsets in all HSPs.
+        offsetAdjuster = OffsetAdjuster(readIntervals, base=logBase)
+        for hsp in titleAlignments.hsps():
+            offsetAdjuster.adjustHSP(hsp)
+        # A function for adjusting other offsets, below.
+        adjustOffset = offsetAdjuster.adjustOffset
+    else:
+        def adjustOffset(offset):
+            return offset
+
+    # It would be more efficient to only walk through all HSPs once and
+    # compute these values all at once, but for now this is simple and clear.
+    maxY = int(ceil(titleAlignments.bestHsp().score.score))
+    minY = int(titleAlignments.worstHsp().score.score)
+    maxX = max(hsp.readEndInSubject for hsp in titleAlignments.hsps())
+    minX = min(hsp.readStartInSubject for hsp in titleAlignments.hsps())
+
+    if xRange == 'subject':
+        # We'll display a graph for the full subject range. Adjust X axis
+        # min/max to make sure we cover at least zero to the sequence length.
+        maxX = max(titleAlignments.subjectLength, maxX)
+        minX = min(0, minX)
+
+    # Swap min & max Y values, if needed, as it's possible we are dealing
+    # with LSPs but that the score adjuster made numerically greater values
+    # for those that were small.
+    if maxY < minY:
+        (maxY, minY) = (minY, maxY)
+
+    if logLinearXAxis:
+        # Adjust minX and maxX if we have gaps at the subject start or end.
+        gaps = list(readIntervals.walk())
+        if gaps:
+            # Check start of first gap:
+            intervalType, (start, stop) = gaps[0]
+            if intervalType == ReadIntervals.EMPTY:
+                adjustedStart = adjustOffset(start)
+                if adjustedStart < minX:
+                    minX = adjustedStart
+            # Check stop of last gap:
+            intervalType, (start, stop) = gaps[-1]
+            if intervalType == ReadIntervals.EMPTY:
+                adjustedStop = adjustOffset(stop)
+                if adjustedStop > maxX:
+                    maxX = adjustedStop
+
+    # We're all set up to start plotting the graph.
+
+    # Add light grey vertical rectangles to show the logarithmic gaps. Add
+    # these first so that reads will be plotted on top of them. Only draw
+    # gaps that are more than SMALLEST_LOGGED_GAP_TO_DISPLAY pixels wide as
+    # we could have millions of tiny gaps for a bacteria and drawing them
+    # all will be slow and only serves to make the entire background grey.
+    if logLinearXAxis and len(offsetAdjuster.adjustments()) < 100:
+        for (intervalType, interval) in readIntervals.walk():
+            if intervalType == ReadIntervals.EMPTY:
+                adjustedStart = adjustOffset(interval[0])
+                adjustedStop = adjustOffset(interval[1])
+                width = adjustedStop - adjustedStart
+                if width >= SMALLEST_LOGGED_GAP_TO_DISPLAY:
+                    readsAx.axvspan(adjustedStart, adjustedStop,
+                                    color='#f4f4f4')
+    else:
+        # Add horizontal lines for all the query sequences. These will be the
+        # grey 'whiskers' in the plots once we (below) draw the matched part
+        # on top of part of them.
+        if addQueryLines:
+            for hsp in titleAlignments.hsps():
+                y = hsp.score.score
+                line = Line2D([hsp.readStartInSubject, hsp.readEndInSubject],
+                              [y, y], color='#aaaaaa')
+                readsAx.add_line(line)
+
+        # Add the horizontal BLAST alignment lines.
+
+        # If an idList is given set things up to look up read colors.
+        readColor = {}
+        if idList:
+            for color, reads in idList.items():
+                for read in reads:
+                    if read in readColor:
+                        raise ValueError('Read %s is specified multiple '
+                                         'times in idList' % read)
+                    else:
+                        readColor[read] = color
+
+        # Draw the matched region.
+        for titleAlignment in titleAlignments:
+            readId = titleAlignment.read.id
+            for hsp in titleAlignment.hsps:
+                y = hsp.score.score
+                line = Line2D([hsp.subjectStart, hsp.subjectEnd], [y, y],
+                              color=readColor.get(readId, 'blue'))
+                readsAx.add_line(line)
+
+    if showFeatures:
+        if subjectIsNucleotides:
+            featureAdder = NucleotideFeatureAdder()
+        else:
+            featureAdder = ProteinFeatureAdder()
+
+        features = featureAdder.add(featureAx, title, minX, maxX,
+                                    adjustOffset)
+
+        # If there are features and there weren't too many of them, add
+        # vertical feature lines to the reads and ORF axes.
+        if features and not featureAdder.tooManyFeaturesToPlot:
+            for feature in features:
+                start = feature.start
+                end = feature.end
+                color = feature.color
+                readsAx.axvline(x=start, color=color)
+                readsAx.axvline(x=end, color='#cccccc')
+    else:
+        features = None
+
+    # We'll return some information we've gathered.
+    result = {
+        'adjustOffset': adjustOffset,
+        'features': features,
+        'minX': minX,
+        'maxX': maxX,
+        'minY': minY,
+        'maxY': maxY,
+    }
+
+    # Allow the class of titlesAlignments to add to the plot, if it has a
+    # method for doing so.
+    try:
+        adjuster = readsAlignments.adjustPlot
+    except AttributeError:
+        pass
+    else:
+        adjuster(readsAx)
+
+    # Titles, axis, etc.
+    if createFigure:
+        readCount = titleAlignments.readCount()
+        hspCount = titleAlignments.hspCount()
+        figure.suptitle(
+            '%s (%s)\nLength %d %s, %d read%s, %d HSP%s.' %
+            (
+                fill(titleAlignments.subjectTitle, 80), accession,
+                titleAlignments.subjectLength,
+                'nt' if subjectIsNucleotides else 'aa',
+                readCount, '' if readCount == 1 else 's',
+                hspCount, '' if hspCount == 1 else 's'
+            ),
+            fontsize=20)
+
+    # Add a title and y-axis label, but only if we made the reads axes.
+    if createdReadsAx:
+        readsAx.set_title('Read alignments', fontsize=20)
+        ylabel = readsAlignments.params.scoreTitle
+        if rankScores:
+            ylabel += ' rank'
+        plt.ylabel(ylabel, fontsize=17)
+
+    # Set the x-axis limits.
+    readsAx.set_xlim([minX - 1, maxX + 1])
+
+    readsAx.set_ylim([0, int(maxY * Y_AXIS_UPPER_PADDING)])
+    readsAx.grid()
+    if createFigure:
+        if showFigure:
+            plt.show()
+        if imageFile:
+            figure.savefig(imageFile)
+    stop = time()
+    if not quiet:
+        report('Graph generated in %.3f mins.' % ((stop - startTime) / 60.0))
+
+    return result
+
+
+def alignmentPanel(titlesAlignments, sortOn='maxScore', idList=False,
+                   equalizeXAxes=False, xRange='subject', logLinearXAxis=False,
+                   rankScores=False, showFeatures=True,
+                   logBase=DEFAULT_LOG_LINEAR_X_AXIS_BASE):
+    """
+    Produces a rectangular panel of graphs that each contain an alignment graph
+    against a given sequence.
+
+    @param titlesAlignments: A L{dark.titles.TitlesAlignments} instance.
+    @param sortOn: The attribute to sort subplots on. Either "maxScore",
+        "medianScore", "readCount", "length", or "title".
+    @param idList: A dictionary. Keys are colors and values are lists of read
+        ids that should be colored using that color.
+    @param equalizeXAxes: If C{True}, adjust the X axis on each alignment plot
+        to be the same.
+    @param xRange: Set to either 'subject' or 'reads' to indicate the range of
+        the X axis.
+    @param logLinearXAxis: If C{True}, convert read offsets so that empty
+        regions in the plots we're preparing will only be as wide as their
+        logged actual values.
+    @param logBase: The logarithm base to use if logLinearXAxis is C{True}.
+    @param: rankScores: If C{True}, change the scores for the reads for each
+        title to be their rank (worst to best).
+    @param showFeatures: If C{True}, look online for features of the subject
+        sequences.
+    @raise ValueError: If C{outputDir} exists but is not a directory or if
+        C{xRange} is not "subject" or "reads".
+    """
+
+    if xRange not in ('subject', 'reads'):
+        raise ValueError('xRange must be either "subject" or "reads".')
+
+    start = time()
+    titles = titlesAlignments.sortTitles(sortOn)
+    cols = 5
+    rows = int(len(titles) / cols) + (0 if len(titles) % cols == 0 else 1)
+    figure, ax = plt.subplots(rows, cols, squeeze=False)
+    allGraphInfo = {}
+    coords = dimensionalIterator((rows, cols))
+
+    report('Plotting %d titles in %dx%d grid, sorted on %s' %
+           (len(titles), rows, cols, sortOn))
+
+    for i, title in enumerate(titles):
+        titleAlignments = titlesAlignments[title]
+        row, col = next(coords)
+        report('%d: %s %s' % (i, title, NCBISequenceLinkURL(title, '')))
+
+        # Add a small plot to the alignment panel.
+        graphInfo = alignmentGraph(
+            titlesAlignments, title, addQueryLines=True,
+            showFeatures=showFeatures, rankScores=rankScores,
+            logLinearXAxis=logLinearXAxis, logBase=logBase,
+            createFigure=False, showFigure=False,
+            readsAx=ax[row][col], quiet=True, idList=idList, xRange=xRange)
+
+        allGraphInfo[title] = graphInfo
+        readCount = titleAlignments.readCount()
+        hspCount = titleAlignments.hspCount()
+
+        # Make a short title for the small panel blue plot, ignoring any
+        # leading NCBI gi / accession numbers.
+        if title.startswith('gi|') and title.find(' ') > -1:
+            shortTitle = title.split(' ', 1)[1][:40]
+        else:
+            shortTitle = title[:40]
+
+        plotTitle = ('%d: %s\nLength %d, %d read%s, %d HSP%s.' % (
+            i, shortTitle, titleAlignments.subjectLength,
+            readCount, '' if readCount == 1 else 's',
+            hspCount, '' if hspCount == 1 else 's'))
+
+        if hspCount:
+            if rankScores:
+                plotTitle += '\nY axis is ranked score'
+            else:
+                plotTitle += '\nmax %.2f, median %.2f' % (
+                    titleAlignments.bestHsp().score.score,
+                    titleAlignments.medianScore())
+
+        ax[row][col].set_title(plotTitle, fontsize=10)
+
+    maxX = max(graphInfo['maxX'] for graphInfo in allGraphInfo.values())
+    minX = min(graphInfo['minX'] for graphInfo in allGraphInfo.values())
+    maxY = max(graphInfo['maxY'] for graphInfo in allGraphInfo.values())
+    minY = min(graphInfo['minY'] for graphInfo in allGraphInfo.values())
+
+    # Post-process graphs to adjust axes, etc.
+
+    coords = dimensionalIterator((rows, cols))
+    for title in titles:
+        titleAlignments = titlesAlignments[title]
+        row, col = next(coords)
+        a = ax[row][col]
+        a.set_ylim([0, int(maxY * Y_AXIS_UPPER_PADDING)])
+        if equalizeXAxes:
+            a.set_xlim([minX, maxX])
+        a.set_yticks([])
+        a.set_xticks([])
+
+        if xRange == 'subject' and minX < 0:
+            # Add a vertical line at x=0 so we can see the 'whiskers' of
+            # reads that extend to the left of the sequence we're aligning
+            # against.
+            a.axvline(x=0, color='#cccccc')
+
+        # Add a line on the right of each sub-plot so we can see where the
+        # sequence ends (as all panel graphs have the same width and we
+        # otherwise couldn't tell).
+        sequenceLen = titleAlignments.subjectLength
+        if logLinearXAxis:
+            sequenceLen = allGraphInfo[title]['adjustOffset'](sequenceLen)
+        a.axvline(x=sequenceLen, color='#cccccc')
+
+    # Hide the final panel graphs (if any) that have no content. We do this
+    # because the panel is a rectangular grid and some of the plots at the
+    # end of the last row may be unused.
+    for row, col in coords:
+        ax[row][col].axis('off')
+
+    # plt.subplots_adjust(left=0.01, bottom=0.01, right=0.99, top=0.93,
+    # wspace=0.1, hspace=None)
+    plt.subplots_adjust(hspace=0.4)
+    figure.suptitle('X: %d to %d, Y (%s): %d to %d' %
+                    (minX, maxX,
+                     titlesAlignments.readsAlignments.params.scoreTitle,
+                     int(minY), int(maxY)), fontsize=20)
+    figure.set_size_inches(5 * cols, 3 * rows, forward=True)
+    figure.show()
+    stop = time()
+    report('Alignment panel generated in %.3f mins.' % ((stop - start) / 60.0))
+
+
+def alignmentPanelHTML(titlesAlignments, proteinGenomeDB, sortOn='maxScore',
+                       outputDir=None, idList=False, equalizeXAxes=False,
+                       xRange='subject', logLinearXAxis=False,
+                       logBase=DEFAULT_LOG_LINEAR_X_AXIS_BASE,
+                       rankScores=False, showFeatures=True,
+                       subjectType='protein'):
+    """
+    Produces an HTML index file in C{outputDir} and a collection of alignment
+    graphs and FASTA files to summarize the information in C{titlesAlignments}.
+
+    @param titlesAlignments: A L{dark.titles.TitlesAlignments} instance.
+    @param proteinGenomeDB: A L{dark.civ.proteins.SqliteIndex} instance.
+    @param sortOn: The attribute to sort subplots on. Either "maxScore",
+        "medianScore", "readCount", "length", or "title".
+    @param outputDir: Specifies a C{str} directory to write the HTML to. If
+        the directory does not exist it will be created.
+    @param idList: A dictionary. Keys are colors and values are lists of read
+        ids that should be colored using that color.
+    @param equalizeXAxes: If C{True}, adjust the X axis on each alignment plot
+        to be the same.
+    @param xRange: Set to either 'subject' or 'reads' to indicate the range of
+        the X axis.
+    @param logLinearXAxis: If C{True}, convert read offsets so that empty
+        regions in the plots we're preparing will only be as wide as their
+        logged actual values.
+    @param logBase: The logarithm base to use if logLinearXAxis is C{True}.
+    @param: rankScores: If C{True}, change the scores for the reads for each
+        title to be their rank (worst to best).
+    @param showFeatures: If C{True}, look online for features of the subject
+        sequences.
+    @param subjectType: A C{str} indicating whether the matched subjects are
+        'protein' or 'genome'. This is used to determine what accession number
+        to extract from the matched title.
+    @raise TypeError: If C{outputDir} is C{None}.
+    @raise ValueError: If C{outputDir} is None or exists but is not a
+        directory or if C{xRange} is not "subject" or "reads".
+    """
+    if subjectType not in ('protein', 'genome'):
+        raise ValueError('subjectType must be either "protein" or "genome".')
+
+    if xRange not in ('subject', 'reads'):
+        raise ValueError('xRange must be either "subject" or "reads".')
+
+    if equalizeXAxes:
+        raise NotImplementedError('This feature is not yet implemented.')
+
+    titles = titlesAlignments.sortTitles(sortOn)
+
+    if os.access(outputDir, os.F_OK):
+        # outputDir exists. Check it's a directory.
+        if not S_ISDIR(os.stat(outputDir).st_mode):
+            raise ValueError('%r is not a directory.' % outputDir)
+    else:
+        if outputDir is None:
+            raise ValueError('The outputDir needs to be specified.')
+        else:
+            os.mkdir(outputDir)
+
+    htmlWriter = AlignmentPanelHTMLWriter(outputDir, titlesAlignments)
+
+    getAccession = (proteinGenomeDB.proteinAccession
+                    if subjectType == 'protein'
+                    else proteinGenomeDB.genomeAccession)
+
+    for title in titles:
+        accession = getAccession(title)
+        imageBasename = accession + '.png'
+        imageFile = join(outputDir, imageBasename)
+        graphInfo = alignmentGraph(
+            titlesAlignments, title, accession, addQueryLines=True,
+            showFeatures=showFeatures, rankScores=rankScores,
+            logLinearXAxis=logLinearXAxis, logBase=logBase,
+            showFigure=False, imageFile=imageFile,
+            quiet=True, idList=idList, xRange=xRange)
+
+        # Close the image plot to make sure memory is flushed.
+        plt.close()
+        htmlWriter.addImage(imageBasename, accession, title, graphInfo)
+
+    htmlWriter.close()
diff --git a/dark/civ/html.py b/dark/civ/html.py
new file mode 100644
index 0000000..4b43804
--- /dev/null
+++ b/dark/civ/html.py
@@ -0,0 +1,198 @@
+from __future__ import print_function
+
+from os.path import join
+
+from dark.fastq import FastqReads
+from dark.html import NCBISequenceLinkURL
+
+
+class AlignmentPanelHTMLWriter(object):
+    """
+    Produces HTML details of a rectangular panel of graphs that each
+    contain an alignment graph against a given sequence. This is
+    supplementary output info for the AlignmentPanel class in graphics.py.
+
+    @param outputDir: The C{str} directory to write files into.
+    @param titlesAlignments: A L{dark.titles.TitlesAlignments} instance.
+    """
+    def __init__(self, outputDir, titlesAlignments):
+        self._outputDir = outputDir
+        self._titlesAlignments = titlesAlignments
+        self._images = []
+
+    def addImage(self, imageBasename, accession, title, graphInfo):
+        self._images.append({
+            'accession': accession,
+            'graphInfo': graphInfo,
+            'imageBasename': imageBasename,
+            'title': title
+        })
+
+    def close(self):
+        with open(join(self._outputDir, 'index.html'), 'w') as fp:
+            self._writeHeader(fp)
+            self._writeBody(fp)
+            self._writeFooter(fp)
+        with open(join(self._outputDir, 'style.css'), 'w') as fp:
+            self._writeCSS(fp)
+
+    def _writeHeader(self, fp):
+        fp.write("""\
+<html>
+  <head>
+    <title>Read alignments for %d matched subjects</title>
+    <link rel="stylesheet" type="text/css" href="style.css">
+  </head>
+  <body>
+    <div id="content">
+        """ % len(self._images))
+
+    def _writeBody(self, fp):
+        fp.write('<h1>Read alignments for %d matched subjects</h1>\n' %
+                 len(self._images))
+
+        # Write out an alignment panel as a table.
+        cols = 6
+        fp.write('<table><tbody>\n')
+
+        for i, image in enumerate(self._images):
+            title = image['title']
+            accession = image['accession']
+            if i % cols == 0:
+                fp.write('<tr>\n')
+
+            fp.write(
+                '<td><a id="small_%s"></a><a href="#big_%s"><img src="%s" '
+                'class="thumbnail"/></a></td>\n' %
+                (accession, accession, image['imageBasename']))
+
+            if i % cols == cols - 1:
+                fp.write('</tr>')
+
+        # Add empty cells to the final table row, and close the row, if
+        # necessary.
+        if i % cols < cols - 1:
+            while i % cols < cols - 1:
+                fp.write('<td>&nbsp;</td>\n')
+                i += 1
+            fp.write('</tr>\n')
+
+        fp.write('</tbody></table>\n')
+
+        # Write out the full images with additional detail.
+        for i, image in enumerate(self._images):
+            title = image['title']
+            accession = image['accession']
+            titleAlignments = self._titlesAlignments[title]
+            graphInfo = image['graphInfo']
+            readFormat = self._writeReads(image)
+            fp.write("""
+      <a id="big_%s"></a>
+      <h3>%d: %s</h3>
+      <p>
+            Length: %d.
+            Read count: %d.
+            HSP count: %d.
+            <a href="%s.%s">%s</a>.
+            <a href="#small_%s">Top panel.</a>
+"""
+                     % (accession,
+                        i, title,
+                        titleAlignments.subjectLength,
+                        titleAlignments.readCount(),
+                        titleAlignments.hspCount(),
+                        accession, readFormat, readFormat,
+                        accession))
+
+            url = NCBISequenceLinkURL(title)
+            if url:
+                fp.write('<a href="%s" target="_blank">NCBI</a>.' % url)
+
+            # Write out feature information.
+            if graphInfo['features'] is None:
+                # Feature lookup was False (or we were offline).
+                pass
+            elif len(graphInfo['features']) == 0:
+                fp.write('There were no features.')
+            else:
+                fp.write('<a href="%s">Features</a>' %
+                         self._writeFeatures(i, image))
+
+            # Write out the titles that this title invalidated due to its
+            # read set.
+            readSetFilter = self._titlesAlignments.readSetFilter
+            if readSetFilter:
+                invalidated = readSetFilter.invalidates(title)
+                if invalidated:
+                    nInvalidated = len(invalidated)
+                    fp.write('<br/>This title invalidated %d other%s due to '
+                             'its read set:<ul>'
+                             % (nInvalidated,
+                                '' if nInvalidated == 1 else 's'))
+                    for title in invalidated:
+                        fp.write('<li>%s</li>' % title)
+                    fp.write('</ul>')
+
+            fp.write('</p><img src="%s" class="full-size"/>' %
+                     image['imageBasename'])
+
+    def _writeFooter(self, fp):
+        fp.write("""\
+    </div>
+  </body>
+</html>
+""")
+
+    def _writeCSS(self, fp):
+        fp.write("""\
+#content {
+  width: 95%;
+  margin: auto;
+}
+img.thumbnail {
+  height: 300px;
+}
+img.full-size {
+  height: 900px;
+}
+""")
+
+    def _writeReads(self, image):
+        """
+        Write a FASTA or FASTQ file containing the set of reads that hit a
+        sequence.
+
+        @param image: A member of self._images.
+        @return: A C{str}, either 'fasta' or 'fastq' indicating the format
+            of the reads in C{self._titlesAlignments}.
+        """
+        if isinstance(self._titlesAlignments.readsAlignments.reads,
+                      FastqReads):
+            format_ = 'fastq'
+        else:
+            format_ = 'fasta'
+        filename = join(self._outputDir,
+                        '%s.%s' % (image['accession'], format_))
+        titleAlignments = self._titlesAlignments[image['title']]
+        with open(filename, 'w') as fp:
+            for titleAlignment in titleAlignments:
+                fp.write(titleAlignment.read.toString(format_))
+        return format_
+
+    def _writeFeatures(self, image):
+        """
+        Write a text file containing the features as a table.
+
+        @param image: A member of self._images.
+        @return: The C{str} features file name - just the base name, not
+            including the path to the file.
+        """
+        basename = image['accession'] + '-features.txt'
+        filename = join(self._outputDir, basename)
+        featureList = image['graphInfo']['features']
+        # Note that the following (deliberately) creates an empty features
+        # file if there were no features.
+        with open(filename, 'w') as fp:
+            for feature in featureList:
+                fp.write('%s\n\n' % feature.feature)
+        return basename
diff --git a/dark/civ/proteins.py b/dark/civ/proteins.py
new file mode 100644
index 0000000..4426f54
--- /dev/null
+++ b/dark/civ/proteins.py
@@ -0,0 +1,1779 @@
+from __future__ import print_function
+
+import os
+import re
+import sqlite3
+import sys
+import numpy as np
+from Bio import SeqIO
+from cachetools import LRUCache, cachedmethod
+from collections import defaultdict
+from json import load
+from operator import attrgetter, itemgetter
+from os.path import dirname, exists, join
+from six import string_types
+from six.moves.urllib.parse import quote
+from textwrap import fill
+from warnings import warn
+
+from dark.dimension import dimensionalIterator
+from dark.errors import DatabaseDuplicationError
+from dark.fasta import FastaReads
+from dark.fastq import FastqReads
+from dark.filter import TitleFilter
+from dark.genbank import GenomeRanges
+from dark.html import NCBISequenceLinkURL
+from dark.reads import Reads
+from dark.taxonomy import isRNAVirus, formatLineage
+
+
+class PathogenSampleFiles(object):
+    """
+    Maintain a cache of FASTA/FASTQ file names for the samples that contain a
+    given pathogen, create de-duplicated (by read id) FASTA/FASTQ files
+    for each pathogen/sample pair, provide functions to write out index files
+    of samples numbers (which are generated here in C{self.add}),
+    and provide a filename lookup function for pathogen/sample combinations
+    or just pathogen accessions by themselves.
+
+    @param proteinGrouper: An instance of C{ProteinGrouper}.
+    @param format_: A C{str}, either 'fasta' or 'fastq' indicating the format
+        of the files containing the reads matching proteins.
+    @raise ValueError: If C{format_} is unknown.
+    """
+    def __init__(self, proteinGrouper, format_='fasta'):
+        self._proteinGrouper = proteinGrouper
+        if format_ in ('fasta', 'fastq'):
+            self._format = format_
+            self._readsClass = FastaReads if format_ == 'fasta' else FastqReads
+        else:
+            raise ValueError("format_ must be either 'fasta' or 'fastq'.")
+        self._pathogens = {}
+        self._samples = {}
+        self._readsFilenames = {}
+
+    def add(self, genomeAccession, sampleName):
+        """
+        Add a (pathogen accession number, sample name) combination and get its
+        FASTA/FASTQ file name and unique read count. Write the FASTA/FASTQ file
+        if it does not already exist. Save the unique read count into
+        C{self._proteinGrouper}.
+
+        @param genomeAccession: A C{str} pathogen accession number.
+        @param sampleName: A C{str} sample name.
+        @return: A C{str} giving the FASTA/FASTQ file name holding all the
+            reads (without duplicates, by id) from the sample that matched the
+            proteins in the given pathogen.
+        """
+        sampleIndex = self._samples.setdefault(sampleName, len(self._samples))
+
+        try:
+            return self._readsFilenames[(genomeAccession, sampleIndex)]
+        except KeyError:
+            reads = Reads()
+            for proteinMatch in self._proteinGrouper.genomeAccessions[
+                    genomeAccession][sampleName]['proteins'].values():
+                for read in self._readsClass(proteinMatch['readsFilename']):
+                    reads.add(read)
+            saveFilename = join(
+                proteinMatch['outDir'],
+                'pathogen-%s-sample-%d.%s' % (genomeAccession, sampleIndex,
+                                              self._format))
+            reads.filter(removeDuplicatesById=True)
+            nReads = reads.save(saveFilename, format_=self._format)
+            # Save the unique read count into self._proteinGrouper
+            self._proteinGrouper.genomeAccessions[
+                genomeAccession][sampleName]['uniqueReadCount'] = nReads
+            self._readsFilenames[
+                (genomeAccession, sampleIndex)] = saveFilename
+            return saveFilename
+
+    def lookup(self, genomeAccession, sampleName):
+        """
+        Look up a pathogen accession number, sample name combination and get
+        its FASTA/FASTQ file name.
+
+        This method should be used instead of C{add} in situations where
+        you want an exception to be raised if a pathogen/sample combination has
+        not already been passed to C{add}.
+
+        @param genomeAccession: A C{str} pathogen accession number.
+        @param sampleName: A C{str} sample name.
+        @raise KeyError: If the pathogen accession number or sample name have
+            not been seen, either individually or in combination.
+        @return: A C{str} filename retrieved from self._readsFilenames
+        """
+        return self._readsFilenames[
+            (genomeAccession, self._samples[sampleName])]
+
+    def writeSampleIndex(self, fp):
+        """
+        Write a file of sample indices and names, sorted by index.
+
+        @param fp: A file-like object, opened for writing.
+        """
+        print('\n'.join(
+            '%d %s' % (index, name) for (index, name) in
+            sorted((index, name) for (name, index) in self._samples.items())
+        ), file=fp)
+
+
+class ProteinGrouper(object):
+    """
+    Group matched proteins by the pathogen they come from.
+
+    @param proteinGenomeDatabase: A connection to an Sqlite3 database
+        holding protein and genome information, as built by
+        C{make-protein-database.py}.
+    @param assetDir: The C{str} directory name where
+        C{noninteractive-alignment-panel.py} put its HTML, blue plot and
+        alignment panel images, and FASTA or FASTQ files. This must be relative
+        to the filenames that will later be passed to C{addFile}.
+    @param sampleName: A C{str} sample name. This takes precedence over
+        C{sampleNameRegex} (the two cannot be used together, obviously).
+    @param sampleNameRegex: A C{str} regular expression that can be used to
+        extract a short sample name from full file names subsequently passed
+        to C{self.addFile}. The regular expression must have a matching group
+        (delimited by parentheses) to capture the part of the file name that
+        should be used as the sample name.
+    @param format_: A C{str}, either 'fasta' or 'fastq' indicating the format
+        of the files containing the reads matching proteins.
+    @param saveReadLengths: If C{True}, save the lengths of all reads matching
+        proteins.
+    @param titleRegex: A regex that pathogen names must match.
+        Note that this matching is done on the final part of the protein title
+        in square brackets, according to the convention used by the NCBI viral
+        refseq database and RVDB.
+    @param negativeTitleRegex: A regex that pathogen names must not match.
+        Note that this matching is done on the final part of the protein title
+        in square brackets, according to the convention used by the NCBI viral
+        refseq database and RVDB.
+    @param pathogenDataDir: The C{str} directory where per-pathogen information
+        (e.g., collected reads across all samples) should be written. Will be
+        created (in C{self.toHTML}) if it doesn't exist.
+    @raise ValueError: If C{format_} is unknown.
+    """
+
+    VIRALZONE = 'https://viralzone.expasy.org/search?query='
+    ICTV = 'https://talk.ictvonline.org/search-124283882/?q='
+    READCOUNT_MARKER = '*READ-COUNT*'
+    READ_AND_HSP_COUNT_STR_SEP = '/'
+
+    def __init__(self, proteinGenomeDatabase, assetDir='out', sampleName=None,
+                 sampleNameRegex=None, format_='fasta',
+                 saveReadLengths=False, titleRegex=None,
+                 negativeTitleRegex=None, pathogenDataDir='pathogen-data'):
+        self._db = proteinGenomeDatabase
+        self._assetDir = assetDir
+        self._sampleName = sampleName
+        self._sampleNameRegex = (re.compile(sampleNameRegex) if sampleNameRegex
+                                 else None)
+        if format_ in ('fasta', 'fastq'):
+            self._format = format_
+        else:
+            raise ValueError("format_ must be either 'fasta' or 'fastq'.")
+        self._saveReadLengths = saveReadLengths
+
+        if titleRegex or negativeTitleRegex:
+            self.titleFilter = TitleFilter(
+                positiveRegex=titleRegex, negativeRegex=negativeTitleRegex)
+        else:
+            self.titleFilter = None
+
+        self._pathogenDataDir = pathogenDataDir
+
+        # genomeAccessions will be a dict of dicts of dicts. The first
+        # two keys will be a pathogen accession and a sample name. The
+        # final dict will contain 'proteins' (a list of dicts) and
+        # 'uniqueReadCount' (an int).
+        self.genomeAccessions = defaultdict(dict)
+        # sampleNames is keyed by sample name and will have values that hold
+        # the sample's alignment panel index.html file.
+        self.sampleNames = {}
+        self.pathogenSampleFiles = PathogenSampleFiles(self, format_=format_)
+
+    def _title(self, pathogenType):
+        """
+        Create a title summarizing the pathogens and samples.
+
+        @param pathogenType: A C{str}, either 'viral' or 'bacterial'.
+        @return: A C{str} title.
+        """
+
+        assert pathogenType in ('viral', 'bacterial')
+
+        nPathogens = len(self.genomeAccessions)
+        nSamples = len(self.sampleNames)
+
+        if pathogenType == 'bacterial':
+            what = 'bacterium' if nPathogens == 1 else 'bacteria'
+        else:
+            what = 'virus%s' % ('' if nPathogens == 1 else 'es')
+
+        return (
+            'Overall, proteins from %d %s were found in %d sample%s.' %
+            (nPathogens, what, nSamples, '' if nSamples == 1 else 's'))
+
+    def addFile(self, filename, fp):
+        """
+        Read and record protein information for a sample.
+
+        @param filename: A C{str} file name.
+        @param fp: An open file pointer to read the file's data from.
+        @raise ValueError: If information for a pathogen/protein/sample
+            combination is given more than once.
+        """
+        if self._sampleName:
+            sampleName = self._sampleName
+        elif self._sampleNameRegex:
+            match = self._sampleNameRegex.search(filename)
+            if match:
+                sampleName = match.group(1)
+            else:
+                sampleName = filename
+        else:
+            sampleName = filename
+
+        outDir = join(dirname(filename), self._assetDir)
+
+        self.sampleNames[sampleName] = join(outDir, 'index.html')
+
+        for index, proteinLine in enumerate(fp):
+            (coverage, medianScore, bestScore, readCount, hspCount,
+             proteinLength, longName) = proteinLine.split(None, 6)
+
+            proteinInfo = self._db.findProtein(longName)
+            proteinName = (proteinInfo['product'] or proteinInfo['gene'] or
+                           'unknown')
+            proteinAccession = proteinInfo['accession']
+
+            genomeInfo = self._db.findGenome(longName)
+            genomeName = genomeInfo['name']
+            genomeAccession = genomeInfo['accession']
+
+            # Ignore genomes with names we don't want.
+            if (self.titleFilter and self.titleFilter.accept(
+                    genomeName) == TitleFilter.REJECT):
+                continue
+
+            if sampleName not in self.genomeAccessions[genomeAccession]:
+                self.genomeAccessions[genomeAccession][sampleName] = {
+                    'proteins': {},
+                    'uniqueReadCount': None,
+                }
+
+            proteins = self.genomeAccessions[
+                genomeAccession][sampleName]['proteins']
+
+            # We should only receive one line of information for a given
+            # genome/sample/protein combination.
+            if proteinAccession in proteins:
+                raise ValueError(
+                    'Protein %r already seen for genome %r (%s) sample %r.' %
+                    (proteinAccession, genomeName, genomeAccession,
+                     sampleName))
+
+            readsFilename = join(outDir,
+                                 '%s.%s' % (proteinAccession, self._format))
+
+            if longName.startswith(SqliteIndexWriter.SEQUENCE_ID_PREFIX +
+                                   SqliteIndexWriter.SEQUENCE_ID_SEPARATOR):
+                proteinURL = NCBISequenceLinkURL(longName, field=2)
+                genomeURL = NCBISequenceLinkURL(longName, field=4)
+            else:
+                proteinURL = genomeURL = None
+
+            proteinInfo = proteins[proteinAccession] = {
+                'accession': proteinAccession,
+                'bestScore': float(bestScore),
+                'bluePlotFilename': join(outDir, '%s.png' % proteinAccession),
+                'coverage': float(coverage),
+                'readsFilename': readsFilename,
+                'hspCount': int(hspCount),
+                'index': index,
+                'medianScore': float(medianScore),
+                'outDir': outDir,
+                'proteinLength': int(proteinLength),
+                'proteinName': proteinName,
+                'proteinURL': proteinURL,
+                'genomeURL': genomeURL,
+                'readCount': int(readCount),
+            }
+
+            if proteinInfo['readCount'] == proteinInfo['hspCount']:
+                proteinInfo['readAndHspCountStr'] = readCount
+            else:
+                proteinInfo['readAndHspCountStr'] = '%s%s%s' % (
+                    readCount, self.READ_AND_HSP_COUNT_STR_SEP, hspCount)
+
+            if self._saveReadLengths:
+                readsClass = (FastaReads if self._format == 'fasta'
+                              else FastqReads)
+                proteins[proteinName]['readLengths'] = tuple(
+                    len(read) for read in readsClass(readsFilename))
+
+    def _computeUniqueReadCounts(self):
+        """
+        Add all pathogen / sample combinations to self.pathogenSampleFiles.
+
+        This will make all de-duplicated (by id) FASTA/FASTQ files and store
+        the number of de-duplicated reads into C{self.genomeAccessions}.
+        """
+        for genomeAccession, samples in self.genomeAccessions.items():
+            for sampleName in samples:
+                self.pathogenSampleFiles.add(genomeAccession, sampleName)
+
+    def toStr(self, title='Summary of pathogens', preamble=None):
+        """
+        Produce a string representation of the pathogen summary.
+
+        @param title: The C{str} title for the output.
+        @param preamble: The C{str} descriptive preamble, or C{None} if no
+            preamble is needed.
+        @return: A C{str} suitable for printing.
+        """
+        # Note that the string representation contains much less
+        # information than the HTML summary. E.g., it does not contain the
+        # unique (de-duplicated, by id) read count, since that is only computed
+        # when we are making combined FASTA files of reads matching a
+        # pathogen.
+        readCountGetter = itemgetter('readCount')
+        result = []
+        append = result.append
+
+        result.extend((title, ''))
+        if preamble:
+            result.extend((preamble, ''))
+        result.extend((self._title(), ''))
+
+        for pathogenName in sorted(self.pathogenNames):
+            samples = self.pathogenNames[pathogenName]
+            sampleCount = len(samples)
+            append('%s (in %d sample%s)' %
+                   (pathogenName,
+                    sampleCount, '' if sampleCount == 1 else 's'))
+            for sampleName in sorted(samples):
+                proteins = samples[sampleName]['proteins']
+                proteinCount = len(proteins)
+                totalReads = sum(readCountGetter(p) for p in proteins.values())
+                append('  %s (%d protein%s, %d read%s)' %
+                       (sampleName,
+                        proteinCount, '' if proteinCount == 1 else 's',
+                        totalReads, '' if totalReads == 1 else 's'))
+                for proteinName in sorted(proteins):
+                    append(
+                        '    %(coverage).2f\t%(medianScore).2f\t'
+                        '%(bestScore).2f\t%(readAndHspCountStr)11s\t'
+                        '%(proteinName)s'
+                        % proteins[proteinName])
+            append('')
+
+        return '\n'.join(result)
+
+    def toHTML(self, pathogenPanelFilename=None, minProteinFraction=0.0,
+               pathogenType='viral', title='Summary of pathogens',
+               preamble=None, sampleIndexFilename=None,
+               omitVirusLinks=False):
+        """
+        Produce an HTML string representation of the pathogen summary.
+
+        @param pathogenPanelFilename: If not C{None}, a C{str} filename to
+            write a pathogen panel PNG image to.
+        @param minProteinFraction: The C{float} minimum fraction of proteins
+            in a pathogen that must be matched by a sample in order for that
+            pathogen to be displayed for that sample.
+        @param pathogenType: A C{str} giving the type of the pathogen involved,
+            either 'bacterial' or 'viral'.
+        @param title: The C{str} title for the HTML page.
+        @param preamble: The C{str} descriptive preamble for the HTML page, or
+            C{None} if no preamble is needed.
+        @param sampleIndexFilename: A C{str} filename to write a sample index
+            file to. Lines in the file will have an integer index, a space, and
+            then the sample name.
+        @param omitVirusLinks: If C{True}, links to ICTV and ViralZone will be
+            omitted in output.
+        @return: An HTML C{str} suitable for printing.
+        """
+        if pathogenType not in ('bacterial', 'viral'):
+            raise ValueError(
+                "Unrecognized pathogenType argument: %r. Value must be either "
+                "'bacterial' or 'viral'." % pathogenType)
+
+        if not exists(self._pathogenDataDir):
+            os.mkdir(self._pathogenDataDir)
+
+        self._computeUniqueReadCounts()
+
+        if pathogenPanelFilename:
+            self.pathogenPanel(pathogenPanelFilename)
+
+        if sampleIndexFilename:
+            with open(sampleIndexFilename, 'w') as fp:
+                self.pathogenSampleFiles.writeSampleIndex(fp)
+
+        # Figure out if we have to delete some pathogens because the
+        # fraction of their proteins that we have matches for is too low.
+        if minProteinFraction > 0.0:
+            toDelete = defaultdict(list)
+            for genomeAccession in self.genomeAccessions:
+                genomeInfo = self._db.findGenome(genomeAccession)
+                proteinCount = genomeInfo['proteinCount']
+                assert proteinCount > 0
+                for s in self.genomeAccessions[genomeAccession]:
+                    sampleProteinFraction = (
+                        len(self.genomeAccessions[
+                            genomeAccession][s]['proteins']) /
+                        proteinCount)
+                    if sampleProteinFraction < minProteinFraction:
+                        toDelete[genomeAccession].append(s)
+
+            for genomeAccession, samples in toDelete.items():
+                for sample in samples:
+                    del self.genomeAccessions[genomeAccession][sample]
+
+        genomeAccessions = sorted(
+            genomeAccession for genomeAccession in self.genomeAccessions
+            if len(self.genomeAccessions[genomeAccession]) > 0)
+        nPathogenNames = len(genomeAccessions)
+        sampleNames = sorted(self.sampleNames)
+
+        result = [
+            '<html>',
+            '<head>',
+            '<title>',
+            title,
+            '</title>',
+            '<meta charset="UTF-8">',
+            '</head>',
+            '<body>',
+            '<style>',
+            '''\
+            body {
+                margin-left: 2%;
+                margin-right: 2%;
+            }
+            hr {
+                display: block;
+                margin-top: 0.5em;
+                margin-bottom: 0.5em;
+                margin-left: auto;
+                margin-right: auto;
+                border-style: inset;
+                border-width: 1px;
+            }
+            p.pathogen {
+                margin-top: 10px;
+                margin-bottom: 3px;
+            }
+            p.sample {
+                margin-top: 10px;
+                margin-bottom: 3px;
+            }
+            .sample {
+                margin-top: 5px;
+                margin-bottom: 2px;
+            }
+            ul {
+                margin-bottom: 2px;
+            }
+            .indented {
+                margin-left: 2em;
+            }
+            .sample-name {
+                font-size: 125%;
+                font-weight: bold;
+            }
+            .pathogen-name {
+                font-size: 125%;
+                font-weight: bold;
+            }
+            .index-name {
+                font-weight: bold;
+            }
+            .index {
+                font-size: small;
+            }
+            .protein-name {
+            }
+            .stats {
+                font-family: "Courier New", Courier, monospace;
+                white-space: pre;
+            }
+            .protein-list {
+                margin-top: 2px;
+            }''',
+            '</style>',
+            '</head>',
+            '<body>',
+        ]
+
+        proteinFieldsDescription = [
+            '<p>',
+            'In all bullet point protein lists below, there are the following '
+            'fields:',
+            '<ol>',
+            '<li>Coverage fraction.</li>',
+            '<li>Median bit score.</li>',
+            '<li>Best bit score.</li>',
+            '<li>Read count (if the HSP count differs, read and HSP ',
+            ('counts are both given, separated by "%s").</li>' %
+             self.READ_AND_HSP_COUNT_STR_SEP),
+            '<li>Protein length (in amino acids).</li>',
+        ]
+
+        if self._saveReadLengths:
+            proteinFieldsDescription.append(
+                '<li>All read lengths (in parentheses).</li>')
+
+        proteinFieldsDescription.extend([
+            '<li>Protein name.</li>',
+            '</ol>',
+            '</p>',
+        ])
+
+        append = result.append
+
+        append('<h1>%s</h1>' % title)
+        if preamble:
+            append('<p>%s</p>' % preamble)
+        append('<p>')
+        append(self._title(pathogenType))
+
+        if minProteinFraction > 0.0:
+            percent = minProteinFraction * 100.0
+            if nPathogenNames < len(self.genomeAccessions):
+                if nPathogenNames == 1:
+                    append('Pathogen protein fraction filtering has been '
+                           'applied, so information on only 1 pathogen is '
+                           'displayed. This is the only pathogen for which at '
+                           'least one sample matches at least %.2f%% of the '
+                           'pathogen proteins.' % percent)
+                else:
+                    append('Pathogen protein fraction filtering has been '
+                           'applied, so information on only %d pathogens is '
+                           'displayed. These are the only pathogens for which '
+                           'at least one sample matches at least %.2f%% of '
+                           'the pathogen proteins.' % (nPathogenNames,
+                                                       percent))
+            else:
+                append('Pathogen protein fraction filtering was applied, '
+                       'but all pathogens have at least %.2f%% of their '
+                       'proteins matched by at least one sample.' % percent)
+
+        append('</p>')
+
+        if pathogenPanelFilename:
+            append('<p>')
+            append('<a href="%s">Panel showing read count per pathogen, per '
+                   'sample.</a>' % pathogenPanelFilename)
+            append('Red vertical bars indicate samples with an unusually high '
+                   'read count.')
+            append('</p>')
+
+        result.extend(proteinFieldsDescription)
+
+        # Write a linked table of contents by pathogen.
+        append('<p><span class="index-name">Pathogen index:</span>')
+        append('<span class="index">')
+        for genomeAccession in genomeAccessions:
+            genomeInfo = self._db.findGenome(genomeAccession)
+            append('<a href="#pathogen-%s">%s</a>' % (genomeAccession,
+                                                      genomeInfo['organism']))
+            append('&middot;')
+        # Get rid of final middle dot and add a period.
+        result.pop()
+        result[-1] += '.'
+        append('</span></p>')
+
+        # Write a linked table of contents by sample.
+        append('<p><span class="index-name">Sample index:</span>')
+        append('<span class="index">')
+        for sampleName in sampleNames:
+            append('<a href="#sample-%s">%s</a>' % (sampleName, sampleName))
+            append('&middot;')
+        # Get rid of final middle dot and add a period.
+        result.pop()
+        result[-1] += '.'
+        append('</span></p>')
+
+        # Write all pathogens (with samples (with proteins)).
+        append('<hr>')
+        append('<h1>Pathogens by sample</h1>')
+
+        for genomeAccession in genomeAccessions:
+            samples = self.genomeAccessions[genomeAccession]
+            sampleCount = len(samples)
+            genomeInfo = self._db.findGenome(genomeAccession)
+            pathogenProteinCount = genomeInfo['proteinCount']
+            if pathogenType == 'viral' and not omitVirusLinks:
+                quoted = quote(genomeInfo['organism'])
+                pathogenLinksHTML = (
+                    ' (<a href="%s%s">ICTV</a>, <a href="%s%s">ViralZone</a>)'
+                ) % (self.ICTV, quoted, self.VIRALZONE, quoted)
+            else:
+                pathogenLinksHTML = ''
+
+            withStr = (' with %d protein%s' %
+                       (pathogenProteinCount,
+                        '' if pathogenProteinCount == 1 else 's'))
+
+            pathogenReadsFilename = join(
+                self._pathogenDataDir,
+                'pathogen-%s.%s' % (genomeAccession, self._format))
+
+            pathogenReadsFp = open(pathogenReadsFilename, 'w')
+            pathogenReadCount = 0
+
+            append(
+                '<a id="pathogen-%s"></a>'
+                '<p class="pathogen">'
+                '<span class="pathogen-name">%s</span>'
+                '%s %s, '
+                'was matched by %d sample%s '
+                '(<a href="%s">%s</a> in total):'
+                '</p>' %
+                (genomeAccession,
+                 genomeInfo['organism'],
+                 pathogenLinksHTML, withStr,
+                 sampleCount, '' if sampleCount == 1 else 's',
+                 pathogenReadsFilename, self.READCOUNT_MARKER))
+
+            # Remember where we are in the output result so we can fill in
+            # the total read count once we have processed all samples for
+            # this pathogen. Not nice, I know.
+            pathogenReadCountLineIndex = len(result) - 1
+
+            for sampleName in sorted(samples):
+                readsFileName = self.pathogenSampleFiles.lookup(
+                    genomeAccession, sampleName)
+
+                # Copy the read data from the per-sample reads for this
+                # pathogen into the per-pathogen file of reads.
+                with open(readsFileName) as readsFp:
+                    while True:
+                        data = readsFp.read(4096)
+                        if data:
+                            pathogenReadsFp.write(data)
+                        else:
+                            break
+
+                proteins = samples[sampleName]['proteins']
+                proteinCount = len(proteins)
+                uniqueReadCount = samples[sampleName]['uniqueReadCount']
+                pathogenReadCount += uniqueReadCount
+                proteinCountHTML = '%d protein%s, ' % (
+                    proteinCount, '' if proteinCount == 1 else 's')
+
+                append(
+                    '<p class="sample indented">'
+                    'Sample <a href="#sample-%s">%s</a> '
+                    '(%s<a href="%s">%d de-duplicated (by id) '
+                    'read%s</a>, <a href="%s">panel</a>):</p>' %
+                    (sampleName, sampleName,
+                     proteinCountHTML,
+                     readsFileName,
+                     uniqueReadCount, '' if uniqueReadCount == 1 else 's',
+                     self.sampleNames[sampleName]))
+                append('<ul class="protein-list indented">')
+                for proteinName in sorted(proteins):
+                    proteinMatch = proteins[proteinName]
+                    append(
+                        '<li>'
+                        '<span class="stats">'
+                        '%(coverage).2f %(medianScore)6.2f %(bestScore)6.2f '
+                        '%(readAndHspCountStr)11s %(proteinLength)4d'
+                        % proteinMatch
+                    )
+
+                    if self._saveReadLengths:
+                        append(' (%s)' % ', '.join(
+                            map(str, sorted(proteinMatch['readLengths']))))
+
+                    # Add the </span> with no intermediate whitespace
+                    # because the 'stats' CSS class uses 'pre' on
+                    # whitespace, which results in a newline when we use
+                    # '\n'.join(result).
+                    result[-1] += '</span>'
+
+                    append(
+                        '<span class="protein-name">'
+                        '%(proteinName)s'
+                        '</span> '
+                        '(<a href="%(bluePlotFilename)s">blue plot</a>, '
+                        '<a href="%(readsFilename)s">reads</a>'
+                        % proteinMatch)
+
+                    if proteinMatch['proteinURL']:
+                        # Append this directly to the last string in result, to
+                        # avoid introducing whitespace when we join result
+                        # using '\n'.
+                        result[-1] += (', <a href="%s">NCBI protein</a>' %
+                                       proteinMatch['proteinURL'])
+
+                    if proteinMatch['genomeURL']:
+                        # Append this directly to the last string in result, to
+                        # avoid introducing whitespace when we join result
+                        # using '\n'.
+                        result[-1] += (', <a href="%s">NCBI genome</a>' %
+                                       proteinMatch['genomeURL'])
+                    result[-1] += ')'
+
+                    append('</li>')
+
+                append('</ul>')
+
+            pathogenReadsFp.close()
+
+            # Sanity check there's a read count marker text in our output
+            # where we expect it.
+            readCountLine = result[pathogenReadCountLineIndex]
+            if readCountLine.find(self.READCOUNT_MARKER) == -1:
+                raise ValueError(
+                    'Could not find pathogen read count marker (%s) in result '
+                    'index %d text (%s).' %
+                    (self.READCOUNT_MARKER, pathogenReadCountLineIndex,
+                     readCountLine))
+
+            # Put the read count into the pathogen summary line we wrote
+            # earlier, replacing the read count marker with the correct
+            # text.
+            result[pathogenReadCountLineIndex] = readCountLine.replace(
+                self.READCOUNT_MARKER,
+                '%d read%s' % (pathogenReadCount,
+                               '' if pathogenReadCount == 1 else 's'))
+
+        # Write all samples (with pathogens (with proteins)).
+        append('<hr>')
+        append('<h1>Samples by pathogen</h1>')
+
+        for sampleName in sampleNames:
+            samplePathogenAccessions = [
+                accession for accession in self.genomeAccessions
+                if sampleName in self.genomeAccessions[accession]]
+
+            if len(samplePathogenAccessions):
+                append(
+                    '<a id="sample-%s"></a>'
+                    '<p class="sample">Sample '
+                    '<span class="sample-name">%s</span> '
+                    'matched proteins from %d pathogen%s, '
+                    '<a href="%s">panel</a>:</p>' %
+                    (sampleName, sampleName, len(samplePathogenAccessions),
+                     '' if len(samplePathogenAccessions) == 1 else 's',
+                     self.sampleNames[sampleName]))
+            else:
+                append(
+                    '<a id="sample-%s"></a>'
+                    '<p class="sample">Sample '
+                    '<span class="sample-name">%s</span> '
+                    'did not match anything.</p>' %
+                    (sampleName, sampleName))
+                continue
+
+            for genomeAccession in sorted(samplePathogenAccessions):
+                genomeInfo = self._db.findGenome(genomeAccession)
+                readsFileName = self.pathogenSampleFiles.lookup(
+                    genomeAccession, sampleName)
+                proteins = self.genomeAccessions[genomeAccession][sampleName][
+                    'proteins']
+                uniqueReadCount = self.genomeAccessions[
+                    genomeAccession][sampleName]['uniqueReadCount']
+                proteinCount = len(proteins)
+                pathogenProteinCount = genomeInfo['proteinCount']
+                proteinCountStr = '%d/%d protein%s' % (
+                    proteinCount, pathogenProteinCount,
+                    '' if pathogenProteinCount == 1 else 's')
+
+                append(
+                    '<p class="sample indented">'
+                    '<a href="#pathogen-%s">%s</a> %s, '
+                    '<a href="%s">%d de-duplicated (by id) read%s</a>:</p>' %
+                    (genomeAccession, genomeInfo['organism'],
+                     proteinCountStr, readsFileName,
+                     uniqueReadCount, '' if uniqueReadCount == 1 else 's'))
+                append('<ul class="protein-list indented">')
+                for proteinAccession in sorted(proteins):
+                    proteinMatch = proteins[proteinAccession]
+                    append(
+                        '<li>'
+                        '<span class="stats">'
+                        '%(coverage).2f %(medianScore)6.2f %(bestScore)6.2f '
+                        '%(readAndHspCountStr)11s %(proteinLength)4d '
+                        '</span> '
+                        '<span class="protein-name">'
+                        '%(proteinName)s'
+                        '</span> '
+                        '(<a href="%(bluePlotFilename)s">blue plot</a>, '
+                        '<a href="%(readsFilename)s">reads</a>'
+                        % proteinMatch)
+
+                    if proteinMatch['proteinURL']:
+                        # Append this directly to the last string in result, to
+                        # avoid introducing whitespace when we join result
+                        # using '\n'.
+                        result[-1] += (', <a href="%s">NCBI protein</a>' %
+                                       proteinMatch['proteinURL'])
+
+                    if proteinMatch['genomeURL']:
+                        # Append this directly to the last string in result, to
+                        # avoid introducing whitespace when we join result
+                        # using '\n'.
+                        result[-1] += (', <a href="%s">NCBI genome</a>' %
+                                       proteinMatch['genomeURL'])
+
+                    result[-1] += ')'
+
+                    append('</li>')
+
+                append('</ul>')
+
+        append('</body>')
+        append('</html>')
+
+        return '\n'.join(result)
+
+    def _pathogenSamplePlot(self, genomeAccession, sampleNames, ax):
+        """
+        Make an image of a graph giving pathogen read count (Y axis) versus
+        sample id (X axis).
+
+        @param genomeAccession: A C{str} pathogen accession number.
+        @param sampleNames: A sorted C{list} of sample names.
+        @param ax: A matplotlib C{axes} instance.
+        """
+        readCounts = []
+        for sampleName in sampleNames:
+            try:
+                readCount = self.genomeAccessions[genomeAccession][sampleName][
+                    'uniqueReadCount']
+            except KeyError:
+                readCount = 0
+            readCounts.append(readCount)
+
+        highlight = 'r'
+        normal = 'gray'
+        sdMultiple = 2.5
+        minReadsForHighlighting = 10
+        highlighted = []
+
+        if len(readCounts) == 1:
+            if readCounts[0] > minReadsForHighlighting:
+                color = [highlight]
+                highlighted.append(sampleNames[0])
+            else:
+                color = [normal]
+        else:
+            mean = np.mean(readCounts)
+            sd = np.std(readCounts)
+            color = []
+            for readCount, sampleName in zip(readCounts, sampleNames):
+                if (readCount > (sdMultiple * sd) + mean and
+                        readCount >= minReadsForHighlighting):
+                    color.append(highlight)
+                    highlighted.append(sampleName)
+                else:
+                    color.append(normal)
+
+        nSamples = len(sampleNames)
+        x = np.arange(nSamples)
+        yMin = np.zeros(nSamples)
+        ax.set_xticks([])
+        ax.set_xlim((-0.5, nSamples - 0.5))
+        ax.vlines(x, yMin, readCounts, color=color)
+        if highlighted:
+            title = '%s\nIn red: %s' % (
+                genomeAccession, fill(', '.join(highlighted), 50))
+        else:
+            # Add a newline to keep the first line of each title at the
+            # same place as those titles that have an "In red:" second
+            # line.
+            title = genomeAccession + '\n'
+
+        ax.set_title(title, fontsize=10)
+        ax.tick_params(axis='both', which='major', labelsize=8)
+        ax.tick_params(axis='both', which='minor', labelsize=6)
+
+    def pathogenPanel(self, filename):
+        """
+        Make a panel of images, with each image being a graph giving pathogen
+        de-duplicated (by id) read count (Y axis) versus sample id (X axis).
+
+        @param filename: A C{str} file name to write the image to.
+        """
+        import matplotlib
+        matplotlib.use('PDF')
+        import matplotlib.pyplot as plt
+
+        self._computeUniqueReadCounts()
+        genomeAccessions = sorted(self.genomeAccessions)
+        sampleNames = sorted(self.sampleNames)
+
+        cols = 5
+        rows = int(len(genomeAccessions) / cols) + (
+            0 if len(genomeAccessions) % cols == 0 else 1)
+        figure, ax = plt.subplots(rows, cols, squeeze=False)
+
+        coords = dimensionalIterator((rows, cols))
+
+        for genomeAccession in genomeAccessions:
+            row, col = next(coords)
+            self._pathogenSamplePlot(genomeAccession, sampleNames,
+                                     ax[row][col])
+
+        # Hide the final panel graphs (if any) that have no content. We do
+        # this because the panel is a rectangular grid and some of the
+        # plots at the end of the last row may be unused.
+        for row, col in coords:
+            ax[row][col].axis('off')
+
+        figure.suptitle(
+            ('Per-sample read count for %d pathogen%s and %d sample%s.\n\n'
+             'Sample name%s: %s') % (
+                 len(genomeAccessions),
+                 '' if len(genomeAccessions) == 1 else 's',
+                 len(sampleNames),
+                 '' if len(sampleNames) == 1 else 's',
+                 '' if len(sampleNames) == 1 else 's',
+                 fill(', '.join(sampleNames), 50)),
+            fontsize=20)
+        figure.set_size_inches(5.0 * cols, 2.0 * rows, forward=True)
+        plt.subplots_adjust(hspace=0.4)
+
+        figure.savefig(filename)
+
+
+class _Genome(object):
+    """
+    Hold genome information, mirroring the attributes of a BioPython
+    GenBank record.
+
+    @param d: A C{dict} holding genome information (see below).
+    """
+    def __init__(self, d):
+        self.id = d['id']
+        self.description = d['name']
+        self.seq = d['sequence']
+        self.annotations = {
+            'taxonomy': d['taxonomy'],
+        }
+        self.lineage = d.get('lineage')
+        self.features = [_GenomeFeature(f) for f in d['features']]
+
+
+class _GenomeLocation(object):
+    """
+    Hold genome feature location information, mirroring the attributes of a
+    BioPython GenBank record.
+
+    @param start: An C{int} start location.
+    @param end: An C{int} stop location.
+    """
+    def __init__(self, start, end, strand):
+        self.start = start
+        self.end = end
+        self.strand = strand
+
+    def __str__(self):
+        return '[%d:%d](%s)' % (self.start, self.end,
+                                '+' if self.strand == 1 else '-')
+
+
+class _GenomeFeature(object):
+    """
+    Hold genome feature information, mirroring the attributes of a BioPython
+    GenBank record.
+
+    @param d: A C{dict} holding genome feature information.
+    """
+    def __init__(self, d):
+        self.type = d['type']
+        self.qualifiers = d['qualifiers']
+        self.strand = 1
+        location = d['qualifiers']['location']
+        self.location = _GenomeLocation(location['start'], location['stop'],
+                                        self.strand)
+
+
+class SqliteIndexWriter(object):
+    """
+    Create or update an Sqlite3 database holding information about proteins and
+    the genomes they come from.
+
+    @param dbFilename: A C{str} file name containing an sqlite3 database. If
+        the file does not exist it will be created. The special string
+        ":memory:" can be used to create an in-memory database.
+    @param fastaFp: A file-pointer to which the protein FASTA is written.
+    """
+    PROTEIN_ACCESSION_FIELD = 2
+    GENOME_ACCESSION_FIELD = 4
+    TAXONOMY_SEPARATOR = '\t'
+    SEQUENCE_ID_PREFIX = 'civ'
+    SEQUENCE_ID_SEPARATOR = '|'
+
+    def __init__(self, dbFilename, fastaFp=sys.stdout):
+        self._connection = sqlite3.connect(dbFilename)
+        self._fastaFp = fastaFp
+
+        cur = self._connection.cursor()
+        cur.executescript('''
+            CREATE TABLE IF NOT EXISTS proteins (
+                accession VARCHAR UNIQUE PRIMARY KEY,
+                genomeAccession VARCHAR NOT NULL,
+                sequence VARCHAR NOT NULL,
+                length INTEGER NOT NULL,
+                offsets VARCHAR NOT NULL,
+                forward INTEGER NOT NULL,
+                circular INTEGER NOT NULL,
+                rangeCount INTEGER NOT NULL,
+                gene VARCHAR,
+                note VARCHAR,
+                product VARCHAR,
+                FOREIGN KEY (genomeAccession)
+                    REFERENCES genomes (accession)
+            );
+
+            CREATE TABLE IF NOT EXISTS genomes (
+                accession VARCHAR UNIQUE PRIMARY KEY,
+                organism VARCHAR NOT NULL,
+                name VARCHAR NOT NULL,
+                sequence VARCHAR NOT NULL,
+                length INTEGER NOT NULL,
+                proteinCount INTEGER NOT NULL,
+                host VARCHAR,
+                note VARCHAR,
+                taxonomy VARCHAR NOT NULL,
+                databaseName VARCHAR
+            );
+            ''')
+        self._connection.commit()
+
+    def addGenBankFile(self, filename, rnaOnly=False, databaseName=None,
+                       lineageFetcher=None, proteinSource='GENBANK',
+                       genomeSource='GENBANK', duplicationPolicy='error',
+                       logfp=None):
+        """
+        Add proteins from a GenBank file.
+
+        @param filename: A C{str} file name, with the file in GenBank format
+            (see https://www.ncbi.nlm.nih.gov/Sitemap/samplerecord.html).
+        @param rnaOnly: If C{True}, only include RNA viruses.
+        @param databaseName: A C{str} indicating the database the records
+            in C{filename} came from (e.g., 'refseq' or 'RVDB').
+        @param lineageFetcher: A function that returns taxonomy lineage
+            information in the tuple-of-tuples form returned by
+            C{dark.taxonomy.AccessionLineageFetcher.lineage}. Must be given if
+            C{rnaOnly} is C{True}.
+        @param proteinSource: A C{str} giving the source of the protein
+            accession number. This becomes part of the sequence id printed
+            in the protein FASTA output.
+        @param genomeSource: A C{str} giving the source of the genome
+            accession number. This becomes part of the sequence id printed
+            in the protein FASTA output.
+        @param duplicationPolicy: A C{str} indicating what to do if a
+            to-be-inserted accession number is already present in the database.
+            "error" results in a ValueError being raised, "ignore" means ignore
+            the duplicate. It should also be possible to update (i.e., replace)
+            but that is not supported yet.
+        @param logfp: If not C{None}, a file pointer to write verbose
+            progress output to.
+        @raise DatabaseDuplicationError: If a duplicate accession number is
+            encountered and C{duplicationPolicy} is 'error'.
+        @return: A tuple containing two C{int}s: the number of genome sequences
+            in the added file and the total number of proteins found.
+        """
+        if rnaOnly:
+            assert lineageFetcher, (
+                'a lineage fetcher must be passed if rnaOnly is True.')
+
+        assert self.SEQUENCE_ID_SEPARATOR not in proteinSource, (
+            'proteinSource cannot contain %r as that is used as a separator.' %
+            self.SEQUENCE_ID_SEPARATOR)
+
+        assert self.SEQUENCE_ID_SEPARATOR not in genomeSource, (
+            'genomeSource cannot contain %r as that is used as a separator.' %
+            self.SEQUENCE_ID_SEPARATOR)
+
+        genomeCount = totalProteinCount = 0
+
+        with open(filename) as fp:
+            genomes = SeqIO.parse(fp, 'gb')
+            with self._connection:
+                for genome in genomes:
+                    if logfp:
+                        print('\n%s: %s' % (genome.id, genome.description),
+                              file=logfp)
+                        for k, v in genome.annotations.items():
+                            if k not in ('references', 'comment',
+                                         'structured_comment'):
+                                print('  %s = %r' % (k, v), file=logfp)
+                    if rnaOnly:
+                        try:
+                            lineage = lineageFetcher(genome.id)
+                        except ValueError as e:
+                            print('ValueError calling lineage fetcher: %s' % e,
+                                  file=sys.stderr)
+                            lineage = None
+
+                        if lineage:
+                            print('  Lineage:', file=logfp)
+                            print(formatLineage(lineage, prefix='    '),
+                                  file=logfp)
+                            if not isRNAVirus(lineage):
+                                if logfp:
+                                    print('  %s is not an RNA virus. '
+                                          'Skipping.' % genome.description,
+                                          file=logfp)
+                                continue
+                            else:
+                                if logfp:
+                                    print('  %s is an RNA virus.' %
+                                          genome.description, file=logfp)
+                        else:
+                            print('Could not look up taxonomy lineage for %s '
+                                  '(%s). Cannot confirm as RNA. Skipping.' % (
+                                      genome.id, genome.description),
+                                  file=sys.stderr)
+                            continue
+
+                    proteinCount = len(list(self._genomeProteins(genome)))
+
+                    if self.addGenome(
+                            genome, proteinCount, databaseName,
+                            duplicationPolicy=duplicationPolicy, logfp=logfp):
+
+                        self.addProteins(
+                            genome, proteinSource=proteinSource,
+                            genomeSource=genomeSource,
+                            duplicationPolicy=duplicationPolicy, logfp=logfp)
+
+                        totalProteinCount += proteinCount
+                        genomeCount += 1
+
+        return genomeCount, totalProteinCount
+
+    def addJSONFile(self, filename, rnaOnly=False, databaseName=None,
+                    lineageFetcher=None, proteinSource='GENBANK',
+                    genomeSource='GENBANK', duplicationPolicy='error',
+                    logfp=None):
+        """
+        Add proteins from a JSON infor file.
+
+        @param filename: A C{str} file name, in JSON format.
+        @param rnaOnly: If C{True}, only include RNA viruses.
+        @param databaseName: A C{str} indicating the database the records
+            in C{filename} came from (e.g., 'refseq' or 'RVDB').
+        @param lineageFetcher: A function that returns taxonomy lineage
+            information in the tuple-of-tuples form returned by
+            C{dark.taxonomy.AccessionLineageFetcher.lineage}. Must be given if
+            C{rnaOnly} is C{True}.
+        @param proteinSource: A C{str} giving the source of the protein
+            accession number. This becomes part of the sequence id printed
+            in the protein FASTA output.
+        @param genomeSource: A C{str} giving the source of the genome
+            accession number. This becomes part of the sequence id printed
+            in the protein FASTA output.
+        @param duplicationPolicy: A C{str} indicating what to do if a
+            to-be-inserted accession number is already present in the database.
+            "error" results in a ValueError being raised, "ignore" means ignore
+            the duplicate. It should also be possible to update (i.e., replace)
+            but that is not supported yet.
+        @param logfp: If not C{None}, a file pointer to write verbose
+            progress output to.
+        @raise DatabaseDuplicationError: If a duplicate accession number is
+            encountered and C{duplicationPolicy} is 'error'.
+        @return: A tuple containing two C{int}s: the number of genome sequences
+            in the added file and the total number of proteins found.
+        """
+        if rnaOnly:
+            assert lineageFetcher, (
+                'a lineage fetcher must be passed if rnaOnly is True.')
+
+        assert self.SEQUENCE_ID_SEPARATOR not in proteinSource, (
+            'proteinSource cannot contain %r as that is used as a separator.' %
+            self.SEQUENCE_ID_SEPARATOR)
+
+        assert self.SEQUENCE_ID_SEPARATOR not in genomeSource, (
+            'genomeSource cannot contain %r as that is used as a separator.' %
+            self.SEQUENCE_ID_SEPARATOR)
+
+        genomeCount = totalProteinCount = 0
+
+        with open(filename) as fp:
+            genome = _Genome(load(fp))
+
+        with self._connection:
+            if logfp:
+                print('\n%s: %s' % (genome.id, genome.description),
+                      file=logfp)
+                for k, v in genome.annotations.items():
+                    if k not in ('references', 'comment',
+                                 'structured_comment'):
+                        print('  %s = %r' % (k, v), file=logfp)
+            if rnaOnly:
+                if genome.lineage:
+                    lineage = genome.lineage
+                else:
+                    try:
+                        lineage = lineageFetcher(genome.id)
+                    except ValueError as e:
+                        print('ValueError calling lineage fetcher: %s' % e,
+                              file=sys.stderr)
+                        lineage = None
+
+                if lineage:
+                    print('  Lineage:', file=logfp)
+                    print(formatLineage(lineage, prefix='    '), file=logfp)
+                    if not isRNAVirus(lineage):
+                        if logfp:
+                            print('  %s is not an RNA virus. Skipping.' %
+                                  genome.description, file=logfp)
+                        return (0, 0)
+                    else:
+                        if logfp:
+                            print('  %s is an RNA virus.' %
+                                  genome.description, file=logfp)
+                else:
+                    print('Could not find taxonomy lineage for %s '
+                          '(%s). Cannot confirm as RNA. Skipping.' % (
+                              genome.id, genome.description),
+                          file=sys.stderr)
+                    return (0, 0)
+
+            proteinCount = len(list(self._genomeProteins(genome)))
+
+            if self.addGenome(
+                    genome, proteinCount, databaseName,
+                    duplicationPolicy=duplicationPolicy, logfp=logfp):
+
+                self.addProteins(
+                    genome, proteinSource=proteinSource,
+                    genomeSource=genomeSource,
+                    duplicationPolicy=duplicationPolicy, logfp=logfp)
+
+                totalProteinCount += proteinCount
+                genomeCount += 1
+
+        return genomeCount, totalProteinCount
+
+    def addGenome(self, genome, proteinCount, databaseName,
+                  duplicationPolicy='error', logfp=None):
+        """
+        Add information about a genome to the genomes table.
+
+        @param genome: A GenBank genome record, as parsed by SeqIO.parse
+        @param proteinCount: The C{int} number of proteins in the genome.
+        @param databaseName: A C{str} indicating the database the records
+            in C{filename} came from (e.g., 'refseq' or 'RVDB').
+        @param duplicationPolicy: A C{str} indicating what to do if a
+            to-be-inserted accession number is already present in the database.
+            "error" results in a ValueError being raised, "ignore" means ignore
+            the duplicate. It should also be possible to update (i.e., replace)
+            but that is not supported yet.
+        @param logfp: If not C{None}, a file pointer to write verbose
+            progress output to.
+        @raise DatabaseDuplicationError: If a duplicate accession number is
+            encountered and C{duplicationPolicy} is 'error'.
+        @return: C{True} if the genome was added, else C{False}.
+        """
+        sequence = str(genome.seq)
+        taxonomy = self.TAXONOMY_SEPARATOR.join(genome.annotations['taxonomy'])
+        source = self._sourceInfo(genome, logfp=logfp)
+
+        if source is None:
+            # The lack of a source is logged by self._sourceInfo.
+            return False
+
+        try:
+            self._connection.execute(
+                'INSERT INTO genomes(accession, organism, name, sequence, '
+                'length, proteinCount, host, note, taxonomy, databaseName) '
+                'VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?)',
+                (genome.id, source['organism'], genome.description, sequence,
+                 len(sequence), proteinCount, source['host'],
+                 source.get('note'), taxonomy, databaseName))
+        except sqlite3.IntegrityError as e:
+            if str(e).find('UNIQUE constraint failed') > -1:
+                if duplicationPolicy == 'error':
+                    raise DatabaseDuplicationError(
+                        'Genome information for %r already present in '
+                        'database.' % genome.id)
+                elif duplicationPolicy == 'ignore':
+                    if logfp:
+                        print(
+                            'Genome information for %r already present in '
+                            'database. Ignoring.' % genome.id, file=logfp)
+                    return False
+                else:
+                    raise NotImplementedError(
+                        'Unknown duplication policy (%s) found when '
+                        'attempting to insert genome information for %s.' %
+                        (duplicationPolicy, genome.id))
+            else:
+                raise
+        else:
+            return True
+
+    def addProtein(self, accession, genomeAccession, sequence, offsets,
+                   forward, circular, rangeCount, gene=None, note=None,
+                   product=None, duplicationPolicy='error', logfp=None):
+        """
+        Add information about a protein to the proteins table.
+
+        @param accession: A C{str} protein accession id.
+        @param genomeAccession: A C{str} genome accession id (the genome to
+            which this protein belongs).
+        @param sequence: A C{str} protein amino acid sequence.
+        @param offsets: A C{str} describing the offsets of the protein in the
+            genome (as obtained from C{SeqIO.parse} on a GenBank file).
+        @param forward: A C{bool}, C{True} if the protein occurs on the
+            forward strand of the genome, C{False} if on the complement strand.
+            Note that this is converted to an C{int} in the database.
+        @param circular: A C{bool}, C{True} if the protein crosses the genome
+            boundary and is therefore circular, C{False} if not. Note that
+            this is converted to an C{int} in the database.
+        @param rangeCount: The C{int} number of ranges (regions) the protein
+            comes from in the genome.
+        @param gene: A C{str} gene name, or C{None} if no gene is known.
+        @param note: A C{str} note about the protein, or C{None}.
+        @param product: A C{str} description of the protein product (e.g.,
+            "putative replication initiation protein"), or C{None}.
+        @param duplicationPolicy: A C{str} indicating what to do if a
+            to-be-inserted accession number is already present in the database.
+            "error" results in a ValueError being raised, "ignore" means ignore
+            the duplicate. It should also be possible to update (i.e., replace)
+            but that is not supported yet.
+        @param logfp: If not C{None}, a file pointer to write verbose
+            progress output to.
+        @raise DatabaseDuplicationError: If a duplicate accession number is
+            encountered and C{duplicationPolicy} is 'error'.
+        """
+        try:
+            self._connection.execute(
+                'INSERT INTO proteins('
+                'accession, genomeAccession, sequence, length, offsets, '
+                'forward, circular, rangeCount, gene, note, product) '
+                'VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?)',
+                (accession, genomeAccession, sequence, len(sequence), offsets,
+                 int(forward), int(circular), rangeCount, gene, note, product))
+        except sqlite3.IntegrityError as e:
+            if str(e).find('UNIQUE constraint failed') > -1:
+                if duplicationPolicy == 'error':
+                    raise DatabaseDuplicationError(
+                        'Protein information for %r already present in '
+                        'database.' % accession)
+                elif duplicationPolicy == 'ignore':
+                    if logfp:
+                        print(
+                            'Protein information for %r already present in '
+                            'database. Ignoring.' % accession, file=logfp)
+                else:
+                    raise NotImplementedError(
+                        'Unknown duplication policy (%s) found when '
+                        'attempting to insert protein information for %s.' %
+                        (duplicationPolicy, accession))
+            else:
+                raise
+        else:
+            if logfp:
+                print('    Protein %s: genome=%s product=%s' % (
+                    accession, genomeAccession, product), file=logfp)
+
+    def addProteins(self, genome, proteinSource='GENBANK',
+                    genomeSource='GENBANK', duplicationPolicy='error',
+                    logfp=None):
+        """
+        Add proteins from a Genbank genome record to the proteins database and
+        write out their sequences to the proteins FASTA file (in
+        C{self._fastaFp}).
+
+        @param genome: Either a GenBank genome record, as parsed by
+            C{SeqIO.parse} or a C{_Genome} instance (which behaves like the
+            former).
+        @param proteinSource: A C{str} giving the source of the protein
+            accession number. This becomes part of the sequence id printed
+            in the protein FASTA output.
+        @param genomeSource: A C{str} giving the source of the genome
+            accession number. This becomes part of the sequence id printed
+            in the protein FASTA output.
+        @param duplicationPolicy: A C{str} indicating what to do if a
+            to-be-inserted accession number is already present in the database.
+            "error" results in a ValueError being raised, "ignore" means ignore
+            the duplicate. It should also be possible to update (i.e., replace)
+            but that is not supported yet.
+        @param logfp: If not C{None}, a file pointer to write verbose
+            progress output to.
+        @raise DatabaseDuplicationError: If a duplicate accession number is
+            encountered and C{duplicationPolicy} is 'error'.
+        """
+        genomeLen = len(genome.seq)
+        source = self._sourceInfo(genome, logfp=logfp)
+        # source must be present. addGenome would skip this genome otherwise.
+        assert source
+
+        for fInfo in self._genomeProteins(genome, logfp=logfp):
+
+            # Write FASTA for the protein.
+            seqId = self.SEQUENCE_ID_SEPARATOR.join((
+                self.SEQUENCE_ID_PREFIX,
+                proteinSource, fInfo['proteinId'],
+                genomeSource, genome.id,
+                fInfo['product']))
+
+            print('>%s [%s]\n%s' %
+                  (seqId, source['organism'], fInfo['translation']),
+                  file=self._fastaFp)
+
+            self.addProtein(
+                fInfo['proteinId'], genome.id, fInfo['translation'],
+                fInfo['featureLocation'], fInfo['forward'],
+                fInfo['circular'],
+                fInfo['ranges'].distinctRangeCount(genomeLen),
+                gene=fInfo['gene'], note=fInfo['note'],
+                product=fInfo['product'], duplicationPolicy=duplicationPolicy,
+                logfp=logfp)
+
+    def _sourceInfo(self, genome, logfp):
+        """
+        Extract summary information from a genome source feature.
+
+        @param genome: A GenBank genome record, as parsed by SeqIO.parse
+        @param logfp: If not C{None}, a file pointer to write verbose
+            progress output to.
+        @return: A C{dict} with keys for the various pieces of information
+            (if any) found in the source feature (see the return value below
+            for detail). Or C{None} if no source feature is found.
+        """
+        result = {}
+
+        for feature in genome.features:
+            if feature.type == 'source':
+                for key in 'host', 'note', 'organism':
+                    try:
+                        values = feature.qualifiers[key]
+                    except KeyError:
+                        value = None
+                    else:
+                        assert len(values) == 1
+                        value = values[0]
+
+                    result[key] = value
+                break
+        else:
+            warn('Genome %r (accession %s) had no source feature! '
+                 'Skipping.' % (genome.description, genome.id))
+            return
+
+        return result
+
+    def _cdsInfo(self, genome, feature, logfp=None):
+        """
+        Extract summary information from a genome CDS feature.
+
+        @param genome: A GenBank genome record, as parsed by SeqIO.parse
+        @param feature: A feature from a genome, as produced by BioPython's
+            GenBank parser.
+        @param logfp: If not C{None}, a file pointer to write verbose
+            progress output to.
+        @return: A C{dict} with keys for the various pieces of information
+            found in the feature (see the return value below for detail).
+            Or C{None} if the feature is not of interest or otherwise invalid.
+        """
+        qualifiers = feature.qualifiers
+
+        # Check in advance that all feature qualifiers we're interested in
+        # have the right lengths, if they're present.
+        for key in 'gene', 'note', 'product', 'protein_id', 'translation':
+            if key in qualifiers:
+                assert len(qualifiers[key]) == 1, (
+                    'GenBank qualifier key %s is not length one %r' %
+                    (key, qualifiers[key]))
+
+        # A protein id is mandatory.
+        if 'protein_id' in qualifiers:
+            proteinId = qualifiers['protein_id'][0]
+        else:
+            if 'translation' in qualifiers:
+                warn('Genome %r (accession %s) has CDS feature with no '
+                     'protein_id feature but has a translation! '
+                     'Skipping.\nFeature: %s' %
+                     (genome.description, genome.id, feature))
+            return
+
+        # A translated (i.e., amino acid) sequence is mandatory.
+        if 'translation' in qualifiers:
+            translation = qualifiers['translation'][0]
+        else:
+            warn('Genome %r (accession %s) has CDS feature with protein '
+                 '%r with no translated sequence. Skipping.' %
+                 (genome.description, genome.id, proteinId))
+            return
+
+        featureLocation = str(feature.location)
+
+        # Make sure the feature's location string can be parsed.
+        try:
+            ranges = GenomeRanges(featureLocation)
+        except ValueError as e:
+            warn('Genome %r  (accession %s) contains unparseable CDS '
+                 'location for protein %r. Skipping. Error: %s' %
+                 (genome.description, genome.id, proteinId, e))
+            return
+        else:
+            # Does the protein span the end of the genome? This indicates a
+            # circular genome.
+            circular = int(ranges.circular(len(genome.seq)))
+
+        if feature.location.start >= feature.location.end:
+            warn('Genome %r (accession %s) contains feature with start '
+                 '(%d) >= stop (%d). Skipping.\nFeature: %s' %
+                 (genome.description, genome.id, feature.location.start,
+                  feature.location.end, feature))
+            return
+
+        strand = feature.strand
+        if strand is None:
+            # The strands of the protein in the genome are not all the same
+            # (see Bio.SeqFeature.CompoundLocation._get_strand).  The
+            # protein is formed by the combination of reading one strand in
+            # one direction and the other in the other direction.
+            #
+            # This occurs just once in all 1.17M proteins found in all 700K
+            # RVDB (C-RVDBv15.1) genomes, for protein YP_656697.1 on the
+            # Ranid herpesvirus 1 strain McKinnell genome (NC_008211.1).
+            #
+            # This situation makes turning DIAMOND protein output into
+            # SAM very complicated because a match on such a protein
+            # cannot be stored as a SAM linear alignment. It instead
+            # requires a multi-line 'supplementary' alignment. The code
+            # and tests for that are more complex than I want to deal
+            # with at the moment, just for the sake of one protein in a
+            # frog herpesvirus.
+            warn('Genome %s (accession %s) has protein %r with mixed '
+                 'orientation!' % (genome.description, genome.id,
+                                   proteinId))
+            return
+        elif strand == 0:
+            # This never occurs for proteins corresponding to genomes in
+            # the RVDB database C-RVDBv15.1.
+            warn('Genome %r (accession %s) has protein %r with feature '
+                 'with strand of zero!' %
+                 (genome.description, genome.id, proteinId))
+            return
+        else:
+            assert strand in (1, -1)
+            forward = strand == 1
+            # Make sure the strand agrees with the orientations in the
+            # string BioPython makes out of the locations.
+            assert ranges.orientations() == {forward}
+
+        return {
+            'circular': circular,
+            'featureLocation': featureLocation,
+            'forward': forward,
+            'gene': qualifiers.get('gene', [''])[0],
+            'note': qualifiers.get('note', [''])[0],
+            'product': qualifiers.get('product', ['UNKNOWN'])[0],
+            'proteinId': proteinId,
+            'ranges': ranges,
+            'strand': strand,
+            'translation': translation,
+        }
+
+    def _genomeProteins(self, genome, logfp=None):
+        """
+        Get proteins (CDS features) that we can process from a genome, along
+        with information extracted from each.
+
+        @param genome: A GenBank genome record, as parsed by SeqIO.parse
+        @param logfp: If not C{None}, a file pointer to write verbose
+            progress output to.
+        @return: A generator yielding feature info C{dict}s as returned by
+            C{self._cdsInfo}.
+        """
+        for feature in genome.features:
+            if feature.type == 'CDS':
+                featureInfo = self._cdsInfo(genome, feature, logfp=None)
+                if featureInfo:
+                    yield featureInfo
+
+    def close(self):
+        """
+        Create indices on the accesssion ids and close the connection.
+        """
+        cur = self._connection.cursor()
+        cur.execute('CREATE UNIQUE INDEX IF NOT EXISTS protein_idx ON '
+                    'proteins(accession)')
+        cur.execute('CREATE UNIQUE INDEX IF NOT EXISTS genomes_idx ON '
+                    'genomes(accession)')
+        self._connection.commit()
+        self._connection.close()
+        self._connection = None
+
+    def __enter__(self):
+        return self
+
+    def __exit__(self, excType, excValue, traceback):
+        self.close()
+
+
+class SqliteIndex(object):
+    """
+    Provide lookup access to an Sqlite3 database holding information about
+    proteins and the genomes they come from.
+
+    @param dbFilenameOrConnection: Either a C{str} file name containing an
+        sqlite3 database as created by C{SqliteIndexWriter} or an already
+        open connection to such a database. Note that an already open
+        connection will not be closed by self.close().
+    @param lookupCacheSize: The C{int} size of the memoization cache
+        for the protein and genome lookup functions (each has its own
+        memoization cache).
+    """
+    PROTEIN_ACCESSION_FIELD = 2
+    GENOME_ACCESSION_FIELD = 4
+
+    def __init__(self, dbFilenameOrConnection, lookupCacheSize=1024):
+        if isinstance(dbFilenameOrConnection, string_types):
+            self._connection = sqlite3.connect(dbFilenameOrConnection)
+            self._closeConnection = True
+        else:
+            self._connection = dbFilenameOrConnection
+            self._closeConnection = False
+        self._connection.row_factory = sqlite3.Row
+        self._proteinCache = LRUCache(maxsize=lookupCacheSize)
+        self._genomeCache = LRUCache(maxsize=lookupCacheSize)
+
+    def genomeAccession(self, id_):
+        """
+        Get the genome accession info from a sequence id.
+
+        @param id_: A C{str} sequence id in the form
+            'civ|GENBANK|%s|GENBANK|%s|%s [%s]' where the genome accession
+            is in the fifth '|'-separated field.
+        @raise IndexError: If C{id_} does not have enough |-separated fields.
+        @return: The C{str} accession number.
+        """
+        return id_.split('|', self.GENOME_ACCESSION_FIELD + 1)[
+            self.GENOME_ACCESSION_FIELD]
+
+    def proteinAccession(self, id_):
+        """
+        Get the protein accession info from a sequence id.
+
+        @param id_: A C{str} sequence id in the form
+            'civ|GENBANK|%s|GENBANK|%s|%s [%s]' where the protein accession
+            is in the third '|'-separated field.
+        @raise IndexError: If C{id_} does not have enough |-separated fields.
+        @return: The C{str} accession number.
+        """
+        return id_.split('|', self.PROTEIN_ACCESSION_FIELD + 1)[
+            self.PROTEIN_ACCESSION_FIELD]
+
+    @cachedmethod(attrgetter('_genomeCache'))
+    def _findGenome(self, accession):
+        """
+        Find info about a genome, given an accession number.
+
+        @param accession: A C{str} accession number.
+        @return: A C{dict} with keys corresponding to the names of the columns
+            in the genomes database table, else C{None} if C{id_} cannot be
+            found.
+        """
+        cur = self.execute(
+            'SELECT * FROM genomes WHERE accession = ?', (accession,))
+        row = cur.fetchone()
+
+        if row:
+            result = dict(row)
+            result['accession'] = accession
+            result['taxonomy'] = result['taxonomy'].split(
+                SqliteIndexWriter.TAXONOMY_SEPARATOR)
+            return result
+
+    def findGenome(self, id_):
+        """
+        Find info about a genome, given a sequence id.
+
+        @param id_: A C{str} sequence id. This is either of the form
+            'civ|GENBANK|%s|GENBANK|%s|%s [%s]' where the genome id is in the
+            5th '|'-delimited field, or else is the nucleotide sequence
+            accession number as already extracted.
+        @return: A C{dict} with keys corresponding to the names of the columns
+            in the genomes database table, else C{None} if C{id_} cannot be
+            found.
+        """
+        try:
+            accession = self.genomeAccession(id_)
+        except IndexError:
+            accession = id_
+
+        return self._findGenome(accession)
+
+    @cachedmethod(attrgetter('_proteinCache'))
+    def _findProtein(self, accession):
+        """
+        Find info about a protein, given an accession number.
+
+        @param accession: A C{str} accession number.
+        @return: A C{dict} with keys corresponding to the names of the columns
+            in the proteins database table, else C{None} if C{id_} cannot be
+            found.
+        """
+        cur = self.execute(
+            'SELECT * FROM proteins WHERE accession = ?', (accession,))
+        row = cur.fetchone()
+        if row:
+            result = dict(row)
+            result['forward'] = bool(result['forward'])
+            result['circular'] = bool(result['circular'])
+            result['length'] = int(result['length'])
+            result['accession'] = accession
+            return result
+
+    def findProtein(self, id_):
+        """
+        Find info about a protein, given a sequence id.
+
+        @param id_: A C{str} sequence id. This is either of the form
+            'civ|GENBANK|%s|GENBANK|%s|%s [%s]' where the protein id is in the
+            3rd '|'-delimited field, or else is the protein accession number as
+            already extracted.
+        @return: A C{dict} with keys corresponding to the names of the columns
+            in the proteins database table, else C{None} if C{id_} cannot be
+            found.
+        """
+        try:
+            accession = self.proteinAccession(id_)
+        except IndexError:
+            accession = id_
+
+        return self._findProtein(accession)
+
+    def execute(self, query, *args):
+        """
+        Execute an SQL statement. See
+        https://docs.python.org/3.5/library/sqlite3.html#sqlite3.Cursor.execute
+        for full argument details.
+
+        @param query: A C{str} SQL query.
+        @param args: Additional arguments (if any) to pass to the sqlite3
+            execute command.
+        @return: An sqlite3 cursor.
+        """
+        cur = self._connection.cursor()
+        cur.execute(query, *args)
+        return cur
+
+    def close(self):
+        """
+        Close the database connection (if we opened it).
+        """
+        if self._closeConnection:
+            self._connection.close()
+        self._connection = None
+
+    def __enter__(self):
+        return self
+
+    def __exit__(self, excType, excValue, traceback):
+        self.close()
diff --git a/dark/diamond/conversion.py b/dark/diamond/conversion.py
index 8dbca3f..234545a 100644
--- a/dark/diamond/conversion.py
+++ b/dark/diamond/conversion.py
@@ -17,8 +17,8 @@ from dark.diamond.hsp import normalizeHSP
 # command line via --outfmt 6) that must be given to DIAMOND blastx to
 # allow its output to be parsed by convert-diamond-to-sam.py (which uses
 # diamondTabularFormatToDicts (below)).
-FIELDS = ('bitscore btop qframe qend qqual qlen qseq qseqid qstart slen '
-          'sstart stitle')
+FIELDS = ('bitscore btop qframe qend full_qqual qlen full_qseq qseqid '
+          'qstart slen sstart stitle')
 
 # The keys in the following are DIAMOND format 6 field names. The values
 # are one-argument functions that take a string and return an appropriately
@@ -52,55 +52,70 @@ DIAMOND_FIELD_CONVERTER = {
 }
 
 
-def diamondTabularFormatToDicts(filename, fieldNames=None):
+class DiamondTabularFormat(object):
     """
-    Read DIAMOND tabular (--outfmt 6) output and convert lines to dictionaries.
+    Read/convert DIAMOND TAB-separated format (format #6).
 
-    @param filename: Either a C{str} file name or an open file pointer.
     @param fieldNames: A C{list} or C{tuple} of C{str} DIAMOND field names.
         Run 'diamond -help' to see the full list. If C{None}, a default set of
         fields will be used, as compatible with convert-diamond-to-sam.py
-    @raise ValueError: If a line of C{filename} does not have the expected
-        number of TAB-separated fields (i.e., len(fieldNames)). Or if
-        C{fieldNames} is empty or contains duplicates.
-    @return: A generator that yields C{dict}s with keys that are the DIAMOND
-        field names and values as converted by DIAMOND_FIELD_CONVERTER.
     """
-    fieldNames = fieldNames or FIELDS.split()
-    nFields = len(fieldNames)
-    if not nFields:
-        raise ValueError('fieldNames cannot be empty.')
-
-    c = Counter(fieldNames)
-    if c.most_common(1)[0][1] > 1:
-        raise ValueError(
-            'fieldNames contains duplicated names: %s.' %
-            (', '.join(sorted(x[0] for x in c.most_common() if x[1] > 1))))
-
-    def identity(x):
-        return x
-
-    convertFunc = DIAMOND_FIELD_CONVERTER.get
-
-    with as_handle(filename) as fp:
-        for count, line in enumerate(fp, start=1):
-            result = {}
-            line = line[:-1]
-            values = line.split('\t')
-            if len(values) != nFields:
-                raise ValueError(
-                    'Line %d of %s had %d field values (expected %d). '
-                    'To provide input for this function, DIAMOND must be '
-                    'called with "--outfmt 6 %s" (without the quotes). '
-                    'The offending input line was %r.' %
-                    (count,
-                     (filename if isinstance(filename, six.string_types)
-                      else 'input'),
-                     len(values), nFields, FIELDS, line))
-            for fieldName, value in zip(fieldNames, values):
-                value = convertFunc(fieldName, identity)(value)
-                result[fieldName] = value
-            yield result
+    def __init__(self, fieldNames=None):
+        self._fieldNames = fieldNames or FIELDS.split()
+        self._nFields = len(self._fieldNames)
+        if not self._nFields:
+            raise ValueError('fieldNames cannot be empty.')
+
+        c = Counter(self._fieldNames)
+        if c.most_common(1)[0][1] > 1:
+            raise ValueError(
+                'field names contains duplicated names: %s.' %
+                (', '.join(sorted(x[0] for x in c.most_common() if x[1] > 1))))
+
+        def identity(x):
+            return x
+
+        self._convertors = tuple(DIAMOND_FIELD_CONVERTER.get(field, identity)
+                                 for field in self._fieldNames)
+
+    def diamondFieldsToDict(self, line):
+        """
+        Convert a line of DIAMOND tabular (--outfmt 6) output to a dictionary
+
+        @param line: A C{str} line read from DIAMOND output. This may contain
+            a trailing whitespace (including newline and carriage return).
+        @raise ValueError: If a line of C{filename} does not have the expected
+            number of TAB-separated fields (i.e., len(self._fieldNames)).
+        @return: A C{dict} with keys that are the DIAMOND field names, with
+            values as converted by the functions set up in __init__.
+        """
+        values = line.rstrip().split('\t')
+
+        if len(values) != self._nFields:
+            raise ValueError(
+                'DIAMOND output line had %d field values (expected %d). '
+                'The offending input line was %r.' %
+                (len(values), self._nFields, line))
+
+        return dict(
+            (fieldName, func(value))
+            for fieldName, func, value in zip(
+                self._fieldNames, self._convertors, values))
+
+    def diamondTabularFormatToDicts(self, filename):
+        """
+        Read DIAMOND tabular (--outfmt 6) output and convert lines to
+        dictionaries.
+
+        @param filename: Either a C{str} file name or an open file pointer.
+        @raise ValueError: If a line of C{filename} does not have the expected
+            number of TAB-separated fields (i.e., len(fieldNames)).
+        @return: A generator that yields C{dict}s returned by
+            C{diamondFieldsToDict}, above.
+        """
+        with as_handle(filename) as fp:
+            for line in fp:
+                yield self.diamondFieldsToDict(line)
 
 
 class DiamondTabularFormatReader(object):
diff --git a/dark/diamond/run.py b/dark/diamond/run.py
new file mode 100644
index 0000000..80cc5e7
--- /dev/null
+++ b/dark/diamond/run.py
@@ -0,0 +1,109 @@
+from __future__ import print_function
+
+from tempfile import mkdtemp
+from shutil import rmtree
+from os.path import join
+from subprocess import CalledProcessError
+
+from dark.diamond.conversion import FIELDS, DiamondTabularFormat
+from dark.process import Executor
+from dark.utils import cd
+
+
+def diamondInstalled():
+    """
+    Test if DIAMOND is installed.
+
+    @return: A C{bool}, which is C{True} if DIAMOND seems to be installed.
+    """
+    try:
+        Executor().execute('diamond help')
+    except CalledProcessError:
+        return False
+    else:
+        return True
+
+
+class DiamondExecutor(object):
+    """
+
+    @param dryRun: If C{True} do not actually execute the DIAMOND commands.
+    """
+    SUBJECTS_FILENAME = 'subjects.fasta'
+    QUERIES_FILENAME = 'queries.fasta'
+    OUTPUT_FILENAME = 'diamond.tsv'
+
+    def __init__(self, dryRun=False):
+        self._dirty = False
+        self._dir = mkdtemp()
+        self._subjectsFp = None
+        self._subjectsExist = False
+        self._executor = Executor(dryRun)
+
+    def addSubject(self, subject):
+        """
+        Add a subject sequence to the database.
+
+        @param subject: A C{dark.reads.Read} instance.
+        """
+        if self._subjectsFp is None:
+            self._subjectsFp = open(join(self._dir, self.SUBJECTS_FILENAME),
+                                    'a', encoding='utf-8')
+
+        print(subject.toString('fasta'), end='', file=self._subjectsFp)
+        self._subjectsExist = self._dirty = True
+
+    def cleanup(self):
+        """
+        Remove the temporary directory we made.
+        """
+        if self._subjectsFp:
+            self._subjectsFp.close()
+            self._subjectsFp = None
+        rmtree(self._dir)
+
+    def search(self, reads, fieldNames=None):
+        """
+        Match reads against the database.
+
+        @param reads: An instance of C{dark.reads.Reads}.
+        @param fieldNames: An iterable of C{str} field names for DIAMOND
+            tabular output (format 6). See diamond help for the names of all
+            available fields.
+        @return: A generator that yields C{dict}s with keys as in
+            C{fieldNames}.
+        """
+        if not self._subjectsExist:
+            raise ValueError('No subject sequences in the database')
+
+        with cd(self._dir):
+            if self._dirty:
+                self._subjectsFp.close()
+                self._subjectsFp = None
+                self._executor.execute('diamond makedb --db database --in %s' %
+                                       self.SUBJECTS_FILENAME)
+
+            with open(self.QUERIES_FILENAME, 'w') as fp:
+                count = reads.save(fp, format_='fastq')
+
+            if count == 0:
+                raise ValueError('No query sequences were passed')
+
+            fieldNames = fieldNames or FIELDS.split()
+
+            self._executor.execute(
+                'diamond blastx --db database --query %s --outfmt 6 %s > %s' %
+                (self.QUERIES_FILENAME, ' '.join(fieldNames),
+                 self.OUTPUT_FILENAME))
+
+            dtf = DiamondTabularFormat(fieldNames)
+
+            for diamondDict in dtf.diamondTabularFormatToDicts(
+                    self.OUTPUT_FILENAME):
+                yield diamondDict
+
+    def __enter__(self):
+        return self
+
+    def __exit__(self, excType, excValue, traceback):
+        self.cleanup()
diff --git a/dark/diamond/sam.py b/dark/diamond/sam.py
new file mode 100644
index 0000000..89f0c30
--- /dev/null
+++ b/dark/diamond/sam.py
@@ -0,0 +1,377 @@
+from __future__ import print_function
+
+from os import unlink
+from tempfile import TemporaryFile
+from resource import getrlimit, RLIMIT_NOFILE
+
+from dark.btop import btop2cigar
+from dark.diamond.conversion import DiamondTabularFormat
+from dark.fpcache import FilePointerCache
+from dark.genbank import GenomeRanges
+from dark.reads import DNARead
+from dark.utils import asHandle
+
+
+class SAMWriter(object):
+    """
+    A class that collects SAM output lines and writes SAM to a file.
+
+    @param ram: If C{True}, do not use a temporary file to hold the non-header
+        SAM output. This will run faster but use more memory since all
+        non-header SAM output will be stored in RAM and only written out when
+        the full header can be determined.
+    """
+    def __init__(self, ram=False):
+        self._ram = ram
+        self._referenceLengths = {}
+
+    def addMatchLine(self, line, genome, fp):
+        """
+        Add a line of SAM output.
+
+        @param line: A C{str} line of valid pre-formatted SAM output.
+        @param genome: A C{dict} with information about the (nucleotide) genome
+            that the protein in the DIAMOND match comes from. The C{dict} is
+            as returned by C{dark.proteins.SqliteIndex.findGenome}.
+        @param fp: The file pointer to write the line to.
+        """
+        genomeName = genome['accession']
+        genomeLength = genome['length']
+        try:
+            preexistingLength = self._referenceLengths[genomeName]
+        except KeyError:
+            self._referenceLengths[genomeName] = genomeLength
+        else:
+            if preexistingLength != genomeLength:
+                raise ValueError(
+                    'Reference %r passed with length %d but has already been '
+                    'given with a different length (%d).' %
+                    (genomeName, genomeLength, preexistingLength))
+
+        print(line, file=fp)
+
+    def save(self, fp, tmpfp=None, headerLines=None):
+        """
+        Write SAM output to a file.
+
+        @param fp: An open file pointer to write the full SAM to.
+        @param tmpfp: An open file pointer with the non-header SAM output.
+            Only used if C{self._ram} is C{False}
+        @param headerLines: An optional iterable of C{str}s SAM header
+            lines (with no trailing newlines).
+        """
+        # Print SAM headers.
+        print('\n'.join(
+            (headerLines or []) +
+            [
+                '@SQ\tSN:%s\tLN:%d' % (name, self._referenceLengths[name])
+                for name in sorted(self._referenceLengths)
+            ]), file=fp)
+
+        # Print SAM match lines.
+        if self._ram:
+            # Trade speed for memory and print the non-header lines one at a
+            # time instead of making and printing a potentially massive string
+            # with '\n'.join() (although arguably we could just make the single
+            # string, seeing as self._ram is True.)
+            for line in self._nonHeaderLines:
+                print(line, file=fp)
+        else:
+            tmpfp.seek(0)
+            for line in tmpfp:
+                print(line, end='', file=fp)
+
+
+class _DiamondSAMWriter(object):
+    """
+    A base class for converting DIAMOND tabular output to SAM.
+
+    @param genomesProteins: An open sqlite3 database with information on
+        genomes and proteins, as made by the dark-matter
+        make-protein-database.py script.
+    @param mappingQuality: The mapping quality to use for the SAM MAPQ field
+        (#5). The default (255) indicates that mapping quality information is
+        not available.
+    @param ram: If C{True}, do not use a temporary file to hold the non-header
+        SAM output. This will run faster but use more memory since all
+        non-header SAM output will be stored in RAM and only written out when
+        the full header can be determined.
+    @param keepDescriptions: If C{True}, do not discard text after the first
+        space in query or subject sequence ids. Note that this violates the
+        SAM specification, but since SAM files are TAB-separated there may be
+        only a small chance this will cause problems downstream.
+    """
+    def __init__(self, genomesProteins, mappingQuality=255, ram=False,
+                 keepDescriptions=False):
+        self._genomesProteins = genomesProteins
+        self._mappingQuality = mappingQuality
+        self._ram = ram
+        self._keepDescriptions = keepDescriptions
+        self._strToDiamondDict = DiamondTabularFormat().diamondFieldsToDict
+        self._referenceLengths = {}
+
+    def _preprocessMatch(self, diamondStr):
+        """
+        Add information from a row of DIAMOND tabular output.
+
+        @param diamondStr: A C{str} with TAB-separated fields from DIAMOND
+            output format 6.
+        @return: A C{tuple} containing the match C{dict} (from
+            C{self._strToDiamondDict}), the protein C{dict} as looked up
+            by C{self._genomesProteins.findProtein}, and the genome C{dict}
+            as looked up by C{self._genomesProteins.findGenome}.
+        """
+        match = self._strToDiamondDict(diamondStr)
+        stitle = match['stitle']
+
+        genomeAccession = self._genomesProteins.genomeAccession(stitle)
+        genome = self._genomesProteins.findGenome(genomeAccession)
+        if genome is None:
+            raise ValueError(
+                'Could not find accession %r in genomes database table.' %
+                genomeAccession)
+
+        proteinAccession = self._genomesProteins.proteinAccession(stitle)
+        protein = self._genomesProteins.findProtein(proteinAccession)
+        if protein is None:
+            raise ValueError(
+                'Could not find accession %r in proteins database table.' %
+                proteinAccession)
+
+        self._referenceLengths[genomeAccession] = genome['length']
+
+        return match, protein, genome
+
+    def _SAMLine(self, match, protein, genome):
+        """
+        Convert DIAMOND match information to a line of SAM file output.
+
+        @param match: A C{dict} with information about the DIAMOND match, as
+            returned by C{DiamondTabularFormat().diamondFieldsToDict} which
+            has been called for us by C{self._preprocessMatch}.
+        @param protein: A C{dict} with information about the protein the
+            DIAMOND was for. The C{dict} is as returned by
+            C{dark.proteins.SqliteIndex.findProtein}.
+        @param genome: A C{dict} with information about the (nucleotide) genome
+            that the protein in the DIAMOND match comes from. The C{dict} is
+            as returned by C{dark.proteins.SqliteIndex.findGenome}.
+        @return: A TAB-separated C{str} line of SAM.
+        """
+        qseqid = (match['qseqid'] if self._keepDescriptions else
+                  match['qseqid'].split(None, 1)[0])
+
+        ranges = GenomeRanges(protein['offsets'])
+
+        # matchStartInGenome = ranges.startInGenome(match)
+        # queryStartInGenome = matchStartInGenome - match['qstart'] - 1
+        # orientations = ranges.orientations()
+
+        print('ranges', ranges)
+        print('protein', protein)
+        print('match', match)
+        # If the query frame is less than zero, the match was with a
+        # translation of the reverse-complemented query. We'll put the
+        # reverse complement into the SAM output. This seems to be standard
+        # / accepted practice, based on my web searches.  See e.g.,
+        # https://www.biostars.org/p/131891/ for what Bowtie2 does and for
+        # some comments on this issue for SAM/BAM files in general.
+        if match['qframe'] > 0:
+            flag = 0
+            qseq = match['full_qseq']
+            qqual = match['full_qqual'] or '*'
+        else:
+            flag = 16
+            qseq = DNARead('id',
+                           match['full_qseq']).reverseComplement().sequence
+            qqual = match['full_qqual'][::-1] if match['full_qqual'] else '*'
+
+        # Make a CIGAR string, including hard-clipped bases at the start and
+        # end of the query (DIAMOND outputs a hard-clipped query sequence).
+        startClipCount = match['qstart'] - 1
+        endClipCount = match['qlen'] - match['qend']
+
+        assert startClipCount >= 0
+        assert endClipCount >= 0, (
+            'Query sequence %s has length %d but the qend value is %d' %
+            (qseq, len(match['qseq']), match['qend']))
+
+        cigar = (
+            ('%dH' % startClipCount if startClipCount else '') +
+            btop2cigar(match['btop'], concise=False, aa=True) +
+            ('%dH' % endClipCount if endClipCount else ''))
+
+        return '\t'.join(map(str, (
+            # 1. QNAME
+            qseqid,
+            # 2. FLAG
+            flag,
+            # 3. RNAME
+            genome['accession'],
+            # 4. POS. This needs to be a 1-based offset into the
+            # nucleotide-equivalent of the DIAMOND subject sequence (which was
+            # a protein since that is how DIAMOND operates). Because DIAMOND
+            # gives back a 1-based protein location, we adjust to 0-based,
+            # multiply by 3 to get to nucleotides, then adjust to 1-based.
+            3 * (match['sstart'] - 1) + 1,
+            # 5. MAPQ
+            self._mappingQuality,
+            # 6. CIGAR
+            cigar,
+            # 7. RNEXT
+            '*',
+            # 8. PNEXT
+            0,
+            # 9. TLEN
+            0,
+            # 10. SEQ
+            qseq,
+            # 11. QUAL
+            qqual,
+            # 12. Alignment score
+            'AS:i:%d' % int(match['bitscore']))))
+
+    def addMatch(self, diamondStr):
+        """
+        Add information from a row of DIAMOND tabular output.
+
+        @param diamondStr: A C{str} with TAB-separated fields from DIAMOND
+            output format 6.
+        """
+        raise NotImplementedError('addMatch must be implemented by a subclass')
+
+    def save(self, fp):
+        """
+        Write SAM output to a file.
+        """
+        raise NotImplementedError('save must be implemented by a subclass')
+
+
+class SimpleDiamondSAMWriter(_DiamondSAMWriter):
+    """
+    Convert DIAMOND tabular output to SAM.
+
+    @param genomesProteins: An open sqlite3 database with information on
+        genomes and proteins, as made by the dark-matter
+        make-protein-database.py script.
+    @param mappingQuality: The mapping quality to use for the SAM MAPQ field
+        (#5). The default (255) indicates that mapping quality information is
+        not available.
+    @param ram: If C{True}, do not use a temporary file to hold the non-header
+        SAM output. This will run faster but use more memory since all
+        non-header SAM output will be stored in RAM and only written out when
+        the full header can be determined.
+    @param keepDescriptions: If C{True}, do not discard text after the first
+        space in query or subject sequence ids. Note that this violates the
+        SAM specification, but since SAM files are TAB-separated there may be
+        only a small chance this will cause problems downstream.
+    """
+    def __init__(self, genomesProteins, mappingQuality=255, ram=False,
+                 keepDescriptions=False):
+        _DiamondSAMWriter.__init__(
+            self, genomesProteins, mappingQuality=mappingQuality, ram=ram,
+            keepDescriptions=keepDescriptions)
+        self._tf = TemporaryFile(mode='w+t', encoding='utf-8')
+        self._writer = SAMWriter(ram=ram)
+
+    def addMatch(self, diamondStr):
+        """
+        Add information from a row of DIAMOND tabular output.
+
+        @param diamondStr: A C{str} with TAB-separated fields from DIAMOND
+            output format 6.
+        """
+        match, protein, genome = self._preprocessMatch(diamondStr)
+        self._writer.addMatchLine(
+            self._SAMLine(match, protein, genome), genome, self._tf)
+
+    def save(self, filename):
+        """
+        Write SAM output to a file.
+
+        @param filename: A C{str} filename or an already-open file handle.
+        """
+        with asHandle(filename) as fp:
+            self._writer.save(fp, tmpfp=self._tf)
+        self._tf.close()
+        self._tf = self._writer = None
+
+
+class PerReferenceDiamondSAMWriter(_DiamondSAMWriter):
+    """
+    Convert DIAMOND tabular output to SAM and write it to per-reference
+    output files.
+
+    @param genomesProteins: An open sqlite3 database with information on
+        genomes and proteins, as made by the dark-matter
+        make-protein-database.py script.
+    @param baseFilenameFunc: A function that produces a base filename to save
+        SAM output to for a given reference name. A '.sam' suffix will be
+        added to these files.
+    @param mappingQuality: The mapping quality to use for the SAM MAPQ field
+        (#5). The default (255) indicates that mapping quality information is
+        not available.
+    @param ram: If C{True}, do not use a temporary file to hold the non-header
+        SAM output. This will run faster but use more memory since all
+        non-header SAM output will be stored in RAM and only written out when
+        the full header can be determined.
+    @param keepDescriptions: If C{True}, do not discard text after the first
+        space in query or subject sequence ids. Note that this violates the
+        SAM specification, but since SAM files are TAB-separated there may be
+        only a small chance this will cause problems downstream.
+    @param fpcMaxsize: An C{int} maximum size for the file pointer cache used
+        because SAM output may be being written to an arbitrary number of
+        output files so we need to make sure the OS limit on open files is not
+        exceeded and that we are not inefficiently constantly re-opening and
+        writing single lines to files). If C{None} or C{0} is passed, a size
+        of half the maximum number of files will be used.
+    """
+    SAM_SUFFIX = '.sam'
+    TMP_SUFFIX = '.tmp'
+
+    def __init__(self, genomesProteins, baseFilenameFunc=None,
+                 mappingQuality=255, ram=False, keepDescriptions=False,
+                 fpcMaxsize=0):
+        _DiamondSAMWriter.__init__(
+            self, genomesProteins, mappingQuality=mappingQuality, ram=ram,
+            keepDescriptions=keepDescriptions)
+
+        self._writers = {}
+        self._baseFilenameFunc = baseFilenameFunc
+        self._fpc = FilePointerCache(
+            maxsize=(fpcMaxsize or getrlimit(RLIMIT_NOFILE)[0] >> 1),
+            openArgs={'mode': 'wt'}, reopenArgs={'mode': 'at'})
+
+    def addMatch(self, diamondStr):
+        """
+        Add information from a row of DIAMOND tabular output.
+
+        @param diamondStr: A C{str} with TAB-separated fields from DIAMOND
+            output format 6.
+        """
+        match, protein, genome = self._preprocessMatch(diamondStr)
+        genomeAccession = genome['accession']
+        basename = self._baseFilenameFunc(match['stitle'])
+
+        try:
+            writer, _ = self._writers[genomeAccession]
+        except KeyError:
+            writer = SAMWriter()
+            self._writers[genomeAccession] = (writer, basename)
+
+        writer.addMatchLine(self._SAMLine(match, protein, genome), genome,
+                            self._fpc.open(basename + self.TMP_SUFFIX))
+
+    def save(self):
+        """
+        Write SAM output to all files.
+        """
+        # Close the open file descriptor cache.
+        self._fpc.close()
+
+        # Write out SAM (including headers) for each reference.
+        for writer, basename in self._writers.values():
+            tmp = basename + self.TMP_SUFFIX
+            with open(basename + self.SAM_SUFFIX, 'wt') as samfp:
+                with open(tmp) as tmpfp:
+                    writer.save(samfp, tmpfp)
+            unlink(tmp)
diff --git a/dark/errors.py b/dark/errors.py
new file mode 100644
index 0000000..6cc332f
--- /dev/null
+++ b/dark/errors.py
@@ -0,0 +1,2 @@
+class DatabaseDuplicationError(Exception):
+    """An attempt was made to insert a duplicate into a database."""
diff --git a/dark/fpcache.py b/dark/fpcache.py
new file mode 100644
index 0000000..1611afb
--- /dev/null
+++ b/dark/fpcache.py
@@ -0,0 +1,70 @@
+from os.path import exists
+from threading import RLock
+from cachetools import LFUCache, cached
+
+
+class _FilePointerCache(LFUCache):
+    """
+    A least-frequently-used file pointer cache class with an eviction method
+    that closes open files.
+    """
+    def popitem(self):
+        """
+        Evict a cache item.
+
+        @return: The C{str} name of the evicted file.
+        """
+        _, (filename, fp) = super().popitem()
+        fp.close()
+        return filename
+
+    def close(self):
+        """
+        Close all open files in the cache.
+        """
+        while self.currsize:
+            self.popitem()
+
+
+class FilePointerCache(object):
+    """
+    A file pointer cache class for simultaneously opening multiple files,
+    without exceeeding the per-process operating system limit on open file
+    handles.
+
+    @param maxsize: The C{int} maximum size of the cache.
+    @param openArgs: If not C{None}, a C{dict} of keyword arguments to pass to
+        open when opening a new file. If C{None}, the file will be opened
+        with the default mode ('rt').
+    @param reopenArgs: If not C{None}, a C{dict} of keyword arguments to pass
+        to open when opening an already existing file. If C{None}, the file
+        will be opened with the default mode ('rt').
+    """
+    def __init__(self, maxsize=32, openArgs=None, reopenArgs=None):
+        self._openArgs = openArgs or {}
+        self._reopenArgs = reopenArgs or {}
+        self._filePointerCache = _FilePointerCache(maxsize=maxsize)
+        self._lock = RLock()
+
+        @cached(self._filePointerCache, lock=self._lock)
+        def _open(filename):
+            if exists(filename):
+                fp = open(filename, **self._reopenArgs)
+            else:
+                fp = open(filename, **self._openArgs)
+            return filename, fp
+
+        self._open = _open
+
+    def open(self, filename):
+        return self._open(filename)[1]
+
+    def close(self):
+        with self._lock:
+            self._filePointerCache.close()
+
+    def __enter__(self):
+        return self
+
+    def __exit__(self, exc_type, exc_value, traceback):
+        self.close()
diff --git a/dark/genbank.py b/dark/genbank.py
new file mode 100644
index 0000000..dfcfbe6
--- /dev/null
+++ b/dark/genbank.py
@@ -0,0 +1,221 @@
+import warnings
+
+
+class GenomeRanges(object):
+    """
+    Split and manage a GenBank range string as converted by BioPython.
+
+    @param rangeStr: A C{str} indicating the (0-based) genome nucleotide
+        offsets covered by a protein. This can have the following example
+        forms:
+
+        [9462:10137](+)
+        [11969:12575](-)
+        join{[126386:126881](-), [125941:126232](+)}
+        join{[153005:153821](+), [82030:82618](-), [80480:81305](-)}
+
+        Note that this is the string as returned from BioPython SeqIO.parse
+        when it parses a GenBank flat file, not the string that is in such a
+        file. For that use the C{splitRange} function you can find in an
+        earlier version of this file (i.e., look in its git history).
+
+    @raise ValueError: If C{rangeStr} cannot be correctly parsed.
+    @return: A C{list} whose elements are 3-C{tuple}s holding the C{int} start
+        and stop offsets and a C{bool} that indicates whether the offsets
+        corresponds to the forward strand (i.e., will be C{True} when
+        there is a '(+)' in a range and C{False} when there is a '(-)' which
+        indicates the reverse strand).
+
+        Example arguments and their return values:
+
+            '[9462:10137](+)' => ((9462, 10137, True),)
+
+            '[11969:12575](-)' => ((11969, 12575, False),)
+
+            'join{[126386:126881](-), [125941:126232](+)}' =>
+                ((126386, 126881, False),
+                 (125941, 126232, True))
+    """
+    def __init__(self, rangeStr):
+        if rangeStr.startswith('join{') and rangeStr[-1] == '}':
+            join = True
+            inner = rangeStr[5:-1]
+        else:
+            join = False
+            inner = rangeStr
+
+        ranges = []
+        subRanges = inner.split(', ')
+        nRanges = len(subRanges)
+
+        if nRanges == 1 and join:
+            raise ValueError('Could not parse GenBank range string "%s". '
+                             'join{} can only be used with multiple ranges.' %
+                             rangeStr)
+        elif nRanges > 1 and not join:
+            raise ValueError('Could not parse GenBank range string "%s". '
+                             'Multiple ranges must be wrapped in join{}.' %
+                             rangeStr)
+
+        for subRange in subRanges:
+            if subRange.endswith('](+)'):
+                forward = True
+            elif subRange.endswith('](-)'):
+                forward = False
+            else:
+                raise ValueError('Could not parse GenBank range string "%s". '
+                                 'Range "%s" does not end with ](+) or ](-).' %
+                                 (rangeStr, subRange))
+
+            if not subRange.startswith('['):
+                raise ValueError('Could not parse GenBank range string "%s". '
+                                 'Range "%s" does not start with "[".' %
+                                 (rangeStr, subRange))
+            try:
+                start, stop = subRange[1:-4].split(':')
+            except ValueError as e:
+                raise ValueError('Could not parse GenBank range string "%s". '
+                                 'Original parsing ValueError was "%s".' %
+                                 (rangeStr, e))
+            else:
+                if start.startswith('<'):
+                    start = start[1:]
+                if stop.startswith('>'):
+                    stop = stop[1:]
+                start, stop = map(int, (start, stop))
+                if start > stop:
+                    raise ValueError(
+                        'Could not parse GenBank range string "%s". '
+                        'Offset values (%d, %d) cannot decrease.' %
+                        (rangeStr, start, stop))
+
+                ranges.append((start, stop, forward))
+
+        self.ranges = self._mergeContiguousRanges(ranges)
+        self._nRanges = len(self.ranges)
+
+    def _mergeContiguousRanges(self, ranges):
+        """
+        Merge ranges that are contiguous (follow each other immediately on the
+        genome).
+
+        @param ranges: An iterable of (start, stop, forward) tuples.
+        @return: A C{tuple} of (start, stop, forward) tuples with contiguous
+            ranges found in C{ranges} merged.
+        """
+        result = []
+        lastStart = lastStop = lastForward = None
+
+        for index, (start, stop, forward) in enumerate(ranges):
+            if lastStart is None:
+                lastStart, lastStop, lastForward = start, stop, forward
+            else:
+                if start == lastStop and forward == lastForward:
+                    # This range continues the previous one.
+                    warnings.warn(
+                        'Contiguous GenBank ranges detected: [%d:%d] '
+                        'followed by [%d:%d].' %
+                        (lastStart, lastStop, start, stop))
+                    lastStop = stop
+                else:
+                    # Emit the range that just got terminated.
+                    result.append((lastStart, lastStop, lastForward))
+                    # And remember this one.
+                    lastStart, lastStop, lastForward = start, stop, forward
+
+        # Emit the final range.
+        result.append((lastStart, lastStop, lastForward))
+
+        return tuple(result)
+
+    def __str__(self):
+        return '<%s: %s>' % (
+            self.__class__.__name__, ', '.join(map(str, self.ranges)))
+
+    def circular(self, genomeLength):
+        """
+        Determine whether the offset ranges of a protein in a genome span the
+        end of the genome (indicating that the genome may be circular).
+
+        @param genomeLength: The C{int} length of the genome.
+        @return: A C{bool}, C{True} if the ranges overlaps the genome end.
+        """
+        if self._nRanges == 1:
+            # If there is only one range, we simply return False even though it
+            # could be that the single range spans the whole genome. That is
+            # not the same as having two ranges, one of which ends at the end
+            # of the genome and the next starting at the start of the genome.
+            #
+            # I.e., a range tuples like this ((0, 100, False),) on a genome of
+            # length 100 will return False, whereas ((0, 75, False), (75, 100,
+            # False)) on a genome of length 100 will return True.
+            #
+            # The decision about whether you consider the degenerate case of
+            # ((0, 100, False),) to indicate circularity is a bit
+            # arbitrary. Here I've decided that it should not because I intend
+            # to use this function to decide whether to output multiple SAM
+            # lines for the match of a read against a reference. In this
+            # degenerate case just a single SAM line suffices, whereas in a
+            # non-degenerate case (such as ((0, 75, False), (75, 100, False))
+            # on a genome of length 100), a chimeric SAM alignment must be
+            # used, which requires quite different processing than a normal
+            # linear match. YMMV, and if you don't like it you can easily add
+            # an argument to this function that changes this behaviour (leaving
+            # the default to do what the code currently does).
+            return False
+
+        for index, (start, stop, _) in enumerate(self.ranges):
+            if stop == genomeLength:
+                return self.ranges[(index + 1) % self._nRanges][0] == 0
+
+        return False
+
+    def startInGenome(self, match):
+        """
+        Calculate the position in the nucleotide genome where a DIAMOND match
+        begins.
+
+        @param match: A C{dict} with information about the DIAMOND match, as
+            returned by C{self._preprocessMatch}. Must contain a 'sstart' key
+            giving a 1-based C{int} offset of the start of the match in the
+            protein.
+        @return: The C{int} offset of the start of the DIAMOND match in the
+            genome.
+        """
+
+        # Calculate the start of the match in the genome, given its start in
+        # the protein.
+        offsetInGenome = remaining = (match['sstart'] - 1) * 3
+
+        for start, stop, _ in self.ranges:
+            rangeWidth = stop - start
+            if remaining < rangeWidth:
+                # The match starts in this range.
+                return start + remaining
+            else:
+                remaining -= rangeWidth
+        else:
+            raise ValueError(
+                'Starting nucleotide offset %d not found in protein '
+                'nucleotide ranges %s.' %
+                (offsetInGenome,
+                 ', '.join(('(%d, %d)' % (i, j))
+                           for i, j, _ in self.ranges)))
+
+    def orientations(self):
+        """
+        Produce the set of all orientations for our ranges.
+
+        @return: A C{set} of C{True} and C{False} range orientations.
+        """
+        return set(r[2] for r in self.ranges)
+
+    def distinctRangeCount(self, genomeLength):
+        """
+        Determine the number of distinct ranges, given the genome length.
+
+        @param genomeLength: The C{int} length of the genome.
+        @return: An C{int} which is the number of ranges that are not
+            contiguous with one another.
+        """
+        return len(self.ranges) - int(self.circular(genomeLength))
diff --git a/dark/process.py b/dark/process.py
index f94620e..f21aa8b 100644
--- a/dark/process.py
+++ b/dark/process.py
@@ -29,6 +29,8 @@ class Executor(object):
             name), in which case the shell is not used.
         @return: A C{CompletedProcess} instance. This has attributes such as
             C{returncode}, C{stdout}, and C{stderr}. See pydoc subprocess.
+            If C{dryRun} is passed as C{True} to C{__init__}, C{None} is
+            returned.
         """
         if isinstance(command, six.string_types):
             # Can't have newlines in a command given to the shell.
diff --git a/dark/proteins.py b/dark/proteins.py
index d6cefee..ddebdd3 100644
--- a/dark/proteins.py
+++ b/dark/proteins.py
@@ -1,4 +1,4 @@
-from __future__ import division, print_function
+from __future__ import print_function, division
 
 import os
 from collections import defaultdict, Counter
@@ -407,8 +407,8 @@ class ProteinGrouper(object):
         Produce a string representation of the pathogen summary.
 
         @param title: The C{str} title for the output.
-        @param preamble: The C{str} descriptive preamble for the HTML page, or
-            C{None} if no preamble is needed.
+        @param preamble: The C{str} descriptive preamble, or C{None} if no
+            preamble is needed.
         @return: A C{str} suitable for printing.
         """
         # Note that the string representation contains much less
diff --git a/dark/taxonomy.py b/dark/taxonomy.py
index a0c27c7..52ff2d3 100644
--- a/dark/taxonomy.py
+++ b/dark/taxonomy.py
@@ -1,3 +1,8 @@
+import sqlite3
+from six import string_types
+from operator import attrgetter
+from cachetools import LRUCache, cachedmethod
+
 from dark.database import getDatabaseConnection
 
 
@@ -31,7 +36,7 @@ class LineageFetcher(object):
 
         lineage = []
         gi = int(title.split('|')[1])
-        query = 'SELECT taxID from gi_taxid where gi = %d' % gi
+        query = 'SELECT taxID FROM gi_taxid WHERE gi = %d' % gi
 
         try:
             while True:
@@ -39,12 +44,12 @@ class LineageFetcher(object):
                 taxID = self._cursor.fetchone()[0]
                 if taxID == 1:
                     break
-                query = 'SELECT name from names where taxId = %s' % taxID
+                query = 'SELECT name FROM names WHERE taxId = %s' % taxID
                 self._cursor.execute(query)
                 scientificName = self._cursor.fetchone()[0]
                 lineage.append((taxID, scientificName))
                 # Move up to the parent.
-                query = ('SELECT parent_taxID from nodes where taxID = %s' %
+                query = ('SELECT parent_taxID FROM nodes WHERE taxID = %s' %
                          taxID)
         except TypeError:
             lineage = []
@@ -60,3 +65,158 @@ class LineageFetcher(object):
         self._cursor.close()
         self._db.close()
         self._cursor = self._db = self._cache = None
+
+
+class AccessionLineageFetcher(object):
+    """
+    Provide access to the NCBI taxonomy database so we can retrieve the
+    taxonomy lineage corresponding to an accession number.
+
+    @param dbFilenameOrConnection: Either a C{str} database filename or an
+        open sqlite3 database. The database must contain taxonomy information
+        with tables and columns named as in the building scripts used in
+        https://github.com/acorg/ncbi-taxonomy-database
+    """
+    CACHE_SIZE = 10E6
+
+    def __init__(self, dbFilenameOrConnection):
+        if isinstance(dbFilenameOrConnection, string_types):
+            self._db = sqlite3.connect(dbFilenameOrConnection)
+            self._closeConnection = True
+        else:
+            self._db = dbFilenameOrConnection
+            self._closeConnection = False
+        self._cache = LRUCache(maxsize=self.CACHE_SIZE)
+
+    @cachedmethod(attrgetter('_cache'))
+    def lineage(self, accession):
+        """
+        Get lineage information from the taxonomy database for a given title.
+
+        @param accession: A C{str} accession number. This must of coures
+            include the version (e.g., the ".1" part of "DQ011818.1") if the
+            taxonomy database in use uses it (as is the case with the database
+            built by https://github.com/acorg/ncbi-taxonomy-database
+        @raise ValueError: If a taxonomy id cannot be found in the names or
+            nodes table.
+        @raise AttributeError: If C{self.close} has been called.
+        @return: A C{tuple} of the taxonomic categories of the title. Each
+            tuple element is a 3-tuple of (C{int}, C{str}, C{str}) giving a
+            taxonomy id a (scientific) name, and the rank (species, genus,
+            etc). The first element in the list corresponds to the passed
+            C{accession}, number and each successive element is the parent of
+            the preceeding one. The top level of the taxonomy hierarchy
+            (with taxid = 1) is not included in the returned list. If no
+            taxonomy information is found for C{accession}, return C{None}.
+        """
+        cursor = self._db.cursor()
+        execute = cursor.execute
+        fetchone = cursor.fetchone
+
+        execute('SELECT taxid FROM accession_taxid WHERE accession = ?',
+                (accession,))
+        row = fetchone()
+        if row:
+            taxid = int(row[0])
+            lineage = []
+
+            while taxid != 1:
+                execute('SELECT name FROM names WHERE taxid = ?', (taxid,))
+                result = fetchone()
+                if result is None:
+                    raise ValueError(
+                        'Could not find taxonomy id %r in names table' %
+                        (taxid,))
+                name = result[0]
+                execute('SELECT parent_taxid, rank FROM nodes WHERE taxid = ?',
+                        (taxid,))
+                result = fetchone()
+                if result is None:
+                    raise ValueError(
+                        'Could not find taxonomy id %r in nodes table' %
+                        (taxid,))
+                parentTaxid, rank = result
+                lineage.append((taxid, name, rank))
+                taxid = parentTaxid
+
+            return tuple(lineage) or None
+
+    def close(self):
+        """
+        Close the database connection (if we opened it).
+        """
+        if self._closeConnection:
+            self._db.close()
+        # Set self._db to None so self.lineage will raise an AttributeError
+        # exception if it is called again.
+        self._db = None
+
+
+def isRetrovirus(lineage):
+    """
+    Determine whether a lineage corresponds to a retrovirus.
+
+    @param lineage: A C{tuple} of taxonomy id, scientific name, and rank
+        as returned by C{AccessionLineageFetcher.lineage}.
+    @return: C{True} if the lineage corresponds to a retrovirus, C{False}
+        otherwise.
+    """
+    for taxid, name, rank in lineage:
+        if rank == 'family' and name == 'Retroviridae':
+            return True
+
+    return False
+
+
+def isRNAVirus(lineage):
+    """
+    Determine whether a lineage corresponds to an RNA virus.
+
+    @param lineage: A C{tuple} of taxonomy id, scientific name, and rank
+        as returned by C{AccessionLineageFetcher.lineage}.
+    @return: C{True} if the lineage corresponds to an RNA virus, C{False}
+        otherwise.
+    """
+    for taxid, name, rank in lineage:
+        if rank == 'realm' and name == 'Riboviria':
+            return True
+
+    return False
+
+
+def formatLineage(lineage, namesOnly=False, separator=None, prefix=''):
+    """
+    Format a lineage for printing.
+
+    @param lineage: A C{tuple} of taxonomy id, scientific name, and rank
+        as returned by C{AccessionLineageFetcher.lineage}.
+    @param namesOnly: If C{True} only print taxonomic names.
+    @param separator: A C{str} separator to put between fields. If C{None},
+        return a space-padded aligned columns.
+    @param prefix: A C{str} to put at the start of each line.
+    @return: A formatted C{str} for printing.
+    """
+    taxids, names, ranks = [], [], []
+
+    for (taxid, name, rank) in lineage:
+        taxids.append(str(taxid))
+        ranks.append('-' if rank == 'no rank' else rank)
+        names.append(name)
+
+    if namesOnly:
+        # The separator is guaranteed to be set by our caller.
+        return prefix + separator.join(names)
+
+    if separator is not None:
+        return '\n'.join(
+            '%s%s%s%s%s%s' % (prefix, rank, separator, name, separator, taxid)
+            for (taxid, name, rank) in zip(taxids, names, ranks))
+    else:
+        taxidWidth = max(len(x) for x in taxids)
+        nameWidth = max(len(x) for x in names)
+        rankWidth = max(len(x) for x in ranks)
+
+        return '\n'.join(
+            '%s%-*s %-*s %*s' % (
+                prefix, rankWidth, rank, nameWidth, name, taxidWidth, taxid)
+            for (taxid, name, rank) in zip(taxids, names, ranks))
diff --git a/dark/utils.py b/dark/utils.py
index 810088f..0059330 100644
--- a/dark/utils.py
+++ b/dark/utils.py
@@ -1,5 +1,6 @@
 from __future__ import division
 
+import os
 import string
 import six
 import bz2
@@ -103,27 +104,32 @@ def median(l):
 
 
 @contextmanager
-def asHandle(fileNameOrHandle, mode='r'):
+def asHandle(fileNameOrHandle, mode='rt', encoding='UTF-8'):
     """
-    Decorator for file opening that makes it easy to open compressed files.
+    Decorator for file opening that makes it easy to open compressed files
+    and which can be passed an already-open file handle or a file name.
     Based on L{Bio.File.as_handle}.
 
     @param fileNameOrHandle: Either a C{str} or a file handle.
+    @param mode: The C{str} mode to use for opening the file.
+    @param encoding: The C{str} encoding to use when opening the file.
     @return: A generator that can be turned into a context manager via
         L{contextlib.contextmanager}.
     """
     if isinstance(fileNameOrHandle, six.string_types):
         if fileNameOrHandle.endswith('.gz'):
             if six.PY3:
-                yield gzip.open(fileNameOrHandle, mode='rt', encoding='UTF-8')
+                yield gzip.open(fileNameOrHandle, mode=mode, encoding=encoding)
             else:
-                yield gzip.GzipFile(fileNameOrHandle)
+                yield gzip.GzipFile(fileNameOrHandle, mode=mode)
         elif fileNameOrHandle.endswith('.bz2'):
             if six.PY3:
-                yield bz2.open(fileNameOrHandle, mode='rt', encoding='UTF-8')
+                yield bz2.open(fileNameOrHandle, mode=mode, encoding=encoding)
             else:
-                yield bz2.BZ2File(fileNameOrHandle)
+                yield bz2.BZ2File(fileNameOrHandle, mode=mode)
         else:
+            # Putting mode=mode, encoding=encoding into the following
+            # causes a hard-to-understand error from the mocking library.
             with open(fileNameOrHandle) as fp:
                 yield fp
     else:
@@ -251,3 +257,18 @@ def countPrint(mesg, count, len1, len2=None):
                         mesg,
                         count, len1, percentage(count, len1),
                         count, len2, percentage(count, len2)))
+
+
+@contextmanager
+def cd(newdir):
+    """
+    Trivial context manager for temporarily switching directory.
+
+    @param dir: A C{str} directory to cd to.
+    """
+    olddir = os.getcwd()
+    try:
+        os.chdir(newdir)
+        yield newdir
+    finally:
+        os.chdir(olddir)
diff --git a/doc/protein-database.md b/doc/protein-database.md
new file mode 100644
index 0000000..cbea76e
--- /dev/null
+++ b/doc/protein-database.md
@@ -0,0 +1,147 @@
+# THIS IS ALL OUT OF DATE! IGNORE.
+
+# Making a protein database from the RVDB nucleotide database
+
+There is a `Makefile` in the `misc` directory that can perform most of the
+steps described below for you. To use it, try something like:
+
+```sh
+$ mkdir ~/protein-database
+$ cd ~/protein-database
+$ ln -s ~/dark-matter/misc/Makefile-protein-database Makefile
+$ make download
+$ make
+```
+
+and cross your fingers.  Else, read on and do these things one at a time.
+
+## Executables
+
+There are scripts run below called `download-genbank.sh` and
+`make-protein-database.py`. They live in the
+[dark-matter](https://github.com/acorg/dark-matter) `bin` directory, which
+should be in your `PATH`. Try `pip install dark-matter` if you don't
+already have it installed.
+
+## Get the clustered RVDB nucleotide database
+
+Get `C-RVDBv15.1.fasta` (or a later version) from
+[here](https://hive.biochemistry.gwu.edu/rvdb).
+
+`C-RVDBv15.1.fasta` contains 699,223 viral nucleotide genomes. We want to
+make a database of the proteins in the sequences (actually, the coding
+regions).
+
+This is an alternative to using the protein version of the RVDB (available
+[here](https://rvdb-prot.pasteur.fr/)) which doesn't do exactly what we
+need (we need more intermediate information, such as exactly where the
+protein sequences fall in the nucleotide genomes).
+
+## Download the GenBank nucleotide to accession to taxonomy id file
+
+You will need `nucl_gb.accession2taxid.gz` which you can download from
+[here](ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/accession2taxid/nucl_gb.accession2taxid.gz)
+(there's no need to uncompress it).
+
+### Download the GenBank info for all nucleotide genomes
+
+The following uses GNU `parallel` to download files containing 1000
+[GenBank flat-file](https://www.ncbi.nlm.nih.gov/Sitemap/samplerecord.html)
+format records for the accession numbers of the sequences in
+`C-RVDBv15.1.fasta`:
+
+```sh
+$ download-genbank.sh ids C-RVDBv15.1.fasta accession-to-name
+```
+
+If you have an NCBI API key, put its value into your environment in
+`NCBI_API_KEY`. If that variable is set the script will send 10 requests
+per second to NCBI (instead of 3). See
+[here](https://www.ncbi.nlm.nih.gov/books/NBK25497/) for how to obtain an
+NCBI API key.
+
+
+This creates a directory called `ids` (specified by the 3rd argument,
+above), holding files with names like
+
+```
+$ ls -1 ids/ids_aaaa*
+ids/ids_aaaaa.genbank
+ids/ids_aaaab.genbank
+ids/ids_aaaac.genbank
+ids/ids_aaaad.genbank
+ids/ids_aaaae.genbank
+```
+
+### Accession numbers to names
+
+The above also creates a file called `accession-to-name` (specified by its
+3rd argument) that maps accession numbers to nucleotide sequence names (as
+in `C-RVDBv15.1.fasta`). That file will be needed below when we make the
+protein file for DIAMOND and put the nucleotide sequence name into square
+brackets at the end of the line.
+
+## Make the protein FASTA and database
+
+Run
+
+```sh
+$ make-protein-database.py --gb ids/ids_*.genbank \
+      --databaseFile genomes-proteins.db \
+      --nucleotideAccessionToTaxidFile nucl_gb.accession2taxid.gz \
+      --accessionToNameFile accession-to-name \
+      --force > proteins.fasta
+```
+
+This creates `proteins.fasta` which contains an amino acid sequence for all
+coding regions in all the nucleotide sequences in `C-RVDBv15.1.fasta`, as
+extracted from the `ids/ids_*.genbank` files downloaded above.
+
+The `proteins.fasta` file can be given to `DIAMOND` to make its database.
+
+It also writes a `genomes-proteins.db` file. This is an sqlite3 database with
+tables as follows:
+
+```
+CREATE TABLE proteins (
+    accession VARCHAR UNIQUE PRIMARY KEY,
+    sequence VARCHAR NOT NULL,
+    offsets VARCHAR NOT NULL,
+    reversed INTEGER,
+    gene VARCHAR,
+    note VARCHAR
+);
+
+CREATE TABLE genomes (
+    accession VARCHAR UNIQUE PRIMARY KEY,
+    length INTEGER,
+    proteinCount INTEGER,
+    taxid INTEGER
+);
+```
+
+The former holds information about the protein sequences that will be used
+later, principally for determining coverage of the nucleotide sequence the
+proteins came from.  The latter contains the number of proteins (coding
+regions, actually) for a given nucleotide genome.
+
+## Make a nucleotide sequence id to taxonomy id database
+
+Later, when we match (using DIAMOND) protein sequences, we will be able to
+extract their corresponding nucleotide genome accession number from their
+sequence id (it's the 2nd accession number in a name like
+`acc|GENBANK|YP_009137150.1|GENBANK|NC_001798|neurovirulence protein`).
+We'll want to look up the corresponding NCBI taxonomy id, which we can do
+either by the protein accession number or the nucleotide one (since they're
+obviously from the same organism). Because there are fewer nucleotide
+sequences we here make a small database mapping nucleotide accession number
+to taxid.
+
+So:
+
+```sh
+$ ./add-accession-taxids.py
+```
+
+which adds taxonomy ids to all the accession numbers in the `genomes`
+database table above. This is used in later processing and filtering.
diff --git a/requirements.txt b/requirements.txt
index 53ecf5d..125aa85 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,5 +1,6 @@
 biopython>=1.71
 bz2file>=0.98
+cachetools>=3.1.0
 ete3>=3.1.1
 flake8>=3.5.0
 ipython>=3.1.0
@@ -7,7 +8,7 @@ matplotlib>=1.4.3
 mock>=1.0.1; python_version == '2.7'
 mysql-connector-python==8.0.11
 numpy>=1.14.2
-pyfaidx==0.5.0
+pyfaidx>=0.5.0
 pysam>=0.15.2
 pytest>=4.0.2
 pyzmq>=14.3.1
diff --git a/setup.py b/setup.py
index 3e65832..bf4d48a 100644
--- a/setup.py
+++ b/setup.py
@@ -24,8 +24,11 @@ def version():
 # bin files that are not (yet) part of the distribution.
 scripts = [
     'bin/aa-info.py',
+    'bin/aa-to-dna.py',
     'bin/aa-to-properties.py',
+    'bin/accession-to-taxonomy.py',
     'bin/adaptor-distances.py',
+    'bin/alignment-panel-civ.py',
     'bin/alignments-per-read.py',
     'bin/bit-score-to-e-value.py',
     'bin/cat-json-blast-records.py',
@@ -40,6 +43,7 @@ scripts = [
     'bin/create-newick-relabeling-output.py',
     'bin/dark-matter-version.py',
     'bin/dna-to-aa.py',
+    'bin/download-genbank.sh',
     'bin/e-value-to-bit-score.py',
     'bin/extract-ORFs.py',
     'bin/fasta-base-indices.py',
@@ -67,15 +71,18 @@ scripts = [
     'bin/graph-evalues.py',
     'bin/local-align.py',
     'bin/make-fasta-database.py',
+    'bin/make-protein-database.py',
     'bin/ncbi-fetch-id.py',
     'bin/newick-to-ascii.py',
     'bin/noninteractive-alignment-panel.py',
+    'bin/parse-genbank-flat-file.py',
     'bin/position-summary.py',
     'bin/pre-commit.sh',
     'bin/print-blast-xml-for-derek.py',
     'bin/print-blast-xml.py',
     'bin/print-read-lengths.py',
     'bin/proteins-to-pathogens.py',
+    'bin/proteins-to-pathogens-civ.py',
     'bin/randomize-fasta.py',
     'bin/read-blast-json.py',
     'bin/read-blast-xml.py',
@@ -86,6 +93,7 @@ scripts = [
     'bin/sam-references.py',
     'bin/sff-to-fastq.py',
     'bin/split-fasta-by-adaptors.py',
+    'bin/subset-protein-database.py',
     'bin/summarize-fasta-bases.py',
     'bin/summarize-reads.py',
     'bin/trim-primers.py',
@@ -98,7 +106,7 @@ setup(name='dark-matter',
       packages=['dark', 'dark.blast', 'dark.diamond'],
       url='https://github.com/acorg/dark-matter',
       download_url='https://github.com/acorg/dark-matter',
-      author='Terry Jones, Barbara Muehlemann, Sophie Mathias, Tali Veith',
+      author='Terry Jones, Barbara Muehlemann, Tali Veith, Sophie Mathias',
       author_email='tcj25@cam.ac.uk',
       keywords=['virus discovery'],
       classifiers=[
@@ -123,6 +131,7 @@ setup(name='dark-matter',
           'pyfaidx>=0.4.8.4',
           'pyzmq>=14.3.1',
           'requests>=2.18.4',
+          'cachetools>=3.1.0',
           'simplejson>=3.5.3',
           'six>=1.11.0',
       ])
diff --git a/tox.ini b/tox.ini
index b47fa7a..b038fc1 100644
--- a/tox.ini
+++ b/tox.ini
@@ -1,24 +1,29 @@
 [tox]
 # Adding pypy to the envlist doesn't work yet.
-envlist = py27,py35
+envlist = py27,py35,py36,py37
 skip_missing_interpreters = True
 
 [testenv]
 commands =
   {envbindir}/python setup.py develop
-  python -m discover -v
+  pytest
 
 [testenv:py27]
 deps =
-    discover==0.4.0
     -rrequirements.txt
 
 [testenv:py35]
 deps =
-    discover==0.4.0
+    -rrequirements.txt
+
+[testenv:py36]
+deps =
+    -rrequirements.txt
+
+[testenv:py37]
+deps =
     -rrequirements.txt
 
 [testenv:pypy]
 deps =
-    discover==0.4.0
     -rrequirements-pypy.txt


Test patch:
diff --git a/test/diamond/sample_proteins.py b/test/diamond/sample_proteins.py
new file mode 100644
index 0000000..e46fff8
--- /dev/null
+++ b/test/diamond/sample_proteins.py
@@ -0,0 +1,93 @@
+SAMPLE_DATA = {
+    'genomes':
+    {
+        'NC_030446.1':
+        {
+            'id': ('NC_030446.1 Tibetan frog hepatitis B virus isolate '
+                   '243398, complete genome'),
+            'genome': (
+                'CAGTAGTTATTACGGAAACTATTGTTGTTTGATTATCATGTACTGCCTGTAATTCCATG'
+                'GATCCTATTGTTAGATTTTCTTTGCAATTAGTCGAAAACTTACCTGATGATTATTTTCC'
+                'TGACAATGATCAATTTTTGAGGCATGCATTTGATGTGGCACACCTGTATTGGGAAGATG'
+                'CTGAAAATAAAACTTACCATGTTCAAACTTCTCGTTCAATCTATGAAGTTATCCTCATC'
+                'TTAAGAACTGTCCCCCCTTTGTTGGGATTCACTCACAGACAGTTGAGGGAACATTGGGA'
+                'ATACCTCAGAGATCAGGAAAATTGGGGAAACCAAACACCTCCTCCATATCCCCTTCCTG'
+                'AGTTGCCTCATTCTTTCCCTGTGTCCATTAGGGACCCACCATTGCTTACTCAGTTCCAA'
+                'GCTATCACTGATAGACACCGAGAGCTGAGAGGTGCAGTAAATGTGGGTCAACCTAATTC'
+                'TTACAGGGAACTGTGGATGCTCTCAGAATTCAACAGATTCTTGGGGATTATGTTAGAGT'
+                'CACAAACACGGTTTCTAATAAAAAGACTCTGGTATCACATAAACTGTATTGCCTGGGGA'
+                'GAAGCCACTGTCCACAACCACGTAGCTAAGGTAAGAACTTGGTTGGCTACGCCAGCCCC'
+                'ATACAGGGGTCAAAACGCTCCCCTTATTGAAGAGATTGTGGAGCCGCAACGTCAGAGGA'
+                'CGTCCAGGCCCAATTCTAGAAGAGCTAGAGGGGGCTCTCACCCAGGCGGGGACAGAGGA'
+                'CTATATTATCAGAGGACAGTCAGGAGATTTAATAGTCCAAGGGAAAGAGGATCTAGACG'
+                'GAGGAGTCCTTCCCCTGCAGGCAGAACCTGAACGCCTGTTAGGACCATTGAAGGGACTC'
+                'TATCAAAATGACGTTCCCCCCTTCAATCCGTCTTGGGAGATCCCTACCATTCAAGGGGT'
+                'ACACCTTCAAGAAGATATTATCAGAGAGATTCCCGATAGGAGATGGGCTTACATTTACC'
+                'CAATTAGGTTTTGGCCAAACGAGGTTAAATACCGCCCCCGGGAACAAGGAATCACAGAT'
+                'AAATACCCCGAACATACCAGAGCACATTCAGAAAAAGTCTCCCTCTACTTGCAGACATT'
+                'ATATAACAATGGGATCTTATACAAAAGGGAATCCAGGCATGTGCTTTTATTCAATGATA'
+                'ATCCTTACAGCTGGGAGTCTAATAATCTACTGCCCTGCCACTGCCATAAAAGCCAGGAC'
+                'TACGAGGGCCTCGGAAGGCAGCCACCGAGCAGACGGAGCAATGAGGTTCTTCGCCCGAG'
+                'TGTGGAGTTGGATAGGGGGAAACCCCTCACAACACGTCTCAGCCCAAGAGCACAGAGAC'
+                'CGGATGGAAGACCCTCACGCTTTCCAGCAATACTTAGCAACCCTGCCGGACATCAACCG'
+                'AGCGGCATTGTTAAGTTTAAACAAACACTTCCGGCTGGACAATCACCAAGGATTTCTAC'
+                'AGATACTAGACAATGCCACGTTCGCGGAAATAGCAGAACATTATTACAGATTCCGCCAA'
+                'GCATATGCCAACCACAGCATAACATGGACTCCTCCCCACGAGGAGGAATACTGGGAGGA'
+                'ACACGCCCTTTATCCCGAGAAGGCGGCTCCGGTAGAACCTATCATACCGGTTCCCGGGC'
+                'CCAAAAGAGCTCCTCCTCTAAAGCCACGTCCTCAACCCCGAAACAGGACAAAGCTCAGG'
+                'AAGCCATCAGACGCCATACCTGTGCTGCTCGCAATTCCTACCAGCAGCTCAACAGAAGT'
+                'AGGCCGGCTTCCATATCAAATCTACCCCCACGAGTTACTACCTATCCCGAAGCTGGAAG'
+                'AAAAGAGTACAACTGCAAACATGAGCGGCGCTGTTTCTCCCGCTCTGTGGTTCCCAGTG'
+                'GTTGCACTGGTCGTGTATTTCTTGTTAACAAAAGCTCGGGAAATTCTAGCCAAGCTAGA'
+                'TTGGTGGTGGACTTCTCACAGTTTTCGAGGGGCAAACGGCGAGTGCATTTTCCAAGATA'
+                'CTTCTCCCCAAACACAGAAGCACTTAGCAGGATTTTGCCCCCTCACATGTGGAGGATTT'
+                'CTCTGGACCTCGCTGAGGCTTTTTACCATTTGCCTTATCATACTTCTAGCAGTCTGCGG'
+                'CTATCTGTTTCTGATGGAAAAAGGGTTTACCTCTTTAGGCGAGCAGCTATGGGAGTTGG'
+                'GCTCAGCCCCTGGAACCTTACTCTCTTCACTTCTATCATCGCTGAAAGGATCAGGAAAT'
+                'ATTTTTCTATCCATTGCTTTGCTTACATGGATGACTTCATCTTATCTCACCCTAACCCC'
+                'CGCTACCTCATTGCTGCTGCTAACGGGGTTGTGGGCTACTTACAACAATATGGGGTCAG'
+                'GATAAACTTTGACAAATCTACACTAACCCCCACCAGACAAATTGTGTTTATGGGACTTC'
+                'AAATAACACCCACACACATACATGTCACCAAAGAAAAGTTCCAGAAGCTTCATGATTTA'
+                'TTAACTTCAATCAACACCGACGTCTTCTATGATTGGAAAATTCTTCAAAGATTTGCTGG'
+                'ACATTTTAATCAGATTGTTGCTTTCTCATTGTTTGGGTCACATATTCTAGGACCCATAT'
+                'ATTTAGCTATCACTAACAAAGCTGATTTTTGTTTTACTAAATCATATGTAACTGACCTG'
+                'AGATATGTATTGAACCATGTCTATGATTTGAGACTCAGGGCTTCCCCCAAACACATTGT'
+                'ACCATCTGTAGCAGCTGACGCTACACTTACACATGGGGGAATTGCTACTCGACAGGGAG'
+                'AGGAAATTCAATTGACCTTTTCTGAGCAGAGGCCGGTCCATGTTCAGGAATTGTTAATT'
+                'GCCCTGGTGGCTCTTCTCTGTTTTGTACCGAAAGGGTTACTGTCTGACTCGACATTTGT'
+                'GGTGCACAAACAACTTCATACTTTGCCCTATGCGTTGGCTGCAGTAGGCCATTATCTCT'
+                'TAAGCCGTGCAAAGGTCACCTATGTACACACATCGCTTAACCCTGCTGACCCTGTGTCA'
+                'CGCCTCCTGGCCAATCACAGCGCTCCACGCTGGCCTTTACTTGTACACACCTCCTTGGT'
+                'TAAACCACTATATATACCATGGAGCCTACTCCGTTCTCCATAATTACACACCTCCTTCG'
+                'GAGCACCAGGT'),
+        },
+    },
+
+    'proteins':
+    {
+        'YP_009259545.1':
+        {
+            'id': (
+                'acc|GENBANK|YP_009259545.1|GENBANK|NC_030446.1|polymerase '
+                '[Tibetan frog hepatitis B virus isolate 243398, complete '
+                'genome]'),
+            'protein': (
+                'MKRLWSRNVRGRPGPILEELEGALTQAGTEDYIIRGQSGDLIVQ'
+                'GKEDLDGGVLPLQAEPERLLGPLKGLYQNDVPPFNPSWEIPTIQGVHLQEDIIREIPD'
+                'RRWAYIYPIRFWPNEVKYRPREQGITDKYPEHTRAHSEKVSLYLQTLYNNGILYKRES'
+                'RHVLLFNDNPYSWESNNLLPCHCHKSQDYEGLGRQPPSRRSNEVLRPSVELDRGKPLT'
+                'TRLSPRAQRPDGRPSRFPAILSNPAGHQPSGIVKFKQTLPAGQSPRISTDTRQCHVRG'
+                'NSRTLLQIPPSICQPQHNMDSSPRGGILGGTRPLSREGGSGRTYHTGSRAQKSSSSKA'
+                'TSSTPKQDKAQEAIRRHTCAARNSYQQLNRSRPASISNLPPRVTTYPEAGRKEYNCKH'
+                'ERRCFSRSVVPSGCTGRVFLVNKSSGNSSQARLVVDFSQFSRGKRRVHFPRYFSPNTE'
+                'ALSRILPPHMWRISLDLAEAFYHLPYHTSSSLRLSVSDGKRVYLFRRAAMGVGLSPWN'
+                'LTLFTSIIAERIRKYFSIHCFAYMDDFILSHPNPRYLIAAANGVVGYLQQYGVRINFD'
+                'KSTLTPTRQIVFMGLQITPTHIHVTKEKFQKLHDLLTSINTDVFYDWKILQRFAGHFN'
+                'QIVAFSLFGSHILGPIYLAITNKADFCFTKSYVTDLRYVLNHVYDLRLRASPKHIVPS'
+                'VAADATLTHGGIATRQGEEIQLTFSEQRPVHVQELLIALVALLCFVPKGLLSDSTFVV'
+                'HKQLHTLPYALAAVGHYLLSRAKVTYVHTSLNPADPVSRLLANHSAPRWPLLVHTSLV'
+                'KPLYIPWSLLRSP'),
+            'range': '[676:3111](+)',
+        }
+    }
+}
diff --git a/test/diamond/test_conversion.py b/test/diamond/test_conversion.py
index e8e5c01..93cad66 100644
--- a/test/diamond/test_conversion.py
+++ b/test/diamond/test_conversion.py
@@ -15,8 +15,7 @@ from ..mocking import mockOpen
 from json import dumps
 
 from dark.diamond.conversion import (
-    JSONRecordsReader, DiamondTabularFormatReader, diamondTabularFormatToDicts,
-    FIELDS)
+    JSONRecordsReader, DiamondTabularFormatReader, DiamondTabularFormat)
 from dark.reads import Reads, AARead
 
 
@@ -1147,50 +1146,49 @@ class TestDiamondTabularFormatToDicts(TestCase):
     """
     def testDuplicatesInFieldNameList(self):
         """
-        If a field name list that contains duplicates is passed, the function
-        must raise a ValueError.
+        If a field name list that contains duplicates is passed, the
+        DiamondTabularFormat __init__ function must raise a ValueError.
         """
-        error = '^fieldNames contains duplicated names: a, b\\.$'
+        error = '^field names contains duplicated names: a, b\\.$'
         assertRaisesRegex(
-            self, ValueError, error, list,
-            diamondTabularFormatToDicts(None, ['a', 'b', 'a', 'c', 'b']))
+            self, ValueError, error,
+            DiamondTabularFormat, ['a', 'b', 'a', 'c', 'b'])
 
     def testTooFewFields(self):
         """
         If an input line does not have enough fields, a ValueError must be
         raised.
         """
+        dtf = DiamondTabularFormat(['a', 'b', 'c'])
         data = StringIO('a\tb\n')
         error = (
-            "^Line 1 of input had 2 field values \\(expected 3\\)\\. " +
-            ('To provide input for this function, DIAMOND must be called ' +
-             'with "--outfmt 6 %s" \\(without the quotes\\)\\. ' % FIELDS) +
-            "The offending input line was 'a\\\\tb'\\.")
+            r"^DIAMOND output line had 2 field values \(expected 3\)\. "
+            r"The offending input line was 'a\\tb\\n'\.")
         assertRaisesRegex(
             self, ValueError, error, list,
-            diamondTabularFormatToDicts(data, ['a', 'b', 'c']))
+            dtf.diamondTabularFormatToDicts(data))
 
     def testTooManyFields(self):
         """
         If an input line has too many fields, a ValueError must be raised.
         """
+        dtf = DiamondTabularFormat(['a', 'b'])
         data = StringIO('a\tb\tc\n')
         error = (
-            "^Line 1 of input had 3 field values \\(expected 2\\)\\. " +
-            ('To provide input for this function, DIAMOND must be called ' +
-             'with "--outfmt 6 %s" \\(without the quotes\\)\\. ' % FIELDS) +
-            "The offending input line was 'a\\\\tb\\\\tc'\\.")
+            r"^DIAMOND output line had 3 field values \(expected 2\)\. "
+            r"The offending input line was 'a\\tb\\tc\\n'\.")
         assertRaisesRegex(
             self, ValueError, error, list,
-            diamondTabularFormatToDicts(data, ['a', 'b']))
+            dtf.diamondTabularFormatToDicts(data))
 
     def testUnknownField(self):
         """
         An unknown field name must result in a returned field name and value
         that are identical to those in the function call and its input string.
         """
+        dtf = DiamondTabularFormat(['__blah__'])
         data = StringIO('3.5\n')
-        (result,) = list(diamondTabularFormatToDicts(data, ['__blah__']))
+        (result,) = list(dtf.diamondTabularFormatToDicts(data))
         self.assertEqual({'__blah__': '3.5'}, result)
 
     def testConversions(self):
@@ -1214,7 +1212,8 @@ class TestDiamondTabularFormatToDicts(TestCase):
             ('3.5 1.7 1 7 4 10 12 1 2 ACGT\n'
              '3.6 1.8 2 8 5 11 13 2 3 TGCA').replace(' ', '\t') + '\n'
         )
-        (result1, result2) = list(diamondTabularFormatToDicts(data, fields))
+        dtf = DiamondTabularFormat(fields)
+        (result1, result2) = list(dtf.diamondTabularFormatToDicts(data))
 
         self.assertEqual(
             {
diff --git a/test/diamond/test_run.py b/test/diamond/test_run.py
new file mode 100644
index 0000000..8dc0555
--- /dev/null
+++ b/test/diamond/test_run.py
@@ -0,0 +1,112 @@
+from unittest import TestCase, skipUnless
+from six import assertRaisesRegex
+
+from dark.aa import CODONS
+from dark.diamond.run import DiamondExecutor, diamondInstalled
+from dark.reads import Read, Reads
+
+from .sample_proteins import SAMPLE_DATA
+
+
+@skipUnless(diamondInstalled(), 'DIAMOND is not installed')
+class TestDiamondExecutor(TestCase):
+    """
+    Test the execution of DIAMOND.
+    """
+    def testEmptyDatabase(self):
+        """
+        If no subject sequences have been added, doing a search must result in
+        a ValueError.
+        """
+        de = DiamondExecutor()
+        queries = Reads([Read('id', 'ACGT')])
+        error = '^No subject sequences in the database$'
+        assertRaisesRegex(self, ValueError, error, list, de.search(queries))
+        de.cleanup()
+
+    def testNoQueries(self):
+        """
+        If no query sequences are passed, doing a search must result in
+        a ValueError.
+        """
+        de = DiamondExecutor()
+        de.addSubject(Read('id', 'ACGT'))
+        queries = Reads()
+        error = '^No query sequences were passed$'
+        assertRaisesRegex(self, ValueError, error, list, de.search(queries))
+        de.cleanup()
+
+    def testOne(self):
+        """
+        If one query sequence is passed, doing a search must result in the
+        expected result.
+        """
+        qid = 'query'
+        qseq = 'CCAGAGCGCACATACTGAATAGCAAACGGATTCTCCTTCGGGTCACACTCAATTGGG'
+        qqual = 'EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE'
+        sid = 'subject'
+
+        with DiamondExecutor() as de:
+            de.addSubject(Read(
+                sid,
+                'MVMWMEVPCSGVDKYIDSLEKFIEKRIPVKMTGKRNNKLRVKIVAKKATKKKNEVTKLG'
+                'AALRALGGLGGGAVGSLFGNPATGSTIGTGLGAALSRWLGRGDYRVSSNSVVQQSLKGT'
+                'SSIPSMHQDGQSVVVRHKEFVTEVRGNTSFLVRGTFDINPGRAETFPWLAGVASRFQEY'
+                'KIRGLVWHYVPSSGTAVSGTDPALGTVMLQTSYRSNDIPPANKMEVLNEYWSSESVPSE'
+                'AFCHPIECDPKENPFNIQYVRTDKVPDGDSKLLYDLGTTHLCVSGQQANDVVLGDLWCT'
+                'YEIELKKPIVNSNVTSVARSAALAFTGTLDLNSWFNGSVVNFGSLDVTANVKTISFPAR'
+                'LTGRFLVTVTIIASTTFTAADLSGTPLTTNCSLFALPTGSTYLRTVLGGTTPTVNIATL'
+                'QFAVDIQDSAKTASVTCLGSFTGAATRSMVTVTPYLM'))
+            queries = Reads([Read(qid, qseq, qqual)])
+            (result,) = list(de.search(queries))
+
+        self.assertEqual(
+            {
+                'bitscore': 40.8,
+                'btop': '11AN5ST',
+                'qframe': -2,
+                'qend': 3,
+                'full_qqual': qqual,
+                'qlen': 57,
+                'full_qseq': qseq,
+                'qseqid': qid,
+                'qstart': 56,
+                'slen': 450,
+                'sstart': 241,
+                'stitle': sid,
+            },
+            result)
+
+    def testYP_009259545(self):
+        """
+        Test for a match against YP_009259545
+        """
+        proteinAccession = 'YP_009259545.1'
+        proteinSequence = SAMPLE_DATA['proteins'][proteinAccession]['protein']
+        proteinId = SAMPLE_DATA['proteins'][proteinAccession]['id']
+
+        qid = 'query'
+        qseq = ''.join(CODONS[aa][0] for aa in proteinSequence[10:50])
+        qqual = 'E' * len(qseq)
+
+        with DiamondExecutor() as de:
+            de.addSubject(Read(proteinId, proteinSequence))
+            queries = Reads([Read(qid, qseq, qqual)])
+            (result,) = list(de.search(queries))
+
+        self.assertEqual(
+            {
+                'bitscore': 83.6,
+                'btop': '40',
+                'qframe': 1,
+                'qend': 120,
+                'full_qqual': qqual,
+                'qlen': len(qseq),
+                'full_qseq': qseq,
+                'qseqid': 'query',
+                'qstart': 1,
+                'slen': len(proteinSequence),
+                'sstart': 11,
+                'stitle': proteinId,
+            },
+            result)
diff --git a/test/diamond/test_sam.py b/test/diamond/test_sam.py
new file mode 100644
index 0000000..cace660
--- /dev/null
+++ b/test/diamond/test_sam.py
@@ -0,0 +1,125 @@
+from unittest import TestCase, skipUnless, skip
+from six import StringIO
+
+from dark.aa import CODONS
+from dark.civ.proteins import SqliteIndex, SqliteIndexWriter, _Genome
+from dark.diamond.sam import SimpleDiamondSAMWriter
+from dark.diamond.run import DiamondExecutor, diamondInstalled
+from dark.genbank import GenomeRanges
+from dark.reads import Read, Reads
+
+from .sample_proteins import SAMPLE_DATA
+
+
+@skip('SAM tests still under construction')
+@skipUnless(diamondInstalled(), 'DIAMOND is not installed')
+class TestSimpleDiamondSAMWriter(TestCase):
+    """
+    Test the SimpleDiamondSAMWriter class.
+    """
+    def testTibetanFrogHBV(self):
+        """
+        Test that Tibetan frogs can get HBV.
+        """
+        proteinAccession = 'YP_009259545.1'
+        proteinSequence = SAMPLE_DATA['proteins'][proteinAccession]['protein']
+        proteinId = SAMPLE_DATA['proteins'][proteinAccession]['id']
+        proteinRange = SAMPLE_DATA['proteins'][proteinAccession]['range']
+        ranges = GenomeRanges(proteinRange)
+        queryStartInProtein = 10  # This is a 0-based amino acid offset.
+        queryLenInProtein = 40  # This is an amino acid length.
+
+        genomeAccession = 'NC_030446.1'
+        genomeSequence = SAMPLE_DATA['genomes'][genomeAccession]['genome']
+        genomeLen = len(genomeSequence)
+
+        # The query sequence is nucleotides that match the amino acids in the
+        # protein. Here we use the first ([0] in the below) codon for each
+        # amino acid to make the nucleotide sequence.
+        queryId = 'query'
+        querySequence = ''.join(
+            CODONS[aa][0]
+            for aa in proteinSequence[queryStartInProtein:
+                                      queryStartInProtein + queryLenInProtein])
+        queryQuality = 'E' * len(querySequence)
+
+        # Use the protein sequence to make a DIAMOND database and run DIAMOND
+        # on the query. Yes, this really runs DIAMOND, so you need to have it
+        # installed, with its executable somewhere in your shell's PATH.
+        with DiamondExecutor() as de:
+            de.addSubject(Read(proteinId, proteinSequence))
+            queries = Reads([Read(queryId, querySequence, queryQuality)])
+            (diamondResult,) = list(de.search(queries))
+
+        # Make sure DIAMOND gives us back what we expected.
+        self.assertEqual(
+            {
+                'bitscore': 83.6,
+                'btop': str(queryLenInProtein),  # Exact match of all AAs.
+                'qframe': 1,
+                'qend': 3 * queryLenInProtein,
+                'full_qqual': queryQuality,
+                'qlen': len(querySequence),
+                'full_qseq': querySequence,
+                'qseqid': 'query',
+                'qstart': 1,
+                'slen': len(proteinSequence),
+                'sstart': queryStartInProtein + 1,  # DIAMOND is 1-based.
+                'stitle': proteinId,
+            },
+            diamondResult)
+
+        # Make a genomes/proteins sqlite database and add information about
+        # the protein and the nucleotide genome it comes from.
+        db = SqliteIndexWriter(':memory:')
+        db.addProtein(
+            proteinAccession, genomeAccession, proteinSequence,
+            proteinRange, True, ranges.circular(genomeLen),
+            ranges.distinctRangeCount(genomeLen))
+
+        genome = _Genome(
+            # genomeAccession, genomeName, genomeSequence
+        )
+        db.addGenome(genome, 1, 'test-db')
+
+        # Make a DIAMOND-to-SAM writer and give it the DIAMOND output.
+        writer = SimpleDiamondSAMWriter(SqliteIndex(db._connection))
+        writer.addMatch(
+            '\t'.join(map(str, (
+                diamondResult['bitscore'],
+                diamondResult['btop'],
+                diamondResult['qframe'],
+                diamondResult['qend'],
+                diamondResult['full_qqual'],
+                diamondResult['qlen'],
+                diamondResult['full_qseq'],
+                diamondResult['qseqid'],
+                diamondResult['qstart'],
+                diamondResult['slen'],
+                diamondResult['sstart'],
+                diamondResult['stitle']))))
+
+        # Tell the writer to save the matches as SAM and check the result.
+        fp = StringIO()
+        writer.save(filename=fp)
+
+        flags = '31'
+        self.assertEqual(
+            '\n'.join([
+                '@SQ\tSN:%s\tLN:%d' % (genomeAccession, len(genomeSequence)),
+                '\t'.join([
+                    queryId,
+                    '0',
+                    genomeAccession,
+                    flags,
+                    '255',
+                    '120M',  # (Exact) match of 40 AAs.
+                    '*',
+                    '0',
+                    '0',
+                    querySequence,
+                    queryQuality,
+                    'AS:i:%d' % int(diamondResult['bitscore']),
+                ])
+            ]) + '\n',
+            fp.getvalue())
diff --git a/test/test_cigar.py b/test/test_cigar.py
new file mode 100644
index 0000000..4e74143
--- /dev/null
+++ b/test/test_cigar.py
@@ -0,0 +1,147 @@
+from unittest import TestCase
+from six import assertRaisesRegex
+
+from dark.cigar import CINS, CDEL, CMATCH, CEQUAL, CDIFF, dna2cigar
+
+
+class TestConstants(TestCase):
+    """
+    Test constants in dark.cigar
+    """
+    def testOperations(self):
+        """
+        Make sure the CIGAR operations have the expected one-letter codes.
+        """
+        self.assertEqual('I', CINS)
+        self.assertEqual('D', CDEL)
+        self.assertEqual('M', CMATCH)
+        self.assertEqual('=', CEQUAL)
+        self.assertEqual('X', CDIFF)
+
+
+class TestDNA2CIGAR(TestCase):
+    """
+    Test the dna2cigar function.
+    """
+    def testUnequalStrings(self):
+        """
+        If two unequal strings are passed, a ValueError must be raised.
+        """
+        error = r"^Sequences 'hey' and 'there' of unequal length \(3 != 5\)\.$"
+        assertRaisesRegex(self, ValueError, error, dna2cigar, 'hey', 'there')
+
+    def testEmptyStrings(self):
+        """
+        If two empty strings are passed, a ValueError must be raised.
+        """
+        error = r"^Two sequences of zero length were passed\.$"
+        assertRaisesRegex(self, ValueError, error, dna2cigar, '', '')
+
+    def testEqualStringsLength1(self):
+        """
+        If two equal strings of length 1 are passed, the correct CIGAR string
+        must be returned.
+        """
+        self.assertEqual('1=', dna2cigar('A', 'A'))
+
+    def testEqualStringsLength1MismatchedCase(self):
+        """
+        If two strings of length 1 are passed and the're the same but for
+        case, the CIGAR string must indicate that they do not match.
+        """
+        self.assertEqual('1X', dna2cigar('A', 'a'))
+
+    def testEqualStringsLength1Concise(self):
+        """
+        If two equal strings of length 1 are passed and the concise argument
+        is True, the correct CIGAR string must be returned.
+        """
+        self.assertEqual('1M', dna2cigar('A', 'A', concise=True))
+
+    def testUnequalStringsLength1(self):
+        """
+        If two unequal strings of length 1 are passed, the correct CIGAR
+        string must be returned.
+        """
+        self.assertEqual('1X', dna2cigar('A', 'G'))
+
+    def testUnequalStringsLength1Concise(self):
+        """
+        If two unequal strings of length 1 are passed and the concise argument
+        is True, the correct CIGAR string must be returned.
+        """
+        self.assertEqual('1M', dna2cigar('A', 'G', concise=True))
+
+    def testEqualStringsLength2(self):
+        """
+        If two equal strings of length 2 are passed, the correct CIGAR string
+        must be returned.
+        """
+        self.assertEqual('2=', dna2cigar('AT', 'AT'))
+
+    def testEqualStringsLength2Concise(self):
+        """
+        If two equal strings of length 2 are passed and the concise argument
+        is True, the correct CIGAR string must be returned.
+        """
+        self.assertEqual('2M', dna2cigar('AT', 'AT', concise=True))
+
+    def testUnequalStringsLength2(self):
+        """
+        If two unequal strings of length 2 are passed, the correct CIGAR
+        string must be returned.
+        """
+        self.assertEqual('2X', dna2cigar('AT', 'GC'))
+
+    def testUnequalStringsLength2Concise(self):
+        """
+        If two unequal strings of length 2 are passed and the concise argument
+        is True, the correct CIGAR string must be returned.
+        """
+        self.assertEqual('2M', dna2cigar('AT', 'GC', concise=True))
+
+    def testMixedEqualityStringsLength2(self):
+        """
+        If two strings of length 2 are passed and they match in just one
+        place, the correct CIGAR string must be returned.
+        """
+        self.assertEqual('1=1X', dna2cigar('AG', 'AT'))
+
+    def testMixedEqualityStringsLength2Concise(self):
+        """
+        If two equal strings of length 2 are passed and the concise argument
+        is True, the correct CIGAR string must be returned.
+        """
+        self.assertEqual('2M', dna2cigar('AG', 'AT', concise=True))
+
+    def testMixedEqualityStringsLength2Alt(self):
+        """
+        If two strings of length 2 are passed and they match in just one
+        place, the correct CIGAR string must be returned.
+        """
+        self.assertEqual('1X1=', dna2cigar('GA', 'TA'))
+
+    def testMixedEqualityStringsLength2ConciseAlt(self):
+        """
+        If two equal strings of length 2 are passed and the concise argument
+        is True, the correct CIGAR string must be returned.
+        """
+        self.assertEqual('2M', dna2cigar('GA', 'TA', concise=True))
+
+    def testMixedMatch(self):
+        """
+        If two strings with matching and non-matching regions must result
+        in the expected CIGAR string.
+        """
+        self.assertEqual('3=4X5=2X',
+                         dna2cigar('ACGTTTTCCCTTGG',
+                                   'ACGAAAACCCTTTT'))
+
+    def testMixedMatchConcise(self):
+        """
+        If two strings with matching and non-matching regions must result
+        in the expected CIGAR string when the concise argument is True.
+        """
+        self.assertEqual('14M',
+                         dna2cigar('ACGTTTTCCCTTGG',
+                                   'ACGAAAACCCTTTT', concise=True))
diff --git a/test/test_features.py b/test/test_features.py
index d50b1f6..6df6d52 100644
--- a/test/test_features.py
+++ b/test/test_features.py
@@ -17,7 +17,7 @@ except ImportError:
     else:
         raise
 else:
-    matplotlib.use('PDF')
+    matplotlib.use('Agg')
     from matplotlib import pyplot as plt
     havePlt = True
 
@@ -29,10 +29,16 @@ except ImportError:
 from Bio.SeqFeature import FeatureLocation, SeqFeature
 from Bio.SeqRecord import SeqRecord
 
+from random import uniform
+
 from dark.features import (Feature, FeatureList, _FeatureAdder,
                            ProteinFeatureAdder, NucleotideFeatureAdder)
 
 
+def _randomLabel():
+    return str(uniform(1, 1e9))
+
+
 class Test_Feature(TestCase):
     """
     Tests of the C{Feature} class.
@@ -301,7 +307,7 @@ class Test_FeatureAdder(TestCase):
             return None
 
         featureAdder = _FeatureAdder()
-        fig = plt.subplot(111)
+        fig = plt.subplot(111, label=_randomLabel())
         fig.set_title = MagicMock()
         featureAdder.add(fig, 'title', 0, 100, sequenceFetcher=fetcher)
         fig.set_title.assert_called_with('Target sequence features',
@@ -316,7 +322,7 @@ class Test_FeatureAdder(TestCase):
             return None
 
         featureAdder = _FeatureAdder()
-        fig = plt.subplot(111)
+        fig = plt.subplot(111, label=_randomLabel())
         fig.set_yticks = MagicMock()
         featureAdder.add(fig, 'title', 0, 100, sequenceFetcher=fetcher)
         fig.set_yticks.assert_called_with([])
@@ -332,7 +338,7 @@ class Test_FeatureAdder(TestCase):
             return None
 
         featureAdder = _FeatureAdder()
-        fig = plt.subplot(111)
+        fig = plt.subplot(111, label=_randomLabel())
         fig.text = MagicMock()
         fig.axis = MagicMock()
         featureAdder.add(fig, 'title', 0, 300, sequenceFetcher=fetcher)
@@ -353,7 +359,7 @@ class Test_FeatureAdder(TestCase):
 
         featureAdder = _FeatureAdder()
         featureAdder.WANTED_TYPES = ('site',)
-        fig = plt.subplot(111)
+        fig = plt.subplot(111, label=_randomLabel())
         fig.text = MagicMock()
         fig.axis = MagicMock()
         result = featureAdder.add(fig, 'title', 0, 300,
@@ -381,7 +387,7 @@ class Test_FeatureAdder(TestCase):
 
         featureAdder = _FeatureAdder()
         featureAdder.WANTED_TYPES = ('site',)
-        fig = plt.subplot(111)
+        fig = plt.subplot(111, label=_randomLabel())
         self.assertRaises(NotImplementedError, featureAdder.add, fig, 'title',
                           0, 300, sequenceFetcher=fetcher)
 
@@ -400,7 +406,7 @@ class Test_FeatureAdder(TestCase):
 
         featureAdder = _FeatureAdder()
         featureAdder.WANTED_TYPES = ('site',)
-        fig = plt.subplot(111)
+        fig = plt.subplot(111, label=_randomLabel())
         fig.text = MagicMock()
         fig.axis = MagicMock()
         result = featureAdder.add(fig, 'title', 0, 300,
@@ -434,7 +440,7 @@ class TestProteinFeatureAdder(TestCase):
             return SeqRecord(None, features=[feature])
 
         featureAdder = ProteinFeatureAdder()
-        fig = plt.subplot(111)
+        fig = plt.subplot(111, label=_randomLabel())
         fig.plot = MagicMock()
         result = featureAdder.add(fig, 'title', 0, 300,
                                   sequenceFetcher=fetcher)
@@ -455,7 +461,7 @@ class TestProteinFeatureAdder(TestCase):
             return SeqRecord(None, features=[feature])
 
         featureAdder = ProteinFeatureAdder()
-        fig = plt.subplot(111)
+        fig = plt.subplot(111, label=_randomLabel())
         fig.plot = MagicMock()
         fig.axis = MagicMock()
         fig.legend = MagicMock()
@@ -492,7 +498,7 @@ class TestNucleotideFeatureAdder(TestCase):
             return SeqRecord(None, features=[feature])
 
         featureAdder = NucleotideFeatureAdder()
-        fig = plt.subplot(111)
+        fig = plt.subplot(111, label=_randomLabel())
         fig.plot = MagicMock()
         result = featureAdder.add(fig, 'title', 0, 300,
                                   sequenceFetcher=fetcher)
@@ -513,7 +519,7 @@ class TestNucleotideFeatureAdder(TestCase):
             return SeqRecord(None, features=[feature])
 
         featureAdder = NucleotideFeatureAdder()
-        fig = plt.subplot(111)
+        fig = plt.subplot(111, label=_randomLabel())
         fig.plot = MagicMock()
         fig.axis = MagicMock()
         fig.legend = MagicMock()
@@ -548,7 +554,7 @@ class TestNucleotideFeatureAdder(TestCase):
             return SeqRecord(None, features=[feature])
 
         featureAdder = NucleotideFeatureAdder()
-        fig = plt.subplot(111)
+        fig = plt.subplot(111, label=_randomLabel())
         fig.plot = MagicMock()
         fig.axis = MagicMock()
         fig.legend = MagicMock()
@@ -587,7 +593,7 @@ class TestNucleotideFeatureAdder(TestCase):
             return SeqRecord(None, features=[feature1, feature2])
 
         featureAdder = NucleotideFeatureAdder()
-        fig = plt.subplot(111)
+        fig = plt.subplot(111, label=_randomLabel())
         fig.plot = MagicMock()
         fig.axis = MagicMock()
         fig.legend = MagicMock()
diff --git a/test/test_genbank.py b/test/test_genbank.py
new file mode 100644
index 0000000..a0e39a6
--- /dev/null
+++ b/test/test_genbank.py
@@ -0,0 +1,441 @@
+from unittest import TestCase
+from six import assertRaisesRegex
+import warnings
+
+from dark.genbank import GenomeRanges
+
+
+class TestErrors(TestCase):
+    """
+    Test incorrect ways of instantitating the GenomeRanges class.
+    """
+    def testEmptyString(self):
+        """
+        The GenomeRanges class must raise a ValueError when passed an range
+        empty string.
+        """
+        error = (r'^Could not parse GenBank range string ""\. '
+                 r'Range "" does not end with \]\(\+\) or \]\(-\)\.$')
+        assertRaisesRegex(self, ValueError, error, GenomeRanges, '')
+
+    def testUnparseableString(self):
+        """
+        The GenomeRanges class must raise a ValueError when passed an
+        unparseable string.
+        """
+        error = (r'^Could not parse GenBank range string "x"\. '
+                 r'Range "x" does not end with \]\(\+\) or \]\(-\)\.$')
+        assertRaisesRegex(self, ValueError, error, GenomeRanges, 'x')
+
+    def testNoComplementIndicator(self):
+        """
+        The GenomeRanges class must raise a ValueError when passed an a range
+        that does not end with ](+] or ](-).
+        """
+        error = (r'^Could not parse GenBank range string "\[33:40\]"\. '
+                 r'Range "\[33:40\]" does not end with \]\(\+\) or \]\(-\)\.$')
+        assertRaisesRegex(self, ValueError, error, GenomeRanges, '[33:40]')
+
+    def testNoLeadingSquareBracket(self):
+        """
+        The GenomeRanges class must raise a ValueError when passed an a range
+        that does not start with [.
+        """
+        error = (r'^Could not parse GenBank range string "33:40\]\(-\)"\. '
+                 r'Range "33:40\]\(-\)" does not start with "\["\.$')
+        assertRaisesRegex(self, ValueError, error, GenomeRanges, '33:40](-)')
+
+    def testJoinWithNoRange(self):
+        """
+        The GenomeRanges class must raise a ValueError when passed a join
+        directive containing no ranges.
+        """
+        error = (r'^Could not parse GenBank range string "join{}"\. '
+                 r'join{} can only be used with multiple ranges\.$')
+        assertRaisesRegex(self, ValueError, error, GenomeRanges, 'join{}')
+
+    def testDecreasingOrder(self):
+        """
+        The GenomeRanges class must raise a ValueError when passed a string in
+        which the start, stop values decrease.
+        """
+        error = (r'^Could not parse GenBank range string "\[7:4\]\(\+\)"\. '
+                 r'Offset values \(7, 4\) cannot decrease\.$')
+        assertRaisesRegex(self, ValueError, error, GenomeRanges, '[7:4](+)')
+
+    def testContiguousRanges(self):
+        """
+        The GenomeRanges class must raise a ValueError if passed two ranges
+        that are contiguous and which have the same (-) direction.
+        """
+        message = (
+            'Contiguous GenBank ranges detected: [3:5] followed by [5:7].')
+        with warnings.catch_warnings(record=True) as w:
+            warnings.simplefilter('always')
+            GenomeRanges('join{[3:5](-), [5:7](-)}')
+            self.assertEqual(1, len(w))
+            self.assertEqual(message, str(w[0].message))
+
+    def testContiguousRangesComplementStrand(self):
+        """
+        The GenomeRanges class must raise a ValueError if passed two ranges
+        that are contiguous and which have the same (+) direction.
+        """
+        message = (
+            'Contiguous GenBank ranges detected: [3:5] followed by [5:7].')
+        with warnings.catch_warnings(record=True) as w:
+            warnings.simplefilter('always')
+            GenomeRanges('join{[3:5](+), [5:7](+)}')
+            self.assertEqual(1, len(w))
+            self.assertEqual(message, str(w[0].message))
+
+    def testTwoNakedRanges(self):
+        """
+        Two ranges with no join qualifier must result in a ValueError.
+        """
+        error = (r'^Could not parse GenBank range string '
+                 r'"\[3:5\(\+\), 7:9\(-\)\]"\. '
+                 r'Multiple ranges must be wrapped in join{}\.$')
+        assertRaisesRegex(self, ValueError, error, GenomeRanges,
+                          '[3:5(+), 7:9(-)]')
+
+    def testOneJoinedRange(self):
+        """
+        A single range with a join() must result in a ValueError.
+        """
+        error = (r'^Could not parse GenBank range string '
+                 r'"join{\[3:5\]\(-\)}"\. '
+                 r'join{} can only be used with multiple ranges\.$')
+        assertRaisesRegex(self, ValueError, error, GenomeRanges,
+                          'join{[3:5](-)}')
+
+
+class TestRanges(TestCase):
+    """
+    Test the ranges attribute of the GenomeRanges class.
+    """
+    def testOneNakedRangePositive(self):
+        """
+        A single range on the positive strand must result in the the expected
+        ranges value stored.
+        """
+        gr = GenomeRanges('[3:5](+)')
+        self.assertEqual(((3, 5, True),), gr.ranges)
+
+    def testOneNakedRangeNegative(self):
+        """
+        A single range on the negative strand must result in the the expected
+        ranges value stored.
+        """
+        gr = GenomeRanges('[3:5](-)')
+        self.assertEqual(((3, 5, False),), gr.ranges)
+
+    def testTwoJoinedRanges(self):
+        """
+        Two joined ranges must return the expected result.
+        """
+        gr = GenomeRanges('join{[3:5](+), [7:9](-)}')
+        self.assertEqual(((3, 5, True), (7, 9, False)), gr.ranges)
+
+    def testThreeJoinedRanges(self):
+        """
+        Three joined ranges must return the expected result.
+        """
+        gr = GenomeRanges('join{[3:5](+), [7:9](-), [17:19](-)}')
+        self.assertEqual(((3, 5, True), (7, 9, False), (17, 19, False)),
+                         gr.ranges)
+
+    def testTwoJoinedContiguousRangesMismatchedStrands(self):
+        """
+        Two joined ranges that are contiguous but not on the same strand must
+        return the expected unmerged (two-range) result.
+        """
+        gr = GenomeRanges('join{[3:5](-), [5:9](+)}')
+        self.assertEqual(((3, 5, False), (5, 9, True)), gr.ranges)
+
+    def testTwoJoinedContiguousRangesComplementStrand(self):
+        """
+        Two joined ranges that are contiguous on the same (complement) strand
+        must return the expected (single-range) result.
+        """
+        message = (
+            'Contiguous GenBank ranges detected: [3:5] followed by [5:9].')
+        with warnings.catch_warnings(record=True) as w:
+            warnings.simplefilter('always')
+            gr = GenomeRanges('join{[3:5](-), [5:9](-)}')
+            self.assertEqual(((3, 9, False),), gr.ranges)
+            self.assertEqual(1, len(w))
+            self.assertEqual(message, str(w[0].message))
+
+    def testTwoJoinedContiguousRanges(self):
+        """
+        Two joined ranges that are contiguous on the same strand must return
+        the expected (single-range) result.
+        """
+        message = (
+            'Contiguous GenBank ranges detected: [3:5] followed by [5:9].')
+        with warnings.catch_warnings(record=True) as w:
+            warnings.simplefilter('always')
+            gr = GenomeRanges('join{[3:5](+), [5:9](+)}')
+            self.assertEqual(((3, 9, True),), gr.ranges)
+            self.assertEqual(1, len(w))
+            self.assertEqual(message, str(w[0].message))
+
+    def testThreeJoinedContiguousRanges(self):
+        """
+        Three joined ranges that are contiguous must return the expected
+        (single-range) result.
+        """
+        message1 = (
+            'Contiguous GenBank ranges detected: [3:5] followed by [5:7].')
+        # Note that the second warning message has a range that doesn't
+        # correspond to any of the ranges in the GenomeRanges
+        # initialization string. That's because by the time the second
+        # warning is issued the first two ranges ([3:5] and [5:7]) have
+        # been merged into one ([3:7]).
+        message2 = (
+            'Contiguous GenBank ranges detected: [3:7] followed by [7:9].')
+        with warnings.catch_warnings(record=True) as w:
+            warnings.simplefilter('always')
+            gr = GenomeRanges('join{[3:5](-), [5:7](-), [7:9](-)}')
+            self.assertEqual(((3, 9, False),), gr.ranges)
+            self.assertEqual(2, len(w))
+            self.assertEqual(message1, str(w[0].message))
+            self.assertEqual(message2, str(w[1].message))
+
+    def testTwoJoinedContiguousRangesInTheMiddleOfFourRanges(self):
+        """
+        Two joined ranges that are contiguous must return the expected
+        result when there are non-contiguous ranges surrounding them.
+        """
+        message = (
+            'Contiguous GenBank ranges detected: [3:5] followed by [5:9].')
+        with warnings.catch_warnings(record=True) as w:
+            warnings.simplefilter('always')
+            gr = GenomeRanges('join{[0:2](-), [3:5](-), [5:9](-), [12:15](-)}')
+            self.assertEqual(((0, 2, False), (3, 9, False), (12, 15, False)),
+                             gr.ranges)
+            self.assertEqual(1, len(w))
+            self.assertEqual(message, str(w[0].message))
+
+
+class TestCircular(TestCase):
+    """
+    Test the GenomeRanges class ranges method.
+    """
+    def testOneRange(self):
+        """
+        The circular method must return False when given only a single range
+        tuple that is fully contained within the genome.
+        """
+        self.assertFalse(GenomeRanges('[20:40](+)').circular(100))
+
+    def testOneRangeEndingAtGenomeEnd(self):
+        """
+        The circular method must return False when given only a single range
+        tuple that ends at the end of the genome.
+        """
+        self.assertFalse(GenomeRanges('[20:40](+)').circular(40))
+
+    def testOneRangeSpanningTheWholeGenome(self):
+        """
+        The circular method must return False when given only a single range
+        tuple that ends at the end of the genome.
+        """
+        self.assertFalse(GenomeRanges('[0:40](+)').circular(40))
+
+    def testTwoRangesThatAreCircular(self):
+        """
+        The circular method must return True when given two ranges that span
+        the end of the genome.
+        """
+        self.assertTrue(
+            GenomeRanges('join{[20:40](+), [0:10](+)}').circular(40))
+
+    def testThreeRangesThatAreCircular(self):
+        """
+        The circular method must return True when given three ranges that span
+        the end of the genome.
+        """
+        self.assertTrue(
+            GenomeRanges(
+                'join{[20:30](+), [31:40](+), [0:10](+)}').circular(40))
+
+    def testThreeRangesThatAreCircular2(self):
+        """
+        The circular method must return True when given three ranges that span
+        the end of the genome when one of the ranges is the last in the passed
+        BioPython GenBank string and the second range (starting at zero) is the
+        first.
+        """
+        self.assertTrue(
+            GenomeRanges(
+                'join{[0:10](+), [40:60](+), [80:90](+)}').circular(90))
+
+
+class TestStartInGenome(TestCase):
+    """
+    Test the GenomeRanges class startInGenome method.
+    """
+    def testStartOutsideGenomeRange(self):
+        """
+        If the start location in the DIAMOND match is greater than the length
+        of the genome, a ValueError must be raised.
+        """
+        # The starting offset in the genome is 27 because DIAMOND returns
+        # 1-based offsets into the protein subject. So in Python terms the
+        # start offset is amino acid 9, or nucleotide 27.
+        ranges = GenomeRanges('[10:20](+)')
+        error = (r'^Starting nucleotide offset 27 not found in protein '
+                 r'nucleotide ranges \(10, 20\)\.$')
+        assertRaisesRegex(self, ValueError, error, ranges.startInGenome,
+                          {'sstart': 10})
+
+    def testStartOutsideGenomeRangeTwoRanges(self):
+        """
+        If the start location in the DIAMOND match is greater than the length
+        of the genome when two ranges are given, a ValueError must be raised.
+        """
+        # The starting offset in the genome is 27 because DIAMOND returns
+        # 1-based offsets into the protein subject. So in Python terms the
+        # start offset is amino acid 9, or nucleotide 27.
+        ranges = GenomeRanges('join{[5:10](+), [12:18](+)}')
+        error = (r'^Starting nucleotide offset 27 not found in protein '
+                 r'nucleotide ranges \(5, 10\), \(12, 18\)\.$')
+        assertRaisesRegex(self, ValueError, error, ranges.startInGenome,
+                          {'sstart': 10})
+
+    def testMatchAtStartOfFirstRange(self):
+        """
+        If the start location in the DIAMOND match is at the very beginning of
+        the first range, the correct offset must be returned.
+        """
+        # The offset is 5 because the match starts 0 nucleotides into the
+        # protein (0 = (1 - 1 ) * 3) and the protein starts at position 5 in
+        # the genome. So the match begins at nucleotide 0 + 5 = 5.
+        ranges = GenomeRanges('[5:35](+)')
+        self.assertEqual(5, ranges.startInGenome({'sstart': 1}))
+
+    def testMatchAtEndOfFirstRange(self):
+        """
+        If the start location in the DIAMOND match is at the very end of
+        the first range, the correct offset must be returned.
+        """
+        # The offset is 20 because the match starts 15 nucleotides into the
+        # protein (15 = (6 - 1 ) * 3) and the protein starts at position 5 in
+        # the genome. So the match begins at nucleotide 15 + 5 = 20.
+        ranges = GenomeRanges('[5:21](+)')
+        self.assertEqual(20, ranges.startInGenome({'sstart': 6}))
+
+    def testInFirstRange(self):
+        """
+        If the start location in the DIAMOND match is in the first range,
+        the correct offset must be returned.
+        """
+        # The offset is 32 because the match starts 27 nucleotides into the
+        # protein (27 = (10 - 1 ) * 3) and the protein starts at position 5 in
+        # the genome. So the match begins at nucleotide 27 + 5 = 32.
+        ranges = GenomeRanges('[5:35](+)')
+        self.assertEqual(32, ranges.startInGenome({'sstart': 10}))
+
+    def testMatchAtStartOfSecondRange(self):
+        """
+        If the start location in the DIAMOND match is at the very beginning
+        of the second range, the correct offset must be returned.
+        """
+        # The offset is 50 because the match starts 30 nucleotides into the
+        # protein (45 = (16 - 1 ) * 3) and the protein has a range of 30
+        # nucleotides (35 - 5 = 30) and then a range of 35 nucleotides (85 -
+        # 50 = 35). So the match begins 15 (45 - 30 = 15) nucleotides into
+        # the second range (which starts at 50), and 15 + 50 = 65.
+        ranges = GenomeRanges('join{[5:35](+), [50:85](+)}')
+        self.assertEqual(50, ranges.startInGenome({'sstart': 11}))
+
+    def testInSecondRange(self):
+        """
+        If the start location in the DIAMOND match is in the second range,
+        the correct offset must be returned.
+        """
+        # The offset is 65 because the match starts 45 nucleotides into the
+        # protein (45 = (16 - 1 ) * 3) and the protein has a range of 30
+        # nucleotides (35 - 5 = 30) and then a range of 35 nucleotides (85 -
+        # 50 = 35). So the match begins 15 (45 - 30 = 15) nucleotides into
+        # the second range (which starts at 50), and 15 + 50 = 65.
+        ranges = GenomeRanges('join{[5:35](+), [50:85](+)}')
+        self.assertEqual(65, ranges.startInGenome({'sstart': 16}))
+
+    def testInThirdRange(self):
+        """
+        If the start location in the DIAMOND match is in the third range,
+        the correct offset must be returned.
+        """
+        # The offset is 2900 because the match starts 1200 nucleotides into
+        # the protein (1200 = (401 - 1 ) * 3) and the protein has a range of
+        # 100 nucleotides (100 - 0 = 100), then a range of 200 nucleotides
+        # (600 - 400 = 200), then a range of 1000 nucleotides. So the match
+        # begins 900 (1200 - 300 = 900) nucleotides into the third range
+        # (which starts at 2000), and 2000 + 900 = 2900.
+        ranges = GenomeRanges('join{[0:100](+), [400:600](+), [2000:3000](+)}')
+        self.assertEqual(2900, ranges.startInGenome({'sstart': 401}))
+
+
+class TestOrientations(TestCase):
+    """
+    Test the GenomeRanges class orientations method.
+    """
+    def testAllComplement(self):
+        """
+        If all ranges are on the complement strand, the orientations method
+        must return just True.
+        """
+        ranges = GenomeRanges('join{[0:100](-), [400:600](-)}')
+        self.assertEqual({False}, ranges.orientations())
+
+    def testNoComplement(self):
+        """
+        If no ranges are on the complement strand, the orientations method
+        must return just False.
+        """
+        ranges = GenomeRanges('join{[0:100](+), [400:600](+)}')
+        self.assertEqual({True}, ranges.orientations())
+
+    def testMixed(self):
+        """
+        If the ranges are on the both the regular and the complement strand,
+        the orientations method must return the set {True, False}.
+        """
+        ranges = GenomeRanges('join{[0:100](+), [400:600](-)}')
+        self.assertEqual({True, False}, ranges.orientations())
+
+
+class TestDistinctRangeCount(TestCase):
+    """
+    Test the GenomeRanges class distinctRangeCount method.
+    """
+    def testOneRange(self):
+        """
+        If there is only one range, 1 must be returned.
+        """
+        ranges = GenomeRanges('[0:100](+)')
+        self.assertEqual(1, ranges.distinctRangeCount(100))
+
+    def testTwoRange(self):
+        """
+        If there are 2 ranges, 2 must be returned.
+        """
+        ranges = GenomeRanges('join{[0:100](+), [400:600](+)}')
+        self.assertEqual(2, ranges.distinctRangeCount(100))
+
+    def testThreeRange(self):
+        """
+        If there are 3 ranges, 3 must be returned.
+        """
+        ranges = GenomeRanges('join{[0:100](+), [400:600](+), [2000:3000](+)}')
+        self.assertEqual(3, ranges.distinctRangeCount(100))
+
+    def testThreeRangesCircular(self):
+        """
+        If there are 3 ranges but the genome is circular, 2 must be returned.
+        """
+        ranges = GenomeRanges('join{[0:100](+), [400:600](+), [2000:3000](+)}')
+        self.assertEqual(2, ranges.distinctRangeCount(3000))
diff --git a/test/test_taxonomy.py b/test/test_taxonomy.py
index 0909fd8..8c52802 100644
--- a/test/test_taxonomy.py
+++ b/test/test_taxonomy.py
@@ -1,6 +1,9 @@
 from unittest import TestCase
+import six
+import sqlite3
 
-from dark.taxonomy import LineageFetcher
+from dark.taxonomy import (
+    LineageFetcher, AccessionLineageFetcher, isRetrovirus, isRNAVirus)
 
 
 class FakeCursor(object):
@@ -61,3 +64,213 @@ class TestLineageFetcher(TestCase):
                 (2, 'Vira'),
             ],
             lineage)
+
+
+class TestAccessionLineageFetcher(TestCase):
+    """
+    Test the AccessionLineageFetcher class.
+    """
+
+    def _makeFetcher(self, data):
+        """
+        Make a taxonomy database and put it into an AccessionLineageFetcher
+        instance.
+
+        @param data: A C{dict} with 'taxids', 'names', and 'nodes' keys,
+            each of which is a C{list} of C{tuple}s containing values to
+            insert into the respective tables.
+        @return: A C{AccessionLineageFetcher} instance.
+        """
+        db = sqlite3.connect(':memory:')
+        cursor = db.cursor()
+
+        # The tables created here must be identical to those expected by
+        # the AccessionLineageFetcher class. See also
+        # https://github.com/acorg/ncbi-taxonomy-database
+        cursor.executescript('''
+            CREATE TABLE nodes (
+                taxid INTEGER NOT NULL,
+                parent_taxid INTEGER NOT NULL,
+                rank VARCHAR NOT NULL
+            );
+
+            CREATE TABLE accession_taxid (
+                accession VARCHAR UNIQUE PRIMARY KEY,
+                taxid INTEGER NOT NULL
+            );
+
+            CREATE TABLE names (
+                taxid INTEGER NOT NULL,
+                name VARCHAR NOT NULL
+            );
+
+            CREATE INDEX nodes_idx ON nodes(taxid);
+            CREATE INDEX accession_idx ON accession_taxid(accession);
+            CREATE INDEX name_idx ON names(taxid);
+        ''')
+
+        execute = cursor.execute
+
+        for accession, taxid in data.get('taxids', []):
+            execute('INSERT INTO accession_taxid VALUES (?, ?)',
+                    (accession, taxid))
+
+        for taxid, name in data.get('names', []):
+            execute('INSERT INTO names VALUES (?, ?)',
+                    (taxid, name))
+
+        for taxid, parentTaxid, rank in data.get('nodes', []):
+            execute('INSERT INTO nodes VALUES (?, ?, ?)',
+                    (taxid, parentTaxid, rank))
+
+        return AccessionLineageFetcher(db)
+
+    def testUnknownTaxid(self):
+        """
+        If a taxonomy id is not present in the accession_taxid table, the
+        lineage fetcher must return C{None}.
+        """
+        fetcher = self._makeFetcher({})
+        self.assertIsNone(fetcher.lineage('DQ011818.1'))
+        fetcher.close()
+
+    def testReusingClosedFetcher(self):
+        """
+        Trying to re-use a closed fetcher must result in an AttributeError
+        being raised.
+        """
+        fetcher = self._makeFetcher({})
+        fetcher.close()
+        error = "^'NoneType' object has no attribute 'cursor'$"
+        six.assertRaisesRegex(self, AttributeError, error, fetcher.lineage,
+                              'DQ011818.1')
+
+    def testTaxidNotInNamesTable(self):
+        """
+        If a taxonomy id is not present in the names table, a ValueError must
+        be raised.
+        """
+        fetcher = self._makeFetcher({
+            'taxids': (
+                ('DQ011818.1', 500),
+            ),
+        })
+        error = '^Could not find taxonomy id 500 in names table$'
+        six.assertRaisesRegex(self, ValueError, error, fetcher.lineage,
+                              'DQ011818.1')
+        fetcher.close()
+
+    def testTaxidNotInNodesTable(self):
+        """
+        If a taxonomy id is not present in the nodes table, a ValueError must
+        be raised.
+        """
+        fetcher = self._makeFetcher({
+            'taxids': (
+                ('DQ011818.1', 500),
+            ),
+            'names': (
+                (500, 'Porcine endogenous retrovirus A'),
+            ),
+        })
+        error = '^Could not find taxonomy id 500 in nodes table$'
+        six.assertRaisesRegex(self, ValueError, error, fetcher.lineage,
+                              'DQ011818.1')
+        fetcher.close()
+
+    def testFullLookup(self):
+        """
+        The full taxonomy of an accession number must be retrievable.
+        """
+        fetcher = self._makeFetcher({
+            'taxids': (
+                ('DQ011818.1', 500),
+            ),
+            'names': (
+                (500, 'Porcine endogenous retrovirus A'),
+                (400, 'Retroviruses'),
+                (300, 'Viruses'),
+            ),
+            'nodes': (
+                (500, 400, 'species'),
+                (400, 300, 'genus'),
+                (300, 1, 'realm')
+            ),
+        })
+
+        self.assertEqual(
+            (
+                (500, 'Porcine endogenous retrovirus A', 'species'),
+                (400, 'Retroviruses', 'genus'),
+                (300, 'Viruses', 'realm'),
+            ),
+            fetcher.lineage('DQ011818.1'))
+        fetcher.close()
+
+
+class TestIsRetrovirus(TestCase):
+    """
+    Test the isRetrovirus function.
+    """
+    def testYes(self):
+        """
+        A retrovirus must be detected.
+        """
+        self.assertTrue(
+            isRetrovirus((
+                (11676, 'Human immunodeficiency virus 1', 'species'),
+                (11646, 'Lentivirus', 'genus'),
+                (327045, 'Orthoretrovirinae', 'subfamily'),
+                (11632, 'Retroviridae', 'family'),
+                (2169561, 'Ortervirales', 'order'),
+                (10239, 'Viruses', 'superkingdom'),
+            )))
+
+    def testNo(self):
+        """
+        A non-retrovirus must be detected.
+        """
+        self.assertFalse(
+            isRetrovirus((
+                (10310, 'Human alphaherpesvirus 2', 'species'),
+                (10294, 'Simplexvirus', 'genus'),
+                (10293, 'Alphaherpesvirinae', 'subfamily'),
+                (10292, 'Herpesviridae', 'family'),
+                (548681, 'Herpesvirales', 'order'),
+                (10239, 'Viruses', 'superkingdom'),
+            )))
+
+
+class TestIsRNAVirus(TestCase):
+    """
+    Test the isRNAVirus function.
+    """
+    def testByRealm(self):
+        """
+        Anything in the Riboviria realm must be classified as an RNA virus.
+        """
+        self.assertTrue(
+            isRNAVirus((
+                (11234, 'Measles morbillivirus', 'species'),
+                (11229, 'Morbillivirus', 'genus'),
+                (2560076, 'Orthoparamyxovirinae', 'subfamily'),
+                (11158, 'Paramyxoviridae', 'family'),
+                (11157, 'Mononegavirales', 'order'),
+                (2497574, 'Monjiviricetes', 'class'),
+                (2497570, 'Haploviricotina', 'subphylum'),
+                (2497569, 'Negarnaviricota', 'phylum'),
+                (2559587, 'Riboviria', 'realm'),
+                (10239, 'Viruses', 'superkingdom'),
+            )))
+
+    def testHBV(self):
+        """
+        An HBV virus must not be classified as an RNA virus.
+        """
+        self.assertFalse(
+            isRNAVirus((
+                (1508712, 'Tent-making bat hepatitis B virus', 'species'),
+                (10405, 'Orthohepadnavirus', 'genus'),
+                (10404, 'Hepadnaviridae', 'family'),
+                (10239, 'Viruses', 'superkingdom'),
+            )))