summary: further improvements to protein handling in sourmash

[ctb]

This is an update of and replacement for #999, which raised a lot of issues around how we were doing protein k-mer calculations.

This issue is being updated after the release of sourmash 4.1.

---

Over the past year, several of the issues in #999 were resolved by the release of sourmash v4, which introduced sourmash sketch (via #1159)

Taking from @bluegenes excellent summary, here are the remaining unresolved issues from #999.

• Translated k=33 should be comparable with protein k=11 (equivalent ksizes). Differentiating by hash functions per can we just use the hash function to flag incompatible signatures, instead of DNA/protein/etc? #751 could facilitate this. In Change --protein k-size meaning  #574, @luizirber pointed out that can we just use the hash function to flag incompatible signatures, instead of DNA/protein/etc? #751 is already enabled in the rust code?
• Potential translation improvements:
  - translate whole sequence, then emit kmers (with --protein, translate whole sequence and then emit k-mers #664)
  - pad translated seqs (Pad non-multiples-of-three nucleotide sequences for protein translation and deal with codon degeneracy elegantly #659)
  - potentially add 3-frame translation (Add option for strand-specific / 3-frame translation for --protein #657)
• Potentially enhance input checks for add_sequence and add_protein? From Document and/or improve protein MinHash API. #720: "in the provided example in Setting "scaled=1" doesn't retrieve all protein k-mers #701, I'm not sure why add_sequence doesn't complain when it's given a bunch of non-ACTG characters."

I do think differentiating sketches by hash functions #751 is its whole own thing and not specifically protein-esque.

---

Notes and thoughts:

It would be nice to figure out if #1037, which checks the first 100bp of FASTA files, is a good approach. Thoughts from #999 (comment) that seem relevant -

• I think we do need both command line and API level checking. The command line can make use of additional info (filename, aggregated across sequences, etc) while the API has to do the trickier job of working with only the sequence it's given.
• I am leaning towards add_dna_sequence and add_protein_sequence at the API level;
• it's not 100% clear how robust it will be to check that any given k-mer is DNA vs prot;
• one strategy might be to look at what fraction of k-mers are valid alphabet;

#1277 changed the Python layer so that ksize for protein was "correct" (the actual length of the word, not k*3!). This still needs to be changed at the Rust layer, though, which would involve changing the JSON signature formats and version.

Also see "Next steps for sourmash sketch" #1169.

[ctb]

here are some notes I took while working through #1159 - copied over from #999 (comment)

Changes to signature computation

Suggested changes to signature computation:

• translate whole read per with --protein, translate whole sequence and then emit k-mers #664
• support top strand/bottom strand only for translate RNA input Add option for strand-specific / 3-frame translation for --protein #657

Signature JSON format changes

What signature format changes should we do in tandem? see #268 for rollup issue

• ksize change for protein Change --protein k-size meaning  #574
• move hash function into minhash section?
• rename 'signatures' list to 'minhashes'?
• support actual command-line concatenation of signatures? investigate modifications to JSON parsing to allow cat to work directly #1093
• could add bp and input file list per add number of bp & sequences => signature file #246, Add documentation to sourmash signature describe for source file origin #769
• ould also add information about preprocessing per Provide a standard format / provenance info for commands used in data preprocessing #269 on a per-command basis; OR, just add information e.g. input type was genome, input type was reads, etc. based on command executed.

Changes to MinHash API

general issue here, #338
more here, #720
and more here, #885, although that is mostly about docs and tests now.

• consider lifting things out to signature, too, per Better conceptual definition on minhash/signature/query #616
• switch to using moltype enums in creation update MinHash attributes/properties for consistency, and/or switch to using moltype attribute #1136
• add_sequence changes per add_sequence method for MinHash #186? add_sequence(<str>, moltype=...)
  • moltype should be enum
  • moltype could be assumed to be DNA for DNA minhashes, aa|prot|protein for protein minhashse.
  • could have it be dna|rna_top|rna_bottom for protein minhashes to translate incoming?
• do we want to make different MinHash classes for different moltypes?
  • protein MinHashes are more complicated than DNA..
  • the default ksize/scaled/etc could mirror the CLI defaults