5
$\begingroup$

I'm looking for a tool which, given a truth vcf file and a test vcf file, calculates the phase/switch error rate. I performed phasing of a vcf using WhatsHap and want to compare the outcome to some ground truth phased vcf I have. I can't find a tool, and I don't want to write it myself. It seems like a rather common thing... so I'd expect it to exist, but my google-fu is weak today.

Example of test lines:

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  sample

chr10   100588  .       G       A       .       PASS    AF=0.18;AC=1;NS=150;AN=2;EAS_AF=0.14;EUR_AF=0.39;AFR_AF=0.08;AMR_AF=0.31;SAS_AF=0.2;VT=SNP;DP=11907     GT:PS   1|0:48005
chr10   102385  .       A       G       .       PASS    AF=0.07;AC=1;NS=150;AN=2;EAS_AF=0.02;EUR_AF=0.0;AFR_AF=0.1;AMR_AF=0.19;SAS_AF=0.01;VT=SNP;DP=22665      GT:PS   0|1:48005
chr10   105170  .       G       A       .       PASS    AF=0.03;AC=1;NS=150;AN=2;EAS_AF=0.02;EUR_AF=0.0;AFR_AF=0.0;AMR_AF=0.19;SAS_AF=0.0;VT=SNP;DP=17140       GT:PS   0|1:48005
chr10   105365  .       G       T       .       PASS    AF=0.19;AC=1;NS=150;AN=2;EAS_AF=0.14;EUR_AF=0.39;AFR_AF=0.09;AMR_AF=0.31;SAS_AF=0.2;VT=SNP;DP=21845     GT:PS   1|0:48005
chr10   106057  .       T       C       .       PASS    AF=0.13;AC=1;NS=150;AN=2;EAS_AF=0.02;EUR_AF=0.0;AFR_AF=0.32;AMR_AF=0.19;SAS_AF=0.01;VT=SNP;DP=22559     GT:PS   0|1:48005

Example of truth lines

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  sample
chr10   100554  .       C       CT,CTT  30      .       .       GT:AD   1|2:0,1,1
chr10   100588  .       G       A       30      .       .       GT:AD   0|1:1,1
chr10   102385  .       A       G       30      .       .       GT:AD   1|0:1,1
chr10   102636  .       T       C       30      .       .       GT:AD   1|1:0,2
chr10   102757  .       T       C       30      .       .       GT:AD   0|1:1,1
chr10   104815  .       GCC     G,GC    30      .       .       GT:AD   2|1:0,1,1
chr10   105170  .       G       A       30      .       .       GT:AD   1|0:1,1
chr10   105365  .       G       T       30      .       .       GT:AD   0|1:1,1
chr10   106057  .       T       C       30      .       .       GT:AD   1|0:1,1
chr10   106110  .       C       G       30      .       .       GT:AD   1|1:0,2
chr10   106261  .       A       G       30      .       .       GT:AD   1|0:1,1
chr10   108612  .       T       G       30      .       .       GT:AD   1|0:1,1
chr10   108646  .       A       AT      30      .       .       GT:AD   0|1:1,1
chr10   108834  .       G       A       30      .       .       GT:AD   1|0:1,1
chr10   110840  .       C       T       30      .       .       GT:AD   1|1:0,2
chr10   111743  .       A       AT      30      .       .       GT:AD   1|0:1,1
chr10   112016  .       T       A       30      .       .       GT:AD   1|0:1,1
chr10   112262  .       C       T       30      .       .       GT:AD   1|1:0,2
chr10   113006  .       C       T       30      .       .       GT:AD   0|1:1,1
chr10   113031  .       G       A       30      .       .       GT:AD   0|1:1,1
chr10   113136  .       G       A       30      .       .       GT:AD   1|0:1,1
chr10   113359  .       A       T       30      .       .       GT:AD   1|1:0,2
chr10   113583  .       C       G       30      .       .       GT:AD   0|1:1,1
chr10   113997  .       G       C       30      .       .       GT:AD   1|1:0,2
Improve this question
$\endgroup$
2

1 Answer 1

Reset to default
3
$\begingroup$

Assuming your files are as you show, and the genotype field is always the 1st entry of the sample information field (the 10th), then you can just do this directly with a simple awk script:

$ awk -F"\t" -vOFS="\t" 'BEGIN{print "#CHROM","POS","ID","REF","ALT","Truth","Test"}{ if(/^[^#]/){var=$1$2$4$5; gt=substr($10,1,3); if(NR==FNR){a[var]=gt; next;} if(var in a && a[var]!=gt){print $1,$2,$3,$4,$5,gt,a[var]}}}' truth.vcf test.vcf 
#CHROM  POS     ID  REF ALT Truth   Test
chr10   100588  .   G   A   1|0 0|1
chr10   102385  .   A   G   0|1 1|0
chr10   105170  .   G   A   0|1 1|0
chr10   105365  .   G   T   1|0 0|1
chr10   106057  .   T   C   0|1 1|0

Here's the script split into separate lines for clarity (you can still copy paste directly into a terminal):

awk -F"\t" -vOFS="\t" '
    BEGIN{
        print "#CHROM","POS","ID","REF","ALT","Truth","Test"
    }
    { 
        if(/^[^#]/){
            var=$1$2$4$5; 
            gt=substr($10,1,3); 
            if(NR==FNR){
                a[var]=gt; 
                next;   
            } 
            if(var in a && a[var]!=gt){
                print $1,$2,$3,$4,$5,gt,a[var]
            }
        }
    }' truth.vcf test.vcf

The basic idea is:

  • awk -F"\t" -vOFS="\t": set the input and output field separators to a tab.
  • BEGIN{print "#CHROM ... }: print a header.
  • if(/^[^#]/){: skip the header lines of the input vcf files
  • var=$1$2$4$5: get a "signature" for this variant, by concatenating the chromosome, position, ref and alt fields.
  • gt=substr($10,1,3); : get the genotype value.
  • if(NR==FNR){: if this is the 1st file/
  • a[var]=gt; next: save the genotype for this variant and skip to the next line.
  • if(var in a && a[var]!=gt){: this part will only be run for the 2nd file because of the 'next' in the previous block. If this variant was in the truth set and if its genotype there isn't the same as in the test set, print the variant and the two genotypes.

If you need the error rate, you can extend the script thusly:

awk -F"\t" -vOFS="\t" '
BEGIN{
    print "#CHROM","POS","ID","REF","ALT","Truth","Test"
}
{ 
    if(/^[^#]/){
        var=$1$2$4$5; 
        gt=substr($10,1,3); 
        if(NR==FNR){
            a[var]=gt; 
            next;   
        }
        tot++; 
        if(var in a && a[var]!=gt){
            bad++;
            print $1,$2,$3,$4,$5,gt,a[var]
        }
    }
}
END{
    printf "Correct: %s of %s (%s%)\n",tot-bad,tot,(tot-bad)*100/tot;
    printf "Wrong: %s of %s (%s%)\n",bad,tot,bad*100/tot;
}' truth.vcf test.vcf

On your example dataset, that prints:

#CHROM  POS ID  REF ALT Truth   Test
chr10   100588  .   G   A   1|0 0|1
chr10   102385  .   A   G   0|1 1|0
chr10   105170  .   G   A   0|1 1|0
chr10   105365  .   G   T   1|0 0|1
chr10   106057  .   T   C   0|1 1|0
Correct: 0 of 5 (0%)
Wrong: 5 of 5 (100%)
Improve this answer
$\endgroup$
4
  • 2
    $\begingroup$ In the test VCF, "PS" means "phase set". Switches between phase sets are not counted as errors. In general, we do need a proper script for evaluation. There is not a common script probably because the inputs often change a lot. $\endgroup$
    – user172818
    Mar 13, 2019 at 18:21
  • $\begingroup$ @user172818 oh wow, thanks! I hadn't encountered the PS before. Does this mean my simplistic approach is basically useless? $\endgroup$
    – terdon ♦
    Mar 13, 2019 at 18:33
  • $\begingroup$ I think it will be useful when "PS" is not involved. For example, Hi-C phasing is supposed to give whole-chromosome phasing. We don't need "PS". Your script is applicable in that case. $\endgroup$
    – user172818
    Mar 13, 2019 at 20:55
  • $\begingroup$ Yes, while the truth set has "left" and "right" haplotypes (paternal or maternal) the test set is split in haplotype/phase blocks. Within a block, we can figure out which SNPs are on the same haplotype (or aren't), but not across blocks. My example is not perfect, but I believe it's 100% correctly phased here. $\endgroup$
    – Wouter De Coster
    Mar 14, 2019 at 11:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.