ENH: Add iter_regex and find_features

Fixes #612. Adds iter_regex to skbio._sequence.BiologicalSequence and
find_features to skbio._sequence.NucleotideSequence. Supported features are
currently purine_run and pyrimidine run, but more can be added by specifying
additional regular expressions.
I could be mistaken, but isn't aho-corasick intended for exact string
matching, and not as flexible as a regular expression?  If that is the
case, it would be practically impossible, for example, to locate
arbitrary-length runs of purines.

Your point about HMMs is well taken, and perhaps a solution would be to add
a parameter to find_features that takes a function as an argument (but
defaults to None) that is used when the caller needs to locate a type of
feature that's not in the hard-coded set of feature types.

As for the data representation of the actual _sequence attribute: I think
it makes sense to represent it as a numpy array, but this does not preclude
using string operations (e.g. regular expressions from python's builtin re
module) since, as I previously mentioned in response to your comment on the
matter, casting from a million-length numpy array to a python string was on
the order of ~2 microseconds in my tests, which shouldn't be a deal breaker.


On Wed, Aug 13, 2014 at 1:26 PM, Daniel McDonald <notifications@github.com>
wrote:

> In skbio/sequence/_sequence.py:
>
> > @@ -1324,6 +1350,34 @@ def reverse_complement(self):
> >          return self._complement(reversed(self))
> >      rc = reverse_complement
> >
> > +    def find_features(self, feature_type, min_length=1):
> > +        """Search the sequence for features
> > +
> > +        Parameters
> > +        ----------
> > +        feature_type : {'purine_run', 'pyrimidine_run'}
>
> Thanks @ebolyen <https://github.com/ebolyen>, that was fantastic. SCFGs
> have been widely used within RNA, particularly secondary structure
> prediction (you may be interested in consan
> <http://selab.janelia.org/software/consan/README>) and, AFAIK, many
> features of DNA/RNA/Protein can be described using CFGs. But, I do not know
> if all features can be, or if that is as efficient. For instance, HMMER
> uses HMMs to define coding features, and Infernal uses covariance models
> for defining non-coding RNAs for performance reasons, but possibly also due
> to the nature of the features they describe. I'm nearly certain that some
> features of interest, such as genes with alternative splicing or
> trans-regulatory effects cannot be represented with a regular grammar but I
> may be wrong in thought here. With respect to automatons though, we found
> in practice that Aho-Corasick was much faster than regex, while still
> allowing for multiple nodes and amb iguities . Regarding passing a specific
> function to produce a boolean vector for a simple FA, the above example is
> as verbose as:
>
> def pyrimidines_run(chars):
>     return array([c in 'ctu' for c in chars])
>
> But, I think we actually need to step back and decide two things:
>
>    - the level of complexity we want this method to support
>    - if the underlying data representation is or is not going to be a
>    numpy array in the future
>
> On the first point, if we want to be able to search an arbitrary sequence
> for a protein family, then we likely want to instead have support for HMMs
> and we probably should just wrap HMMER and save ourselves the investment to
> make it correct and fast. If we want to just be able to search for simple
> features, such as repeats (e.g., ARGARG...), then an automaton will be the
> most performant and verbose. If we want to just masks (e.g., all
> purine/pyrimidines), then boolean arrays will probably be the most
> performant due to SIMD.
>
> On the second point, if we use a numpy array, it makes boolean array ops
> and potentially other ops in BiologicalSequence very fast. It does not
> preclude using regex or other automatons, but we may need a translation
> layer. I really think that we should resolve the data representation first
> though.
>
> —
> Reply to this email directly or view it on GitHub
> <https://github.com/biocore/scikit-bio/pull/613/files#r16190630>.
>