dhtslib.file.file source code

1 module dhtslib.file.file;
2 
3 import core.stdc.stdio : SEEK_SET;
4 
5 import std.stdio;
6 import std.string : toStringz, fromStringz;
7 import std.utf : toUTFz;
8 import std.traits : isSomeString;
9 import std.algorithm : map;
10 import std.array : array;
11 
12 import dhtslib.memory;
13 import dhtslib.file.iterator;
14 import dhtslib : initKstring;
15 import htslib;
16 
17 enum HtslibFileFormatMode
18 {
19     Binary = 'b',
20     Cram = 'c',
21     Gzip = 'g',
22     Uncompressed = 'u',
23     Bgzf = 'z',
24     ZlibCompress0 = '0',
25     ZlibCompress1 = '1',
26     ZlibCompress2 = '2',
27     ZlibCompress3 = '3',
28     ZlibCompress4 = '4',
29     ZlibCompress5 = '5',
30     ZlibCompress6 = '6',
31     ZlibCompress7 = '7',
32     ZlibCompress8 = '8',
33     ZlibCompress9 = '9',
34 }
35 
36 
37 /** 
38  * Some shortcut mode strings for different file types
39  * As a reminder:
40  *      b  binary format (BAM, BCF, etc) rather than text (SAM, VCF, etc)
41         c  CRAM format
42         g  gzip compressed
43         u  uncompressed
44         z  bgzf compressed
45         [0-9]  zlib compression level
46  */
47 enum HtslibFileWriteMode
48 {
49     Bam             = "wb",
50     Cram            = "wc",
51     Sam             = "w",
52     UncompressedBam = "wb0",
53     GzippedSam      = "wg",
54     BgzippedSam     = "wz",
55 
56     Bcf             = "wb",
57     UncompressedBcf = "wb0",
58     Vcf             = "w",
59     GzippedVcf      = "wg",
60     BgzippedVcf     = "wz",
61 
62     Text            = "w",
63     GzippedText     = "wg",
64     BgzippedText    = "wz",
65     
66 }
67 
68 /** 
69  * HtslibFile is an abstraction for htslib's htsFile using dhtslib.memory for reference counting.
70  * 
71  * Ideally this struct and its methods can replace file io across dhtslib through a 
72  * standard interface. It can also centralize the ability to duplicate a file so that 
73  * it can be iterated several times. It pairs with HtslibIterator.
74  */
75 struct HtslibFile
76 {
77     /// dhtslib.memory htsFile* rc wrapper
78     /// File pointer
79     HtsFile fp; 
80 
81     /// D File if used
82     File f;
83 
84     /// dhtslib.memory kstring_t rc wrapper
85     /// Since we need the filename as a null terminated string anyway
86     /// just use kstring
87     Kstring fn;
88 
89     char[] mode;
90 
91     /// SAM/BAM/CRAM index 
92     HtsIdx idx;
93 
94     /// SAM/BAM/CRAM header
95     BamHdr bamHdr;
96     /// BCF/VCF header
97     BcfHdr bcfHdr;
98     /// text header
99     Kstring textHdr;
100 
101     /// Tabix index 
102     Tbx tbx;
103 
104     /// is file end reached?
105     bool eof;
106 
107     /// offset of file after header
108     /// if it has been loaded 
109     off_t headerOffset = -1;
110 
111     /// allow HtslibFile to be used as 
112     /// underlying ptr type
113     alias getFilePtr this;
114 
115     /// get underlying file pointer wrapper
116     @property nothrow pure @nogc
117     ref inout(HtsFile) getFilePtr() inout return
118     {
119         return this.fp;
120     }
121 
122     /// File or string ctor
123     this(T)(T f)
124     if ((isSomeString!T || is(T == File)))
125     {
126         this(f, "r");
127     }
128 
129     /// File or string ctor with mode
130     this(T1, T2)(T1 f, T2 mode)
131     if ((isSomeString!T1 || is(T1 == File)) && (isSomeString!T2 || is(T2 == HtslibFileWriteMode)))
132     {
133         this.mode = mode.dup ~ '\0';
134         // open file
135         static if (isSomeString!T1)
136         {
137             this.fn = Kstring(initKstring());
138             ks_initialize(this.fn);
139             kputsn(f.ptr, f.length, this.fn);
140             this.fp = HtsFile(hts_open(this.fn.s, this.mode.ptr));
141         }
142         else static if (is(T1 == File))
143         {
144             this.fn = Kstring(initKstring());
145             ks_initialize(this.fn);
146             kputsn(f.name.ptr, f.name.length, this.fn);
147 
148             auto hf = hdopen(f.fileno, m.ptr);
149             this.fp = HtsFile(hts_hopen(hf, this.fn.s, m.ptr));
150         }
151         else static assert(0);
152     }
153 
154     /// Duplicate the file with the same offset
155     HtslibFile dup()
156     {
157         // open a new file
158         auto filename = fromStringz(this.fn.s).idup;
159         auto newFile = HtslibFile(filename, this.mode);
160 
161         // seek to the current offset
162         newFile.seek(this.tell);
163 
164         // copy internal fields
165         newFile.idx = this.idx;
166         newFile.bamHdr = this.bamHdr;
167         newFile.bcfHdr = this.bcfHdr;
168         newFile.textHdr = this.textHdr;
169         newFile.tbx = this.tbx;
170         newFile.eof = this.eof;
171         newFile.headerOffset = this.headerOffset;
172         return newFile;
173     }
174 
175     /// set extra multithreading
176     void setExtraThreads(int extra)
177     {
178         hts_set_threads(this.fp, extra);
179     }
180 
181     /// get file offset
182     off_t tell()
183     {
184         if(this.fp.is_bgzf) return bgzf_tell(this.fp.fp.bgzf);
185         else return htell(this.fp.fp.hfile);
186     }
187 
188     /// seek to offset in file
189     void seek(long loc)
190     {
191         long err;
192         if(this.fp.is_bgzf) err = bgzf_seek(this.fp.fp.bgzf, loc, SEEK_SET);
193         else err = hseek(this.fp.fp.hfile, loc, SEEK_SET);
194         if(err < 0) hts_log_error(__FUNCTION__, "Error seeking htsFile");
195     }
196 
197     /// reset file reading
198     void reset()
199     {
200         this.seek(0);
201     }
202 
203     /// reset file reading
204     void resetToFirstRecord()
205     {
206         if(headerOffset != -1)
207             this.seek(headerOffset);
208         else {
209             hts_log_error(__FUNCTION__, "Cannot resetToFirstRecord, header offset unknown");
210             hts_log_error(__FUNCTION__, "Has the header been loaded?");
211         }
212     }
213 
214     /// Read a SAM/BAM header, VCF/BCF header, or text header and internally store it
215     void loadHeader(char commentChar = '#')
216     {
217         switch(this.fp.format.format){
218             case htsExactFormat.sam:
219             case htsExactFormat.bam:
220             case htsExactFormat.cram:
221                 this.loadHeader!BamHdr(commentChar);
222                 break;
223             case htsExactFormat.vcf:
224             case htsExactFormat.bcf:
225                 this.loadHeader!BcfHdr(commentChar);    
226                 break;
227             default:
228                 this.loadHeader!Kstring(commentChar);
229         }
230     }
231 
232     /// Read a SAM/BAM header, VCF/BCF header, or text header and internally store it
233     void loadHeader(T)(char commentChar)
234     if(is(T == BamHdr) || is(T == BcfHdr) || is(T == Kstring))
235     {
236         static if(is(T == BamHdr)){
237             this.bamHdr = BamHdr(sam_hdr_read(this.fp));
238         }
239         else static if(is(T == BcfHdr)){
240             this.bcfHdr = BcfHdr(bcf_hdr_read(this.fp));
241         }
242         else static if(is(T == Kstring)){
243             // parse header from a text file
244             this.textHdr = Kstring(initKstring);
245             ks_initialize(this.textHdr);
246 
247             auto ks = Kstring(initKstring);
248             ks_initialize(ks);
249 
250             // peek at first char, if commentChar (#)
251             // save the line and repeat
252             if(this.fp.is_bgzf){
253                 
254                 while(true){
255                     if(cast(char) bgzf_peek(this.fp.fp.bgzf) == commentChar){
256                         hts_getline(this.fp, cast(int)'\n', ks);
257                         kputs(ks.s, this.textHdr);
258                         kputc(cast(int)'\n', this.textHdr);
259                     } else break;
260                 }
261             }else{
262                 while(true){
263                     char c;
264                     hpeek(this.fp.fp.hfile, &c, 1);
265                     if(c  == commentChar){
266                         hts_getline(this.fp, cast(int)'\n', ks);
267                         kputs(ks.s, this.textHdr);
268                         kputc(cast(int)'\n', this.textHdr);
269                     } else break;
270                 }
271                 this.textHdr.l--;
272                 kputc(cast(int)'\0', this.textHdr);
273             }
274         }
275         // set header offset
276         this.headerOffset = this.tell;
277     }
278 
279     /// set header for a HtlsibFile 
280     void setHeader(T)(T hdr)
281     if(is(T == BamHdr) || is(T == BcfHdr))
282     {
283         static if(is(T == BamHdr))
284             this.bamHdr = hdr;
285         else static if(is(T == BcfHdr))
286             this.bcfHdr = hdr;
287         else static assert(0);
288     }
289 
290     /// write internally stored header
291     void writeHeader()
292     {
293         assert(this.bamHdr.initialized || this.bcfHdr.initialized);
294         assert(!(this.bamHdr.initialized && this.bcfHdr.initialized));
295         assert((this.bamHdr.initialized && this.bamHdr.getRef != null) 
296             || (this.bcfHdr.initialized && this.bcfHdr.getRef != null));
297         int err;
298         if(this.bamHdr.initialized)
299             err = sam_hdr_write(this.fp, this.bamHdr);
300         else
301             err = bcf_hdr_write(this.fp, this.bcfHdr);
302         if(err < 0) hts_log_error(__FUNCTION__, "Error writing SAM/BAM/VCF/BCF header");
303     }
304 
305     /// load BAM/BCF index
306     HtsIdx loadHtsIndex()
307     {
308         if(this.fp.format.format == htsExactFormat.bam || this.fp.format.format == htsExactFormat.cram)
309             this.idx = HtsIdx(sam_index_load(this.fp, cast(const(char)*)this.fn.s));
310         else if(this.fp.format.format == htsExactFormat.bcf)
311             this.idx = HtsIdx(bcf_index_load(cast(const(char)*)this.fn.s));
312         return this.idx;
313     }
314 
315     /// load SAM/BAM index from filename
316     HtsIdx loadHtsIndex(string idxFile)
317     {
318         if(this.fp.format.format == htsExactFormat.bam || this.fp.format.format == htsExactFormat.cram)
319             this.idx = HtsIdx(sam_index_load2(this.fp, this.fn.s, toStringz(idxFile)));
320         else if(this.fp.format.format == htsExactFormat.bcf)
321             this.idx = HtsIdx(bcf_index_load2(this.fn.s, toStringz(idxFile)));
322         return this.idx;
323     }
324 
325     /// load tabix index
326     Tbx loadTabixIndex()
327     {
328         this.tbx = Tbx(tbx_index_load(this.fn.s));
329         return this.tbx;
330     }
331 
332     /// load tabix index from file
333     Tbx loadTabixIndex(T)(T idxFile)
334     if(isSomeString!T)
335     {
336         this.tbx = Tbx(tbx_index_load2(this.fn.s, toStringz(idxFile)));
337         return this.tbx;
338     }
339 
340     /// Query htsFile with tid, start, and end
341     /// returns an HtslibIterator that has type T
342     /// requires index be loaded first
343     /// can use Tabix index or BAI/CSI index
344     auto query(T)(int tid, long beg, long end)
345     if(is(T == Bam1) || is(T == Bcf1) || is(T == Kstring))
346     {
347         // if T == Kstring
348         // we assume using tabix index
349         static if(is(T == Kstring)){
350             if(this.tbx != null){
351                 auto itr = HtsItr(tbx_itr_queryi(this.tbx, tid, beg, end));
352                 return HtslibIterator!Kstring(this, itr);
353             }else{
354                 throw new Exception("No TABIX index is loaded");
355             }
356         }else{
357             // BAI/CSI Index
358             if(this.idx != null){
359                 static if(is(T == Bam1)){
360                     auto itr = HtsItr(sam_itr_queryi(this.idx, tid, beg, end));
361                     return HtslibIterator!Bam1(this, itr);
362                 }else static if(is(T == Bcf1)){
363                     auto itr = HtsItr(hts_itr_query(this.idx, tid, beg, end, &bcf_readrec));
364                     return HtslibIterator!Bcf1(this, itr);
365                 }
366             // Tabix Index
367             }else if(this.tbx != null){
368                 static if(is(T == Bam1)){
369                     auto itr = HtsItr(tbx_itr_queryi(this.tbx, tid, beg, end));
370                     return HtslibIterator!Bam1(this, itr);
371                 }else static if(is(T == Bcf1)){
372                     auto itr = HtsItr(tbx_itr_queryi(this.tbx, tid, beg, end));
373                     return HtslibIterator!Bcf1(this, itr);
374                 }
375             }else{
376                 throw new Exception("No TABIX/BAI/CSI index is loaded");
377             }
378         }
379         
380         
381     }
382 
383     /// Query htsFile with string region
384     /// returns an HtslibIterator that has type T
385     /// requires index and header be loaded first
386     /// can use Tabix index or BAI/CSI index
387     auto query(T)(string region)
388     if(is(T == Bam1) || is(T == Bcf1))
389     {
390         // BAI/CSI Index
391         if(this.idx != null){
392             static if(is(T == Bam1)){
393                 auto itr = HtsItr(sam_itr_querys(this.idx, this.bamHdr, toStringz(region)));
394                 return HtslibIterator!Bam1(this, itr);
395             }else static if(is(T == Bcf1)){
396                 auto itr = HtsItr(bcf_itr_querys(this.idx, this.bcfHdr, toStringz(region)));
397                 return HtslibIterator!Bcf1(this, itr);
398             }
399         // Tabix Index
400         }else if(this.tbx != null){
401             static if(is(T == Bam1)){
402                 auto itr = HtsItr(tbx_itr_querys(this.tbx, toStringz(region)));
403                 return HtslibIterator!Bam1(this, itr);
404             }else static if(is(T == Bcf1)){
405                 auto itr = HtsItr(tbx_itr_querys(this.tbx, toStringz(region)));
406                 return HtslibIterator!Bcf1(this, itr);
407             }else static if(is(T == Kstring)){
408                 auto itr = HtsItr(tbx_itr_querys(this.tbx, toStringz(region)));
409                 return HtslibIterator!Kstring(this, itr);
410             }
411         }else{
412             throw new Exception("No TABIX/BAI/CSI index is loaded");
413         }
414     }
415 
416     /// Query htsFile with array of regions
417     /// returns an HtslibIterator that has type Bam1
418     /// requires index and header be loaded first
419     auto query(string[] regions)
420     {
421         auto cQueries = regions.map!(toUTFz!(char *)).array;
422 
423         assert(this.idx.getRef != null);
424         auto itr = HtsItr(sam_itr_regarray(this.idx, this.bamHdr, cQueries.ptr, cast(int)regions.length));
425         return HtslibIterator!Bam1(this, itr);
426     }
427 
428     /// iterate over all records in a file
429     /// returns a HtslibIterator
430     auto byRecord(T)()
431     if(is(T == Bam1) || is(T == Bcf1) || is(T == Kstring))
432     {
433         return HtslibIterator!T(this);
434     }
435 
436     /// write SAM/BAM/VCF/BCF record, string, or ubyte data
437     /// requires the header be loaded if writing SAM/BAM/VCF/BCF record
438     void write(T)(T rec)
439     if(isSomeString!T || is(T == Bam1) || is(T == Bcf1) || is(T == Kstring) || is(T == ubyte[]))
440     {
441         long err;
442         static if(is(T == Bam1)){
443             err = sam_write1(this.fp, this.bamHdr, rec);
444         }else static if(is(T == Bcf1)){
445             err = bcf_write(this.fp, this.bcfHdr, rec);
446         }else static if(isSomeString!T || is(T == ubyte[])){
447             if(this.fp.is_bgzf) err = bgzf_write(this.fp.fp.bgzf, rec.ptr, rec.length);
448             else err = hwrite(this.fp.fp.hfile, rec.ptr, rec.length);
449         }else static if(is(T == Kstring)){
450             if(this.fp.is_bgzf) err = bgzf_write(this.fp.fp.bgzf, rec.s, rec.l);
451             else err = hwrite(this.fp.fp.hfile, rec.s, rec.l);
452         }else static assert(0);
453         if(err < 0) hts_log_error(__FUNCTION__, "Error writing SAM/BAM/VCF/BCF record, string, or ubyte data");
454     }
455 
456     /// write a string with a newline appended
457     void writeln(T)(T rec)
458     if(isSomeString!T)
459     {
460         rec = rec ~ '\n';
461         write(rec);
462     }
463 
464     /// read a BAM/SAM/BCF/VCF record
465     auto readRecord(T)()
466     if(is(T == Bam1) || is(T == Bcf1) || is(T == Kstring))
467     {
468         long err;
469         static if(is(T == Bam1)){
470             assert(this.bamHdr.initialized && this.bamHdr.getRef != null);
471             auto rec = bam_init1;
472             err = sam_read1(this.fp, this.bamHdr, rec);
473         }
474         else static if(is(T == Bcf1)){
475             assert(this.bcfHdr.initialized && this.bcfHdr.getRef != null);
476             auto rec = bcf_init;
477             err = bcf_read(this.fp, this.bcfHdr, rec);
478         }else static if(is(T == Kstring)){
479             auto rec = initKstring;
480             ks_initialize(rec);
481             err = hts_getline(this.fp, cast(int)'\n', rec);
482         }
483         if(err < -1) hts_log_error(__FUNCTION__, "Error reading Bam/Bcf record");
484         else if(err == -1) eof = true;
485         return T(rec);
486     }
487 
488 }
489 
490 debug(dhtslib_unittest) unittest
491 {
492     hts_log_info(__FUNCTION__, "Testing HtslibFile text reading and writing");
493     {
494         auto f = HtslibFile("/tmp/htslibfile.test.txt", HtslibFileWriteMode.Text);
495         f.writeln("hello");
496         f.writeln("test");
497 
498         f = HtslibFile("/tmp/htslibfile.test.txt.gz", HtslibFileWriteMode.GzippedText);
499         f.writeln("hello");
500         f.writeln("test");
501 
502         f = HtslibFile("/tmp/htslibfile.test.txt.bgz", HtslibFileWriteMode.BgzippedText);
503         f.writeln("hello");
504         f.writeln("test");
505 
506         f = HtslibFile("/tmp/htslibfile.test.txt.9gz", "w9");
507         f.writeln("hello");
508         f.writeln("test");
509     }
510     {
511         auto f = HtslibFile("/tmp/htslibfile.test.txt");
512         assert(fromStringz(f.readRecord!Kstring.s) == "hello");
513         assert(fromStringz(f.readRecord!Kstring.s) == "test");
514 
515         f = HtslibFile("/tmp/htslibfile.test.txt.gz");
516         assert(fromStringz(f.readRecord!Kstring.s) == "hello");
517         assert(fromStringz(f.readRecord!Kstring.s) == "test");
518 
519         f = HtslibFile("/tmp/htslibfile.test.txt.bgz");
520         assert(fromStringz(f.readRecord!Kstring.s) == "hello");
521         assert(fromStringz(f.readRecord!Kstring.s) == "test");
522 
523         f = HtslibFile("/tmp/htslibfile.test.txt.9gz");
524         assert(fromStringz(f.readRecord!Kstring.s) == "hello");
525         assert(fromStringz(f.readRecord!Kstring.s) == "test");
526 
527     }
528     {
529         auto f = HtslibFile("/tmp/htslibfile.test.txt.gz");
530         assert(fromStringz(f.readRecord!Kstring.s) == "hello");
531         assert(fromStringz(f.readRecord!Kstring.s) == "test");
532     }
533     {
534         auto f = HtslibFile("/tmp/htslibfile.test.txt");
535         assert(f.byRecord!Kstring.map!(x=> fromStringz(x.s).idup).array == ["hello", "test"] );
536     }
537 }
538 
539 debug(dhtslib_unittest) unittest
540 {
541     hts_log_info(__FUNCTION__, "Testing HtslibFile SAM/BAM reading and writing");
542     import std.path:buildPath,dirName;
543     auto fn = buildPath(dirName(dirName(dirName(dirName(__FILE__)))),"htslib","test","range.bam");
544     {
545         auto f = HtslibFile(fn);
546         f.loadHeader;
547         auto h = f.bamHdr;
548         auto read = f.readRecord!Bam1();
549 
550         f = HtslibFile("/tmp/htslibfile.test.sam", HtslibFileWriteMode.Sam);
551         f.setHeader(h);
552         f.writeHeader;
553         f.write(read);
554         
555         f = HtslibFile("/tmp/htslibfile.test.bam", HtslibFileWriteMode.Bam);
556         f.setHeader(h);
557         f.writeHeader;
558         f.write(read);
559 
560         f = HtslibFile("/tmp/htslibfile.test.ubam", HtslibFileWriteMode.UncompressedBam);
561         f.setHeader(h);
562         f.writeHeader;
563         f.write(read);
564 
565         f = HtslibFile("/tmp/htslibfile.test.sam.gz", HtslibFileWriteMode.GzippedSam);
566         f.setHeader(h);
567         f.writeHeader;
568         f.write(read);
569 
570         f = HtslibFile("/tmp/htslibfile.test.sam.bgz", HtslibFileWriteMode.BgzippedSam);
571         f.setHeader(h);
572         f.writeHeader;
573         f.write(read);
574 
575     }
576     {
577 
578         auto f = HtslibFile("/tmp/htslibfile.test.sam");
579         f.loadHeader;
580         auto h = f.bamHdr;
581         auto read = f.readRecord!Bam1();
582         assert(fromStringz(bam_get_qname(read)) == "HS18_09653:4:1315:19857:61712");
583         
584         f = HtslibFile("/tmp/htslibfile.test.bam");
585         f.loadHeader;
586         h = f.bamHdr;
587         read = f.readRecord!Bam1();
588         assert(fromStringz(bam_get_qname(read)) == "HS18_09653:4:1315:19857:61712");
589 
590         f = HtslibFile("/tmp/htslibfile.test.ubam");
591         f.loadHeader;
592         h = f.bamHdr;
593         read = f.readRecord!Bam1();
594         assert(fromStringz(bam_get_qname(read)) == "HS18_09653:4:1315:19857:61712");
595 
596         f = HtslibFile("/tmp/htslibfile.test.sam.gz");
597         f.loadHeader;
598         h = f.bamHdr;
599         read = f.readRecord!Bam1();
600         assert(fromStringz(bam_get_qname(read)) == "HS18_09653:4:1315:19857:61712");
601 
602         f = HtslibFile("/tmp/htslibfile.test.sam.bgz");
603         f.loadHeader;
604         h = f.bamHdr;
605         read = f.readRecord!Bam1();
606         assert(fromStringz(bam_get_qname(read)) == "HS18_09653:4:1315:19857:61712");
607     }
608 }
609 
610 debug(dhtslib_unittest) unittest
611 {
612     hts_log_info(__FUNCTION__, "Testing HtslibFile BCF reading and writing");
613     import std.path:buildPath,dirName;
614     auto fn = buildPath(dirName(dirName(dirName(dirName(__FILE__)))),"htslib","test","tabix","vcf_file.vcf.gz");
615     {
616         auto f = HtslibFile(fn);
617         f.loadHeader;
618         auto h = f.bcfHdr;
619         auto read = f.readRecord!Bcf1();
620 
621         f = HtslibFile("/tmp/htslibfile.test.vcf", HtslibFileWriteMode.Vcf);
622         f.setHeader(h);
623         f.writeHeader;
624         f.write(read);
625         
626         f = HtslibFile("/tmp/htslibfile.test.bcf", HtslibFileWriteMode.Bcf);
627         f.setHeader(h);
628         f.writeHeader;
629         f.write(read);
630 
631         f = HtslibFile("/tmp/htslibfile.test.ubcf", HtslibFileWriteMode.UncompressedBcf);
632         f.setHeader(h);
633         f.writeHeader;
634         f.write(read);
635 
636         f = HtslibFile("/tmp/htslibfile.test.vcf.gz", HtslibFileWriteMode.GzippedVcf);
637         f.setHeader(h);
638         f.writeHeader;
639         f.write(read);
640 
641         f = HtslibFile("/tmp/htslibfile.test.vcf.bgz", HtslibFileWriteMode.BgzippedVcf);
642         f.setHeader(h);
643         f.writeHeader;
644         f.write(read);
645 
646     }
647     {
648 
649         auto f = HtslibFile("/tmp/htslibfile.test.vcf");
650         f.loadHeader;
651         auto h = f.bcfHdr;
652         auto read = f.readRecord!Bcf1();
653         assert(read.pos == 3000149);
654         
655         f = HtslibFile("/tmp/htslibfile.test.bcf");
656         f.loadHeader;
657         h = f.bcfHdr;
658         read = f.readRecord!Bcf1();
659         assert(read.pos == 3000149);
660 
661         f = HtslibFile("/tmp/htslibfile.test.ubcf");
662         f.loadHeader;
663         h = f.bcfHdr;
664         read = f.readRecord!Bcf1();
665         assert(read.pos == 3000149);
666 
667         f = HtslibFile("/tmp/htslibfile.test.vcf.gz");
668         f.loadHeader;
669         h = f.bcfHdr;
670         read = f.readRecord!Bcf1();
671         assert(read.pos == 3000149);
672 
673         f = HtslibFile("/tmp/htslibfile.test.vcf.bgz");
674         f.loadHeader;
675         h = f.bcfHdr;
676         read = f.readRecord!Bcf1();
677         assert(read.pos == 3000149);
678     }
679 }