Interpreting short indel calls in 1000 Genomes Project VCFs

Question

Consider the following short indel polymorphism rs59679400 on chr7.

$ rsidx search dbSNP/dbSNP_GRCh38.{vcf.gz,rsidx} rs59679400 | cut -f 1-7
[rsidx] running version 0.3
NC_000007.14    12510038        rs59679400      GAAAA   G,GA,GAA,GAAA,GAAAAA    .       .

In the Phase III (2015) 1000 Genomes VCFs (GRCh37), I see the following record for the variant.

$ rsidx search 1000GenomesPhaseIII/ALL.chr7.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.{vcf.gz,rsidx} rs59679400 | cut -f 1-7
[rsidx] running version 0.3
7       12549664        rs59679400      GA      G       100     PASS

That's easy enough to interpret, although it looks like only two of the several possible alleles at that locus are represented.

When I look at the VCFs from the latest NYGC update to the 1000 Genomes data, I see the following. (I'll note that it's really inconvenient that RS numbers are no longer included for any variants.)

$ tabix 1000GenomesNYGC/1kGP_high_coverage_Illumina.chr7.filtered.SNV_INDEL_SV_phased_panel.vcf.gz chr7:12510038-12510038 | cut -f 1-7
chr7    12474780        HGSV_106040     G       <DUP>   .       .
chr7    12510038        7:12510038:GA:G GA      G       .       .
chr7    12510038        7:12510039:GAA:G        GAA     G       .       .

Ignoring the first record, which is a larger structural variant, we see two records at the position in question. My initial question is why they've used two records to encode variation for a single indel variant? As shown in the dbSNP record in the first snippet, VCF allows one to encode an arbitrary number of alternate alleles in a single record.

My second—and more important—question is how one is to interpret these records and resolve the sequence for any particular individual genotyped at the locus? Selecting an arbitrary sample HG00103, I get the following allele calls for the two records? What am I to make of this?

>>> from pysam import VariantFile
>>> vcf = VariantFile("1000GenomesNYGC/1kGP_high_coverage_Illumina.chr7.filtered.SNV_INDEL_SV_phased_panel.vcf.gz")
>>> records = list(vcf.fetch("chr7", 12510037, 12510038))
>>> len(records)
3
>>> records[1].samples["HG00103"].alleles
('G', 'G')
>>> records[2].samples["HG00103"].alleles
('GAA', 'GAA')

Answer 1

this is Marta from the NYGC. Thanks for reaching out with your questions.

The full SNV/INDEL callset that we released (listed as "SNV/INDEL VCFs" in the key resource table in our paper) includes multiallelic representation of all multiallelic sites (i.e. multiple ALT alleles are listed together in column 5 separated by commas) and the record you showed as an example above is listed there in the following way:

chr7 12510038 . GAAA GAA,GA,G

The VCF file you were looking at is the phased panel that we generated and released alongside the full unphased callset. The phasing software tools that we used in the paper (SHAPEIT2, SHAPEIT4, Eagle2) as well as most other phasing tools do not support multiallelic sites and require that such sites are split into separate lines and represented as biallelic for the purposes of phasing. As you noted, this is not ideal and can lead to difficulties in interpretation (as others explained above) but it's the best we can do with the existing approaches. We mention this in the Methods section of the paper as well as in the README that's included with the phased panel on the IGSR ftp (all linked in the paper). Please let us know if you have any further questions related to the high-coverage 1kGP resources. Thanks!

Answer 2

Older GATK was using the so-called variant model where overlapping variants were combined during calling. More recent GATK is using the allelic model where one non-reference allele is evaluated at a time.

With the allelic model, a reference allele doesn't mean the sample possesses a reference allele; it means the sample doesn't have the alternate allele. For this particular example, records[1] means a homozygous -A deletion. records[2] means the sample doesn't have the -AA deletion. The two records are consistent with each other.

It is still possible to have conflictive records (for example, homozygous -A deletion and homozygous -AA deletion at the same location). At least this used to be a problem with GATK. Not sure if/how much this has been improved recently.

Answer 3

Often multi-allelic sites would be split into several lines (one per alternate allele) to make filtering easier for specific types of alternate alleles (e.g. filter out MNV vs SNV).

VCF-level, one thing you can do to identify the actual genotype over the locus easier is to join the split alleles into a single multi-allelic record. You can easily do this using bcftools:

bcftools view --samples HG00103 --regions chr7:12510038-12510038 1kGP_high_coverage_Illumina.chr7.filtered.SNV_INDEL_SV_phased_panel.vcf.gz | bcftools norm -m + | bcftools view -H

From here is clear that starting on position chr7:12510038, instead of "GAA", you have homozygous "GA" (1|1), therefore a hmozygous "A" deletion, which could be either the first or second "A" at positions 12510039 and 12510040.