Preserve phasing when mergin half-missing genotypes

Fixes #2331
samtools · Dec 11, 2024 · 6559a12 · 6559a12
1 parent 313c597
commit 6559a12
Show file tree

Hide file tree

Showing 6 changed files with 34 additions and 2 deletions.
diff --git a/NEWS b/NEWS
@@ -21,6 +21,10 @@ Changes affecting specific commands:
     - Check the input GFF for features outside transcript boundaries and extend the transcript
       to contain the feature fully (#2323)
 
+* bcftools merge
+
+    - Preserve phasing in half-missing genotypes (#2331)
+
 * bcftools query
 
     - The functions used in -i/-e filtering expressions (such as SUM, MEDIAN, etc) can be

diff --git a/test/merge.phased.1.1.out b/test/merge.phased.1.1.out
@@ -0,0 +1,7 @@
+##fileformat=VCFv4.2
+##FILTER=<ID=PASS,Description="All filters passed">
+##reference=ref.fa
+##contig=<ID=1,length=307041717>
+##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	a	b	c	d
+1	5695120	.	C	T	.	.	.	GT	0|.	.|0	1|.	.|1
diff --git a/test/merge.phased.1.a.vcf b/test/merge.phased.1.a.vcf
@@ -0,0 +1,6 @@
+##fileformat=VCFv4.2
+##reference=ref.fa
+##contig=<ID=1,length=307041717>
+##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	a	b
+1	5695120	.	C	T	.	.	.	GT	0|.	.|0
diff --git a/test/merge.phased.1.b.vcf b/test/merge.phased.1.b.vcf
@@ -0,0 +1,6 @@
+##fileformat=VCFv4.2
+##reference=ref.fa
+##contig=<ID=1,length=307041717>
+##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	c	d
+1	5695120	.	C	T	.	.	.	GT	1|.	.|1
diff --git a/test/test.pl b/test/test.pl
@@ -61,6 +61,7 @@
 run_test(\&test_vcf_isec,$opts,in=>['isec-miss.1.1','isec-miss.1.2','isec-miss.1.3'],out=>'isec-miss.1.1.out',args=>'-R {PATH}/isec-miss.1.regs.txt -n +1');
 run_test(\&test_vcf_isec,$opts,in=>['isec-miss.2.1','isec-miss.2.2','isec-miss.2.3'],out=>'isec-miss.2.1.out',args=>'-n +1 -r 20:100,20:140,12:55,20:140,20:100');
 run_test(\&test_vcf_isec,$opts,in=>['isec-miss.2.1','isec-miss.2.2','isec-miss.2.3'],out=>'isec-miss.2.1.out',args=>'-R {PATH}/isec-miss.1.regs.txt -n +1');
+run_test(\&test_vcf_merge,$opts,in=>['merge.phased.1.a','merge.phased.1.b'],out=>'merge.phased.1.1.out',args=>'');
 run_test(\&test_vcf_merge,$opts,in=>['merge.11.a','merge.11.b'],out=>'merge.11.1.out',args=>'');
 run_test(\&test_vcf_merge,$opts,in=>['merge.join.a','merge.join.b'],out=>'merge.join.1.out',args=>'-i AF:join');
 run_test(\&test_vcf_merge,$opts,in=>['merge.LPL.a'],out=>'merge.LPL.0.out',args=>'--force-single');

diff --git a/vcfmerge.c b/vcfmerge.c
@@ -1759,7 +1759,11 @@ void merge_GT(args_t *args, bcf_fmt_t **fmt_map, bcf1_t *out)
                         type_t val = convert(&p_ori[k * sizeof(type_t)]); \
                         if ( val==vector_end ) break; /* smaller ploidy */ \
                         ma->smpl_ploidy[ismpl+j]++; \
-                        if ( bcf_gt_is_missing(val) ) tmp[k] = 0; /* missing allele */ \
+                        if ( bcf_gt_is_missing(val) ) \
+                        { \
+                            if ( bcf_gt_is_phased(val) ) tmp[k] = 1; /* missing allele, phased */ \
+                            else tmp[k] = 0; /* missing allele, unphased */ \
+                        } \
                         else tmp[k] = val; \
                     } \
                     for (; k<nsize; k++) tmp[k] = bcf_int32_vector_end; \
@@ -1777,7 +1781,11 @@ void merge_GT(args_t *args, bcf_fmt_t **fmt_map, bcf1_t *out)
                     type_t val = convert(&p_ori[k * sizeof(type_t)]); \
                     if ( val==vector_end ) break; /* smaller ploidy */ \
                     ma->smpl_ploidy[ismpl+j]++; \
-                    if ( bcf_gt_is_missing(val) ) tmp[k] = 0; /* missing allele */ \
+                    if ( bcf_gt_is_missing(val) ) \
+                    { \
+                        if ( bcf_gt_is_phased(val) ) tmp[k] = 1; /* missing allele, phased */ \
+                        else tmp[k] = 0; /* missing allele, unphased */ \
+                    } \
                     else \
                     { \
                         int al = (val>>1) - 1; \