feat(mask): Add --mask-gaps {all,terminals} option
#1286
Conversation
This allows`augur mask` to allow masking of gaps at terminals via `--mask-gaps terminals` or all gaps via `--mask-gaps all`. Masking all gaps symmetrizes treatments of gaps and insertions. It is useful when gaps are jumpy due to sequencing artefacts and allows gaps to be entirely ignored Further discussion: https://bedfordlab.slack.com/archives/C01LCTT7JNN/p1692443660080129
Codecov ReportPatch coverage:
Additional details and impacted files@@ Coverage Diff @@
## master #1286 +/- ##
==========================================
+ Coverage 69.74% 69.77% +0.03%
==========================================
Files 67 67
Lines 7146 7154 +8
Branches 1742 1745 +3
==========================================
+ Hits 4984 4992 +8
Misses 1855 1855
Partials 307 307
☔ View full report in Codecov by Sentry. |
|
See #1048 for some good discussion around this. In the end I didn't implement it because If this PR goes in can you please update https://docs.nextstrain.org/en/latest/guides/bioinformatics/missing-sequence-data.html? |
|
Thanks @jameshadfield! I missed the discussion back then. It seems like my use case (gaps often wrong, sometimes better to strip them) wasn't discussed back then. I don't think I've ever used augur align. It seems like masking gaps belongs most cleanly into mask. I'll check out your prior art! |
|
Nice tests there! I've pretty much reinvented the wheel only worse 😅 |
| @@ -78,7 +78,7 @@ def mask_vcf(mask_sites, in_file, out_file, cleanup=True): | |||
| pass | |||
|
|
|||
|
|
|||
| def mask_sequence(sequence, mask_sites, mask_from_beginning, mask_from_end, mask_invalid): | |||
| def mask_sequence(sequence, mask_sites, mask_from_beginning, mask_from_end, mask_invalid, mask_gaps): | |||
There was a problem hiding this comment.
| def mask_sequence(sequence, mask_sites, mask_from_beginning, mask_from_end, mask_invalid, mask_gaps): | |
| def mask_sequence(sequence, mask_sites, mask_from_beginning, mask_from_end, mask_invalid, mask_gaps=None): |
To match the default noted in docstring.
| if beginning + end > sequence_length: | ||
| beginning, end = sequence_length, 0 | ||
|
|
||
| seq = str(sequence.seq)[beginning:-end or None] | ||
|
|
||
| if mask_invalid: | ||
| seq = "".join(nuc if nuc in VALID_NUCLEOTIDES else "N" for nuc in seq) | ||
|
|
||
| if mask_gaps == "all": | ||
| seq = "".join(nuc if nuc != "-" else "N" for nuc in seq) |
There was a problem hiding this comment.
(non-blocking)
| seq = "".join(nuc if nuc != "-" else "N" for nuc in seq) | |
| seq = "".join("N" if nuc == "-" else nuc for nuc in seq) |
This reads slightly better to me - it pairs the mask character with the gap character.
This allows
augur maskto allow masking of gaps at terminals via--mask-gaps terminalsor all gaps via--mask-gaps all.Masking all gaps symmetrizes treatments of gaps and insertions.
It is useful when gaps are jumpy due to sequencing artefacts and allows
gaps to be entirely ignored
Further discussion: https://bedfordlab.slack.com/archives/C01LCTT7JNN/p1692443660080129
Testing
I tested all options in a workflow. As this is a new option that is a noop unless the option is explicitly set, it shouldn't break any existing code even if the implementation happens to have a bug.
If this PR gets approved, we might want to add a simple cram test, but it'd be a waste of time to do this now before feedback on whether this PR is at all welcome.
Checklist