Panelapp Australia upgrading Gencode/Ensembl version

[drtconway]

Hi Seqr-people,

We're working on a major upgrade for Panelapp Australia to migrate the underlying gene annotation from Gencode 27 (Ensembl 90) to Gencode 49 (Ensembl 115). This entails quite a few genes changing symbol (name), and Ensembl ID, and in some cases both.

This will lead to some problems - some fairly direct, and some more subtle.

I'm still investigating, and will add more details as I discover them, but let me start with what we have determined so far.

Where the blob of JSON from the Panelapp API is parsed, the hard coded "90" here, will break:

seqr/panelapp/panelapp_utils.py

Line 23 in 7e611fc

gene['gene_data']['ensembl_genes'].get('GRch38', {}).get('90', {}).get('ensembl_id'): gene

But there's a more subtle and pervasive problem, which is that there will be gene symbols and Ensembl IDs that will be inconsistent with GeneInfo. GeneInfo can be updated with the relevant gencode version, I'm still trying to elucidate what that "fixes" and what remains problematic.

One area of concern is in gene_utils.py

seqr/seqr/utils/gene_utils.py

Line 198 in 7e611fc

gene_ids.update({gene_ids[0] for gene_ids in gene_symbols_to_ids.values() if len(gene_ids)})

where gene symbols are resolved into gene IDs for query resolution. The code currently finds the most "recent" match only, which means existing projects/samples stored in the backend will not match the query without being reloaded.

We thought we should give a heads-up, since this will affect any Seqr instance using Panelapp Australia for gene lists. I understand that Genomics England are also working toward an update of their underlying genome annotation, but I can't speak to any detail about that.

We would be keen to work with you to make the transition as smooth as possible for everyone!

Tom, for the Bioinformatics Methods Group at MCRI.

[ssadedin]

@hanars - just paging you to alert to this.

We are proposing to implement the above change on the PanelApp API this coming Friday, which will stop panels updating from PanelApp Australia until at least a basic patch is implemented along the lines above.

We will propose a patch later today to address the problem, along with a couple of other changes that improve the safety when gene symbols have multiple definitions across different versions of Ensembl / Gencode.

Apologies for the short timeline for this change!

[hanars]

Hi all, really appreciate the proactive conversation and support you all are providing for this change! To address some of the specific issues raised above:

The code currently finds the most "recent" match only, which means existing projects/samples stored in the backend will not match the query without being reloaded.

Luckily this is not how seqr works anymore. While historically we would load new projects with newer annotations while old projects would retain older annotations, this is not the case for any supported version of seqr. All data in seqr is now called with the same version of VEP (110) and gene annotations are consistent across all loaded data. In fact, VEP 110 uses gencode v44 which indicates that we should be updating the latest gencode version of seqr regardless of this coming PanelApp change

GeneInfo can be updated with the relevant gencode version, I'm still trying to elucidate what that "fixes" and what remains problematic.

We have bene intentionally slow to update our gencode version because changing gene symbols can potentially break certain ongoing reports we need to provide for grant compliance. As such, we will be unable to update the gencode version with less than a weeks notice, as we need to take the time to quantify and assess the impacts to these reports. While going forward we will update the version, for the purposes of assessing the impact of this coming panel app change we need to assess it relative to the current available gencode versions in seqr.

We are proposing to implement the above change on the PanelApp API this coming Friday

Currently, we update PanelApp over the weekend so I think our plan will be to do a PanelApp update on Thursday, disable our usual weekend update, and then reenable it sometime next week when we have had time to test out and release the patch. Let me know if that makes sense to you, or if you would recommend a different approach?

Thanks again for all your work on helping support our team through this transition!

[drtconway]

We completely appreciate the complex issues relating to upgrading the annotation! Gene symbols (names) change, Ensembl IDs change, new genes appear, or change from non-coding to coding, and so on. The upgrade in Panelapp Australia is from Gencode 27 (Ensembl 90 - circa 2018), to Gencode 49 - the most recent, and we have be quite paranoid about making sure we don't make errors that result in genes being dropped, etc. All of which is to say, we're on board with being cautious! In the case of Panelapp Australia, one of the important drivers was the fact that there are a steady trickle of genes being recognised with important disease association which are not in the old Gencode, or are incorrectly annotated there, so upgrading the annotation has become imperative to keeping gene lists accurate and up to date with the state of our knowledge.

After some analysis, here are 3 small changes which we believe will address our issues. From what you write, you may not want to incorporate them into the main source tree, but we offer them at least as insight in to what we're doing. All three are independent of the back-end (ES, Hail, CH), unless I am mistaken, though with the new CH back-end (which I have only had the most cursory look at myself) other strategies to address the issues may be apparent.

1. Parsing the Panelapp Australia JSON correctly

The following change just makes the parsing of the JSON find the most recent version, rather than being hard coded to '90'. I noticed the original code in our repo was formatted differently (with line breaks), but is otherwise the same as yours.

Without this change (or a similar one), seqR instances will simply fail to update any panels.

bioinfomethods@0cba0d5

2. Hard-fail rather than drop genes from panels

Currently, the code to update panels emits a log warning if any genes in a panel do not have an Ensembl ID present in GeneInfo, and drops the gene from the gene list. (There's also a latent issue that we have not tried to address here: it doesn't check if the gene symbol in the panel is consistent with what is in GeneInfo, which is potentially a problem.)

Because it is a log message in the Panelapp Australia update code, and not at all user-visible, we are not comfortable with dropping genes, so this change turns the warning into an exception which will cause the panel update to fail outright, which we consider preferable.

bioinfomethods@3e0c91f

3. Search for all Ensembl IDs, not just the latest one.

We will need to support multiple annotation versions across projects at least in the short term in our instance. (We are still using the ES backend, but will hopefully be upgrading to the CH backend soon!)

Further, we noticed in our GeneInfo table there are some instances where more than one Ensembl ID is present for a named gene within a Gencode version. The Gencode folks have gradually been eliminating these cases, and converging on names mapping to a single Ensembl ID, but even the latest Gencode still has instances where a name points to more than one Ensembl ID. I'm not absolutely certain, but I think the ordering of these duplicates, when sorted by Gencode release is not guaranteed to be stable, in which case the behaviour is not deterministic, but even if I am wrong about that, and it is deterministic, the ambiguity is obviously concerning.

The following change addresses both problems by populating the set of Ensembl IDs to search with all the results rather than just the first. In the case of multiple annotation versions, at most one of the alternatives will match for a given family, though in the ambiguous case you may get both (or all) IDs in results. We consider a few extra results preferable to missing results.

This change is also defensive in the presence of gene symbol changes in the sense that there are instances we see in the upgrade from Gencode 27 to Gencode 49 where a symbol "moves" from one gene to another - so assuming both the old and new annotations are in the GeneInfo table, you would see results from both.

bioinfomethods@a4d954f

[hanars]

Hi Tom,

This all makes sense. We would definitely want to incorporate change 1 if you are open to putting a pull request up against the main branch.

For change 2, I would be concerned that by hard failing we would lose the ability to update PanelApp in the main seqr instance entirely, as it seems almost certain to me that there will be genes in panel app with no corresponding ensembl ID in seqr. We would much rather ingest all the PanelApp data that corresponds to genes in seqr and skips information for unsupported genes, rather than failing the updates for all the genes we do have just because there are also genes we do not have. Does that make sense? Or is there a piece to this that I am missing?

For change 3, since multiple annotation versions are not an issue in the main seqr anymore theres no need for us to incorporate this update in the main branch. I would recommend that once you update to the latest seqr backend you deprecate this change as well, as it will have negative performance implication for Clickhouse searches.

Best
Hana

[hanars]

For the duplication you are managing in change 3, I do want to flag that the majority of cases where there are multiple ensembl IDs associated with a gene symbol within a gencode release are cases where different genome reference builds have different ensembl IDs (i.e. one ensembl ID only has data for hg19 and one only has hg38). In build 39 there are over 400 genes that have duplicate ids for a given symbol, but only ~50 of those are duplicated across both supported genome builds. Since you note that gencode seems to be trying to reduce the frequency of this issue in later releases, it might be with doing a quick check to see if different refernece builds accounts for all the remaining cases in the latest release, or if there are still some of these left for real.

For the cases where its not a build specific issue (assuming there are any of those left in the latest release), I am also not convinced that choosing all the genes is necessarily accurate. In gencode 39 for example MKKS (a gene associated with several panels) has 3 corresponding ensembl IDs, one of which is 1.7kb and the other two of which are < 0.2kb. To me, its not clear that PanelApp would really be making assertions about all of these genes. To be clear, arbitrarily picking the first time we see a gene symbol the way we do now is nondeterministic and not necessarily going to pick the correct ensembl ID either, so I am not opposed to the changes you are making, but its something to keep in mind

[drtconway]

With respect to the duplication issue, if I understand the code correctly, while Panelapp Australia panels contain Ensembl IDs and most certainly takes care to point to the correct Ensembl ID where symbols are non-unique, when the search form is populated in the UI, we basically round-trip through gene symbols. This means rather than respecting the Ensembl ID in the panel, you get whatever the query on GeneInfo gives back. This means already, you don't necessarily end up querying on the Ensembl ID that was in panelapp.

I agree that having extra Ensembl IDs in the query in the back end has a performance cost whatever the backend, but we are also very concerned about false negatives - where a variant may exist in a sample and the search quietly fails to return it.

In fact this is the motivation for the second change - the hard fail. We would rather fail to update a panel, over quietly (from the user's point of view) dropping genes from panels. The user is much more likely to notice that (e.g.) Mendeliome hasn't updates for a surprising amount of time, than noticing that that a few of the 5000ish genes are missing.

The possibility of false negatives keeps us awake at night!

[hanars]

Regarding the duplication issue, I totally agree that the current behavior of arbitrarily picking one ensembl ID is wrong! I was just suggesting that rather than choosing all possible ensembl IDs for a given gene symbol in PanelApp, which I think may be inaccurately adding certain genes to gene lists they do not belong in, you may want to investigate if theres a different deterministic way to disambiguate these genes. Since its a very small number of impacted genes in v39 I would hope its an even smaller list in v49, which means doing an upfront manual review of these duplicates is actually feasible and its possible theres a clear pattern in these genes that could be used to deterministically choose the correct ensembl ID(s) going forward. It may even be worth asking a PanelApp curator their thoughts on a specific example if there is one. The example I gave above is v39 specific and I think fixed in v49, but in that case for instance I would assume that the small secondary genes have never been evaluated for gene/disease association and that all the PanelApp curations for that gene refer only to the single lager gene - in that case surfacing to users that the smaller genes are specifically associated with a disease would be false. I agree that false negatives are an important concern, but false positives with misleading evidence presented to the user may not be a better solution. Again, this may be moot in v49 depending on what duplicate genes actually are left once you control for genome build, its possible that there are no examples like the one I found in v39, but it feels worth investigating.

Regarding the hard failure, I agree that that behavior makes sense once seqr has a gencode version that corresponds to the PanelApp gencode version, as in that scenario there is no reason why you would expect to get a gene from PanelApp with no corresponding gene in seqr and would want to investigate that unforeseen circumstance. However, as long seqr is several gencode releases behind PanelApp it is a certainty that there will be genes in some panels absent from seqr. In the scenario where 100 new genes are added to Mendeliome and 97 of them have a corresponding v39 gene and 3 do not it would be a far better user experience that we surface this new information for 97 genes then that we fully disable the update (especially since there can be no variants annotated in those 3 missing v49 genes in the first place)

[drtconway]

Well, Panelapp Australia has migrated, so if people do nothing, then panels will stop updating due to the "90" issue.

With respect to the duplication issue, I note that our general position is that we would much rather include a few extra (spurious) variants in search results, rather than miss variants, which appears different to your desire to avoid spurious variants. Both are defensible in context!

One thing that I don't understand, which maybe you could spell out a bit for me is about the CH backend. I know in the ES backend, the gene ID is in the ES document for each variant, and that is a key for normal searches with gene lists. Your response earlier sounded like the CH backend behaves differently, so changing the gene annotation wouldn't require reloading samples. Did I misunderstand? Can you explain that?

Thanks Tom.

[hanars]

Regarding the "90" issue, I understood that someone from your team would be putting up a Pull Request against the main seqr repo to handle these breaking changes sometime this week. If that is not the case please let me know so we can try and allocate some bandwidth toward this on our end.

Regarding the CH backend, theres a known performance hit in this backend as gene lists get longer (so a search in 5000 genes runs a bit faster than one in 5500 genes). So its not a huge deal, but including hundreds of extra genes in a search that ultimately are incapable of returning any hits because they correspond to gene IDs that have no loaded data in the backend is going to cause a performance hit with no actual benefit to search. In the CH backend gene annotations are consistent across all projects, so theres no advantage to including older gene IDs in searches. The data model is fundamentally different in this backend, so variant-level information such as VEP is stored once per variant in a relational model to the sample-specific data for a variant (gentoype, GQ, etc), so regardless of when a sample is loaded it will have the same annotations as all other data

[drtconway]

That's great news, re the annotaitons in the CH back end - thanks for clarifying.

I've now tested the changes in our local source tree (which has some divergence from yours), so I can now prepare a PR for that change. We just need to make sure that people running seqR instances understand they need to look at the update logs when they update their panels, and may need to update their Gencode release.

Tom.