# ---------------------------------------------------------------------------- # Copyright (c) 2013--, scikit-bio development team. # # Distributed under the terms of the Modified BSD License. # # The full license is in the file COPYING.txt, distributed with this software. # ---------------------------------------------------------------------------- import re import collections import numbers from contextlib import contextmanager import numpy as np import pandas as pd import skbio.sequence.distance from skbio._base import SkbioObject from skbio.metadata._mixin import (MetadataMixin, PositionalMetadataMixin, IntervalMetadataMixin) from skbio.metadata import IntervalMetadata from skbio.sequence._repr import _SequenceReprBuilder from skbio.util._decorator import (stable, experimental, classonlymethod, overrides) class Sequence(MetadataMixin, PositionalMetadataMixin, IntervalMetadataMixin, collections.abc.Sequence, SkbioObject): """Store generic sequence data and optional associated metadata. ``Sequence`` objects do not enforce an alphabet or grammar and are thus the most generic objects for storing sequence data. ``Sequence`` objects do not necessarily represent biological sequences. For example, ``Sequence`` can be used to represent a position in a multiple sequence alignment. Subclasses ``DNA``, ``RNA``, and ``Protein`` enforce the IUPAC character set [1]_ for, and provide operations specific to, each respective molecule type. ``Sequence`` objects consist of the underlying sequence data, as well as optional metadata and positional metadata. The underlying sequence is immutable, while the metdata and positional metadata are mutable. Parameters ---------- sequence : str, Sequence, or 1D np.ndarray (np.uint8 or '\\|S1') Characters representing the sequence itself. metadata : dict, optional Arbitrary metadata which applies to the entire sequence. A shallow copy of the ``dict`` will be made (see Examples section below for details). positional_metadata : pd.DataFrame consumable, optional Arbitrary per-character metadata (e.g., sequence read quality scores). Must be able to be passed directly to ``pd.DataFrame`` constructor. Each column of metadata must be the same length as `sequence`. A shallow copy of the positional metadata will be made if necessary (see Examples section below for details). interval_metadata : IntervalMetadata Arbitrary metadata which applies to intervals within a sequence to store interval features (such as genes, ncRNA on the sequence). lowercase : bool or str, optional If ``True``, lowercase sequence characters will be converted to uppercase characters. If ``False``, no characters will be converted. If a str, it will be treated as a key into the positional metadata of the object. All lowercase characters will be converted to uppercase, and a ``True`` value will be stored in a boolean array in the positional metadata under the key. See Also -------- DNA RNA Protein References ---------- .. [1] Nomenclature for incompletely specified bases in nucleic acid sequences: recommendations 1984. Nucleic Acids Res. May 10, 1985; 13(9): 3021-3030. A Cornish-Bowden Examples -------- >>> from pprint import pprint >>> from skbio import Sequence >>> from skbio.metadata import IntervalMetadata **Creating sequences:** Create a sequence without any metadata: >>> seq = Sequence('GGUCGUGAAGGA') >>> seq Sequence --------------- Stats: length: 12 --------------- 0 GGUCGUGAAG GA Create a sequence with metadata and positional metadata: >>> metadata = {'authors': ['Alice'], 'desc':'seq desc', 'id':'seq-id'} >>> positional_metadata = {'exons': [True, True, False, True], ... 'quality': [3, 3, 4, 10]} >>> interval_metadata = IntervalMetadata(4) >>> interval = interval_metadata.add([(1, 3)], metadata={'gene': 'sagA'}) >>> seq = Sequence('ACGT', metadata=metadata, ... positional_metadata=positional_metadata, ... interval_metadata=interval_metadata) >>> seq Sequence ----------------------------- Metadata: 'authors': 'desc': 'seq desc' 'id': 'seq-id' Positional metadata: 'exons': 'quality': Interval metadata: 1 interval feature Stats: length: 4 ----------------------------- 0 ACGT **Retrieving underlying sequence data:** Retrieve underlying sequence: >>> seq.values # doctest: +NORMALIZE_WHITESPACE array([b'A', b'C', b'G', b'T'], dtype='|S1') Underlying sequence immutable: >>> values = np.array([b'T', b'C', b'G', b'A'], dtype='|S1') >>> seq.values = values # doctest: +SKIP Traceback (most recent call last): ... AttributeError: can't set attribute >>> seq.values[0] = b'T' Traceback (most recent call last): ... ValueError: assignment destination is read-only **Retrieving sequence metadata:** Retrieve metadata: >>> pprint(seq.metadata) # using pprint to display dict in sorted order {'authors': ['Alice'], 'desc': 'seq desc', 'id': 'seq-id'} Retrieve positional metadata: >>> seq.positional_metadata exons quality 0 True 3 1 True 3 2 False 4 3 True 10 Retrieve interval metadata: >>> seq.interval_metadata # doctest: +ELLIPSIS 1 interval feature ------------------ Interval(interval_metadata=<...>, bounds=[(1, 3)], \ fuzzy=[(False, False)], metadata={'gene': 'sagA'}) **Updating sequence metadata:** .. warning:: Be aware that a shallow copy of ``metadata`` and ``positional_metadata`` is made for performance. Since a deep copy is not made, changes made to mutable Python objects stored as metadata may affect the metadata of other ``Sequence`` objects or anything else that shares a reference to the object. The following examples illustrate this behavior. First, let's create a sequence and update its metadata: >>> metadata = {'id':'seq-id', 'desc':'seq desc', 'authors': ['Alice']} >>> seq = Sequence('ACGT', metadata=metadata) >>> seq.metadata['id'] = 'new-id' >>> seq.metadata['pubmed'] = 12345 >>> pprint(seq.metadata) {'authors': ['Alice'], 'desc': 'seq desc', 'id': 'new-id', 'pubmed': 12345} Note that the original metadata dictionary (stored in variable ``metadata``) hasn't changed because a shallow copy was made: >>> pprint(metadata) {'authors': ['Alice'], 'desc': 'seq desc', 'id': 'seq-id'} >>> seq.metadata == metadata False Note however that since only a *shallow* copy was made, updates to mutable objects will also change the original metadata dictionary: >>> seq.metadata['authors'].append('Bob') >>> seq.metadata['authors'] ['Alice', 'Bob'] >>> metadata['authors'] ['Alice', 'Bob'] This behavior can also occur when manipulating a sequence that has been derived from another sequence: >>> subseq = seq[1:3] >>> subseq Sequence ----------------------------- Metadata: 'authors': 'desc': 'seq desc' 'id': 'new-id' 'pubmed': 12345 Stats: length: 2 ----------------------------- 0 CG >>> pprint(subseq.metadata) {'authors': ['Alice', 'Bob'], 'desc': 'seq desc', 'id': 'new-id', 'pubmed': 12345} The subsequence has inherited the metadata of its parent sequence. If we update the subsequence's author list, we see the changes propagated in the parent sequence and original metadata dictionary: >>> subseq.metadata['authors'].append('Carol') >>> subseq.metadata['authors'] ['Alice', 'Bob', 'Carol'] >>> seq.metadata['authors'] ['Alice', 'Bob', 'Carol'] >>> metadata['authors'] ['Alice', 'Bob', 'Carol'] The behavior for updating positional metadata is similar. Let's create a new sequence with positional metadata that is already stored in a ``pd.DataFrame``: >>> positional_metadata = pd.DataFrame( ... {'list': [[], [], [], []], 'quality': [3, 3, 4, 10]}) >>> seq = Sequence('ACGT', positional_metadata=positional_metadata) >>> seq Sequence ----------------------------- Positional metadata: 'list': 'quality': Stats: length: 4 ----------------------------- 0 ACGT >>> seq.positional_metadata list quality 0 [] 3 1 [] 3 2 [] 4 3 [] 10 Now let's update the sequence's positional metadata by adding a new column and changing a value in another column: >>> seq.positional_metadata['gaps'] = [False, False, False, False] >>> seq.positional_metadata.loc[0, 'quality'] = 999 >>> seq.positional_metadata list quality gaps 0 [] 999 False 1 [] 3 False 2 [] 4 False 3 [] 10 False Note that the original positional metadata (stored in variable ``positional_metadata``) hasn't changed because a shallow copy was made: >>> positional_metadata list quality 0 [] 3 1 [] 3 2 [] 4 3 [] 10 >>> seq.positional_metadata.equals(positional_metadata) False Next let's create a sequence that has been derived from another sequence: >>> subseq = seq[1:3] >>> subseq Sequence ----------------------------- Positional metadata: 'list': 'quality': 'gaps': Stats: length: 2 ----------------------------- 0 CG >>> subseq.positional_metadata list quality gaps 0 [] 3 False 1 [] 4 False As described above for metadata, since only a *shallow* copy was made of the positional metadata, updates to mutable objects will also change the parent sequence's positional metadata and the original positional metadata ``pd.DataFrame``: >>> subseq.positional_metadata.loc[0, 'list'].append('item') >>> subseq.positional_metadata list quality gaps 0 [item] 3 False 1 [] 4 False >>> seq.positional_metadata list quality gaps 0 [] 999 False 1 [item] 3 False 2 [] 4 False 3 [] 10 False >>> positional_metadata list quality 0 [] 3 1 [item] 3 2 [] 4 3 [] 10 You can also update the interval metadata. Let's re-create a ``Sequence`` object with interval metadata at first: >>> seq = Sequence('ACGT') >>> interval = seq.interval_metadata.add( ... [(1, 3)], metadata={'gene': 'foo'}) You can update directly on the ``Interval`` object: >>> interval # doctest: +ELLIPSIS Interval(interval_metadata=<...>, bounds=[(1, 3)], \ fuzzy=[(False, False)], metadata={'gene': 'foo'}) >>> interval.bounds = [(0, 2)] >>> interval # doctest: +ELLIPSIS Interval(interval_metadata=<...>, bounds=[(0, 2)], \ fuzzy=[(False, False)], metadata={'gene': 'foo'}) You can also query and obtain the interval features you are interested and then modify them: >>> intervals = list(seq.interval_metadata.query(metadata={'gene': 'foo'})) >>> intervals[0].fuzzy = [(True, False)] >>> print(intervals[0]) # doctest: +ELLIPSIS Interval(interval_metadata=<...>, bounds=[(0, 2)], \ fuzzy=[(True, False)], metadata={'gene': 'foo'}) """ _number_of_extended_ascii_codes = 256 # ASCII is built such that the difference between uppercase and lowercase # is the 6th bit. _ascii_invert_case_bit_offset = 32 _ascii_lowercase_boundary = 90 default_write_format = 'fasta' __hash__ = None @property @stable(as_of="0.4.0") def values(self): r"""Array containing underlying sequence characters. Notes ----- This property is not writeable. Examples -------- >>> from skbio import Sequence >>> s = Sequence('AACGA') >>> s.values # doctest: +NORMALIZE_WHITESPACE array([b'A', b'A', b'C', b'G', b'A'], dtype='|S1') """ return self._bytes.view('|S1') @property def __array_interface__(self): r"""Array interface for compatibility with numpy. This property allows a ``Sequence`` object to share its underlying data buffer (``Sequence.values``) with numpy. See [1]_ for more details. References ---------- .. [1] http://docs.scipy.org/doc/numpy/reference/arrays.interface.html Examples -------- >>> import numpy as np >>> from skbio import Sequence >>> seq = Sequence('ABC123') >>> np.asarray(seq) # doctest: +NORMALIZE_WHITESPACE array([b'A', b'B', b'C', b'1', b'2', b'3'], dtype='|S1') """ return self.values.__array_interface__ @property @experimental(as_of="0.4.1") def observed_chars(self): r"""Set of observed characters in the sequence. Notes ----- This property is not writeable. Examples -------- >>> from skbio import Sequence >>> s = Sequence('AACGAC') >>> s.observed_chars == {'G', 'A', 'C'} True """ return set(str(self)) @property def _string(self): return self._bytes.tobytes() @classonlymethod @experimental(as_of="0.4.1") def concat(cls, sequences, how='strict'): r"""Concatenate an iterable of ``Sequence`` objects. Parameters ---------- sequences : iterable (Sequence) An iterable of ``Sequence`` objects or appropriate subclasses. how : {'strict', 'inner', 'outer'}, optional How to intersect the `positional_metadata` of the sequences. If 'strict': the `positional_metadata` must have the exact same columns; 'inner': an inner-join of the columns (only the shared set of columns are used); 'outer': an outer-join of the columns (all columns are used: missing values will be padded with NaN). Returns ------- Sequence The returned sequence will be an instance of the class which called this class-method. Raises ------ ValueError If `how` is not one of: 'strict', 'inner', or 'outer'. ValueError If `how` is 'strict' and the `positional_metadata` of each sequence does not have the same columns. TypeError If the sequences cannot be cast as the calling class. Notes ----- The sequence-wide metadata (``Sequence.metadata``) is not retained during concatenation. Sequence objects can be cast to a different type only when the new type is an ancestor or child of the original type. Casting between sibling types is not allowed, e.g. ``DNA`` -> ``RNA`` is not allowed, but ``DNA`` -> ``Sequence`` or ``Sequence`` -> ``DNA`` would be. Examples -------- Concatenate two DNA sequences into a new DNA object: >>> from skbio import DNA, Sequence >>> s1 = DNA("ACGT") >>> s2 = DNA("GGAA") >>> DNA.concat([s1, s2]) DNA -------------------------- Stats: length: 8 has gaps: False has degenerates: False has definites: True GC-content: 50.00% -------------------------- 0 ACGTGGAA Concatenate DNA sequences into a Sequence object (type coercion): >>> Sequence.concat([s1, s2]) Sequence ------------- Stats: length: 8 ------------- 0 ACGTGGAA Positional metadata is conserved: >>> s1 = DNA('AcgT', lowercase='one') >>> s2 = DNA('GGaA', lowercase='one', ... positional_metadata={'two': [1, 2, 3, 4]}) >>> result = DNA.concat([s1, s2], how='outer') >>> result DNA --------------------------- Positional metadata: 'one': 'two': Stats: length: 8 has gaps: False has degenerates: False has definites: True GC-content: 50.00% --------------------------- 0 ACGTGGAA >>> result.positional_metadata one two 0 False NaN 1 True NaN 2 True NaN 3 False NaN 4 False 1.0 5 False 2.0 6 True 3.0 7 False 4.0 """ if how not in {'strict', 'inner', 'outer'}: raise ValueError("`how` must be 'strict', 'inner', or 'outer'.") seqs = list(sequences) if len(seqs) == 0: return cls("") for seq in seqs: seq._assert_can_cast_to(cls) if how == 'strict': how = 'inner' cols = set() for s in seqs: if s.has_positional_metadata(): cols.add(frozenset(s.positional_metadata)) else: cols.add(frozenset()) if len(cols) > 1: raise ValueError("The positional metadata of the sequences do" " not have matching columns. Consider setting" " how='inner' or how='outer'") seq_data = [] pm_data = [] for seq in seqs: seq_data.append(seq._bytes) pm_data.append(seq.positional_metadata) if not seq.has_positional_metadata(): del seq.positional_metadata pm = pd.concat(pm_data, join=how, ignore_index=True, sort=True) bytes_ = np.concatenate(seq_data) im = IntervalMetadata.concat(i.interval_metadata for i in seqs) return cls(bytes_, positional_metadata=pm, interval_metadata=im) @classmethod def _assert_can_cast_to(cls, target): if not (issubclass(cls, target) or issubclass(target, cls)): raise TypeError("Cannot cast %r as %r." % (cls.__name__, target.__name__)) @overrides(PositionalMetadataMixin) def _positional_metadata_axis_len_(self): return len(self) @overrides(IntervalMetadataMixin) def _interval_metadata_axis_len_(self): return len(self) @stable(as_of="0.4.0") def __init__(self, sequence, metadata=None, positional_metadata=None, interval_metadata=None, lowercase=False): if isinstance(sequence, np.ndarray): if sequence.dtype == np.uint8: self._set_bytes_contiguous(sequence) elif sequence.dtype == '|S1': sequence = sequence.view(np.uint8) # Guarantee the sequence is an array (might be scalar before # this). if sequence.shape == (): sequence = np.array([sequence], dtype=np.uint8) self._set_bytes_contiguous(sequence) else: raise TypeError( "Can only create sequence from numpy.ndarray of dtype " "np.uint8 or '|S1'. Invalid dtype: %s" % sequence.dtype) elif isinstance(sequence, Sequence): # Sequence casting is acceptable between direct # decendants/ancestors sequence._assert_can_cast_to(type(self)) if metadata is None and sequence.has_metadata(): metadata = sequence.metadata if (positional_metadata is None and sequence.has_positional_metadata()): positional_metadata = sequence.positional_metadata if (interval_metadata is None and sequence.has_interval_metadata()): interval_metadata = sequence.interval_metadata sequence = sequence._bytes self._owns_bytes = False self._set_bytes(sequence) else: # Encode as ascii to raise UnicodeEncodeError if necessary. if isinstance(sequence, str): sequence = sequence.encode("ascii") s = np.frombuffer(sequence, dtype=np.uint8) # There are two possibilities (to our knowledge) at this point: # Either the sequence we were given was something string-like, # (else it would not have made it past frombuffer), or it was a # numpy scalar, and so our length must be 1. if isinstance(sequence, np.generic) and len(s) != 1: raise TypeError("Can cannot create a sequence with %r" % type(sequence).__name__) sequence = s self._owns_bytes = False self._set_bytes(sequence) MetadataMixin._init_(self, metadata=metadata) PositionalMetadataMixin._init_( self, positional_metadata=positional_metadata) IntervalMetadataMixin._init_( self, interval_metadata=interval_metadata) if lowercase is False: pass elif lowercase is True or isinstance(lowercase, str): lowercase_mask = self._bytes > self._ascii_lowercase_boundary self._convert_to_uppercase(lowercase_mask) # If it isn't True, it must be a string_type if not (lowercase is True): self.positional_metadata[lowercase] = lowercase_mask else: raise TypeError("lowercase keyword argument expected a bool or " "string, but got %s" % type(lowercase)) def _set_bytes_contiguous(self, sequence): r"""Munge the sequence data into a numpy array of dtype uint8.""" if not sequence.flags['C_CONTIGUOUS']: # numpy doesn't support views of non-contiguous arrays. Since we're # making heavy use of views internally, and users may also supply # us with a view, make sure we *always* store a contiguous array to # avoid hard-to-track bugs. See # https://github.com/numpy/numpy/issues/5716 sequence = np.ascontiguousarray(sequence) self._owns_bytes = True else: self._owns_bytes = False self._set_bytes(sequence) def _set_bytes(self, sequence): sequence.flags.writeable = False self._bytes = sequence def _convert_to_uppercase(self, lowercase): if np.any(lowercase): with self._byte_ownership(): self._bytes[lowercase] ^= self._ascii_invert_case_bit_offset @stable(as_of="0.4.0") def __contains__(self, subsequence): r"""Determine if a subsequence is contained in this sequence. Parameters ---------- subsequence : str, Sequence, or 1D np.ndarray (np.uint8 or '\|S1') The putative subsequence. Returns ------- bool Indicates whether `subsequence` is contained in this sequence. Raises ------ TypeError If `subsequence` is a ``Sequence`` object with a different type than this sequence. Examples -------- >>> from skbio import Sequence >>> s = Sequence('GGUCGUGAAGGA') >>> 'GGU' in s True >>> 'CCC' in s False """ return self._munge_to_bytestring(subsequence, "in") in self._string @stable(as_of="0.4.0") def __eq__(self, other): r"""Determine if this sequence is equal to another. Sequences are equal if they are *exactly* the same type and their sequence characters, metadata, and positional metadata are the same. Parameters ---------- other : Sequence Sequence to test for equality against. Returns ------- bool Indicates whether this sequence is equal to `other`. Examples -------- Define two ``Sequence`` objects that have the same underlying sequence of characters: >>> from skbio import Sequence >>> s = Sequence('ACGT') >>> t = Sequence('ACGT') The two sequences are considered equal because they are the same type, their underlying sequence of characters are the same, and their optional metadata attributes (``metadata`` and ``positional_metadata``) were not provided: >>> s == t True >>> t == s True Define another sequence object with a different sequence of characters than the previous two sequence objects: >>> u = Sequence('ACGA') >>> u == t False Define a sequence with the same sequence of characters as ``u`` but with different metadata, positional metadata, and interval metadata: >>> v = Sequence('ACGA', metadata={'id': 'abc'}, ... positional_metadata={'quality':[1, 5, 3, 3]}) >>> _ = v.interval_metadata.add([(0, 1)]) The two sequences are not considered equal because their metadata, positional metadata, and interval metadata do not match: >>> u == v False """ # checks ordered from least to most expensive if self.__class__ != other.__class__: return False if not MetadataMixin._eq_(self, other): return False if self._string != other._string: return False if not PositionalMetadataMixin._eq_(self, other): return False if not IntervalMetadataMixin._eq_(self, other): return False return True @stable(as_of="0.4.0") def __ne__(self, other): r"""Determine if this sequence is not equal to another. Sequences are not equal if they are not *exactly* the same type, or their sequence characters, metadata, or positional metadata differ. Parameters ---------- other : Sequence Sequence to test for inequality against. Returns ------- bool Indicates whether this sequence is not equal to `other`. Examples -------- >>> from skbio import Sequence >>> s = Sequence('ACGT') >>> t = Sequence('ACGT') >>> s != t False >>> u = Sequence('ACGA') >>> u != t True >>> v = Sequence('ACGA', metadata={'id': 'v'}) >>> u != v True """ return not (self == other) @stable(as_of="0.4.0") def __getitem__(self, indexable): r"""Slice this sequence. Notes ----- This drops the ``self.interval_metadata`` from the returned new ``Sequence`` object. Parameters ---------- indexable : int, slice, iterable (int and slice), 1D array_like (bool) The position(s) to return from this sequence. If `indexable` is an iterable of integers, these are assumed to be indices in the sequence to keep. If `indexable` is a 1D ``array_like`` of booleans, these are assumed to be the positions in the sequence to keep. Returns ------- Sequence New sequence containing the position(s) specified by `indexable` in this sequence. Positional metadata will be sliced in the same manner and included in the returned sequence. `metadata` is included in the returned sequence. Examples -------- >>> from skbio import Sequence >>> s = Sequence('GGUCGUGAAGGA') Obtain a single character from the sequence: >>> s[1] Sequence ------------- Stats: length: 1 ------------- 0 G Obtain a slice: >>> s[7:] Sequence ------------- Stats: length: 5 ------------- 0 AAGGA Obtain characters at the following indices: >>> s[[3, 4, 7, 0, 3]] Sequence ------------- Stats: length: 5 ------------- 0 CGAGC Obtain characters at positions evaluating to `True`: >>> s = Sequence('GGUCG') >>> index = [True, False, True, 'a' is 'a', False] >>> s[index] Sequence ------------- Stats: length: 3 ------------- 0 GUC """ if (not isinstance(indexable, np.ndarray) and ((not isinstance(indexable, str)) and hasattr(indexable, '__iter__'))): indexable_ = indexable indexable = np.asarray(indexable) if indexable.dtype == object: indexable = list(indexable_) # TODO: Don't blow out memory if len(indexable) == 0: # indexing with an empty list, so convert to ndarray and # fall through to ndarray slicing below indexable = np.asarray(indexable) else: seq = np.concatenate( list(_slices_from_iter(self._bytes, indexable))) index = _as_slice_if_single_index(indexable) positional_metadata = None if self.has_positional_metadata(): pos_md_slices = list(_slices_from_iter( self.positional_metadata, index)) positional_metadata = pd.concat(pos_md_slices, sort=True) metadata = None if self.has_metadata(): metadata = self.metadata return self._constructor( sequence=seq, metadata=metadata, positional_metadata=positional_metadata) elif (isinstance(indexable, str) or isinstance(indexable, bool)): raise IndexError("Cannot index with %s type: %r" % (type(indexable).__name__, indexable)) if (isinstance(indexable, np.ndarray) and indexable.dtype == bool and len(indexable) != len(self)): raise IndexError("An boolean vector index must be the same length" " as the sequence (%d, not %d)." % (len(self), len(indexable))) if isinstance(indexable, np.ndarray) and indexable.size == 0: # convert an empty ndarray to a supported dtype for slicing a numpy # array indexable = indexable.astype(int) seq = self._bytes[indexable] positional_metadata = self._slice_positional_metadata(indexable) metadata = None if self.has_metadata(): metadata = self.metadata return self._constructor( sequence=seq, metadata=metadata, positional_metadata=positional_metadata) def _slice_positional_metadata(self, indexable): if self.has_positional_metadata(): if _is_single_index(indexable): index = _single_index_to_slice(indexable) else: index = indexable return self.positional_metadata.iloc[index] else: return None @stable(as_of="0.4.0") def __len__(self): r"""Return the number of characters in this sequence. Returns ------- int The length of this sequence. Examples -------- >>> from skbio import Sequence >>> s = Sequence('GGUC') >>> len(s) 4 """ return self._bytes.size @stable(as_of="0.4.0") def __bool__(self): r"""Returns truth value (truthiness) of sequence. Returns ------- bool True if length of sequence is greater than 0, else False. Examples -------- >>> from skbio import Sequence >>> bool(Sequence('')) False >>> bool(Sequence('ACGT')) True """ return len(self) > 0 @stable(as_of="0.4.0") def __iter__(self): r"""Iterate over positions in this sequence. Yields ------ Sequence Single character subsequence, one for each position in the sequence. Examples -------- >>> from skbio import Sequence >>> s = Sequence('GGUC') >>> for c in s: ... str(c) 'G' 'G' 'U' 'C' """ for i in range(len(self)): yield self[i] @stable(as_of="0.4.0") def __reversed__(self): r"""Iterate over positions in this sequence in reverse order. Yields ------ Sequence Single character subsequence, one for each position in the sequence. Examples -------- >>> from skbio import Sequence >>> s = Sequence('GGUC') >>> for c in reversed(s): ... str(c) 'C' 'U' 'G' 'G' """ return iter(self[::-1]) @stable(as_of="0.4.0") def __str__(self): r"""Return sequence characters as a string. Returns ------- str Sequence characters as a string. No metadata or positional metadata will be included. See Also -------- sequence Examples -------- >>> from skbio import Sequence >>> s = Sequence('GGUCGUAAAGGA', metadata={'id':'hello'}) >>> str(s) 'GGUCGUAAAGGA' """ return str(self._string.decode("ascii")) @stable(as_of="0.4.0") def __repr__(self): r"""Return a string representation of this sequence object. Representation includes: * sequence type * metadata keys and values: will display key/value if it is an understood type, otherwise just the type will be displayed. If it is an understood type whose representation is too long, just the type will be displayed * positional metadata: column names and column dtypes will be displayed in the order they appear in the positional metadata ``pd.DataFrame``. Column names (i.e., keys) follow the same display rules as metadata keys * interval metadata: the number of interval features will be displayed. * sequence stats (e.g., length) * up to five lines of chunked sequence data. Each line of chunked sequence data displays the current position in the sequence Returns ------- str String representation of this sequence object. Notes ----- Subclasses can override Sequence._repr_stats to provide custom statistics. Examples -------- Short sequence without metadata: >>> from skbio import Sequence >>> from skbio.metadata._interval import IntervalMetadata >>> Sequence('ACGTAATGGATACGTAATGCA') Sequence ------------------------- Stats: length: 21 ------------------------- 0 ACGTAATGGA TACGTAATGC A Longer sequence displays first two lines and last two lines: >>> Sequence('ACGT' * 100) Sequence --------------------------------------------------------------------- Stats: length: 400 --------------------------------------------------------------------- 0 ACGTACGTAC GTACGTACGT ACGTACGTAC GTACGTACGT ACGTACGTAC GTACGTACGT 60 ACGTACGTAC GTACGTACGT ACGTACGTAC GTACGTACGT ACGTACGTAC GTACGTACGT ... 300 ACGTACGTAC GTACGTACGT ACGTACGTAC GTACGTACGT ACGTACGTAC GTACGTACGT 360 ACGTACGTAC GTACGTACGT ACGTACGTAC GTACGTACGT Sequence with metadata, positional metadata, and interval metadata: >>> metadata = { ... 'id': 'seq-id', ... 'description': 'description of the sequence, wrapping across ' ... 'lines if it\'s too long', ... 'authors': ['Alice', 'Bob', 'Carol'], ... 'year': 2015, ... 'published': True ... } >>> positional_metadata = { ... 'exons': [True, True, False, True], ... 'quality': [3, 10, 11, 10] ... } >>> seq = Sequence('ACGT', metadata=metadata, ... positional_metadata=positional_metadata) >>> _ = seq.interval_metadata.add([(0, 2)], metadata={'gene': 'sagA'}) >>> seq Sequence ---------------------------------------------------------------------- Metadata: 'authors': 'description': "description of the sequence, wrapping across lines if it's too long" 'id': 'seq-id' 'published': True 'year': 2015 Positional metadata: 'exons': 'quality': Interval metadata: 1 interval feature Stats: length: 4 ---------------------------------------------------------------------- 0 ACGT """ return _SequenceReprBuilder( seq=self, width=71, # 79 for pep8, 8 space indent for docstrings indent=4, chunk_size=10).build() def _repr_stats(self): r"""Define statistics to display in the sequence's repr. Subclasses can override this method to provide type-specific statistics. This method computes a single statistic: length. Returns ------- list List of tuples where each tuple represents a statistic. Each tuple contains exactly two ``str`` elements: the statistic's name/label, and the str-formatted value of the statistic. Ordering of statistics (i.e., list order) determines display order in the sequence repr. """ return [('length', '%d' % len(self))] @stable(as_of="0.4.0") def __copy__(self): r"""Return a shallow copy of this sequence. See Also -------- copy Notes ----- This method is equivalent to ``seq.copy(deep=False)``. """ return self._copy(False, {}) @stable(as_of="0.4.0") def __deepcopy__(self, memo): r"""Return a deep copy of this sequence. See Also -------- copy Notes ----- This method is equivalent to ``seq.copy(deep=True)``. """ return self._copy(True, memo) def _copy(self, deep, memo): # strategy: copy the sequence without metadata first, then set metadata # attributes with copies. we take this approach instead of simply # passing the metadata through the Sequence constructor because we # don't want to copy twice (this could happen when deep=True, where we # deep copy here and then shallow copy in the Sequence constructor). we # also directly set the private metadata attributes instead of using # their public setters to avoid an unnecessary copy # we don't make a distinction between deep vs. shallow copy of bytes # because dtype=np.uint8. we only need to make the distinction when # dealing with object dtype bytes_ = np.copy(self._bytes) seq_copy = self._constructor(sequence=bytes_, metadata=None, positional_metadata=None, interval_metadata=None) if deep: seq_copy._metadata = MetadataMixin._deepcopy_(self, memo) seq_copy._positional_metadata = \ PositionalMetadataMixin._deepcopy_(self, memo) seq_copy._interval_metadata = IntervalMetadataMixin._deepcopy_( self, memo) else: seq_copy._metadata = MetadataMixin._copy_(self) seq_copy._positional_metadata = \ PositionalMetadataMixin._copy_(self) seq_copy._interval_metadata = IntervalMetadataMixin._copy_( self) return seq_copy @stable(as_of='0.4.0') def lowercase(self, lowercase): r"""Return a case-sensitive string representation of the sequence. Parameters ---------- lowercase: str or boolean vector If lowercase is a boolean vector, it is used to set sequence characters to lowercase in the output string. True values in the boolean vector correspond to lowercase characters. If lowercase is a str, it is treated like a key into the positional metadata, pointing to a column which must be a boolean vector. That boolean vector is then used as described previously. Returns ------- str String representation of sequence with specified characters set to lowercase. Examples -------- >>> from skbio import Sequence >>> s = Sequence('ACGT') >>> s.lowercase([True, True, False, False]) 'acGT' >>> s = Sequence('ACGT', ... positional_metadata={ ... 'exons': [True, False, False, True]}) >>> s.lowercase('exons') 'aCGt' Constructor automatically populates a column in positional metadata when the ``lowercase`` keyword argument is provided with a column name: >>> s = Sequence('ACgt', lowercase='introns') >>> s.lowercase('introns') 'ACgt' >>> s = Sequence('ACGT', lowercase='introns') >>> s.lowercase('introns') 'ACGT' """ index = self._munge_to_index_array(lowercase) outbytes = self._bytes.copy() outbytes[index] ^= self._ascii_invert_case_bit_offset return str(outbytes.tobytes().decode('ascii')) @stable(as_of="0.4.0") def count(self, subsequence, start=None, end=None): r"""Count occurrences of a subsequence in this sequence. Parameters ---------- subsequence : str, Sequence, or 1D np.ndarray (np.uint8 or '\|S1') Subsequence to count occurrences of. start : int, optional The position at which to start counting (inclusive). end : int, optional The position at which to stop counting (exclusive). Returns ------- int Number of occurrences of `subsequence` in this sequence. Raises ------ ValueError If `subsequence` is of length 0. TypeError If `subsequence` is a ``Sequence`` object with a different type than this sequence. Examples -------- >>> from skbio import Sequence >>> s = Sequence('GGUCG') >>> s.count('G') 3 >>> s.count('GG') 1 >>> s.count('T') 0 >>> s.count('G', 2, 5) 1 """ if len(subsequence) == 0: raise ValueError("`count` is not defined for empty subsequences.") return self._string.count( self._munge_to_bytestring(subsequence, "count"), start, end) @experimental(as_of="0.5.0") def replace(self, where, character): r"""Replace values in this sequence with a different character. Parameters ---------- where : 1D array_like (bool) or iterable (slices or ints) or str Indicates positions in the sequence to replace with `character`. Can be a boolean vector, an iterable of indices/slices, or a string that is a key in `positional_metadata` pointing to a boolean vector. character : str or bytes Character that will replace chosen items in this sequence. Returns ------- Sequence Copy of this sequence, with chosen items replaced with chosen character. All metadata is retained. Examples -------- Let's create and display a Sequence: >>> from skbio import Sequence >>> sequence = Sequence('GGTACCAACG') >>> str(sequence) 'GGTACCAACG' Let's call ``replace`` on the Sequence using a boolean vector for ``where`` and assign it to a new variable: >>> seq = sequence.replace([False, False, False, True, False, False, ... True, True, False, False], '-') Let's take a look at the new Sequence: >>> str(seq) 'GGT-CC--CG' Other types of input are accepted by the ``where`` parameter. Let's pass in a list of indices and slices that is equivalent to the boolean vector we used previously: >>> str(seq) == str(sequence.replace([3, slice(6, 8)], '-')) True ``where`` also accepts a boolean vector contained in ``Sequence.positional_metadata``: >>> sequence.positional_metadata = {'where': ... [False, False, False, True, False, ... False, True, True, False, False]} Let's pass in the key ``'where'`` and compare to ``seq``: >>> str(seq) == str(sequence.replace('where', '-')) True """ if type(character) is not bytes: character = character.encode('ascii') character = ord(character) index = self._munge_to_index_array(where) seq_bytes = self._bytes.copy() seq_bytes[index] = character metadata = None if self.has_metadata(): metadata = self.metadata positional_metadata = None if self.has_positional_metadata(): positional_metadata = self.positional_metadata interval_metadata = None if self.has_interval_metadata(): interval_metadata = self.interval_metadata # Use __class__ instead of _constructor so that validations are # performed for subclasses (the user could have introduced invalid # characters). return self.__class__(seq_bytes, metadata=metadata, positional_metadata=positional_metadata, interval_metadata=interval_metadata) @stable(as_of="0.4.0") def index(self, subsequence, start=None, end=None): r"""Find position where subsequence first occurs in the sequence. Parameters ---------- subsequence : str, Sequence, or 1D np.ndarray (np.uint8 or '\|S1') Subsequence to search for in this sequence. start : int, optional The position at which to start searching (inclusive). end : int, optional The position at which to stop searching (exclusive). Returns ------- int Position where `subsequence` first occurs in this sequence. Raises ------ ValueError If `subsequence` is not present in this sequence. TypeError If `subsequence` is a ``Sequence`` object with a different type than this sequence. Examples -------- >>> from skbio import Sequence >>> s = Sequence('ACACGACGTT-') >>> s.index('ACG') 2 """ try: return self._string.index( self._munge_to_bytestring(subsequence, "index"), start, end) except ValueError: raise ValueError( "%r is not present in %r." % (subsequence, self)) @experimental(as_of="0.4.0") def distance(self, other, metric=None): r"""Compute the distance to another sequence. Parameters ---------- other : str, Sequence, or 1D np.ndarray (np.uint8 or '\|S1') Sequence to compute the distance to. If `other` is a ``Sequence`` object, it must be the same type as this sequence. Other input types will be converted into a ``Sequence`` object of the same type as this sequence. metric : function, optional Function used to compute the distance between this sequence and `other`. If ``None`` (the default), Hamming distance will be used (:func:`skbio.sequence.distance.hamming`). `metric` should take two ``skbio.Sequence`` objects and return a ``float``. The sequence objects passed to `metric` will be the same type as this sequence. See :mod:`skbio.sequence.distance` for other predefined metrics that can be supplied via `metric`. Returns ------- float Distance between this sequence and `other` as defined by `metric`. Raises ------ TypeError If `other` is a ``Sequence`` object with a different type than this sequence. See Also -------- skbio.sequence.distance fraction_diff fraction_same Examples -------- >>> from skbio import Sequence >>> s = Sequence('GGUC') >>> t = Sequence('AGUC') Compute Hamming distance (the default metric): >>> s.distance(t) 0.25 Use a custom metric: >>> def custom_metric(s1, s2): return 0.42 >>> s.distance(t, custom_metric) 0.42 """ # TODO refactor this method to accept a name (string) of the distance # metric to apply and accept **kwargs other = self._munge_to_self_type(other, 'distance') if metric is None: metric = skbio.sequence.distance.hamming return float(metric(self, other)) @stable(as_of="0.4.0") def matches(self, other): r"""Find positions that match with another sequence. Parameters ---------- other : str, Sequence, or 1D np.ndarray (np.uint8 or '\|S1') Sequence to compare to. Returns ------- 1D np.ndarray (bool) Boolean vector where ``True`` at position ``i`` indicates a match between the sequences at their positions ``i``. Raises ------ ValueError If the sequences are not the same length. TypeError If `other` is a ``Sequence`` object with a different type than this sequence. See Also -------- mismatches Examples -------- >>> from skbio import Sequence >>> s = Sequence('GGUC') >>> t = Sequence('GAUU') >>> s.matches(t) array([ True, False, True, False], dtype=bool) """ other = self._munge_to_sequence(other, 'matches/mismatches') if len(self) != len(other): raise ValueError("Match and mismatch vectors can only be " "generated from equal length sequences.") return self._bytes == other._bytes @stable(as_of="0.4.0") def mismatches(self, other): r"""Find positions that do not match with another sequence. Parameters ---------- other : str, Sequence, or 1D np.ndarray (np.uint8 or '\|S1') Sequence to compare to. Returns ------- 1D np.ndarray (bool) Boolean vector where ``True`` at position ``i`` indicates a mismatch between the sequences at their positions ``i``. Raises ------ ValueError If the sequences are not the same length. TypeError If `other` is a ``Sequence`` object with a different type than this sequence. See Also -------- matches Examples -------- >>> from skbio import Sequence >>> s = Sequence('GGUC') >>> t = Sequence('GAUU') >>> s.mismatches(t) array([False, True, False, True], dtype=bool) """ return np.invert(self.matches(other)) @stable(as_of="0.4.0") def match_frequency(self, other, relative=False): r"""Return count of positions that are the same between two sequences. Parameters ---------- other : str, Sequence, or 1D np.ndarray (np.uint8 or '\|S1') Sequence to compare to. relative : bool, optional If ``True``, return the relative frequency of matches instead of the count. Returns ------- int or float Number of positions that are the same between the sequences. This will be an ``int`` if `relative` is ``False`` and a ``float`` if `relative` is ``True``. Raises ------ ValueError If the sequences are not the same length. TypeError If `other` is a ``Sequence`` object with a different type than this sequence. See Also -------- mismatch_frequency matches mismatches distance Examples -------- >>> from skbio import Sequence >>> s = Sequence('GGUC') >>> t = Sequence('AGUC') >>> s.match_frequency(t) 3 >>> s.match_frequency(t, relative=True) 0.75 """ if relative: return float(self.matches(other).mean()) else: return int(self.matches(other).sum()) @stable(as_of="0.4.0") def mismatch_frequency(self, other, relative=False): r"""Return count of positions that differ between two sequences. Parameters ---------- other : str, Sequence, or 1D np.ndarray (np.uint8 or '\|S1') Sequence to compare to. relative : bool, optional If ``True``, return the relative frequency of mismatches instead of the count. Returns ------- int or float Number of positions that differ between the sequences. This will be an ``int`` if `relative` is ``False`` and a ``float`` if `relative` is ``True``. Raises ------ ValueError If the sequences are not the same length. TypeError If `other` is a ``Sequence`` object with a different type than this sequence. See Also -------- match_frequency matches mismatches distance Examples -------- >>> from skbio import Sequence >>> s = Sequence('GGUC') >>> t = Sequence('AGUC') >>> s.mismatch_frequency(t) 1 >>> s.mismatch_frequency(t, relative=True) 0.25 """ if relative: return float(self.mismatches(other).mean()) else: return int(self.mismatches(other).sum()) @experimental(as_of="0.4.1") def frequencies(self, chars=None, relative=False): r"""Compute frequencies of characters in the sequence. Parameters ---------- chars : str or set of str, optional Characters to compute the frequencies of. May be a ``str`` containing a single character or a ``set`` of single-character strings. If ``None``, frequencies will be computed for all characters present in the sequence. relative : bool, optional If ``True``, return the relative frequency of each character instead of its count. If `chars` is provided, relative frequencies will be computed with respect to the number of characters in the sequence, **not** the total count of characters observed in `chars`. Thus, the relative frequencies will not necessarily sum to 1.0 if `chars` is provided. Returns ------- dict Frequencies of characters in the sequence. Raises ------ TypeError If `chars` is not a ``str`` or ``set`` of ``str``. ValueError If `chars` is not a single-character ``str`` or a ``set`` of single-character strings. ValueError If `chars` contains characters outside the allowable range of characters in a ``Sequence`` object. See Also -------- kmer_frequencies iter_kmers Notes ----- If the sequence is empty (i.e., length zero), ``relative=True``, **and** `chars` is provided, the relative frequency of each specified character will be ``np.nan``. If `chars` is not provided, this method is equivalent to, but faster than, ``seq.kmer_frequencies(k=1)``. If `chars` is not provided, it is equivalent to, but faster than, passing ``chars=seq.observed_chars``. Examples -------- Compute character frequencies of a sequence: >>> from pprint import pprint >>> from skbio import Sequence >>> seq = Sequence('AGAAGACC') >>> freqs = seq.frequencies() >>> pprint(freqs) # using pprint to display dict in sorted order {'A': 4, 'C': 2, 'G': 2} Compute relative character frequencies: >>> freqs = seq.frequencies(relative=True) >>> pprint(freqs) {'A': 0.5, 'C': 0.25, 'G': 0.25} Compute relative frequencies of characters A, C, and T: >>> freqs = seq.frequencies(chars={'A', 'C', 'T'}, relative=True) >>> pprint(freqs) {'A': 0.5, 'C': 0.25, 'T': 0.0} Note that since character T is not in the sequence we receive a relative frequency of 0.0. The relative frequencies of A and C are relative to the number of characters in the sequence (8), **not** the number of A and C characters (4 + 2 = 6). """ freqs = np.bincount(self._bytes, minlength=self._number_of_extended_ascii_codes) if chars is not None: chars, indices = self._chars_to_indices(chars) else: indices, = np.nonzero(freqs) # Downcast from int64 to uint8 then convert to str. This is safe # because we are guaranteed to have indices in the range 0 to 255 # inclusive. chars = indices.astype(np.uint8).tobytes().decode('ascii') obs_counts = freqs[indices] if relative: obs_counts = obs_counts / len(self) # Use tolist() for minor performance gain. return dict(zip(chars, obs_counts.tolist())) def _chars_to_indices(self, chars): r"""Helper for Sequence.frequencies.""" if isinstance(chars, (str, bytes)): chars = set([chars]) elif not isinstance(chars, set): raise TypeError( "`chars` must be of type `set`, not %r" % type(chars).__name__) # Impose an (arbitrary) ordering to `chars` so that we can return # `indices` in that same order. chars = list(chars) indices = [] for char in chars: if not isinstance(char, (str, bytes)): raise TypeError( "Each element of `chars` must be string-like, not %r" % type(char).__name__) if len(char) != 1: raise ValueError( "Each element of `chars` must contain a single " "character (found %d characters)" % len(char)) index = ord(char) if index >= self._number_of_extended_ascii_codes: raise ValueError( "Character %r in `chars` is outside the range of " "allowable characters in a `Sequence` object." % char) indices.append(index) return chars, indices @stable(as_of="0.4.0") def iter_kmers(self, k, overlap=True): r"""Generate kmers of length `k` from this sequence. Parameters ---------- k : int The kmer length. overlap : bool, optional Defines whether the kmers should be overlapping or not. Yields ------ Sequence kmer of length `k` contained in this sequence. Raises ------ ValueError If `k` is less than 1. Examples -------- >>> from skbio import Sequence >>> s = Sequence('ACACGACGTT') >>> for kmer in s.iter_kmers(4, overlap=False): ... str(kmer) 'ACAC' 'GACG' >>> for kmer in s.iter_kmers(3, overlap=True): ... str(kmer) 'ACA' 'CAC' 'ACG' 'CGA' 'GAC' 'ACG' 'CGT' 'GTT' """ if k < 1: raise ValueError("k must be greater than 0.") if overlap: step = 1 count = len(self) - k + 1 else: step = k count = len(self) // k if len(self) == 0 or self.has_positional_metadata(): # Slower path when sequence is empty or positional metadata needs # to be sliced. for i in range(0, len(self) - k + 1, step): yield self[i:i+k] else: # Optimized path when positional metadata doesn't need slicing. kmers = np.lib.stride_tricks.as_strided( self._bytes, shape=(k, count), strides=(1, step)).T metadata = None if self.has_metadata(): metadata = self.metadata for s in kmers: yield self._constructor( sequence=s, metadata=metadata, positional_metadata=None) @stable(as_of="0.4.0") def kmer_frequencies(self, k, overlap=True, relative=False): r"""Return counts of words of length `k` from this sequence. Parameters ---------- k : int The word length. overlap : bool, optional Defines whether the kmers should be overlapping or not. relative : bool, optional If ``True``, return the relative frequency of each kmer instead of its count. Returns ------- dict Frequencies of words of length `k` contained in this sequence. Raises ------ ValueError If `k` is less than 1. Examples -------- >>> from pprint import pprint >>> from skbio import Sequence >>> s = Sequence('ACACATTTATTA') >>> freqs = s.kmer_frequencies(3, overlap=False) >>> pprint(freqs) # using pprint to display dict in sorted order {'ACA': 1, 'CAT': 1, 'TTA': 2} >>> freqs = s.kmer_frequencies(3, relative=True, overlap=False) >>> pprint(freqs) {'ACA': 0.25, 'CAT': 0.25, 'TTA': 0.5} """ kmers = self.iter_kmers(k, overlap=overlap) freqs = dict(collections.Counter((str(seq) for seq in kmers))) if relative: if overlap: num_kmers = len(self) - k + 1 else: num_kmers = len(self) // k relative_freqs = {} for kmer, count in freqs.items(): relative_freqs[kmer] = count / num_kmers freqs = relative_freqs return freqs @stable(as_of="0.4.0") def find_with_regex(self, regex, ignore=None): r"""Generate slices for patterns matched by a regular expression. Parameters ---------- regex : str or regular expression object String to be compiled into a regular expression, or a pre- compiled regular expression object (e.g., from calling ``re.compile``). ignore : 1D array_like (bool) or iterable (slices or ints), optional Indicate the positions to ignore when matching. Yields ------ slice Location where the regular expression matched. Examples -------- >>> from skbio import Sequence >>> s = Sequence('AATATACCGGTTATAA') >>> for match in s.find_with_regex('(TATA+)'): ... match ... str(s[match]) slice(2, 6, None) 'TATA' slice(11, 16, None) 'TATAA' """ if isinstance(regex, str): regex = re.compile(regex) lookup = np.arange(len(self)) if ignore is None: string = str(self) else: ignore = self._munge_to_index_array(ignore) lookup = np.delete(lookup, ignore) string = str(self[lookup]) for match in regex.finditer(string): # We start at 1 because we don't want the group that contains all # other groups. for g in range(1, len(match.groups())+1): yield slice(lookup[match.start(g)], lookup[match.end(g) - 1] + 1) @stable(as_of="0.4.0") def iter_contiguous(self, included, min_length=1, invert=False): r"""Yield contiguous subsequences based on `included`. Parameters ---------- included : 1D array_like (bool) or iterable (slices or ints) `included` is transformed into a flat boolean vector where each position will either be included or skipped. All contiguous included positions will be yielded as a single region. min_length : int, optional The minimum length of a subsequence for it to be yielded. Default is 1. invert : bool, optional Whether to invert `included` such that it describes what should be skipped instead of included. Default is False. Yields ------ Sequence Contiguous subsequence as indicated by `included`. Notes ----- If slices provide adjacent ranges, then they will be considered the same contiguous subsequence. Examples -------- Here we use `iter_contiguous` to find all of the contiguous ungapped sequences using a boolean vector derived from our DNA sequence. >>> from skbio import DNA >>> s = DNA('AAA--TT-CCCC-G-') >>> no_gaps = ~s.gaps() >>> for ungapped_subsequence in s.iter_contiguous(no_gaps, ... min_length=2): ... print(ungapped_subsequence) AAA TT CCCC Note how the last potential subsequence was skipped because it would have been smaller than our `min_length` which was set to 2. We can also use `iter_contiguous` on a generator of slices as is produced by `find_motifs` (and `find_with_regex`). >>> from skbio import Protein >>> s = Protein('ACDFNASANFTACGNPNRTESL') >>> for subseq in s.iter_contiguous(s.find_motifs('N-glycosylation')): ... print(subseq) NASANFTA NRTE Note how the first subsequence contains two N-glycosylation sites. This happened because they were contiguous. """ idx = self._munge_to_index_array(included) if invert: idx = np.delete(np.arange(len(self)), idx) # Adapted from http://stackoverflow.com/a/7353335/579416 for contig in np.split(idx, np.where(np.diff(idx) != 1)[0] + 1): r = self[contig] if len(r) >= min_length: yield r def _constructor(self, **kwargs): return self.__class__(**kwargs) def _munge_to_index_array(self, sliceable): r"""Return an index array from something isomorphic to a boolean vector. """ if isinstance(sliceable, str): if sliceable in self.positional_metadata: if self.positional_metadata[sliceable].dtype == bool: sliceable = self.positional_metadata[sliceable] else: raise TypeError("Column '%s' in positional metadata does " "not correspond to a boolean vector" % sliceable) else: raise ValueError("No positional metadata associated with key " "'%s'" % sliceable) if not hasattr(sliceable, 'dtype') or (hasattr(sliceable, 'dtype') and sliceable.dtype == 'object'): sliceable = tuple(sliceable) bool_mode = False int_mode = False for s in sliceable: if isinstance(s, (bool, np.bool_)): bool_mode = True elif isinstance(s, (slice, int, np.signedinteger)) or ( hasattr(s, 'dtype') and s.dtype != bool): int_mode = True else: raise TypeError("Invalid type in iterable: %s, must be one" " of {bool, int, slice, np.signedinteger}" % s.__class__.__name__) if bool_mode and int_mode: raise TypeError("Cannot provide iterable of both bool and" " int.") sliceable = np.r_[sliceable] if sliceable.dtype == bool: if sliceable.size != len(self): raise ValueError("Boolean array (%d) does not match length of" " sequence (%d)." % (sliceable.size, len(self))) normalized, = np.where(sliceable) else: normalized = np.bincount(sliceable) if np.any(normalized > 1): raise ValueError("Overlapping index regions are not allowed.") normalized, = np.where(normalized) if np.any(normalized != sliceable): raise ValueError("Index regions are out of order.") return normalized def _munge_to_self_type(self, other, method): if isinstance(other, Sequence): if type(other) != type(self): raise TypeError("Cannot use %s and %s together with `%s`" % (self.__class__.__name__, other.__class__.__name__, method)) else: return other return self.__class__(other) def _munge_to_sequence(self, other, method): if isinstance(other, Sequence): if type(other) != type(self): raise TypeError("Cannot use %s and %s together with `%s`" % (self.__class__.__name__, other.__class__.__name__, method)) else: return other # We don't use self.__class__ or self._constructor here because we want # to construct the most general type of Sequence object in order to # avoid validation errors. return Sequence(other) def _munge_to_bytestring(self, other, method): if type(other) is bytes: return other elif isinstance(other, str): return other.encode('ascii') else: return self._munge_to_sequence(other, method)._string @contextmanager def _byte_ownership(self): if not self._owns_bytes: self._bytes = self._bytes.copy() self._owns_bytes = True self._bytes.flags.writeable = True yield self._bytes.flags.writeable = False def _single_index_to_slice(start_index): end_index = None if start_index == -1 else start_index+1 return slice(start_index, end_index) def _is_single_index(index): return (isinstance(index, numbers.Integral) and not isinstance(index, bool)) def _as_slice_if_single_index(indexable): if _is_single_index(indexable): return _single_index_to_slice(indexable) else: return indexable def _slices_from_iter(array, indexables): for i in indexables: if isinstance(i, slice): pass elif _is_single_index(i): i = _single_index_to_slice(i) else: raise IndexError("Cannot slice sequence from iterable " "containing %r." % i) yield array[i]