Add stats for the number of sites matched in the GT-vs-GT, GT-vs-PL, etc modes.

This information is important for interpretation of the discordance score, because
only the GT-vs-GT matching can be interpreted as the number of mismatching genotypes.

Author: pd3

NEWS

@@ -45,6 +45,12 @@ Changes affecting specific commands:
       is newly printed as
             3_prime_utr&NMD_transcript|PCGF3|ENST00000430644|NMD
 
+* bcftools gtcheck
+
+    - Add stats for the number of sites matched in the GT-vs-GT, GT-vs-PL, etc modes. This
+      information is important for interpretation of the discordance score, as only the
+      GT-vs-GT matching can be interpreted as the number of mismatching genotypes.
+
 * bcftools +mendelian2
 
     - Fix in command line argument parsing, the `-p` and `-P` options were not

doc/bcftools.1

@@ -2,12 +2,12 @@
 .\"     Title: bcftools
 .\"    Author: [see the "AUTHOR(S)" section]
 .\" Generator: Asciidoctor 2.0.15.dev
-.\"      Date: 2023-05-30
+.\"      Date: 2023-06-02
 .\"    Manual: \ \&
 .\"    Source: \ \&
 .\"  Language: English
 .\"
-.TH "BCFTOOLS" "1" "2023-05-30" "\ \&" "\ \&"
+.TH "BCFTOOLS" "1" "2023-06-02" "\ \&" "\ \&"
 .ie \n(.g .ds Aq \(aq
 .el       .ds Aq '
 .ss \n[.ss] 0
@@ -51,7 +51,7 @@ standard input (stdin) and outputs to the standard output (stdout). Several
 commands can thus be  combined  with  Unix pipes.
 .SS "VERSION"
 .sp
-This manual page was last updated \fB2023\-05\-30 09:18 BST\fP and refers to bcftools git version \fB1.17\-50\-ga8249495+\fP.
+This manual page was last updated \fB2023\-06\-02 11:27 BST\fP and refers to bcftools git version \fB1.17\-52\-g0773541c+\fP.
 .SS "BCF1"
 .sp
 The obsolete BCF1 format output by versions of samtools <= 0.1.19 is \fBnot\fP
@@ -2478,6 +2478,10 @@ option is given, the identity of samples from \fIquery.vcf.gz\fP
 is checked against the samples in the \fB\-g\fP file.
 Without the \fB\-g\fP option, multi\-sample cross\-check of samples in \fIquery.vcf.gz\fP is performed.
 .sp
+Note that the interpretation of the discordance score depends on the options provided (specifically \fB\-e\fP and
+\fB\-u\fP) and on the available annotations (FORMAT/PL vs FORMAT/GT).
+The discordance score can be interpreted as the number of mismatching genotypes if only GT\-vs\-GT matching is performed.
+.sp
 \fB\-\-distinctive\-sites\fP \fINUM[,MEM[,DIR]]\fP
 .RS 4
 Find sites that can distinguish between at least NUM sample pairs. If the number is smaller or equal to 1,
@@ -2496,11 +2500,18 @@ Stop after first record to estimate required time.
 Interpret genotypes and genotype likelihoods probabilistically. The value of \fIINT\fP
 represents genotype quality when GT tag is used (e.g. Q=30 represents one error in 1,000 genotypes and
 Q=40 one error in 10,000 genotypes) and is ignored when PL tag is used (in that case an arbitrary
-non\-zero integer can be provided). See also the \fB\-u, \-\-use\fP option below. If set to 0,
-the discordance equals to the number of mismatching genotypes when GT vs GT is compared.
-Note that the values with and without \fB\-e\fP are not comparable, only values generated
-with \fB\-e 0\fP correspond to mismatching genotypes.
-If performance is an issue, set to 0 for faster run but less accurate results.
+non\-zero integer can be provided).
+\~
+.br
+\~
+.br
+If \fB\-e\fP is set to 0, the discordance score can be interpreted as the number of mismatching genotypes,
+but only in the GT\-vs\-GT matching mode. See the \fB\-u, \-\-use\fP option below for additional notes and caveats.
+\~
+.br
+\~
+.br
+If performance is an issue, set \fB\-e 0\fP for faster run times but less accurate results.
 .RE
 .sp
 \fB\-g, \-\-genotypes\fP \fIFILE\fP
@@ -2581,8 +2592,15 @@ see \fBCommon Options\fP
 \fB\-u, \-\-use\fP \fITAG1\fP[,\fITAG2\fP]
 .RS 4
 specifies which tag to use in the query file (\fITAG1\fP) and the \fB\-g\fP (\fITAG2\fP) file.
-By default, the PL tag is used in the query file and GT in the \fB\-g\fP file when
-available.
+By default, the PL tag is used in the query file and, when available, the GT tags in the
+\fB\-g\fP file.
+\~
+.br
+\~
+.br
+Note that when the requested tag is not available, the program will attempt to use
+the other tag. The output includes the number of sites that were matched by the four
+possible mode (for example GT\-vs\-GT or GT\-vs\-PL).
 .RE
 .sp
 \fBExamples:\fP

doc/bcftools.html

@@ -50,7 +50,7 @@ <h2 id="_description">DESCRIPTION</h2>
 <div class="sect2">
 <h3 id="_version">VERSION</h3>
 <div class="paragraph">
-<p>This manual page was last updated <strong>2023-05-30 09:18 BST</strong> and refers to bcftools git version <strong>1.17-50-ga8249495+</strong>.</p>
+<p>This manual page was last updated <strong>2023-06-02 11:27 BST</strong> and refers to bcftools git version <strong>1.17-52-g0773541c+</strong>.</p>
 </div>
 </div>
 <div class="sect2">
@@ -2178,6 +2178,11 @@ <h3 id="gtcheck">bcftools gtcheck [<em>OPTIONS</em>] [<strong>-g</strong> <em>ge
 is checked against the samples in the <strong>-g</strong> file.
 Without the <strong>-g</strong> option, multi-sample cross-check of samples in <em>query.vcf.gz</em> is performed.</p>
 </div>
+<div class="paragraph">
+<p>Note that the interpretation of the discordance score depends on the options provided (specifically <strong>-e</strong> and
+<strong>-u</strong>) and on the available annotations (FORMAT/PL vs FORMAT/GT).
+The discordance score can be interpreted as the number of mismatching genotypes if only GT-vs-GT matching is performed.</p>
+</div>
 <div class="dlist">
 <dl>
 <dt class="hdlist1"><strong>--distinctive-sites</strong> <em>NUM[,MEM[,DIR]]</em></dt>
@@ -2196,11 +2201,14 @@ <h3 id="gtcheck">bcftools gtcheck [<em>OPTIONS</em>] [<strong>-g</strong> <em>ge
 <p>Interpret genotypes and genotype likelihoods probabilistically. The value of <em>INT</em>
 represents genotype quality when GT tag is used (e.g. Q=30 represents one error in 1,000 genotypes and
 Q=40 one error in 10,000 genotypes) and is ignored when PL tag is used (in that case an arbitrary
-non-zero integer can be provided). See also the <strong>-u, --use</strong> option below. If set to 0,
-the discordance equals to the number of mismatching genotypes when GT vs GT is compared.
-Note that the values with and without <strong>-e</strong> are not comparable, only values generated
-with <strong>-e 0</strong> correspond to mismatching genotypes.
-If performance is an issue, set to 0 for faster run but less accurate results.</p>
+non-zero integer can be provided).
+&#160;<br>
+&#160;<br>
+If <strong>-e</strong> is set to 0, the discordance score can be interpreted as the number of mismatching genotypes,
+but only in the GT-vs-GT matching mode. See the <strong>-u, --use</strong> option below for additional notes and caveats.
+&#160;<br>
+&#160;<br>
+If performance is an issue, set <strong>-e 0</strong> for faster run times but less accurate results.</p>
 </dd>
 <dt class="hdlist1"><strong>-g, --genotypes</strong> <em>FILE</em></dt>
 <dd>
@@ -2274,8 +2282,13 @@ <h3 id="gtcheck">bcftools gtcheck [<em>OPTIONS</em>] [<strong>-g</strong> <em>ge
 <dt class="hdlist1"><strong>-u, --use</strong> <em>TAG1</em>[,<em>TAG2</em>]</dt>
 <dd>
 <p>specifies which tag to use in the query file (<em>TAG1</em>) and the <strong>-g</strong> (<em>TAG2</em>) file.
-By default, the PL tag is used in the query file and GT in the <strong>-g</strong> file when
-available.</p>
+By default, the PL tag is used in the query file and, when available, the GT tags in the
+<strong>-g</strong> file.
+&#160;<br>
+&#160;<br>
+Note that when the requested tag is not available, the program will attempt to use
+the other tag. The output includes the number of sites that were matched by the four
+possible mode (for example GT-vs-GT or GT-vs-PL).</p>
 </dd>
 </dl>
 </div>
@@ -5257,7 +5270,7 @@ <h2 id="_copying">COPYING</h2>
 </div>
 <div id="footer">
 <div id="footer-text">
-Last updated 2023-05-30 09:18:06 +0100
+Last updated 2023-06-02 11:27:10 +0100
 </div>
 </div>
 </body>

doc/bcftools.txt

@@ -1624,6 +1624,10 @@ option is given, the identity of samples from 'query.vcf.gz'
 is checked against the samples in the *-g* file.
 Without the *-g* option, multi-sample cross-check of samples in 'query.vcf.gz' is performed.
 
+Note that the interpretation of the discordance score depends on the options provided (specifically *-e* and
+*-u*) and on the available annotations (FORMAT/PL vs FORMAT/GT).
+The discordance score can be interpreted as the number of mismatching genotypes if only GT-vs-GT matching is performed.
+
 *--distinctive-sites* 'NUM[,MEM[,DIR]]'::
     Find sites that can distinguish between at least NUM sample pairs. If the number is smaller or equal to 1,
     it is interpreted as the fraction of pairs.  The optional MEM string sets the maximum memory used for
@@ -1637,11 +1641,14 @@ Without the *-g* option, multi-sample cross-check of samples in 'query.vcf.gz' i
     Interpret genotypes and genotype likelihoods probabilistically. The value of 'INT'
     represents genotype quality when GT tag is used (e.g. Q=30 represents one error in 1,000 genotypes and
     Q=40 one error in 10,000 genotypes) and is ignored when PL tag is used (in that case an arbitrary
-    non-zero integer can be provided). See also the *-u, --use* option below. If set to 0,
-    the discordance equals to the number of mismatching genotypes when GT vs GT is compared.
-    Note that the values with and without *-e* are not comparable, only values generated
-    with *-e 0* correspond to mismatching genotypes.
-    If performance is an issue, set to 0 for faster run but less accurate results.
+    non-zero integer can be provided).
+    {nbsp} +
+    {nbsp} +
+    If *-e* is set to 0, the discordance score can be interpreted as the number of mismatching genotypes,
+    but only in the GT-vs-GT matching mode. See the *-u, --use* option below for additional notes and caveats.
+    {nbsp} +
+    {nbsp} +
+    If performance is an issue, set *-e 0* for faster run times but less accurate results.
 
 *-g, --genotypes* 'FILE'::
     VCF/BCF file with reference genotypes to compare against
@@ -1696,8 +1703,13 @@ Without the *-g* option, multi-sample cross-check of samples in 'query.vcf.gz' i
 
 *-u, --use* 'TAG1'[,'TAG2']::
     specifies which tag to use in the query file ('TAG1') and the *-g* ('TAG2') file.
-    By default, the PL tag is used in the query file and GT in the *-g* file when
-    available.
+    By default, the PL tag is used in the query file and, when available, the GT tags in the
+    *-g* file.
+    {nbsp} +
+    {nbsp} +
+    Note that when the requested tag is not available, the program will attempt to use
+    the other tag. The output includes the number of sites that were matched by the four
+    possible mode (for example GT-vs-GT or GT-vs-PL).
 
 *Examples:*
 ----

test/gtcheck.5.1.out

@@ -5,4 +5,8 @@ INFO	sites-skipped-monoallelic	1
 INFO	sites-skipped-no-data	1
 INFO	sites-skipped-GT-not-diploid	1
 INFO	sites-skipped-PL-not-diploid	1
+INFO	sites-used-PL-vs-PL	0
+INFO	sites-used-PL-vs-GT	1
+INFO	sites-used-GT-vs-PL	0
+INFO	sites-used-GT-vs-GT	1
 DC	A	A	3.000150e-04	4.605170e+01	2

vcfgtcheck.c

@@ -1,6 +1,6 @@
 /*  vcfgtcheck.c -- Check sample identity.
 
-    Copyright (C) 2013-2021 Genome Research Ltd.
+    Copyright (C) 2013-2023 Genome Research Ltd.
 
     Author: Petr Danecek <[email protected]>
 
@@ -59,6 +59,7 @@ typedef struct
     int argc, gt_samples_is_file, qry_samples_is_file, regions_is_file, targets_is_file, pair_samples_is_file;
     int regions_overlap, targets_overlap;
     int qry_use_GT,gt_use_GT, nqry_smpl,ngt_smpl, *qry_smpl,*gt_smpl;
+    int nused[2][2];
     double *pdiff, *qry_prob, *gt_prob;
     uint32_t *ndiff,*ncnt,ncmp, npairs;
     int32_t *qry_arr,*gt_arr, nqry_arr,ngt_arr;
@@ -309,7 +310,7 @@ static void init_data(args_t *args)
         init_samples(args->qry_samples, args->qry_samples_is_file, &args->qry_smpl, &args->nqry_smpl, args->qry_hdr, args->qry_fname);
     }
     if ( args->gt_samples )
-    {   
+    {
         init_samples(args->gt_samples, args->gt_samples_is_file, &args->gt_smpl, &args->ngt_smpl,
             args->gt_hdr ? args->gt_hdr : args->qry_hdr,
             args->gt_fname ? args->gt_fname : args->qry_fname);
@@ -377,7 +378,7 @@ static void init_data(args_t *args)
         args->gt_prob  = args->cross_check ? args->qry_prob : (double*) malloc(3*args->ngt_smpl*sizeof(*args->gt_prob));
 
         // dsg2prob: the first index is bitmask of 8 possible dsg combinations (only 1<<0,1<<2,1<<3 are set, accessing
-        // anything else indicated an error, this is just to reuse gt_to_dsg()); the second index are the corresponding 
+        // anything else indicated an error, this is just to reuse gt_to_dsg()); the second index are the corresponding
         // probabilities of 0/0, 0/1, and 1/1 genotypes
         for (i=0; i<8; i++)
             for (j=0; j<3; j++)
@@ -555,7 +556,9 @@ static void process_line(args_t *args)
         args->gt_arr = args->qry_arr;
     }
 
+    // stats: number of compared sites, and used tags
     args->ncmp++;
+    args->nused[qry_use_GT][gt_use_GT]++;
 
     double af,hwe_dsg[8];
     if ( args->calc_hwe_prob )
@@ -636,7 +639,7 @@ static void process_line(args_t *args)
                 gt_dsg = gt_use_GT ? gt_to_prob(args,ptr,gt_prob) : pl_to_prob(args,ptr,gt_prob);
                 if ( !gt_dsg ) continue;                        // missing value
                 if ( args->hom_only && !(gt_dsg&5) ) continue;  // not a hom
-               
+
                 ptr = args->qry_arr + args->pairs[i].iqry*nqry1;
                 qry_dsg = qry_use_GT ? gt_to_prob(args,ptr,qry_prob) : pl_to_prob(args,ptr,qry_prob);
                 if ( !qry_dsg ) continue;                       // missing value
@@ -797,11 +800,15 @@ static void report(args_t *args)
     fprintf(args->fp,"INFO\tsites-skipped-no-data\t%u\n",args->nskip_no_data);
     fprintf(args->fp,"INFO\tsites-skipped-GT-not-diploid\t%u\n",args->nskip_dip_GT);
     fprintf(args->fp,"INFO\tsites-skipped-PL-not-diploid\t%u\n",args->nskip_dip_PL);
+    fprintf(args->fp,"INFO\tsites-used-PL-vs-PL\t%u\n",args->nused[0][0]);
+    fprintf(args->fp,"INFO\tsites-used-PL-vs-GT\t%u\n",args->nused[0][1]);
+    fprintf(args->fp,"INFO\tsites-used-GT-vs-PL\t%u\n",args->nused[1][0]);
+    fprintf(args->fp,"INFO\tsites-used-GT-vs-GT\t%u\n",args->nused[1][1]);
     fprintf(args->fp,"# DC, discordance:\n");
     fprintf(args->fp,"#     - query sample\n");
     fprintf(args->fp,"#     - genotyped sample\n");
-    fprintf(args->fp,"#     - discordance (number of mismatches; smaller is better)\n");
-    fprintf(args->fp,"#     - negative log of HWE probability at matching sites (rare genotypes mataches are more informative, bigger is better)\n");
+    fprintf(args->fp,"#     - discordance (either an abstract score or number of mismatches, see -e/-u in the man page for details; smaller is better)\n");
+    fprintf(args->fp,"#     - negative log of HWE probability at matching sites (rare genotypes matches are more informative, bigger is better)\n");
     fprintf(args->fp,"#     - number of sites compared (bigger is better)\n");
     fprintf(args->fp,"#DC\t[2]Query Sample\t[3]Genotyped Sample\t[4]Discordance\t[5]-log P(HWE)\t[6]Number of sites compared\n");
 
@@ -1023,7 +1030,7 @@ static int is_input_okay(args_t *args, int nmatch)
     return 1;
 
 not_okay:
-    fprintf(stderr,"INFO: skipping %s:%"PRIhts_pos", %s. (This is printed only once.)\n", 
+    fprintf(stderr,"INFO: skipping %s:%"PRIhts_pos", %s. (This is printed only once.)\n",
         bcf_seqname(hdr,rec),rec->pos+1,msg);
     return 0;
 }
@@ -1097,7 +1104,7 @@ int main_vcfgtcheck(int argc, char *argv[])
     args->es_max_mem = strdup("500M");
 
     // In simulated sample swaps the minimum error was 0.3 and maximum intra-sample error was 0.23
-    //    - min_inter: pairs with smaller err value will be considered identical 
+    //    - min_inter: pairs with smaller err value will be considered identical
     //    - max_intra: pairs with err value bigger than abs(max_intra_err) will be considered
     //                  different. If negative, the cutoff may be heuristically lowered
     args->min_inter_err =  0.23;
@@ -1169,7 +1176,7 @@ int main_vcfgtcheck(int argc, char *argv[])
             case 3 : args->calc_hwe_prob = 0; break;
             case 4 : error("The option -S, --target-sample has been deprecated\n"); break;
             case 5 : args->dry_run = 1; break;
-            case 6 : 
+            case 6 :
                 args->distinctive_sites = strtod(optarg,&tmp);
                 if ( *tmp )
                 {
@@ -1202,7 +1209,7 @@ int main_vcfgtcheck(int argc, char *argv[])
                 else if ( !strncasecmp("qry:",optarg,4) ) args->qry_samples = optarg+4;
                 else error("Which one? Query samples (qry:%s) or genotype samples (gt:%s)?\n",optarg,optarg);
                 break;
-            case 'S': 
+            case 'S':
                 if ( !strncasecmp("gt:",optarg,3) ) args->gt_samples = optarg+3, args->gt_samples_is_file = 1;
                 else if ( !strncasecmp("qry:",optarg,4) ) args->qry_samples = optarg+4, args->qry_samples_is_file = 1;
                 else error("Which one? Query samples (qry:%s) or genotype samples (gt:%s)?\n",optarg,optarg);