Fix SnpEff exonic property

[wm75]

A Gemini user on usegalaxy.eu reported yesterday that gemini's is_exonic column is wrongly populated for SnpEff-annotated composite effects. I think this simple typo here is the reason.

[wm75]

In addition, I think the consequence property of OldSnpEff may be wrong, too. The old EFF annotations are using +, not &, right?

[wm75]

Ah, OldSnpEff does not even have a self.effects['Transcript_BioType'] (hence the now failing test). Also, SnpEff.coding seems to suffer from the same consequence vs consequences mix-up.

I could fix those, but maybe tell me if I'm completely wrong with my assumptions first.

[brentp]

yes. looks like you're right. not sure how this has escaped noticed for so long. that'd be great if you updated to include OldSnpEff.

[wm75]

Ok, so this was a bit more work than I intended originally, but I cleaned up the class hierarchy of Effect and its subclasses quite a bit and fixed a number of bugs, mostly around the detection of exonic and coding effects of variants:

• the Effect base class now defines exonic, coding, splicing and lof methods, and their is* counterparts. Most subclasses simply use the inherited methods, only OldSnpEff overwrites most of them. This makes it easier to ensure a consistent meaning of these methods where the is* methods now always return the answer for the top consequence, while the others scan all consequences.
• in particular, coding and is_coding profit from this because they have been defined really inconsistently before for the different subclasses (and some of the previous definitions shouldn't have worked properly).
• OldSnpEff now tries to use so-mapped effects in all those methods falling back to assuming new-style effects only if that fails. (I don't think that coding, exonic and their is_* counterparts did work before for OldSnpEff
• to deal with bcftools output, the truncated (without variant) terms in the exonic impacts set are now also lower-cased, just like they are for the IMPACTS_SEVERITY list. Please check if this is the correct thing to do since I haven't looked at any bcftools output, but just relied on the lowercasing for IMPACTS_SEVERITY to be the correct thing. -I've added a couple of additional effects to EXONIC_IMPACTS and added a new CODING_IMPACTS set that bundles the logic of what is considered coding in one place. Note that I chose not to consider gained start codons coding. I did so by "majority vote" since this exact behavior differed between Effect subclasses before - I think the BCFT subclass considered them coding. If you disagree on this choice just change the behavior as you see fit.
• I've also made the code for the severity property a bit clearer (I hope you agree)
• I've moved the consequences property definition to the base class making its exact behavior (especially the caching) uniform across the subclasses.
• Finally, I've also changed the caching mechanism in top_consequence to reflect that of consequences.

[wm75]

The thing is though that after all this cleanup work it becomes apparent that most of it won't help Gemini because from https://github.com/arq5x/gemini/blob/5db2e52ae4aee413a8780df538565c90a38e4b11/gemini/gemini_load_chunk.py#L498 onwards Gemini only uses the top impact effect and its top consequence to populate is_exonic, is_coding, is_lof and is_splicing. At least for is_exonic and is_coding this is clearly against what the documentation of the Gemini database scheme says. I'm just not sure which behavior, the documented or the implemented one, is the desired one.

If you want I can raise this as a separate issue in the Gemini repo.

[brentp]

hi, thanks for looking into this.

re the gemini note. I think it's doing the right thing. the top_impact is stored in the variants table. and it will have the correct is_exonic, etc for that chosen impact. all other impacts are stored in the variant_impacts table with is_exonic et al set according to each impact.

[wm75]

yeah, I guess you can argue like that, which would mean that the docs need an update.

[wm75]

The counter argument would be one of practicality: what's documented - is_exonic = 1 if any effect on any transcript is an exonic one - is just very hard to derive from the variant_impacts table. Haven't actually verified that it's possible at all.