VCF info and format fields #471

tomwhite · 2021-02-23T17:42:42Z

This changes the VCF reader to handle INFO and FORMAT fields of type Flag, Integer, or Float, and where number is 1, as suggested in #464.

A few things for discussion:

Unlike scikit-allel which has a single fields list, I've created two, info_fields and format_fields to avoid ambiguities (like DP, which can apply to either), Alternatively, we could have a single list, and namespace them, e.g. FORMAT/DP.
Naming. As it stands INFO fields are prefixed with variant_, and FORMAT fields with call_. So call_DP is an example. Do we want to add a renaming keyword (like scikit-allel has) so it could be renamed to e.g. call_depth?
The GT FORMAT field is special since it is always included, even if format_fields doesn't include it. Should we make it optional?
The vcf_to_zarr_sequential is getting a bit unwieldy and should probably be refactored into a class or set of classes to handle various field types.

jeromekelleher · 2021-02-23T18:29:10Z

Naming. As it stands INFO fields are prefixed with variant_, and FORMAT fields with call_. So call_DP is an example. Do we want to add a renaming keyword (like scikit-allel has) so it could be renamed to e.g. call_depth?

One quick comment here: I'd tend to avoid renaming stuff from VCF. Ultimately we're just representing the VCF data so we should keep as close to the original as we can, IMO. By using call_DP it's clear we're not doing anything clever, just showing the user the value of the DP field in the VCF.

jeromekelleher · 2021-02-23T18:39:36Z

This is great @tomwhite! I've had a quick skim through and it looks very nice.

I think it would be worthwhile building up a library of weird and wonderful VCF examples at this point, I suspect some of the assumptions being made won't hold up in practise (that was my repeated experience with the VCF parser in wormtable. Using htslib/cyvcf2 is a huge help, but I suspect people play pretty fast and loose with INFO field types and so on, and we'll need to provide the user with ways of defining the numpy dtype, etc, themselves.

timothymillar · 2021-02-24T08:43:51Z

sgkit/io/vcf/vcf_reader.py

+                for info_field, arr in info_field_arrays.items():
+                    if info_field in [k for k, _ in variant.INFO]:
+                        if isinstance(variant.INFO[info_field], tuple):
+                            # for the case Number='.', just use the first value


It might be better to also throw an error for the '.' case as this is commonly used indicate that the field is variable in length, see this old htslib issue

Agreed - let's not support this until it's done properly.

tomwhite · 2021-02-24T12:00:35Z

I think it would be worthwhile building up a library of weird and wonderful VCF examples at this point

I've used test data from GATK for this in the past. It also looks like scikit-allel has some test VCFs that have been chosen/altered to test issues such as these.

tomwhite · 2021-02-24T12:34:32Z

Here's an update which better encapsulates the INFO and FORMAT field parsing, and removes support for the Number '.' case.

timothymillar · 2021-02-24T21:55:18Z

sgkit/io/vcf/vcf_reader.py

+    def add_variant(self, i: int, variant: Any) -> None:
+        key = self.key
+        if self.category == "INFO":
+            if key in [k for k, _ in variant.INFO]:


I think you can use variant.INFO.get(key) here

jeromekelleher · 2021-02-25T08:18:09Z

Here's a ready-made set of example VCFs @tomwhite: https://github.com/ga4gh/ga4gh-server/tree/master/tests/data/datasets/dataset1/variants

This is from the old GA4GH REST API reference implementation, where we did try to support all the INFO/GT fields in VCFs.

tomwhite · 2021-03-02T16:41:04Z

I've overhauled this quite a bit with the following changes:

Support all VCF Type fields (Integer, Float, String, Flag)
Support the following VCF Number fields: any integer, A, R (but not . or G).
Allow field definitions to be overridden using the field_defs argument. (This is useful to override a . Number field with a fixed number, for example.)
Allow fields to be specified by a wildcard (e.g. INFO/*).
Convert INFO and FORMAT fields from scikit-allel format, and add roundtrip tests to check equivalence. This is a useful check since it's effectively testing two VCF parsers against one another.

I've checked equivalence for all the test files in sgkit's test data, and it would be fairly straightforward to extend this to other sets of VCF files. That could be done separately though.

timothymillar · 2021-03-02T19:22:38Z

This looks great @tomwhite. The size of "G" can actually be calculated from n_alleles and ploidy using n-choose-k with replacement

from math import comb
comb(n_alleles + ploidy - 1, ploidy)

hammer · 2021-03-02T19:46:26Z

@timothymillar do you have any VCFs with data from variable ploidy samples that you might be able to contribute to our testing efforts?

timothymillar · 2021-03-03T06:31:28Z

@hammer I'm looking into this, we can definitely share some data for (a chromosome of) a 4x potato population which is already public. I'm not aware of a true mixed ploidy example but at the very least I can call the 4x material with a range of ploidy values to make something similar.

jeromekelleher

This is fantastic @tomwhite, minor comments above! Probably some of these should be pushed out to follow-up PRs.

jeromekelleher · 2021-03-03T09:01:53Z

sgkit/io/vcf/vcf_reader.py

+        return "bool", False
+    elif vcf_type == "Integer":
+        return "i4", -1
+    elif vcf_type == "Float":


Worth putting a note in here that BCF stores these values as single precision, so there's no point in using doubles.

jeromekelleher · 2021-03-03T09:07:27Z

sgkit/io/vcf/vcf_reader.py

+    def update_dataset(
+        self, ds: xr.Dataset, add_str_max_length_attrs: bool = False
+    ) -> None:
+        # cyvcf2 represents missing Integer values as the minimum int32 value


I'm slightly worried about this - could we be silently altering data? Suppose someone uses INT32_MIN in their data, won't we silently change this to fill_value then? I guess this is an upstream problem with cyvcf2.

I opened brentp/cyvcf2#195 for this.

I don't think they could use INT32_MIN for anything else as this value is explicitly interpreted as a 'missing' integer in the VCF spec (page 30)

jeromekelleher · 2021-03-03T09:10:44Z

sgkit/io/vcf/vcf_reader.py

@@ -261,6 +467,24 @@ def vcf_to_zarrs(
        If True, genotype calls with more alleles than the specified (maximum) ploidy value
        will be truncated to size ploidy. If false, calls with more alleles than the
        specified ploidy will raise an exception.
+    max_alt_alleles
+        The (maximum) number of alternate alleles in the VCF file. Any records with more than
+        this number of alternate alleles will have the extra alleles dropped.


Perhaps not for this PR, but we should probably raise a warning if this happens. Users should know if their ALTs are being dropped, and ideally it should be discoverable later if they have been truncated. Could we add some metadata tracking when this has happened?

Opened #487

jeromekelleher · 2021-03-03T09:11:04Z

sgkit/io/vcf/vcf_reader.py

+        this number of alternate alleles will have the extra alleles dropped.
+    fields
+        Extra fields to extract data for. A list of strings, with ``INFO`` or ``FORMAT`` prefixes.
+        Wildcards are permitted too, for example: ``["INFO/*", "FORMAT/DP"]``.


jeromekelleher · 2021-03-03T09:15:00Z

sgkit/io/vcf/vcf_reader.py

+        ``dimension``. The first three correspond to VCF header values, and ``dimension`` is
+        the name of the final dimension in the array for the case where ``Number`` is a fixed
+        integer larger than 1. For example,
+        ``{"INFO/AC": {"Number": "A"}, "FORMAT/HQ": {"dimension": "haplotypes"}}``


Good to explain what "A" means here; maybe we should say something like "see section 1.2.2 of the VCF specification <url>__ for further details on the available Number specifications."

mergify · 2021-03-08T09:01:48Z

This PR has conflicts, @tomwhite please rebase and push updated version 🙏

dtype) for arbitrary dataset variables.

tomwhite added the data representation Issues related to how data is represented: data types, data structures, indexes, access methods, etc label Feb 23, 2021

tomwhite mentioned this pull request Feb 23, 2021

Replicating Hail's GWAS Tutorial #463

Closed

timothymillar reviewed Feb 24, 2021

View reviewed changes

tomwhite force-pushed the vcf-info-and-format-fields branch from 9c17620 to f711779 Compare March 2, 2021 16:40

jeromekelleher approved these changes Mar 3, 2021

View reviewed changes

mergify bot added the conflict PR conflict label Mar 8, 2021

tomwhite added 4 commits March 9, 2021 09:58

Generalize enforcement of chunk sizes

3a55f79

Generalize mechanism for computing max length of strings (for setting

b0883f9

dtype) for arbitrary dataset variables.

Support reading INFO and FORMAT fields from VCF files.

0c215e5

vcfzarr roundtrip

71d3352

tomwhite mentioned this pull request Mar 9, 2021

Issue warning if more alternate alleles than max_alt_alleles when reading VCF #487

Closed

Add links to VCF spec for reading INFO/FORMAT fields

f3a919c

tomwhite force-pushed the vcf-info-and-format-fields branch from f711779 to f3a919c Compare March 9, 2021 10:34

mergify bot removed the conflict PR conflict label Mar 9, 2021

Fix mypy error

b5e72ac

tomwhite added the auto-merge Auto merge label for mergify test flight label Mar 9, 2021

mergify bot merged commit 10e3e01 into sgkit-dev:master Mar 9, 2021

tomwhite mentioned this pull request Mar 16, 2021

Support Number=G VCF fields #493

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

VCF info and format fields #471

VCF info and format fields #471

tomwhite commented Feb 23, 2021

jeromekelleher commented Feb 23, 2021

jeromekelleher commented Feb 23, 2021

timothymillar Feb 24, 2021

tomwhite Feb 24, 2021

tomwhite commented Feb 24, 2021

tomwhite commented Feb 24, 2021

timothymillar Feb 24, 2021

jeromekelleher commented Feb 25, 2021

tomwhite commented Mar 2, 2021

timothymillar commented Mar 2, 2021

hammer commented Mar 2, 2021

timothymillar commented Mar 3, 2021

jeromekelleher left a comment

jeromekelleher Mar 3, 2021

jeromekelleher Mar 3, 2021

tomwhite Mar 9, 2021

timothymillar Mar 9, 2021

jeromekelleher Mar 3, 2021

tomwhite Mar 9, 2021

jeromekelleher Mar 3, 2021

jeromekelleher Mar 3, 2021

mergify bot commented Mar 8, 2021

VCF info and format fields #471

VCF info and format fields #471

Conversation

tomwhite commented Feb 23, 2021

jeromekelleher commented Feb 23, 2021

jeromekelleher commented Feb 23, 2021

Choose a reason for hiding this comment

Choose a reason for hiding this comment

tomwhite commented Feb 24, 2021

tomwhite commented Feb 24, 2021

Choose a reason for hiding this comment

jeromekelleher commented Feb 25, 2021

tomwhite commented Mar 2, 2021

timothymillar commented Mar 2, 2021

hammer commented Mar 2, 2021

timothymillar commented Mar 3, 2021

jeromekelleher left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mergify bot commented Mar 8, 2021