Make +check-ploidy optionally use missing genotypes. Replaces #1531

plugins/check-ploidy.c

@@ -1,5 +1,5 @@
 /* 
-    Copyright (C) 2017 Genome Research Ltd.
+    Copyright (C) 2017-2021 Genome Research Ltd.
 
     Author: Petr Danecek <[email protected]>
 
@@ -47,7 +47,7 @@ typedef struct
 {
     int argc;
     char **argv;
-    int rid, gt_id, ndat;
+    int rid, gt_id, ndat, ignore_missing;
     dat_t *dat;
     bcf_hdr_t *hdr;
 }
@@ -69,8 +69,15 @@ const char *usage(void)
         "Options:\n"
         "   run \"bcftools plugin\" for a list of common options\n"
         "\n"
+        "Plugin options:\n"
+        "   -m, --use-missing           use missing and half-missing genotypes such as ., ./., 0/1\n"
+        "\n"
         "Example:\n"
+        "   # report ploidy, ignore missing genotypes\n"
         "   bcftools +check-ploidy file.bcf\n"
+        "\n"
+        "   # use missing genotypes\n"
+        "   bcftools +check-ploidy file.bcf -- -m\n"
         "\n";
 }
 
@@ -79,6 +86,25 @@ int init(int argc, char **argv, bcf_hdr_t *in, bcf_hdr_t *out)
     args = (args_t*) calloc(1,sizeof(args_t));
     args->argc = argc; args->argv = argv;
     if ( !in ) error("%s",usage());
+
+    args->ignore_missing = 1;
+    static struct option loptions[] =
+    {
+        {"use-missing",0,0,'m'},
+        {0,0,0,0}
+    };
+    int c;
+    while ((c = getopt_long(argc, argv, "?hm",loptions,NULL)) >= 0)
+    {
+        switch (c)
+        {
+            case 'm': args->ignore_missing = 0; break;
+            case 'h':
+            case '?':
+            default: error("%s", usage()); break;
+        }
+    }
+
     args->hdr  = in;
     args->ndat = bcf_hdr_nsamples(args->hdr);
     args->dat  = (dat_t*) calloc(args->ndat,sizeof(dat_t));
@@ -124,7 +150,7 @@ bcf1_t *process(bcf1_t *rec)
             for (nal=0; nal<fmt_gt->n; nal++) \
             { \
                 if ( p[nal]==vector_end ) break; /* smaller ploidy */ \
-                if ( bcf_gt_is_missing(p[nal]) ) { missing=1; break; } /* missing allele */ \
+                if ( bcf_gt_is_missing(p[nal]) && args->ignore_missing ) { missing=1; break; } /* missing allele */ \
             } \
             if ( !nal || missing ) continue; /* missing genotype */ \
             dat_t *dat = &args->dat[i]; \

test/checkploidy.2.out

@@ -0,0 +1,2 @@
+# [1]Sample	[2]Chromosome	[3]Region Start	[4]Region End	[5]Ploidy
+S1	20	1	1	2

test/checkploidy.2.vcf

@@ -0,0 +1,11 @@
+##fileformat=VCFv4.1
+##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
+##contig=<ID=20,assembly=b37>
+##contig=<ID=21,assembly=b37>
+##contig=<ID=22,assembly=b37>
+##contig=<ID=X,assembly=b37>
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	S1
+20	1	1	C	G	.	PASS	.	GT	0/0
+21	1	2	C	G	.	PASS	.	GT	./0
+22	1	3	C	G	.	PASS	.	GT	./.
+X	1	4	C	G	.	PASS	.	GT	.

test/checkploidy.3.out

@@ -0,0 +1,5 @@
+# [1]Sample	[2]Chromosome	[3]Region Start	[4]Region End	[5]Ploidy
+S1	20	1	1	2
+S1	21	1	1	2
+S1	22	1	1	2
+S1	X	1	1	1

test/test.pl

@@ -436,6 +436,8 @@
 test_vcf_annotate($opts,in=>'annotate22',vcf=>'annotate22',out=>'annotate30.out',args=>'-c FMT/XX,INFO/XX -x FMT/XX,INFO/XX');
 test_vcf_annotate($opts,in=>'annotate23',tab=>'annotate23',out=>'annotate31.out',args=>'-c CHROM,POS,~ID,REF,ALT,INFO/END');
 test_vcf_plugin($opts,in=>'checkploidy',out=>'checkploidy.out',cmd=>'+check-ploidy --no-version');
+test_vcf_plugin($opts,in=>'checkploidy.2',out=>'checkploidy.2.out',cmd=>'+check-ploidy --no-version');
+test_vcf_plugin($opts,in=>'checkploidy.2',out=>'checkploidy.3.out',cmd=>'+check-ploidy --no-version',args=>'-- -m');
 test_vcf_plugin($opts,in=>'plugin1',out=>'missing2ref.out',cmd=>'+missing2ref --no-version');
 test_vcf_plugin($opts,in=>'plugin1',out=>'missing2ref.out',cmd=>'+setGT --no-version',args=>'-- -t . -n 0');
 test_vcf_plugin($opts,in=>'setGT',out=>'setGT.1.out',cmd=>'+setGT --no-version',args=>'-- -t q -n 0 -i \'GT~"." && FMT/DP=30 && GQ=150\'');