Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

fix: translate - strand annotation #966

Closed
wants to merge 3 commits into from
Closed

fix: translate - strand annotation #966

Changes from 1 commit
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
5afabe0
fix: translate `-` strand annotation now correct
corneliusroemer Jun 7, 2022
64e12a1
fix: translate strand - make `+` default, only `-` iff `-1`
corneliusroemer Jun 7, 2022
ea4abc1
chore: add fix to CHANGES.md
corneliusroemer Jun 7, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion augur/translate.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -382,7 +382,7 @@ def run(args):
'type':feat.type,
'start':int(feat.location.start)+1,
'end':int(feat.location.end),
'strand': '+' if feat.location.strand else '-'}
'strand': '+' if feat.location.strand == 1 else '-'}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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}

Copy link
Member

@jameshadfield jameshadfield Sep 15, 2022 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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 🐛

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more. 

'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.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Member

@tsibley tsibley Aug 18, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

. 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".

if is_vcf: #need to add our own nuc
annotations['nuc'] = {'seqid':args.reference_sequence,
'type':feat.type,
Expand Down