This PR is ready for review.

In this PR, we add in Mutation and MutationSet classes. The intent is to represent mutations made to a SeqLike. I've taken inspiration from multiple places, but the biggest one has been the discussion on #57.

The use cases for this Mutation class are primarily two-fold:

1. Adding them to a SeqLike returns another mutated SeqLike.
2. Subtracting (diffing) a SeqLike from another yields a MutationSet.

While working through use case 2, I noticed how it's a bit tricky:

• adding a Mutation or MutationSet to a SeqLike results in a new SeqLike,
• but subtracting the new SeqLike from the original SeqLike might not necessarily result in the original MutationSet,
• yet adding the new MutationSet to the original SeqLike will give back a SeqLike identical to the new SeqLike.

An example of this phenomena is shown in a new notebook, docs/notebooks/mutations.ipynb. For reviewing purposes, that is probably the go-to notebook for understanding what Mutation and MutationSet classes can do; the rest are implementation details.

Would love to get feedback on this PR -- especially if there are semantics that I haven't yet thought of.

TODO list of what's left:

• [x] FIX TEST: SeqLike classes have two more class methods. Therefore, 77 is correct. Reference here.
• [x] Add .positions class method, which returns a list of positions to mutate.
• [x] Switch out magical_parse() for __new__() under Mutation.

I was playing around with the library (it was really nice and great job!) and hit this bug: instantiation like SeqLike('ATCGATC') or SeqLike('ATCGATC', None) will fail. From the doc: 1) seq_like is supposedly optional but is required in the implementation (code) and 2) logic regarding seq_like==None only occurs when sequence is a SeqLike (code) as opposed to a native type and thus seems not very useful? (There's a def determine__type() (code) that I think is intended for the job but it's not used.)

In addition, some general comments regarding dispatch, feel free to ignore them if they are out of context: def _construct_seqlike share the same signature for alphabet and codon_map, and def determine__type_and_alphabet share same signature for seq_type and sequence. Thus it seems clearer to extract the only signature that varies? E.g. for the latter case we only dispatch functions based on alphabet (or even directly use if alphabet==None: alphabet=determine_alphabet(_type, sequence) rather than dispatching). Current setup for these two functions seems to create duplicate logics and is error prone (types for sequence here should be the same?)

Again, happy to file a PR for anything specific above.

Hi team, thanks for creating this tool, it looks really nice!

Just a question regarding #noqa: DAR201: it's generated in the doc in quite a few places (e.g. https://modernatx.github.io/seqlike/reference/seqlike/#seqlike.SeqLike.SeqLike.deepcopy--noqa-dar101). I've tested that it can successfully be ommited via wrapping with <!-- --> (ref). Curious if the team likes the fix? Happy to submit a MR.

As mentioned by @ndousis in an earlier email thread, some of @dmnfarrell's code made it into this library. We would like to acknowledge @dmnfarrell as a contributor in the codebase.

@dmnfarrell is this something you would be amenable to? I would essentially ask the all-contributors bot to add your contribution in.

Prepending a string to SeqLike results unexpectedly in an appended version:

In [1]: from seqlike import SeqLike

In [2]: "ACTG" + SeqLike("TTTT", "nt", id="test")
Out[2]: 
*** NT: SeqRecord(seq=Seq('TTTTACTG'), id='test', name='<unknown name>', description='<unknown description>', dbxrefs=[]) 

This should either (1) result in an error ("cannot prepend a string to a SeqLike"), or (2) yield a new SeqLike with the metadata of the parent SeqLike and the correctly-ordered sequence string ("ACTGTTTT"). If the latter, should this new SeqLike be renumbered?

As discussed in #49. The naming here is because we inherit from the collections abstract base class of Sequence. Seqs are coupled to Biopython, Sequences are more generic.

Needs new tests and some documentation, just want to pause here for general feedback and workflow testing.

All current tests pass.

edit for usage:

In [1]: from seqlike.SequenceLike import SequenceLike

In [2]: s = SequenceLike(["abc", "abc", "qwe", "asd"])

In [3]: s
Out[3]: ['abc', 'abc', 'qwe', 'asd']

In [4]: s.to_index()
Out[4]: array([0., 0., 2., 1.])

In [5]: s.to_onehot()
Out[5]: 
array([[1., 0., 0.],
       [1., 0., 0.],
       [0., 0., 1.],
       [0., 1., 0.]])

In [6]: s.alphabet
Out[6]: ['abc', 'asd', 'qwe']

In [9]: s2 = SequenceLike([0, 0, 2, 1], alphabet=["abc", "asd", "qwe"], encoding='index')

In [10]: s2
Out[10]: array(['abc', 'abc', 'qwe', 'asd'], dtype=object)

Also a quick note: when dispatching for SeqLike, I guess we can use SeqLikeType.__args__ as the type instead of hardcoding a list of potential types, otherwise it's hard to keep them consistent (since I saw there's a plan to include torch.tensor)

Originally posted by @pagpires in https://github.com/modernatx/seqlike/issues/41#issuecomment-1000942209

Removed **kwargs from SeqLike aa, and hard-coded arguments in the call to self.translate so as to return original attributes like id and name. This is the behavior I would expect from aa().

I left the interface to SeqLike translate as is to maintain flexibility.

Closes #25

I had a couple import errors (see below) when testing this package out that are related to undeclared dependencies in setup.py. Some were in requirements.txt, but not declared in setup.py's install_requires, so I moved them all to setup.py and removed requirements.txt (so this doesn't happen again :crossed_fingers:).

I also added 2 extras:

• seqlike[test]: this has the pytest* deps (which I removed from the main install_requires)
• seqlike[notebook]: adds bokeh, which is (currently) required when using ipython/jupyter (I think there were some refactors to make bokeh a lazy import, but it is still imported eagerly/required in notebooks)

Here are the errors I saw:

.../lib/python3.7/site-packages/seqlike/SeqLike.py in <module>
      9 from typing import Callable, Optional, Union
     10
---> 11 import lazy_loader as lazy
     12
     13 from Bio.Seq import Seq

ModuleNotFoundError: No module named 'lazy_loader'

.../lib/python3.7/site-packages/seqlike/draw_utils.py in <module>
     19 try:
     20     get_ipython
---> 21     from bokeh.io import output_notebook
     22
     23     output_notebook()

ModuleNotFoundError: No module named 'bokeh'

The wrapper function muscle_alignment permits alignment using Muscle 3.8:

from seqlike.alignment_commands import muscle_alignment

sequences.seq.align(aligner=muscle_alignment, muscle_arg1=something, muscle_arg2=something)

and addresses #24. All tests pass except tests/test_assets.py::test_free_mono_font_exists. Two notes:

1. the latest version of Muscle is v5.1 and has a different interface than v3.8; MuscleCommandline is compatible with v3.8.
2. the preserve_order parameter (preserves original sequence order, as aligner may try to group sequences by similarity) may still be buggy.

It'd be nice to support arbitrary alphabets for sequences that are not necessarily string-type. For e.g. we may want to do sequence of codons, or sequence of other entities.

Doing so would allow us to access the to_onehot() or to_index() capabilities of SeqLike objects without necessarily being bound to BioPython SeqRecord/Seq objects.

Potential challenges:

1. We would break the "default to SeqRecords pair" that we assume in SeqLike. A list of codons is neither!
2. We may need to rearchitect the SeqLike object such that there is a .sequence and .alphabet, which the encoder functions expect (?). .to_*() functions.
3. We may need a more generic SeqLike object from which our current SeqLikes inherit.

A good concrete first step here is to create an Abstract Base Class for discussion purposes.

I've encountered the situation where we need to represent mutations of a sequence. Having written essentially the same code over and over, I thought it might be good to talk about some of these ideas here.

I was thinking of something along the lines of two classes: a Mutation and a MutationSet. Defining them as classes allows for certain semantics:

s = SeqLike('MKAIL')
mut = Mutation('A', 2, 'C')
# mut's repr would look like A3C
mut2 = Mutation('K', 1, 'R')

We may consider 'addition' to be an application of mutations to a reference sequence:

s2 = s + mut
# s2 <-- SeqLike('MKCIL')

s2 = s + mut2
# s2 <-- SeqLike('MRAIL')

Mutations can also be offset by position:

mut3 = mut - 1
# mut3 <-- A1C, but would raise an error if wt sequence does not have A at index 1.

If we need to hold multiple mutations together, we might use a MutationSet:

mutations = MutationSet([mut, mut2])
# mutations <-- [K1R, A2C] (automatically sorted by position, then by letter)
s2 = s + mutations
# s2 <-- SeqLike('MRCIL')

MutationSets could also be offset by position:

mutations + 1  # would give us [K2R, A3C]

We could also consider subtraction of two SeqLikes to give us the 'diff' as a MutationSet:

mutset = s2 - s1  # (left is canonically considered the 'wt')
# mutset <-- [R1K, C2A]
mutset = s1 - s2 
# mutset <-- [K1R, A2C]

I'm not sure what other semantics I might have missed here. Any thoughts?

UPDATE: I changed the position numbers above to reflect Python indexing rules, not canonical indexing. I am sure we could magically handle both, but IMO because SeqLike uses Python indexing rules, Mutations and MutationSets should also use Python indexing rules for positions.

I have a need to adjust alignment parameters; for example, I have encountered something akin to this issue, and the proposed solution from the author of MAFFT is to adjust one of the MAFFT parameters.

Adjusting alignment parameters via the .seq.align() API might be helpful. A few designs for the user-facing API that I can think of include:

# default aligner is MAFFT, so we can pass through the command line options via kwargs.
sequences.seq.align(ep=1.59, op=0.0)

# want to use MUSCLE instead of MAFFT
from seqlike.AlignCommandLine import MuscleCommandLine as muscle
sequences.seq.align(aligner=muscle, muscle_arg1=something, muscle_arg2=something)

We should ensure that all contributions are recognized. Following the all contributors spec, we should use the bot to help recognize contributors of all kinds to the project.

https://allcontributors.org/