VCFReader

1 import std.exception: assertThrown;
2 import std.stdio: writeln, writefln;
3 
4 hts_set_log_level(htsLogLevel.HTS_LOG_TRACE);
5 
6 
7 auto vw = VCFWriter("/dev/null");
8 
9 vw.addHeaderLineRaw("##INFO=<ID=NS,Number=1,Type=Integer,Description=\"Number of Samples With Data\">");
10 vw.addHeaderLineKV("INFO", "<ID=DP,Number=1,Type=Integer,Description=\"Total Depth\">");
11 // ##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency">
12 vw.addTag!"INFO"("AF", "A", "Integer", "Number of Samples With Data");
13 vw.addHeaderLineRaw("##contig=<ID=20,length=62435964,assembly=B36,md5=f126cdf8a6e0c7f379d618ff66beb2da,species=\"Homo sapiens\",taxonomy=x>"); // @suppress(dscanner.style.long_line)
14 vw.addHeaderLineRaw("##FILTER=<ID=q10,Description=\"Quality below 10\">");
15 
16 // Exercise header
17 assert(vw.vcfhdr.nsamples == 0);
18 vw.addSample("NA12878");
19 assert(vw.vcfhdr.nsamples == 1);
20 
21 auto r = new VCFRecord(vw.vcfhdr, bcf_init1());
22 
23 r.chrom = "20";
24 assert(r.chrom == "20");
25 assertThrown(r.chrom = "chr20");   // headerline chromosome is "20" not "chr20"
26 
27 r.pos = 999_999;
28 assert(r.pos == 999_999);
29 r.pos = 62_435_964 + 1;     // Exceeds contig length
30 
31 // Test ID field
32 // note that in an empty freshly initialized bcf1_t, the ID field is "" rather than "."
33 // Will consider adding initialization code
34 //writefln("ID: %s", r.id);
35 r.id = "";
36 assert(r.id == ".");
37 r.id = ".";
38 assert(r.id == ".");
39 r.id = "rs001";
40 assert(r.id == "rs001");
41 r.addID("rs999");
42 assert(r.id == "rs001;rs999");
43 r.addID("rs001");   // test duplicate checking
44 assert(r.id == "rs001;rs999");
45 
46 // Test REF/ALT allele setters and getters
47 // many overloads to test for good coverage
48 // Test matrix: Set {zero, one, two, three} alleles * {set by string, set by array} = 8 test cases
49 //  * also, there is setAlleles(ref, alt) and setAlleles(ref, alt1, ...) 
50 //  * TODO: will also need to retreive alt alleles as array and string
51 
52 string[] alleles_array;
53 
54 // Zero alleles
55 r.alleles("");
56 const auto zs = r.allelesAsArray;
57 const auto zsr = r.refAllele;
58 assert(zs.length == 0);
59 assert(zsr == "");
60 
61 r.alleles(alleles_array);
62 const auto za = r.allelesAsArray;
63 const auto zar = r.refAllele;
64 assert(za.length == 0);
65 assert(zar == "");
66 
67 // One allele
68 r.alleles("C");
69 const auto os = r.allelesAsArray;
70 const auto osr = r.refAllele;
71 assert(os.length == 1 && os[0] == "C");
72 assert(osr == "C");
73 
74 r.alleles(["C"]);
75 const auto oa = r.allelesAsArray;
76 const auto oar = r.refAllele;
77 assert(oa.length == 1 && oa[0] == "C");
78 assert(oar == "C");
79 
80 // Two alleles
81 r.alleles("C,T");
82 const auto ts = r.allelesAsArray;
83 const auto tsr= r.refAllele;
84 assert(ts == ["C", "T"]);
85 assert(tsr== "C");
86 
87 r.alleles(["C", "T"]);
88 const auto ta = r.allelesAsArray;
89 const auto tar= r.refAllele;
90 assert(ta == ["C", "T"]);
91 assert(tar== "C");
92 
93 const taaa = r.altAllelesAsArray;
94 const taas = r.altAllelesAsString;
95 assert(taaa == ["T"]);
96 assert(taas == "T");
97 
98 // alternate setter for >= 2 alleles
99 r.setAlleles("A", "G");
100 assert(r.allelesAsArray == ["A", "G"]);
101 
102 // Three alleles
103 r.alleles("A,C,T");
104 assert(r.allelesAsArray == ["A", "C", "T"]);
105 assert(r.refAllele == "A");
106 assert(r.altAllelesAsString == "C,T");
107 
108 // alternate the alleles for testing purposes
109 r.alleles(["G", "A", "C"]);
110 assert(r.allelesAsArray == ["G", "A", "C"]);
111 assert(r.refAllele == "G");
112 assert(r.altAllelesAsString == "A,C");
113 
114 r.setAlleles("A", "C", "T");
115 assert(r.allelesAsArray == ["A", "C", "T"]);
116 assert(r.refAllele == "A");
117 assert(r.altAllelesAsString == "C,T");
118 
119 
120 // Test QUAL
121 r.qual = 1.0;
122 assert(r.qual == 1.0);
123 // now test setting qual without unpacking
124 // TODO: see https://forum.dlang.org/post/hebouvswxlslqhovzaia@forum.dlang.org, once resolved (or once factory function written),
125 //  add template value param BCF_UN_STR
126 auto rr = new VCFRecord(vw.vcfhdr, "20\t17330\t.\tT\tA\t3\t.\tNS=3;DP=11;AF=0.017\n"); // @suppress(dscanner.style.long_line)
127 rr.qual = 3.0;
128 assert(rr.qual == 3.0);
129 
130 
131 // Test FILTER
132 assert(r.hasFilter("PASS"));
133 assert(r.hasFilter("."));
134 // add q10 filter (is in header)
135 r.filter = "q10";
136 assert(r.filter == "q10");
137 assert(!r.hasFilter("PASS"));   // if has another filter, no longer PASS (unless added explicitly)
138 assert(!r.hasFilter("."));      // if has another filter, no longer has . (which is PASS, or no filter [spec conflict?])
139 
140 // add q30 filteR (not in header)
141 r.filter = "q30";
142 assert(r.filter == "q10");  // i.e., unchanged
143 
144 assert(r.hasFilter("q10"));
145 assert(!r.hasFilter("q30"));
146 
147 r.removeAllFilters();
148 assert(r.hasFilter("."));
149 r.addFilter("q10");
150 assert(r.hasFilter("q10"));
151 r.removeFilter("q10");
152 assert(r.hasFilter("PASS"));
153 
154 
155 // Finally, print the records:
156 writefln("VCF records via toString:\n%s%s", r, rr);
157 
158 writeln("\ndhtslib.vcf: all tests passed\n");