<gh_stars>0
# coding: utf-8
"""
Adapters
"""
import itertools
import logging
import os
import pickle
import re
import sys
from urllib.error import URLError
from urllib.request import urlopen
from atropos import align
from atropos.align import Match
from atropos.io.seqio import ColorspaceSequence, FastaReader
from atropos.util import (
    IUPAC_BASES, GC_BASES, MergingDict, NestedDict, CountingDict, Const, 
    reverse_complement, ALPHABETS)
from atropos.util import colorspace as cs

class AdapterType(object):
    """Struct for adapter type information.
    
    Args:
        name: Adapter type name.
        desc: Adapter type description.
        flags: Alignment flags.
    """
    def __init__(self, name, desc, *flags):
        self.name = name
        self.desc = desc
        self.flags = flags[0]
        for i in range(1, len(flags)):
            self.flags |= flags[i]
    
    def asdict(self):
        """Returns AdapterType fields in a dict.
        """
        return dict(name=self.name, desc=self.desc, flags=Const(self.flags))

ADAPTER_TYPES = dict(
    back=AdapterType(
        'back', "regular 3'",
        align.START_WITHIN_SEQ2,
        align.STOP_WITHIN_SEQ2,
        align.STOP_WITHIN_SEQ1),
    front=AdapterType(
        'front', "regular 5'",
        align.START_WITHIN_SEQ2,
        align.STOP_WITHIN_SEQ2,
        align.START_WITHIN_SEQ1),
    prefix=AdapterType('prefix', "anchored 5'", align.STOP_WITHIN_SEQ2),
    suffix=AdapterType('suffix', "anchored 3'", align.START_WITHIN_SEQ2),
    anywhere=AdapterType('anywhere', "variable 5'/3'", align.SEMIGLOBAL),
    linked=AdapterType('linked', 'linked', 'linked')
)

def where_int_to_dict(where):
    """Convert a "where" flags integer to a dictionary of values from the
    corresponding AdapterType.
    """
    for adapter_type in ADAPTER_TYPES.values():
        if where == adapter_type.flags:
            return adapter_type.asdict()
    raise ValueError("Invalid WHERE value: {}".format(where))

# TODO get rid of these constants
BACK = ADAPTER_TYPES['back'].flags
FRONT = ADAPTER_TYPES['front'].flags
PREFIX = ADAPTER_TYPES['prefix'].flags
SUFFIX = ADAPTER_TYPES['suffix'].flags
ANYWHERE = ADAPTER_TYPES['anywhere'].flags
LINKED = ADAPTER_TYPES['linked'].flags

# TODO: specify this externally rather than hard-coding
DEFAULT_ADAPTERS_URL = "https://raw.githubusercontent.com/jdidion/atropos/master/atropos/adapters/sequencing_adapters.fa"
DEFAULT_ADAPTERS_PATH = os.path.join(
    os.path.dirname(__file__), 'sequencing_adapters.fa')

class AdapterParser(object):
    """Factory for Adapter classes that all use the same parameters (error rate,
    indels etc.). The given **kwargs will be passed to the Adapter constructors.
    
    Args:
        colorspace: Whether the adapter sequences are in colorspace.
        cache: The AdapterCache to use for fetching/storing named adapters.
    """
    def __init__(self, colorspace=False, cache=None, **kwargs):
        self.colorspace = colorspace
        self.cache = cache
        self.constructor_args = kwargs
        self.adapter_class = ColorspaceAdapter if colorspace else Adapter
    
    def parse(self, spec, cmdline_type='back'):
        """Parse an adapter specification not using ``file:`` notation and return
        an object of an appropriate Adapter class. The notation for anchored
        5' and 3' adapters is supported. If the name parameter is None, then
        an attempt is made to extract the name from the specification
        (If spec is 'name=ADAPTER', name will be 'name'.)
        
        Args:
            spec: The adapter spec.
            name: The adapter name. If not provided, one is automatically
                generated.
            cmdline_type: describes which commandline parameter was used (``-a``
                is 'back', ``-b`` is 'anywhere', and ``-g`` is 'front').
        
        TODO: describe the adapter spec format
        
        Returns:
            An :class:`Adapter` instance.
        """
        if spec.startswith('file:'):
            # read adapter sequences from a file
            with FastaReader(spec[5:]) as fasta:
                for record in fasta:
                    name = record.name.split(None, 1)[0]
                    yield self.parse_from_spec(
                        record.sequence, cmdline_type, name)
        else:
            yield self.parse_from_spec(spec, cmdline_type)
    
    def parse_from_spec(self, spec, cmdline_type='back', name=None):
        if cmdline_type not in ADAPTER_TYPES:
            raise ValueError(
                'cmdline_type cannot be {0!r}'.format(cmdline_type))
        orig_spec = spec
        where = ADAPTER_TYPES[cmdline_type].flags
        
        if name is None and spec is None:
            raise ValueError('Either name or spec must be given')
        elif name is None:
            if self.cache and self.cache.has_name(spec):
                name = spec
                spec = self.cache.get_for_name(name)
        elif spec is None:
            if self.cache and self.cache.has_name(name):
                spec = self.cache.get_for_name(name)
        
        if spec is None:
            raise ValueError('Name not found: {}'.format(name))
        elif name is None:
            name, spec = _extract_name_from_spec(spec)
        
        if self.cache and name is not None:
            self.cache.add(name, spec)
        
        front_anchored, back_anchored = False, False
        if spec.startswith('^'):
            spec = spec[1:]
            front_anchored = True
        if spec.endswith('$'):
            spec = spec[:-1]
            back_anchored = True
        
        sequence1, middle, sequence2 = spec.partition('...')
        
        if where == ANYWHERE:
            if front_anchored or back_anchored:
                raise ValueError("'anywhere' (-b) adapters may not be anchored")
            if middle == '...':
                raise ValueError("'anywhere' (-b) adapters may not be linked")
            return self.adapter_class(
                sequence=spec, where=where, name=name, **self.constructor_args)
        
        assert where == FRONT or where == BACK
        if middle == '...':
            if not sequence1:
                if where == BACK:  # -a ...ADAPTER
                    spec = sequence2
                else:  # -g ...ADAPTER
                    raise ValueError('Invalid adapter specification')
            elif not sequence2:
                if where == BACK:  # -a ADAPTER...
                    spec = sequence1
                    where = FRONT
                    front_anchored = True
                else:  # -g ADAPTER...
                    spec = sequence1
            else:
                # linked adapter
                if self.colorspace:
                    raise NotImplementedError(
                        'Using linked adapters in colorspace is not supported')
                # automatically anchor 5' adapter if -a is used
                if where == BACK:
                    front_anchored = True
                
                return LinkedAdapter(sequence1, sequence2, name=name,
                    front_anchored=front_anchored, back_anchored=back_anchored,
                    **self.constructor_args)
        
        if front_anchored and back_anchored:
            raise ValueError(
                'Trying to use both "^" and "$" in adapter specification '
                '{!r}'.format(orig_spec))
        if front_anchored:
            if where == BACK:
                raise ValueError("Cannot anchor the 3' adapter at its 5' end")
            where = PREFIX
        elif back_anchored:
            if where == FRONT:
                raise ValueError("Cannot anchor 5' adapter at 3' end")
            where = SUFFIX
        
        return self.adapter_class(
            sequence=spec, where=where, name=name, **self.constructor_args)
    
    def parse_multi(self, back=None, anywhere=None, front=None):
        """Parse all three types of commandline options that can be used to
        specify adapters. back, anywhere and front are lists of strings,
        corresponding to the respective commandline types (-a, -b, -g).
        
        Args:
            back: Back-adapter specs.
            anywhere: Anywhere-adapter specs.
            front: Front-adapter specs.
        
        Return:
            A list of appropriate Adapter classes.
        """
        adapters = []
        for specs, cmdline_type in (
                (back, 'back'), (anywhere, 'anywhere'), (front, 'front')):
            if not specs:
                continue
            for spec in specs:
                adapters.extend(self.parse(spec, cmdline_type))
        return adapters

class Adapter(object):
    """An adapter knows how to match itself to a read. In particular, it knows
    where it should be within the read and how to interpret wildcard characters.
    
    Args:
        sequence: The adapter sequence as string. Will be converted to
            uppercase. Also, Us will be converted to Ts.
        where: One of the BACK, FRONT, PREFIX, SUFFIX or ANYWHERE constants.
            This influences where the adapter is allowed to appear within in the
            read and also which part of the read is removed.
        max_error_rate: Maximum allowed error rate. The error rate is the number
            of errors in the alignment divided by the length of the part of the
            alignment that matches the adapter.
        min_overlap: Minimum length of the part of the alignment that
            matches the adapter.
        read_wildcards: Whether IUPAC wildcards in the read are allowed.
        adapter_wildcards: Whether IUPAC wildcards in the adapter are allowed.
        name: optional name of the adapter. If not provided, the name is set to
            a unique number.
        indels: Whether indels are allowed in adapter matches.
        indel_cost: Cost of indel bases in alignment scoring.
        match_probability: Callabale that computes the probability that a match
            could happen by random chance.
        max_rmp: Maximum random-match probability below which a match is
            considered genuine.
        gc_content: Expected GC content of sequences.
        alphabet: The Alphabet to use to validate the adapter sequence.
    """
    def __init__(
            self, sequence, where, max_error_rate=0.1, min_overlap=3,
            read_wildcards=False, adapter_wildcards=True, name=None,
            indels=True, indel_cost=1, match_probability=None, max_rmp=None,
            gc_content=0.5, alphabet=None):
        if len(sequence) == 0:
            raise ValueError("Empty adapter sequence")
        # TODO: all of this validation code should be wrapped up in Alphabet
        sequence = parse_braces(sequence.upper().replace('U', 'T'))
        seq_set = set(sequence)
        if seq_set <= set('ACGT'):
            adapter_wildcards = False
        if adapter_wildcards and not seq_set <= IUPAC_BASES:
            raise ValueError(
                "Invalid character(s) in adapter sequence: {}".format(
                    ','.join(seq_set - IUPAC_BASES)))
        if alphabet:
            if isinstance(alphabet, str):
                alphabet = ALPHABETS[alphabet]
            alphabet.validate_string(sequence)
        
        self.debug = False
        self.name = _generate_adapter_name() if name is None else name
        self.sequence = sequence
        self.where = where
        self.max_error_rate = max_error_rate
        self.min_overlap = min(min_overlap, len(self.sequence))
        self.match_probability = match_probability
        self.max_rmp = max_rmp
        self.gc_content = gc_content
        self.indels = indels
        self.adapter_wildcards = adapter_wildcards
        self.read_wildcards = read_wildcards
        # redirect trimmed() to appropriate function depending on adapter type
        trimmers = {
            FRONT: self._trimmed_front,
            PREFIX: self._trimmed_front,
            BACK: self._trimmed_back,
            SUFFIX: self._trimmed_back,
            ANYWHERE: self._trimmed_anywhere
        }
        self.trimmed = trimmers[where]
        if where == ANYWHERE:
            self._front_flag = None  # means: guess
        else:
            self._front_flag = where not in (BACK, SUFFIX)
        # statistics about length of removed sequences
        self.lengths_front = CountingDict()
        self.lengths_back = CountingDict()
        self.errors_front = NestedDict()
        self.errors_back = NestedDict()
        self.adjacent_bases = { 'A': 0, 'C': 0, 'G': 0, 'T': 0, '': 0 }
        self.aligner = align.Aligner(
            self.sequence, self.max_error_rate, flags=self.where,
            wildcard_ref=self.adapter_wildcards,
            wildcard_query=self.read_wildcards)
        self.aligner.min_overlap = self.min_overlap
        if self.indels:
            self.aligner.indel_cost = indel_cost
        else:
            # TODO
            # When indels are disallowed, an entirely different algorithm
            # should be used.
            self.aligner.indel_cost = 100000
    
    def __repr__(self):
        return '<Adapter(name="{name}", sequence="{sequence}", where={where}, '\
            'max_error_rate={max_error_rate}, min_overlap={min_overlap}, '\
            'read_wildcards={read_wildcards}, '\
            'adapter_wildcards={adapter_wildcards}, '\
            'indels={indels})>'.format(**vars(self))
    
    def enable_debug(self):
        """Print out the dynamic programming matrix after matching a read to an
        adapter.
        """
        self.debug = True
        self.aligner.enable_debug()
    
    def match_to(self, read):
        """Attempt to match this adapter to the given read.
        
        Args:
            read: A :class:`Sequence` instance.
        
        Returns:
            A :class:`Match` instance if a match was found; return None if no
            match was found given the matching criteria (minimum overlap length,
            maximum error rate).
        """
        read_seq = read.sequence.upper()
        
        # try to find an exact match first unless wildcards are allowed
        pos = -1
        if not self.adapter_wildcards:
            if self.where == PREFIX:
                if read_seq.startswith(self.sequence):
                    pos = 0
            elif self.where == SUFFIX:
                if read_seq.endswith(self.sequence):
                    pos = (len(read_seq) - len(self.sequence))
            else:
                pos = read_seq.find(self.sequence)
        
        if pos >= 0:
            seqlen = len(self.sequence)
            return Match(
                0, seqlen, pos, pos + seqlen, seqlen, 0, self._front_flag,
                self, read)
        
        # try approximate matching
        if not self.indels and self.where in (PREFIX, SUFFIX):
            if self.where == PREFIX:
                alignment = align.compare_prefixes(
                    self.sequence, read_seq,
                    wildcard_ref=self.adapter_wildcards,
                    wildcard_query=self.read_wildcards)
            else:
                alignment = align.compare_suffixes(
                    self.sequence, read_seq,
                    wildcard_ref=self.adapter_wildcards,
                    wildcard_query=self.read_wildcards)
        else:
            alignment = self.aligner.locate(read_seq)
            if self.debug:
                print(self.aligner.dpmatrix)  # pragma: no cover
        
        if alignment:
            astart, astop, rstart, rstop, matches, errors = alignment
            size = astop - astart
            if ((
                    size >=
                    self.min_overlap and errors / size <=
                    self.max_error_rate
                ) and (
                    self.max_rmp is None or
                    self.match_probability(matches, size) <= self.max_rmp)):
                return Match(
                    astart, astop, rstart, rstop, matches, errors,
                    self._front_flag, self, read)
        
        return None
    
    def _trimmed_anywhere(self, match):
        """Trims an adapter from either the front or back of sequence.
        
        Returns:
            A :class:`Sequence` instance: the trimmed read.
        """
        if match.front:
            return self._trimmed_front(match)
        else:
            return self._trimmed_back(match)
    
    def _trimmed_front(self, match):
        """Trims an adapter from the front of sequence.
        
        Returns:
            A :class:`Sequence` instance: the trimmed read.
        """
        # TODO move away
        self.lengths_front[match.rstop] += 1
        self.errors_front[match.rstop][match.errors] += 1
        return match.read[match.rstop:]
    
    def _trimmed_back(self, match):
        """Trims an adapter from the back of sequence.
        
        Returns:
            A :class:`Sequence` instance: the trimmed read.
        """
        self.lengths_back[len(match.read) - match.rstart] += 1
        self.errors_back[len(match.read) - match.rstart][match.errors] += 1
        adjacent_base = match.read.sequence[match.rstart-1:match.rstart]
        if adjacent_base not in 'ACGT':
            adjacent_base = ''
        self.adjacent_bases[adjacent_base] += 1
        return match.read[:match.rstart]
    
    def __len__(self):
        return len(self.sequence)
    
    def random_match_probabilities(self):
        """Estimate probabilities that this adapter matches a random sequence. 
        Indels are not taken into account.

        Returns:
            A list of probabilities the same length as this adapter's sequence,
            where the value at position 'i' is the probability that i bases of
            this adapter match a random sequence.
        """
        if self._front_flag:
            seq = self.sequence[::-1]
        else:
            seq = self.sequence
        
        #matches = 0
        base_probs = (self.gc_content / 2.0, (1 - self.gc_content) / 2.0)
        probabilities = [1.0] + ([0] * len(seq))
        c_bases = frozenset(GC_BASES if self.adapter_wildcards else 'GC')
        
        # TODO: this doesn't work; need to figure out if RandomMatchProbability
        # can be used for this.
        #for idx, base in enumerate(seq, 1):
        #    if base in c_bases:
        #        matches += 1
        #    probabilities[idx] = self.match_probability(
        #        matches, idx, *base_probs)
        
        cur_p = 1.0
        for idx, base in enumerate(seq, 1):
            cur_p *= base_probs[0 if base in c_bases else 1]
            probabilities[idx] = cur_p
        return probabilities
    
    def summarize(self):
        """Summarize the activities of this :class:`Adapter`.
        """
        total_front = sum(self.lengths_front.values())
        total_back = sum(self.lengths_back.values())
        
        stats = MergingDict(
            adapter_class=self.__class__.__name__,
            total_front=total_front,
            total_back=total_back,
            total=total_front + total_back,
            match_probabilities=Const(self.random_match_probabilities()))
        
        where = self.where
        assert (
            where in (ANYWHERE, LINKED) or
            (where in (BACK, SUFFIX) and total_front == 0) or
            (where in (FRONT, PREFIX) and total_back == 0))
        
        stats["where"] = where_int_to_dict(where)
        stats["sequence"] = Const(self.sequence)
        stats["max_error_rate"] = Const(self.max_error_rate)
        if where in (ANYWHERE, FRONT, PREFIX):
            stats["lengths_front"] = self.lengths_front
            stats["errors_front"] = self.errors_front
        if where in (ANYWHERE, BACK, SUFFIX):
            stats["lengths_back"] = self.lengths_back
            stats["errors_back"] = self.errors_back
        if where in (BACK, SUFFIX):
            stats["adjacent_bases"] = self.adjacent_bases
        
        return stats

class ColorspaceAdapter(Adapter):
    """An adapter for a colorspace sequence.
    
    Args:
        args, kwargs: Arguments to pass to :class:`Adapter` constructor.
    """
    def __init__(self, *args, **kwargs):
        if kwargs.get('adapter_wildcards', False):
            raise ValueError("Wildcards not supported for colorspace adapters")
        kwargs['adapter_wildcards'] = False
        super().__init__(*args, **kwargs)
        has_nucleotide_seq = False
        if set(self.sequence) <= set('ACGT'):
            # adapter was given in basespace
            self.nucleotide_sequence = self.sequence
            has_nucleotide_seq = True
            self.sequence = cs.encode(self.sequence)[1:]
        if self.where in (PREFIX, FRONT) and not has_nucleotide_seq:
            raise ValueError(
                "A 5' colorspace adapter needs to be given in nucleotide space")
        self.aligner.reference = self.sequence
    
    def match_to(self, read):
        """Attempt to match this adapter to the given read.
        
        Args:
            read: A :class:`Sequence` instance.
        
        Returns:
            A :class:`Match` instance if a match was found; return None if no
            match was found given the matching criteria (minimum overlap length,
            maximum error rate).
        """
        if self.where != PREFIX:
            return super().match_to(read)
        # create artificial adapter that includes a first color that encodes the
        # transition from primer base into adapter
        asequence = (
            cs.ENCODE[read.primer + self.nucleotide_sequence[0:1]] +
            self.sequence)
        
        pos = 0 if read.sequence.startswith(asequence) else -1
        if pos >= 0:
            match = Match(
                0, len(asequence), pos, pos + len(asequence),
                len(asequence), 0, self._front_flag, self, read)
        else:
            # try approximate matching
            self.aligner.reference = asequence
            alignment = self.aligner.locate(read.sequence)
            if self.debug:
                print(self.aligner.dpmatrix)  # pragma: no cover
            if alignment is not None:
                match = Match(*(alignment + (self._front_flag, self, read)))
            else:
                match = None
        
        if match is None:
            return None
        assert (
            match.length > 0 and
            match.errors / match.length <= self.max_error_rate)
        assert match.length >= self.min_overlap
        return match

    def _trimmed_front(self, match):
        """Trims an adapter from the front of sequence.
        
        Returns:
            A :class:`Sequence` instance: the trimmed read.
        """
        read = match.read
        self.lengths_front[match.rstop] += 1
        self.errors_front[match.rstop][match.errors] += 1
        # to remove a front adapter, we need to re-encode the first color
        # following the adapter match
        color_after_adapter = read.sequence[match.rstop:match.rstop + 1]
        if not color_after_adapter:
            # the read is empty
            return read[match.rstop:]
        base_after_adapter = cs.DECODE[
            self.nucleotide_sequence[-1:] + color_after_adapter]
        new_first_color = cs.ENCODE[read.primer + base_after_adapter]
        new_read = read[:]
        new_read.sequence = new_first_color + read.sequence[(match.rstop + 1):]
        new_read.qualities = None
        if read.qualities:
            new_read.qualities = read.qualities[match.rstop:]
        return new_read

    def _trimmed_back(self, match):
        """Trims an adapter from the back of sequence.
        
        Returns:
            A :class:`Sequence` instance: the trimmed read.
        """
        # trim one more color if long enough
        adjusted_rstart = max(match.rstart - 1, 0)
        self.lengths_back[len(match.read) - adjusted_rstart] += 1
        self.errors_back[len(match.read) - adjusted_rstart][match.errors] += 1
        return match.read[:adjusted_rstart]

    def __repr__(self):
        return '<ColorspaceAdapter(sequence={0!r}, where={1})>'.format(
            self.sequence, self.where)

# TODO Consolidate Match and a LinkedMatch class.

class LinkedMatch(object):
    """Represent a match of a LinkedAdapter.
    
    Args:
        front_match: The match to the front of the sequence.
        back_match: The match to the back of the sequence.
        adapter: The matched adapter.
    """
    def __init__(self, front_match, back_match, adapter):
        self.front_match = front_match
        self.back_match = back_match
        self.adapter = adapter
        assert front_match is not None
    
    def get_info_record(self):
        """Returns the info record for the either the back or forward match.
        """
        if self.back_match:
            return self.back_match.get_info_record()
        else:
            return self.front_match.get_info_record()

class LinkedAdapter(object):
    """An adapter with linked front and back sequences.
    
    Args:
        front_sequence: Front adapter sequence.
        back_sequence: Back adapter sequence.
        front_anchored: Whether the front adapter is anchored.
        back_anchored: Whether the back adapter is anchored.
        name: Adapter name.
        kwargs: passed on to individual :class:`Adapter` constructors.
    """
    def __init__(
            self, front_sequence, back_sequence, front_anchored=True,
            back_anchored=False, name=None, **kwargs):
        assert front_anchored and not back_anchored
        where1 = PREFIX if front_anchored else FRONT
        where2 = SUFFIX if back_anchored else BACK
        self.front_anchored = front_anchored
        self.back_anchored = back_anchored
        
        # The following attributes are needed for the report
        self.where = LINKED
        self.name = _generate_adapter_name() if name is None else name
        self.front_adapter = Adapter(
            front_sequence, where=where1, name=None, **kwargs)
        self.back_adapter = Adapter(
            back_sequence, where=where2, name=None, **kwargs)
    
    def enable_debug(self):
        """Enable debug on adapters.
        """
        self.front_adapter.enable_debug()
        self.back_adapter.enable_debug()
    
    def match_to(self, read):
        """Match the linked adapters against the given read. If the 'front'
        adapter is not found, the 'back' adapter is not searched for.
        
        Args:
            read: A :class:`Sequence` instance.
        
        Returns:
            A :class:`Match` instance if a match was found; return None if no
            match was found given the matching criteria (minimum overlap length,
            maximum error rate).
        """
        front_match = self.front_adapter.match_to(read)
        if front_match is None:
            return None
        # TODO use match.trimmed() instead as soon as that does not update
        # statistics anymore
        read = read[front_match.rstop:]
        back_match = self.back_adapter.match_to(read)
        return LinkedMatch(front_match, back_match, self)

    def trimmed(self, match):
        """Returns the read trimmed with the front and/or back adapter
        trimmer(s).
        
        Args:
            match: The match to trim.
        
        Returns:
            The trimmed read.
        """
        front_trimmed = self.front_adapter.trimmed(match.front_match)
        if match.back_match:
            return self.back_adapter.trimmed(match.back_match)
        else:
            return front_trimmed
    
    def summarize(self):
        """Returns the summary dict for this adapter.
        """
        total_front = sum(self.front_adapter.lengths_front.values())
        total_back = sum(self.back_adapter.lengths_back.values())
        
        stats = MergingDict(
            total_front=total_front,
            total_back=total_back,
            total=total_front + total_back)
        
        where = self.where
        assert (
            where in (ANYWHERE, LINKED) or
            (where in (BACK, SUFFIX) and total_front == 0) or
            (where in (FRONT, PREFIX) and total_back == 0))
        
        stats["where"] = where_int_to_dict(where)
        stats["front_sequence"] = Const(self.front_adapter.sequence)
        stats["front_match_probabilities"] = Const(
            self.front_adapter.random_match_probabilities())
        stats["back_sequence"] = Const(self.back_adapter.sequence)
        stats["back_match_probabilities"] = Const(
            self.back_adapter.random_match_probabilities())
        stats["front_max_error_rate"] = Const(self.front_adapter.max_error_rate)
        stats["back_max_error_rate"] = Const(self.back_adapter.max_error_rate)
        stats["front_lengths_front"] = self.front_adapter.lengths_front
        stats["front_lengths_back"] = self.front_adapter.lengths_back
        stats["back_lengths_front"] = self.back_adapter.lengths_front
        stats["back_lengths_back"] = self.back_adapter.lengths_back
        # have to clone these nested dicts and set them
        # up with a custom merge function
        stats["front_errors_front"] = self.front_adapter.errors_front
        stats["front_errors_back"] = self.front_adapter.errors_back
        stats["back_errors_front"] = self.back_adapter.errors_front
        stats["back_errors_back"] = self.back_adapter.errors_back
        
        return stats

class AdapterCache(object):
    """Cache for known adapters.
    
    Args:
        path: Path to file where cache is written.
        auto_reverse_complement: Whether adapter reverse-complements should
            automatically be added to the cache.
    """
    def __init__(self, path=".adapters", auto_reverse_complement=False):
        self.path = path
        self.auto_reverse_complement = auto_reverse_complement
        if path and os.path.exists(path):
            with open(path, "rb") as cache:
                self.seq_to_name, self.name_to_seq = pickle.load(cache)
        else:
            self.seq_to_name = {}
            self.name_to_seq = {}
    
    @property
    def empty(self):
        """Whether the cache is empty.
        """
        return len(self.seq_to_name) == 0
    
    def save(self):
        """Save the cache to file.
        """
        if self.path is not None:
            with open(self.path, "wb") as cache:
                pickle.dump((self.seq_to_name, self.name_to_seq), cache)
    
    def add(self, name, seq):
        """Add a sequence to the cache.
        
        Args:
            name: Adapter name.
            seq: Adapter sequence.
        """
        self._add(name, seq)
        if self.auto_reverse_complement:
            self._add("{}_rc".format(name), reverse_complement(seq))
    
    def _add(self, name, seq):
        if seq not in self.seq_to_name:
            self.seq_to_name[seq] = set()
        self.seq_to_name[seq].add(name)
        self.name_to_seq[name] = seq
    
    def load_from_file(self, path=DEFAULT_ADAPTERS_PATH):
        """Load cached data from a file.
        
        Args:
            path: Path from which to load.
        """
        with open(path, "rt") as i:
            return self.load_from_fasta(i)

    def load_from_url(self, url=DEFAULT_ADAPTERS_URL):
        """Load adapter data from a URL.
        
        Args:
            url: URL from which to load.
        """
        logging.getLogger().info(
            "Loading list of known contaminants from %s", url)
        try:
            fasta = urlopen(url).read().decode().split("\n")
            return self.load_from_fasta(fasta)
        except URLError:
            if url.startswith("file:"):
                url = url[5:]
            return self.load_from_file(url)
    
    def load_from_fasta(self, fasta):
        """Load adapter data from a FASTA file.
        
        Args:
            fasta: FASTA file.
        """
        close = False
        if isinstance(fasta, str):
            fasta = open(fasta, 'rt')
            close = True
        num_records = None
        with FastaReader(fasta) as fasta:
            for num_records, record in enumerate(fasta, 1):
                name = record.name.split(None, 1)[0]
                seq = record.sequence
                self.add(name, seq)
        if close:
            fasta.close()
        return num_records
    
    def load_default(self):
        """Tries to load from default URL first, then from default path.
        """
        try:
            return self.load_from_url()
        except (OSError, IOError):
            logging.getLogger().warning(
                "Error loading adapters from URL %s; loading from file",
                DEFAULT_ADAPTERS_URL)
        try:
            return self.load_from_file()
        except IOError:
            logging.getLogger().warning(
                "Error loading adapters from file %s; loading from file",
                DEFAULT_ADAPTERS_PATH)
    
    @property
    def names(self):
        """Sequence of adapter names.
        """
        return list(self.name_to_seq.keys())
    
    @property
    def sequences(self):
        """Sequence of adapter sequences.
        """
        return list(self.seq_to_name.keys())
    
    def iter_names(self):
        """Returns an iterator over adapter names.
        """
        return self.name_to_seq.items()
    
    def iter_sequences(self):
        """Returns an iterator over adapter sequences.
        """
        return self.seq_to_name.items()
    
    def has_name(self, name):
        """Returns whether this cache contains the specified name.
        
        Args:
            name: The adapter name.
        """
        return name in self.name_to_seq
    
    def get_for_name(self, name):
        """Returns the sequence associated with a name.
        
        Args:
            name: The name to fetch.
        
        Returns:
            The sequence.
        """
        return self.name_to_seq[name]
    
    def has_seq(self, seq):
        """Tests whether a sequence is in the cache.
        
        Args:
            seq: The sequence to check.
        
        Returns:
            True if the sequence is in the cache.
        """
        return seq in self.seq_to_name
    
    def get_for_seq(self, seq):
        """Returns the name associated with a given sequence.
        
        Args:
            seq: The sequence to fetch.
        
        Returns:
            The name associated with the sequence.
        """
        return list(self.seq_to_name[seq])
    
    def summarize(self):
        """Returns a summary dict. Does *not* add sequence info.
        """
        return dict(
            path=self.path,
            auto_reverse_complement=self.auto_reverse_complement,
            num_adapter_names=len(self.name_to_seq),
            num_adapter_seqs=len(self.seq_to_name))

def parse_braces(sequence):
    """Replace all occurrences of ``x{n}`` (where x is any character) with n
    occurrences of x. Raise ValueError if the expression cannot be parsed.
    
    Examples:
        >>> parse_braces('TGA{5}CT')
        TGAAAAACT
    """
    # Simple DFA with four states, encoded in prev
    result = ''
    prev = None
    for char in re.split(r'(\{|\})', sequence):
        if char == '':
            continue
        if prev is None:
            if char == '{':
                raise ValueError('"{" must be used after a character')
            if char == '}':
                raise ValueError('"}" cannot be used here')
            prev = char
            result += char
        elif prev == '{':
            prev = int(char)
            if not 0 <= prev <= 10000:
                raise ValueError('Value {} invalid'.format(prev))
        elif isinstance(prev, int):
            if char != '}':
                raise ValueError('"}" expected')
            result = result[:-1] + result[-1] * prev
            prev = None
        else:
            if char != '{':
                raise ValueError('Expected "{"')
            prev = '{'
    # Check if we are in a non-terminating state
    if isinstance(prev, int) or prev == '{':
        raise ValueError("Unterminated expression")
    return result

def _extract_name_from_spec(spec):
    """Parse an adapter specification given as 'name=adapt' into 'name'
    and 'adapt'.
    
    Args:
        spec: Adapter spec.
    
    Returns:
        (name, spec)
    """
    fields = spec.split('=', 1)
    if len(fields) > 1:
        name, spec = fields
        name = name.strip()
    else:
        name = None
    spec = spec.strip()
    return name, spec

ADAPTER_ID_GENERATOR = itertools.count(1)

def _generate_adapter_name():
    return str(next(ADAPTER_ID_GENERATOR))

# coding: utf-8
"""
Adapters
"""
import itertools
import logging
import os
import pickle
import re
import sys
from urllib.error import URLError
from urllib.request import urlopen
from atropos import align
from atropos.align import Match
from atropos.io.seqio import ColorspaceSequence, FastaReader
from atropos.util import (
    IUPAC_BASES, GC_BASES, MergingDict, NestedDict, CountingDict, Const, 
    reverse_complement, ALPHABETS)
from atropos.util import colorspace as cs

class AdapterType(object):
    """Struct for adapter type information.
    
    Args:
        name: Adapter type name.
        desc: Adapter type description.
        flags: Alignment flags.
    """
    def __init__(self, name, desc, *flags):
        self.name = name
        self.desc = desc
        self.flags = flags[0]
        for i in range(1, len(flags)):
            self.flags |= flags[i]
    
    def asdict(self):
        """Returns AdapterType fields in a dict.
        """
        return dict(name=self.name, desc=self.desc, flags=Const(self.flags))

ADAPTER_TYPES = dict(
    back=AdapterType(
        'back', "regular 3'",
        align.START_WITHIN_SEQ2,
        align.STOP_WITHIN_SEQ2,
        align.STOP_WITHIN_SEQ1),
    front=AdapterType(
        'front', "regular 5'",
        align.START_WITHIN_SEQ2,
        align.STOP_WITHIN_SEQ2,
        align.START_WITHIN_SEQ1),
    prefix=AdapterType('prefix', "anchored 5'", align.STOP_WITHIN_SEQ2),
    suffix=AdapterType('suffix', "anchored 3'", align.START_WITHIN_SEQ2),
    anywhere=AdapterType('anywhere', "variable 5'/3'", align.SEMIGLOBAL),
    linked=AdapterType('linked', 'linked', 'linked')
)

def where_int_to_dict(where):
    """Convert a "where" flags integer to a dictionary of values from the
    corresponding AdapterType.
    """
    for adapter_type in ADAPTER_TYPES.values():
        if where == adapter_type.flags:
            return adapter_type.asdict()
    raise ValueError("Invalid WHERE value: {}".format(where))

# TODO get rid of these constants
BACK = ADAPTER_TYPES['back'].flags
FRONT = ADAPTER_TYPES['front'].flags
PREFIX = ADAPTER_TYPES['prefix'].flags
SUFFIX = ADAPTER_TYPES['suffix'].flags
ANYWHERE = ADAPTER_TYPES['anywhere'].flags
LINKED = ADAPTER_TYPES['linked'].flags

# TODO: specify this externally rather than hard-coding
DEFAULT_ADAPTERS_URL = "https://raw.githubusercontent.com/jdidion/atropos/master/atropos/adapters/sequencing_adapters.fa"
DEFAULT_ADAPTERS_PATH = os.path.join(
    os.path.dirname(__file__), 'sequencing_adapters.fa')

class AdapterParser(object):
    """Factory for Adapter classes that all use the same parameters (error rate,
    indels etc.). The given **kwargs will be passed to the Adapter constructors.
    
    Args:
        colorspace: Whether the adapter sequences are in colorspace.
        cache: The AdapterCache to use for fetching/storing named adapters.
    """
    def __init__(self, colorspace=False, cache=None, **kwargs):
        self.colorspace = colorspace
        self.cache = cache
        self.constructor_args = kwargs
        self.adapter_class = ColorspaceAdapter if colorspace else Adapter
    
    def parse(self, spec, cmdline_type='back'):
        """Parse an adapter specification not using ``file:`` notation and return
        an object of an appropriate Adapter class. The notation for anchored
        5' and 3' adapters is supported. If the name parameter is None, then
        an attempt is made to extract the name from the specification
        (If spec is 'name=ADAPTER', name will be 'name'.)
        
        Args:
            spec: The adapter spec.
            name: The adapter name. If not provided, one is automatically
                generated.
            cmdline_type: describes which commandline parameter was used (``-a``
                is 'back', ``-b`` is 'anywhere', and ``-g`` is 'front').
        
        TODO: describe the adapter spec format
        
        Returns:
            An :class:`Adapter` instance.
        """
        if spec.startswith('file:'):
            # read adapter sequences from a file
            with FastaReader(spec[5:]) as fasta:
                for record in fasta:
                    name = record.name.split(None, 1)[0]
                    yield self.parse_from_spec(
                        record.sequence, cmdline_type, name)
        else:
            yield self.parse_from_spec(spec, cmdline_type)
    
    def parse_from_spec(self, spec, cmdline_type='back', name=None):
        if cmdline_type not in ADAPTER_TYPES:
            raise ValueError(
                'cmdline_type cannot be {0!r}'.format(cmdline_type))
        orig_spec = spec
        where = ADAPTER_TYPES[cmdline_type].flags
        
        if name is None and spec is None:
            raise ValueError('Either name or spec must be given')
        elif name is None:
            if self.cache and self.cache.has_name(spec):
                name = spec
                spec = self.cache.get_for_name(name)
        elif spec is None:
            if self.cache and self.cache.has_name(name):
                spec = self.cache.get_for_name(name)
        
        if spec is None:
            raise ValueError('Name not found: {}'.format(name))
        elif name is None:
            name, spec = _extract_name_from_spec(spec)
        
        if self.cache and name is not None:
            self.cache.add(name, spec)
        
        front_anchored, back_anchored = False, False
        if spec.startswith('^'):
            spec = spec[1:]
            front_anchored = True
        if spec.endswith('$'):
            spec = spec[:-1]
            back_anchored = True
        
        sequence1, middle, sequence2 = spec.partition('...')
        
        if where == ANYWHERE:
            if front_anchored or back_anchored:
                raise ValueError("'anywhere' (-b) adapters may not be anchored")
            if middle == '...':
                raise ValueError("'anywhere' (-b) adapters may not be linked")
            return self.adapter_class(
                sequence=spec, where=where, name=name, **self.constructor_args)
        
        assert where == FRONT or where == BACK
        if middle == '...':
            if not sequence1:
                if where == BACK:  # -a ...ADAPTER
                    spec = sequence2
                else:  # -g ...ADAPTER
                    raise ValueError('Invalid adapter specification')
            elif not sequence2:
                if where == BACK:  # -a ADAPTER...
                    spec = sequence1
                    where = FRONT
                    front_anchored = True
                else:  # -g ADAPTER...
                    spec = sequence1
            else:
                # linked adapter
                if self.colorspace:
                    raise NotImplementedError(
                        'Using linked adapters in colorspace is not supported')
                # automatically anchor 5' adapter if -a is used
                if where == BACK:
                    front_anchored = True
                
                return LinkedAdapter(sequence1, sequence2, name=name,
                    front_anchored=front_anchored, back_anchored=back_anchored,
                    **self.constructor_args)
        
        if front_anchored and back_anchored:
            raise ValueError(
                'Trying to use both "^" and "$" in adapter specification '
                '{!r}'.format(orig_spec))
        if front_anchored:
            if where == BACK:
                raise ValueError("Cannot anchor the 3' adapter at its 5' end")
            where = PREFIX
        elif back_anchored:
            if where == FRONT:
                raise ValueError("Cannot anchor 5' adapter at 3' end")
            where = SUFFIX
        
        return self.adapter_class(
            sequence=spec, where=where, name=name, **self.constructor_args)
    
    def parse_multi(self, back=None, anywhere=None, front=None):
        """Parse all three types of commandline options that can be used to
        specify adapters. back, anywhere and front are lists of strings,
        corresponding to the respective commandline types (-a, -b, -g).
        
        Args:
            back: Back-adapter specs.
            anywhere: Anywhere-adapter specs.
            front: Front-adapter specs.
        
        Return:
            A list of appropriate Adapter classes.
        """
        adapters = []
        for specs, cmdline_type in (
                (back, 'back'), (anywhere, 'anywhere'), (front, 'front')):
            if not specs:
                continue
            for spec in specs:
                adapters.extend(self.parse(spec, cmdline_type))
        return adapters

class Adapter(object):
    """An adapter knows how to match itself to a read. In particular, it knows
    where it should be within the read and how to interpret wildcard characters.
    
    Args:
        sequence: The adapter sequence as string. Will be converted to
            uppercase. Also, Us will be converted to Ts.
        where: One of the BACK, FRONT, PREFIX, SUFFIX or ANYWHERE constants.
            This influences where the adapter is allowed to appear within in the
            read and also which part of the read is removed.
        max_error_rate: Maximum allowed error rate. The error rate is the number
            of errors in the alignment divided by the length of the part of the
            alignment that matches the adapter.
        min_overlap: Minimum length of the part of the alignment that
            matches the adapter.
        read_wildcards: Whether IUPAC wildcards in the read are allowed.
        adapter_wildcards: Whether IUPAC wildcards in the adapter are allowed.
        name: optional name of the adapter. If not provided, the name is set to
            a unique number.
        indels: Whether indels are allowed in adapter matches.
        indel_cost: Cost of indel bases in alignment scoring.
        match_probability: Callabale that computes the probability that a match
            could happen by random chance.
        max_rmp: Maximum random-match probability below which a match is
            considered genuine.
        gc_content: Expected GC content of sequences.
        alphabet: The Alphabet to use to validate the adapter sequence.
    """
    def __init__(
            self, sequence, where, max_error_rate=0.1, min_overlap=3,
            read_wildcards=False, adapter_wildcards=True, name=None,
            indels=True, indel_cost=1, match_probability=None, max_rmp=None,
            gc_content=0.5, alphabet=None):
        if len(sequence) == 0:
            raise ValueError("Empty adapter sequence")
        # TODO: all of this validation code should be wrapped up in Alphabet
        sequence = parse_braces(sequence.upper().replace('U', 'T'))
        seq_set = set(sequence)
        if seq_set <= set('ACGT'):
            adapter_wildcards = False
        if adapter_wildcards and not seq_set <= IUPAC_BASES:
            raise ValueError(
                "Invalid character(s) in adapter sequence: {}".format(
                    ','.join(seq_set - IUPAC_BASES)))
        if alphabet:
            if isinstance(alphabet, str):
                alphabet = ALPHABETS[alphabet]
            alphabet.validate_string(sequence)
        
        self.debug = False
        self.name = _generate_adapter_name() if name is None else name
        self.sequence = sequence
        self.where = where
        self.max_error_rate = max_error_rate
        self.min_overlap = min(min_overlap, len(self.sequence))
        self.match_probability = match_probability
        self.max_rmp = max_rmp
        self.gc_content = gc_content
        self.indels = indels
        self.adapter_wildcards = adapter_wildcards
        self.read_wildcards = read_wildcards
        # redirect trimmed() to appropriate function depending on adapter type
        trimmers = {
            FRONT: self._trimmed_front,
            PREFIX: self._trimmed_front,
            BACK: self._trimmed_back,
            SUFFIX: self._trimmed_back,
            ANYWHERE: self._trimmed_anywhere
        }
        self.trimmed = trimmers[where]
        if where == ANYWHERE:
            self._front_flag = None  # means: guess
        else:
            self._front_flag = where not in (BACK, SUFFIX)
        # statistics about length of removed sequences
        self.lengths_front = CountingDict()
        self.lengths_back = CountingDict()
        self.errors_front = NestedDict()
        self.errors_back = NestedDict()
        self.adjacent_bases = { 'A': 0, 'C': 0, 'G': 0, 'T': 0, '': 0 }
        self.aligner = align.Aligner(
            self.sequence, self.max_error_rate, flags=self.where,
            wildcard_ref=self.adapter_wildcards,
            wildcard_query=self.read_wildcards)
        self.aligner.min_overlap = self.min_overlap
        if self.indels:
            self.aligner.indel_cost = indel_cost
        else:
            # TODO
            # When indels are disallowed, an entirely different algorithm
            # should be used.
            self.aligner.indel_cost = 100000
    
    def __repr__(self):
        return '<Adapter(name="{name}", sequence="{sequence}", where={where}, '\
            'max_error_rate={max_error_rate}, min_overlap={min_overlap}, '\
            'read_wildcards={read_wildcards}, '\
            'adapter_wildcards={adapter_wildcards}, '\
            'indels={indels})>'.format(**vars(self))
    
    def enable_debug(self):
        """Print out the dynamic programming matrix after matching a read to an
        adapter.
        """
        self.debug = True
        self.aligner.enable_debug()
    
    def match_to(self, read):
        """Attempt to match this adapter to the given read.
        
        Args:
            read: A :class:`Sequence` instance.
        
        Returns:
            A :class:`Match` instance if a match was found; return None if no
            match was found given the matching criteria (minimum overlap length,
            maximum error rate).
        """
        read_seq = read.sequence.upper()
        
        # try to find an exact match first unless wildcards are allowed
        pos = -1
        if not self.adapter_wildcards:
            if self.where == PREFIX:
                if read_seq.startswith(self.sequence):
                    pos = 0
            elif self.where == SUFFIX:
                if read_seq.endswith(self.sequence):
                    pos = (len(read_seq) - len(self.sequence))
            else:
                pos = read_seq.find(self.sequence)
        
        if pos >= 0:
            seqlen = len(self.sequence)
            return Match(
                0, seqlen, pos, pos + seqlen, seqlen, 0, self._front_flag,
                self, read)
        
        # try approximate matching
        if not self.indels and self.where in (PREFIX, SUFFIX):
            if self.where == PREFIX:
                alignment = align.compare_prefixes(
                    self.sequence, read_seq,
                    wildcard_ref=self.adapter_wildcards,
                    wildcard_query=self.read_wildcards)
            else:
                alignment = align.compare_suffixes(
                    self.sequence, read_seq,
                    wildcard_ref=self.adapter_wildcards,
                    wildcard_query=self.read_wildcards)
        else:
            alignment = self.aligner.locate(read_seq)
            if self.debug:
                print(self.aligner.dpmatrix)  # pragma: no cover
        
        if alignment:
            astart, astop, rstart, rstop, matches, errors = alignment
            size = astop - astart
            if ((
                    size >=
                    self.min_overlap and errors / size <=
                    self.max_error_rate
                ) and (
                    self.max_rmp is None or
                    self.match_probability(matches, size) <= self.max_rmp)):
                return Match(
                    astart, astop, rstart, rstop, matches, errors,
                    self._front_flag, self, read)
        
        return None
    
    def _trimmed_anywhere(self, match):
        """Trims an adapter from either the front or back of sequence.
        
        Returns:
            A :class:`Sequence` instance: the trimmed read.
        """
        if match.front:
            return self._trimmed_front(match)
        else:
            return self._trimmed_back(match)
    
    def _trimmed_front(self, match):
        """Trims an adapter from the front of sequence.
        
        Returns:
            A :class:`Sequence` instance: the trimmed read.
        """
        # TODO move away
        self.lengths_front[match.rstop] += 1
        self.errors_front[match.rstop][match.errors] += 1
        return match.read[match.rstop:]
    
    def _trimmed_back(self, match):
        """Trims an adapter from the back of sequence.
        
        Returns:
            A :class:`Sequence` instance: the trimmed read.
        """
        self.lengths_back[len(match.read) - match.rstart] += 1
        self.errors_back[len(match.read) - match.rstart][match.errors] += 1
        adjacent_base = match.read.sequence[match.rstart-1:match.rstart]
        if adjacent_base not in 'ACGT':
            adjacent_base = ''
        self.adjacent_bases[adjacent_base] += 1
        return match.read[:match.rstart]
    
    def __len__(self):
        return len(self.sequence)
    
    def random_match_probabilities(self):
        """Estimate probabilities that this adapter matches a random sequence. 
        Indels are not taken into account.

        Returns:
            A list of probabilities the same length as this adapter's sequence,
            where the value at position 'i' is the probability that i bases of
            this adapter match a random sequence.
        """
        if self._front_flag:
            seq = self.sequence[::-1]
        else:
            seq = self.sequence
        
        #matches = 0
        base_probs = (self.gc_content / 2.0, (1 - self.gc_content) / 2.0)
        probabilities = [1.0] + ([0] * len(seq))
        c_bases = frozenset(GC_BASES if self.adapter_wildcards else 'GC')
        
        # TODO: this doesn't work; need to figure out if RandomMatchProbability
        # can be used for this.
        #for idx, base in enumerate(seq, 1):
        #    if base in c_bases:
        #        matches += 1
        #    probabilities[idx] = self.match_probability(
        #        matches, idx, *base_probs)
        
        cur_p = 1.0
        for idx, base in enumerate(seq, 1):
            cur_p *= base_probs[0 if base in c_bases else 1]
            probabilities[idx] = cur_p
        return probabilities
    
    def summarize(self):
        """Summarize the activities of this :class:`Adapter`.
        """
        total_front = sum(self.lengths_front.values())
        total_back = sum(self.lengths_back.values())
        
        stats = MergingDict(
            adapter_class=self.__class__.__name__,
            total_front=total_front,
            total_back=total_back,
            total=total_front + total_back,
            match_probabilities=Const(self.random_match_probabilities()))
        
        where = self.where
        assert (
            where in (ANYWHERE, LINKED) or
            (where in (BACK, SUFFIX) and total_front == 0) or
            (where in (FRONT, PREFIX) and total_back == 0))
        
        stats["where"] = where_int_to_dict(where)
        stats["sequence"] = Const(self.sequence)
        stats["max_error_rate"] = Const(self.max_error_rate)
        if where in (ANYWHERE, FRONT, PREFIX):
            stats["lengths_front"] = self.lengths_front
            stats["errors_front"] = self.errors_front
        if where in (ANYWHERE, BACK, SUFFIX):
            stats["lengths_back"] = self.lengths_back
            stats["errors_back"] = self.errors_back
        if where in (BACK, SUFFIX):
            stats["adjacent_bases"] = self.adjacent_bases
        
        return stats

class ColorspaceAdapter(Adapter):
    """An adapter for a colorspace sequence.
    
    Args:
        args, kwargs: Arguments to pass to :class:`Adapter` constructor.
    """
    def __init__(self, *args, **kwargs):
        if kwargs.get('adapter_wildcards', False):
            raise ValueError("Wildcards not supported for colorspace adapters")
        kwargs['adapter_wildcards'] = False
        super().__init__(*args, **kwargs)
        has_nucleotide_seq = False
        if set(self.sequence) <= set('ACGT'):
            # adapter was given in basespace
            self.nucleotide_sequence = self.sequence
            has_nucleotide_seq = True
            self.sequence = cs.encode(self.sequence)[1:]
        if self.where in (PREFIX, FRONT) and not has_nucleotide_seq:
            raise ValueError(
                "A 5' colorspace adapter needs to be given in nucleotide space")
        self.aligner.reference = self.sequence
    
    def match_to(self, read):
        """Attempt to match this adapter to the given read.
        
        Args:
            read: A :class:`Sequence` instance.
        
        Returns:
            A :class:`Match` instance if a match was found; return None if no
            match was found given the matching criteria (minimum overlap length,
            maximum error rate).
        """
        if self.where != PREFIX:
            return super().match_to(read)
        # create artificial adapter that includes a first color that encodes the
        # transition from primer base into adapter
        asequence = (
            cs.ENCODE[read.primer + self.nucleotide_sequence[0:1]] +
            self.sequence)
        
        pos = 0 if read.sequence.startswith(asequence) else -1
        if pos >= 0:
            match = Match(
                0, len(asequence), pos, pos + len(asequence),
                len(asequence), 0, self._front_flag, self, read)
        else:
            # try approximate matching
            self.aligner.reference = asequence
            alignment = self.aligner.locate(read.sequence)
            if self.debug:
                print(self.aligner.dpmatrix)  # pragma: no cover
            if alignment is not None:
                match = Match(*(alignment + (self._front_flag, self, read)))
            else:
                match = None
        
        if match is None:
            return None
        assert (
            match.length > 0 and
            match.errors / match.length <= self.max_error_rate)
        assert match.length >= self.min_overlap
        return match

    def _trimmed_front(self, match):
        """Trims an adapter from the front of sequence.
        
        Returns:
            A :class:`Sequence` instance: the trimmed read.
        """
        read = match.read
        self.lengths_front[match.rstop] += 1
        self.errors_front[match.rstop][match.errors] += 1
        # to remove a front adapter, we need to re-encode the first color
        # following the adapter match
        color_after_adapter = read.sequence[match.rstop:match.rstop + 1]
        if not color_after_adapter:
            # the read is empty
            return read[match.rstop:]
        base_after_adapter = cs.DECODE[
            self.nucleotide_sequence[-1:] + color_after_adapter]
        new_first_color = cs.ENCODE[read.primer + base_after_adapter]
        new_read = read[:]
        new_read.sequence = new_first_color + read.sequence[(match.rstop + 1):]
        new_read.qualities = None
        if read.qualities:
            new_read.qualities = read.qualities[match.rstop:]
        return new_read

    def _trimmed_back(self, match):
        """Trims an adapter from the back of sequence.
        
        Returns:
            A :class:`Sequence` instance: the trimmed read.
        """
        # trim one more color if long enough
        adjusted_rstart = max(match.rstart - 1, 0)
        self.lengths_back[len(match.read) - adjusted_rstart] += 1
        self.errors_back[len(match.read) - adjusted_rstart][match.errors] += 1
        return match.read[:adjusted_rstart]

    def __repr__(self):
        return '<ColorspaceAdapter(sequence={0!r}, where={1})>'.format(
            self.sequence, self.where)

# TODO Consolidate Match and a LinkedMatch class.

class LinkedMatch(object):
    """Represent a match of a LinkedAdapter.
    
    Args:
        front_match: The match to the front of the sequence.
        back_match: The match to the back of the sequence.
        adapter: The matched adapter.
    """
    def __init__(self, front_match, back_match, adapter):
        self.front_match = front_match
        self.back_match = back_match
        self.adapter = adapter
        assert front_match is not None
    
    def get_info_record(self):
        """Returns the info record for the either the back or forward match.
        """
        if self.back_match:
            return self.back_match.get_info_record()
        else:
            return self.front_match.get_info_record()

class LinkedAdapter(object):
    """An adapter with linked front and back sequences.
    
    Args:
        front_sequence: Front adapter sequence.
        back_sequence: Back adapter sequence.
        front_anchored: Whether the front adapter is anchored.
        back_anchored: Whether the back adapter is anchored.
        name: Adapter name.
        kwargs: passed on to individual :class:`Adapter` constructors.
    """
    def __init__(
            self, front_sequence, back_sequence, front_anchored=True,
            back_anchored=False, name=None, **kwargs):
        assert front_anchored and not back_anchored
        where1 = PREFIX if front_anchored else FRONT
        where2 = SUFFIX if back_anchored else BACK
        self.front_anchored = front_anchored
        self.back_anchored = back_anchored
        
        # The following attributes are needed for the report
        self.where = LINKED
        self.name = _generate_adapter_name() if name is None else name
        self.front_adapter = Adapter(
            front_sequence, where=where1, name=None, **kwargs)
        self.back_adapter = Adapter(
            back_sequence, where=where2, name=None, **kwargs)
    
    def enable_debug(self):
        """Enable debug on adapters.
        """
        self.front_adapter.enable_debug()
        self.back_adapter.enable_debug()
    
    def match_to(self, read):
        """Match the linked adapters against the given read. If the 'front'
        adapter is not found, the 'back' adapter is not searched for.
        
        Args:
            read: A :class:`Sequence` instance.
        
        Returns:
            A :class:`Match` instance if a match was found; return None if no
            match was found given the matching criteria (minimum overlap length,
            maximum error rate).
        """
        front_match = self.front_adapter.match_to(read)
        if front_match is None:
            return None
        # TODO use match.trimmed() instead as soon as that does not update
        # statistics anymore
        read = read[front_match.rstop:]
        back_match = self.back_adapter.match_to(read)
        return LinkedMatch(front_match, back_match, self)

    def trimmed(self, match):
        """Returns the read trimmed with the front and/or back adapter
        trimmer(s).
        
        Args:
            match: The match to trim.
        
        Returns:
            The trimmed read.
        """
        front_trimmed = self.front_adapter.trimmed(match.front_match)
        if match.back_match:
            return self.back_adapter.trimmed(match.back_match)
        else:
            return front_trimmed
    
    def summarize(self):
        """Returns the summary dict for this adapter.
        """
        total_front = sum(self.front_adapter.lengths_front.values())
        total_back = sum(self.back_adapter.lengths_back.values())
        
        stats = MergingDict(
            total_front=total_front,
            total_back=total_back,
            total=total_front + total_back)
        
        where = self.where
        assert (
            where in (ANYWHERE, LINKED) or
            (where in (BACK, SUFFIX) and total_front == 0) or
            (where in (FRONT, PREFIX) and total_back == 0))
        
        stats["where"] = where_int_to_dict(where)
        stats["front_sequence"] = Const(self.front_adapter.sequence)
        stats["front_match_probabilities"] = Const(
            self.front_adapter.random_match_probabilities())
        stats["back_sequence"] = Const(self.back_adapter.sequence)
        stats["back_match_probabilities"] = Const(
            self.back_adapter.random_match_probabilities())
        stats["front_max_error_rate"] = Const(self.front_adapter.max_error_rate)
        stats["back_max_error_rate"] = Const(self.back_adapter.max_error_rate)
        stats["front_lengths_front"] = self.front_adapter.lengths_front
        stats["front_lengths_back"] = self.front_adapter.lengths_back
        stats["back_lengths_front"] = self.back_adapter.lengths_front
        stats["back_lengths_back"] = self.back_adapter.lengths_back
        # have to clone these nested dicts and set them
        # up with a custom merge function
        stats["front_errors_front"] = self.front_adapter.errors_front
        stats["front_errors_back"] = self.front_adapter.errors_back
        stats["back_errors_front"] = self.back_adapter.errors_front
        stats["back_errors_back"] = self.back_adapter.errors_back
        
        return stats

class AdapterCache(object):
    """Cache for known adapters.
    
    Args:
        path: Path to file where cache is written.
        auto_reverse_complement: Whether adapter reverse-complements should
            automatically be added to the cache.
    """
    def __init__(self, path=".adapters", auto_reverse_complement=False):
        self.path = path
        self.auto_reverse_complement = auto_reverse_complement
        if path and os.path.exists(path):
            with open(path, "rb") as cache:
                self.seq_to_name, self.name_to_seq = pickle.load(cache)
        else:
            self.seq_to_name = {}
            self.name_to_seq = {}
    
    @property
    def empty(self):
        """Whether the cache is empty.
        """
        return len(self.seq_to_name) == 0
    
    def save(self):
        """Save the cache to file.
        """
        if self.path is not None:
            with open(self.path, "wb") as cache:
                pickle.dump((self.seq_to_name, self.name_to_seq), cache)
    
    def add(self, name, seq):
        """Add a sequence to the cache.
        
        Args:
            name: Adapter name.
            seq: Adapter sequence.
        """
        self._add(name, seq)
        if self.auto_reverse_complement:
            self._add("{}_rc".format(name), reverse_complement(seq))
    
    def _add(self, name, seq):
        if seq not in self.seq_to_name:
            self.seq_to_name[seq] = set()
        self.seq_to_name[seq].add(name)
        self.name_to_seq[name] = seq
    
    def load_from_file(self, path=DEFAULT_ADAPTERS_PATH):
        """Load cached data from a file.
        
        Args:
            path: Path from which to load.
        """
        with open(path, "rt") as i:
            return self.load_from_fasta(i)

    def load_from_url(self, url=DEFAULT_ADAPTERS_URL):
        """Load adapter data from a URL.
        
        Args:
            url: URL from which to load.
        """
        logging.getLogger().info(
            "Loading list of known contaminants from %s", url)
        try:
            fasta = urlopen(url).read().decode().split("\n")
            return self.load_from_fasta(fasta)
        except URLError:
            if url.startswith("file:"):
                url = url[5:]
            return self.load_from_file(url)
    
    def load_from_fasta(self, fasta):
        """Load adapter data from a FASTA file.
        
        Args:
            fasta: FASTA file.
        """
        close = False
        if isinstance(fasta, str):
            fasta = open(fasta, 'rt')
            close = True
        num_records = None
        with FastaReader(fasta) as fasta:
            for num_records, record in enumerate(fasta, 1):
                name = record.name.split(None, 1)[0]
                seq = record.sequence
                self.add(name, seq)
        if close:
            fasta.close()
        return num_records
    
    def load_default(self):
        """Tries to load from default URL first, then from default path.
        """
        try:
            return self.load_from_url()
        except (OSError, IOError):
            logging.getLogger().warning(
                "Error loading adapters from URL %s; loading from file",
                DEFAULT_ADAPTERS_URL)
        try:
            return self.load_from_file()
        except IOError:
            logging.getLogger().warning(
                "Error loading adapters from file %s; loading from file",
                DEFAULT_ADAPTERS_PATH)
    
    @property
    def names(self):
        """Sequence of adapter names.
        """
        return list(self.name_to_seq.keys())
    
    @property
    def sequences(self):
        """Sequence of adapter sequences.
        """
        return list(self.seq_to_name.keys())
    
    def iter_names(self):
        """Returns an iterator over adapter names.
        """
        return self.name_to_seq.items()
    
    def iter_sequences(self):
        """Returns an iterator over adapter sequences.
        """
        return self.seq_to_name.items()
    
    def has_name(self, name):
        """Returns whether this cache contains the specified name.
        
        Args:
            name: The adapter name.
        """
        return name in self.name_to_seq
    
    def get_for_name(self, name):
        """Returns the sequence associated with a name.
        
        Args:
            name: The name to fetch.
        
        Returns:
            The sequence.
        """
        return self.name_to_seq[name]
    
    def has_seq(self, seq):
        """Tests whether a sequence is in the cache.
        
        Args:
            seq: The sequence to check.
        
        Returns:
            True if the sequence is in the cache.
        """
        return seq in self.seq_to_name
    
    def get_for_seq(self, seq):
        """Returns the name associated with a given sequence.
        
        Args:
            seq: The sequence to fetch.
        
        Returns:
            The name associated with the sequence.
        """
        return list(self.seq_to_name[seq])
    
    def summarize(self):
        """Returns a summary dict. Does *not* add sequence info.
        """
        return dict(
            path=self.path,
            auto_reverse_complement=self.auto_reverse_complement,
            num_adapter_names=len(self.name_to_seq),
            num_adapter_seqs=len(self.seq_to_name))

def parse_braces(sequence):
    """Replace all occurrences of ``x{n}`` (where x is any character) with n
    occurrences of x. Raise ValueError if the expression cannot be parsed.
    
    Examples:
        >>> parse_braces('TGA{5}CT')
        TGAAAAACT
    """
    # Simple DFA with four states, encoded in prev
    result = ''
    prev = None
    for char in re.split(r'(\{|\})', sequence):
        if char == '':
            continue
        if prev is None:
            if char == '{':
                raise ValueError('"{" must be used after a character')
            if char == '}':
                raise ValueError('"}" cannot be used here')
            prev = char
            result += char
        elif prev == '{':
            prev = int(char)
            if not 0 <= prev <= 10000:
                raise ValueError('Value {} invalid'.format(prev))
        elif isinstance(prev, int):
            if char != '}':
                raise ValueError('"}" expected')
            result = result[:-1] + result[-1] * prev
            prev = None
        else:
            if char != '{':
                raise ValueError('Expected "{"')
            prev = '{'
    # Check if we are in a non-terminating state
    if isinstance(prev, int) or prev == '{':
        raise ValueError("Unterminated expression")
    return result

def _extract_name_from_spec(spec):
    """Parse an adapter specification given as 'name=adapt' into 'name'
    and 'adapt'.
    
    Args:
        spec: Adapter spec.
    
    Returns:
        (name, spec)
    """
    fields = spec.split('=', 1)
    if len(fields) > 1:
        name, spec = fields
        name = name.strip()
    else:
        name = None
    spec = spec.strip()
    return name, spec

ADAPTER_ID_GENERATOR = itertools.count(1)

def _generate_adapter_name():
    return str(next(ADAPTER_ID_GENERATOR))
# coding: utf-8 Adapters Struct for adapter type information.          Args:         name: Adapter type name.         desc: Adapter type description.         flags: Alignment flags. Returns AdapterType fields in a dict. Convert a "where" flags integer to a dictionary of values from the     corresponding AdapterType. # TODO get rid of these constants # TODO: specify this externally rather than hard-coding Factory for Adapter classes that all use the same parameters (error rate,     indels etc.). The given **kwargs will be passed to the Adapter constructors.          Args:         colorspace: Whether the adapter sequences are in colorspace.         cache: The AdapterCache to use for fetching/storing named adapters. Parse an adapter specification not using ``file:`` notation and return         an object of an appropriate Adapter class. The notation for anchored         5' and 3' adapters is supported. If the name parameter is None, then         an attempt is made to extract the name from the specification         (If spec is 'name=ADAPTER', name will be 'name'.)                  Args:             spec: The adapter spec.             name: The adapter name. If not provided, one is automatically                 generated.             cmdline_type: describes which commandline parameter was used (``-a``                 is 'back', ``-b`` is 'anywhere', and ``-g`` is 'front').                  TODO: describe the adapter spec format                  Returns:             An :class:`Adapter` instance. # read adapter sequences from a file # -a ...ADAPTER # -g ...ADAPTER # -a ADAPTER... # -g ADAPTER... # linked adapter # automatically anchor 5' adapter if -a is used Parse all three types of commandline options that can be used to         specify adapters. back, anywhere and front are lists of strings,         corresponding to the respective commandline types (-a, -b, -g).                  Args:             back: Back-adapter specs.             anywhere: Anywhere-adapter specs.             front: Front-adapter specs.                  Return:             A list of appropriate Adapter classes. An adapter knows how to match itself to a read. In particular, it knows     where it should be within the read and how to interpret wildcard characters.          Args:         sequence: The adapter sequence as string. Will be converted to             uppercase. Also, Us will be converted to Ts.         where: One of the BACK, FRONT, PREFIX, SUFFIX or ANYWHERE constants.             This influences where the adapter is allowed to appear within in the             read and also which part of the read is removed.         max_error_rate: Maximum allowed error rate. The error rate is the number             of errors in the alignment divided by the length of the part of the             alignment that matches the adapter.         min_overlap: Minimum length of the part of the alignment that             matches the adapter.         read_wildcards: Whether IUPAC wildcards in the read are allowed.         adapter_wildcards: Whether IUPAC wildcards in the adapter are allowed.         name: optional name of the adapter. If not provided, the name is set to             a unique number.         indels: Whether indels are allowed in adapter matches.         indel_cost: Cost of indel bases in alignment scoring.         match_probability: Callabale that computes the probability that a match             could happen by random chance.         max_rmp: Maximum random-match probability below which a match is             considered genuine.         gc_content: Expected GC content of sequences.         alphabet: The Alphabet to use to validate the adapter sequence. # TODO: all of this validation code should be wrapped up in Alphabet # redirect trimmed() to appropriate function depending on adapter type # means: guess # statistics about length of removed sequences # TODO # When indels are disallowed, an entirely different algorithm # should be used. Print out the dynamic programming matrix after matching a read to an         adapter. Attempt to match this adapter to the given read.                  Args:             read: A :class:`Sequence` instance.                  Returns:             A :class:`Match` instance if a match was found; return None if no             match was found given the matching criteria (minimum overlap length,             maximum error rate). # try to find an exact match first unless wildcards are allowed # try approximate matching # pragma: no cover Trims an adapter from either the front or back of sequence.                  Returns:             A :class:`Sequence` instance: the trimmed read. Trims an adapter from the front of sequence.                  Returns:             A :class:`Sequence` instance: the trimmed read. # TODO move away Trims an adapter from the back of sequence.                  Returns:             A :class:`Sequence` instance: the trimmed read. Estimate probabilities that this adapter matches a random sequence.          Indels are not taken into account.          Returns:             A list of probabilities the same length as this adapter's sequence,             where the value at position 'i' is the probability that i bases of             this adapter match a random sequence. #matches = 0 # TODO: this doesn't work; need to figure out if RandomMatchProbability # can be used for this. #for idx, base in enumerate(seq, 1): #    if base in c_bases: #        matches += 1 #    probabilities[idx] = self.match_probability( #        matches, idx, *base_probs) Summarize the activities of this :class:`Adapter`. An adapter for a colorspace sequence.          Args:         args, kwargs: Arguments to pass to :class:`Adapter` constructor. # adapter was given in basespace Attempt to match this adapter to the given read.                  Args:             read: A :class:`Sequence` instance.                  Returns:             A :class:`Match` instance if a match was found; return None if no             match was found given the matching criteria (minimum overlap length,             maximum error rate). # create artificial adapter that includes a first color that encodes the # transition from primer base into adapter # try approximate matching # pragma: no cover Trims an adapter from the front of sequence.                  Returns:             A :class:`Sequence` instance: the trimmed read. # to remove a front adapter, we need to re-encode the first color # following the adapter match # the read is empty Trims an adapter from the back of sequence.                  Returns:             A :class:`Sequence` instance: the trimmed read. # trim one more color if long enough # TODO Consolidate Match and a LinkedMatch class. Represent a match of a LinkedAdapter.          Args:         front_match: The match to the front of the sequence.         back_match: The match to the back of the sequence.         adapter: The matched adapter. Returns the info record for the either the back or forward match. An adapter with linked front and back sequences.          Args:         front_sequence: Front adapter sequence.         back_sequence: Back adapter sequence.         front_anchored: Whether the front adapter is anchored.         back_anchored: Whether the back adapter is anchored.         name: Adapter name.         kwargs: passed on to individual :class:`Adapter` constructors. # The following attributes are needed for the report Enable debug on adapters. Match the linked adapters against the given read. If the 'front'         adapter is not found, the 'back' adapter is not searched for.                  Args:             read: A :class:`Sequence` instance.                  Returns:             A :class:`Match` instance if a match was found; return None if no             match was found given the matching criteria (minimum overlap length,             maximum error rate). # TODO use match.trimmed() instead as soon as that does not update # statistics anymore Returns the read trimmed with the front and/or back adapter         trimmer(s).                  Args:             match: The match to trim.                  Returns:             The trimmed read. Returns the summary dict for this adapter. # have to clone these nested dicts and set them # up with a custom merge function Cache for known adapters.          Args:         path: Path to file where cache is written.         auto_reverse_complement: Whether adapter reverse-complements should             automatically be added to the cache. Whether the cache is empty. Save the cache to file. Add a sequence to the cache.                  Args:             name: Adapter name.             seq: Adapter sequence. Load cached data from a file.                  Args:             path: Path from which to load. Load adapter data from a URL.                  Args:             url: URL from which to load. Load adapter data from a FASTA file.                  Args:             fasta: FASTA file. Tries to load from default URL first, then from default path. Sequence of adapter names. Sequence of adapter sequences. Returns an iterator over adapter names. Returns an iterator over adapter sequences. Returns whether this cache contains the specified name.                  Args:             name: The adapter name. Returns the sequence associated with a name.                  Args:             name: The name to fetch.                  Returns:             The sequence. Tests whether a sequence is in the cache.                  Args:             seq: The sequence to check.                  Returns:             True if the sequence is in the cache. Returns the name associated with a given sequence.                  Args:             seq: The sequence to fetch.                  Returns:             The name associated with the sequence. Returns a summary dict. Does *not* add sequence info. Replace all occurrences of ``x{n}`` (where x is any character) with n     occurrences of x. Raise ValueError if the expression cannot be parsed.          Examples:         >>> parse_braces('TGA{5}CT')         TGAAAAACT # Simple DFA with four states, encoded in prev # Check if we are in a non-terminating state Parse an adapter specification given as 'name=adapt' into 'name'     and 'adapt'.          Args:         spec: Adapter spec.          Returns:         (name, spec)