Skip to content

Add support for VCFv4.4 / VCFv4.5 Number= fields #1874

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 17 commits into
base: develop
Choose a base branch
from
Open

Add support for VCFv4.4 / VCFv4.5 Number= fields #1874

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
17 commits
Select commit Hold shift + click to select a range
408d66b
Add initial support for VCFv4.4 / VCFv4.5 Number= fields
daviesrob Jan 16, 2025
0ddccae
Add tests for VCF Number= header parsing
daviesrob Jul 8, 2025
f9a23c8
Beef up bcf_hdr_check_sanity with a more complete set of Number=types.
jkbonfield Apr 17, 2025
e914656
Add more nuance for header Number=? checking.
jkbonfield Apr 17, 2025
63feb78
Make bcf_remove_allele_set() work with local alleles
daviesrob Jun 25, 2025
8a2f9f8
Add bcf_remove_allele_set() tests
daviesrob Jul 1, 2025
6f1fafc
Fix up bcf_remove_allele_set() type=String handling
daviesrob Jul 9, 2025
0a39f1e
Add tests for string types
daviesrob Jul 9, 2025
c954b48
Make bcf_remove_allele_set() work for INFO tags with multiples of A
daviesrob Jul 25, 2025
1288e9f
Add tests for CICN, CIEND, CILEN, CIPOS, MEINFO, and METRANS tags
daviesrob Jul 25, 2025
5ceb231
Make bcf_remove_allele_set() work for <CNV:TR> related tags
daviesrob Jul 31, 2025
0841858
Add tests for <CNV:TR> related tags
daviesrob Jul 31, 2025
55c819c
Add MQ, MQ0 and NS checking; missed in earlier commit
jkbonfield Sep 11, 2025
ce53372
Fix check on valid GTs
daviesrob Sep 23, 2025
7d4d067
Use nR_ori when checking for valid LAA values
daviesrob Sep 23, 2025
8b6689c
Preserve LAA tag when all referenced alleles are removed
daviesrob Sep 25, 2025
082160d
Remove CNV:TR-related INFO tags when empty
daviesrob Sep 25, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Makefile
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -872,7 +872,7 @@ test/test_time_funcs.o: test/test_time_funcs.c config.h $(hts_time_funcs_h)
test/test_view.o: test/test_view.c config.h $(cram_h) $(htslib_sam_h) $(htslib_vcf_h) $(htslib_hts_log_h)
test/test_faidx.o: test/test_faidx.c config.h $(htslib_faidx_h)
test/test_index.o: test/test_index.c config.h $(htslib_sam_h) $(htslib_vcf_h)
test/test-vcf-api.o: test/test-vcf-api.c config.h $(htslib_hts_h) $(htslib_vcf_h) $(htslib_kstring_h) $(htslib_kseq_h)
test/test-vcf-api.o: test/test-vcf-api.c config.h $(htslib_hts_h) $(htslib_vcf_h) $(htslib_vcfutils_h) $(htslib_kbitset_h) $(htslib_kstring_h) $(htslib_kseq_h)
test/test-vcf-sweep.o: test/test-vcf-sweep.c config.h $(htslib_vcf_sweep_h)
test/test-bcf-sr.o: test/test-bcf-sr.c config.h $(htslib_hts_defs_h) $(htslib_synced_bcf_reader_h) $(htslib_hts_h) $(htslib_vcf_h)
test/test-bcf-translate.o: test/test-bcf-translate.c config.h $(htslib_vcf_h)
Expand Down
20 changes: 15 additions & 5 deletions htslib/vcf.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,11 +65,21 @@ extern "C" {
#define BCF_HT_STR 3
#define BCF_HT_LONG (BCF_HT_INT | 0x100) // BCF_HT_INT, but for int64_t values; VCF only!

#define BCF_VL_FIXED 0 // variable length
#define BCF_VL_VAR 1
#define BCF_VL_A 2
#define BCF_VL_G 3
#define BCF_VL_R 4
// Values for INFO / FORMAT "Number=" fields
#define BCF_VL_FIXED 0 ///< Integer defining a fixed number of items
#define BCF_VL_VAR 1 ///< Generic variable length ("Number=.")
#define BCF_VL_A 2 ///< One value per alternate allele
#define BCF_VL_G 3 ///< One value for each possible genotype
#define BCF_VL_R 4 ///< One value for each allele, including the reference

// The following was introduced in VCFv4.4 and is only valid for FORMAT header lines
#define BCF_VL_P 5 ///< One value for each allele value defined in GT

// The following were introduced in VCFv4.5 and are only valid for FORMAT header lines
#define BCF_VL_LA 6 ///< As BCF_VL_A, but only alt alleles listed in LAA are considered present
#define BCF_VL_LG 7 ///< As BCF_VL_G, but only alt alleles listed in LAA are considered present
#define BCF_VL_LR 8 ///< As BCF_VL_R, but only alt alleles listed in LAA are considered present
#define BCF_VL_M 9 ///< One value for each possible base modification for the corresponding ChEBI ID

/* === Dictionary ===

Expand Down
244 changes: 244 additions & 0 deletions test/test-vcf-api.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,8 @@ DEALINGS IN THE SOFTWARE. */

#include "../htslib/hts.h"
#include "../htslib/vcf.h"
#include "../htslib/vcfutils.h"
#include "../htslib/kbitset.h"
#include "../htslib/kstring.h"
#include "../htslib/kseq.h"

Expand Down Expand Up @@ -802,6 +804,246 @@ void test_rlen_values(void)
hts_set_log_level(logging);
}

static void test_vl_types(void)
{
const char *test_vcf = "data:,"
"##fileformat=VCFv4.5\n"
"##FILTER=<ID=PASS,Description=\"All filters passed\">\n"
"##contig=<ID=5,length=243199373>\n"
"##INFO=<ID=FIXED_1_INFO,Number=1,Type=Integer,Description=\"Fixed number 1\">\n"
"##INFO=<ID=FIXED_4_INFO,Number=4,Type=Float,Description=\"Fixed number 4\">\n"
"##INFO=<ID=VL_DOT_INFO,Number=.,Type=Integer,Description=\"Variable number\">\n"
"##INFO=<ID=VL_A_INFO,Number=A,Type=Integer,Description=\"One value for each ALT allele\">\n"
"##INFO=<ID=VL_G_INFO,Number=G,Type=Integer,Description=\"One value for each possible genotype\">\n"
"##INFO=<ID=VL_R_INFO,Number=R,Type=Integer,Description=\"One value for each allele including REF\">\n"
"##FORMAT=<ID=FIXED_1_FMT,Number=1,Type=String,Description=\"Fixed number 1\">\n"
"##FORMAT=<ID=FIXED_4_FMT,Number=4,Type=String,Description=\"Fixed number 4\">\n"
"##FORMAT=<ID=VL_DOT_FMT,Number=.,Type=String,Description=\"Variable number\">\n"
"##FORMAT=<ID=VL_A_FMT,Number=A,Type=Integer,Description=\"One value for each ALT allele\">\n"
"##FORMAT=<ID=VL_G_FMT,Number=G,Type=Integer,Description=\"One value for each possible genotype\">\n"
"##FORMAT=<ID=VL_R_FMT,Number=R,Type=Integer,Description=\"One value for each allele including REF\">\n"
"##FORMAT=<ID=VL_P_FMT,Number=P,Type=String,Description=\"One value for each allele value defined in GT\">\n"
"##FORMAT=<ID=VL_LA_FMT,Number=LA,Type=Integer,Description=\"One value for each local ALT allele\">\n"
"##FORMAT=<ID=VL_LG_FMT,Number=LG,Type=Integer,Description=\"One value for each local genotype\">\n"
"##FORMAT=<ID=VL_LR_FMT,Number=LR,Type=Integer,Description=\"One value for each local allele including REF\">\n"
"##FORMAT=<ID=VL_M_FMT,Number=M,Type=Integer,Description=\"One value for each posible base modification of the given type\">\n"
"#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tAAA\n";

typedef struct {
const char *id;
int type;
int expected_vl_code;
int expected_number;
} expected_types_t;

expected_types_t expected[] = {
{ "FIXED_1_INFO", BCF_HL_INFO, BCF_VL_FIXED, 1 },
{ "FIXED_4_INFO", BCF_HL_INFO, BCF_VL_FIXED, 4 },
{ "VL_DOT_INFO", BCF_HL_INFO, BCF_VL_VAR, 0xfffff },
{ "VL_A_INFO", BCF_HL_INFO, BCF_VL_A, 0xfffff },
{ "VL_G_INFO", BCF_HL_INFO, BCF_VL_G, 0xfffff },
{ "VL_R_INFO", BCF_HL_INFO, BCF_VL_R, 0xfffff },
{ "FIXED_1_FMT", BCF_HL_FMT, BCF_VL_FIXED, 1 },
{ "FIXED_4_FMT", BCF_HL_FMT, BCF_VL_FIXED, 4 },
{ "VL_DOT_FMT", BCF_HL_FMT, BCF_VL_VAR, 0xfffff },
{ "VL_A_FMT", BCF_HL_FMT, BCF_VL_A, 0xfffff },
{ "VL_G_FMT", BCF_HL_FMT, BCF_VL_G, 0xfffff },
{ "VL_R_FMT", BCF_HL_FMT, BCF_VL_R, 0xfffff },
{ "VL_P_FMT", BCF_HL_FMT, BCF_VL_P, 0xfffff },
{ "VL_LA_FMT", BCF_HL_FMT, BCF_VL_LA, 0xfffff },
{ "VL_LG_FMT", BCF_HL_FMT, BCF_VL_LG, 0xfffff },
{ "VL_LR_FMT", BCF_HL_FMT, BCF_VL_LR, 0xfffff },
{ "VL_M_FMT", BCF_HL_FMT, BCF_VL_M, 0xfffff },
};

htsFile *fp = hts_open(test_vcf, "r");
if (!fp) error("Failed to open test data\n");
bcf_hdr_t *hdr = bcf_hdr_read(fp);
if (!hdr) error("Failed to read BCF header\n");
size_t i;
for (i = 0; i < sizeof(expected)/sizeof(expected[0]); i++)
{
int id_num = bcf_hdr_id2int(hdr, BCF_DT_ID, expected[i].id);
if (id_num < 0)
{
error("Couldn't look up VCF header ID %s\n", expected[i].id);
}
int vl_code = bcf_hdr_id2length(hdr, expected[i].type, id_num);
if (vl_code != expected[i].expected_vl_code)
{
const char *length_types[] = {
"BCF_VL_FIXED", "BCF_VL_VAR", "BCF_VL_A", "BCF_VL_G",
"BCF_VL_R", "BCF_VL_P", "BCF_VL_LA", "BCF_VL_LG",
"BCF_VL_LR", "BCF_VL_M"
};
if (vl_code >= 0 && vl_code < sizeof(length_types)/sizeof(length_types[0]))
{
error("Unexpected length code for %s: expected %s got %s\n",
expected[i].id,
length_types[expected[i].expected_vl_code],
length_types[vl_code]);
}
else
{
error("Unexpected length code for %s: expected %s got %d\n",
expected[i].id,
length_types[expected[i].expected_vl_code], vl_code);
}
}
int num = bcf_hdr_id2number(hdr, expected[i].type, id_num);
if (num != expected[i].expected_number)
{
error("Unexpected number for %s: expected %d%s got %d%s\n",
expected[i].id,
expected[i].expected_number,
(expected[i].expected_number == 0xfffff
? " (= code for not fixed)" : ""),
num,
num == 0xfffff ? " (= code for not fixed)" : "");
}
}
bcf_hdr_destroy(hdr);
bcf_close(fp);
}

static int read_vcf_line(const char *line, bcf_hdr_t *hdr, bcf1_t *rec,
kstring_t *kstr)
{
int ret;
if (kputsn(line, strlen(line), ks_clear(kstr)) < 0)
return -1;

ret = vcf_parse(kstr, hdr, rec);
if (ret != 0) {
fprintf(stderr, "vcf_parse(\"%s\", hdr, rec) returned %d\n",
ks_c_str(kstr), ret);
}
return ret;
}

static void chomp(kstring_t *kstr)
{
if (kstr->l < 1)
return;
if (kstr->s[kstr->l - 1] == '\n')
kstr->s[--kstr->l] = '\0';
}

// Test bcf_remove_allele_set()
void test_bcf_remove_allele_set(void)
{
char header[] = "##fileformat=VCFv4.5\n"
"##FILTER=<ID=PASS,Description=\"All filters passed\">\n"
"##contig=<ID=5,length=243199373>\n"
"##INFO=<ID=AC,Number=A,Type=Integer,Description=\"Allele count\">\n"
"##INFO=<ID=AD,Number=R,Type=Integer,Description=\"Allele depth\">\n"
"##INFO=<ID=AF,Number=A,Type=Float,Description=\"Allele frequency\">\n"
"##INFO=<ID=CN,Number=A,Type=Float,Description=\"Copy number of CNV/breakpoint\">\n"
"##INFO=<ID=CICN,Number=.,Type=Float,Description=\"Confidence interval around copy number\">\n"
"##INFO=<ID=CIEND,Number=.,Type=Integer,Description=\"Confidence interval around the inferred END for symbolic structural variants\">\n"
"##INFO=<ID=CILEN,Number=.,Type=Integer,Description=\"Confidence interval for the SVLEN field\">\n"
"##INFO=<ID=CIPOS,Number=.,Type=Integer,Description=\"Confidence interval around POS for symbolic structural variants\">\n"
"##INFO=<ID=CIRB,Number=.,Type=Integer,Description=\"Confidence interval around RB\">\n"
"##INFO=<ID=CIRUC,Number=.,Type=Float,Description=\"Confidence interval around RUC\">\n"
"##INFO=<ID=IMPRECISE,Number=0,Type=Flag,Description=\"Imprecise structural variant\">\n"
"##INFO=<ID=MEINFO,Number=.,Type=String,Description=\"Mobile element info of the form NAME,START,END,POLARITY\">\n"
"##INFO=<ID=METRANS,Number=.,Type=String,Description=\"Mobile element transduction info of the form CHR,START,END,POLARITY\">\n"
"##INFO=<ID=RB,Number=.,Type=Integer,Description=\"Total number of bases in the corresponding repeat sequence\">\n"
"##INFO=<ID=RN,Number=A,Type=Integer,Description=\"Total number of repeat sequences in this allele\">\n"
"##INFO=<ID=RUB,Number=.,Type=Integer,Description=\"Number of bases in each individual repeat unit\">\n"
"##INFO=<ID=RUC,Number=.,Type=Float,Description=\"Repeat unit count of corresponding repeat sequence\">\n"
"##INFO=<ID=RUL,Number=.,Type=Integer,Description=\"Repeat unit length of the corresponding repeat sequence\">\n"
"##INFO=<ID=RUS,Number=.,Type=String,Description=\"Repeat unit sequence of the corresponding repeat sequence\">\n"
"##INFO=<ID=SVLEN,Number=A,Type=Integer,Description=\"Length of structural variant\">\n"
"##INFO=<ID=VL_A_STR_INFO,Number=A,Type=String,Description=\"INFO string Number=A\">\n"
"##INFO=<ID=VL_R_STR_INFO,Number=R,Type=String,Description=\"INFO string Number=R\">\n"
"##FORMAT=<ID=AD,Number=R,Type=Integer,Description=\"Allele depth\">\n"
"##FORMAT=<ID=EC,Number=A,Type=Integer,Description=\"Expected allele count\">\n"
"##FORMAT=<ID=GT,Number=1,Type=String,Description=\"Genotype\">\n"
"##FORMAT=<ID=PL,Number=G,Type=Integer,Description=\"List of Phred-scaled genotype likelihoods\">\n"
"##FORMAT=<ID=LAA,Number=.,Type=Integer,Description=\"Local alleles\">\n"
"##FORMAT=<ID=LAD,Number=LR,Type=Integer,Description=\"Local allele depth\">\n"
"##FORMAT=<ID=LEC,Number=LA,Type=Integer,Description=\"Local expected allele count\">\n"
"##FORMAT=<ID=LPL,Number=LG,Type=Integer,Description=\"List of Phred-scaled local genotype likelihoods\">\n"
"##FORMAT=<ID=VL_A_STR_FMT,Number=A,Type=String,Description=\"FMT string Number=A\">\n"
"##FORMAT=<ID=VL_G_STR_FMT,Number=G,Type=String,Description=\"FMT string Number=G\">\n"
"##FORMAT=<ID=VL_LA_STR_FMT,Number=LA,Type=String,Description=\"FMT string Number=LA\">\n"
"##FORMAT=<ID=VL_LG_STR_FMT,Number=LG,Type=String,Description=\"FMT string Number=LG\">\n"
"##FORMAT=<ID=VL_LR_STR_FMT,Number=LR,Type=String,Description=\"FMT string Number=LR\">\n"
"##FORMAT=<ID=VL_R_STR_FMT,Number=R,Type=String,Description=\"FMT string Number=R\">\n"
"##ALT=<ID=CNV:TR,Description=\"Tandem repeat determined based on DNA abundance\">\n"
"#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tAAA\tBBB\tCCC\n";
const char * inputs[] = {
"5\t110285\t.\tT\tC,<*>\t.\tPASS\tAC=1,0;AD=6,5,0;AF=0.99,0.01;VL_A_STR_INFO=alt_c,alt_nonref;VL_R_STR_INFO=ref,alt_c,alt_nonref\tGT:AD:EC:PL:VL_A_STR_FMT:VL_G_STR_FMT:VL_R_STR_FMT\t.:.:.:.:.:.:.\t0/1:6,5,0:4,0:114,0,15,35,73,113:alt_c,alt_nonref:gt_00,gt_01,gt_11,gt_02,gt_12,gt_22:ref,alt_c,alt_nonref\t.:.:.:.:.:.:.",
"5\t110290\t.\tT\tC,A\t.\tPASS\tAC=90,1;AD=6,5,6;AF=0.009,0.0001;VL_A_STR_INFO=alt_c,alt_a;VL_R_STR_INFO=ref,alt_c,alt_a\tGT:LAA:LAD:LEC:LPL:VL_LA_STR_FMT:VL_LG_STR_FMT:VL_LR_STR_FMT\t0/0:.:3:.:0:.:gt_00:ref\t0/1:1,2:3,2,0:44,27:114,0,15,35,73,113:alt_c,alt_a:gt_00,gt_01,gt_11,gt_02,gt_12,gt_22:ref,alt_c,alt_a\t1/1:1:0,3:46:110,15,0:alt_c:gt_00,gt_01,gt_11:ref,alt_c",
"5\t110350\t.\tT\t<INS>,<INS>\t.\tPASS\tIMPRECISE;SVLEN=100,200;CIEND=-50,50,-25,25;CIPOS=-10,10,-20,20\tGT\t0/1\t0/1\t0/1",
"5\t110500\t.\tT\t<CNV>,<CNV>\t.\tPASS\tIMPRECISE;SVLEN=50,100;CILEN=0,25,-25,25;CN=2,4;CICN=-0.5,1,-1.5,1.5\tGT\t0/1\t0/1\t0/1",
"5\t110700\t.\tA\t<INS:ME>,<INS:ME>\t.\tPASS\tMEINFO=AluY,1,260,+,FLAM_C,1,110,-;METRANS=1,94820,95080,+,1,129678,129788,-\tGT\t0/1\t0/1\t0/1",
"5\t112000\t.\tC\t<CNV:TR>,<CNV:TR>\t.\tPASS\tRN=2,1;RUS=CAG,TTG,CA;RUL=3,3,2;RB=12,6,6;RUC=4,2,3;RUB=3,3,3,3,3,3,2,2,2;SVLEN=18,6",
"5\t113000\t.\tT\tC,A\t.\tPASS\tAC=90,1;AD=6,5,6;AF=0.009,0.0001;VL_A_STR_INFO=alt_c,alt_a;VL_R_STR_INFO=ref,alt_c,alt_a\tGT:LAA:LAD:LEC:LPL:VL_LA_STR_FMT:VL_LG_STR_FMT:VL_LR_STR_FMT\t0/0:.:3:.:0:.:gt_00:ref\t0/1:1,2:3,2,0:44,27:114,0,15,35,73,113:alt_c,alt_a:gt_00,gt_01,gt_11,gt_02,gt_12,gt_22:ref,alt_c,alt_a\t1/1:1:0,3:46:110,15,0:alt_c:gt_00,gt_01,gt_11:ref,alt_c",
"5\t114000\t.\tT\tC,A\t.\tPASS\tAC=90,1;AD=6,5,6;AF=0.009,0.0001;VL_A_STR_INFO=alt_c,alt_a;VL_R_STR_INFO=ref,alt_c,alt_a\tGT:LAA:LAD:LEC:LPL:VL_LA_STR_FMT:VL_LG_STR_FMT:VL_LR_STR_FMT\t0/0:.:3:.:0:.:gt_00:ref\t0/1:1,2:3,2,0:44,27:114,0,15,35,73,113:alt_c,alt_a:gt_00,gt_01,gt_11,gt_02,gt_12,gt_22:ref,alt_c,alt_a\t1/1:1:0,3:46:110,15,0:alt_c:gt_00,gt_01,gt_11:ref,alt_c",
"5\t115000\t.\tC\t<CNV:TR>,<CNV:TR>\t.\tPASS\tRN=2,1;RUS=CAG,TTG,CA;RUL=3,3,2;RB=12,6,6;RUC=4,2,3;RUB=3,3,3,3,3,3,2,2,2;SVLEN=18,6",
};
const char * expected[] = {
// 2nd allele removed
"5\t110285\t.\tT\tC\t.\tPASS\tAC=1;AD=6,5;AF=0.99;VL_A_STR_INFO=alt_c;VL_R_STR_INFO=ref,alt_c\tGT:AD:EC:PL:VL_A_STR_FMT:VL_G_STR_FMT:VL_R_STR_FMT\t.:.:.:.:.:.:.\t0/1:6,5:4:114,0,15:alt_c:gt_00,gt_01,gt_11:ref,alt_c\t.:.:.:.:.:.:.",
"5\t110290\t.\tT\tC\t.\tPASS\tAC=90;AD=6,5;AF=0.009;VL_A_STR_INFO=alt_c;VL_R_STR_INFO=ref,alt_c\tGT:LAA:LAD:LEC:LPL:VL_LA_STR_FMT:VL_LG_STR_FMT:VL_LR_STR_FMT\t0/0:.:3:.:0:.:gt_00:ref\t0/1:1:3,2:44:114,0,15:alt_c:gt_00,gt_01,gt_11:ref,alt_c\t1/1:1:0,3:46:110,15,0:alt_c:gt_00,gt_01,gt_11:ref,alt_c",
"5\t110350\t.\tT\t<INS>\t.\tPASS\tIMPRECISE;SVLEN=100;CIEND=-50,50;CIPOS=-10,10\tGT\t0/1\t0/1\t0/1",
"5\t110500\t.\tT\t<CNV>\t.\tPASS\tIMPRECISE;SVLEN=50;CILEN=0,25;CN=2;CICN=-0.5,1\tGT\t0/1\t0/1\t0/1",
"5\t110700\t.\tA\t<INS:ME>\t.\tPASS\tMEINFO=AluY,1,260,+;METRANS=1,94820,95080,+\tGT\t0/1\t0/1\t0/1",
"5\t112000\t.\tC\t<CNV:TR>\t.\tPASS\tRN=2;RUS=CAG,TTG;RUL=3,3;RB=12,6;RUC=4,2;RUB=3,3,3,3,3,3;SVLEN=18",
// 1st allele removed
"5\t113000\t.\tT\tA\t.\tPASS\tAC=1;AD=6,6;AF=0.0001;VL_A_STR_INFO=alt_a;VL_R_STR_INFO=ref,alt_a\tGT:LAA:LAD:LEC:LPL:VL_LA_STR_FMT:VL_LG_STR_FMT:VL_LR_STR_FMT\t0/0:.:3:.:0:.:gt_00:ref\t0/.:1:3,0:27:114,35,113:alt_a:gt_00,gt_02,gt_22:ref,alt_a\t./.:.:0:.:110:.:gt_00:ref",
// Both alleles removed
"5\t114000\t.\tT\t.\t.\tPASS\tAD=6;VL_R_STR_INFO=ref\tGT:LAA:LAD:LEC:LPL:VL_LA_STR_FMT:VL_LG_STR_FMT:VL_LR_STR_FMT\t0/0:.:3:.:0:.:gt_00:ref\t0/.:.:3:.:114:.:gt_00:ref\t./.:.:0:.:110:.:gt_00:ref",
"5\t115000\t.\tC\t.\t.\tPASS\t.",
};

kstring_t kstr = KS_INITIALIZE;

bcf_hdr_t *hdr = bcf_hdr_init("r");
bcf1_t *rec = bcf_init();
struct kbitset_t *rm_set = kbs_init(3);
size_t i;

if (!hdr)
error("bcf_hdr_init() failed\n");

if (!rec)
error("bcf_init() failed\n");

if (!rm_set)
error("kbs_init() failed\n");

check0(ks_resize(&kstr, 1000));
check0(bcf_hdr_parse(hdr, header));
for (i = 0; i < sizeof(inputs)/sizeof(inputs[0]); i++)
{
check0(read_vcf_line(inputs[i], hdr, rec, &kstr));
kbs_clear(rm_set);
if (rec->pos == 113000 - 1) {
kbs_insert(rm_set, 1);
} else if (rec->pos >= 114000 - 1) {
kbs_insert(rm_set, 1);
kbs_insert(rm_set, 2);
} else {
kbs_insert(rm_set, 2);
}
check0(bcf_remove_allele_set(hdr, rec, rm_set));
check0(vcf_format(hdr, rec, ks_clear(&kstr)));
chomp(&kstr);
if (strcmp(expected[i], ks_c_str(&kstr)) != 0)
{
error("bcf_remove_allele_set() output differs\nExpected:\n%s\nGot:\n%s\n",
expected[i], ks_c_str(&kstr));
}
}
bcf_destroy(rec);
bcf_hdr_destroy(hdr);
ks_free(&kstr);
kbs_destroy(rm_set);
}

int main(int argc, char **argv)
{
char *fname = argc>1 ? argv[1] : "rmme.bcf";
Expand All @@ -815,9 +1057,11 @@ int main(int argc, char **argv)
iterator(fname);

// additional tests. quiet unless there's a failure.
test_vl_types();
test_get_info_values(fname);
test_invalid_end_tag();
test_open_format();
test_rlen_values();
test_bcf_remove_allele_set();
return 0;
}
Loading
Loading