Meaning of the FORMAT fields of the VCF file coming from GIAB project

Question

After reading the GIAB paper in https://www.biorxiv.org/content/early/2018/05/25/281006 and its Figure 1, I am still having trouble understanding the data inside the GIAB VCF file for HG001 (HG001_GRCh38_GIAB_highconf_CG-IllFB-IllGATKHC-Ion-10X-SOLID_CHROM1-X_v.3.3.2_all.vcf). In particular, I need to understand the information that led to a FILTER value of PASS for some variants.

The first case that I would like to discuss is the following variant:

chr1    110176467   .   A   G   50  PASS    platforms=3;platformnames=Illumina,CG,10X;datasets=3;datasetnames=HiSeqPE300x,CGnormal,10XChromium;callsets=4;callsetnames=HiSeqPE300xGATK,CGnormal,HiSeqPE300xfreebayes,10XGATKhaplo;datasetsmissingcall=IonExome,SolidPE50x50bp,SolidSE75bp;callable=CS_HiSeqPE300xGATK_callable,CS_CGnormal_callable,CS_HiSeqPE300xfreebayes_callable    GT:PS:DP:ADALL:AD:GQ    0/1:.:512:114,106:145,148:760

Translating the information in the line to English, this is what I get:

• We considered data from the following sequencing technologies: Illumina,CG,10X,IonExome,SolidPE50x50bp, and SolidSE75bp.
• Of those, some of them (IonExome,SolidPE50x50bp,SolidSE75bp) were missing the call, so we discarded them.
• Therefore, the data for the variant is going to come from 3 platforms (Illumina,CG,10X), 3 datasets (HiSeqPE300x,CGnormal,10XChromium), and 4 callsets (HiSeqPE300xGATK,CGnormal,HiSeqPE300xfreebayes,10XGATKhaplo), because we analyzed the data from HiSeqPE300x with GATK and also with Freebayes.
• However, the callsets from CS_HiSeqPE300xGATK,CS_CGnormal, and CS_HiSeqPE300xfreebayes had this call in a region with low coverage of high MQ reads.
• Across all the platforms, DP=512
• ADALL (Values 114,106) is coming from “all" the datasets. Are those “all” the set (HiSeqPE300x + Cgnormal + 10XChromium), or the set (Illumina + CG + 10X + IonExome + SolidPE50x50bp + SolidSE75bp).
• AD (Values 145,148) is coming from the same “all” datasets as ADALL, but they are unfiltered, so 145>114 and 148>106.

All this make sense, but I would like to know what case of Figure 1c was applied to arrive at the PASS value.

The second variant that I would like to mention is this one:

chr1    5705293 .   T   C   50  PASS    platforms=3;platformnames=Illumina,CG,10X;datasets=3;datasetnames=HiSeqPE300x,CGnormal,10XChromium;callsets=4;callsetnames=HiSeqPE300xGATK,CGnormal,HiSeqPE300xfreebayes,10XGATKhaplo;datasetsmissingcall=IonExome,SolidPE50x50bp,SolidSE75bp;callable=CS_HiSeqPE300xGATK_callable,CS_CGnormal_callable,CS_HiSeqPE300xfreebayes_callable;filt=CS_HiSeqPE300xGATK_filt   GT:PS:DP:ADALL:AD:GQ    0/1:.:599:121,139:49,28:627

Again, my translation would be:

• We considered data from the following sequencing technologies: Illumina,CG,10X,IonExome,SolidPE50x50bp, and SolidSE75bp.

• Of those, some of them (IonExome,SolidPE50x50bp,SolidSE75bp) were missing the call, so we discarded them.

• Therefore, the data for the variant is going to come from 3 platforms (Illumina,CG,10X), 3 datasets (HiSeqPE300x,CGnormal,10XChromium), and 4 callsets (HiSeqPE300xGATK,CGnormal,HiSeqPE300xfreebayes,10XGATKhaplo), because we analyzed the data from HiSeqPE300x with GATK and also with Freebayes.

• However, the callsets from CS_HiSeqPE300xGATK,CS_CGnormal, and CS_HiSeqPE300xfreebayes had this call in a region with low coverage of high MQ reads.
• Also, the data from CS_HiSeqPE300xGATK_filt was filtered/discarded
• Across all the platforms, DP=599
• ADALL (Values 121,139) is coming from “all" the datasets. I assume (Cgnormal + 10XChromium), because CS_HiSeqPE300xGATK_filt was filtered.
• AD (Values 49,28) are unfiltered.

But notice that 49<121 and 28<139. How come that there are less unfiltered reads (AD)? Should not be the case that AD > ADALL for all instances?

I don’t know how to reconcile the two variants. Aren’t they providing opposite information in the AD and ADALL fields? Which specific datasets and callsets are involved in getting AD and ADALL for each variant?

I would also appreciate answers to a couple of other questions:

• Why are these two variants a PASS if there are 3 “callable” datasets in regions of low coverage?

• Did the variants follow an arbitration process? I guess not, if the INFO flag “arbitrated” is missing.

For reference, here are the descriptions in the header of the VCF file that I have been using:

##INFO=<ID=platforms,Number=1,Type=Integer,Description="Number of different platforms for which at least one callset called this genotype, whether filtered or not">
##INFO=<ID=platformnames,Number=.,Type=String,Description="Names of platforms for which at least one callset called this genotype, whether filtered or not">
##INFO=<ID=datasets,Number=1,Type=Integer,Description="Number of different datasets for which at least one callset called this genotype, whether filtered or not">
##INFO=<ID=datasetnames,Number=.,Type=String,Description="Names of datasets for which at least one callset called this genotype, whether filtered or not">
##FORMAT=<ID=ADALL,Number=R,Type=Integer,Description="Net allele depths across all datasets">
##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Net allele depths across all unfiltered datasets with called genotype">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Total read depth summed across all datasets, excluding MQ0 reads">
##INFO=<ID=datasetsmissingcall,Number=.,Type=Integer,Description="Names of datasets that are missing a call or have an incorrect call at this location, and the high-confidence call is a variant">
##INFO=<ID=callable,Number=.,Type=String,Description="List of callsets that had this call in a region with low coverage of high MQ reads.">
##INFO=<ID=filt,Number=.,Type=String,Description="List of callsets that had this call filtered.">
##INFO=<ID=arbitrated,Number=1,Type=String,Description="TRUE if callsets had discordant calls so that arbitration was needed.">

Answer 1

The filtering criteria should be defined in the header of the VCF file. Can you include the header from your VCF file?

VCF Format and example header:

FILTERs that have been applied to the data should be described as follows:

##FILTER=<ID=ID,Description=”description”>

header:

##fileformat=VCFv4.0
##fileDate=20090805
##source=myImputationProgramV3.1
##reference=1000GenomesPilot-NCBI36
##phasing=partial
##INFO=<ID=NS,Number=1,Type=Integer,Description="Number of Samples With Data">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
##INFO=<ID=AF,Number=.,Type=Float,Description="Allele Frequency">
##INFO=<ID=AA,Number=1,Type=String,Description="Ancestral Allele">
##INFO=<ID=DB,Number=0,Type=Flag,Description="dbSNP membership, build 129">
##INFO=<ID=H2,Number=0,Type=Flag,Description="HapMap2 membership">
##FILTER=<ID=q10,Description="Quality below 10">
##FILTER=<ID=s50,Description="Less than 50% of samples have data">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth">
##FORMAT=<ID=HQ,Number=2,Type=Integer,Description="Haplotype Quality">

Answer 2

For your second question, I don't think that variant matches any of the examples shown in Figure 1c. Note that that panel's caption says it shows four examples of arbitration. Per you last sub-question, I don't think this variant was subjected to arbitration.

For your second question, I believe the subtlety is in the field descriptions:

ADALL = Net allele depths across all datasets

AD = Net allele depths across all unfiltered datasets with called genotype

My understanding is that ADALL is a count of allele depths that includes samples where the variant was filtered out or the genotype was not concordant; something to that effect. So generally I would expect ADALL > AD.

I'm honestly not sure exactly what they mean by "unfiltered datasets". That could mean "counting all reads without filtering on quality" or "callsets to which variant filtering was not yet applied". Maybe there's another subtlety there that explains why the two variants have opposite AD vs ADALL proportions.

For your additional questions:

Why are these two variants a PASS if there are 3 “callable” datasets in regions of low coverage?

According to the paper, having the variant in low-coverage regions isn't disqualifying if the genotype calls are concordant, which makes sense to me. Low coverage makes me more suspicious, but if multiple callsets are in agreement, I'm still willing to consider the call good (unless of course they all exhibit the same error mode, which is the main potential weakness of the consensus approach used by GiaB).

Did the variants follow an arbitration process? I guess not, if the INFO flag “arbitrated” is missing.

Right, if that flag is not present, it means arbitration was not needed/applied.