Skip to content

Commit

Permalink
Revamped line matching code to fix problems in gVCF merging
Browse files Browse the repository at this point in the history
Previously gVCF blocks were correctly split but in some cases FORMAT
fields would not be updated at split sites.

This commit revamps the line matching code, adds comments and makes it
more understandable.

Fixes #1891 and revisits #1164
  • Loading branch information
pd3 committed Apr 4, 2023
1 parent 796594f commit 60c65eb
Show file tree
Hide file tree
Showing 19 changed files with 346 additions and 177 deletions.
3 changes: 3 additions & 0 deletions NEWS
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,9 @@ Changes affecting specific commands:
- Use values pertaining to the unknown allele (<*> or <NON_REF>) when available
to prevent mixtures of known and missing values (#1888)

- Revamped line matching code to fix problems in gVCF merging where split gVCF blocks
would not update genotypes (#1891, #1164).

* bcftools norm

- The `-m, --multiallelics +` mode now preserves phasing (#1893)
Expand Down
19 changes: 10 additions & 9 deletions test/merge.10.1.out
Original file line number Diff line number Diff line change
@@ -1,13 +1,14 @@
##fileformat=VCFv4.3
##FILTER=<ID=PASS,Description="All filters passed">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##contig=<ID=1,assembly=b37,length=249250621>
##reference=file:ref.fa
#CHROM POS ID REF ALT QUAL FILTER INFO
1 3000000 . C CCG . . .
1 3000000 . C CAA . . .
1 3000150 . CC C . . .
1 3000150 . C CTT . . .
1 3000152 . C A . . .
1 3000152 . C CA . . .
1 3000154 . C A . . .
1 3000154 . C T . . .
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT a b
1 3000000 . C CCG . . . GT 0/1 ./.
1 3000000 . C CAA . . . GT ./. 0/1
1 3000150 . CC C . . . GT 0/1 ./.
1 3000150 . C CTT . . . GT ./. 0/1
1 3000152 . C A . . . GT 0/1 ./.
1 3000152 . C CA . . . GT ./. 0/1
1 3000154 . C A . . . GT 0/1 ./.
1 3000154 . C T . . . GT ./. 0/1
13 changes: 7 additions & 6 deletions test/merge.10.2.out
Original file line number Diff line number Diff line change
@@ -1,10 +1,11 @@
##fileformat=VCFv4.3
##FILTER=<ID=PASS,Description="All filters passed">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##contig=<ID=1,assembly=b37,length=249250621>
##reference=file:ref.fa
#CHROM POS ID REF ALT QUAL FILTER INFO
1 3000000 . C CCG,CAA . . .
1 3000150 . CC C,CTTC . . .
1 3000152 . C A . . .
1 3000152 . C CA . . .
1 3000154 . C A,T . . .
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT a b
1 3000000 . C CCG,CAA . . . GT 0/1 0/2
1 3000150 . CC C,CTTC . . . GT 0/1 0/2
1 3000152 . C A . . . GT 0/1 ./.
1 3000152 . C CA . . . GT ./. 0/1
1 3000154 . C A,T . . . GT 0/1 0/2
15 changes: 8 additions & 7 deletions test/merge.10.3.out
Original file line number Diff line number Diff line change
@@ -1,11 +1,12 @@
##fileformat=VCFv4.3
##FILTER=<ID=PASS,Description="All filters passed">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##contig=<ID=1,assembly=b37,length=249250621>
##reference=file:ref.fa
#CHROM POS ID REF ALT QUAL FILTER INFO
1 3000000 . C CCG,CAA . . .
1 3000150 . CC C,CTTC . . .
1 3000152 . C A . . .
1 3000152 . C CA . . .
1 3000154 . C A . . .
1 3000154 . C T . . .
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT a b
1 3000000 . C CCG,CAA . . . GT 0/1 0/2
1 3000150 . CC C . . . GT 0/1 ./.
1 3000150 . C CTT . . . GT ./. 0/1
1 3000152 . C A . . . GT 0/1 ./.
1 3000152 . C CA . . . GT ./. 0/1
1 3000154 . C A,T . . . GT 0/1 0/2
11 changes: 6 additions & 5 deletions test/merge.10.a.vcf
Original file line number Diff line number Diff line change
@@ -1,8 +1,9 @@
##fileformat=VCFv4.3
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##contig=<ID=1,assembly=b37,length=249250621>
##reference=file:ref.fa
#CHROM POS ID REF ALT QUAL FILTER INFO
1 3000000 . C CCG . . .
1 3000150 . CC C . . .
1 3000152 . C A . . .
1 3000154 . C A . . .
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT a
1 3000000 . C CCG . . . GT 0/1
1 3000150 . CC C . . . GT 0/1
1 3000152 . C A . . . GT 0/1
1 3000154 . C A . . . GT 0/1
11 changes: 6 additions & 5 deletions test/merge.10.b.vcf
Original file line number Diff line number Diff line change
@@ -1,8 +1,9 @@
##fileformat=VCFv4.3
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##contig=<ID=1,assembly=b37,length=249250621>
##reference=file:ref.fa
#CHROM POS ID REF ALT QUAL FILTER INFO
1 3000000 . C CAA . . .
1 3000150 . C CTT . . .
1 3000152 . C CA . . .
1 3000154 . C T . . .
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT b
1 3000000 . C CAA . . . GT 0/1
1 3000150 . C CTT . . . GT 0/1
1 3000152 . C CA . . . GT 0/1
1 3000154 . C T . . . GT 0/1
9 changes: 7 additions & 2 deletions test/merge.gvcf.10.1.out
Original file line number Diff line number Diff line change
@@ -1,6 +1,11 @@
##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##contig=<ID=chr1,length=248956422>
##contig=<ID=chr2,length=248956422>
##INFO=<ID=END,Number=1,Type=Integer,Description="Stop position of the interval">
#CHROM POS ID REF ALT QUAL FILTER INFO
chr1 1 . A <*>,C . . END=2
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT a b
chr1 1 . A <*>,C . . END=2 GT 0/0 0/2
chr2 1 . A <*> . . END=2 GT 0/0 ./.
chr2 3 . G C,<*>,A . . . GT 0/1 0/3
chr2 4 . T <*> . . END=6 GT 0/0 0/0
11 changes: 8 additions & 3 deletions test/merge.gvcf.10.2.out
Original file line number Diff line number Diff line change
@@ -1,7 +1,12 @@
##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##contig=<ID=chr1,length=248956422>
##contig=<ID=chr2,length=248956422>
##INFO=<ID=END,Number=1,Type=Integer,Description="Stop position of the interval">
#CHROM POS ID REF ALT QUAL FILTER INFO
chr1 1 . A <*> . . END=2
chr1 1 . A C . . .
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT a b
chr1 1 . A <*>,C . . END=2 GT 0/0 0/2
chr2 1 . A <*> . . END=2 GT 0/0 ./.
chr2 3 . G C,<*> . . . GT 0/1 ./.
chr2 3 . G A,<*> . . . GT ./. 0/1
chr2 4 . T <*> . . END=6 GT 0/0 0/0
12 changes: 9 additions & 3 deletions test/merge.gvcf.10.3.out
Original file line number Diff line number Diff line change
@@ -1,7 +1,13 @@
##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##contig=<ID=chr1,length=248956422>
##contig=<ID=chr2,length=248956422>
##INFO=<ID=END,Number=1,Type=Integer,Description="Stop position of the interval">
#CHROM POS ID REF ALT QUAL FILTER INFO
chr1 1 . A <*>,C . . .
chr1 2 . C <*> . . .
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT a b
chr1 1 . A <*>,C . . . GT 0/0 0/2
chr1 2 . C <*> . . . GT 0/0 ./.
chr2 1 . A <*> . . END=2 GT 0/0 ./.
chr2 3 . G C,<*>,A . . . GT 0/1 0/3
chr2 4 . T <*> . . END=6 GT 0/0 0/0
chr2 7 . G <*> . . END=8 GT ./. 0/0
13 changes: 10 additions & 3 deletions test/merge.gvcf.10.4.out
Original file line number Diff line number Diff line change
@@ -1,7 +1,14 @@
##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##contig=<ID=chr1,length=248956422>
##contig=<ID=chr2,length=248956422>
##INFO=<ID=END,Number=1,Type=Integer,Description="Stop position of the interval">
#CHROM POS ID REF ALT QUAL FILTER INFO
chr1 1 . A C . . .
chr1 2 . C <*> . . .
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT a b
chr1 1 . A <*>,C . . . GT 0/0 0/2
chr1 2 . C <*> . . . GT 0/0 ./.
chr2 1 . A <*> . . END=2 GT 0/0 ./.
chr2 3 . G C,<*> . . . GT 0/1 ./.
chr2 3 . G A,<*> . . . GT ./. 0/1
chr2 4 . T <*> . . END=6 GT 0/0 0/0
chr2 7 . G <*> . . END=8 GT ./. 0/0
9 changes: 7 additions & 2 deletions test/merge.gvcf.10.a.vcf
Original file line number Diff line number Diff line change
@@ -1,6 +1,11 @@
##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##contig=<ID=chr1,length=248956422>
##contig=<ID=chr2,length=248956422>
##INFO=<ID=END,Number=1,Type=Integer,Description="Stop position of the interval">
#CHROM POS ID REF ALT QUAL FILTER INFO
chr1 1 . A <*> . . END=2
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT a
chr1 1 . A <*> . . END=2 GT 0/0
chr2 1 . A <*> . . END=2 GT 0/0
chr2 3 . G C,<*> . . . GT 0/1
chr2 4 . T <*> . . END=6 GT 0/0
9 changes: 7 additions & 2 deletions test/merge.gvcf.10.b.vcf
Original file line number Diff line number Diff line change
@@ -1,5 +1,10 @@
##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##contig=<ID=chr1,length=248956422>
#CHROM POS ID REF ALT QUAL FILTER INFO
chr1 1 . A C . . .
##contig=<ID=chr2,length=248956422>
##INFO=<ID=END,Number=1,Type=Integer,Description="Stop position of the interval">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT b
chr1 1 . A C . . . GT 0/1
chr2 3 . G A,<*> . . . GT 0/1
chr2 4 . T <*> . . END=8 GT 0/0
2 changes: 2 additions & 0 deletions test/merge.gvcf.10.fa
Original file line number Diff line number Diff line change
@@ -1,2 +1,4 @@
>chr1
ACGT
>chr2
ACGTACGT
2 changes: 1 addition & 1 deletion test/merge.gvcf.2.out
Original file line number Diff line number Diff line change
Expand Up @@ -19,8 +19,8 @@
2 21444416 . G <*> . . END=21444427;MinDP=5;QS=1,0 PL:DP 0,15,125:5 .:. .:.
2 21444428 . C <*> . . END=21444429;MinDP=2;QS=2,0 PL:DP 0,15,125:5 0,6,51:2 .:.
2 21444430 . TCAA T,TAA,<*> 0 . MinDP=2;QS=1.60366,0.304878,0.0914634,0 PL:DP:DV 37,0,79,35,73,113,113,113,113,113:5:2 0,51,51,51,51,51,6,51,51,51:2:. .:.:.
2 21444431 . C <*> . . MinDP=4;QS=1,0 PL:DP 0,12,110:4 .:. .:.
2 21444431 . CA C,CAAACAAAAAA 0 . QS=0.75,0.25,1 PL:DP:DV 0,4,10,10,10,10:4:1 28,28,28,3,28,0:1:1 .:.:.
2 21444431 . C <*> . . MinDP=4;QS=1,0 PL:DP 0,12,110:4 .:. .:.
2 21444433 . C <*> 0 . END=21444444;QS=0.75,0.25 PL:DP:DV 0,4,10:4:1 .:.:. .:.:.
3 1 . C <*> 0 . END=4;MinDP=33;QS=0.75,0.25 PL:DP:DV 0,4,10:4:1 .:.:. .:.:.
3 5 . C <*>,T 0 . MinDP=33;QS=1.5,0.25,0.25 PL:DP:DV 0,4,10,10,10,10:4:1 0,10,10,4,10,10:4:1 .:.:.
Expand Down
10 changes: 10 additions & 0 deletions test/merge.gvcf.5.1.out
Original file line number Diff line number Diff line change
@@ -0,0 +1,10 @@
##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##ALT=<ID=NON_REF,Description="Represents any possible alternative allele at this location">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##INFO=<ID=END,Number=1,Type=Integer,Description="Stop position of the interval">
##contig=<ID=chr22,length=10514891>
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT 11 13
chr22 10510106 . T <NON_REF> . . END=10510112 GT 0/0 ./.
chr22 10510113 . C <NON_REF>,T . . . GT 0/0 2/2
chr22 10510114 . N <NON_REF> . . END=10510117 GT 0/0 ./.
7 changes: 7 additions & 0 deletions test/merge.gvcf.5.a.vcf
Original file line number Diff line number Diff line change
@@ -0,0 +1,7 @@
##fileformat=VCFv4.2
##ALT=<ID=NON_REF,Description="Represents any possible alternative allele at this location">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##INFO=<ID=END,Number=1,Type=Integer,Description="Stop position of the interval">
##contig=<ID=chr22,length=10514891>
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT 11
chr22 10510106 . T <NON_REF> . . END=10510117 GT 0/0
7 changes: 7 additions & 0 deletions test/merge.gvcf.5.b.vcf
Original file line number Diff line number Diff line change
@@ -0,0 +1,7 @@
##fileformat=VCFv4.2
##ALT=<ID=NON_REF,Description="Represents any possible alternative allele at this location">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##INFO=<ID=END,Number=1,Type=Integer,Description="Stop position of the interval">
##contig=<ID=chr22,length=10514891>
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT 13
chr22 10510113 . C T,<NON_REF> . . . GT 1/1
6 changes: 3 additions & 3 deletions test/test.pl
Original file line number Diff line number Diff line change
Expand Up @@ -91,8 +91,6 @@
run_test(\&test_vcf_merge,$opts,in=>['merge.gvcf.10.a','merge.gvcf.10.b'],out=>'merge.gvcf.10.2.out',args=>'-m none');
run_test(\&test_vcf_merge,$opts,in=>['merge.gvcf.10.a','merge.gvcf.10.b'],out=>'merge.gvcf.10.3.out',args=>'-g {PATH}/merge.gvcf.10.fa');
run_test(\&test_vcf_merge,$opts,in=>['merge.gvcf.10.a','merge.gvcf.10.b'],out=>'merge.gvcf.10.4.out',args=>'-g {PATH}/merge.gvcf.10.fa -m none');
run_test(\&test_vcf_merge,$opts,in=>['merge.gvcf.10.b','merge.gvcf.10.a'],out=>'merge.gvcf.10.5.out',args=>'-g {PATH}/merge.gvcf.10.fa');
run_test(\&test_vcf_merge,$opts,in=>['merge.gvcf.10.b','merge.gvcf.10.a'],out=>'merge.gvcf.10.4.out',args=>'-g {PATH}/merge.gvcf.10.fa -m none');
run_test(\&test_vcf_merge,$opts,in=>['merge.noidx.a','merge.noidx.b','merge.noidx.c'],out=>'merge.noidx.abc.out',args=>'');
run_test(\&test_vcf_merge,$opts,in=>['merge.noidx.a','merge.noidx.b','merge.noidx.c'],out=>'merge.noidx.abc.out',args=>'--no-index',noidx=>1);
run_test(\&test_vcf_merge,$opts,in=>['merge.8.a','merge.8.b'],out=>'merge.8.out',args=>'');
Expand All @@ -104,6 +102,8 @@
run_test(\&test_vcf_merge,$opts,in=>['merge.10.a','merge.10.b'],out=>'merge.10.3.out',args=>'-m snp-ins-del');
run_test(\&test_vcf_merge,$opts,in=>['merge.mrules.1.a','merge.mrules.1.b'],out=>'merge.mrules.1.1.out',args=>'--gvcf -');
run_test(\&test_vcf_merge,$opts,in=>['merge.mrules.1.a','merge.mrules.1.b'],out=>'merge.mrules.1.2.out',args=>'--gvcf - -M AD:.,PL:.');
run_test(\&test_vcf_merge,$opts,in=>['merge.gvcf.5.a','merge.gvcf.5.b'],out=>'merge.gvcf.5.1.out',args=>'--gvcf -');
run_test(\&test_vcf_merge,$opts,in=>['merge.gvcf.5.a','merge.gvcf.5.b'],out=>'merge.gvcf.5.1.out',args=>'--gvcf - --merge none');
# run_test(\&test_vcf_merge_big,$opts,in=>'merge_big.1',out=>'merge_big.1.1',nsmpl=>79000,nfiles=>79,nalts=>486,args=>''); # commented out for speed
run_test(\&test_vcf_query,$opts,in=>'query.string',out=>'query.string.1.out',args=>q[-f '%CHROM\\t%POS\\t%CLNREVSTAT\\n' -i'CLNREVSTAT="criteria_provided,_conflicting_interpretations"']);
run_test(\&test_vcf_query,$opts,in=>'query.string',out=>'query.string.1.out',args=>q[-f '%CHROM\\t%POS\\t%CLNREVSTAT\\n' -i'CLNREVSTAT="criteria_provided" || CLNREVSTAT="_conflicting_interpretations"']);
Expand Down Expand Up @@ -1284,7 +1284,7 @@ sub test_vcf_merge
$args =~ s/{PATH}/$$opts{path}/g;
my $files = join(' ',@files);
test_cmd($opts,%args,cmd=>"$$opts{bin}/bcftools merge --no-version $args $files", exp_fix=>1);
test_cmd($opts,%args,cmd=>"$$opts{bin}/bcftools merge -Ob $args $files | $$opts{bin}/bcftools view | grep -v ^##bcftools_", exp_fix => 1);
test_cmd($opts,%args,cmd=>"$$opts{bin}/bcftools merge --no-version -Ob $args $files | $$opts{bin}/bcftools view --no-version | grep -v ^##bcftools_", exp_fix => 1);
}
}
sub test_vcf_isec
Expand Down
Loading

0 comments on commit 60c65eb

Please sign in to comment.