fix: translate - strand annotation

[corneliusroemer]

Annotations output by augur translate always contained + irrespective
of what the input gff or gb file specified and biopython returned

This bug was down to the assumption that biopython feat.location.strand
returns a boolean, when it in fact returns numeric 1 or -1 for pos/neg strand

In order to be backwards compatible and output + by default, for example when no strand directionality is given, I reversed the test, so it puts - iff strand is -1, otherwise it's +

Testing

The change is minimal, it's a very localized bug fix.

There is almost no testing for augur translate, so this would be good to set up at some point, but it shouldn't block this important fix (it's directly relevant for all MPX builds). Maybe @victorlin can add this to the backlog of Augur work?

This makes Auspice display pos and neg strands correctly as shown below:

[rneher]

there is something wrong with the tests though...

[ZarulHanifah]

Sorry for coming in, just want to ask, does augur translate currently assumes that all CDS in the genbank file are forward strand?

I am working on orf virus, and some of the genes are actually in reverse strand, and I do notice that the amino acids are different that what is present in the genbank translations.

[jameshadfield]

Sorry for coming in, just want to ask, does augur translate currently assumes that all CDS in the genbank file are forward strand?

Hey @ZarulHanifah -- this is purely a visualisation fix (whether the gene was displayed above/below the line in auspice), it doesn't change the actual AA translations at all.

[jameshadfield]

Ouch, let's get this functionality fixed! (The failing test logs have expired, so I can't see what's wrong there.)

I am confused as to how we do this correctly for TB 👇 as running the (VCF) example in ./tests/builds/tb has feat.location.strand of -1 or 1 (and thus all are exported as being on the positive strand).

[tsibley]

The docs for FeatureLocation.strand say it can be -1, 1, 0, or None, described as:

For nucleotide features you will also want to specify the strand, use 1 for the forward (plus) strand, -1 for the reverse (negative) strand, 0 for stranded but strand unknown (? in GFF3), or None for when the strand does not apply (dot in GFF3), e.g. features on proteins.

We should probably account for all of those. One way is something like this:

'strand': {+1:'+', -1:'-', 0:'?', None:None}[feat.location.strand]

We'd then want to also update the annotation schema definition to match the new values.

[tsibley]

Such a fix is a good time to systematically review strand handling and make sure there aren't similar bugs elsewhere. For example, this conditional on strand should also probably account for all values instead of only testing for -1?

augur/augur/translate.py

Lines 167 to 176 in ea4abc1

	if feature.strand == -1:
	aaRepLocs = {(end-start-i-1)//3:safe_translate( str_reverse_comp( "".join([sequences[seqk][key+start]
	if key+start in sequences[seqk].keys() else ref[key+start]
	for key in range(i-i%3,i+3-i%3)]) ))
	for i in genNucSites}
	else:
	aaRepLocs = {i//3:safe_translate( "".join([sequences[seqk][key+start]
	if key+start in sequences[seqk].keys() else ref[key+start]
	for key in range(i-i%3,i+3-i%3)]) )
	for i in genNucSites}

[jameshadfield]

We should probably account for all of those.

💯 Yeah let's use this opportunity to add some tests and move to a style where we ensure the value is one of the expected values, and error if it's not (rather than a "if this then ... else assume this" approach)

However if the value is not '+1' or '-1' then how do we translate it? It seems to me we either assume it's positive strand, and export it as such, or we don't translate the gene and therefore don't export it in the annotation. I believe the current code intends to do the former, although there's a 🐛

[tsibley]

However if the value is not '+1' or '-1' then how do we translate it?

Relatedly, if we do make an assumption for unknown (0/?) or unstranded (None) annotations and translate them, and thus export them, do we pass that information along to Auspice? Or do we keep Auspice unaware? Basically, does strand in the export v2 schema stay only +/- or gain ?/null too?

If we do pass thru from Augur strand values that aren't +1/-1, then we'll need to update at least a few spots in Auspice.

[jameshadfield]

do we pass that information along to Auspice? Or do we keep Auspice unaware?

We should pass the information to auspice, IMO, whilst being aware that our current visualisation is not going to show them any differently than a positive-strand gene (for now). This issue suggests adding a tooltip where we could surface that information.

[victorlin]

'strand': {+1:'+', -1:'-', 0:'?', None:None}[feat.location.strand]

This was implemented in #1211.

I started #1279 to follow up with updating the schema and adding tests.

[jameshadfield]

Following up on this, why did we choose to keep the BioPython representation of None when "the strand does not apply", rather than the GFF representation of "."? We use GFF spec for the other strand values.

[tsibley]

. is GFF3's explicit null value so I suggested Python's None because it translates to JSON's null, which is of course JSON's explicit null value. That is, I chose to preserve semantics of the value over the syntax of it.

I believe GFF3 used . because it wanted an explicit null value to avoid the pitfalls of empty fields in a TSV row, e.g. it wants the clearness of \t.\t instead of \t\t. JSON's null serves a similar purpose.

That all said, I think it would be equally reasonable to choose to preserve the GFF3 syntax so that we can say "strand is as defined in GFF3".

[tsibley]

I am confused as to how we do this correctly for TB point_down

We did, apparently, but don't any more. That TB build hasn't been updated since 2018, and it contains both -1 and +1 strand values:

$ curl -fsSL --compressed https://data.nextstrain.org/tb_global_meta.json | jq -c '.annotations | to_entries | .[] | {key, strand: .value.strand}'
{"key":"gyrB","strand":1}
{"key":"gyrA","strand":1}
{"key":"Rv0010c","strand":-1}
{"key":"Rv0011c","strand":-1}
{"key":"Rv0026","strand":1}
{"key":"Rv0039c","strand":-1}
{"key":"ponA1","strand":1}
{"key":"Rv0147","strand":1}
{"key":"mce1R","strand":-1}
{"key":"lprO","strand":-1}
…

[tsibley]

Looks like 74125f5 is what broke that, first released in Augur 6.0.0, when we moved from +1/-1 in v1 JSONs to +/- in v2.

[victorlin]

The described bug has been addressed in #1211.