Stdpopsim: Associate species with upstream genome assemblies and IDs

[jeromekelleher]

As discussed in the meeting 2020-01-14, currently species genome data is hand-curated in a haphazard manner. It would be much better if we pulled information about chromosomes and their lengths from an upstream resource like UCSC or Ensembl.

The current proposal is that each species should have a list of associated genome assemblies, which have the ID of the species assembly and derive the genome data directly from the upstream resource. Thus, we want to decouple the genome assembly from the species definition, and also make the species class explicitly link up with upstream genomic resources on this species. This will be useful down the line when we want to pull in gene annotations, etc.

Ideally, we would like to have a local development script which will take an (e.g.) NCBI taxon ID as a parameter, and write out a stub species definition file in the catalog. This would include the most recent genome assembly data. I'm not sure that this is feasible with the upstream APIs, but it would be really nice if we could do it.

[grahamgower]

How about including a copy of the assembly report, and extracting chromsome lengths, etc. from that? The reports are nice tab separated files. I didn't see an obvious API for finding the file from an assembly name, so it might remain a manual process to navigate to the NCBI assembly (e.g. for CanFam3.1: https://www.ncbi.nlm.nih.gov/assembly/GCF_000002285.3/), then
click through to the assembly report (ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/002/285/GCF_000002285.3_CanFam3.1/GCF_000002285.3_CanFam3.1_assembly_report.txt). The QC person would then add a checksum for the assembly report file, to be confirmed in a unit test.

[jeromekelleher]

That could be a good option.

[grahamgower]

There is a tab-separated table for mapping species names/ids to assembly accessions in eukaryotes.txt here:
ftp://ftp.ncbi.nlm.nih.gov/genomes/GENOME_REPORTS/

There are tables with more complete info, which also lists the "assembly name" (eg. TAIR10.1) and ftp path for the assemblies. But these files are much bigger. 50 mb for refseq and 160 mb for genbank.
ftp://ftp.ncbi.nlm.nih.gov/genomes/ASSEMBLY_REPORTS/

[jeromekelleher]

Putting this in to the next release milestone for now as it would be nice to cover. We can remove it if we end up with too much to do.

[grahamgower]

Definitely one for after the next release. Be great to get @andrewkern working on it!

[jeromekelleher]

Agreed - I've created the 0.1.3 milestone so we can kick these cans in there.

[ndukler]

Copied over from duplicate issue:

As we add annotations and recombination maps to stdpopsim some will inevitably be built on different assemblies. Currently we have no ability to handle this.

I think that we should implement an assembly object that contains an:

id
genome definition
citation info
Then modify the existing recombination maps to be tied to a specific assembly id. Furthermore, people should be able to list the assemblies and which recombination maps and annotations are tied to each one. Lastly there will have to be some error checking when people start using some of these things together to ensure that they are all from the same assembly. We'll also have defaults so this is backwards compatible.

This will also require a bit of reorganization of the AWS storage and the addition of the genome definitions there.

[ndukler]

I would suggest going with the Ensembl REST API for handling both genome maps and annotations. It's language agnostic and well maintained.

[andrewkern]

we would want this to propagate over to the annotations in #560

[jeromekelleher]

I've been figuring out the Ensembl API and I think I'm making pretty good progress. We can link up with Emsembl to automatically pull down information about chromosomes etc for the current default assembly, and I can see a way to using this to autogenerate some code for us (I don't think we should pull basic stuff like this down from Ensembl on demand - it would be slow and brittle).

Managing multiple assemblies is super painful, though. So, a fundamental question we need to answer is, do we need to know about multiple assemblies, or do we just fix at what Ensembl has for a particular release? I think fixing at a particular release will work for everything except recombination maps. And this is very tedious, because we will want to mix recombination maps like HapMap GRCh37 with annotations from Emsembl which is on GRCh38.

So, what to do? I hate to say it, but maybe it's simplest to lift the recombination maps over to the latest genome build?

[andrewkern]

oh no liftover!!! have you given it a shot to see how much we lose? https://genome.ucsc.edu/cgi-bin/hgLiftOver

[apragsdale]

Using lifted over recombination maps has caused me some serious headaches before, since the more recent genome build will have segments that are inverted or transposed relative to the previous builds. It can unrealistically stretch and compress recombination rates over regions, and you can easily get things like regions with negative recombination rates (would msprime just error out because of invalid recombination map input in that case?). Worth keeping in mind..

I suppose if that's the only option, we might have to do our own liftover. But if there are published resources that re-infer recombination maps for more recent builds, I think those should be preferred.

[ndukler]

I was thinking of lifting over the recombination maps in 100bp bins or so, and not lift over bins that fail certain criteria. Then re-compute the recombination maps at the same resolution they were at on the source genome (e.g. 50kb) and handle the missing regions by either dropping them from the denominator or replacing them with the chromosome-wide average. Would this would avoid many of the headaches you mentioned @apragsdale ?

[ndukler]

The main issue with re-inferred recombination maps is that they may use slightly different data so they would probably have to be called something different so we would essentially just lose some recombination maps over time as we gained others. I suppose that's kind of already happening but if possible the lifted over maps might make things a bit more apples-to-apples between genome versions.

[reedacartwright]

It would probably help development to reach out the the Ensembl REST developers to get them involved.

[jeromekelleher]

Good idea @reedacartwright, I'll reach out to them.

[ndukler]

Here are the mapping files. They have 12 of the more common species so not as many as I though but they have everything we currently use except for Pongo.

[ndukler]

I like the format too. Has all the information we need. I think we're just going to have to make building a table like this part of adding a species. Don't think we can really automate it.

[reedacartwright]

We can't automate it, but maybe we can ask Ensembl to provide the information we need at every release so stdpopsim can download it from the REST api.

[jeromekelleher]

I just had a go at pulling in the mapping files @ndukler - it works, but it's messy. I'll email the Ensembl devs tomorrow morning, they must have to deal with this all the time.

[jeromekelleher]

I've email some folks I know at EBI. Will report back.

[jeromekelleher]

So it turns out Ensembl do provide synonyms in the REST API. Doh, it was right there in the docs all the time. It's not that extensive, but it'll do.

We don't have sex chr information from Ensembl, but I figure that's an easy enough thing for us to curate manually. It's not like that sort of thing will ever change.

[ndukler]

Ok, we just need to standardize the file format for this kind of thing. Also are circular chromosomes supported in msprime or slim?

[jeromekelleher]

Also are circular chromosomes supported in msprime or slim?

No, unlikely to be implemented in msprime either, I'd imagine.

[grahamgower]

What's the problem with circular chromosomes? Is it just recombination?

[reedacartwright]

You probably don't need to do anything special for them except disable recombination. Even if you are simulating indels modeling the ends likely doesn't add too much to your analysis.

[jeromekelleher]

Oh, OK. I'd imagined these were circular chromosomes with recombination. If there's no recomb, then I agree with @reedacartwright

[jeromekelleher]

Closing - we've settled on using Ensembl for the moment.