--fill-gaps eventually adds a furin cleavage site to everyone

[babarlelephant]

Hi, in augur align you are calling --fill-gaps which is just doing

 _seq = _seq.replace('-', 'N')

so that in the nextstrain tree everyone gets a FCS.

I don't have access to filtered_H5Nx_HA.fasta and metadata-with-clade_H5Nx_HA.tsv so I can't make a proposal of modification of the Snakemake with a python script in order to leave the FCS indels but if I did I would be happy to make some tries.

The only experiment I did was to download a few thousands of H5Nx sequences from genbank, to label with LPAI those with a deletion at the FCS, then to remove from the alignment the FCS region, to run mega, and I got that the LPAI sequences where quite clustered together.

[jameshadfield]

I think this is due to the default behavior of augur translate being to infer ambiguous nucleotides, and could be remedied by adding --keep-ambiguous to the ancestral rule.

[lmoncla]

Thanks for flagging this! This is definitely an error caused by the combination of --fill-gaps in augur align and not having --keep-ambiguous in augur translate. This has been fixed with this commit, which does a few things:

1. Removes --fill-gaps in augur align
2. Adds --keep-ambiguous in augur translate
3. Adds a new rule, rule cleavage-site, which reads in the HA alignment, translates it, and infers whether the 4 amino acids immediately preceding HA2 encode a furin cleavage motif. Here, I've defined a furin cleavage motif as the sequence R-X-K/R-R, where X is any amino acid. This results in 2 new colorby options: furin cleavage motif, which can be present or absent, and cleavage site sequence, which displays the sequence of the 4 amino acids immediately preceding HA2.

As expected, we now display that some high-path H5N1 and H5Nx lineages have furin cleavage sites, but that the vast majority of H9N2 and H7N9 sequences do not.