text: user681: 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
  
  
  
  
  user298: 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%)
  
  
original_text: 681: <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>
   298: <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>
  