Repository: acorg/dark-matter

Issue:
Factor the common filtering code out of sam.py and filter-sam.py
I'm writing another SAM utility that also needs filtering, so time to factor out the common stuff in `sam.py` and `filter-sam.py`

Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 1d79db9..c91e688 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,10 @@
+## 3.0.34 July 18, 2018
+
+Factored common SAM filtering code out into `dark.sam.SAMFilter`. Split
+common FASTA command-line options into those for filtering (this is just
+inclusion/exclusion) and those for editing (which change the FASTA records
+in some way).
+
 ## 3.0.33 July 14, 2018
 
 Added `compare-consensuses.py` script.
diff --git a/bin/fasta-identity-table.py b/bin/fasta-identity-table.py
index c67cc1e..5bc7d01 100755
--- a/bin/fasta-identity-table.py
+++ b/bin/fasta-identity-table.py
@@ -9,7 +9,8 @@ from operator import itemgetter
 
 from dark.dna import compareDNAReads
 from dark.filter import (
-    addFASTAFilteringCommandLineOptions, parseFASTAFilteringCommandLineOptions)
+    addFASTAFilteringCommandLineOptions, parseFASTAFilteringCommandLineOptions,
+    addFASTAEditingCommandLineOptions, parseFASTAEditingCommandLineOptions)
 from dark.reads import addFASTACommandLineOptions, parseFASTACommandLineOptions
 
 
@@ -504,16 +505,20 @@ if __name__ == '__main__':
 
     addFASTACommandLineOptions(parser)
     addFASTAFilteringCommandLineOptions(parser)
+    addFASTAEditingCommandLineOptions(parser)
     args = parser.parse_args()
 
     colors = parseColors(args.color, args.defaultColor)
     # Sanity check - the last threshold must be zero.
     assert colors[-1][0] == 0.0
 
+    reads = parseFASTAEditingCommandLineOptions(
+        args, parseFASTAFilteringCommandLineOptions(
+            args, parseFASTACommandLineOptions(args)))
+
     # Collect the reads into a dict, keeping the insertion order.
     reads1 = OrderedDict()
-    for read in parseFASTAFilteringCommandLineOptions(
-            args, parseFASTACommandLineOptions(args)):
+    for read in reads:
         reads1[read.id] = read
 
     if args.fastaFile2:
diff --git a/bin/filter-fasta.py b/bin/filter-fasta.py
index cc9aa05..1c7d0f1 100755
--- a/bin/filter-fasta.py
+++ b/bin/filter-fasta.py
@@ -5,7 +5,8 @@ from __future__ import print_function, division
 import sys
 
 from dark.filter import (
-    addFASTAFilteringCommandLineOptions, parseFASTAFilteringCommandLineOptions)
+    addFASTAFilteringCommandLineOptions, parseFASTAFilteringCommandLineOptions,
+    addFASTAEditingCommandLineOptions, parseFASTAEditingCommandLineOptions)
 from dark.reads import addFASTACommandLineOptions, parseFASTACommandLineOptions
 
 
@@ -35,9 +36,13 @@ if __name__ == '__main__':
 
     addFASTACommandLineOptions(parser)
     addFASTAFilteringCommandLineOptions(parser)
+    addFASTAEditingCommandLineOptions(parser)
+
     args = parser.parse_args()
-    reads = parseFASTAFilteringCommandLineOptions(
-        args, parseFASTACommandLineOptions(args))
+
+    reads = parseFASTAEditingCommandLineOptions(
+        args, parseFASTAFilteringCommandLineOptions(
+            args, parseFASTACommandLineOptions(args)))
 
     saveAs = (
         args.saveAs or
diff --git a/bin/filter-sam.py b/bin/filter-sam.py
index 363fed1..ddc34a6 100755
--- a/bin/filter-sam.py
+++ b/bin/filter-sam.py
@@ -7,8 +7,8 @@ from pysam import AlignmentFile
 
 from dark.filter import (
     addFASTAFilteringCommandLineOptions, parseFASTAFilteringCommandLineOptions)
-from dark.reads import Read, Reads
-from dark.sam import samfile
+from dark.reads import Reads
+from dark.sam import samfile, SAMFilter
 
 
 if __name__ == '__main__':
@@ -19,10 +19,6 @@ if __name__ == '__main__':
         description=('Given a SAM/BAM file and a set of filtering criteria '
                      'write filtered SAM/BAM to stdout.'))
 
-    parser.add_argument(
-        'samfile',
-        help='The SAM/BAM file to filter.')
-
     parser.add_argument(
         '--quiet', action='store_true', default=False,
         help='If True, do not print the final summary.')
@@ -35,102 +31,67 @@ if __name__ == '__main__':
         '--checkResultCount', type=int,
         help=('The number of alignments expected in the output. If this '
               'number is not seen, the script exits with status 1 (and an '
-              'error message is printed unless --quiet is used).'))
+              'error message is also printed, unless --quiet was used).'))
 
-    parser.add_argument(
-        '--dropUnmapped', default=False, action='store_true',
-        help='If given, unmapped matches will not be output.')
-
-    parser.add_argument(
-        '--dropSecondary', default=False, action='store_true',
-        help='If given, secondary matches will not be output.')
-
-    parser.add_argument(
-        '--dropSupplementary', default=False, action='store_true',
-        help='If given, supplementary matches will not be output.')
-
-    parser.add_argument(
-        '--dropDuplicates', default=False, action='store_true',
-        help=('If given, matches flagged as optical or PCR duplicates will '
-              'not be output.'))
-
-    parser.add_argument(
-        '--keepQCFailures', default=False, action='store_true',
-        help=('If given, reads that are considered quality control failures '
-              'will be included in the output.'))
-
-    parser.add_argument(
-        '--referenceWhitelist', metavar='NAME', action='append',
-        help=('A reference sequence id whose alignments should be output. '
-              'If omitted, alignments against all references will be output. '
-              'May be repeated.'))
-
-    parser.add_argument(
-        '--referenceBlacklist', metavar='NAME', action='append',
-        help=('A reference sequence id whose alignments should not be output. '
-              'If omitted, alignments against all references will be output. '
-              'May be repeated.'))
-
-    addFASTAFilteringCommandLineOptions(parser, filteringSAM=True)
+    addFASTAFilteringCommandLineOptions(parser)
+    SAMFilter.addFilteringOptions(parser)
 
     args = parser.parse_args()
-
-    referenceWhitelist = (set(args.referenceWhitelist)
-                          if args.referenceWhitelist else None)
-    referenceBlacklist = (set(args.referenceBlacklist)
-                          if args.referenceBlacklist else None)
-
-    if referenceWhitelist and referenceBlacklist and (
-            referenceWhitelist & referenceBlacklist):
-        raise ValueError(
-            'The reference whitelist and blacklist sets are not '
-            'disjoint. The intersection contains %s.' %
-            ', '.join(sorted(referenceWhitelist & referenceBlacklist)))
-
-    reads = parseFASTAFilteringCommandLineOptions(args, Reads(),
-                                                  filteringSAM=True)
-    filterRead = reads.filterRead
-    dropUnmapped = args.dropUnmapped
-    dropSecondary = args.dropSecondary
-    dropSupplementary = args.dropSupplementary
-    dropDuplicates = args.dropDuplicates
-    keepQCFailures = args.keepQCFailures
-
-    kept = total = 0
-    with samfile(args.samfile) as samAlignment:
-        out = AlignmentFile(sys.stdout, mode='w' + args.bam,
-                            template=samAlignment)
-        save = out.write
-        for alignment in samAlignment.fetch():
-            total += 1
-
-            if (filterRead(Read(alignment.query_name,
-                                alignment.query_sequence,
-                                alignment.qual)) and
-                not (
-                    (alignment.is_unmapped and dropUnmapped) or
-                    (alignment.is_secondary and dropSecondary) or
-                    (alignment.is_supplementary and dropSupplementary) or
-                    (alignment.is_duplicate and dropDuplicates) or
-                    (alignment.is_qcfail and not keepQCFailures) or
-                    (referenceWhitelist is not None and
-                     alignment.reference_name not in referenceWhitelist) or
-                    (referenceBlacklist is not None and
-                     alignment.reference_name in referenceBlacklist))):
-                kept += 1
-                save(alignment)
-        out.close()
+    reads = parseFASTAFilteringCommandLineOptions(args, Reads())
+    samFilter = SAMFilter.parseFilteringOptions(args, reads.filterRead,
+                                                storeQueryIds=True)
+
+    if samFilter.referenceIds:
+        # Make a header that only includes the wanted reference ids (if
+        # any).
+        print('If you get a segmentation violation, this is (I think) '
+              'due to the following pysam issue '
+              'https://github.com/pysam-developers/pysam/issues/716 It is '
+              'safer to run without using --referenceId.', file=sys.stderr)
+        with samfile(args.samfile) as sam:
+            header = sam.header.as_dict()
+        referenceIds = samFilter.referenceIds
+        sequences = [sequence for sequence in header['SQ']
+                     if sequence['SN'] in referenceIds]
+        referenceIdsFound = set(sequence['SN'] for sequence in sequences)
+        notFound = referenceIds - referenceIdsFound
+        if notFound:
+            if len(notFound) == len(referenceIds):
+                raise ValueError(
+                    'None of the provided --referenceId ids were found in %s' %
+                    args.samfile)
+            else:
+                # Just warn about the missing references.
+                print(
+                    'WARNING: reference%s (%s) given with --referenceId not '
+                    'found in %s.' % (
+                        '' if len(notFound) == 1 else 's',
+                        ', '.join(sorted(notFound)), args.samfile),
+                    file=sys.stderr)
+
+        header['SQ'] = sequences
+        out = AlignmentFile(sys.stdout, mode='w' + args.bam, header=header)
+    else:
+        with samfile(args.samfile) as sam:
+            out = AlignmentFile(sys.stdout, mode='w' + args.bam, template=sam)
+
+    save = out.write
+    kept = 0
+    for kept, alignment in enumerate(samFilter.alignments(), start=1):
+        save(alignment)
+
+    out.close()
 
     if not args.quiet:
+        total = samFilter.alignmentCount
         print('Read %d alignment%s, kept %d (%.2f%%).' %
               (total, '' if total == 1 else 's', kept,
                0.0 if total == 0 else kept / total * 100.0), file=sys.stderr)
 
-    if args.checkResultCount is not None:
-        if kept != args.checkResultCount:
-            if not args.quiet:
-                print('Did not write the expected %d alignment%s (wrote %d).' %
-                      (args.checkResultCount,
-                       '' if args.checkResultCount == 1 else 's', kept),
-                      file=sys.stderr)
-            sys.exit(1)
+    if args.checkResultCount is not None and kept != args.checkResultCount:
+        if not args.quiet:
+            print('Did not write the expected %d alignment%s (wrote %d).' %
+                  (args.checkResultCount,
+                   '' if args.checkResultCount == 1 else 's', kept),
+                  file=sys.stderr)
+        sys.exit(1)
diff --git a/bin/sam-references.py b/bin/sam-references.py
new file mode 100755
index 0000000..e11ed60
--- /dev/null
+++ b/bin/sam-references.py
@@ -0,0 +1,17 @@
+#!/usr/bin/env python
+
+from __future__ import print_function
+
+import argparse
+
+from dark.sam import samReferencesToStr
+
+parser = argparse.ArgumentParser(
+    formatter_class=argparse.ArgumentDefaultsHelpFormatter,
+    description='Print SAM/BAM file reference names and lengths.')
+
+parser.add_argument(
+    'samfile', metavar='FILENAME',
+    help='The name of a SAM/BAM alignment file.')
+
+print(samReferencesToStr(parser.parse_args().samfile))
diff --git a/bin/sam-to-fasta-alignment.py b/bin/sam-to-fasta-alignment.py
index 8c34984..b7735b9 100755
--- a/bin/sam-to-fasta-alignment.py
+++ b/bin/sam-to-fasta-alignment.py
@@ -9,29 +9,16 @@ from __future__ import division, print_function
 import sys
 import argparse
 
-from dark.sam import (
-    PaddedSAM, UnequalReferenceLengthError, UnknownReference)
+from dark.filter import (
+    addFASTAFilteringCommandLineOptions, parseFASTAFilteringCommandLineOptions)
+from dark.reads import Reads
+from dark.sam import SAMFilter, PaddedSAM
+from dark.utils import nucleotidesToStr
 
 parser = argparse.ArgumentParser(
     formatter_class=argparse.ArgumentDefaultsHelpFormatter,
     description='Produce aligned FASTA queries from a SAM/BAM file.')
 
-parser.add_argument(
-    'samFile', help='The SAM/BAM file to read.')
-
-parser.add_argument(
-    '--minLength', type=int, default=0,
-    help='Do not output anything for queries shorter than this.')
-
-parser.add_argument(
-    '--referenceName',
-    help=('The name of the reference sequence to print alignments for. '
-          'This is only needed if the SAM/BAM alignment was made against '
-          'multiple references *and* they are not all of the same length. '
-          'If there is only one reference sequence or if all reference '
-          'sequences are of the same length, there is no need to provide a '
-          'reference name.'))
-
 parser.add_argument(
     '--rcSuffix', default='',
     help=('A string to add to the end of query names that are reverse '
@@ -40,33 +27,15 @@ parser.add_argument(
           '--allowDuplicateIds is not used)'))
 
 parser.add_argument(
-    '--dropSecondary', default=False, action='store_true',
-    help='If given, secondary matches will not be output.')
-
-parser.add_argument(
-    '--dropSupplementary', default=False, action='store_true',
-    help='If given, supplementary matches will not be output.')
-
-parser.add_argument(
-    '--dropDuplicates', default=False, action='store_true',
-    help=('If given, matches flagged as optical or PCR duplicates will '
-          'not be output.'))
-
-parser.add_argument(
-    '--allowDuplicateIds', default=False, action='store_true',
-    help=('If given, repeated query ids (due to secondary or supplemental '
-          'matches) will not have /1, /2, etc. appended to their ids. So '
-          'repeated ids may appear in the output FASTA.'))
-
-parser.add_argument(
-    '--keepQCFailures', default=False, action='store_true',
-    help=('If given, reads that are considered quality control failures will '
-          'be included in the output.'))
-
-parser.add_argument(
-    '--listReferenceNames', default=False, action='store_true',
-    help=('If given, the names of all reference sequences will be '
-          'printed to standard output and the program will exit.'))
+    '--rcNeeded', default=False, action='store_true',
+    help=('If given, queries that are flagged as matching when reverse '
+          'complemented will be reverse complemented in the output. This '
+          'must be used if the program that created the SAM/BAM input '
+          'flags reversed matches but does not also store the reverse '
+          'complemented query. The bwa program (mem and aln followed by '
+          'samse) stores the queries reversed complemented if the match '
+          'was, so this option is not needed for bwa. If in doubt, test the '
+          'output of your matching program as this is very important!'))
 
 parser.add_argument(
     '--listReferenceInsertions', default=False, action='store_true',
@@ -85,86 +54,23 @@ parser.add_argument(
           'indicates that 13 query (3 with T and 10 with G) matches would '
           'insert a nucleotide into the reference at offset 27.'))
 
-parser.add_argument(
-    '--rcNeeded', default=False, action='store_true',
-    help=('If given, queries that are flagged as matching when reverse '
-          'complemented will be reverse complemented in the output. This '
-          'must be used if the program that created the SAM/BAM input '
-          'flags reversed matches but does not also store the reverse '
-          'complemented query. The bwa program (mem and aln followed by '
-          'samse) stores the queries reversed complemented if the match '
-          'was, so this option is not needed for bwa. If in doubt, test the '
-          'output of your matching program as this is very important!'))
-
-
-def baseCountsToStr(counts):
-    """
-    Convert base counts to a string.
-
-    @param counts: A C{Counter} instance.
-    @return: A C{str} representation of nucleotide counts at an offset.
-    """
-    return ' '.join([
-        ('%s:%d' % (base, counts[base])) for base in sorted(counts)])
-
-
-def nucleotidesToStr(nucleotides, prefix=''):
-    """
-    Convert offsets and base counts to a string.
-
-    @param nucleotides: A C{defaultdict(Counter)} instance, keyed
-        by C{int} offset, with nucleotides keying the Counters.
-    @param prefix: A C{str} to put at the start of each line.
-    @return: A C{str} representation of the offsets and nucleotide
-        counts for each.
-    """
-    result = []
-    for offset in sorted(nucleotides):
-        result.append(
-            '%s%d: %s' % (prefix, offset,
-                          baseCountsToStr(nucleotides[offset])))
-    return '\n'.join(result)
-
+SAMFilter.addFilteringOptions(parser)
+addFASTAFilteringCommandLineOptions(parser)
 
 args = parser.parse_args()
-
-paddedSAM = PaddedSAM(args.samFile)
-
-try:
-    if args.listReferenceNames:
-        print(paddedSAM.referencesToStr(), file=sys.stdout)
+reads = parseFASTAFilteringCommandLineOptions(args, Reads())
+samFilter = SAMFilter.parseFilteringOptions(
+    args, filterRead=reads.filterRead)
+paddedSAM = PaddedSAM(samFilter)
+
+for read in paddedSAM.queries(rcSuffix=args.rcSuffix, rcNeeded=args.rcNeeded):
+    print(read.toString('fasta'), end='')
+
+if args.listReferenceInsertions:
+    if paddedSAM.referenceInsertions:
+        print('(0-based) insertions into the reference:\n%s' %
+              nucleotidesToStr(paddedSAM.referenceInsertions, '  '),
+              file=sys.stderr)
     else:
-        try:
-            for read in paddedSAM.queries(
-                    referenceName=args.referenceName,
-                    minLength=args.minLength,
-                    rcSuffix=args.rcSuffix,
-                    dropSecondary=args.dropSecondary,
-                    dropSupplementary=args.dropSupplementary,
-                    dropDuplicates=args.dropDuplicates,
-                    allowDuplicateIds=args.allowDuplicateIds,
-                    keepQCFailures=args.keepQCFailures,
-                    rcNeeded=args.rcNeeded):
-                print(read.toString('fasta'), end='')
-        except UnequalReferenceLengthError as e:
-            raise ValueError(
-                str(e) + ' So it is not clear how long the padded output '
-                'FASTA sequences should be. Use --referenceName to specify '
-                'which reference sequence is the one whose aligned reads you '
-                'want printed. Use --listReferenceNames to see a list of '
-                'reference sequence names and lengths.')
-        except UnknownReference as e:
-            raise ValueError(
-                str(e) + ' Use --listReferenceNames to see a list of '
-                'reference sequence names.')
-
-        if args.listReferenceInsertions:
-            if paddedSAM.referenceInsertions:
-                print('(0-based) insertions into the reference:\n%s' %
-                      nucleotidesToStr(paddedSAM.referenceInsertions, '  '),
-                      file=sys.stderr)
-            else:
-                print('No matches required an insertion into the reference.',
-                      file=sys.stderr)
-finally:
-    paddedSAM.close()
+        print('No matches required an insertion into the reference.',
+              file=sys.stderr)
diff --git a/dark/__init__.py b/dark/__init__.py
index 68c249c..b490eac 100644
--- a/dark/__init__.py
+++ b/dark/__init__.py
@@ -7,4 +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.33'
+__version__ = '3.0.34'
diff --git a/dark/filter.py b/dark/filter.py
index cd42063..05d36e8 100644
--- a/dark/filter.py
+++ b/dark/filter.py
@@ -169,15 +169,15 @@ class ReadSetFilter(object):
             return []
 
 
-def addFASTAFilteringCommandLineOptions(parser, filteringSAM=False):
+def addFASTAFilteringCommandLineOptions(parser):
     """
     Add standard FASTA filtering command-line options to an argparse parser.
 
+    These are options that can be used to select or omit entire FASTA records,
+    NOT options that change them (for that see
+    addFASTAEditingCommandLineOptions).
+
     @param parser: An C{argparse.ArgumentParser} instance.
-    @param filteringSAM: If C{True} only add options that can be used when
-        the result of filtering will be used to filter a SAM file. Because
-        we cannot yet filter out (for example) certain sites from a SAM file
-        we do not support some filtering options for SAM.
     """
     parser.add_argument(
         '--minLength', type=int,
@@ -267,27 +267,50 @@ def addFASTAFilteringCommandLineOptions(parser, filteringSAM=False):
         help=('A file of (1-based) sequence numbers to retain. Numbers must '
               'be one per line.'))
 
-    parser.add_argument(
-        '--idLambda', metavar='LAMBDA-FUNCTION',
-        help=('A one-argument function taking and returning a read id. '
-              'E.g., --idLambda "lambda id: id.split(\'_\')[0]" or '
-              '--idLambda "lambda id: id[:10]". If the function returns None, '
-              'the read will be filtered out.'))
 
-    parser.add_argument(
-        '--readLambda', metavar='LAMBDA-FUNCTION',
-        help=('A one-argument function taking and returning a read. '
-              'E.g., --readLambda "lambda r: Read(r.id.split(\'_\')[0], '
-              'r.sequence.strip(\'-\')". Make sure to also modify the quality '
-              'string if you change the length of a FASTQ sequence. If the '
-              'function returns None, the read will be filtered out. The '
-              'function will be passed to eval with the dark.reads classes '
-              'Read, DNARead, AARead, etc. all in scope.'))
+def parseFASTAFilteringCommandLineOptions(args, reads):
+    """
+    Examine parsed FASTA filtering command-line options and return filtered
+    reads.
+
+    @param args: An argparse namespace, as returned by the argparse
+        C{parse_args} function.
+    @param reads: A C{Reads} instance to filter.
+    @return: The filtered C{Reads} instance.
+    """
+    keepSequences = (
+        parseRangeString(args.keepSequences, convertToZeroBased=True)
+        if args.keepSequences else None)
+
+    removeSequences = (
+        parseRangeString(args.removeSequences, convertToZeroBased=True)
+        if args.removeSequences else None)
+
+    return reads.filter(
+        minLength=args.minLength, maxLength=args.maxLength,
+        whitelist=set(args.whitelist) if args.whitelist else None,
+        blacklist=set(args.blacklist) if args.blacklist else None,
+        whitelistFile=args.whitelistFile, blacklistFile=args.blacklistFile,
+        titleRegex=args.titleRegex,
+        negativeTitleRegex=args.negativeTitleRegex,
+        keepSequences=keepSequences, removeSequences=removeSequences,
+        head=args.head, removeDuplicates=args.removeDuplicates,
+        removeDuplicatesById=args.removeDuplicatesById,
+        randomSubset=args.randomSubset, trueLength=args.trueLength,
+        sampleFraction=args.sampleFraction,
+        sequenceNumbersFile=args.sequenceNumbersFile)
+
+
+def addFASTAEditingCommandLineOptions(parser):
+    """
+    Add standard FASTA editing command-line options to an argparse parser.
 
-    # Options below this point cannot be used to filter SAM.
-    if filteringSAM:
-        return
+    These are options that can be used to alter FASTA records, NOT options
+    that simply select or reject those things (for those see
+    addFASTAFilteringCommandLineOptions).
 
+    @param parser: An C{argparse.ArgumentParser} instance.
+    """
     # A mutually exclusive group for --keepSites, --keepSitesFile,
     # --removeSites, and --removeSitesFile.
     group = parser.add_mutually_exclusive_group()
@@ -339,83 +362,72 @@ def addFASTAFilteringCommandLineOptions(parser, filteringSAM=False):
               'description is the part of a sequence id after the '
               'first whitespace (if any).'))
 
+    parser.add_argument(
+        '--idLambda', metavar='LAMBDA-FUNCTION',
+        help=('A one-argument function taking and returning a read id. '
+              'E.g., --idLambda "lambda id: id.split(\'_\')[0]" or '
+              '--idLambda "lambda id: id[:10]". If the function returns None, '
+              'the read will be filtered out.'))
+
+    parser.add_argument(
+        '--readLambda', metavar='LAMBDA-FUNCTION',
+        help=('A one-argument function taking and returning a read. '
+              'E.g., --readLambda "lambda r: Read(r.id.split(\'_\')[0], '
+              'r.sequence.strip(\'-\')". Make sure to also modify the quality '
+              'string if you change the length of a FASTQ sequence. If the '
+              'function returns None, the read will be filtered out. The '
+              'function will be passed to eval with the dark.reads classes '
+              'Read, DNARead, AARead, etc. all in scope.'))
+
 
-def parseFASTAFilteringCommandLineOptions(args, reads, filteringSAM=False):
+def parseFASTAEditingCommandLineOptions(args, reads):
     """
-    Examine parsed command-line options and return information about kept
-    sites and sequences.
+    Examine parsed FASTA editing command-line options and return information
+    about kept sites and sequences.
 
     @param args: An argparse namespace, as returned by the argparse
         C{parse_args} function.
     @param reads: A C{Reads} instance to filter.
-    @param filteringSAM: If C{True} only examine options that can be used when
-        the result of filtering will be used to filter a SAM file. Because
-        we cannot yet filter out (for example) certain sites from a SAM file
-        we do not support some filtering options for SAM.
     @return: The filtered C{Reads} instance.
     """
-    keepSequences = (
-        parseRangeString(args.keepSequences, convertToZeroBased=True)
-        if args.keepSequences else None)
-
-    removeSequences = (
-        parseRangeString(args.removeSequences, convertToZeroBased=True)
-        if args.removeSequences else None)
-
-    if filteringSAM:
-        removeGaps = removeDescriptions = False
-        truncateTitlesAfter = keepSites = removeSites = None
-    else:
-        removeGaps = args.removeGaps
-        removeDescriptions = args.removeDescriptions
-        truncateTitlesAfter = args.truncateTitlesAfter
-        keepSites = (
-            parseRangeString(args.keepSites, convertToZeroBased=True)
-            if args.keepSites else None)
-
-        if args.keepSitesFile:
-            keepSites = keepSites or set()
-            with open(args.keepSitesFile) as fp:
-                for lineNumber, line in enumerate(fp):
-                    try:
-                        keepSites.update(
-                            parseRangeString(line, convertToZeroBased=True))
-                    except ValueError as e:
-                        raise ValueError(
-                            'Keep sites file %r line %d could not be parsed: '
-                            '%s' % (args.keepSitesFile, lineNumber, e))
-
-        removeSites = (
-            parseRangeString(args.removeSites, convertToZeroBased=True)
-            if args.removeSites else None)
-
-        if args.removeSitesFile:
-            removeSites = removeSites or set()
-            with open(args.removeSitesFile) as fp:
-                for lineNumber, line in enumerate(fp):
-                    try:
-                        removeSites.update(
-                            parseRangeString(line, convertToZeroBased=True))
-                    except ValueError as e:
-                        raise ValueError(
-                            'Remove sites file %r line %d parse error: %s'
-                            % (args.removeSitesFile, lineNumber, e))
+    removeGaps = args.removeGaps
+    removeDescriptions = args.removeDescriptions
+    truncateTitlesAfter = args.truncateTitlesAfter
+    keepSites = (
+        parseRangeString(args.keepSites, convertToZeroBased=True)
+        if args.keepSites else None)
+
+    if args.keepSitesFile:
+        keepSites = keepSites or set()
+        with open(args.keepSitesFile) as fp:
+            for lineNumber, line in enumerate(fp):
+                try:
+                    keepSites.update(
+                        parseRangeString(line, convertToZeroBased=True))
+                except ValueError as e:
+                    raise ValueError(
+                        'Keep sites file %r line %d could not be parsed: '
+                        '%s' % (args.keepSitesFile, lineNumber, e))
+
+    removeSites = (
+        parseRangeString(args.removeSites, convertToZeroBased=True)
+        if args.removeSites else None)
+
+    if args.removeSitesFile:
+        removeSites = removeSites or set()
+        with open(args.removeSitesFile) as fp:
+            for lineNumber, line in enumerate(fp):
+                try:
+                    removeSites.update(
+                        parseRangeString(line, convertToZeroBased=True))
+                except ValueError as e:
+                    raise ValueError(
+                        'Remove sites file %r line %d parse error: %s'
+                        % (args.removeSitesFile, lineNumber, e))
 
     return reads.filter(
-        minLength=args.minLength, maxLength=args.maxLength,
         removeGaps=removeGaps,
-        whitelist=set(args.whitelist) if args.whitelist else None,
-        blacklist=set(args.blacklist) if args.blacklist else None,
-        whitelistFile=args.whitelistFile, blacklistFile=args.blacklistFile,
-        titleRegex=args.titleRegex,
-        negativeTitleRegex=args.negativeTitleRegex,
         truncateTitlesAfter=truncateTitlesAfter,
-        keepSequences=keepSequences, removeSequences=removeSequences,
-        head=args.head, removeDuplicates=args.removeDuplicates,
-        removeDuplicatesById=args.removeDuplicatesById,
         removeDescriptions=removeDescriptions,
-        randomSubset=args.randomSubset, trueLength=args.trueLength,
-        sampleFraction=args.sampleFraction,
-        sequenceNumbersFile=args.sequenceNumbersFile,
         idLambda=args.idLambda, readLambda=args.readLambda,
         keepSites=keepSites, removeSites=removeSites)
diff --git a/dark/sam.py b/dark/sam.py
index 4ea07a1..eef7854 100644
--- a/dark/sam.py
+++ b/dark/sam.py
@@ -1,3 +1,6 @@
+import six
+
+from itertools import chain
 from contextlib import contextmanager
 from collections import Counter, defaultdict
 
@@ -9,93 +12,319 @@ from dark.reads import Read, DNARead
 
 
 class UnequalReferenceLengthError(Exception):
-    "The reference sequences in a SAM/BAM file are not all of the same length."
+    "The references of interest in a SAM/BAM file are not of the same length."
 
 
 class UnknownReference(Exception):
     "Reference sequence not found in SAM/BAM file."
 
 
+class InvalidSAM(Exception):
+    "SAM/BAM file has unexpected/invalid content."
+
+
 # From https://samtools.github.io/hts-specs/SAMv1.pdf
 _CONSUMES_QUERY = {CMATCH, CINS, CSOFT_CLIP, CEQUAL, CDIFF}
 _CONSUMES_REFERENCE = {CMATCH, CDEL, CREF_SKIP, CEQUAL, CDIFF}
 
 
+@contextmanager
+def samfile(filename):
+    """
+    A context manager to open and close a SAM/BAM file.
+
+    @param filename: A C{str} file name to open.
+    """
+    f = AlignmentFile(filename)
+    yield f
+    f.close()
+
+
+def samReferencesToStr(filenameOrSamfile, indent=0):
+    """
+    List SAM/BAM file reference names and lengths.
+
+    @param filenameOrSamfile: Either a C{str} SAM/BAM file name or an
+        instance of C{pysam.AlignmentFile}.
+    @param indent: An C{int} number of spaces to indent each line.
+    @return: A C{str} describing known reference names and their lengths.
+    """
+    indent = ' ' * indent
+
+    def _references(sam):
+        result = []
+        for i in range(sam.nreferences):
+            result.append('%s%s (length %d)' % (
+                indent, sam.get_reference_name(i), sam.lengths[i]))
+        return '\n'.join(result)
+
+    if isinstance(filenameOrSamfile, six.string_types):
+        with samfile(filenameOrSamfile) as sam:
+            return _references(sam)
+    else:
+        return _references(sam)
+
+
+class SAMFilter(object):
+    """
+    Filter a SAM/BAM file.
+
+    @param filename: The C{str} name of a BAM/SAM file to filter.
+    @param filterRead: A function that takes a C{dark.reads.Read) instance
+        and returns either C{None} or a C{Read} instance, according to
+        whether the passed read should be omitted or not. If C{None} is
+        passed, no filtering on reads (i.e., queries) is done.
+    @param referenceIds: Either C{None} or a set of C{str} reference ids
+        that should be kept (other references will be dropped).
+    @param storeQueryIds: If C{True}, query ids will be stored (in
+        C{self.queryIds}) as the SAM/BAM file is read.
+    @param dropUnmapped: If C{True}, unmapped matches will be excluded.
+    @param dropSecondary: If C{True}, secondary matches will be excluded.
+    @param dropSupplementary: If C{True}, supplementary matches will be
+        excluded.
+    @param dropDuplicates: If C{True}, matches flagged as optical or PCR
+        duplicates will be excluded.
+    @param keepQCFailures: If C{True}, reads that are considered quality
+        control failures will be included.
+    @param minScore: If not C{None}, alignments with score tag values less
+        than this value will not be output. If given, alignments that do not
+        have a score will not be output.
+    @param maxScore: If not C{None}, alignments with score tag values greater
+        than this value will not be output. If given, alignments that do not
+        have a score will not be output.
+    @param scoreTag: The alignment tag to extract for minScore and maxScore
+        comparisons.
+    """
+    def __init__(self, filename, filterRead=None, referenceIds=None,
+                 storeQueryIds=False, dropUnmapped=False,
+                 dropSecondary=False, dropSupplementary=False,
+                 dropDuplicates=False, keepQCFailures=False, minScore=None,
+                 maxScore=None, scoreTag='AS'):
+        self.filename = filename
+        self.filterRead = filterRead
+        self.referenceIds = referenceIds
+        self.dropUnmapped = dropUnmapped
+        self.dropSecondary = dropSecondary
+        self.dropSupplementary = dropSupplementary
+        self.dropDuplicates = dropDuplicates
+        self.keepQCFailures = keepQCFailures
+        self.storeQueryIds = storeQueryIds
+        self.minScore = minScore
+        self.maxScore = maxScore
+        self.scoreTag = scoreTag
+
+    @staticmethod
+    def addFilteringOptions(parser, samfileIsPositionalArg=False):
+        """
+        Add options to an argument parser for filtering SAM/BAM.
+
+        @param samfileIsPositionalArg: If C{True} the SAM/BAM file must
+            be given as the final argument on the command line (without
+            being preceded by --samfile).
+        @param parser: An C{argparse.ArgumentParser} instance.
+        """
+        parser.add_argument(
+            '%ssamfile' % ('' if samfileIsPositionalArg else '--'),
+            required=True,
+            help='The SAM/BAM file to filter.')
+
+        parser.add_argument(
+            '--referenceId', metavar='ID', nargs='+', action='append',
+            help=('A reference sequence id whose alignments should be kept '
+                  '(alignments against other references will be dropped). '
+                  'If omitted, alignments against all references will be '
+                  'kept. May be repeated.'))
+
+        parser.add_argument(
+            '--dropUnmapped', default=False, action='store_true',
+            help='If given, unmapped matches will not be output.')
+
+        parser.add_argument(
+            '--dropSecondary', default=False, action='store_true',
+            help='If given, secondary matches will not be output.')
+
+        parser.add_argument(
+            '--dropSupplementary', default=False, action='store_true',
+            help='If given, supplementary matches will not be output.')
+
+        parser.add_argument(
+            '--dropDuplicates', default=False, action='store_true',
+            help=('If given, matches flagged as optical or PCR duplicates '
+                  'will not be output.'))
+
+        parser.add_argument(
+            '--keepQCFailures', default=False, action='store_true',
+            help=('If given, reads that are considered quality control '
+                  'failures will be included in the output.'))
+
+        parser.add_argument(
+            '--minScore', type=float,
+            help=('If given, alignments with --scoreTag (default AS) values '
+                  'less than this value will not be output. If given, '
+                  'alignments that do not have a score will not be output.'))
+
+        parser.add_argument(
+            '--maxScore', type=float,
+            help=('If given, alignments with --scoreTag (default AS) values '
+                  'greater than this value will not be output. If given, '
+                  'alignments that do not have a score will not be output.'))
+
+        parser.add_argument(
+            '--scoreTag', default='AS',
+            help=('The alignment tag to extract for --minScore and --maxScore '
+                  'comparisons.'))
+
+    @classmethod
+    def parseFilteringOptions(cls, args, filterRead=None, storeQueryIds=False):
+        """
+        Parse command line options (added in C{addSAMFilteringOptions}.
+
+        @param args: The command line arguments, as returned by
+            C{argparse.parse_args}.
+        @param filterRead: A one-argument function that accepts a read
+            and returns C{None} if the read should be omitted in filtering
+            or else a C{Read} instance.
+        @param storeQueryIds: If C{True}, query ids will be stored as the
+            SAM/BAM file is read.
+        @return: A C{SAMFilter} instance.
+        """
+        referenceIds = (set(chain.from_iterable(args.referenceId))
+                        if args.referenceId else None)
+
+        return cls(
+            args.samfile,
+            filterRead=filterRead,
+            referenceIds=referenceIds,
+            storeQueryIds=storeQueryIds,
+            dropUnmapped=args.dropUnmapped,
+            dropSecondary=args.dropSecondary,
+            dropSupplementary=args.dropSupplementary,
+            dropDuplicates=args.dropDuplicates,
+            keepQCFailures=args.keepQCFailures,
+            minScore=args.minScore,
+            maxScore=args.maxScore)
+
+    def alignments(self):
+        """
+        Get alignments from the SAM/BAM file, subject to filtering.
+        """
+        referenceIds = self.referenceIds
+        dropUnmapped = self.dropUnmapped
+        dropSecondary = self.dropSecondary
+        dropSupplementary = self.dropSupplementary
+        dropDuplicates = self.dropDuplicates
+        keepQCFailures = self.keepQCFailures
+        storeQueryIds = self.storeQueryIds
+        filterRead = self.filterRead
+        minScore = self.minScore
+        maxScore = self.maxScore
+        scoreTag = self.scoreTag
+
+        if storeQueryIds:
+            self.queryIds = queryIds = set()
+
+        count = 0
+        with samfile(self.filename) as samAlignment:
+            for count, alignment in enumerate(samAlignment.fetch(), start=1):
+                if storeQueryIds:
+                    queryIds.add(alignment.query_name)
+
+                if minScore is not None or maxScore is not None:
+                    try:
+                        score = alignment.get_tag(scoreTag)
+                    except KeyError:
+                        continue
+                    else:
+                        if ((minScore is not None and score < minScore) or
+                                (maxScore is not None and score > maxScore)):
+                            continue
+
+                if ((filterRead is None or
+                     filterRead(Read(alignment.query_name,
+                                     alignment.query_sequence,
+                                     alignment.qual))) and
+                    not (
+                        (referenceIds and
+                         alignment.reference_name not in referenceIds) or
+                        (alignment.is_unmapped and dropUnmapped) or
+                        (alignment.is_secondary and dropSecondary) or
+                        (alignment.is_supplementary and dropSupplementary) or
+                        (alignment.is_duplicate and dropDuplicates) or
+                        (alignment.is_qcfail and not keepQCFailures))):
+                    yield alignment
+
+        self.alignmentCount = count
+
+    def referenceLengths(self):
+        """
+        Get the lengths of wanted references.
+
+        @raise UnknownReference: If a reference id is not present in the
+            SAM/BAM file.
+        @return: A C{dict} of C{str} reference id to C{int} length with a key
+            for each reference id in C{self.referenceIds} or for all references
+            if C{self.referenceIds} is C{None}.
+        """
+        result = {}
+        with samfile(self.filename) as sam:
+            if self.referenceIds:
+                for referenceId in self.referenceIds:
+                    tid = sam.get_tid(referenceId)
+                    if tid == -1:
+                        raise UnknownReference(
+                            'Reference %r is not present in the SAM/BAM file.'
+                            % referenceId)
+                    else:
+                        result[referenceId] = sam.lengths[tid]
+            else:
+                result = dict(zip(sam.references, sam.lengths))
+
+        return result
+
+
 class PaddedSAM(object):
     """
     Obtain aligned (padded) queries from a SAM/BAM file.
 
-    @param filename: The C{str} name of the SAM/BAM file.
+    @param samFilter: A C{SAMFilter} instance.
+    @raises UnequalReferenceLengthError: If C{referenceName} is C{None}
+        and the reference sequence lengths in the SAM/BAM file are not all
+        identical.
+    @raises UnknownReference: If C{referenceName} does not exist.
     """
-    def __init__(self, filename):
-        self.samfile = AlignmentFile(filename)
+    def __init__(self, samFilter):
+        referenceLengths = samFilter.referenceLengths()
+
+        if len(set(referenceLengths.values())) != 1:
+            raise UnequalReferenceLengthError(
+                'Your %d SAM/BAM file reference sequence '
+                'lengths (%s) are not all identical.' % (
+                    len(referenceLengths),
+                    ', '.join(
+                        '%s=%d' % (id_, referenceLengths[id_])
+                        for id_ in sorted(referenceLengths))))
+
+        # Get the length of any of the sequences (they are all identical).
+        self.referenceLength = referenceLengths.popitem()[1]
+        self.samFilter = samFilter
         # self.referenceInsertions will be keyed by query id (the query
         # that would cause a reference insertion). The values will be lists
         # of 2-tuples, with each 2-tuple containing an offset into the
         # reference sequence and the C{str} of nucleotides that would be
         # inserted starting at that offset.
         self.referenceInsertions = defaultdict(list)
-        self.alignmentCount = 0
-
-    def close(self):
-        """
-        Close the opened SAM/BAM file.
-        """
-        self.samfile.close()
-
-    def referencesToStr(self, indent=0):
-        """
-        List the reference names and their lengths.
 
-        @param indent: An C{int} number of spaces to indent each line.
-        @return: A C{str} describing known reference names and their lengths.
-        """
-        samfile = self.samfile
-        result = []
-        indent = ' ' * indent
-        for i in range(samfile.nreferences):
-            result.append('%s%s (length %d)' % (
-                indent, samfile.get_reference_name(i), samfile.lengths[i]))
-        return '\n'.join(result)
-
-    def queries(self, referenceName=None, minLength=0, rcSuffix='',
-                dropSecondary=False, dropSupplementary=False,
-                dropDuplicates=False, allowDuplicateIds=False,
-                keepQCFailures=False, rcNeeded=False, padChar='-',
-                queryInsertionChar='N', addAlignment=False,
-                storeQueryIds=False):
+    def queries(self, rcSuffix='', rcNeeded=False, padChar='-',
+                queryInsertionChar='N', allowDuplicateIds=False,
+                addAlignment=False):
         """
         Produce padded (with gaps) queries according to the CIGAR string and
         reference sequence length for each matching query sequence.
 
-        @param referenceName: The C{str} name of the reference sequence to
-            print alignments for. This is only needed if the SAM/BAM alignment
-            was made against multiple references *and* they are not all of the
-            same length. If there is only one reference sequence or if all
-            reference sequences are of the same length, there is no need to
-            provide a reference name (i.e., pass C{None}).
-        @param minLength: Ignore queries shorter than this C{int} value. Note
-            that this refers to the length of the query sequence once it has
-            been aligned to the reference. The alignment may introduce
-            C{queryInsertionChar} characters into the read, and these are
-            counted towards its length because the alignment is assuming the
-            query is missing a base at those locations.
         @param rcSuffix: A C{str} to add to the end of query names that are
             reverse complemented. This is added before the /1, /2, etc., that
             are added for duplicated ids (if there are duplicates and
             C{allowDuplicateIds} is C{False}.
-        @param dropSecondary: If C{True}, secondary matches will not be
-            yielded.
-        @param dropSupplementary: If C{True}, supplementary matches will not be
-            yielded.
-        @param dropDuplicates: If C{True}, matches flagged as optical or PCR
-            duplicates will not be yielded.
-        @param allowDuplicateIds: If C{True}, repeated query ids (due to
-            secondary or supplemental matches) will not have /1, /2, etc.
-            appended to their ids. So repeated ids may appear in the yielded
-            FASTA.
-        @param keepQCFailures: If C{True}, reads that are marked as quality
-            control failures will be included in the output.
         @param rcNeeded: If C{True}, queries that are flagged as matching when
             reverse complemented should have reverse complementing when
             preparing the output sequences. This must be used if the program
@@ -109,39 +338,22 @@ class PaddedSAM(object):
             is inserted as a 'missing' query character (i.e., a base that can
             be assumed to have been lost due to an error) whose existence is
             necessary for the match to continue.
+        @param allowDuplicateIds: If C{True}, repeated query ids (due to
+            secondary or supplemental matches) will not have /1, /2, etc.
+            appended to their ids. So repeated ids may appear in the yielded
+            FASTA.
         @param addAlignment: If C{True} the reads yielded by the returned
             generator will also have an C{alignment} attribute, being the
-            C{pysam.AlignedSegment} for each query.
-        @param storeQueryIds: If C{True}, query ids will be stored in a
-            C{set} attribute called C{queryIds}.
-        @raises UnequalReferenceLengthError: If C{referenceName} is C{None}
-            and the reference sequence lengths in the SAM/BAM file are not all
-            identical.
-        @raises UnknownReference: If C{referenceName} does not exist.
+            C{pysam.AlignedSegment} for the query.
+        @raises InvalidSAM: If a query has an empty SEQ field and either there
+            is no previous alignment or the alignment is not marked as
+            secondary or supplementary.
         @return: A generator that yields C{Read} instances that are padded
             with gap characters to align them to the length of the reference
             sequence. See C{addAlignment}, above, to yield reads with the
             corresponding C{pysam.AlignedSegment}.
         """
-        samfile = self.samfile
-
-        if referenceName:
-            referenceId = samfile.get_tid(referenceName)
-            if referenceId == -1:
-                raise UnknownReference(
-                    'Reference %r is not present in the SAM/BAM file.'
-                    % referenceName)
-            referenceLength = samfile.lengths[referenceId]
-        else:
-            # No reference given. All references must have the same length.
-            if len(set(samfile.lengths)) != 1:
-                raise UnequalReferenceLengthError(
-                    'Your SAM/BAM file has %d reference sequences, and their '
-                    'lengths (%s) are not all identical.' % (
-                        samfile.nreferences,
-                        ', '.join(map(str, sorted(samfile.lengths)))))
-            referenceId = None
-            referenceLength = samfile.lengths[0]
+        referenceLength = self.referenceLength
 
         # Hold the count for each id so we can add /1, /2 etc to duplicate
         # ids (unless --allowDuplicateIds was given).
@@ -149,23 +361,9 @@ class PaddedSAM(object):
 
         MATCH_OPERATIONS = {CMATCH, CEQUAL, CDIFF}
         lastQuery = None
-        if storeQueryIds:
-            self.queryIds = queryIds = set()
 
-        for lineNumber, alignment in enumerate(samfile.fetch(), start=1):
-            self.alignmentCount += 1
-
-            if storeQueryIds:
-                queryIds.add(alignment.query_name)
-
-            if (alignment.is_unmapped or
-                    (alignment.is_secondary and dropSecondary) or
-                    (alignment.is_supplementary and dropSupplementary) or
-                    (alignment.is_duplicate and dropDuplicates) or
-                    (alignment.is_qcfail and not keepQCFailures) or
-                    (referenceId is not None and
-                     alignment.reference_id != referenceId)):
-                continue
+        for lineNumber, alignment in enumerate(
+                self.samFilter.alignments(), start=1):
 
             query = alignment.query_sequence
 
@@ -176,13 +374,13 @@ class PaddedSAM(object):
             if query is None:
                 if alignment.is_secondary or alignment.is_supplementary:
                     if lastQuery is None:
-                        raise ValueError(
+                        raise InvalidSAM(
                             'Query line %d has an empty SEQ field, but no '
                             'previous alignment is present.' % lineNumber)
                     else:
                         query = lastQuery
                 else:
-                    raise ValueError(
+                    raise InvalidSAM(
                         'Query line %d has an empty SEQ field, but the '
                         'alignment is not marked as secondary or '
                         'supplementary.' % lineNumber)
@@ -197,8 +395,8 @@ class PaddedSAM(object):
                 if rcSuffix:
                     alignment.query_name += rcSuffix
 
-            # Adjust the query id if it's a duplicate and we're not
-            # allowing duplicates.
+            # Adjust the query id if it's a duplicate and we're not allowing
+            # duplicates.
             if allowDuplicateIds:
                 queryId = alignment.query_name
             else:
@@ -301,14 +499,9 @@ class PaddedSAM(object):
             assert queryIndex == len(query)
 
             # We cannot test we consumed the entire reference.  The CIGAR
-            # string applies to (i.e., exhausts) the query and is silent about
-            # the part of the reference that to the right of the aligned query.
-
-            # Check the length restriction now that we have (possibly) added
-            # queryInsertionChar characters to pad the query out to the length
-            # it requires to match the reference.
-            if len(alignedSequence) < minLength:
-                continue
+            # string applies to (and exhausts) the query but is silent
+            # about the part of the reference that lies to the right of the
+            # aligned query.
 
             # Put gap characters before and after the aligned sequence so that
             # it is offset properly and matches the length of the reference.
@@ -323,15 +516,3 @@ class PaddedSAM(object):
                 read.alignment = alignment
 
             yield read
-
-
-@contextmanager
-def samfile(filename):
-    """
-    A context manager to open and close a SAM/BAM file.
-
-    @param filename: A C{str} file name to open.
-    """
-    f = AlignmentFile(filename)
-    yield f
-    f.close()
diff --git a/dark/utils.py b/dark/utils.py
index 98912eb..4862cba 100644
--- a/dark/utils.py
+++ b/dark/utils.py
@@ -185,3 +185,32 @@ else:
         def __exit__(self, exc_type, exc_val, exc_tb):
             self.close()
             return False
+
+
+def baseCountsToStr(counts):
+    """
+    Convert base counts to a string.
+
+    @param counts: A C{Counter} instance.
+    @return: A C{str} representation of nucleotide counts at an offset.
+    """
+    return ' '.join([
+        ('%s:%d' % (base, counts[base])) for base in sorted(counts)])
+
+
+def nucleotidesToStr(nucleotides, prefix=''):
+    """
+    Convert offsets and base counts to a string.
+
+    @param nucleotides: A C{defaultdict(Counter)} instance, keyed
+        by C{int} offset, with nucleotides keying the Counters.
+    @param prefix: A C{str} to put at the start of each line.
+    @return: A C{str} representation of the offsets and nucleotide
+        counts for each.
+    """
+    result = []
+    for offset in sorted(nucleotides):
+        result.append(
+            '%s%d: %s' % (prefix, offset,
+                          baseCountsToStr(nucleotides[offset])))
+    return '\n'.join(result)
diff --git a/setup.py b/setup.py
index 1bcc868..baf9d30 100644
--- a/setup.py
+++ b/setup.py
@@ -77,6 +77,7 @@ scripts = [
     'bin/read-blast-xml.py',
     'bin/sam-to-fasta-alignment.py',
     'bin/sam-reference-read-counts.py',
+    'bin/sam-references.py',
     'bin/sff-to-fastq.py',
     'bin/split-fasta-by-adaptors.py',
     'bin/summarize-fasta-bases.py',


Test patch:
diff --git a/test/test_sam.py b/test/test_sam.py
index 6900a13..28c640b 100644
--- a/test/test_sam.py
+++ b/test/test_sam.py
@@ -4,8 +4,10 @@ from tempfile import mkstemp
 from os import close, unlink, write
 from contextlib import contextmanager
 
-from dark.reads import Read
-from dark.sam import PaddedSAM, UnequalReferenceLengthError, UnknownReference
+from dark.reads import Read, ReadFilter
+from dark.sam import (
+    PaddedSAM, SAMFilter, UnequalReferenceLengthError, UnknownReference,
+    InvalidSAM, samReferencesToStr)
 
 
 # These tests actually use the filesystem to read files. That's due to the API
@@ -24,55 +26,241 @@ def dataFile(data):
     unlink(filename)
 
 
-class TestPaddedSAM(TestCase):
+class TestSAMFilter(TestCase):
     """
-    Test the PaddedSAM class.
+    Test the SAMFilter class.
     """
+    def testUnknownReferences(self):
+        """
+        Passing an unknown reference id to the referenceLengths method must
+        result in an UnknownReference exception.
+        """
+        data = '\n'.join([
+            '@SQ SN:id1 LN:90',
+            '@SQ SN:id2 LN:90',
+        ]).replace(' ', '\t')
 
-    # In reading the tests below, it is important to remember that the start
-    # position (in the reference) of the match in SAM format is 1-based.  This
-    # is the 4th field in the non-header SAM lines (i.e., those that don't
-    # start with @). If you look at the code in ../dark/sam.py, pysam provides
-    # a 'reference_start' attribute that is 0-based.
+        with dataFile(data) as filename:
+            sam = SAMFilter(filename, referenceIds={'unknown'})
+            error = ("^Reference 'unknown' is not present in the "
+                     "SAM/BAM file\\.$")
+            assertRaisesRegex(self, UnknownReference, error,
+                              sam.referenceLengths)
 
-    def testReferencesToStr(self):
+    def testStoreQueryIds(self):
         """
-        The referencesToStr method must return the expected string.
+        If we request that query ids are saved, they must be.
         """
         data = '\n'.join([
-            '@SQ SN:id1 LN:90',
-            '@SQ SN:id2 LN:91',
+            '@SQ SN:ref1 LN:10',
+            'query1 0 ref1 2 60 2=2X2M * 0 0 TCTAGG 123456',
+            'query2 0 ref1 2 60 2= * 0 0 TC XY',
+            'query2 0 ref1 2 60 2= * 0 0 TC XY',
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
-            self.assertEqual('id1 (length 90)\nid2 (length 91)',
-                             ps.referencesToStr())
-            ps.close()
+            sf = SAMFilter(filename, storeQueryIds=True)
+            list(sf.alignments())
+            self.assertEqual({'query1', 'query2'}, sf.queryIds)
 
-    def testUnknownReferences(self):
+    def testAlignmentCount(self):
         """
-        Passing an unknown reference name to 'queries' must result in an
-        UnknownReference exception.
+        When all queries have been yielded, the alignment count must be
+        as expected.
         """
         data = '\n'.join([
-            '@SQ SN:id1 LN:90',
-            '@SQ SN:id2 LN:91',
+            '@SQ SN:ref1 LN:10',
+            'query1 0 ref1 2 60 2=2X2M * 0 0 TCTAGG 123456',
+            'query2 0 ref1 2 60 2= * 0 0 TC XY',
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
-            error = ("^Reference 'unknown' is not present in the "
-                     "SAM/BAM file\\.$")
-            queries = ps.queries(referenceName='unknown')
-            assertRaisesRegex(self, UnknownReference, error, list, queries)
-            ps.close()
+            sf = SAMFilter(filename)
+            list(sf.alignments())
+            self.assertEqual(2, sf.alignmentCount)
+
+    def testMinLength(self):
+        """
+        A request for reads that are only longer than a certain value should
+        result in the expected result.
+        """
+        data = '\n'.join([
+            '@SQ SN:ref1 LN:10',
+            'query1 0 ref1 2 60 2=2X2M * 0 0 TCTAGG ZZZZZZ',
+            'query2 0 ref1 2 60 2= * 0 0 TC ZZ',
+        ]).replace(' ', '\t')
+
+        with dataFile(data) as filename:
+            filterRead = ReadFilter(minLength=6).filter
+            sf = SAMFilter(filename, filterRead=filterRead)
+            (alignment,) = list(sf.alignments())
+            self.assertEqual('query1', alignment.query_name)
+
+    def testDropSecondary(self):
+        """
+        Dropping matches flagged as secondary must give the expected result.
+        """
+        data = '\n'.join([
+            '@SQ SN:ref1 LN:10',
+            'query1 0 ref1 2 60 2=2X2M * 0 0 TCTAGG ZZZZZZ',
+            'query2 256 ref1 2 60 2= * 0 0 TC ZZ',
+        ]).replace(' ', '\t')
+
+        with dataFile(data) as filename:
+            sf = SAMFilter(filename, dropSecondary=True)
+            (alignment,) = list(sf.alignments())
+            self.assertEqual('query1', alignment.query_name)
+
+    def testDropSupplementary(self):
+        """
+        Dropping matches flagged as supplementary must give the expected
+        result.
+        """
+        data = '\n'.join([
+            '@SQ SN:ref1 LN:10',
+            'query1 0 ref1 2 60 2=2X2M * 0 0 TCTAGG ZZZZZZ',
+            'query2 2048 ref1 2 60 2= * 0 0 TC ZZ',
+        ]).replace(' ', '\t')
+
+        with dataFile(data) as filename:
+            sf = SAMFilter(filename, dropSupplementary=True)
+            (alignment,) = list(sf.alignments())
+            self.assertEqual('query1', alignment.query_name)
+
+    def testDropDuplicates(self):
+        """
+        Dropping matches flagged as optical or PCR duplicates must give the
+        expected result.
+        """
+        data = '\n'.join([
+            '@SQ SN:ref1 LN:10',
+            'query1 0 ref1 2 60 2=2X2M * 0 0 TCTAGG ZZZZZZ',
+            'query2 1024 ref1 2 60 2= * 0 0 TC ZZ',
+        ]).replace(' ', '\t')
+
+        with dataFile(data) as filename:
+            sf = SAMFilter(filename, dropDuplicates=True)
+            (alignment,) = list(sf.alignments())
+            self.assertEqual('query1', alignment.query_name)
+
+    def testKeepQualityControlFailures(self):
+        """
+        Keeping matches flagged as quality control failures must give the
+        expected result.
+        """
+        data = '\n'.join([
+            '@SQ SN:ref1 LN:10',
+            'query1 0 ref1 2 60 2=2X2M * 0 0 TCTAGG ZZZZZZ',
+            'query2 512 ref1 4 60 2= * 0 0 TC ZZ',
+        ]).replace(' ', '\t')
+
+        with dataFile(data) as filename:
+            sf = SAMFilter(filename, keepQCFailures=True)
+            (alignment1, alignment2) = list(sf.alignments())
+            self.assertEqual('query1', alignment1.query_name)
+            self.assertEqual('query2', alignment2.query_name)
+
+    def testMinScoreNoScores(self):
+        """
+        A request for reads with alignment scores no lower than a given value
+        must produce an empty result when no alignments have scores.
+        """
+        data = '\n'.join([
+            '@SQ SN:ref1 LN:10',
+            'query1 0 ref1 2 60 2=2X2M * 0 0 TCTAGG ZZZZZZ',
+            'query2 0 ref1 2 60 2= * 0 0 TC ZZ',
+        ]).replace(' ', '\t')
+
+        with dataFile(data) as filename:
+            sf = SAMFilter(filename, minScore=6)
+            self.assertEqual([], list(sf.alignments()))
+
+    def testMinScore(self):
+        """
+        A request for reads with alignment scores no lower than a given value
+        must produce the expected result when some alignments have scores.
+        """
+        data = '\n'.join([
+            '@SQ SN:ref1 LN:10',
+            'query1 0 ref1 2 60 2=2X2M * 0 0 TCTAGG ZZZZZZ AS:i:10',
+            'query2 0 ref1 2 60 2= * 0 0 TC ZZ',
+            'query3 0 ref1 2 60 2=2X2M * 0 0 TCTAGG ZZZZZZ AS:i:3',
+        ]).replace(' ', '\t')
+
+        with dataFile(data) as filename:
+            sf = SAMFilter(filename, minScore=6)
+            (alignment,) = list(sf.alignments())
+            self.assertEqual('query1', alignment.query_name)
+
+    def testMaxScoreNoScores(self):
+        """
+        A request for reads with alignment scores no higher than a given value
+        must produce an empty result when no alignments have scores.
+        """
+        data = '\n'.join([
+            '@SQ SN:ref1 LN:10',
+            'query1 0 ref1 2 60 2=2X2M * 0 0 TCTAGG ZZZZZZ',
+            'query2 0 ref1 2 60 2= * 0 0 TC ZZ',
+            'query3 0 ref1 2 60 2=2X2M * 0 0 TCTAGG ZZZZZZ',
+        ]).replace(' ', '\t')
+
+        with dataFile(data) as filename:
+            sf = SAMFilter(filename, maxScore=6)
+            self.assertEqual([], list(sf.alignments()))
+
+    def testMaxScore(self):
+        """
+        A request for reads with alignment scores no higher than a given value
+        must produce the expected result when some alignments have scores.
+        """
+        data = '\n'.join([
+            '@SQ SN:ref1 LN:10',
+            'query1 0 ref1 2 60 2=2X2M * 0 0 TCTAGG ZZZZZZ AS:i:10',
+            'query2 0 ref1 2 60 2= * 0 0 TC ZZ',
+            'query3 0 ref1 2 60 2=2X2M * 0 0 TCTAGG ZZZZZZ AS:i:3',
+        ]).replace(' ', '\t')
+
+        with dataFile(data) as filename:
+            sf = SAMFilter(filename, maxScore=6)
+            (alignment,) = list(sf.alignments())
+            self.assertEqual('query3', alignment.query_name)
+
+    def testMinAndMaxScore(self):
+        """
+        A request for reads with alignment scores no lower or higher than
+        given values must produce the expected result.
+        """
+        data = '\n'.join([
+            '@SQ SN:ref1 LN:10',
+            'query1 0 ref1 2 60 2=2X2M * 0 0 TCTAGG ZZZZZZ AS:i:10',
+            'query2 0 ref1 2 60 2=2X2M * 0 0 TCTAGG ZZZZZZ AS:i:12',
+            'query3 0 ref1 2 60 2= * 0 0 TC ZZ',
+            'query4 0 ref1 2 60 2=2X2M * 0 0 TCTAGG ZZZZZZ AS:i:3',
+            'query5 0 ref1 2 60 2=2X2M * 0 0 TCTAGG ZZZZZZ AS:i:2',
+        ]).replace(' ', '\t')
+
+        with dataFile(data) as filename:
+            sf = SAMFilter(filename, minScore=3, maxScore=10)
+            (alignment1, alignment2) = list(sf.alignments())
+            self.assertEqual('query1', alignment1.query_name)
+            self.assertEqual('query4', alignment2.query_name)
+
+
+class TestPaddedSAM(TestCase):
+    """
+    Test the PaddedSAM class.
+    """
+
+    # In reading the tests below, it is important to remember that the start
+    # position (in the reference) of the match in SAM format is 1-based.  This
+    # is the 4th field in the non-header SAM lines (i.e., those that don't
+    # start with @). If you look at the code in ../dark/sam.py, pysam provides
+    # a 'reference_start' attribute that is 0-based.
 
     def testUnequalReferenceLengths(self):
         """
-        Passing no reference name to 'queries' when the references have
-        different lengths must result in an UnequalReferenceLengthError
-        exception.
+        Passing no reference ids when the references have different lengths
+        must result in an UnequalReferenceLengthError exception.
         """
         data = '\n'.join([
             '@SQ SN:id1 LN:90',
@@ -80,13 +268,10 @@ class TestPaddedSAM(TestCase):
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
-            error = ('^Your SAM/BAM file has 2 reference sequences, and their '
-                     'lengths \(90, 91\) are not all identical\.$')
-            queries = ps.queries()
-            assertRaisesRegex(self, UnequalReferenceLengthError, error, list,
-                              queries)
-            ps.close()
+            error = ('^Your 2 SAM/BAM file reference sequence lengths '
+                     '\\(id1=90, id2=91\\) are not all identical\\.$')
+            assertRaisesRegex(self, UnequalReferenceLengthError, error,
+                              PaddedSAM, SAMFilter(filename))
 
     def testAllMMatch(self):
         """
@@ -98,10 +283,9 @@ class TestPaddedSAM(TestCase):
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
+            ps = PaddedSAM(SAMFilter(filename))
             (read,) = list(ps.queries())
             self.assertEqual(Read('query1', '-TCTAGG---'), read)
-            ps.close()
 
     def testMixedMatch(self):
         """
@@ -114,10 +298,9 @@ class TestPaddedSAM(TestCase):
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
+            ps = PaddedSAM(SAMFilter(filename))
             (read,) = list(ps.queries())
             self.assertEqual(Read('query1', '-TCTAGG---'), read)
-            ps.close()
 
     def testHardClipLeft(self):
         """
@@ -130,10 +313,9 @@ class TestPaddedSAM(TestCase):
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
+            ps = PaddedSAM(SAMFilter(filename))
             (read,) = list(ps.queries())
             self.assertEqual(Read('query1', '-TCTAGG---'), read)
-            ps.close()
 
     def testHardClipRight(self):
         """
@@ -146,10 +328,9 @@ class TestPaddedSAM(TestCase):
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
+            ps = PaddedSAM(SAMFilter(filename))
             (read,) = list(ps.queries())
             self.assertEqual(Read('query1', '-TCTAGG---'), read)
-            ps.close()
 
     def testRcNeeded(self):
         """
@@ -162,10 +343,9 @@ class TestPaddedSAM(TestCase):
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
+            ps = PaddedSAM(SAMFilter(filename))
             (read,) = list(ps.queries(rcNeeded=True))
             self.assertEqual(Read('query1', '-CCTAGA---'), read)
-            ps.close()
 
     def testRcSuffix(self):
         """
@@ -178,10 +358,9 @@ class TestPaddedSAM(TestCase):
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
+            ps = PaddedSAM(SAMFilter(filename))
             (read,) = list(ps.queries(rcSuffix='-rc'))
             self.assertEqual(Read('query1-rc', '-TCTAGG---'), read)
-            ps.close()
 
     def testQuerySoftClipLeft(self):
         """
@@ -194,10 +373,9 @@ class TestPaddedSAM(TestCase):
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
+            ps = PaddedSAM(SAMFilter(filename))
             (read,) = list(ps.queries())
             self.assertEqual(Read('query1', '-TCTAGG---'), read)
-            ps.close()
 
     def testQuerySoftClipReachesLeftEdge(self):
         """
@@ -210,10 +388,9 @@ class TestPaddedSAM(TestCase):
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
+            ps = PaddedSAM(SAMFilter(filename))
             (read,) = list(ps.queries())
             self.assertEqual(Read('query1', 'TCTAGG----'), read)
-            ps.close()
 
     def testQuerySoftClipProtrudesLeft(self):
         """
@@ -226,10 +403,9 @@ class TestPaddedSAM(TestCase):
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
+            ps = PaddedSAM(SAMFilter(filename))
             (read,) = list(ps.queries())
             self.assertEqual(Read('query1', 'AGG-------'), read)
-            ps.close()
 
     def testKF414679SoftClipLeft(self):
         """
@@ -245,10 +421,9 @@ class TestPaddedSAM(TestCase):
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
+            ps = PaddedSAM(SAMFilter(filename))
             (read,) = list(ps.queries())
             self.assertEqual(Read('query1', seq[14:]), read)
-            ps.close()
 
     def testQuerySoftClipRight(self):
         """
@@ -261,10 +436,9 @@ class TestPaddedSAM(TestCase):
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
+            ps = PaddedSAM(SAMFilter(filename))
             (read,) = list(ps.queries())
             self.assertEqual(Read('query1', '---TCTAGG-'), read)
-            ps.close()
 
     def testQuerySoftClipReachesRightEdge(self):
         """
@@ -277,10 +451,9 @@ class TestPaddedSAM(TestCase):
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
+            ps = PaddedSAM(SAMFilter(filename))
             (read,) = list(ps.queries())
             self.assertEqual(Read('query1', '----TCTAGG'), read)
-            ps.close()
 
     def testQuerySoftClipProtrudesRight(self):
         """
@@ -293,10 +466,9 @@ class TestPaddedSAM(TestCase):
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
+            ps = PaddedSAM(SAMFilter(filename))
             (read,) = list(ps.queries())
             self.assertEqual(Read('query1', '-----TCTAG'), read)
-            ps.close()
 
     def testQuerySoftClipProtrudesBothSides(self):
         """
@@ -309,10 +481,9 @@ class TestPaddedSAM(TestCase):
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
+            ps = PaddedSAM(SAMFilter(filename))
             (read,) = list(ps.queries())
             self.assertEqual(Read('query1', 'TAGGCTGACT'), read)
-            ps.close()
 
     def testQueryHardClipAndSoftClipProtrudesBothSides(self):
         """
@@ -326,10 +497,9 @@ class TestPaddedSAM(TestCase):
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
+            ps = PaddedSAM(SAMFilter(filename))
             (read,) = list(ps.queries())
             self.assertEqual(Read('query1', 'TAGGCTGACT'), read)
-            ps.close()
 
     def testReferenceInsertion(self):
         """
@@ -342,7 +512,7 @@ class TestPaddedSAM(TestCase):
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
+            ps = PaddedSAM(SAMFilter(filename))
             (read,) = list(ps.queries())
             self.assertEqual(Read('query1', '-TCGG-----'), read)
             self.assertEqual(
@@ -350,7 +520,6 @@ class TestPaddedSAM(TestCase):
                     'query1': [(3, 'TA')],
                 },
                 ps.referenceInsertions)
-            ps.close()
 
     def testPrimaryAndSecondaryReferenceInsertion(self):
         """
@@ -365,7 +534,7 @@ class TestPaddedSAM(TestCase):
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
+            ps = PaddedSAM(SAMFilter(filename))
             (read1, read2) = list(ps.queries())
             self.assertEqual(Read('query1', '-TCGG-----'), read1)
             self.assertEqual(Read('query1/1', '---TCG----'), read2)
@@ -375,7 +544,6 @@ class TestPaddedSAM(TestCase):
                     'query1/1': [(5, 'TAG')],
                 },
                 ps.referenceInsertions)
-            ps.close()
 
     def testReferenceDeletion(self):
         """
@@ -388,10 +556,9 @@ class TestPaddedSAM(TestCase):
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
+            ps = PaddedSAM(SAMFilter(filename))
             (read,) = list(ps.queries())
             self.assertEqual(Read('query1', '-TCNNTAGG-'), read)
-            ps.close()
 
     def testReferenceDeletionAlternateChar(self):
         """
@@ -404,10 +571,9 @@ class TestPaddedSAM(TestCase):
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
+            ps = PaddedSAM(SAMFilter(filename))
             (read,) = list(ps.queries(queryInsertionChar='?'))
             self.assertEqual(Read('query1', '-TC??TAGG-'), read)
-            ps.close()
 
     def testReferenceSkip(self):
         """
@@ -420,10 +586,9 @@ class TestPaddedSAM(TestCase):
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
+            ps = PaddedSAM(SAMFilter(filename))
             (read,) = list(ps.queries())
             self.assertEqual(Read('query1', '-TCNNTAGG-'), read)
-            ps.close()
 
     def testReferenceSkipAlternateChar(self):
         """
@@ -436,10 +601,9 @@ class TestPaddedSAM(TestCase):
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
+            ps = PaddedSAM(SAMFilter(filename))
             (read,) = list(ps.queries(queryInsertionChar='X'))
             self.assertEqual(Read('query1', '-TCXXTAGG-'), read)
-            ps.close()
 
     def testMixedMatchSpecificReferenceButNoMatches(self):
         """
@@ -447,15 +611,14 @@ class TestPaddedSAM(TestCase):
         has no matches must result in an empty list.
         """
         data = '\n'.join([
-            '@SQ SN:ref1 LN:10',
+            '@SQ SN:ref1 LN:15',
             '@SQ SN:ref2 LN:15',
             'query1 0 ref1 2 60 2=2X2M * 0 0 TCTAGG ZZZZZZ',
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
-            self.assertEqual([], list(ps.queries(referenceName='ref2')))
-            ps.close()
+            ps = PaddedSAM(SAMFilter(filename, referenceIds={'ref2'}))
+            self.assertEqual([], list(ps.queries()))
 
     def testMixedMatchSpecificReference(self):
         """
@@ -464,20 +627,19 @@ class TestPaddedSAM(TestCase):
         """
         data = '\n'.join([
             '@SQ SN:ref1 LN:10',
-            '@SQ SN:ref2 LN:15',
+            '@SQ SN:ref2 LN:10',
             'query1 0 ref1 2 60 2=2X2M * 0 0 TCTAGG ZZZZZZ',
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
-            (read,) = list(ps.queries(referenceName='ref1'))
+            ps = PaddedSAM(SAMFilter(filename, referenceIds={'ref1'}))
+            (read,) = list(ps.queries())
             self.assertEqual(Read('query1', '-TCTAGG---'), read)
-            ps.close()
 
     def testMinLength(self):
         """
-        A request for reads longer than a certain value should result
-        in the expected result.
+        A request for reads that are only longer than a certain value should
+        result in the expected result.
         """
         data = '\n'.join([
             '@SQ SN:ref1 LN:10',
@@ -486,27 +648,10 @@ class TestPaddedSAM(TestCase):
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
-            (read,) = list(ps.queries(minLength=6))
+            filterRead = ReadFilter(minLength=6).filter
+            ps = PaddedSAM(SAMFilter(filename, filterRead=filterRead))
+            (read,) = list(ps.queries())
             self.assertEqual(Read('query1', '-TCTAGG---'), read)
-            ps.close()
-
-    def testMinLengthWithReferenceDeletion(self):
-        """
-        The minLength specification must be applied after deletion of
-        reference bases (which results in the query being lengthened to
-        continue the match).
-        """
-        data = '\n'.join([
-            '@SQ SN:ref1 LN:10',
-            'query1 0 ref1 2 60 2M2D4M * 0 0 TCTAGG ZZZZZZ',
-        ]).replace(' ', '\t')
-
-        with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
-            (read,) = list(ps.queries(minLength=7))
-            self.assertEqual(Read('query1', '-TCNNTAGG-'), read)
-            ps.close()
 
     def testDropSecondary(self):
         """
@@ -519,10 +664,9 @@ class TestPaddedSAM(TestCase):
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
-            (read,) = list(ps.queries(dropSecondary=True))
+            ps = PaddedSAM(SAMFilter(filename, dropSecondary=True))
+            (read,) = list(ps.queries())
             self.assertEqual(Read('query1', '-TCTAGG---'), read)
-            ps.close()
 
     def testDropSupplementary(self):
         """
@@ -536,10 +680,9 @@ class TestPaddedSAM(TestCase):
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
-            (read,) = list(ps.queries(dropSupplementary=True))
+            ps = PaddedSAM(SAMFilter(filename, dropSupplementary=True))
+            (read,) = list(ps.queries())
             self.assertEqual(Read('query1', '-TCTAGG---'), read)
-            ps.close()
 
     def testDropDuplicates(self):
         """
@@ -553,10 +696,9 @@ class TestPaddedSAM(TestCase):
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
-            (read,) = list(ps.queries(dropDuplicates=True))
+            ps = PaddedSAM(SAMFilter(filename, dropDuplicates=True))
+            (read,) = list(ps.queries())
             self.assertEqual(Read('query1', '-TCTAGG---'), read)
-            ps.close()
 
     def testAllowDuplicateIds(self):
         """
@@ -570,11 +712,10 @@ class TestPaddedSAM(TestCase):
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
+            ps = PaddedSAM(SAMFilter(filename))
             (read1, read2) = list(ps.queries(allowDuplicateIds=True))
             self.assertEqual(Read('query1', '-TCTAGG---'), read1)
             self.assertEqual(Read('query1', '--TC------'), read2)
-            ps.close()
 
     def testDuplicateIdDisambiguation(self):
         """
@@ -588,12 +729,11 @@ class TestPaddedSAM(TestCase):
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
+            ps = PaddedSAM(SAMFilter(filename))
             (read1, read2, read3) = list(ps.queries())
             self.assertEqual(Read('query1', '-TCTAGG---'), read1)
             self.assertEqual(Read('query1/1', '--TC------'), read2)
             self.assertEqual(Read('query1/2', '--TCGA----'), read3)
-            ps.close()
 
     def testKeepQualityControlFailures(self):
         """
@@ -607,11 +747,10 @@ class TestPaddedSAM(TestCase):
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
-            (read1, read2) = list(ps.queries(keepQCFailures=True))
+            ps = PaddedSAM(SAMFilter(filename, keepQCFailures=True))
+            (read1, read2) = list(ps.queries())
             self.assertEqual(Read('query1', '-TCTAGG---'), read1)
             self.assertEqual(Read('query2', '---TC-----'), read2)
-            ps.close()
 
     def testSecondaryWithNoPreviousSequence(self):
         """
@@ -624,12 +763,11 @@ class TestPaddedSAM(TestCase):
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
+            ps = PaddedSAM(SAMFilter(filename))
             error = ('^Query line 1 has an empty SEQ field, but no previous '
                      'alignment is present\\.$')
             queries = ps.queries()
-            assertRaisesRegex(self, ValueError, error, list, queries)
-            ps.close()
+            assertRaisesRegex(self, InvalidSAM, error, list, queries)
 
     def testSecondaryWithNoSequence(self):
         """
@@ -644,12 +782,11 @@ class TestPaddedSAM(TestCase):
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
+            ps = PaddedSAM(SAMFilter(filename))
             (read1, read2, read3) = list(ps.queries())
             self.assertEqual(Read('query1', '-TCT------'), read1)
             self.assertEqual(Read('query2', '-TCTA-----'), read2)
             self.assertEqual(Read('query2/1', '-----TCTA-'), read3)
-            ps.close()
 
     def testSupplementaryWithNoPreviousSequence(self):
         """
@@ -662,12 +799,11 @@ class TestPaddedSAM(TestCase):
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
+            ps = PaddedSAM(SAMFilter(filename))
             error = ('^Query line 1 has an empty SEQ field, but no previous '
                      'alignment is present\\.$')
             queries = ps.queries()
-            assertRaisesRegex(self, ValueError, error, list, queries)
-            ps.close()
+            assertRaisesRegex(self, InvalidSAM, error, list, queries)
 
     def testSupplementaryWithNoSequence(self):
         """
@@ -682,12 +818,11 @@ class TestPaddedSAM(TestCase):
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
+            ps = PaddedSAM(SAMFilter(filename))
             (read1, read2, read3) = list(ps.queries())
             self.assertEqual(Read('query1', '-TCT------'), read1)
             self.assertEqual(Read('query2', '-TCTA-----'), read2)
             self.assertEqual(Read('query2/1', '-----TCTA-'), read3)
-            ps.close()
 
     def testNotSecondaryAndNotSupplementaryWithNoSequence(self):
         """
@@ -700,12 +835,11 @@ class TestPaddedSAM(TestCase):
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
+            ps = PaddedSAM(SAMFilter(filename))
             error = ('^Query line 1 has an empty SEQ field, but the alignment '
                      'is not marked as secondary or supplementary\\.$')
             queries = ps.queries()
-            assertRaisesRegex(self, ValueError, error, list, queries)
-            ps.close()
+            assertRaisesRegex(self, InvalidSAM, error, list, queries)
 
     def testAlsoYieldAlignments(self):
         """
@@ -719,7 +853,7 @@ class TestPaddedSAM(TestCase):
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
+            ps = PaddedSAM(SAMFilter(filename))
             (read1, read2) = list(ps.queries(addAlignment=True))
 
             self.assertEqual(Read('query1', '-TCTAGG---'), read1)
@@ -732,38 +866,34 @@ class TestPaddedSAM(TestCase):
             self.assertEqual('XY', ''.join(
                 map(lambda x: chr(x + 33), read2.alignment.query_qualities)))
 
-            ps.close()
 
-    def testAlignmentCount(self):
+class TestSamReferencesToStr(TestCase):
+    """
+    Test the samReferencesToStr function.
+    """
+    def testSimple(self):
         """
-        When all queries have been yielded, the alignment count must be
-        as expected.
+        The referencesToStr method must return the expected string.
         """
         data = '\n'.join([
-            '@SQ SN:ref1 LN:10',
-            'query1 0 ref1 2 60 2=2X2M * 0 0 TCTAGG 123456',
-            'query2 0 ref1 2 60 2= * 0 0 TC XY',
+            '@SQ SN:id1 LN:90',
+            '@SQ SN:id2 LN:91',
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
-            list(ps.queries(addAlignment=True))
-            self.assertEqual(2, ps.alignmentCount)
-            ps.close()
+            self.assertEqual('id1 (length 90)\nid2 (length 91)',
+                             samReferencesToStr(filename))
 
-    def testStoreQueryIds(self):
+    def testIndent(self):
         """
-        If we request that query ids are saved, they must be.
+        The referencesToStr method must return the expected string when
+        passed an indent.
         """
         data = '\n'.join([
-            '@SQ SN:ref1 LN:10',
-            'query1 0 ref1 2 60 2=2X2M * 0 0 TCTAGG 123456',
-            'query2 0 ref1 2 60 2= * 0 0 TC XY',
-            'query2 0 ref1 2 60 2= * 0 0 TC XY',
+            '@SQ SN:id1 LN:90',
+            '@SQ SN:id2 LN:91',
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
-            ps = PaddedSAM(filename)
-            list(ps.queries(storeQueryIds=True))
-            self.assertEqual({'query1', 'query2'}, ps.queryIds)
-            ps.close()
+            self.assertEqual('  id1 (length 90)\n  id2 (length 91)',
+                             samReferencesToStr(filename, indent=2))
diff --git a/test/test_utils.py b/test/test_utils.py
index 298d15f..f1db1b1 100644
--- a/test/test_utils.py
+++ b/test/test_utils.py
@@ -4,6 +4,7 @@ import gzip
 from six.moves import builtins
 from unittest import TestCase
 from six import assertRaisesRegex
+from collections import Counter
 
 try:
     from unittest.mock import patch
@@ -13,7 +14,8 @@ except ImportError:
 from .mocking import mockOpen, File
 
 from dark.utils import (
-    numericallySortFilenames, median, asHandle, parseRangeString, StringIO)
+    numericallySortFilenames, median, asHandle, parseRangeString, StringIO,
+    baseCountsToStr, nucleotidesToStr)
 
 
 class TestNumericallySortFilenames(TestCase):
@@ -313,3 +315,43 @@ class TestStringIO(TestCase):
         with StringIO() as s:
             s.write('hey')
             self.assertEqual('hey', s.getvalue())
+
+
+class TestBaseCountsToStr(TestCase):
+    """
+    Test the baseCountsToStr function.
+    """
+    def testSimple(self):
+        """
+        A simple example must work as expected.
+        """
+        counts = Counter()
+        counts['A'] += 1
+        counts['G'] += 2
+        self.assertEqual('A:1 G:2',
+                         baseCountsToStr(counts))
+
+
+class TestNucleotidesToStr(TestCase):
+    """
+    Test the nucleotidesToStr function.
+    """
+    def testSimple(self):
+        """
+        A simple example must work as expected.
+        """
+        counts1 = Counter()
+        counts1['A'] += 1
+        counts1['G'] += 2
+        counts2 = Counter()
+        counts2['C'] += 1
+        counts2['T'] += 3
+        self.assertEqual(
+            '0: A:1 G:2\n7: C:1 T:3',
+            nucleotidesToStr(
+                {
+                    0: counts1,
+                    7: counts2,
+                }
+            )
+        )