forked from samtools/bcftools
-
Notifications
You must be signed in to change notification settings - Fork 0
/
csq.c
3525 lines (3256 loc) · 142 KB
/
csq.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
/* The MIT License
Copyright (c) 2016-2023 Genome Research Ltd.
Author: Petr Danecek <[email protected]>
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
*/
/*
Things that would be nice to have
- dynamic N_REF_PAD
- for stop-lost events (also in frameshifts) report the number of truncated aa's
- memory could be greatly reduced by indexing gff (but it is quite compact already)
- deletions that go beyond transcript boundaries are not checked at sequence level
- alloc tscript->ref in hap_finalize, introduce fa_off_beg:16,fa_off_end:16
- see test/csq/ENST00000573314/insertion-overlap.vcf #1476288882
Read about transcript types here
http://vega.sanger.ac.uk/info/about/gene_and_transcript_types.html
http://www.ensembl.org/info/genome/variation/predicted_data.html
https://www.gencodegenes.org/pages/biotypes.html
List of supported biotypes
antisense
IG_C_gene
IG_D_gene
IG_J_gene
IG_LV_gene
IG_V_gene
lincRNA
lncRNA .. generic term for 3prime_overlapping_ncRNA, antisense, bidirectional_promoter_lncRNA, lincRNA, macro_lncRNA, non_coding, processed_transcript, sense_intronic, sense_overlapping
macro_lncRNA
miRNA
misc_RNA
Mt_rRNA
Mt_tRNA
polymorphic_pseudogene
processed_transcript
protein_coding, mRNA
ribozyme
rRNA
sRNA
scRNA
scaRNA
sense_intronic
sense_overlapping
snRNA
snoRNA
TR_C_gene
TR_D_gene
TR_J_gene
TR_V_gene
The gff parsing logic
We collect features such by combining gff lines A,B,C as follows:
A .. gene line with a supported biotype
A.ID=~/^gene:/
B .. transcript line referencing A with supported biotype
B.ID=~/^transcript:/ && B.Parent=~/^gene:A.ID/
C .. corresponding CDS, exon, and UTR lines:
C[3] in {"CDS","exon","three_prime_UTR","five_prime_UTR"} && C.Parent=~/^transcript:B.ID/
For coding biotypes ("protein_coding" or "polymorphic_pseudogene") the
complete chain link C -> B -> A is required. For the rest, link B -> A suffices.
The supported consequence types, sorted by impact:
splice_acceptor_variant .. end region of an intron changed (2bp at the 3' end of an intron)
splice_donor_variant .. start region of an intron changed (2bp at the 5' end of an intron)
stop_gained .. DNA sequence variant resulting in a stop codon
frameshift_variant .. number of inserted/deleted bases not a multiple of three, disrupted translational frame
stop_lost .. elongated transcript, stop codon changed
start_lost .. the first codon changed
inframe_altering .. combination of indels leading to unchanged reading frame and length
inframe_insertion .. inserted coding sequence, unchanged reading frame
inframe_deletion .. deleted coding sequence, unchanged reading frame
missense_variant .. amino acid (aa) change, unchanged length
splice_region_variant .. change within 1-3 bases of the exon or 3-8 bases of the intron
synonymous_variant .. DNA sequence variant resulting in no amino acid change
stop_retained_variant .. different stop codon
start_retained_variant .. start codon retained by indel realignment
non_coding_variant .. variant in non-coding sequence, such as RNA gene
5_prime_UTR_variant
3_prime_UTR_variant
intron_variant .. reported only if none of the above
intergenic_variant .. reported only if none of the above
The annotation algorithm.
The algorithm checks if the variant falls in a region of a supported type. The
search is performed in the following order, until a match is found:
1. idx_cds(gf_cds_t) - lookup CDS by position, create haplotypes, call consequences
2. idx_utr(gf_utr_t) - check UTR hits
3. idx_exon(gf_exon_t) - check for splice variants
4. idx_tscript(tscript_t) - check for intronic variants, RNAs, etc.
These regidx indexes are created by parsing a gff3 file as follows:
1. create the array "ftr" of all UTR, CDS, exons. This will be
processed later and pruned based on transcript types we want to keep.
In the same go, create the hash "id2tr" of transcripts to keep
(based on biotype) which maps from transcript_id to a transcript. At
the same time also build the hash "gid2gene" which maps from gene_id to
gf_gene_t pointer.
2. build "idx_cds", "idx_tscript", "idx_utr" and "idx_exon" indexes.
Use only features from "ftr" which are present in "id2tr".
3. clean data that won't be needed anymore: ftr, id2tr, gid2gene.
Data structures.
idx_cds, idx_utr, idx_exon, idx_tscript:
as described above, regidx structures for fast lookup of exons/transcripts
overlapping a region, the payload is a pointer to tscript.cds
*/
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <getopt.h>
#include <math.h>
#include <inttypes.h>
#include <htslib/hts.h>
#include <htslib/vcf.h>
#include <htslib/synced_bcf_reader.h>
#include <htslib/khash.h>
#include <htslib/khash_str2int.h>
#include <htslib/kseq.h>
#include <htslib/faidx.h>
#include <htslib/bgzf.h>
#include <errno.h>
#include <unistd.h>
#include <ctype.h>
#include "bcftools.h"
#include "filter.h"
#include "regidx.h"
#include "kheap.h"
#include "smpl_ilist.h"
#include "rbuf.h"
#include "gff.h"
#ifndef __FUNCTION__
# define __FUNCTION__ __func__
#endif
// Logic of the filters: include or exclude sites which match the filters?
#define FLT_INCLUDE 1
#define FLT_EXCLUDE 2
#define N_REF_PAD 10 // number of bases to avoid boundary effects
// How to treat phased/unphased genotypes
#define PHASE_REQUIRE 0 // --phase r
#define PHASE_MERGE 1 // --phase m
#define PHASE_AS_IS 2 // --phase a
#define PHASE_SKIP 3 // --phase s
#define PHASE_NON_REF 4 // --phase R
#define PHASE_DROP_GT 5 // --samples -
// Node types in the haplotype tree
#define HAP_CDS 0
#define HAP_ROOT 1
#define HAP_SSS 2 // start/stop/splice
#define CSQ_PRINTED_UPSTREAM (1<<0)
#define CSQ_SYNONYMOUS_VARIANT (1<<1)
#define CSQ_MISSENSE_VARIANT (1<<2)
#define CSQ_STOP_LOST (1<<3)
#define CSQ_STOP_GAINED (1<<4)
#define CSQ_INFRAME_DELETION (1<<5)
#define CSQ_INFRAME_INSERTION (1<<6)
#define CSQ_FRAMESHIFT_VARIANT (1<<7)
#define CSQ_SPLICE_ACCEPTOR (1<<8)
#define CSQ_SPLICE_DONOR (1<<9)
#define CSQ_START_LOST (1<<10)
#define CSQ_SPLICE_REGION (1<<11)
#define CSQ_STOP_RETAINED (1<<12)
#define CSQ_UTR5 (1<<13)
#define CSQ_UTR3 (1<<14)
#define CSQ_NON_CODING (1<<15)
#define CSQ_INTRON (1<<16)
//#define CSQ_INTERGENIC (1<<17)
#define CSQ_INFRAME_ALTERING (1<<18)
#define CSQ_UPSTREAM_STOP (1<<19) // adds * in front of the csq string
#define CSQ_INCOMPLETE_CDS (1<<20) // to remove START/STOP in incomplete CDS, see ENSG00000173376/synon.vcf
#define CSQ_CODING_SEQUENCE (1<<21) // cannot tell exactly what it is, but it does affect the coding sequence
#define CSQ_ELONGATION (1<<22) // symbolic insertion
#define CSQ_START_RETAINED (1<<23)
// Haplotype-aware consequences, printed in one vcf record only, the rest has a reference @12345
#define CSQ_COMPOUND (CSQ_SYNONYMOUS_VARIANT|CSQ_MISSENSE_VARIANT|CSQ_STOP_LOST|CSQ_STOP_GAINED| \
CSQ_INFRAME_DELETION|CSQ_INFRAME_INSERTION|CSQ_FRAMESHIFT_VARIANT| \
CSQ_START_LOST|CSQ_STOP_RETAINED|CSQ_INFRAME_ALTERING|CSQ_INCOMPLETE_CDS| \
CSQ_UPSTREAM_STOP|CSQ_START_RETAINED)
#define CSQ_START_STOP (CSQ_STOP_LOST|CSQ_STOP_GAINED|CSQ_STOP_RETAINED|CSQ_START_LOST|CSQ_START_RETAINED)
#define CSQ_PRN_STRAND(csq) ((csq)&CSQ_COMPOUND && !((csq)&(CSQ_SPLICE_ACCEPTOR|CSQ_SPLICE_DONOR|CSQ_SPLICE_REGION)))
#define CSQ_PRN_TSCRIPT (~(CSQ_INTRON|CSQ_NON_CODING))
#define CSQ_PRN_NMD (~(CSQ_INTRON|CSQ_NON_CODING))
#define CSQ_PRN_BIOTYPE CSQ_NON_CODING
// see kput_vcsq()
const char *csq_strings[] =
{
NULL,
"synonymous",
"missense",
"stop_lost",
"stop_gained",
"inframe_deletion",
"inframe_insertion",
"frameshift",
"splice_acceptor",
"splice_donor",
"start_lost",
"splice_region",
"stop_retained",
"5_prime_utr",
"3_prime_utr",
"non_coding",
"intron",
"intergenic",
"inframe_altering",
NULL,
NULL,
"coding_sequence",
"feature_elongation",
"start_retained"
};
/*
Structures related to VCF output:
vcsq_t
information required to assemble consequence lines such as "inframe_deletion|XYZ|ENST01|+|5TY>5I|121ACG>A+124TA>T"
vrec_t
single VCF record and csq tied to this record. (Haplotype can have multiple
consequences in several VCF records. Each record can have multiple consequences
from multiple haplotypes.)
csq_t
a top-level consequence tied to a haplotype
vbuf_t
pos2vbuf
VCF records with the same position clustered together for a fast lookup via pos2vbuf
*/
typedef struct _vbuf_t vbuf_t;
typedef struct _vcsq_t vcsq_t;
struct _vcsq_t
{
uint32_t strand:1,
type:31; // one of CSQ_* types
uint32_t trid;
uint32_t vcf_ial;
uint32_t biotype; // one of GF_* types
char *gene; // gene name
bcf1_t *ref; // if type&CSQ_PRINTED_UPSTREAM, ref consequence "@1234"
kstring_t vstr; // variant string, eg 5TY>5I|121ACG>A+124TA>T
};
typedef struct
{
bcf1_t *line;
uint32_t *fmt_bm; // bitmask of sample consequences with first/second haplotype interleaved
uint32_t nfmt:4, // the bitmask size (the number of integers per sample)
nvcsq:28, mvcsq;
vcsq_t *vcsq; // there can be multiple consequences for a single VCF record
}
vrec_t;
typedef struct
{
uint32_t pos;
vrec_t *vrec; // vcf line that this csq is tied to; needed when printing haplotypes (hap_stage_vcf)
int idx; // 0-based index of the csq at the VCF line, for FMT/BCSQ
vcsq_t type;
}
csq_t;
struct _vbuf_t
{
vrec_t **vrec; // buffer of VCF lines with the same position
int n, m;
uint32_t keep_until; // the maximum transcript end position
};
KHASH_MAP_INIT_INT(pos2vbuf, vbuf_t*)
/*
Structures related to haplotype-aware consequences in coding regions
hap_node_t
node of a haplotype tree. Each transcript has one tree
tscript_t
despite its general name, it is intended for coding transcripts only
hap_t
hstack_t
for traversal of the haplotype tree and braking combined
consequences into independent parts
*/
typedef struct _hap_node_t hap_node_t;
struct _hap_node_t
{
char *seq; // cds segment [parent_node,this_node)
char *var; // variant "ref>alt"
uint32_t type:2, // HAP_ROOT or HAP_CDS
csq:30; // this node's consequence
int dlen; // alt minus ref length: <0 del, >0 ins, 0 substitution
uint32_t rbeg; // variant's VCF position (0-based, inclusive)
int32_t rlen; // variant's rlen; alen=rlen+dlen; fake for non CDS types
uint32_t sbeg; // variant's position on the spliced reference transcript (0-based, inclusive, N_REF_PAD not included)
uint32_t icds; // which exon does this node's variant overlaps
hap_node_t **child, *prev; // children haplotypes and previous coding node
int nchild, mchild;
bcf1_t *cur_rec, *rec; // current VCF record and node's VCF record
int vcf_ial; // which VCF allele generated this node
uint32_t nend; // number of haplotypes ending in this node
int *cur_child, mcur_child; // mapping from the allele to the currently active child
csq_t *csq_list; // list of haplotype's consequences, broken by position (each corresponds to a VCF record)
int ncsq_list, mcsq_list;
};
#define TSCRIPT_AUX(x) ((tscript_t*)(x)->aux)
typedef struct
{
char *ref; // reference sequence, padded with N_REF_PAD bases on both ends
char *sref; // spliced reference sequence, padded with N_REF_PAD bases on both ends
hap_node_t *root; // root of the haplotype tree
hap_node_t **hap; // pointer to haplotype leaves, two for each sample
int nhap, nsref; // number of haplotypes and length of sref, including 2*N_REF_PAD
}
tscript_t;
static inline int cmp_tscript(gf_tscript_t **a, gf_tscript_t **b)
{
return ( (*a)->end < (*b)->end ) ? 1 : 0;
}
KHEAP_INIT(trhp, gf_tscript_t*, cmp_tscript)
typedef khp_trhp_t tr_heap_t;
typedef struct
{
hap_node_t *node; // current node
int ichild; // current child in the active node
int dlen; // total dlen, from the root to the active node
size_t slen; // total sequence length, from the root to the active node
}
hstack_t;
typedef struct
{
int mstack;
hstack_t *stack;
gf_tscript_t *tr; // tr->ref: spliced transcript on ref strand
kstring_t sseq; // spliced haplotype sequence on ref strand
kstring_t tseq; // the variable part of translated haplotype transcript, coding strand
kstring_t tref; // the variable part of translated reference transcript, coding strand
uint32_t sbeg; // stack's sbeg, for cases first node's type is HAP_SSS
int upstream_stop;
}
hap_t;
typedef struct _args_t
{
// the main regidx lookups, from chr:beg-end to overlapping features and
// index iterator
gff_t *gff;
regidx_t *idx_cds, *idx_utr, *idx_exon, *idx_tscript;
regitr_t *itr;
// text tab-delimited output (out) or vcf/bcf output (out_fh)
FILE *out;
htsFile *out_fh;
char *index_fn;
int write_index;
char *dump_gff;
// vcf
bcf_srs_t *sr;
bcf_hdr_t *hdr;
int hdr_nsmpl; // actual number of samples in the vcf, for bcf_update_format_values()
// include or exclude sites which match the filters
filter_t *filter;
char *filter_str;
int filter_logic; // FLT_INCLUDE or FLT_EXCLUDE
// samples to process
int sample_is_file;
char *sample_list;
smpl_ilist_t *smpl;
char *outdir, **argv, *fa_fname, *gff_fname, *output_fname;
char *bcsq_tag;
int argc, output_type, clevel;
int phase, verbosity, local_csq, record_cmd_line;
int ncsq2_max, nfmt_bcsq; // maximum number of csq per site that can be accessed from FORMAT/BCSQ (*2 and 1 bit skipped to avoid BCF missing values)
int ncsq2_small_warned;
int brief_predictions;
int unify_chr_names;
char *chr_name;
int mchr_name;
struct {
int unknown_chr,unknown_tscript_biotype,unknown_strand,unknown_phase,duplicate_id;
int unknown_cds_phase,incomplete_cds,wrong_phase,overlapping_cds;
} warned;
int rid; // current chromosome
tr_heap_t *active_tr; // heap of active transcripts for quick flushing
hap_t *hap; // transcript haplotype recursion
vbuf_t **vcf_buf; // buffered VCF lines to annotate with CSQ and flush
rbuf_t vcf_rbuf; // round buffer indexes to vcf_buf
kh_pos2vbuf_t *pos2vbuf; // fast lookup of buffered lines by position
gf_tscript_t **rm_tr; // buffer of transcripts to clean
int nrm_tr, mrm_tr;
csq_t *csq_buf; // pool of csq not managed by hap_node_t, i.e. non-CDS csqs
int ncsq_buf, mcsq_buf;
int force; // force run under various conditions. Currently only to skip out-of-phase transcripts
int n_threads; // extra compression/decompression threads
faidx_t *fai;
kstring_t str, str2;
int32_t *gt_arr, mgt_arr;
}
args_t;
// AAA, AAC, ...
const char *gencode = "KNKNTTTTRSRSIIMIQHQHPPPPRRRRLLLLEDEDAAAAGGGGVVVV*Y*YSSSS*CWCLFLF";
const uint8_t nt4[] =
{
4,4,4,4, 4,4,4,4, 4,4,4,4, 4,4,4,4,
4,4,4,4, 4,4,4,4, 4,4,4,4, 4,4,4,4,
4,4,4,4, 4,4,4,4, 4,4,4,4, 4,4,4,4,
4,4,4,4, 4,4,4,4, 4,4,4,4, 4,4,4,4,
4,0,4,1, 4,4,4,2, 4,4,4,4, 4,4,4,4,
4,4,4,4, 3,4,4,4, 4,4,4,4, 4,4,4,4,
4,0,4,1, 4,4,4,2, 4,4,4,4, 4,4,4,4,
4,4,4,4, 3
};
const uint8_t cnt4[] =
{
4,4,4,4, 4,4,4,4, 4,4,4,4, 4,4,4,4,
4,4,4,4, 4,4,4,4, 4,4,4,4, 4,4,4,4,
4,4,4,4, 4,4,4,4, 4,4,4,4, 4,4,4,4,
4,4,4,4, 4,4,4,4, 4,4,4,4, 4,4,4,4,
4,3,4,2, 4,4,4,1, 4,4,4,4, 4,4,4,4,
4,4,4,4, 0,4,4,4, 4,4,4,4, 4,4,4,4,
4,3,4,2, 4,4,4,1, 4,4,4,4, 4,4,4,4,
4,4,4,4, 0
};
#define dna2aa(x) gencode[ nt4[(uint8_t)(x)[0]]<<4 | nt4[(uint8_t)(x)[1]]<<2 | nt4[(uint8_t)(x)[2]] ]
#define cdna2aa(x) gencode[ cnt4[(uint8_t)(x)[2]]<<4 | cnt4[(uint8_t)(x)[1]]<<2 | cnt4[(uint8_t)(x)[0]] ]
static inline int ncsq2_to_nfmt(int ncsq2)
{
return 1 + (ncsq2 - 1) / 30;
}
static inline void icsq2_to_bit(int icsq2, int *ival, int *ibit)
{
*ival = icsq2 / 30;
*ibit = icsq2 % 30;
}
void init_data(args_t *args)
{
args->nfmt_bcsq = ncsq2_to_nfmt(args->ncsq2_max);
args->fai = fai_load(args->fa_fname);
if ( !args->fai ) error("Failed to load the fai index: %s\n", args->fa_fname);
args->gff = gff_init(args->gff_fname);
gff_set(args->gff,verbosity,args->verbosity);
gff_set(args->gff,strip_chr_names,args->unify_chr_names);
gff_set(args->gff,force_out_of_phase,args->force);
gff_set(args->gff,dump_fname,args->dump_gff);
gff_parse(args->gff);
args->idx_cds = gff_get(args->gff,idx_cds);
args->idx_utr = gff_get(args->gff,idx_utr);
args->idx_exon = gff_get(args->gff,idx_exon);
args->idx_tscript = gff_get(args->gff,idx_tscript);
args->itr = regitr_init(NULL);
args->rid = -1;
if ( args->filter_str )
args->filter = filter_init(args->hdr, args->filter_str);
args->pos2vbuf = kh_init(pos2vbuf);
args->active_tr = khp_init(trhp);
args->hap = (hap_t*) calloc(1,sizeof(hap_t));
// init samples
if ( !bcf_hdr_nsamples(args->hdr) ) args->phase = PHASE_DROP_GT;
if ( args->sample_list && !strcmp("-",args->sample_list) )
{
// ignore all samples
if ( args->output_type==FT_TAB_TEXT )
{
// significant speedup for plain VCFs
if (bcf_hdr_set_samples(args->hdr,NULL,0) < 0)
error_errno("[%s] Couldn't build sample filter", __func__);
}
args->phase = PHASE_DROP_GT;
}
else
args->smpl = smpl_ilist_init(args->hdr, args->sample_list, args->sample_is_file, SMPL_STRICT);
args->hdr_nsmpl = args->phase==PHASE_DROP_GT ? 0 : bcf_hdr_nsamples(args->hdr);
if ( args->output_type==FT_TAB_TEXT )
{
args->out = args->output_fname ? fopen(args->output_fname,"w") : stdout;
if ( !args->out ) error("Failed to write to %s: %s\n", !strcmp("-",args->output_fname)?"standard output":args->output_fname,strerror(errno));
fprintf(args->out,"# This file was produced by: bcftools +csq(%s+htslib-%s)\n", bcftools_version(),hts_version());
fprintf(args->out,"# The command line was:\tbcftools +%s", args->argv[0]);
int i;
for (i=1; i<args->argc; i++)
fprintf(args->out," %s",args->argv[i]);
fprintf(args->out,"\n");
fprintf(args->out,"# LOG\t[2]Message\n");
fprintf(args->out,"# CSQ"); i = 1;
fprintf(args->out,"\t[%d]Sample", ++i);
fprintf(args->out,"\t[%d]Haplotype", ++i);
fprintf(args->out,"\t[%d]Chromosome", ++i);
fprintf(args->out,"\t[%d]Position", ++i);
fprintf(args->out,"\t[%d]Consequence", ++i);
fprintf(args->out,"\n");
}
else
{
char wmode[8];
set_wmode(wmode,args->output_type,args->output_fname,args->clevel);
args->out_fh = hts_open(args->output_fname ? args->output_fname : "-", wmode);
if ( args->out_fh == NULL ) error("[%s] Error: cannot write to %s: %s\n", __func__,args->output_fname? args->output_fname : "standard output", strerror(errno));
if ( args->n_threads > 0)
hts_set_opt(args->out_fh, HTS_OPT_THREAD_POOL, args->sr->p);
if ( args->record_cmd_line ) bcf_hdr_append_version(args->hdr,args->argc,args->argv,"bcftools/csq");
bcf_hdr_printf(args->hdr,"##INFO=<ID=%s,Number=.,Type=String,Description=\"%s consequence annotation from BCFtools/csq, see http://samtools.github.io/bcftools/howtos/csq-calling.html for details. Format: Consequence|gene|transcript|biotype|strand|amino_acid_change|dna_change\">",args->bcsq_tag, args->local_csq ? "Local" : "Haplotype-aware");
if ( args->hdr_nsmpl )
bcf_hdr_printf(args->hdr,"##FORMAT=<ID=%s,Number=.,Type=Integer,Description=\"Bitmask of indexes to INFO/BCSQ, with interleaved first/second haplotype. Use \\\"bcftools query -f'[%%CHROM\\t%%POS\\t%%SAMPLE\\t%%TBCSQ\\n]'\\\" to translate.\">",args->bcsq_tag);
if ( bcf_hdr_write(args->out_fh, args->hdr)!=0 ) error("[%s] Error: cannot write the header to %s\n", __func__,args->output_fname?args->output_fname:"standard output");
if ( args->write_index && init_index(args->out_fh,args->hdr,args->output_fname,&args->index_fn)<0 ) error("Error: failed to initialise index for %s\n",args->output_fname);
}
if ( args->verbosity > 0 ) fprintf(stderr,"Calling...\n");
}
void destroy_data(args_t *args)
{
if ( args->ncsq2_small_warned )
fprintf(stderr,
"Note: Some samples had too many consequences to be represented in %d bytes. If you need to record them all,\n"
" the limit can be increased by running with `--ncsq %d`.\n",ncsq2_to_nfmt(args->ncsq2_max)/8,1+args->ncsq2_small_warned/2);
regitr_destroy(args->itr);
gff_destroy(args->gff);
if ( args->filter )
filter_destroy(args->filter);
khp_destroy(trhp,args->active_tr);
kh_destroy(pos2vbuf,args->pos2vbuf);
if ( args->smpl ) smpl_ilist_destroy(args->smpl);
int i,j,ret;
if ( args->out_fh )
{
if ( args->write_index )
{
if ( bcf_idx_save(args->out_fh)<0 )
{
if ( hts_close(args->out_fh)!=0 ) error("Error: close failed .. %s\n", args->output_fname?args->output_fname:"stdout");
error("Error: cannot write to index %s\n", args->index_fn);
}
free(args->index_fn);
}
ret = hts_close(args->out_fh);
}
else
ret = fclose(args->out);
if ( ret ) error("Error: close failed .. %s\n", args->output_fname?args->output_fname:"stdout");
for (i=0; i<args->vcf_rbuf.m; i++)
{
vbuf_t *vbuf = args->vcf_buf[i];
if ( !vbuf ) continue;
for (j=0; j<vbuf->m; j++)
{
if ( !vbuf->vrec[j] ) continue;
if ( vbuf->vrec[j]->line ) bcf_destroy(vbuf->vrec[j]->line);
free(vbuf->vrec[j]->fmt_bm);
free(vbuf->vrec[j]->vcsq);
free(vbuf->vrec[j]);
}
free(vbuf->vrec);
free(vbuf);
}
free(args->vcf_buf);
free(args->rm_tr);
free(args->csq_buf);
free(args->hap->stack);
free(args->hap->sseq.s);
free(args->hap->tseq.s);
free(args->hap->tref.s);
free(args->hap);
fai_destroy(args->fai);
free(args->gt_arr);
free(args->str.s);
free(args->str2.s);
free(args->chr_name);
}
/*
The splice_* functions are for consquences around splice sites: start,stop,splice_*
*/
#define SPLICE_VAR_REF 0 // ref: ACGT>ACGT, csq not applicable, skip completely
#define SPLICE_OUTSIDE 1 // splice acceptor or similar; csq set and is done, does not overlap the region
#define SPLICE_INSIDE 2 // overlaps coding region; csq can be set but coding prediction is needed
#define SPLICE_OVERLAP 3 // indel overlaps region boundary, csq set but could not determine csq
typedef struct
{
gf_tscript_t *tr;
struct {
int32_t pos, rlen, alen, ial;
char *ref, *alt;
bcf1_t *rec;
} vcf;
uint16_t check_acceptor:1, // check distance from exon start (fwd) or end (rev)
check_start:1, // this is the first coding exon (relative to transcript orientation), check first (fwd) or last (rev) codon
check_stop:1, // this is the last coding exon (relative to transcript orientation), check last (fwd) or first (rev) codon
check_donor:1, // as with check_acceptor
check_region_beg:1, // do/don't check for splices at this end, eg. in the first or last exon
check_region_end:1, //
check_utr:1, // check splice sites (acceptor/donor/region_*) only if not in utr
set_refalt:1; // set kref,kalt, if set, check also for synonymous events
uint32_t csq;
int tbeg, tend; // number of trimmed bases from beg and end of ref,alt allele
uint32_t ref_beg, // ref coordinates with spurious bases removed, ACC>AC can become AC>A or CC>C, whichever gives
ref_end; // a more conservative csq (the first and last base in kref.s)
kstring_t kref, kalt; // trimmed alleles, set only with SPLICE_OLAP
}
splice_t;
void splice_init(splice_t *splice, bcf1_t *rec)
{
memset(splice,0,sizeof(*splice));
splice->vcf.rec = rec;
splice->vcf.pos = rec->pos;
splice->vcf.rlen = rec->rlen;
splice->vcf.ref = rec->d.allele[0];
splice->csq = 0;
}
static inline void splice_build_hap(splice_t *splice, uint32_t beg, int len)
{
// len>0 .. beg is the first base, del filled from right
// len<0 .. beg is the last base, del filled from left
int rlen, alen, rbeg, abeg; // first base to include (ref coordinates)
if ( len<0 )
{
rlen = alen = -len;
rbeg = beg - rlen + 1;
int dlen = splice->vcf.alen - splice->vcf.rlen;
if ( dlen<0 && beg < splice->ref_end ) // incomplete del, beg is in the middle
dlen += splice->ref_end - beg;
abeg = rbeg + dlen;
}
else
{
rbeg = abeg = beg;
rlen = alen = len;
// check for incomplete del as above??
}
#define XDBG 0
#if XDBG
fprintf(stderr,"build_hap: rbeg=%d + %d abeg=%d \n",rbeg,rlen,abeg);
#endif
splice->kref.l = 0;
splice->kalt.l = 0;
// add the part before vcf.ref, in the vcf.ref and after vcf.ref
int roff; // how many vcf.ref bases already used
if ( rbeg < splice->vcf.pos )
{
assert( splice->tr->beg <= rbeg ); // this can be extended thanks to N_REF_PAD
kputsn(TSCRIPT_AUX(splice->tr)->ref + N_REF_PAD + rbeg - splice->tr->beg, splice->vcf.pos - rbeg, &splice->kref);
roff = 0;
}
else
roff = rbeg - splice->vcf.pos;
#if XDBG
fprintf(stderr,"r1: %s roff=%d\n",splice->kref.s,roff);
#endif
if ( roff < splice->vcf.rlen && splice->kref.l < rlen )
{
int len = splice->vcf.rlen - roff; // len still available in vcf.ref
if ( len > rlen - splice->kref.l ) len = rlen - splice->kref.l; // how much of ref allele is still needed
kputsn(splice->vcf.ref + roff, len, &splice->kref);
}
#if XDBG
fprintf(stderr,"r2: %s\n",splice->kref.s);
#endif
uint32_t end = splice->vcf.pos + splice->vcf.rlen; // position just after the ref allele
if ( splice->kref.l < rlen )
{
if ( end + rlen - splice->kref.l - 1 > splice->tr->end ) // trim, the requested sequence is too long (could be extended, see N_REF_PAD)
rlen -= end + rlen - splice->kref.l - 1 - splice->tr->end;
if ( splice->kref.l < rlen )
kputsn(TSCRIPT_AUX(splice->tr)->ref + N_REF_PAD + end - splice->tr->beg, rlen - splice->kref.l, &splice->kref);
}
#if XDBG
fprintf(stderr,"r3: %s\n",splice->kref.s);
#endif
int aoff;
if ( abeg < splice->vcf.pos )
{
assert( splice->tr->beg <= abeg );
kputsn(TSCRIPT_AUX(splice->tr)->ref + N_REF_PAD + abeg - splice->tr->beg, splice->vcf.pos - abeg, &splice->kalt);
aoff = 0;
}
else
aoff = abeg - splice->vcf.pos;
#if XDBG
fprintf(stderr,"a1: %s aoff=%d\n",splice->kalt.s,aoff);
#endif
if ( aoff < splice->vcf.alen && splice->kalt.l < alen )
{
int len = splice->vcf.alen - aoff; // len still available in vcf.alt
if ( len > alen - splice->kalt.l ) len = alen - splice->kalt.l; // how much of alt allele is still needed
kputsn(splice->vcf.alt + aoff, len, &splice->kalt);
aoff -= len;
}
if ( aoff < 0 ) aoff = 0;
else aoff--;
#if XDBG
fprintf(stderr,"a2: %s aoff=%d\n",splice->kalt.s,aoff);
#endif
end = splice->vcf.pos + splice->vcf.rlen; // position just after the ref allele
if ( splice->kalt.l < alen )
{
if ( end + alen + aoff - splice->kalt.l - 1 > splice->tr->end ) // trim, the requested sequence is too long
alen -= end + alen + aoff - splice->kalt.l - 1 - splice->tr->end;
if ( alen > 0 && alen > splice->kalt.l )
kputsn(TSCRIPT_AUX(splice->tr)->ref + aoff + N_REF_PAD + end - splice->tr->beg, alen - splice->kalt.l, &splice->kalt);
}
#if XDBG
fprintf(stderr,"a3: %s\n",splice->kalt.s);
fprintf(stderr," [%s]\n [%s]\n\n",splice->kref.s,splice->kalt.s);
#endif
}
void csq_stage(args_t *args, csq_t *csq, bcf1_t *rec);
static inline int csq_stage_utr(args_t *args, regitr_t *itr, bcf1_t *rec, uint32_t trid, uint32_t type, int ial)
{
while ( regitr_overlap(itr) )
{
gf_utr_t *utr = regitr_payload(itr, gf_utr_t*);
gf_tscript_t *tr = utr->tr;
if ( tr->id != trid ) continue;
csq_t csq;
memset(&csq, 0, sizeof(csq_t));
csq.pos = rec->pos;
csq.type.type = (utr->which==prime5 ? CSQ_UTR5 : CSQ_UTR3) | type;
csq.type.biotype = tr->type;
csq.type.strand = tr->strand;
csq.type.trid = tr->id;
csq.type.vcf_ial = ial;
csq.type.gene = tr->gene->name;
csq_stage(args, &csq, rec);
return csq.type.type;
}
return 0;
}
static inline void csq_stage_splice(args_t *args, bcf1_t *rec, gf_tscript_t *tr, uint32_t type, int ial)
{
#if XDBG
fprintf(stderr,"csq_stage_splice %d: type=%d\n",rec->pos+1,type);
#endif
if ( !type ) return;
csq_t csq;
memset(&csq, 0, sizeof(csq_t));
csq.pos = rec->pos;
csq.type.type = type;
csq.type.biotype = tr->type;
csq.type.strand = tr->strand;
csq.type.trid = tr->id;
csq.type.vcf_ial = ial;
csq.type.gene = tr->gene->name;
csq_stage(args, &csq, rec);
}
static inline const char *drop_chr_prefix(args_t *args, const char *chr)
{
if ( !args->unify_chr_names ) return chr;
if ( !strncasecmp("chr",chr,3) ) return chr+3;
return chr;
}
static inline const char *add_chr_prefix(args_t *args, const char *chr)
{
if ( !args->unify_chr_names ) return chr;
int len = strlen(chr);
hts_expand(char,len+4,args->mchr_name,args->chr_name);
memcpy(args->chr_name,"chr",3);
memcpy(args->chr_name+3,chr,len+1);
return args->chr_name;
}
static inline int splice_csq_ins(args_t *args, splice_t *splice, uint32_t ex_beg, uint32_t ex_end)
{
// coordinates that matter for consequences, eg AC>ACG trimmed to C>CG, 1bp
// before and after the inserted bases
if ( splice->tbeg || splice->vcf.ref[0]!=splice->vcf.alt[0] )
{
splice->ref_beg = splice->vcf.pos + splice->tbeg - 1;
splice->ref_end = splice->vcf.pos + splice->vcf.rlen - splice->tend;
}
else
{
if ( splice->tend ) splice->tend--;
splice->ref_beg = splice->vcf.pos;
splice->ref_end = splice->vcf.pos + splice->vcf.rlen - splice->tend;
}
#if XDBG
fprintf(stderr,"ins: %s>%s .. ex=%d,%d beg,end=%d,%d tbeg,tend=%d,%d check_utr=%d start,stop,beg,end=%d,%d,%d,%d\n", splice->vcf.ref,splice->vcf.alt,ex_beg,ex_end,splice->ref_beg,splice->ref_end,splice->tbeg,splice->tend,splice->check_utr,splice->check_start,splice->check_stop,splice->check_region_beg,splice->check_region_end);
#endif
int ret;
if ( splice->ref_beg >= ex_end ) // fully outside, beyond the exon
{
if ( splice->check_utr )
{
regitr_t *itr = regitr_init(NULL);
const char *chr = drop_chr_prefix(args, bcf_seqname(args->hdr,splice->vcf.rec));
if ( regidx_overlap(args->idx_utr,chr,splice->ref_beg+1,splice->ref_beg+1, itr) ) // adjacent utr
{
ret = csq_stage_utr(args, itr, splice->vcf.rec, splice->tr->id, splice->csq, splice->vcf.ial);
if ( ret!=0 )
{
regitr_destroy(itr);
return SPLICE_OUTSIDE; // overlaps utr
}
}
regitr_destroy(itr);
}
if ( !splice->check_region_end ) return SPLICE_OUTSIDE;
char *ref = NULL, *alt = NULL;
if ( splice->set_refalt ) // seq identity is checked only when tr->ref is available
{
splice_build_hap(splice, ex_end+1, N_SPLICE_REGION_INTRON);
ref = splice->kref.s, alt = splice->kalt.s;
}
if ( splice->ref_beg < ex_end + N_SPLICE_REGION_INTRON && splice->ref_end > ex_end + N_SPLICE_DONOR )
{
splice->csq |= CSQ_SPLICE_REGION;
if ( ref && !strncmp(ref,alt,N_SPLICE_REGION_INTRON) ) splice->csq |= CSQ_SYNONYMOUS_VARIANT;
}
if ( splice->ref_beg < ex_end + N_SPLICE_DONOR )
{
if ( splice->check_donor && splice->tr->strand==STRAND_FWD ) splice->csq |= CSQ_SPLICE_DONOR;
if ( splice->check_acceptor && splice->tr->strand==STRAND_REV ) splice->csq |= CSQ_SPLICE_ACCEPTOR;
if ( ref && !strncmp(ref,alt,N_SPLICE_DONOR) ) splice->csq |= CSQ_SYNONYMOUS_VARIANT;
}
csq_stage_splice(args, splice->vcf.rec, splice->tr, splice->csq, splice->vcf.ial);
return SPLICE_OUTSIDE;
}
if ( splice->ref_end < ex_beg || (splice->ref_end == ex_beg && !splice->check_region_beg) ) // fully outside, before the exon
{
if ( splice->check_utr )
{
regitr_t *itr = regitr_init(NULL);
const char *chr = drop_chr_prefix(args, bcf_seqname(args->hdr,splice->vcf.rec));
if ( regidx_overlap(args->idx_utr,chr,splice->ref_end-1,splice->ref_end-1, itr) ) // adjacent utr
{
ret = csq_stage_utr(args, itr, splice->vcf.rec, splice->tr->id, splice->csq, splice->vcf.ial);
if ( ret!=0 )
{
regitr_destroy(itr);
return SPLICE_OUTSIDE; // overlaps utr
}
}
regitr_destroy(itr);
}
if ( !splice->check_region_beg ) return SPLICE_OUTSIDE;
char *ref = NULL, *alt = NULL;
if ( splice->set_refalt ) // seq identity is checked only when tr->ref is available
{
splice_build_hap(splice, ex_beg - N_SPLICE_REGION_INTRON, N_SPLICE_REGION_INTRON);
ref = splice->kref.s, alt = splice->kalt.s;
}
if ( splice->ref_end > ex_beg - N_SPLICE_REGION_INTRON && splice->ref_beg < ex_beg - N_SPLICE_DONOR )
{
splice->csq |= CSQ_SPLICE_REGION;
if ( ref && !strncmp(ref,alt,N_SPLICE_REGION_INTRON) ) splice->csq |= CSQ_SYNONYMOUS_VARIANT;
}
if ( splice->ref_end > ex_beg - N_SPLICE_DONOR )
{
if ( splice->check_donor && splice->tr->strand==STRAND_REV ) splice->csq |= CSQ_SPLICE_DONOR;
if ( splice->check_acceptor && splice->tr->strand==STRAND_FWD ) splice->csq |= CSQ_SPLICE_ACCEPTOR;
if ( ref && !strncmp(ref+N_SPLICE_REGION_INTRON-N_SPLICE_DONOR,alt+N_SPLICE_REGION_INTRON-N_SPLICE_DONOR,N_SPLICE_DONOR) ) splice->csq |= CSQ_SYNONYMOUS_VARIANT;
}
csq_stage_splice(args, splice->vcf.rec, splice->tr, splice->csq, splice->vcf.ial);
return SPLICE_OUTSIDE;
}
// overlaps the exon or inside the exon
// possible todo: find better alignment for frameshifting variants?
if ( splice->ref_beg <= ex_beg + 2 ) // in the first 3bp
{
if ( splice->check_region_beg ) splice->csq |= CSQ_SPLICE_REGION;
if ( splice->tr->strand==STRAND_FWD ) { if ( splice->check_start ) splice->csq |= CSQ_START_LOST; }
else { if ( splice->check_stop ) splice->csq |= CSQ_STOP_LOST; }
}
if ( splice->ref_end > ex_end - 2 )
{
if ( splice->check_region_end ) splice->csq |= CSQ_SPLICE_REGION;
if ( splice->tr->strand==STRAND_REV ) { if ( splice->check_start ) splice->csq |= CSQ_START_LOST; }
else { if ( splice->check_stop ) splice->csq |= CSQ_STOP_LOST; }
}
if ( splice->set_refalt )
{
// Make sure the variant will not end up left aligned to avoid overlapping vcf records
// splice_build_hap(splice, splice->ref_beg, splice->vcf.alen - splice->tend - splice->tbeg + 1);
// splice->vcf.rlen -= splice->tbeg + splice->tend - 1;
// if ( splice->kref.l > splice->vcf.rlen ) { splice->kref.l = splice->vcf.rlen; splice->kref.s[splice->kref.l] = 0; }
if ( splice->ref_beg < splice->vcf.pos ) // this must have been caused by too much trimming from right
{
int dlen = splice->vcf.pos - splice->ref_beg;
assert( dlen==1 );
splice->tbeg += dlen;
if ( splice->tbeg + splice->tend == splice->vcf.rlen ) splice->tend -= dlen;
splice->ref_beg = splice->vcf.pos;
}
if ( splice->ref_end==ex_beg ) splice->tend--; // prevent zero-length ref allele
splice_build_hap(splice, splice->ref_beg, splice->vcf.alen - splice->tend - splice->tbeg + 1);
splice->vcf.rlen -= splice->tbeg + splice->tend - 1;
if ( splice->kref.l > splice->vcf.rlen ) { splice->kref.l = splice->vcf.rlen; splice->kref.s[splice->kref.l] = 0; }
}
csq_stage_splice(args, splice->vcf.rec, splice->tr, splice->csq, splice->vcf.ial);
return SPLICE_INSIDE;
}
int shifted_del_synonymous(args_t *args, splice_t *splice, uint32_t ex_beg, uint32_t ex_end)
{
static int small_ref_padding_warned = 0;
gf_tscript_t *tr = splice->tr;
// We know the VCF record overlaps the exon, but does it overlap the start codon?
if ( tr->strand==STRAND_REV && splice->vcf.pos + splice->vcf.rlen + 2 <= ex_end ) return 0;
if ( tr->strand==STRAND_FWD && splice->vcf.pos >= ex_beg + 3 ) return 0;
#if XDBG
fprintf(stderr,"shifted_del_synonymous: %d-%d %s\n",ex_beg,ex_end, tr->strand==STRAND_FWD?"fwd":"rev");
fprintf(stderr," %d .. %s > %s\n",splice->vcf.pos+1,splice->vcf.ref,splice->vcf.alt);
#endif
// is there enough ref sequence for the extension? All coordinates are 0-based
int ref_len = strlen(splice->vcf.ref);
int alt_len = strlen(splice->vcf.alt);
assert( ref_len > alt_len );
int ndel = ref_len - alt_len;
if ( tr->strand==STRAND_REV )
{
int32_t vcf_ref_end = splice->vcf.pos + ref_len - 1; // end pos of the VCF REF allele
int32_t tr_ref_end = splice->tr->end + N_REF_PAD; // the end pos of accessible cached ref seq
if ( vcf_ref_end + ndel > tr_ref_end )
{
if ( !small_ref_padding_warned )
{
fprintf(stderr,"Warning: Could not verify synonymous start/stop at %s:%d due to small N_REF_PAD. (Improve me?)\n",bcf_seqname(args->hdr,splice->vcf.rec),splice->vcf.pos+1);
small_ref_padding_warned = 1;
}
return 0;
}
char *ptr_vcf = splice->vcf.ref + alt_len; // the first deleted base in the VCF REF allele
char *ptr_ref = TSCRIPT_AUX(splice->tr)->ref + N_REF_PAD + (vcf_ref_end + 1 - splice->tr->beg); // the first ref base after the ndel bases deleted
#if XDBG
fprintf(stderr,"vcf: %s\nref: %s\n",ptr_vcf,ptr_ref);
#endif
int i = 0;
while ( ptr_vcf[i] && ptr_vcf[i]==ptr_ref[i] ) i++;
if ( ptr_vcf[i] ) return 0; // the deleted sequence cannot be replaced
}
else
{