Fix overflow for indels longer than 512bp. Fixes #1837

samtools · Jan 4, 2023 · c380ab0 · c380ab0
1 parent be18428
commit c380ab0
Show file tree

Hide file tree

Showing 4 changed files with 16 additions and 5 deletions.
diff --git a/NEWS b/NEWS
@@ -16,6 +16,11 @@ Changes affecting the whole of bcftools, or multiple commands:
 
 Changes affecting specific commands:
 
+* bcftools +allele-length
+
+    - fix overflow for indels longer than 512bp and aggregate alleles equal or larger than
+      that in the same bin (#1837)
+
 * bcftools annotate
 
     - Support sample reordering of annotation file (#1785)

diff --git a/plugins/allele-length.c b/plugins/allele-length.c
@@ -1,6 +1,7 @@
 /*  plugins/allele-length.c -- Calculate stats about the length of alleles
 
     Copyright (C) 2017-2018 GENOMICS plc.
+    Copyright (C) 2023 Genome Research Ltd.
 
     Author: Nicola Asuni <nicola.asuni@genomicsplc.com>
 
@@ -89,12 +90,16 @@ bcf1_t *process(bcf1_t *rec)
 {
     int rl = strlen(rec->d.allele[0]);
     int al = strlen(rec->d.allele[1]);
+    int ral = rl + al;
+    if ( rl >= MAXLEN ) rl = MAXLEN - 1;
+    if ( al >= MAXLEN ) al = MAXLEN - 1;
+    if ( ral >= MAXLEN ) ral = MAXLEN - 1;
     reflen[rl] += 1;
     altlen[al] += 1;
-    refaltlen[(rl + al)] += 1;
+    refaltlen[ral] += 1;
     if ((contain_non_base(rec->d.allele[0])) || (contain_non_base(rec->d.allele[1])))
     {
-        xrefaltlen[(rl + al)] += 1;
+        xrefaltlen[ral] += 1;
         numxvar++;
     }
     numvar++;

diff --git a/test/query.allele-length.tsv b/test/query.allele-length.tsv
@@ -1,6 +1,6 @@
 LENGTH	REF	ALT	REF+ALT	REF+ALT WITH NON-BASE NUCLEOTIDES
 0	0	0	0	0
-1	53	56	0	0
+1	53	57	0	0
 2	9	10	43	0
 3	1	0	19	1
 4	0	0	1	0
@@ -510,5 +510,5 @@ LENGTH	REF	ALT	REF+ALT	REF+ALT WITH NON-BASE NUCLEOTIDES
 508	0	0	0	0
 509	0	0	0	0
 510	0	0	0	0
-511	0	0	0	0
-			66	1
+511	1	0	1	0
+			67	1
diff --git a/test/query.nucleotide.vcf b/test/query.nucleotide.vcf
@@ -148,3 +148,4 @@
 1	11014	rs28484712	G	A	.	.	RS=28484712;RSPOS=11014;dbSNPBuildID=125;SSR=0;SAO=0;VP=0x050000020005000002000100;GENEINFO=DDX11L1:100287102;WGT=1;VC=SNV;R5;ASP
 1	11022	rs28775022	G	A	.	.	RS=28775022;RSPOS=11022;dbSNPBuildID=125;SSR=0;SAO=0;VP=0x050000020005000002000100;GENEINFO=DDX11L1:100287102;WGT=1;VC=SNV;R5;ASP
 1	11063	rs561109771	T	G	.	.	RS=561109771;RSPOS=11063;dbSNPBuildID=142;SSR=0;SAO=0;VP=0x050000020005040024000100;GENEINFO=DDX11L1:100287102;WGT=1;VC=SNV;R5;ASP;VLD;KGPhase3;CAF=0.997,0.002995;COMMON=1
+1	22222	.	AGTTCAGGGCCTGGAGAGAAGGCTGGGAGGCAGGAGCTGGGTCTAAAGAGGCCATTGTAACGATGGAGCTGTGCCTGTGGAGGCTGTTGTGAGGCAGTAGCCTCATCTGCGGAGGCTGCCGTGACATAGGGTATGGGCCTAAATAGGCCATTGTGAGTCATGAGCTTGGTCTGTAGAGGCTGACTGGAGAAAGTTCTGGGCCTGGAGAGGCTGCCGGGAGGTAGGAGCTGGGCCAAAAGATGTAAGCACATTTGCATTTATTAGGCACTTTATTTCCATTATTACACTGTAATATATAATAAAATAATTATAGAACTCACCATAATGTAGAATTAGTGGGCGTGTTAAGCTTGTTTTCCTGCAACTGAATGGTCCCACCTGAGCGTGATGGGAGAAAGTGACAGATCAATAGGTATTAGATTCTCATAAGGACAGCGCAACCTCGATCCCTCACATGCACGGTTCACAACAGGGTGCGTTCTCCTATGAGCATCTAATGCTGCTGCTCATCTGAGAAGGTGGAGCTCAGGTGGGAATGTGAGCAAAGGGGAGTGGCTCTAAATACAGACGAAGCTTCCCTCACTCCCTCACTCGACACCGCTCACCTCCTGCTGTGTGGCTCCTTGCGGCTCCATGGCTCAGGGGTTGGGGACCCCTGCTCAAGTGCATCCAAAGCGACCCTTCCCACACCAGTCTTCACAGTGGTCAAGGGCAGCAACCACTTAGCTCCCAAGGCATGTGCCTCAGCTGGCATTTCGTCACAATCAACAGTAAGTGGTAGCTTGAGTCACTGTGAGGTCACCTACTGGAAATCACCAGCACCCCATTTCCCACTGGCAAAGAGCTCAGCACTGCCCCCTGGGAAACCAAACCTATGCCCAAATCCCATCTGTGTGGGTCTACCTCCTGGGACCCTTCCTAACATATTAGTCAGAGTCCAATCAGGAAGCATAAACCACTCAAAAGTTTAAAGTGGTAAAATTTAATACAGAGAATTATTCATTATAACAGGTGAACAGCATAATGAGAGATTGGCTAGCACAGAGTAAAGAGAACTCTAGAGAATATGGGACTAGCCCAGGCCAGGCATGGTGGCTCATGCCTGAAATTCCAGCCATTACAGAAGCTAATGCAGGAGGATTGCTTAAGGCCAGGAGCTAGAGACCGGTCTGGACGACACAGTGAGACCCTGTCTCTATCCAAAAGAAGAAAAAAGTTAGCTGGGGGTGGTGGTGCACACTTGTAGTCCCAGCTACTCAGAATGCTGAAGTTTGAGCCTGGGAGGTCAAGGCTGCAGTGAGGCATGATTATGCCACTACAGTCCAGCCGGATGACAGAGCAAGACCCTGTCTCAAAGAACAAAACCACAACAACCATTTACAGACAGAAAAGAAATAGAGCTAATAAGCTGAGGAAAGATGTTGAAATGTGACAAGTAAAGTAATATGAGGTCTTTTGTCTATTTAAAATAATCAAACAAAAAATGGCTTACGAAATTATAATACCCTGTGCTGGCAAAGGTGCAGTGAAATGGGCACTTTCTTATACTATGAGGGGTGGTTAAATTGTGTATAAGCCTTCCAGGGTAAAGCCTGTCAATTTTTTAAAATAATGGAGACAGGGTCTCACCATACTGCCATACTGCCTCCTCCAACTCTTGGCCTCAAGCAATCCTCCTCTCTTAGCCTCCCAAAGTGCTAAGATTATAGCTGGGAGGCACCCAAAACCCTGTCAATTTACATCAAGGGTAAGGAGAATGTCCATTCACCATGACTCACAGTAATCTTACTTCTGGGGAGACAATTCAATCTAAGCAAAAGGTCATCTGTACACACACAGTAAAAATCTGGGAGTAACTGAAGACAGAGTTGGTAAGTGAAATAAGAAACAGTTATAAGAAATTAAACTATGGTATCAATAGGCACCTGGTAAAAGGTCAGTTGATGTTAGCTGCTACTTTTTTGTTGTTTTGAGACAGGGTCTCACTCTGTCACCCAGGCTGGAGTGCAGAGGCCTGATCATGACTCACTGCAGTCTCAGCCTCCCTGGGCTCAAGTGATCCTCCCACCTCAGCCTCCCAAGTAGCTGGGACTACAGGAACATGCCACCACACTAGGCTAATTCATGTATTTTTCTGTAGGGATGGTGACTCCCCCTTTGTTTCCAAGGCCTATCGCAAACTCTTGGCCTCGAGCCATCCTCCTGCCTCAGCCTCCCAAAGTGTTGTGATTACCAGTGTGAGCCACCACACCTGGCCAGCTGCTACTTTTATCAATATTATTCTTATTCCACTCAATTAAAAATTATTATTTTCAAGGCTATGCAACAGTATGTATCCTACAGCGTAATTGTAAAAACATACACAGTCGTCATCCCTCAGTATACAGAATTAGTTCCAGCCCCCCATCTCTGCATATACCAAAATCCATGCTTACTCACGTTTCGCTGTCACCCCTCTAGAATCCACGTATACGAAAATTCCAAATGTTAGTTGGGCATAGTGGCAAGCACCTGTAGTCTCAGCCACGTGGGAGGTTGAGGTGGGAGGATCGCTTCAGCCTGGAAGGTTGAGGCTGCAGTCAGCTGCGATAGCACTACTACACTCCAGCCTTGGACAACAGAGGGAGACCCTGTCTCAGAAAAAAAAACAAAATAAAACAGGTTAGAAATTGTAATGAGGTCTGCTGGGCAAAATTCCATATAAGCAATGTATAAATTAATAAAGCAAATCGTGATAAATTAGTACGATTGACTTTCTGGAGTTTCTGACAATAAAAGTAAGGAAAATGCAGAACACAAAGACAGAGAGTAAAAAGAGAAATTAGGAAAGCATTCTACATGTTGAATAGGAAGACACTGGCCATGTTCGTGCAGCGGCAGTATGTCGTGACATGACATACCTTGGAGAGAAGTTAACAGATGAGGAAGTTGATAAAAATCATCAGAGAAGCAAAATACTGGTAGCGACACTCAAGTAAACCATGAAATTTCCATAACTTATGTCAGCAAAGTGGGAATATTGTACAGTGTGTGTTGAAGTTCCTATACAACATTGTTTATCTGCCTTTTGTTTGTTTGTAAGGAATGTACATACTAAAAGTTCTTCTTGCTGTCAAAAGAATATGCGTGAATAAGTCATTTTAACTTATTCTTCTGTTTTTCTTTTATCTTCCTGCCATCATCCCACAGCCTTACTTTAGAAATTTCTTTTTTAGAAAATTGAACAAGTGCTCCCTGTGGTGGCACATACCTCGAGGATGGGAGGCAGGGGTGGAAGGGTCACTTGAGGCCATTAGTTTGACACCAGCCTGGCCAACAAAGTGAGACCCCGTGTCTACAAAACAATTTAAAAATTAGCCAAGTATCGTCATGTATACCTACAGTCCCAGCTATCTGAACTTACTGAGAAT	A	.	.	.