<p>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:</p>

<pre><code><span class="math-container">$ 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 &amp;&amp; a[var]!=gt){print $1,$2,$3,$4,$</span>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
</code></pre>

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

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

<p>The basic idea is:</p>

<ul>
<li><code>awk -F"\t" -vOFS="\t"</code>: set the input and output field separators to a tab.</li>
<li><code>BEGIN{print "#CHROM ... }</code>: print a header.</li>
<li><code>if(/^[^#]/){</code>: skip the header lines of the input vcf files</li>
<li><code>var=$1$2$4$5</code>: get a "signature" for this variant, by concatenating the chromosome, position, ref and alt fields.</li>
<li><code>gt=substr($10,1,3);</code> : get the genotype value.</li>
<li><code>if(NR==FNR){</code>: if this is the 1st file/</li>
<li><code>a[var]=gt; next</code>: save the genotype for this variant and skip to the next line.</li>
<li><code>if(var in a &amp;&amp; a[var]!=gt){</code>: 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.</li>
</ul>

<hr>

<p>If you need the error rate, you can extend the script thusly:</p>

<pre><code>awk -F"\t" -vOFS="\t" '
BEGIN{
    print "#CHROM","POS","ID","REF","ALT","Truth","Test"
}
{ 
    if(/^[^#]/){
        var=<span class="math-container">$1$</span>2<span class="math-container">$4$</span>5; 
        gt=substr(<span class="math-container">$10,1,3); 
        if(NR==FNR){
            a[var]=gt; 
            next;   
        }
        tot++; 
        if(var in a &amp;&amp; a[var]!=gt){
            bad++;
            print $1,$2,$3,$4,$</span>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 
</code></pre>

<p>On your example dataset, that prints:</p>

<pre><code>#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%)
</code></pre>

CC BY-SA 4.0
2019-03-13T18:21:10.230Z
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.
user172818
CC BY-SA 4.0
2019-03-13T18:33:11.587Z
@user172818 oh wow, thanks! I hadn't encountered the `PS` before. Does this mean my simplistic approach is basically useless?
terdon
CC BY-SA 4.0
2019-03-13T20:55:09.677Z
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.
user172818
CC BY-SA 4.0
2019-03-14T11:59:32.257Z
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.
Wouter De Coster
CC BY-SA 4.0
2019-03-13T13:32:40.913Z
2019-03-13T18:19:08.557Z
2019-03-13T18:19:08.557Z
terdon
<p>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 <a href="https://whatshap.readthedocs.io/en/latest/" rel="nofollow noreferrer">WhatsHap</a> 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.</p>

<p>Example of test lines:</p>

<pre><code>#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
</code></pre>

<p>Example of truth lines</p>

<pre><code>#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


</code></pre>

CC BY-SA 4.0
2019-03-13T12:00:11.097Z
Could you give us a few lines of each file to test on? Do you just want to compare the genotype field for the same variant in the two vcf files?
terdon
CC BY-SA 4.0
2019-03-13T13:10:32.733Z
I've edited my question to include example lines. The truth file has every variant phased (it is created based on the synthetic diploid method https://www.nature.com/articles/s41592-018-0054-7 and https://github.com/lh3/CHM-eval/tree/master/dip-call).
Wouter De Coster
<vcf><haplotypes><phasing>
Searching tool to calculate phase/switch error rate
# Searching tool to calculate phase/switch error rate

(Asked by: username_0 on 2019-03-12)

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

**Tags:** vcf, haplotypes, phasing
> **Comments:**
> **username_1**: Could you give us a few lines of each file to test on? Do you just want to compare the genotype field for the same variant in the two vcf files?  
> **username_0**: I've edited my question to include example lines. The truth file has every variant phased (it is created based on the synthetic diploid method https://www.nature.com/articles/s41592-018-0054-7 and https://github.com/lh3/CHM-eval/tree/master/dip-call).  

---

## 1 answers

### Answer by username_1 on 2019-03-13. Score: 3

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%)

> **Comments:**
> **username_2**: 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.  
> **username_1**: @user172818 oh wow, thanks! I hadn't encountered the `PS` before. Does this mean my simplistic approach is basically useless?  
> **username_2**: 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.  
> **username_0**: 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.  

---

