search

dib-lab / khmer

In-memory nucleotide sequence k-mer counting, filtering, graph traversal and more
http://khmer.readthedocs.io/
Other
754 stars 296 forks source link

unique-kmers.py ignores the entirety of any sequence with non-ACGTN in it #1540

Open ctb opened 7 years ago

ctb commented 7 years ago

e.g. if a chromosome from a genome sequence has a single non-ACGTN in it, the entire sequence will be ignored. This is unexpected.

This should probably be changed so that such k-mers can be skipped over, or optionally raise an exception - e.g. in the sbt_search branch of sourmash compute, I now convert non-ACGT bases into 'N', unless --check-sequence is specified. The khmer-consistent approach would be to convert them into 'A' (see khmer.utils.clean_input_reads).

ctb commented 7 years ago

Note also that there's a divide by zero if no k-mers are consumed ;).

standage commented 7 years ago

Discussion from back in September is relevant: #1036. Since then, whenever I've been consuming genomic sequences, I just preprocess to split on non-ACGTs and filter the resulting fragments by some minimum length. But it would be nice to have khmer handle this nicely.

standage commented 7 years ago

Also, some genomic sequences have long internal or terminal stretches of Ns. We'd have to consider whether we want to consume all of these as converted As or ignore (and how exactly that would work).

ctb commented 7 years ago

I'm fine with N->A conversion generally; history/experience suggests that most sequence has 'em anyway so no harm done. But could easily trash any collection of Ns that is longer than k-mer size, for example.

Also, a useful snippet of code that I just used in sourmash:

seq = re.sub('[^ACGT]', 'A', seq)
ctb commented 7 years ago

Regarding the comment above by @standage and my code block: I think we should provide a function in khmer/utils.py that does both the split that @standage does, as well as the 're.sub' above.This would give people an arsenal of tools with which to pre-treat their sequence. We could also put this in screed instead of khmer. Thoughts as to which, or if this is even a good idea?

standage commented 7 years ago

See https://github.com/dib-lab/khmer/issues/1541#issuecomment-275559588.

standage commented 7 years ago

I mean, I think eventually we might want all of the approaches implemented at the C++ level, but...tradeoffs.