Improvements and changes in gtcheck

pd3 · pd3 · commit 983e610ab5d6 · 2023-11-03T09:32:58.000Z
- Add output options `-o, --output` and `-O, --output-type`

- Add filtering options `-i, --include` and `-e, --exclude`

- Rename the short option `-e, --error-probability` from lower case to upper
  case `-E, --error-probability`. This is a backward compatible change, the
  program can tell when the lower case `-e` is used in the previous meaning

- Changes to the output format, replace the DC section with DCv2:

    - adds a new column for the number of matching genotypes

    - fixes in HWE score calculation plus output average HWE score rather
      than absolute HWE score

    - better description of fields
diff --git a/Makefile b/Makefile
@@ -245,7 +245,7 @@ vcfcall.o: vcfcall.c $(htslib_vcf_h) $(htslib_kfunc_h) $(htslib_synced_bcf_reade
 vcfconcat.o: vcfconcat.c $(htslib_vcf_h) $(htslib_synced_bcf_reader_h) $(htslib_kseq_h) $(htslib_bgzf_h) $(htslib_tbx_h) $(htslib_thread_pool_h) $(bcftools_h)
 vcfconvert.o: vcfconvert.c $(htslib_faidx_h) $(htslib_vcf_h) $(htslib_bgzf_h) $(htslib_synced_bcf_reader_h) $(htslib_vcfutils_h) $(htslib_kseq_h) $(bcftools_h) $(filter_h) $(convert_h) $(tsv2vcf_h)
 vcffilter.o: vcffilter.c $(htslib_vcf_h) $(htslib_synced_bcf_reader_h) $(htslib_vcfutils_h) $(bcftools_h) $(filter_h) rbuf.h regidx.h
-vcfgtcheck.o: vcfgtcheck.c $(htslib_vcf_h) $(htslib_synced_bcf_reader_h) $(htslib_vcfutils_h) $(htslib_kbitset_h) $(htslib_hts_os_h) $(bcftools_h) extsort.h
+vcfgtcheck.o: vcfgtcheck.c $(htslib_vcf_h) $(htslib_synced_bcf_reader_h) $(htslib_vcfutils_h) $(htslib_kbitset_h) $(htslib_hts_os_h) $(htslib_bgzf_h) $(bcftools_h) extsort.h filter.h
 vcfindex.o: vcfindex.c $(htslib_vcf_h) $(htslib_tbx_h) $(htslib_kstring_h) $(htslib_bgzf_h) $(bcftools_h)
 vcfisec.o: vcfisec.c $(htslib_vcf_h) $(htslib_synced_bcf_reader_h) $(htslib_vcfutils_h) $(htslib_hts_os_h) $(bcftools_h) $(filter_h)
 vcfmerge.o: vcfmerge.c $(htslib_vcf_h) $(htslib_synced_bcf_reader_h) $(htslib_vcfutils_h) $(htslib_faidx_h) regidx.h $(bcftools_h) vcmp.h $(htslib_khash_h)
diff --git a/NEWS b/NEWS
@@ -28,6 +28,25 @@ Changes affecting specific commands:
     - New `-*, --keep-unseen-allele` option to output the unobserved allele <*>,
       intended for gVCF.
 
+* bcftools gtcheck
+
+    - Add output options `-o, --output` and `-O, --output-type`
+
+    - Add filtering options `-i, --include` and `-e, --exclude`
+
+    - Rename the short option `-e, --error-probability` from lower case to upper
+      case `-E, --error-probability`
+
+    - Changes to the output format, replace the DC section with DCv2:
+
+        - adds a new column for the number of matching genotypes
+
+        - fixes in HWE score calculation plus output average HWE score rather
+          than absolute HWE score
+
+        - better description of fields
+
+
 * bcftools mpileup
 
     - Output MIN_DP rather than MinDP in gVCF mode
diff --git a/doc/bcftools.txt b/doc/bcftools.txt
@@ -1640,25 +1640,38 @@ The discordance score can be interpreted as the number of mismatching genotypes
 *--dry-run*::
     Stop after first record to estimate required time.
 
-*-e, --error-probability* 'INT'::
+*-e, --exclude* ['qry'|'gt']:'EXPRESSION'::
+    Exclude sites from query file ('qry:') or genotype file ('gt:') for which 'EXPRESSION' is true.
+    For valid expressions see *<<expressions,EXPRESSIONS>>*.
+
+*-E, --error-probability* 'INT'::
     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).
     {nbsp} +
     {nbsp} +
-    If *-e* is set to 0, the discordance score can be interpreted as the number of mismatching genotypes,
+    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.
+    If performance is an issue, set *-E 0* for faster run times but less accurate results.
+    {nbsp} +
+    {nbsp} +
+    Note that in previous versions of bcftools (<=1.18), this option used to be a smaller case *-e*. It
+    changed to make room for the filtering option *-e, --exclude* to stay consistent across other
+    commands.
 
 *-g, --genotypes* 'FILE'::
     VCF/BCF file with reference genotypes to compare against
 
 *-H, --homs-only*::
     Homozygous genotypes only, useful with low coverage data (requires *-g, --genotypes*)
 
+*-i, --include* ['qry'|'gt']:'EXPRESSION'::
+    Include sites from query file ('qry:') or genotype file ('gt:') for which 'EXPRESSION' is true.
+    For valid expressions see *<<expressions,EXPRESSIONS>>*.
+
 *--n-matches* 'INT'::
     Print only top INT matches for each sample, 0 for unlimited. Use negative value
     to sort by HWE probability rather than the number of discordant sites. Note
@@ -1668,6 +1681,12 @@ The discordance score can be interpreted as the number of mismatching genotypes
     Disable calculation of HWE probability to reduce memory requirements with
     comparisons between very large number of sample pairs.
 
+*-o, --output* 'FILE'::
+    Write to 'FILE' rather than to standard output, where it is written by default.
+
+*-O, --output-type* 't'|'z'::
+    Write a plain ('t') or compressed ('z') text tab-delimited output.
+
 *-p, --pairs* 'LIST'::
     A comma-separated list of sample pairs to compare. When the *-g* option is given, the first
     sample must be from the query file, the second from the *-g* file, third from the query file
diff --git a/test/gtcheck.1.out b/test/gtcheck.1.out
@@ -1 +1 @@
-DC	s1	s1	0	2.371900e+01	2
+DCv2	s1	s1	0	3.465736e-01	2	2
diff --git a/test/gtcheck.10.out b/test/gtcheck.10.out
@@ -1 +1 @@
-DC	s1	s1	4.002001e-03	2.371900e+01	2
+DCv2	s1	s1	4.002001e-03	3.465736e-01	2	2
diff --git a/test/gtcheck.11.out b/test/gtcheck.11.out
@@ -1 +1 @@
-DC	s1	s1	0.000000e+00	2.371900e+01	2
+DCv2	s1	s1	0.000000e+00	3.465736e-01	2	2
diff --git a/test/gtcheck.12.out b/test/gtcheck.12.out
@@ -1,10 +1,10 @@
-DC	B	A	5.733631e-01	2.253795e+00	2
-DC	C	A	4.938053e+00	8.675006e-01	2
-DC	C	B	2.791391e+00	8.675006e-01	2
-DC	D	A	5.022610e+00	0.000000e+00	2
-DC	D	B	5.178533e+00	0.000000e+00	2
-DC	D	C	4.938053e+00	0.000000e+00	2
-DC	E	A	7.325195e+00	0.000000e+00	2
-DC	E	B	5.178533e+00	0.000000e+00	2
-DC	E	C	2.635468e+00	1.386294e+00	2
-DC	E	D	2.720025e+00	2.407946e+00	2
+DCv2	B	A	5.733631e-01	1.126897e+00	2	2
+DCv2	C	A	4.938053e+00	4.337503e-01	2	2
+DCv2	C	B	2.791391e+00	4.337503e-01	2	2
+DCv2	D	A	5.022610e+00	0.000000e+00	2	2
+DCv2	D	B	5.178533e+00	0.000000e+00	2	2
+DCv2	D	C	4.938053e+00	0.000000e+00	2	2
+DCv2	E	A	7.325195e+00	0.000000e+00	2	2
+DCv2	E	B	5.178533e+00	0.000000e+00	2	2
+DCv2	E	C	2.635468e+00	6.931472e-01	2	2
+DCv2	E	D	2.720025e+00	3.566749e-01	2	2
diff --git a/test/gtcheck.2.out b/test/gtcheck.2.out
@@ -1 +1 @@
-DC	s1	s1	0	6.931472e-01	1
+DCv2	s1	s1	0	6.931472e-01	1	1
diff --git a/test/gtcheck.3.out b/test/gtcheck.3.out
@@ -1,10 +1,10 @@
-DC	B	A	0	2.253795e+00	2
-DC	C	A	1	8.675006e-01	2
-DC	C	B	1	8.675006e-01	2
-DC	D	A	2	0.000000e+00	2
-DC	D	B	2	0.000000e+00	2
-DC	D	C	2	0.000000e+00	2
-DC	E	A	2	0.000000e+00	2
-DC	E	B	2	0.000000e+00	2
-DC	E	C	1	1.386294e+00	2
-DC	E	D	1	2.407946e+00	2
+DCv2	B	A	0	1.126897e+00	2	2
+DCv2	C	A	1	8.675006e-01	2	1
+DCv2	C	B	1	8.675006e-01	2	1
+DCv2	D	A	2	0.000000e+00	2	0
+DCv2	D	B	2	0.000000e+00	2	0
+DCv2	D	C	2	0.000000e+00	2	0
+DCv2	E	A	2	0.000000e+00	2	0
+DCv2	E	B	2	0.000000e+00	2	0
+DCv2	E	C	1	1.386294e+00	2	1
+DCv2	E	D	1	7.133499e-01	2	1
diff --git a/test/gtcheck.4.out b/test/gtcheck.4.out
@@ -1,3 +1,3 @@
-DC	D	C	2	0.000000e+00	2
-DC	E	C	1	1.386294e+00	2
-DC	E	D	1	2.407946e+00	2
+DCv2	D	C	2	0.000000e+00	2	0
+DCv2	E	C	1	1.386294e+00	2	1
+DCv2	E	D	1	7.133499e-01	2	1
diff --git a/test/gtcheck.5.1.out b/test/gtcheck.5.1.out
@@ -5,8 +5,9 @@ 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-skipped-filtering-expression	0
 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
+DCv2	A	A	3.000150e-04	0.000000e+00	2	2
diff --git a/test/gtcheck.5.out b/test/gtcheck.5.out
@@ -1 +1 @@
-DC	B	D	2	0.000000e+00	2
+DCv2	B	D	2	0.000000e+00	2	0
diff --git a/test/gtcheck.6.1.out b/test/gtcheck.6.1.out
@@ -0,0 +1,2 @@
+DCv2	A	B	2.000100e-04	5.753641e-01	1	1
+DCv2	B	C	9.903588e+00	0.000000e+00	1	1
diff --git a/test/gtcheck.6.out b/test/gtcheck.6.out
@@ -1,2 +1,2 @@
-DC	B	C	1	8.675006e-01	2
-DC	B	D	2	0.000000e+00	2
+DCv2	B	C	1	8.675006e-01	2	1
+DCv2	B	D	2	0.000000e+00	2	0
diff --git a/test/gtcheck.6.vcf b/test/gtcheck.6.vcf
@@ -0,0 +1,6 @@
+##fileformat=VCFv4.3
+##reference=ref.fa
+##contig=<ID=1,assembly=b37,length=107349540>
+##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	A	B	C	D
+1	1	.	C	T	.	PASS	.	GT	0/0	0/0	0/1	0/1
diff --git a/test/gtcheck.7.out b/test/gtcheck.7.out
@@ -1 +1 @@
-DC	s1	s1	0	2.302585e+01	1
+DCv2	s1	s1	0	0.000000e+00	1	1
diff --git a/test/gtcheck.8.out b/test/gtcheck.8.out
@@ -1,6 +1,6 @@
-DC	A	B	3	1.734223e+01	9
-DC	C	D	6	1.075056e+01	9
-DC	E	F	9	0.000000e+00	9
+DCv2	A	B	3	5.877867e-01	9	6
+DCv2	C	D	6	3.646431e-01	9	3
+DCv2	E	F	9	0.000000e+00	9	0
 DS	1	3	3	0
 DS	2	3	3	1
 DS	3	3	3	2
diff --git a/test/gtcheck.9.out b/test/gtcheck.9.out
@@ -1,10 +1,10 @@
-DC	B	A	0	2.253795e+00	2
-DC	C	A	1	8.675006e-01	2
-DC	C	B	1	8.675006e-01	2
-DC	D	A	2	0.000000e+00	2
-DC	D	B	2	0.000000e+00	2
-DC	D	C	2	0.000000e+00	2
-DC	E	A	2	0.000000e+00	2
-DC	E	B	2	0.000000e+00	2
-DC	E	C	1	1.386294e+00	2
-DC	E	D	1	2.407946e+00	2
+DCv2	B	A	0	1.126897e+00	2	2
+DCv2	C	A	1	8.675006e-01	2	1
+DCv2	C	B	1	8.675006e-01	2	1
+DCv2	D	A	2	0.000000e+00	2	0
+DCv2	D	B	2	0.000000e+00	2	0
+DCv2	D	C	2	0.000000e+00	2	0
+DCv2	E	A	2	0.000000e+00	2	0
+DCv2	E	B	2	0.000000e+00	2	0
+DCv2	E	C	1	1.386294e+00	2	1
+DCv2	E	D	1	7.133499e-01	2	1
diff --git a/test/gtcheck.ntop.1.out b/test/gtcheck.ntop.1.out
@@ -1,6 +1,6 @@
-DC	smpl	x1	4.951814e+01	2.197225e+00	6
-DC	smpl	x2	9.904588e+00	1.075056e+01	6
-DC	smpl	x3	1.000050e-03	1.075056e+01	5
-DC	smpl	x4	2.971136e+01	7.613325e+00	6
-DC	smpl	x5	1.200060e-03	1.258314e+01	6
-DC	smpl	x6	3.961475e+01	5.416100e+00	6
+DCv2	smpl	x1	4.951814e+01	1.351550e-01	6	6
+DCv2	smpl	x2	9.904588e+00	7.931820e-01	6	6
+DCv2	smpl	x3	1.000050e-03	9.518185e-01	5	5
+DCv2	smpl	x4	2.971136e+01	3.662041e-01	6	6
+DCv2	smpl	x5	1.200060e-03	9.634573e-01	6	6
+DCv2	smpl	x6	3.961475e+01	2.310491e-01	6	6
diff --git a/test/gtcheck.ntop.2.out b/test/gtcheck.ntop.2.out
@@ -1,2 +1,2 @@
-DC	smpl	x5	1.200060e-03	1.258314e+01	6
-DC	smpl	x3	1.000050e-03	1.075056e+01	5
+DCv2	smpl	x5	1.200060e-03	9.634573e-01	6	6
+DCv2	smpl	x3	1.000050e-03	9.518185e-01	5	5
diff --git a/test/test.pl b/test/test.pl
@@ -917,6 +917,7 @@
 run_test(\&test_gtcheck,$opts,in=>'gtcheck.ntop',gts=>'gtcheck.ntop.gts',out=>'gtcheck.ntop.1.out',args=>q[]);
 run_test(\&test_gtcheck,$opts,in=>'gtcheck.ntop',gts=>'gtcheck.ntop.gts',out=>'gtcheck.ntop.2.out',args=>q[--n-matches 2]);
 run_test(\&test_gtcheck,$opts,in=>'gtcheck.5',gts=>'gtcheck.5.gts',out=>'gtcheck.5.1.out',args=>q[],grep=>'grep -v Time');
+run_test(\&test_gtcheck,$opts,in=>'gtcheck.6',out=>'gtcheck.6.1.out',args=>q[-p A,B,B,C]);
 
 print "\nNumber of tests:\n";
 printf "    total   .. %d\n", $$opts{nok}+$$opts{nfailed};
diff --git a/vcfgtcheck.c b/vcfgtcheck.c

Original file line number	Diff line number	Diff line change
`@@ -1 +1 @@`
`1`		`-DC s1 s1 0 2.371900e+01 2`
	`1`	`+DCv2 s1 s1 0 3.465736e-01 2 2`
Original file line number	Diff line number	Diff line change
`@@ -1 +1 @@`
`1`		`-DC s1 s1 4.002001e-03 2.371900e+01 2`
	`1`	`+DCv2 s1 s1 4.002001e-03 3.465736e-01 2 2`
Original file line number	Diff line number	Diff line change
`@@ -1 +1 @@`
`1`		`-DC s1 s1 0.000000e+00 2.371900e+01 2`
	`1`	`+DCv2 s1 s1 0.000000e+00 3.465736e-01 2 2`
Original file line number	Diff line number	Diff line change
`@@ -1 +1 @@`
`1`		`-DC s1 s1 0 6.931472e-01 1`
	`1`	`+DCv2 s1 s1 0 6.931472e-01 1 1`