class Read(object):
    """
    Hold information about a single read.

    @param id: A C{str} describing the read.
    @param sequence: A C{str} of sequence information (might be
        nucleotides or proteins).
    @param quality: An optional C{str} of phred quality scores. If not C{None},
        it must be the same length as C{sequence}.
    @raise ValueError: if the length of the quality string (if any) does not
        match the length of the sequence.
    """
    ALPHABET = None

    def __init__(self, id, sequence, quality=None):
        if quality is not None and len(quality) != len(sequence):
            raise ValueError(
                'Invalid read: sequence length (%d) != quality length (%d)' %
                (len(sequence), len(quality)))

        self.id = id
        self.sequence = sequence
        self.quality = quality

    def __eq__(self, other):
        return (self.id == other.id and
                self.sequence == other.sequence and
                self.quality == other.quality)

    def __ne__(self, other):
        return not self == other

    def __lt__(self, other):
        return ((self.id, self.sequence, self.quality) <
                (other.id, other.sequence, other.quality))

    def __len__(self):
        return len(self.sequence)

    def __hash__(self):
        """
        Calculate a hash key for a read.

        @return: The C{int} hash key for the read.
        """
        if self.quality is None:
            return hash(md5(self.id.encode('UTF-8') + b'\0' +
                            self.sequence.encode('UTF-8')).digest())
        else:
            return hash(md5(self.id.encode('UTF-8') + b'\0' +
                            self.sequence.encode('UTF-8') + b'\0' +
                            self.quality.encode('UTF-8')).digest())

    def __getitem__(self, item):
        sequence = self.sequence[item]
        quality = None if self.quality is None else self.quality[item]
        return self.__class__(self.id, sequence, quality)

    def toString(self, format_):
        """
        Convert the read to a string format.

        @param format_: Either 'fasta', 'fastq' or 'fasta-ss'.
        @raise ValueError: if C{format_} is 'fastq' and the read has no quality
            information, if C{format_} is 'fasta-ss' and the read has no
            structure information, or if an unknown format is requested.
        @return: A C{str} representing the read in the requested format.
        """
        if format_ == 'fasta':
            return '>%s\n%s\n' % (self.id, self.sequence)
        elif format_ == 'fastq':
            if self.quality is None:
                raise ValueError('Read %r has no quality information' %
                                 self.id)
            else:
                return '@%s\n%s\n+%s\n%s\n' % (
                    self.id, self.sequence, self.id, self.quality)
        else:
            raise ValueError("Format must be either 'fasta', 'fastq' or "
                             "'fasta-ss'.")

    def toDict(self):
        """
        Get information about this read in a dictionary.

        @return: A C{dict} with keys/values for the attributes of self.
        """
        return {
            'id': self.id,
            'sequence': self.sequence,
            'quality': self.quality,
        }

    def reverse(self):
        """
        Reverse a read (note that this is NOT a reverse complement).

        @return: The reversed sequence as an instance of the current class.
        """
        return self.__class__(
            self.id,
            self.sequence[::-1],
            None if self.quality is None else self.quality[::-1])

    @classmethod
    def fromDict(cls, d):
        """
        Create a new instance from attribute values provided in a dictionary.

        @param d: A C{dict} with keys/values for the attributes of a new
            instance of this class. Keys 'id' and 'sequence' with C{str} values
            must be provided. A 'quality' C{str} key is optional.
        @return: A new instance of this class, with values taken from C{d}.
        """
        return cls(d['id'], d['sequence'], d.get('quality'))

    def lowComplexityFraction(self):
        """
        What fraction of a read's bases are in low-complexity regions?
        By convention, a region of low complexity is indicated by lowercase
        base letters.

        @return: The C{float} representing the fraction of bases in the
            read that are in regions of low complexity.
        """
        length = len(self)
        if length:
            lowerCount = len(list(filter(str.islower, self.sequence)))
            return float(lowerCount) / length
        else:
            return 0.0

    def walkHSP(self, hsp, includeWhiskers=True):
        """
        Provide information about exactly how a read matches a subject, as
        specified by C{hsp}.

        @param includeWhiskers: If C{True} yield information from the
            (possibly empty) non-matching ends of the read.
        @return: A generator that yields (offset, residue, inMatch) tuples.
            The offset is the offset into the matched subject. The residue is
            the base in the read (which might be '-' to indicate a gap in the
            read was aligned with the subject at this offset). inMatch will be
            C{True} for residues that are part of the HSP match, and C{False}
            for the (possibly non-existent) parts of the read that fall outside
            the HSP (aka, the "whiskers" in an alignment graph).
        """

        # It will be easier to understand the following implementation if
        # you refer to the ASCII art illustration of an HSP in the
        # dark.hsp._Base class in hsp.py

        # Left whisker.
        if includeWhiskers:
            readOffset = 0
            subjectOffset = hsp.readStartInSubject
            while subjectOffset < hsp.subjectStart:
                yield (subjectOffset, self.sequence[readOffset], False)
                readOffset += 1
                subjectOffset += 1

        # Match.
        for subjectOffset, residue in enumerate(hsp.readMatchedSequence,
                                                start=hsp.subjectStart):
            yield (subjectOffset, residue, True)

        # Right whisker.
        if includeWhiskers:
            readOffset = hsp.readEnd
            subjectOffset = hsp.subjectEnd
            while subjectOffset < hsp.readEndInSubject:
                yield (subjectOffset, self.sequence[readOffset], False)
                readOffset += 1
                subjectOffset += 1

    def checkAlphabet(self, count=10):
        """
        A function which checks whether the sequence in a L{dark.Read} object
        corresponds to its readClass. For AA reads, more testing is done in
        dark.Read.AARead.checkAlphabet.

        @param count: An C{int}, indicating how many bases or amino acids at
            the start of the sequence should be considered. If C{None}, all
            bases are checked.
        @return: C{True} if the alphabet characters in the first C{count}
            positions of sequence is a subset of the allowed alphabet for this
            read class, or if the read class has a C{None} alphabet.
        @raise ValueError: If the sequence alphabet is not a subset of the read
            class alphabet.
        """
        if count is None:
            readLetters = set(self.sequence.upper())
        else:
            readLetters = set(self.sequence.upper()[:count])
        # Check if readLetters is a subset of self.ALPHABET.
        if self.ALPHABET is None or readLetters.issubset(self.ALPHABET):
            return readLetters
        raise ValueError("Read alphabet (%r) is not a subset of expected "
                         "alphabet (%r) for read class %s." % (
                             ''.join(sorted(readLetters)),
                             ''.join(sorted(self.ALPHABET)),
                             str(self.__class__.__name__)))

    def newFromSites(self, sites, exclude=False):
        """
        Create a new read from self, with only certain sites.

        @param sites: A set of C{int} 0-based sites (i.e., indices) in
            sequences that should be kept. If C{None} (the default), all sites
            are kept.
        @param exclude: If C{True} the C{sites} will be excluded, not
            included.
        """
        if exclude:
            sites = set(range(len(self))) - sites

        newSequence = []
        if self.quality:
            newQuality = []
            for index, (base, quality) in enumerate(zip(self.sequence,
                                                        self.quality)):
                if index in sites:
                    newSequence.append(base)
                    newQuality.append(quality)
            read = self.__class__(self.id, ''.join(newSequence),
                                  ''.join(newQuality))
        else:
            for index, base in enumerate(self.sequence):
                if index in sites:
                    newSequence.append(base)
            read = self.__class__(self.id, ''.join(newSequence))

        return read class ReadFilter(object):
    """
    Create a function that can be used to filter a set of reads to produce a
    desired subset.

    Note: there are many additional filtering options that could be added,
    e.g., on complexity fraction, on GC %, on quality, etc.

    @param minLength: The minimum acceptable length.
    @param maxLength: The maximum acceptable length.
    @param maxNFraction: The maximum fraction of Ns that can be present in the
        sequence.
    @param removeGaps: If C{True} remove all gaps ('-' characters) from the
        read sequences.
    @param whitelist: If not C{None}, a set of exact read ids that are
        always acceptable (though other characteristics, such as length,
        of a whitelisted id may rule it out).
    @param blacklist: If not C{None}, a set of exact read ids that are
        never acceptable.
    @param whitelistFile: If not C{None}, a C{str} filename containing lines
        that give exact ids that are always acceptable.
    @param blacklistFile: If not C{None}, a C{str} filename containing lines
        that give exact ids that are never acceptable.
    @param titleRegex: A regex that read ids must match.
    @param negativeTitleRegex: A regex that read ids must not match.
    @param truncateTitlesAfter: A string that read ids will be truncated
        beyond. If the truncated version of an id has already been seen,
        that sequence will be skipped.
    @param keepSequences: Either C{None} or a set of C{int}s corresponding
        to reads that should be kept. Indexing starts at zero. The numbers
        refer to the sequential number of the sequence in the list of reads.
    @param removeSequences: Either C{None} or a set of C{int}s corresponding
        to reads that should be removed. Indexing starts at zero. The numbers
        refer to the sequential number of the sequence in the list of reads.
    @param head: If not C{None}, the C{int} number of sequences at the
        start of the reads to return. Later sequences are skipped.
    @param removeDuplicates: If C{True} remove duplicated reads based only on
        sequence identity.
    @param removeDuplicatesById: If C{True} remove duplicated reads based
        only on read id.
    @param removeDuplicatesUseMD5: If C{True}, use MD5 sums instead of the
        full sequence or read when either C{removeDuplicates} or
        C{removeDuplicatesById} are C{True}.
    @param removeDescriptions: If C{True} remove the description (the part
        following the first whitespace) from read ids. The description is
        removed after applying the function specified by --idLambda (if any).
    @param modifier: If not C{None}, a function that is passed a read
        and which either returns a read or C{None}. If it returns a read,
        that read is passed through the filter. If it returns C{None},
        the read is omitted. Such a function can be used to do customized
        filtering, to change sequence ids, etc.
    @param randomSubset: If not C{None}, an C{int} giving the number of
        sequences that should be returned. These will be selected at
        random, using an algorithm that does a single pass over the data
        and which needs to know (in advance) the total number of reads.
        Because a C{Reads} instance does not always know its length, it is
        possible to specify the number of reads by passing a C{trueLength}
        argument. Note that the random selection is done before any other
        filtering. Due to this, if you want to extract a random subset of
        the reads filtered in another way, it will be best to call filter
        twice rather than doing both types of filtering in one step. E.g.,
        you very likely should do this:
            reads.filter(maxLength=100).filter(randomSubset=20)
        rather than this:
            reads.filter(maxLength=100, randomSubset=20)
        The second version will extract a random subset of 20 reads and
        only return those that have length <= 100, so your result may have
        less than 20 reads. The former version extracts reads of the
        desired length and then takes 20 reads at random from that set, so
        you'll always get 20 reads in your result, assuming there are at
        least that many reads satisfying the length filter.
    @param trueLength: The C{int} number of reads in this C{Reads} instance.
        Under normal circumstances it will not be necessary to pass this
        argument. However in some cases a subclass (e.g., with
        C{dark.fasta.FastaReads}) does not know its length until its data
        has been read from disk. In such cases, it is not possible to
        choose a random subset without keeping the subset in memory (which
        is undesirable). See
        https://en.wikipedia.org/wiki/Reservoir_sampling for one approach
        when the set size is unknown. However, it is possible to filter a
        random subset in a single pass over the data without keeping the
        set in memory if the set size is known. C{trueLength} makes it
        possible to pass the actual number of reads (this will obviously
        need to be obtained via some other mechanism).
    @param sampleFraction: If not C{None}, a [0.0, 1.0] C{float} indicating
        a fraction of the reads that should be allowed to pass through the
        filter. The sample size will only be approximately the product of
        the C{sampleFraction} and the number of reads. The sample is taken
        at random. If you try to combine this filter with C{randomSubset}
        a C{ValueError} will be raised. If you need both filters, run them
        one after another.
    @param sequenceNumbersFile: If not C{None}, gives the C{str} name of a
        file containing (1-based) sequence numbers, in ascending order,
        one per line. Only those sequences matching the given numbers will
        be kept.
    @param idLambda: If not C{None}, a C{str} Python lambda function
        specification to use to modify read ids. The function is applied
        before removing the description (if --removeDescriptions is also
        specified).
    @param readLambda: If not C{None}, a C{str} Python lambda function
        specification to use to modify reads. The function will be passed,
        and must return, a single Read (or one of its subclasses). This
        function is called after the --idLambda function, if any.
    @param keepSites: A set of C{int} 0-based sites (i.e., indices) in
        sequences that should be kept. If C{None} (the default), all sites are
        kept.
    @param removeSites: A set of C{int} 0-based sites (i.e., indices) in
        sequences that should be removed. If C{None} (the default), no sites
        are removed.
    @param reverse: If C{True}, reverse the sequences. Reversing happens
        at a very late stage (i.e., after sites are altered via keepSites
        and removeSites).
    @param reverseComplement: If C{True}, replace seqeunces with their reverse
        complements. Reversing happens at a very late stage (i.e., after sites
        are altered via keepSites and removeSites).
    @raises ValueError: If C{randomSubset} and C{sampleFraction} are both
        specified, or if C{randomSubset} is specified but C{trueLength} is not,
        or if the sequence numbers in C{sequenceNumbersFile} are
        non-positive or not ascending, or if both C{keepSites} and
        C{removeSites} are given, or if both C{keepSequences} and
        C{removeSequences} are given.
    @raise AttributeError: If C{reverseComplement} is C{True} but the read
        type does not allow for reverse complementing.
    """

    # TODO, when/if needed: make it possible to pass a seed for the RNG
    # when randomSubset or sampleFraction are used. Also possible is to
    # save and restore the state of the RNG and/or to optionally add
    # 'seed=XXX' to the end of the id of the first read, etc.

    def __init__(self, minLength=None, maxLength=None, maxNFraction=None,
                 removeGaps=False, whitelist=None, blacklist=None,
                 whitelistFile=None, blacklistFile=None,
                 titleRegex=None, negativeTitleRegex=None,
                 truncateTitlesAfter=None, keepSequences=None,
                 removeSequences=None, head=None,
                 removeDuplicates=False, removeDuplicatesById=False,
                 removeDuplicatesUseMD5=False, removeDescriptions=False,
                 modifier=None, randomSubset=None, trueLength=None,
                 sampleFraction=None, sequenceNumbersFile=None, idLambda=None,
                 readLambda=None, keepSites=None, removeSites=None,
                 reverse=False, reverseComplement=False):

        if randomSubset is not None:
            if sampleFraction is not None:
                raise ValueError('randomSubset and sampleFraction cannot be '
                                 'used simultaneously in a filter. Make two '
                                 'read filters instead.')

            if trueLength is None:
                raise ValueError('trueLength must be supplied if randomSubset '
                                 'is specified.')

        self.minLength = minLength
        self.maxLength = maxLength
        self.maxNFraction = maxNFraction
        self.removeGaps = removeGaps
        self.head = head
        self.removeDescriptions = removeDescriptions
        self.modifier = modifier
        self.randomSubset = randomSubset
        self.trueLength = trueLength

        if removeDuplicatesUseMD5 and not (
                removeDuplicates or removeDuplicatesById):
            raise ValueError(
                'If you specify removeDuplicatesUseMD5, you need to also use '
                'one of removeDuplicates or removeDuplicatesById.')
        self.removeDuplicates = removeDuplicates
        self.removeDuplicatesById = removeDuplicatesById
        self.removeDuplicatesUseMD5 = removeDuplicatesUseMD5

        if keepSequences and removeSequences:
            raise ValueError(
                'Cannot simultaneously filter using keepSequences and '
                'removeSequences. Call filter twice in succession instead.')
        self.keepSequences = keepSequences
        self.removeSequences = removeSequences

        if keepSites and removeSites:
            raise ValueError(
                'Cannot simultaneously filter using keepSites and '
                'removeSites. Call filter twice in succession instead.')
        self.keepSites = keepSites
        self.removeSites = removeSites

        if reverseComplement:
            # Make sure reverse is not also set.
            reverse = False
        self.reverse = reverse
        self.reverseComplement = reverseComplement

        self.alwaysFalse = False
        self.yieldCount = 0
        self.readIndex = -1

        def _wantedSequences(filename):
            """
            Read and yield integer sequence numbers from a file.

            @raise ValueError: If the sequence numbers are not all positive or
                are not ascending.
            @return: A generator that yields C{int} sequence numbers.
            """
            with open(filename) as fp:
                lastNumber = None
                for line in fp:
                    n = int(line)
                    if lastNumber is None:
                        if n < 1:
                            raise ValueError(
                                'First line of sequence number file %r must '
                                'be at least 1.' % filename)
                        lastNumber = n
                        yield n
                    else:
                        if n > lastNumber:
                            lastNumber = n
                            yield n
                        else:
                            raise ValueError(
                                'Line number file %r contains non-ascending '
                                'numbers %d and %d.' %
                                (filename, lastNumber, n))

        self.wantedSequenceNumberGeneratorExhausted = False
        self.nextWantedSequenceNumber = None

        if sequenceNumbersFile is not None:
            self.wantedSequenceNumberGenerator = _wantedSequences(
                sequenceNumbersFile)
            try:
                self.nextWantedSequenceNumber = next(
                    self.wantedSequenceNumberGenerator)
            except StopIteration:
                # There was a sequence number file, but it was empty. So no
                # reads will ever be accepted.
                self.alwaysFalse = True

        if (whitelist or blacklist or whitelistFile or blacklistFile or
                titleRegex or negativeTitleRegex or truncateTitlesAfter):
            self.titleFilter = TitleFilter(
                whitelist=whitelist, blacklist=blacklist,
                whitelistFile=whitelistFile, blacklistFile=blacklistFile,
                positiveRegex=titleRegex, negativeRegex=negativeTitleRegex,
                truncateAfter=truncateTitlesAfter)
        else:
            self.titleFilter = None

        if removeDuplicates:
            self.sequencesSeen = set()

        if removeDuplicatesById:
            self.idsSeen = set()

        if sampleFraction is not None:
            if sampleFraction == 0.0:
                # The filter method should always return False.
                self.alwaysFalse = True
            elif sampleFraction == 1.0:
                # Passing 1.0 can be treated the same as passing no value.
                # This makes the filter code below simpler.
                sampleFraction = None
        self.sampleFraction = sampleFraction

        self.idLambda = eval(idLambda) if idLambda else None
        self.readLambda = eval(readLambda) if readLambda else None

    def filter(self, read):
        """
        Check if a read passes the filter.

        @param read: A C{Read} instance.
        @return: C{read} if C{read} passes the filter, C{False} if not.
        """
        self.readIndex += 1

        if self.alwaysFalse:
            return False

        if self.wantedSequenceNumberGeneratorExhausted:
            return False

        if self.nextWantedSequenceNumber is not None:
            if self.readIndex + 1 == self.nextWantedSequenceNumber:
                # We want this sequence.
                try:
                    self.nextWantedSequenceNumber = next(
                        self.wantedSequenceNumberGenerator)
                except StopIteration:
                    # The sequence number iterator ran out of sequence
                    # numbers.  We must let the rest of the filtering
                    # continue for the current sequence in case we
                    # throw it out for other reasons (as we might have
                    # done for any of the earlier wanted sequence
                    # numbers).
                    self.wantedSequenceNumberGeneratorExhausted = True
            else:
                # This sequence isn't one of the ones that's wanted.
                return False

        if (self.sampleFraction is not None and
                uniform(0.0, 1.0) > self.sampleFraction):
            # Note that we don't have to worry about the 0.0 or 1.0
            # cases in the above 'if', as they have been dealt with
            # in self.__init__.
            return False

        if self.randomSubset is not None:
            if self.yieldCount == self.randomSubset:
                # The random subset has already been fully returned.
                # There's no point in going any further through the input.
                self.alwaysFalse = True
                return False
            elif uniform(0.0, 1.0) > ((self.randomSubset - self.yieldCount) /
                                      (self.trueLength - self.readIndex)):
                return False

        if self.head is not None and self.readIndex == self.head:
            # We're completely done.
            self.alwaysFalse = True
            return False

        readLen = len(read)
        if ((self.minLength is not None and readLen < self.minLength) or
                (self.maxLength is not None and readLen > self.maxLength)):
            return False

        if self.maxNFraction is not None:
            nFraction = read.sequence.count('N') / readLen
            if self.maxNFraction < nFraction:
                return False

        if self.removeGaps:
            if read.quality is None:
                read = read.__class__(read.id, read.sequence.replace('-', ''))
            else:
                newSequence = []
                newQuality = []
                for base, quality in zip(read.sequence, read.quality):
                    if base != '-':
                        newSequence.append(base)
                        newQuality.append(quality)
                read = read.__class__(
                    read.id, ''.join(newSequence), ''.join(newQuality))

        if (self.titleFilter and
                self.titleFilter.accept(read.id) == TitleFilter.REJECT):
            return False

        if (self.keepSequences is not None and
                self.readIndex not in self.keepSequences):
            return False

        if (self.removeSequences is not None and
                self.readIndex in self.removeSequences):
            return False

        if self.removeDuplicates:
            if self.removeDuplicatesUseMD5:
                sequence = md5(read.sequence.encode('UTF-8')).digest()
            else:
                sequence = read.sequence
            if sequence in self.sequencesSeen:
                return False
            self.sequencesSeen.add(sequence)

        if self.removeDuplicatesById:
            if self.removeDuplicatesUseMD5:
                id_ = md5(read.id.encode('UTF-8')).digest()
            else:
                id_ = read.id
            if id_ in self.idsSeen:
                return False
            self.idsSeen.add(id_)

        if self.modifier:
            modified = self.modifier(read)
            if modified is None:
                return False
            else:
                read = modified

        # We have to use 'is not None' in the following tests so the empty set
        # is processed properly.
        if self.keepSites is not None:
            read = read.newFromSites(self.keepSites)
        elif self.removeSites is not None:
            read = read.newFromSites(self.removeSites, exclude=True)

        if self.idLambda:
            newId = self.idLambda(read.id)
            if newId is None:
                return False
            else:
                read.id = newId

        if self.readLambda:
            newRead = self.readLambda(read)
            if newRead is None:
                return False
            else:
                read = newRead

        if self.removeDescriptions:
            read.id = read.id.split()[0]

        if self.reverse:
            read = read.reverse()
        elif self.reverseComplement:
            read = read.reverseComplement()

        self.yieldCount += 1
        return read class PaddedSAM(object):
    """
    Obtain aligned (padded) queries from a 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, 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)

    def queries(self, rcSuffix='', rcNeeded=False, padChar='-',
                queryInsertionChar='N', unknownQualityChar='!',
                allowDuplicateIds=False, addAlignment=False):
        """
        Produce padded (with gaps) queries according to the CIGAR string and
        reference sequence length for each matching query sequence.

        @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 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
            that created the SAM/BAM input flags reversed matches but does not
            also store the reverse complemented query.
        @param padChar: A C{str} of length one to use to pad queries with to
            make them the same length as the reference sequence.
        @param queryInsertionChar:  A C{str} of length one to use to insert
            into queries when the CIGAR string indicates that the alignment
            of a query would cause a deletion in the reference. This character
            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 unknownQualityChar: The character to put into the quality
            string when unknown bases are inserted in the query or the query
            is padded on the left/right with gaps.
        @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 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}.
        """
        referenceLength = self.referenceLength

        # Hold the count for each id so we can add /1, /2 etc to duplicate
        # ids (unless --allowDuplicateIds was given).
        idCount = Counter()

        MATCH_OPERATIONS = {CMATCH, CEQUAL, CDIFF}

        for lineNumber, alignment in enumerate(
                self.samFilter.alignments(), start=1):

            query = alignment.query_sequence
            quality = ''.join(chr(q + 33) for q in alignment.query_qualities)

            if alignment.is_reverse:
                if rcNeeded:
                    query = DNARead('id', query).reverseComplement().sequence
                    quality = quality[::-1]
                if rcSuffix:
                    alignment.query_name += rcSuffix

            # Adjust the query id if it's a duplicate and we're not allowing
            # duplicates.
            if allowDuplicateIds:
                queryId = alignment.query_name
            else:
                count = idCount[alignment.query_name]
                idCount[alignment.query_name] += 1
                queryId = alignment.query_name + (
                    '' if count == 0 else '/%d' % count)

            referenceStart = alignment.reference_start
            atStart = True
            queryIndex = 0
            referenceIndex = referenceStart
            alignedSequence = ''
            alignedQuality = ''

            for operation, length in alignment.cigartuples:

                # The operations are tested in the order they appear in
                # https://samtools.github.io/hts-specs/SAMv1.pdf It would be
                # more efficient to test them in order of frequency of
                # occurrence.
                if operation in MATCH_OPERATIONS:
                    atStart = False
                    alignedSequence += query[queryIndex:queryIndex + length]
                    alignedQuality += quality[queryIndex:queryIndex + length]
                elif operation == CINS:
                    # Insertion to the reference. This consumes query bases but
                    # we don't output them because the reference cannot be
                    # changed.  I.e., these bases in the query would need to be
                    # inserted into the reference.  Remove these bases from the
                    # query but record what would have been inserted into the
                    # reference.
                    atStart = False
                    self.referenceInsertions[queryId].append(
                        (referenceIndex,
                         query[queryIndex:queryIndex + length]))
                elif operation == CDEL:
                    # Delete from the reference. Some bases from the reference
                    # would need to be deleted to continue the match. So we put
                    # an insertion into the query to compensate.
                    atStart = False
                    alignedSequence += queryInsertionChar * length
                    alignedQuality += unknownQualityChar * length
                elif operation == CREF_SKIP:
                    # Skipped reference. Opens a gap in the query. For
                    # mRNA-to-genome alignment, an N operation represents an
                    # intron.  For other types of alignments, the
                    # interpretation of N is not defined. So this is unlikely
                    # to occur.
                    atStart = False
                    alignedSequence += queryInsertionChar * length
                    alignedQuality += unknownQualityChar * length
                elif operation == CSOFT_CLIP:
                    # Bases in the query that are not part of the match. We
                    # remove these from the query if they protrude before the
                    # start or after the end of the reference. According to the
                    # SAM docs, 'S' operations may only have 'H' operations
                    # between them and the ends of the CIGAR string.
                    if atStart:
                        # Don't set atStart=False, in case there's another 'S'
                        # operation.
                        unwantedLeft = length - referenceStart
                        if unwantedLeft > 0:
                            # The query protrudes left. Copy its right part.
                            alignedSequence += query[
                                queryIndex + unwantedLeft:queryIndex + length]
                            alignedQuality += quality[
                                queryIndex + unwantedLeft:queryIndex + length]
                            referenceStart = 0
                        else:
                            referenceStart -= length
                            alignedSequence += query[
                                queryIndex:queryIndex + length]
                            alignedQuality += quality[
                                queryIndex:queryIndex + length]
                    else:
                        unwantedRight = (
                            (referenceStart + len(alignedSequence) + length) -
                            referenceLength)

                        if unwantedRight > 0:
                            # The query protrudes right. Copy its left part.
                            alignedSequence += query[
                                queryIndex:queryIndex + length - unwantedRight]
                            alignedQuality += quality[
                                queryIndex:queryIndex + length - unwantedRight]
                        else:
                            alignedSequence += query[
                                queryIndex:queryIndex + length]
                            alignedQuality += quality[
                                queryIndex:queryIndex + length]
                elif operation == CHARD_CLIP:
                    # Some bases have been completely removed from the query.
                    # This (H) can only be present as the first and/or last
                    # operation. There is nothing to do as the bases are simply
                    # not present in the query string in the SAM/BAM file.
                    pass
                elif operation == CPAD:
                    # This is "silent deletion from the padded reference",
                    # which consumes neither query nor reference.
                    atStart = False
                else:
                    raise ValueError('Unknown CIGAR operation:', operation)

                if operation in CONSUMES_QUERY:
                    queryIndex += length

                if operation in CONSUMES_REFERENCE:
                    referenceIndex += length

            if queryIndex != len(query):
                # Oops, we did not consume the entire query.
                raise ValueError(
                    'Query %r not fully consumed when parsing CIGAR string. '
                    'Query %r (len %d), final query index %d, CIGAR: %r' %
                    (alignment.query_name, query, len(query), queryIndex,
                     alignment.cigartuples))

            # We cannot test we consumed the entire reference.  The CIGAR
            # 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.
            padRightLength = (referenceLength -
                              (referenceStart + len(alignedSequence)))
            paddedSequence = (padChar * referenceStart +
                              alignedSequence +
                              padChar * padRightLength)
            paddedQuality = (unknownQualityChar * referenceStart +
                             alignedQuality +
                             unknownQualityChar * padRightLength)

            read = Read(queryId, paddedSequence, paddedQuality)

            if addAlignment:
                read.alignment = alignment

            yield read 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.storeQueryIds = storeQueryIds
        self.dropUnmapped = dropUnmapped
        self.dropSecondary = dropSecondary
        self.dropSupplementary = dropSupplementary
        self.dropDuplicates = dropDuplicates
        self.keepQCFailures = keepQCFailures
        self.minScore = minScore
        self.maxScore = maxScore
        self.scoreTag = scoreTag

        # Detect when there are no filtering criteria, in which case
        # self.filterAlignment can return immediately.
        self.noFiltering = all((
            filterRead is None, referenceIds is None,
            not any((storeQueryIds, dropUnmapped, dropSecondary,
                     dropSupplementary, dropDuplicates, keepQCFailures)),
            minScore is None, maxScore is None))

    @staticmethod
    def addFilteringOptions(parser, samfileIsPositional=False,
                            samfileAction='store', samfileRequired=True,
                            samfileNargs=None, referenceIdRequired=False):
        """
        Add options to an argument parser for filtering SAM/BAM.

        @param parser: An C{argparse.ArgumentParser} instance.
        @param samfileIsPositional: If C{True} the SAM/BAM file must
            be given as the final argument on the command line (without
            being preceded by --samfile).
        @param samfileAction: A C{str} action to take when 'samfile' arguments
            are found on the command line. Pass 'append' to allow multiple SAM
            files.
        @param samfileRequired: If C{True}, --samfile must be given on the
            command line. This is only relevant when samfileIsPositional is
            C{False} because positional arguments are always required. It may
            seem strange to allow no --samfile argument in a class that filters
            SAM files, but this option can be used in conjuction with scripts
            that can optionally take a SAM file (e.g.,
            genome-protein-summary.py).
        @param samfileNargs: The value to pass for 'nargs' in adding the
            samfile option.
        @param referenceIdRequired: If C{True}, make the --referenceId option
            required.
        """
        if samfileIsPositional:
            # Positional arguments are always required.
            assert samfileRequired, ('samfileIsPositional is True, so '
                                     'samfileRequired must also be True.')
            if samfileNargs is None:
                parser.add_argument(
                    'samfile', action=samfileAction,
                    help='The SAM/BAM file to filter.')
            else:
                parser.add_argument(
                    'samfile', action=samfileAction, nargs=samfileNargs,
                    help='The SAM/BAM file to filter.')
        else:
            if samfileNargs is None:
                parser.add_argument(
                    '--samfile', required=samfileRequired,
                    action=samfileAction, help='The SAM/BAM file to filter.')
            else:
                parser.add_argument(
                    '--samfile', required=samfileRequired, nargs=samfileNargs,
                    action=samfileAction, help='The SAM/BAM file to filter.')

        parser.add_argument(
            '--referenceId', metavar='ID', action='append',
            required=referenceIdRequired,
            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, metavar='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, metavar='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', metavar='TAG',
            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.
        """
        return cls(
            args.samfile,
            filterRead=filterRead,
            referenceIds=set(args.referenceId) if args.referenceId else None,
            storeQueryIds=storeQueryIds,
            dropUnmapped=args.dropUnmapped,
            dropSecondary=args.dropSecondary,
            dropSupplementary=args.dropSupplementary,
            dropDuplicates=args.dropDuplicates,
            keepQCFailures=args.keepQCFailures,
            minScore=args.minScore,
            maxScore=args.maxScore)

    def filterAlignment(self, alignment):
        """
        Test an alignment to see if it passes all filtering criteria.

        @param alignment: A pysam alignment instance.
        @return: A C{bool}, C{True} if the alignment passes our filtering,
            C{False} if it should be discarded.
        """
        if self.noFiltering:
            return True

        if self.minScore is not None or self.maxScore is not None:
            try:
                score = alignment.get_tag(self.scoreTag)
            except KeyError:
                return False
            else:
                if ((self.minScore is not None and score < self.minScore) or
                        (self.maxScore is not None and score > self.maxScore)):
                    return False

        return ((self.filterRead is None or
                 self.filterRead(Read(alignment.query_name,
                                      alignment.query_sequence,
                                      alignment.qual))) and
                not (
                    (self.referenceIds and
                     alignment.reference_name not in self.referenceIds) or
                    (alignment.is_unmapped and self.dropUnmapped) or
                    (alignment.is_secondary and self.dropSecondary) or
                    (alignment.is_supplementary and self.dropSupplementary) or
                    (alignment.is_duplicate and self.dropDuplicates) or
                    (alignment.is_qcfail and not self.keepQCFailures)))

    def alignments(self):
        """
        Get alignments from the SAM/BAM file, subject to filtering.

        @return: A generator that yields pysam alignment instances that pass
            our filtering criteria.
        """
        if self.storeQueryIds:
            queryIds = self.queryIds = set()
        else:
            queryIds = None

        lastAlignment = None
        count = 0
        with samfile(self.filename) as samAlignment:
            for count, alignment in enumerate(samAlignment.fetch(), start=1):
                if queryIds is not None:
                    queryIds.add(alignment.query_name)

                # Secondary and supplementary alignments may have a '*'
                # (pysam returns this as None) SEQ field, indicating that
                # the previous sequence should be used. This is best
                # practice according to section 2.5.2 of
                # https://samtools.github.io/hts-specs/SAMv1.pdf So we use
                # the last alignment query and quality strings if we get
                # None as a query sequence.
                if alignment.query_sequence is None:
                    if lastAlignment is None:
                        raise InvalidSAM(
                            'pysam produced an alignment (number %d) with no '
                            'query sequence without previously giving an '
                            'alignment with a sequence.' % count)
                    # Use the previous query sequence and quality. I'm not
                    # making the call to _hardClip dependent on
                    # alignment.cigartuples (as in the else clause below)
                    # because I don't think it's possible for
                    # alignment.cigartuples to be None in this case. If we
                    # have a second match on a query, then it must be
                    # aligned to something (i.e., it cannot be unmapped
                    # with no CIGAR string). The assertion will tell us if
                    # this is ever not the case.
                    assert alignment.cigartuples
                    (alignment.query_sequence,
                     alignment.query_qualities, _) = _hardClip(
                         lastAlignment.query_sequence,
                         lastAlignment.query_qualities,
                         alignment.cigartuples)
                else:
                    lastAlignment = alignment
                    if alignment.cigartuples:
                        (alignment.query_sequence,
                         alignment.query_qualities, _) = _hardClip(
                             alignment.query_sequence,
                             alignment.query_qualities,
                             alignment.cigartuples)

                if self.filterAlignment(alignment):
                    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 UnequalReferenceLengthError(SamError):
    "The references of interest in a SAM/BAM file are not of the same length." class UnknownReference(SamError):
    "Reference sequence not found in SAM/BAM file." class InvalidSAM(SamError):
    "SAM/BAM file has unexpected/invalid content." 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) def _hardClip(sequence, quality, cigartuples):
    """
    Hard clip (if necessary) a sequence.

    @param sequence: A C{str} nucleotide sequence.
    @param quality: A C{str} quality string, or a C{list} of C{int} quality
        values as returned by pysam, or C{None} if the SAM file had a '*'
        for the quality string (which pysam converts to C{None}).
    @param cigartuples: An iterable of (operation, length) tuples, detailing
        the alignment, as per the SAM specification.
    @return: A 3-tuple consisting of
            1) a hard-clipped C{str} sequence if hard-clipping is indicated by
               the CIGAR operations.
            2) a hard-clipped quality C{str} or C{list} (depending on what
               type we were passed) if hard-clipping is indicated by the CIGAR
               operations.
            3) a Boolean, C{True} if hard clipping was performed by this
               function or C{False} if the hard clipping had already been
               done.
    """
    hardClipCount = cigarLength = 0
    for (operation, length) in cigartuples:
        hardClipCount += operation == CHARD_CLIP
        cigarLength += length if operation in CONSUMES_QUERY else 0

    sequenceLength = len(sequence)
    if quality is not None:
        assert sequenceLength == len(quality)
    clipLeft = clipRight = 0
    clippedSequence = sequence
    clippedQuality = quality

    if sequenceLength > cigarLength:
        alreadyClipped = False
    else:
        assert sequenceLength == cigarLength
        alreadyClipped = True

    if hardClipCount == 0:
        pass
    elif hardClipCount == 1:
        # Hard clip either at the start or the end.
        if cigartuples[0][0] == CHARD_CLIP:
            if not alreadyClipped:
                clipLeft = cigartuples[0][1]
                clippedSequence = sequence[clipLeft:]
                if quality is not None:
                    clippedQuality = quality[clipLeft:]
        elif cigartuples[-1][0] == CHARD_CLIP:
            if not alreadyClipped:
                clipRight = cigartuples[-1][1]
                clippedSequence = sequence[:-clipRight]
                if quality is not None:
                    clippedQuality = quality[:-clipRight]
        else:
            raise ValueError(
                'Invalid CIGAR tuples (%s) contains hard-clipping operation '
                'that is neither at the start nor the end of the sequence.' %
                (cigartuples,))
    elif hardClipCount == 2:
        # Hard clip at both the start and end.
        assert cigartuples[0][0] == cigartuples[-1][0] == CHARD_CLIP
        if not alreadyClipped:
            clipLeft, clipRight = cigartuples[0][1], cigartuples[-1][1]
            clippedSequence = sequence[clipLeft:-clipRight]
            if quality is not None:
                clippedQuality = quality[clipLeft:-clipRight]
    else:
        raise ValueError(
            'Invalid CIGAR tuples (%s) specifies hard-clipping %d times (2 '
            'is the maximum).' % (cigartuples, hardClipCount))

    weClipped = bool(clipLeft or clipRight)

    if weClipped:
        assert not alreadyClipped
        if len(clippedSequence) + clipLeft + clipRight != sequenceLength:
            raise ValueError(
                'Sequence %r (length %d) clipped to %r (length %d), but the '
                'difference between these two lengths (%d) is not equal to '
                'the sum (%d) of the left and right clip lengths (%d and %d '
                'respectively). CIGAR tuples: %s' %
                (sequence, len(sequence),
                 clippedSequence, len(clippedSequence),
                 abs(len(sequence) - len(clippedSequence)),
                 clipLeft + clipRight, clipLeft, clipRight, cigartuples))
    else:
        assert len(clippedSequence) == sequenceLength
        if quality is not None:
            assert len(clippedQuality) == sequenceLength

    return clippedSequence, clippedQuality, weClipped class DistanceMatrix:
    """
    Maintain score (matching) information between reads and reference sequences
    as read from a SAM file (or files) and provide methods for calculating
    distance between references.
    """
    def __init__(self):
        # self.scores is a similarity matrix between references and
        # queries from a SAM file. The values are 1.0 or 0.0 when queries
        # match / do not match references unless we are using alignment
        # scores, as generated by tools like bowtie2 and bwa and stored in
        # a SAM tag. In that case, the alignment score value is stored (see
        # https://en.wikipedia.org/wiki/SAM_(file_format)#Optional_fields
        # for SAM tag info).
        #
        # The dictionary is keyed by str reference ids and the values are
        # dicts keyed by str query ids. The values of the sub-dictionaries
        # are the float scores (similarities). E.g.
        #
        # {
        #     'reference-1': {
        #         'query-1': 77.5,
        #         'query-2': 37.0,
        #     },
        #     'reference-2': {
        #         'query-1': 78.0,
        #     },
        # }
        #
        # If a reference/query id pair is not present in self.scores, it
        # means that query did not match that reference. The score method
        # (below) will assign such a pair a similarity score of 0.0.
        self.scores = {}

    def addFile(self, filename, scoreTag=None):
        """
        Add information from a SAM file.

        @param filename: A C{str} SAM file name to read.
        @param scoreTag: A C{str} alignment score tag to use instead of binary
            match / non-match values that are used if C{scoreTag} is C{None}.
            The score tag must be present in the optional SAM fields.
        """
        if scoreTag:
            def getScore(alignment, filename, count):
                """
                Get the alignment score from a SAM alignment.

                @param alignment: A samtools alignment.
                @param filename: The C{str} SAM file name being read.
                @param count: The C{int} line number in the SAM file.
                """
                try:
                    score = alignment.get_tag(scoreTag)
                except KeyError:
                    raise ValueError(
                        f'Alignment {count} in {filename!r} has no '
                        f'{scoreTag!r} score tag.')
                else:
                    if score < 0.0:
                        raise ValueError(
                            f'Alignment {count} in {filename!r} has tag '
                            f'{scoreTag!r} with negative value ({score}).')
                    return score
        else:
            def getScore(*args):
                """
                Return a binary score of 1.0 seeing as we know the query was
                mapped and no score tag is in use.

                @param args: ignored.
                """
                return 1.0

        scores = self.scores

        with samfile(filename) as samAlignment:
            for count, alignment in enumerate(samAlignment.fetch(), start=1):
                if alignment.is_unmapped:
                    continue

                queryId = alignment.query_name
                referenceId = alignment.reference_name
                score = getScore(alignment, filename, count)

                try:
                    preExisting = scores[referenceId][queryId]
                except KeyError:
                    if referenceId not in scores:
                        scores[referenceId] = {}

                    scores[referenceId][queryId] = score
                else:
                    if score > preExisting:
                        scores[referenceId][queryId] = score

    def score(self, referenceId, queryId):
        """
        Get the score for a reference and query.

        Although this method name does not start with an underscore, it's not
        likely to be used directly when this class is used to calculate
        distances between references according to query matching. The scores
        between references and queries (returned by this function) are just
        a component of the overall distance.

        @param referenceId: A C{str} reference id.
        @param queryId: A C{str} query id.
        @return: A C{float} score (zero if there was no match between the query
            and reference) or if the reference id or query id have never been
            seen. This could arguably cause a KeyError, but I decided against
            that to keep the code simpler / faster. If you pass in unknown
            query or reference ids, that's your problem (this should be very
            unlikely because our caller is operating on a SAM file, so the
            place they're most likely to be getting query and reference ids
            from to pass us is from that same SAM file.
        """
        try:
            return self.scores[referenceId][queryId]
        except KeyError:
            return 0.0

    def jaccardDistance(self, referenceId1, referenceId2):
        """
        Get the Jaccard distance between two references.  See
        https://en.wikipedia.org/wiki/Jaccard_index

        @param referenceId1: A C{str} reference id.
        @param referenceId2: A C{str} reference id.
        @return: A C{float} distance.
        """
        queryIds1 = set(self.scores.get(referenceId1, []))
        queryIds2 = set(self.scores.get(referenceId2, []))

        denominator = len(queryIds1 | queryIds2)

        if denominator:
            numerator = len(queryIds1 & queryIds2)
            similarity = numerator / denominator
            distance = 1.0 - similarity
            return distance
        else:
            # The references have no matching queries in common.
            return 1.0

    def soergelDistance(self, referenceId1, referenceId2):
        """
        Get the Soergel (i.e., weighted Jaccard) distance between two
        references. See https://en.wikipedia.org/wiki/Jaccard_index

        Note that this will be the same as the jaccardDistance unless a
        scoreTag was used in reading the SAM file.

        @param referenceId1: A C{str} reference id.
        @param referenceId2: A C{str} reference id.
        @return: A C{float} distance.
        """
        queryIds1 = set(self.scores.get(referenceId1, []))
        queryIds2 = set(self.scores.get(referenceId2, []))

        score = self.score
        scores = [
            (score(referenceId1, queryId), score(referenceId2, queryId))
            for queryId in (queryIds1 | queryIds2)]

        denominator = sum(max(a, b) for a, b in scores)

        if denominator:
            numerator = sum(min(a, b) for a, b in scores)
            similarity = numerator / denominator
            distance = 1.0 - similarity
            return distance
        else:
            # Neither reference had a positive match against any query.
            return 1.0

    def matrix(self, referenceIds=None, metric='soergel', similarity=False,
               returnDict=False):
        """
        Compute a distance matrix.

        @param referenceIds: An iterable of C{str} reference ids. If C{None},
            all known reference ids will be used.
        @param metric: A C{str}, either 'soergel' or 'jaccard'.
        @param similarity: If C{True}, return a similarity matrix.
        @param returnDict: If C{True}, return a C{dict} (see @return, below).
        @return: If C{returnDict} is C{False}, return a symmetric square
            C{np.array} of C{float} inter-reference distances (or similarities
            if C{similarity} is C{True}) in the range (0.0, 1.0). The order of
            the entries in the rows and columns of the returned array matches
            that of C{referenceIds} (if given) or the order in C{self.scores}
            otherwise.  If C{returnDict} is C{True}, return a C{dict} of
            C{dict}s, indexed by the two reference ids, with values as above.
        """
        assert metric in {'jaccard', 'soergel'}
        referenceIds = tuple(referenceIds or self.scores)
        nIds = len(referenceIds)
        matrix = defaultdict(dict) if returnDict else np.empty((nIds, nIds))
        identityValue = 1.0 if similarity else 0.0
        func = (self.jaccardDistance if metric == 'jaccard' else
                self.soergelDistance)

        if similarity:
            def metricFunc(ref1, ref2):
                return 1.0 - func(ref1, ref2)
        else:
            def metricFunc(ref1, ref2):
                return func(ref1, ref2)

        for index1, ref1 in enumerate(referenceIds):
            for index2, ref2 in enumerate(referenceIds):
                i, j = (ref1, ref2) if returnDict else (index1, index2)
                value = identityValue if i == j else metricFunc(ref1, ref2)
                matrix[i][j] = matrix[j][i] = value

        return matrix

    def save(self, fp):
        """
        Save the similarity matrix as JSON.

        @param fp: An open file pointer to write to.
        """
        dump(self.scores, fp, sort_keys=True, indent=4)

    def load(self, fp):
        """
        Load a similarity matrix from a JSON file.

        @param fp: An open file pointer to read from.
        """
        self.scores = load(fp)
from six import assertRaisesRegex
from unittest import TestCase
from tempfile import mkstemp
from os import close, unlink, write
from contextlib import contextmanager
from io import StringIO
from json import dumps, loads
from pysam import CHARD_CLIP, CMATCH
from dark.reads import Read, ReadFilter
from dark.sam import (
    PaddedSAM, SAMFilter, UnequalReferenceLengthError, UnknownReference,
    InvalidSAM, samReferencesToStr, _hardClip, DistanceMatrix)
import numpy as np
        """
        data = '\n'.join([
            '@SQ SN:id1 LN:90',
            '@SQ SN:id2 LN:91',
        ]).replace(' ', '\t')

        with dataFile(data) as filename:
            self.assertEqual('  id1 (length 90)\n  id2 (length 91)',
                             samReferencesToStr(filename, indent=2))


class TestHardClip(TestCase):
    """
    Test the _hardClip function.
    """
    def testHardClipInMiddle(self):
        """
        If hard clipping is given as an operation not at the beginning or end
        of the sequence, a ValueError must be raised.
        """
        error = (
            '^Invalid CIGAR tuples .* contains hard-clipping operation '
            'that is neither at the start nor the end of the sequence\\.$')
        assertRaisesRegex(
            self, ValueError, error, _hardClip, 'CGT', '123',
            ((CMATCH, 1), (CHARD_CLIP, 1), (CMATCH, 1),))

    def testThreeHardClips(self):
        """
        If hard clipping is specified more than twice, a ValueError must be
        raised.
        """
        error = ('^Invalid CIGAR tuples .* specifies hard-clipping 3 times '
                 '\\(2 is the maximum\\).$')
        assertRaisesRegex(
            self, ValueError, error, _hardClip, 'CGT', '123',
            ((CHARD_CLIP, 1), (CHARD_CLIP, 1), (CHARD_CLIP, 1),))

    def testNoClip(self):
        """
        If no hard clipping is indicated, the function must return the
        original sequence.
        """
        self.assertEqual(('CGT', '123', False),
                         _hardClip('CGT', '123', ((CMATCH, 3),)))

    def testClipLeft(self):
        """
        If hard clipping on the left is indicated, and has not been done,
        the function must return the expected sequence.
        """
        self.assertEqual(
            ('CGT', '456', True),
            _hardClip('CAACGT', '123456', ((CHARD_CLIP, 3), (CMATCH, 3),)))

    def testClipRight(self):
        """
        If hard clipping on the right is indicated, and has not been done,
        the function must return the expected sequence.
        """
        self.assertEqual(
            ('CA', '12', True),
            _hardClip('CAACGT', '123456', ((CMATCH, 2), (CHARD_CLIP, 4),)))

    def testClipBoth(self):
        """
        If hard clipping on the left and right is indicated, and has not been
        done, the function must return the expected sequence.
        """
        self.assertEqual(
            ('AA', '23', True),
            _hardClip('CAACGT', '123456',
                      ((CHARD_CLIP, 1), (CMATCH, 2), (CHARD_CLIP, 3),)))

    def testClipLeftAlreadyDone(self):
        """
        If hard clipping on the left is indicated, and has already been done,
        the function must return the expected sequence.
        """
        self.assertEqual(
            ('CGT', '123', False),
            _hardClip('CGT', '123', ((CHARD_CLIP, 3), (CMATCH, 3),)))

    def testClipRightAlreadyDone(self):
        """
        If hard clipping on the right is indicated, and has already been done,
        the function must return the expected sequence.
        """
        self.assertEqual(
            ('CA', '12', False),
            _hardClip('CA', '12', ((CMATCH, 2), (CHARD_CLIP, 4),)))

    def testClipBothAlreadyDone(self):
        """
        If hard clipping on the left and right is indicated, and has already
        been done, the function must return the expected sequence.
        """
        self.assertEqual(
            ('AA', '12', False),
            _hardClip('AA', '12',
                      ((CHARD_CLIP, 1), (CMATCH, 2), (CHARD_CLIP, 3),)))


class TestDistanceMatrix(TestCase):
    """
    Test the DistanceMatrix class.
    """
    def testEmpty(self):
        """
        The similarity (scores) matrix must be empty after intitialization.
        """
        dm = DistanceMatrix()
        self.assertEqual({}, dm.scores)

    def testNoMatchesHasNoMatchedReferenceIds(self):
        """
        If a SAM file with no query records is added, there should be no
        matched reference ids (even though there is a reference in the SAM
        header, it has no matches).
        """

        dm = DistanceMatrix()