Make long overlapping deletions work again after 8ee2a03

caot · Aug 21, 2018 · cc52852 · cc52852
1 parent 983f7da
commit cc52852
Show file tree

Hide file tree

Showing 6 changed files with 17 additions and 3 deletions.
diff --git a/csq.c b/csq.c
@@ -3159,23 +3159,25 @@ void tscript_init_ref(args_t *args, tscript_t *tr, const char *chr)
     int pad_end = len - (tr->end - tr->beg + 1 + pad_beg);
     if ( pad_beg + pad_end != 2*N_REF_PAD )
     {
-        char *ref = (char*) malloc(tr->end - tr->beg + 1 + 2*N_REF_PAD);
+        char *ref = (char*) malloc(tr->end - tr->beg + 1 + 2*N_REF_PAD + 1);
         for (i=0; i < N_REF_PAD - pad_beg; i++) ref[i] = 'N';
         memcpy(ref+i, tr->ref, len);
         len += i;
         for (i=0; i < N_REF_PAD - pad_end; i++) ref[i+len] = 'N';
+        ref[i+len] = 0;
         free(tr->ref);
         tr->ref = ref;
     }
 }
 
 static void sanity_check_ref(args_t *args, tscript_t *tr, bcf1_t *rec)
 {
-    int vbeg = rec->pos >= tr->beg ? 0 : tr->beg - rec->pos;
+    int vbeg = 0;
     int rbeg = rec->pos - tr->beg + N_REF_PAD;
+    if ( rbeg < 0 ) { vbeg += abs(rbeg); rbeg = 0; }
     char *ref = tr->ref + rbeg;
     char *vcf = rec->d.allele[0] + vbeg;
-    assert( vcf - rec->d.allele[0] < strlen(rec->d.allele[0]) );
+    assert( vcf - rec->d.allele[0] < strlen(rec->d.allele[0]) && ref - tr->ref < tr->end - tr->beg + 2*N_REF_PAD );
     int i = 0;
     while ( ref[i] && vcf[i] )
     {

diff --git a/test/csq/EDUMMY0002/EDUMMY0002.fa b/test/csq/EDUMMY0002/EDUMMY0002.fa
@@ -24,3 +24,6 @@ CCGTTTTATCTCTCGTTGCCGACGCTGGTGGTCAGTTATGCGCTGGGCATTGTGCTGAAA
 GGGCAACTGGCGCTGGGGATGCTGCTGGCGGTGCAAATAGCCGCGGGCGTGCTGGCGTTT
 GTGGTGATGATTGTGCTGTCACGCCATCCGCTGGTGGTGGAAGTGAAGCGTCAGTTTTGT
 CGCAGCGAAAAAATGAAAATGCTTTTACGGGCGGGGTGA
+>2
+ATGAGCTTACGTGAAAAAACCATGAGCGGCGCGAAGTGGTCGGCTATTGCCACGGTGATC
+
diff --git a/test/csq/EDUMMY0002/EDUMMY0002.fa.fai b/test/csq/EDUMMY0002/EDUMMY0002.fa.fai
@@ -1 +1,2 @@
 1	1479	3	60	61
+2	60	1510	60	61
diff --git a/test/csq/EDUMMY0002/EDUMMY0002.gff b/test/csq/EDUMMY0002/EDUMMY0002.gff
@@ -1,3 +1,6 @@
 1	.	gene	1	1479	.	+	.	ID=gene:wzx1_5;Name=wzx1_5;locus_tag=locus0000;gene_biotype=protein_coding;
 1	.	transcript	1	1479	.	+	.	ID=transcript:wzx1_5;Parent=gene:wzx1_5;locus_tag=locus0000;gene_biotype=protein_coding;
 1	.	CDS	1	1479	.	+	.	Parent=transcript:wzx1_5;gene=gene:wzx1_5
+2	.	gene	22	30	.	+	.	ID=gene:g2;Name=dummy2;gene_biotype=protein_coding;
+2	.	transcript	22	30	.	+	.	ID=transcript:t2;Parent=gene:g2;gene_biotype=protein_coding;
+2	.	CDS	22	30	.	+	.	Parent=transcript:t2
diff --git a/test/csq/EDUMMY0002/boundaries.txt b/test/csq/EDUMMY0002/boundaries.txt
@@ -7,3 +7,6 @@
 1479	A	C	stop_lost|wzx1_5|wzx1_5|protein_coding|+|493*>493C|1479A>C
 1479	A	C	stop_lost|wzx1_5|wzx1_5|protein_coding|+|493*>493C|1479A>C
 
+11	GTGAAAAAACCATGA	G	coding_sequence|dummy2|t2|protein_coding
+11	GTGAAAAAACCATGA	G	coding_sequence|dummy2|t2|protein_coding
+
diff --git a/test/csq/EDUMMY0002/boundaries.vcf b/test/csq/EDUMMY0002/boundaries.vcf
@@ -1,8 +1,10 @@
 ##fileformat=VCFv4.2
 ##reference=EDUMMY0002.fa
 ##contig=<ID=1,length=1479>
+##contig=<ID=2,length=1479>
 ##INFO=<ID=EXP,Number=1,Type=String,Description="Expected consequence">
 #CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
 1	6	.	C	A	.	.	EXP=missense|wzx1_5|wzx1_5|protein_coding|+|2S>2R|6C>A
 1	1474	.	G	C	.	.	EXP=missense|wzx1_5|wzx1_5|protein_coding|+|492G>492R|1474G>C
 1	1479	.	A	C	.	.	EXP=stop_lost|wzx1_5|wzx1_5|protein_coding|+|493*>493C|1479A>C
+2	11	.	GTGAAAAAACCATGA	G	.	.	EXP=coding_sequence|dummy2|t2|protein_coding