Mask gap characters

[jameshadfield]

This extends augur mask to be able to
(a) mask all gaps in all sequences
(b) mask all gaps in sequences which have gaps but no Ns.

A quick scan of 100k SARS-CoV-2 samples (samples over time & geography)
reveals that around a third have gap characters but no Ns. This
indicates bioinformatics pipelines using gap characters to represent
missing data are widespread. Option (b) above is intended to fix this.

The current implementation results in no noticeable slowdown. Testing on
the above 100k data set with our common nCoV pipeline masking parameters
takes ~13s. Masking gap characters doesn't change this.

Closes #1043

Note that our nCoV pipeline uses a script instead of augur mask. Looking at this I think the only extra functionality is the ability to mask terminal gaps. This would be trivial to add to augur mask (but if we are going to mask all gaps in all sequences then we don't need it). I would suggest adding this as another flag --mask-terminal-gaps as it doesn't fit well with --mask-gaps in this PR. An alternative would be to change this PR to use --mask-all-gaps + --mask-gaps-if-no-Ns flags. 🤔

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 61.81%. Comparing base (862bbb6) to head (a14d495).
Report is 1221 commits behind head on master.

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1048      +/-   ##
==========================================
+ Coverage   61.68%   61.81%   +0.13%     
==========================================
  Files          52       52              
  Lines        6287     6309      +22     
  Branches     1583     1587       +4     
==========================================
+ Hits         3878     3900      +22     
  Misses       2138     2138              
  Partials      271      271              

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[jameshadfield]

Revisiting this today, I think there are two different ways to expose the (command line) arguments:

1. --mask-gaps {all, if-only-gaps, terminal} and allow >1 values, e.g. --mask-gaps if-only-gaps terminal
2. Three flags: --mask-all-gaps, --mask-gaps-if-no-ns, --mask-terminal-gaps

I prefer (1). Thoughts?

[tsibley]

--mask-gaps if-only-gaps terminal to me reads as the terminal option conditioned by if-only-gaps or something. It doesn't parse to me as the union of their effects (e.g. always mask terminal gaps, and mask non-terminal gaps if the sequence doesn't contain any Ns).

I'd suggest (2) (or similarly worded options) as the more conventional approach here.

--mask-gaps-if-no-ns → --mask-all-gaps-if-no-Ns?

[tsibley]

Maybe the verbose but very clear --mask-all-gaps-in-sequences-without-Ns?

[emmahodcroft]

In theory I like 1 (simpler), but I agree with @tsibley that it's not necessarily easy to understand how the arguments interplay. For me I think the options are clear enough - but --mask-all-gaps-if-no-Ns is more confusing for me (I think because it references the 'all gaps' option) - I'd say it needs to be either the short original or the very verbose third option.

[jameshadfield]

Update: I've implemented the ability to mask terminal gaps by reusing the code from https://github.com/nextstrain/ncov/blob/master/scripts/mask-alignment.py & added tests. This will allow our ncov pipeline to use augur mask again.

The argument interface remains as (1) from above, however I'm more than happy to change to (2).

[emmahodcroft]

@rneher points out that:

augur ancestral treats terminal gaps as missing data by default.

So do we need the script that we currently use in ncov to mask terminal gaps (and/or the propose new flag)?

Also I felt like I understood this better earlier this week but now I've got myself a bit mixed up - the --mask-gaps-if-no-ns flag would turn a column of gaps to Ns only if it consists only of gaps? (So not if it has Ns or bases?) Or only not if it has Ns?
If the latter, is there a risk here that we fill in true but well-defined deletions (where nobody is calling them as N)? Or is that near-impossible?

[rneher]

note that augur align has this functionality.

augur align -h 
[...]
 --fill-gaps           If gaps represent missing data rather than true indels, replace by N after aligning. (default: False)

so the question is whether this should be in mask rather than align?

[emmahodcroft]

Sorry, I was behind on Slack but just caught up and I think I have been interpreting this wrongly: I assumed this was a column-based approach but now I think it might be a 'per sequence' approach. I'm a little less sure about this - it seems like outside of nCoV, how confident would one be that a sequence that uses only gaps is 'wrong'? But on the other hand, not sure I see another good solution. However, may raise the question of whether this should be a script for nCoV or something more widely useable in augur?

[jameshadfield]

Revisiting this with fresh eyes I'm going to close this PR -- the introduced functionality is either available in existing augur commands (augur align --fill-gaps) or is pathogen-specific and can be implemented in script. I've instead documented our current behaviour in this PR.

I though about re-writing this PR to only implement augur mask --mask-gaps, which would do largely the same thing as augur align --fill-gaps but would be useful in the cases where we are calling the aligner directly (e.g. using nextalign). Happy to do this if there's interest.