Repository: bionumpy/bionumpy (language: python)

Issue:
Parse vcf header
Implementation of VCFHeader (take string as input, no need to split by new line before passing into the function), including a simple test (INFO is not tested yet)



Discussion / hints:

Fix patch:
diff --git a/benchmarks/config.yml b/benchmarks/config.yml
index ceea95f6b..410928826 100644
--- a/benchmarks/config.yml
+++ b/benchmarks/config.yml
@@ -87,7 +87,7 @@ analysis:
 
   vcf_filtering:
     name: "Filtering VCF on allele frequency"
-    tools: ["bionumpy", "bcftools", "pyvcf"]
+    tools: ["bionumpy", "bcftools"]
     validate_equal: false
     output_file_ending: ".filtered.csv"
     runs:
diff --git a/benchmarks/rules/interval_intersection.smk b/benchmarks/rules/interval_intersection.smk
index d271acfe8..3489d33b0 100644
--- a/benchmarks/rules/interval_intersection.smk
+++ b/benchmarks/rules/interval_intersection.smk
@@ -23,8 +23,9 @@ rule unique_intersect_bionumpy:
         "results/bionumpy/unique_intersect/{a}-vs-{b}.bed"
     benchmark:
         "benchmarks/unique_intersect/bionumpy/{a}-vs-{b}.txt"
-    script:
-        "../scripts/bionumpy_unique_intersect.py"
+    shell:
+        "python ../scripts/unique_intersect_example.py {input} {output}"
+        #        "../scripts/bionumpy_unique_intersect.py"
 
 
 rule intersect_bionumpy:
diff --git a/benchmarks/scripts/bionumpy_unique_intersect.py b/benchmarks/scripts/bionumpy_unique_intersect.py
index cf169f764..a5a9e33b8 100644
--- a/benchmarks/scripts/bionumpy_unique_intersect.py
+++ b/benchmarks/scripts/bionumpy_unique_intersect.py
@@ -33,9 +33,6 @@ class Snakefile:
 a = global_offset.from_local_interval(a)
 b = global_offset.from_local_interval(b)
 
-#chromosome_encoding = bnp.encodings.string_encodings.StringEncoding(chrom_sizes.name)
-#a = dataclasses.replace(a, chromosome=chromosome_encoding.encode(a.chromosome))
-# b = dataclasses.replace(b, chromosome=chromosome_encoding.encode(b.chromosome))
 result = bnp.arithmetics.unique_intersect(a, b, sum(chrom_sizes.size))
 result = global_offset.to_local_interval(result)
 with bnp.open(snakemake.output[0], "w", buffer_type=bnp.Bed6Buffer) as f:
diff --git a/bionumpy/__init__.py b/bionumpy/__init__.py
index cd819ed6a..d65728ffd 100644
--- a/bionumpy/__init__.py
+++ b/bionumpy/__init__.py
@@ -8,17 +8,17 @@
 
 from .io import (count_entries, open_indexed, MultiLineFastaBuffer, bnp_open,
                  TwoLineFastaBuffer, FastQBuffer, Bed6Buffer, NarrowPeakBuffer,
-                 BedBuffer, VCFBuffer, PhasedVCFMatrixBuffer, VCFMatrixBuffer,
-                 GfaSequenceBuffer, get_bufferclass_for_datatype)
+                 BedBuffer, GfaSequenceBuffer, get_bufferclass_for_datatype)
 from .encodings.alphabet_encoding import (DNAEncoding, RNAENcoding, AminoAcidEncoding)
 from .encoded_array import EncodedArray, EncodedRaggedArray, as_encoded_array, OneToOneEncoding, BaseEncoding, change_encoding
 from .sequence import (get_kmers, get_minimizers, get_motif_scores, count_encoded, match_string)
 from .streams import mean, bincount, histogram, streamable, quantile, MultiStream, groupby
 from .datatypes import SAMEntry, GFFEntry, Bed6, Interval, LocationEntry
+from .bnpdataclass import replace
 from .io.strops import str_equal
 from .util.cli import run_as_commandline
 from .computation_graph import compute
-from . import simulate
+# from . import simulate, VCFBuffer, VCFMatrixBuffer, PhasedVCFMatrixBuffer
 from . import arithmetics
 from . import alignments
 from .genomic_data import Genome, GenomicArray, GenomicIntervals
@@ -36,8 +36,7 @@
 __all__ = ["EncodedArray", "EncodedRaggedArray",
            "KmerEncoder", "Minimizers", "PositionWeightMatrix", "mean",
            "bincount", "streamable", "histogram", "count_entries", "quantile",
-           "BedBuffer", "VCFBuffer", "PhasedVCFMatrixBuffer", "VCFMatrixBuffer",
-           "GfaSequenceBuffer",
+           "BedBuffer", "GfaSequenceBuffer",
            "TwoLineFastaBuffer", "FastQBuffer", "open_indexed", "groupby",
            "SAMBuffer", "GFFBuffer", "Bed6Buffer", "MultiLineFastaBuffer",
            "count_encoded", "DNAEncoding", "RNAENcoding", "AminoAcidEncoding"]
diff --git a/bionumpy/bnpdataclass/__init__.py b/bionumpy/bnpdataclass/__init__.py
index 9833f02f9..fa764e760 100644
--- a/bionumpy/bnpdataclass/__init__.py
+++ b/bionumpy/bnpdataclass/__init__.py
@@ -1,9 +1,6 @@
 import dataclasses
 from .bnpdataclass import bnpdataclass, BNPDataClass, make_dataclass
+from .bnpdataclassfunction import replace
 
-def replace(obj, **kwargs):
-    if hasattr(obj, '__replace__'):
-        return obj.__replace__(**kwargs)
-    return dataclasses.replace(obj, **kwargs)
 
-__all__ = ['bnpdataclass', "BNPDataClass", "make_dataclass"]
+__all__ = ['bnpdataclass', "BNPDataClass", "make_dataclass", 'replace']
diff --git a/bionumpy/bnpdataclass/bnpdataclass.py b/bionumpy/bnpdataclass/bnpdataclass.py
index ba03f0657..0f0be24a4 100644
--- a/bionumpy/bnpdataclass/bnpdataclass.py
+++ b/bionumpy/bnpdataclass/bnpdataclass.py
@@ -173,7 +173,7 @@ def _implicit_format_conversion(cls, obj: npdataclass):
             """
             for field in dataclasses.fields(obj):
                 pre_val = getattr(obj, field.name)
-                if field.type in (int, float):
+                if field.type in (int, float, bool):
                     val = np.asanyarray(pre_val)
                 elif field.type == str:
                     assert isinstance(pre_val, (str, list, EncodedArray, EncodedRaggedArray, RaggedArray, np.ndarray)) or hasattr(pre_val, 'to_numpy'), (field, pre_val, type(pre_val))
@@ -195,10 +195,11 @@ def _implicit_format_conversion(cls, obj: npdataclass):
                     else:
                         val = pre_val
                 elif issubclass(field.type, BNPDataClass):
-                    assert isinstance(pre_val, (field.type, field.type._single_entry))
+                    assert isinstance(pre_val, (field.type, field.type._single_entry)), (field.type, type(pre_val))
                     val = pre_val
                 else:
                     assert False, field.type
+
                 setattr(obj, field.name, val)
 
     NewClass.__name__ = base_class.__name__
diff --git a/bionumpy/bnpdataclass/bnpdataclassfunction.py b/bionumpy/bnpdataclass/bnpdataclassfunction.py
index e5bbe939f..6cb188c69 100644
--- a/bionumpy/bnpdataclass/bnpdataclassfunction.py
+++ b/bionumpy/bnpdataclass/bnpdataclassfunction.py
@@ -10,6 +10,10 @@ def __call__(self, func):
         def new_func(data_object, *args, **kwargs):
             pass
 
+def replace(obj, **kwargs):
+    if hasattr(obj, '__replace__'):
+        return obj.__replace__(**kwargs)
+    return dataclasses.replace(obj, **kwargs)
 
 def apply_to_npdataclass(attribute_name):
     def decorator(func):
@@ -17,7 +21,7 @@ def new_func(np_dataclass, *args, **kwargs):
             if not isinstance(np_dataclass, BNPDataClass):
                 return func(np_dataclass, *args, **kwargs)
             result = func(getattr(np_dataclass, attribute_name), *args, **kwargs)
-            return dataclasses.replace(np_dataclass, **{attribute_name: result})
+            return replace(np_dataclass, **{attribute_name: result})
         return new_func
     return decorator
 
diff --git a/bionumpy/bnpdataclass/lazybnpdataclass.py b/bionumpy/bnpdataclass/lazybnpdataclass.py
new file mode 100644
index 000000000..17791d2e4
--- /dev/null
+++ b/bionumpy/bnpdataclass/lazybnpdataclass.py
@@ -0,0 +1,176 @@
+import dataclasses
+import types
+from functools import lru_cache
+from numbers import Number
+
+import numpy as np
+
+from bionumpy.io.dump_csv import get_column, join_columns
+from bionumpy.io.exceptions import FormatException
+
+
+# from bionumpy import EncodedRaggedArray
+
+
+def translate_types(input_type):
+    if input_type == str:
+        return EncodedRaggedArray
+    elif input_type == int:
+        return np.ndarray
+
+
+def buffer_backed_bnp(old_cls):
+    cls = types.new_class(old_cls.__name__, old_cls.__bases__, {})
+    for i, (var_name, var_type) in enumerate(old_cls.__annotations__.items()):
+        setattr(cls, var_name, BufferBackedDescriptor(i, var_type))
+        setattr(cls, '__init__', lambda self, buffer: setattr(self, '_buffer', buffer))
+    return cls
+
+
+class BufferBackedDescriptor:
+    '''
+    This class is made to access and parse parts of a text buffer lazily.v
+    '''
+
+    def __init__(self, buffer, index, dtype):
+        self._buffer = buffer
+        self._index = index
+        self._dtype = dtype
+
+    def __get__(self, obj, objtype):
+        return self._dtype(obj._buffer.get_field_by_number(self._index, self._dtype))
+
+
+class LazyBNPDataClass:
+    pass
+
+class BaseClass:
+    def __init__(self, buffer):
+        self._buffer = buffer
+
+    def __getattr__(self, var_name):
+        if var_name in self._buffer:
+            return self._buffer[var_name]
+        return super().__getattr__(var_name)
+
+
+class ItemGetter:
+    def __init__(self, buffer: 'FileBuffer', dataclass: dataclasses.dataclass, start_line=0):
+        self._buffer = buffer
+        self._dataclass = dataclass
+        self._field_dict = {field.name: (i, field.type) for i, field in enumerate(dataclasses.fields(dataclass))}
+        self._buffer.validate_if_not()
+        self._start_line = start_line
+
+    @lru_cache()
+    def n_entries(self):
+        return self._buffer.count_entries()
+
+    def __call__(self, name):
+        try:
+            return self._buffer.get_field_by_number(self._field_dict[name][0], self._field_dict[name][1])
+        except FormatException as e:
+            e.line_number += self._start_line
+            raise e
+
+    def __getitem__(self, idx):
+        return self.__class__(self._buffer[idx], self._dataclass)
+
+    @property
+    def buffer(self):
+        return self._buffer
+
+
+def create_lazy_class(dataclass, header=None):
+    field_names = [field.name for field in dataclasses.fields(dataclass)]
+
+
+    class NewClass(dataclass, LazyBNPDataClass):
+        def __init__(self, item_getter, set_values=None, computed_values=None):
+            self._itemgetter = item_getter
+            self._set_values = set_values or {}
+            self._computed_values = computed_values or {}
+            self._computed = False
+            self._data = None
+            if header is not None:
+                self.set_context('header', header)
+
+        def __len__(self):
+            return self._itemgetter.n_entries()
+
+        def __repr__(self):
+            return self[:10].get_data_object().__repr__().replace('with 10 entries', f'with {len(self)} entries')
+
+        def __str__(self):
+            return self[:10].get_data_object().__str__().replace('with 10 entries', f'with {len(self)} entries')
+
+        def __getattr__(self, var_name):
+            if var_name in self._set_values:
+                return self._set_values[var_name]
+            if var_name in field_names:
+                if var_name not in self._computed_values:
+                    self._computed_values[var_name] = self._itemgetter(var_name)
+                return self._computed_values[var_name]
+            return super().__getattr__(var_name)
+
+        def __setattr__(self, key, value):
+            if key in ['_itemgetter', '_set_values', '_computed', '_data', '_computed_values']:
+                return super().__setattr__(key, value)
+            self._set_values[key] = value
+            if key in self._computed_values:
+                del self._computed_values[key]
+
+        def __getitem__(self, idx):
+            if isinstance(idx, Number):
+                idx = [idx]
+                return self[idx].get_data_object()[0]
+            new_dict = {key: value[idx] for key, value in self._set_values.items()}
+            new_computed = {key: value[idx] for key, value in self._computed_values.items()}
+            return self.__class__(self._itemgetter[idx], new_dict, new_computed)
+
+        def __replace__(self, **kwargs):
+            new_dict = {key: value for key, value in self._set_values.items()}
+            new_dict.update(kwargs)
+            return self.__class__(self._itemgetter, new_dict)
+
+        def __iter__(self):
+            return iter(self.get_data_object())
+
+        def get_data_object(self):
+            if not self._computed:
+                self._data = dataclass(*(getattr(self, field.name) for field in dataclasses.fields(dataclass)))
+                self._computed = True
+            return self._data
+
+        @classmethod
+        def empty(cls):
+            return dataclass.empty()
+
+        def __array_function__(self, func, types, args, kwargs):
+            assert all(issubclass(t, LazyBNPDataClass) for t in types), types
+            if func == np.concatenate:
+                objects = [a.get_data_object() for a in args[0]]
+                args = (objects, ) + args[1:]
+                return func(*args, **kwargs)
+            return NotImplemented
+
+        def get_buffer(self, buffer_class=None):
+            if buffer_class is None:
+                buffer_class = self._itemgetter.buffer.__class__
+            if not hasattr(self._itemgetter.buffer, 'get_field_range_as_text') or hasattr(self._itemgetter.buffer, 'SKIP_LAZY') or hasattr(buffer_class, 'SKIP_LAZY'):
+                return self._itemgetter.buffer.from_data(self.get_data_object())
+            columns = []
+            if not self._set_values and self._itemgetter.buffer.__class__ == buffer_class:
+                return self._itemgetter.buffer.data.ravel()
+            for i, field in enumerate(dataclasses.fields(dataclass)):
+                if field.name in self._set_values:
+                    columns.append(get_column(self._set_values[field.name], field.type))
+                else:
+                    columns.append(self._itemgetter.buffer.get_field_range_as_text(i, i+1))
+            return buffer_class.join_fields(columns)
+            # return join_columns(columns, self._itemgetter.buffer.DELIMITER).ravel()
+
+
+    NewClass.__name__ = dataclass.__name__
+    NewClass.__qualname__ = dataclass.__qualname__
+    return NewClass
diff --git a/bionumpy/datatypes/__init__.py b/bionumpy/datatypes/__init__.py
index bb809e9e6..377a1bb53 100644
--- a/bionumpy/datatypes/__init__.py
+++ b/bionumpy/datatypes/__init__.py
@@ -2,7 +2,7 @@
 from typing import List
 from ..encodings import (CigarOpEncoding, BamEncoding, QualityEncoding,
                          CigarEncoding, StrandEncoding)
-from ..encodings.vcf_encoding import PhasedGenotypeRowEncoding, GenotypeRowEncoding
+from ..encodings.vcf_encoding import PhasedGenotypeRowEncoding, GenotypeRowEncoding, PhasedHaplotypeRowEncoding
 from ..bnpdataclass import bnpdataclass
 from .gtf import GFFEntry, GFFExonEntry, GFFGeneEntry, GFFTranscriptEntry, GTFEntry
 
@@ -117,6 +117,11 @@ class PhasedVCFGenotypeEntry(VCFEntry):
     genotypes: PhasedGenotypeRowEncoding
 
 
+@bnpdataclass
+class PhasedVCFHaplotypeEntry(VCFEntry):
+    genotypes: PhasedHaplotypeRowEncoding
+
+
 class SNP(Variant):
     pass
 
diff --git a/bionumpy/encoded_array.py b/bionumpy/encoded_array.py
index db04b688d..69e485670 100644
--- a/bionumpy/encoded_array.py
+++ b/bionumpy/encoded_array.py
@@ -11,6 +11,7 @@
 import numpy as np
 from abc import abstractmethod
 
+
 class Encoding:
     @abstractmethod
     def encode(self, *args, **kwargs):
@@ -22,7 +23,7 @@ def get_labels(self):
 
     def __call__(self, *args, **kwargs):
         return self.encode(*args, **kwargs)
-    
+
     def is_base_encoding(self):
         return False
 
@@ -133,7 +134,6 @@ def __eq__(self, other):
         return isinstance(other, ASCIIEncoding)
 
 
-
 class NumericEncoding(OneToOneEncoding):
     def is_numeric(self):
         return True
@@ -145,6 +145,7 @@ def is_numeric(self):
 def get_base_encodings():
     return [BaseEncoding]  # TODO: add other encodings
 
+
 class EncodingException(Exception):
     pass
 
@@ -163,12 +164,12 @@ def __init__(self, *args, **kwargs):
     def __repr__(self) -> str:
         if len(self) == 0:
             return ''
-        if self.size>1000:
+        if self.size > 1000:
             rows = [str(row) for row in self[:5]]
         else:
             rows = [f"{row}" for row in self]
-        encoding_info = f", {self.encoding}" if not self.encoding.is_base_encoding()  else ""
-        indent = " "*len("encoded_ragged_array([")
+        encoding_info = f", {self.encoding}" if not self.encoding.is_base_encoding() else ""
+        indent = " " * len("encoded_ragged_array([")
         quotes = "'" if self.encoding.is_one_to_one_encoding() else ""
         lines = [f"{indent}{quotes}{row}{quotes}," for row in rows]
         lines[0] = lines[0].replace(indent, "encoded_ragged_array([")
@@ -189,7 +190,7 @@ def __array_ufunc__(self, ufunc, method, *inputs, **kwargs):
         """ Convert any data to `EncodedArray` before calling the `ufunc` on them """
         assert isinstance(self.ravel(), EncodedArray), self.ravel()
         inputs = [as_encoded_array(i, self.ravel().encoding).raw() for i in inputs]
-        #inputs = _parse_ufunc_inputs(inputs, self.ravel().encoding)
+        # inputs = _parse_ufunc_inputs(inputs, self.ravel().encoding)
         kwargs = {key: as_encoded_array(val, self.ravel().encoding).raw() for key, val in kwargs.items()}
 
         ret = super().__array_ufunc__(ufunc, method, *inputs, **kwargs)
@@ -233,10 +234,10 @@ def __init__(self, data: np.ndarray, encoding):
             data = data.data
         self.encoding = encoding
         dtype = None if hasattr(data, "dtype") else np.uint8
-        #self.data = np.asarray(data, dtype=dtype).view(NPSArray)
+        # self.data = np.asarray(data, dtype=dtype).view(NPSArray)
         self.data = np.asarray(data, dtype=dtype)
         self.data = get_NPSArray(self.data)
-        #assert isinstance(self.data, np.ndarray)
+        # assert isinstance(self.data, np.ndarray)
 
     def copy(self):
         return self.__class__(self.data.copy(), self.encoding)
@@ -307,11 +308,11 @@ def __str__(self) -> str:
             if n_dims == 0:
                 return self.encoding.to_string(self.data)
             elif n_dims == 1:
-                data_to_show = data_to_show[0:20]
+                data_to_show = data_to_show  # [0:20]
             elif n_dims == 2:
                 # show first columns and rows
-                #raise NotImplemented("Str for n_dims=2 not implemented")
-                return self.encoding.to_string(self.data[0:10, 0:10]) + "...."
+                # raise NotImplemented("Str for n_dims=2 not implemented")
+                return self.encoding.to_string(self.data[0:10])  # + "...."
                 # data_to_show = None
             text = "[" + ", ".join(self.encoding.to_string(e).strip() for e in data_to_show) + "]"
             return text
@@ -327,13 +328,16 @@ def __str__(self) -> str:
             if len(self.data.shape) == 0:
                 return chr(int(text))
             if len(self.data.shape) == 1:
-                return "".join(chr(n) for n in text[:20]) + "..."*(len(text)>20)
+                return "".join(chr(n) for n in text)  # + "..."*(len(text)>20)
 
             a = np.array([str(self.__class__(seq, self.encoding))
                           for seq in self.data.reshape(-1, self.data.shape[-1])]
                          ).reshape(self.data.shape[:-1])[:20]
             return str(a)
 
+    def __hash__(self):
+        if len(self.shape) <= 1:
+            return hash(self.to_string())
 
     def __getitem__(self, idx) -> "EncodedArray":
         """Delegate the indexing to the underlying numpy array
@@ -403,7 +407,7 @@ def __array_function__(self, func, types, args, kwargs):
         if func == np.concatenate:
             return self.__class__(func([e.data for e in args[0]]), self.encoding)
         if func == np.where:
-            return self.__class__(func(args[0], args[1].data, args[2].data), encoding = self.encoding)
+            return self.__class__(func(args[0], args[1].data, args[2].data), encoding=self.encoding)
         if func == np.zeros_like:
             return self.__class__(func(args[0].data, *args[1:], **kwargs), encoding=self.encoding)
         if func == np.append:
@@ -417,7 +421,8 @@ def __array_function__(self, func, types, args, kwargs):
         if func == np.full_like:
             return full_like(*args, **kwargs)
         if func == np.insert:
-            return self.__class__(func(args[0].data, args[1], args[2].data, *args[3:], **kwargs), encoding = self.encoding)
+            return self.__class__(func(args[0].data, args[1], args[2].data, *args[3:], **kwargs),
+                                  encoding=self.encoding)
         elif func in (np.lib.stride_tricks.sliding_window_view, np.lib.stride_tricks.as_strided):
             return self.__class__(func(args[0].data, *args[1:], **kwargs), self.encoding)
 
@@ -441,7 +446,9 @@ def as_strided(self, *args, **kwargs):
 
 def _parse_ufunc_inputs(inputs, target_encoding):
     for a in inputs:
-        assert isinstance(a, (str, list, EncodedArray, EncodedRaggedArray)), "%s is not str, list, EncodedArray or EncodedRaggedArray" % type(a)
+        assert isinstance(a, (
+        str, list, EncodedArray, EncodedRaggedArray)), "%s is not str, list, EncodedArray or EncodedRaggedArray" % type(
+            a)
         yield as_encoded_array(a, target_encoding).raw()
         # if isinstance(a, str):
         #     a = target_encoding.encode(a)
@@ -480,6 +487,7 @@ def new_func(array, *args, **kwargs):
                     raise Exception
         return result
 
+
 def as_encoded_array(s: object, target_encoding: "Encoding" = None) -> EncodedArray:
     """Main function used to create encoded arrays from e.g. strings orl lists.
     Can be called with already encoded arrays, and will then do nothing.
@@ -519,15 +527,16 @@ def as_encoded_array(s: object, target_encoding: "Encoding" = None) -> EncodedAr
                 if hasattr(s.encoding, 'get_alphabet') and hasattr(target_encoding, 'get_alphabet'):
                     m = s.raw().max()
                     if (s.encoding.get_alphabet()[:m] == target_encoding.get_alphabet()[:m]):
-                        if not m<len(target_encoding.get_alphabet()):
-                            raise EncodingException(f"Trying to encode already encoded array with encoding {s.encoding} to encoding {target_encoding}. However {s} contains {EncodedArray(m, s.encoding)}")
+                        if not m < len(target_encoding.get_alphabet()):
+                            raise EncodingException(
+                                f"Trying to encode already encoded array with encoding {s.encoding} to encoding {target_encoding}. However {s} contains {EncodedArray(m, s.encoding)}")
                         if isinstance(s, EncodedArray):
                             return s.__class__(s.raw(), target_encoding)
                         elif isinstance(s, EncodedRaggedArray):
                             return s.__class__(EncodedArray(s.ravel().raw(), target_encoding), s.shape)
                 raise EncodingException("Trying to encode already encoded array with encoding %s to encoding %s. "
                                         "This is not supported. Use the change_encoding function." % (
-                    s.encoding, target_encoding))
+                                            s.encoding, target_encoding))
     elif target_encoding is None:
         target_encoding = BaseEncoding
 
@@ -544,6 +553,7 @@ def as_encoded_array(s: object, target_encoding: "Encoding" = None) -> EncodedAr
 
     return target_encoding.encode(s)
 
+
 def full_like(a, fill_value, dtype=None, order='K', subok=True, shape=None):
     assert dtype is None
     assert order == 'K'
@@ -551,6 +561,7 @@ def full_like(a, fill_value, dtype=None, order='K', subok=True, shape=None):
     fill_value = a.encoding.encode(fill_value)
     return EncodedArray(np.full_like(a.raw(), fill_value, shape=shape), a.encoding)
 
+
 def from_encoded_array(encoded_array: EncodedArray) -> str:
     """Convert data in an `EncodedArray`/`EncodedRaggedArray into `str`/`List[str]`
 
diff --git a/bionumpy/encodings/__init__.py b/bionumpy/encodings/__init__.py
index 17414662a..3287caf5c 100644
--- a/bionumpy/encodings/__init__.py
+++ b/bionumpy/encodings/__init__.py
@@ -1,4 +1,3 @@
-import numpy as np
 from ..encoded_array import BaseEncoding, Encoding, NumericEncoding
 from .alphabet_encoding import (AlphabetEncoding, DNAEncoding, RNAENcoding,
                                 AminoAcidEncoding, ACGTnEncoding,
@@ -6,7 +5,7 @@
 
 
 __all__ = ["BaseEncoding", "Encoding",
-           "AlphabetEncoding", "ACTGEncoding", "AminoAcidEncoding"]# , "ACTGTwoBitEncoding"]
+           "AlphabetEncoding",  "AminoAcidEncoding"]
 
 
 # class StrandEncoding(Encoding):
diff --git a/bionumpy/encodings/alphabet_encoding.py b/bionumpy/encodings/alphabet_encoding.py
index 77a25aee0..3a578859b 100644
--- a/bionumpy/encodings/alphabet_encoding.py
+++ b/bionumpy/encodings/alphabet_encoding.py
@@ -78,7 +78,7 @@ def __hash__(self):
 DNAEncoding = ACGTEncoding
 ACUGEncoding = AlphabetEncoding("ACUG")
 RNAENcoding = ACUGEncoding
-AminoAcidEncoding = AlphabetEncoding('ACDEFGHIKLMNPQRSTVWY')
+AminoAcidEncoding = AlphabetEncoding('ACDEFGHIKLMNPQRSTVWY*')
 BamEncoding = AlphabetEncoding("=ACMGRSVTWYHKDBN")
 CigarOpEncoding = AlphabetEncoding("MIDNSHP=X")
 
diff --git a/bionumpy/encodings/integer_encoding.py b/bionumpy/encodings/integer_encoding.py
new file mode 100644
index 000000000..9594254ce
--- /dev/null
+++ b/bionumpy/encodings/integer_encoding.py
@@ -0,0 +1,15 @@
+import numpy as np
+
+
+class IntegerEncoding:
+    def __init__(self, nan_symbol: str = None):
+        self._nan_symbol = nan_symbol
+
+    def _encode(self, data: np.ndarray):
+        encoded = data-ord('0')
+        if self._nan_symbol is not None:
+            encoded[data==ord(self._nan_symbol)] = np.nan
+        return encoded
+
+    def _decode(self, data: np.ndarray):
+        return
diff --git a/bionumpy/encodings/vcf_encoding.py b/bionumpy/encodings/vcf_encoding.py
index 3893e9f22..4b0da23ab 100644
--- a/bionumpy/encodings/vcf_encoding.py
+++ b/bionumpy/encodings/vcf_encoding.py
@@ -84,6 +84,47 @@ def to_string(self, e):
     def __repr__(self):
         return "GenotypeRowEncoding"
 
+
+class GenotypeBuffer:
+    _lookup = np.zeros(256, dtype=np.uint8)
+    _lookup[[ord(c) for c in ('0', '1', '.')]] = np.array([0, 1, np.nan])
+
+    def __init__(self, data, field_starts):
+        self._data = data
+        self._field_starts = field_starts
+        self._genotype = None
+        self._phasing_status = None
+
+    def __getitem__(self, index):
+         return self.__class__(self._data, self._field_starts[index])
+
+    @property
+    def genotype(self):
+        if self._genotype is None:
+            self._genotype = self._lookup[self._field_starts[..., None]+[0, 2]]
+        return self._genotype
+
+    @property
+    def phasing_status(self):
+        if self._phasing_status is None:
+            self._phasing_status = self._lookup
+
+        return
+
+
+    def _preprocess_data_for_encoding(self, genotype_rows):
+        # split the row of genotype data
+        from ..io.strops import split, replace_inplace
+        data = genotype_rows.ravel()
+        # hack because the row sometime ends with \n and sometimes with \t
+        replace_inplace(data, "\n", "\t")
+        indices = np.flatnonzero(data == "\t")
+        indices = np.insert(indices, 0, -1)
+        return data[indices[:-1, np.newaxis] + np.array([1, 2, 3])]
+        #data = split(data.ravel(), "\t")[:-1, 0:3]  # don't include last element which is empty
+        #return data
+
+
 class _PhasedGenotypeRowEncoding(_GenotypeRowEncoding):
     """Encoding that can be used when all records are phased
      and there is no missing data, i.e. every genotype is
@@ -110,5 +151,25 @@ def __repr__(self):
         return "PhasedGenotypeRowEncoding"
 
 
+class _PhasedHaplotypeRowEncoding(_GenotypeRowEncoding):
+    """Encoding that encodes each haplotype (not two haplotypes together as _PhasdGenotypeRowEncoding"""
+    _alleles = [str(i) for i in range(10)] + ["."]
+    _alphabet = _alleles
+    _reverse_alphabet_lookup = np.array([ord(c) for c in _alphabet], dtype=np.uint8)
+    _alphabet_lookup = np.zeros(256, dtype=np.uint8)
+    _alphabet_lookup[_reverse_alphabet_lookup] = np.arange(len(_reverse_alphabet_lookup))
+
+    def encode(self, genotype_rows):
+        data = self._preprocess_data_for_encoding(genotype_rows)
+        n_rows = len(genotype_rows)
+        first_haplotypes = self._alphabet_lookup[data[:, 0].raw()]
+        second_haplotypes = self._alphabet_lookup[data[:, 2].raw()]
+        out = np.zeros(len(first_haplotypes)*2, dtype=np.int8)
+        out[::2] = first_haplotypes
+        out[1::2] = second_haplotypes
+        return out.reshape(n_rows, len(out)//n_rows)
+
+
 PhasedGenotypeRowEncoding = _PhasedGenotypeRowEncoding()
+PhasedHaplotypeRowEncoding = _PhasedHaplotypeRowEncoding()
 GenotypeRowEncoding = _GenotypeRowEncoding()
\ No newline at end of file
diff --git a/bionumpy/genomic_data/__init__.py b/bionumpy/genomic_data/__init__.py
index c5410a9d7..78058f876 100644
--- a/bionumpy/genomic_data/__init__.py
+++ b/bionumpy/genomic_data/__init__.py
@@ -4,3 +4,4 @@
 from .annotation import GenomicAnnotation
 from .genomic_sequence import GenomicSequence
 from .genome_context import GenomeContext, keep_all, ignore_underscores
+from .binned_genome import BinnedGenome
\ No newline at end of file
diff --git a/bionumpy/genomic_data/binned_genome.py b/bionumpy/genomic_data/binned_genome.py
new file mode 100644
index 000000000..800a3cf81
--- /dev/null
+++ b/bionumpy/genomic_data/binned_genome.py
@@ -0,0 +1,43 @@
+import numpy as np
+from bionumpy import LocationEntry, bnp_open
+from bionumpy.genomic_data.genome import Genome
+from bionumpy.genomic_data.genome_context_base import GenomeContextBase
+
+
+class BinnedGenome:
+    def __init__(self, genome_context: GenomeContextBase, bin_size: int = 1000):
+        self._genome_context = genome_context
+        self._bin_size = bin_size
+        chrom_sizes = np.array(list(genome_context.chrom_sizes.values()))
+        self._n_bins = (chrom_sizes+bin_size-1)//bin_size
+        self._bin_offsets = np.insert(np.cumsum(self._n_bins), 0, 0)
+        self._counts = np.zeros(self._bin_offsets[-1], dtype=int)
+
+    @classmethod
+    def from_file(cls, filename: str, bin_size: int = 1000):
+        genome = Genome.from_file(filename)
+        return cls(genome.get_genome_context(), bin_size)
+
+    @property
+    def genome_context(self):
+        return self._genome_context
+
+    def count(self, entries: LocationEntry, position_field='position'):
+        chrom_nrs = self._genome_context.encoding.encode(entries.chromosome).raw()
+        offsets = getattr(entries, position_field)//self._bin_size
+        bin_nr = self._bin_offsets[chrom_nrs] + offsets
+        self._counts += np.bincount(bin_nr, minlength=self._bin_offsets[-1])
+
+    def count_file(self, filename, position_field='position'):
+        chunks = bnp_open(filename, 'r').read_chunks()
+        for chunk in chunks:
+            self.count(chunk, position_field=position_field)
+
+    @property
+    def count_dict(self):
+        return {chrom: self._counts[self._bin_offsets[i]:self._bin_offsets[i+1]]
+                for i, chrom in enumerate(self._genome_context.chrom_sizes)}
+
+    def __getitem__(self, chromosome):
+        i = self._genome_context.encoding.encode(chromosome).raw()
+        return self._counts[self._bin_offsets[i]:self._bin_offsets[i+1]]
\ No newline at end of file
diff --git a/bionumpy/genomic_data/genome.py b/bionumpy/genomic_data/genome.py
index c3617ca00..d97bec419 100644
--- a/bionumpy/genomic_data/genome.py
+++ b/bionumpy/genomic_data/genome.py
@@ -194,11 +194,16 @@ def read_locations(self, filename: str, stranded: bool = False, stream: bool = F
         f = bnp_open(filename, buffer_type=buffer_type)
         data = f.read_chunks()
         if not stream:
-            data = np.concatenate(list(data))
+            data_list = list(data)
+            if len(data_list) == 0:
+                data = LocationEntry.empty()
+            else:
+                data = np.concatenate(data_list)
+
         return self.get_locations(data, has_numeric_chromosomes=has_numeric_chromosomes)
 
     def _mask_data_on_extra_chromosomes(self, data, chromosome_field_name='chromosome'):
-        if not isinstance(data, BNPDataClass):
+        if (not isinstance(data, BNPDataClass)) or len(data)==0:
             return data
         encoded_chromosomes = self._genome_context.encoding.encode(getattr(data, chromosome_field_name))
         data = replace(data, **{chromosome_field_name: encoded_chromosomes})
diff --git a/bionumpy/genomic_data/genome_context.py b/bionumpy/genomic_data/genome_context.py
index 0d9f059e9..6eb6cafb7 100644
--- a/bionumpy/genomic_data/genome_context.py
+++ b/bionumpy/genomic_data/genome_context.py
@@ -68,6 +68,8 @@ def mask_data(self, data, chromosome_field_name='chromosome'):
         encoded_chromosomes = as_encoded_array(getattr(data, chromosome_field_name), self.encoding)
         data = replace(data, **{chromosome_field_name: encoded_chromosomes})
         mask = self.is_included(encoded_chromosomes)
+        if np.all(mask):
+            return data
         return data[mask]
 
     @classmethod
diff --git a/bionumpy/genomic_data/genomic_intervals.py b/bionumpy/genomic_data/genomic_intervals.py
index 6b7c9b3f6..29d1c217e 100644
--- a/bionumpy/genomic_data/genomic_intervals.py
+++ b/bionumpy/genomic_data/genomic_intervals.py
@@ -1,4 +1,6 @@
 from abc import ABC, abstractmethod, abstractproperty, abstractclassmethod
+from collections import namedtuple
+
 import numpy as np
 from typing import List, Iterable, Tuple, Dict
 from ..bnpdataclass import BNPDataClass, replace, bnpdataclass
@@ -6,7 +8,7 @@
 from .genomic_track import GenomicArray, GenomicArrayNode
 from .genome_context_base import GenomeContextBase
 from ..datatypes import Interval, Bed6, StrandedInterval, LocationEntry, StrandedLocationEntry
-from ..arithmetics.intervals import get_pileup, merge_intervals, extend_to_size, clip, get_boolean_mask
+from ..arithmetics.intervals import get_pileup, merge_intervals, extend_to_size, clip, get_boolean_mask, RawInterval
 from ..computation_graph import StreamNode, Node, ComputationNode, compute
 import dataclasses
 
@@ -130,7 +132,7 @@ def data(self):
 
     def __replace__(self, **kwargs):
         kwargs = {self._field_dict[kw]: value for kw, value in kwargs.items()}
-        return self.__class__(dataclasses.replace(self._locations, **kwargs), self._genome_context, self._is_stranded, self._field_dict)
+        return self.__class__(replace(self._locations, **kwargs), self._genome_context, self._is_stranded, self._field_dict)
 
     @property
     def chromosome(self):
@@ -614,8 +616,9 @@ def get_mask(self) -> GenomicArray:
         GenomicArray
             Genomic mask
         """
-        go = self._genome_context.global_offset.from_local_interval(self._intervals)
-        global_mask = get_boolean_mask(go, self._genome_context.size)
+        I = RawInterval
+        starts, stops = self._genome_context.global_offset.start_ends_from_intervals(self._intervals)
+        global_mask = get_boolean_mask(I(starts, stops), self._genome_context.size)
         return GenomicArray.from_global_data(global_mask, self._genome_context)
 
     def clip(self) -> 'GenomicIntervalsFull':
@@ -632,7 +635,7 @@ def clip(self) -> 'GenomicIntervalsFull':
                        stop=np.minimum(chrom_sizes, self.stop))
 
     def __replace__(self, **kwargs):
-        return self.__class__(dataclasses.replace(self._intervals, **kwargs), self._genome_context, self._is_stranded)
+        return self.__class__(replace(self._intervals, **kwargs), self._genome_context, self._is_stranded)
 
     def compute(self):
         return self
@@ -763,8 +766,8 @@ def clip(self) -> 'GenomicIntervals':
         return self.__class__(ComputationNode(clip, [self._intervals_node, self._chrom_size_node]), self._genome_context)
 
     def __replace__(self, **kwargs):
-        return self.__class__(ComputationNode(dataclasses.replace, [self._intervals_node], kwargs), self._genome_context)
-        return self.__class__(dataclasses.replace(self._intervals, **kwargs), self._genome_context)
+        return self.__class__(ComputationNode(replace, [self._intervals_node], kwargs), self._genome_context)
+        # return self.__class__(dataclasses.replace(self._intervals, **kwargs), self._genome_context)
 
     def compute(self):
         chromosome, start, stop = compute((self.chromosome, self.start, self.stop))
diff --git a/bionumpy/genomic_data/geometry.py b/bionumpy/genomic_data/geometry.py
index ae71c6d7f..0263ead7e 100644
--- a/bionumpy/genomic_data/geometry.py
+++ b/bionumpy/genomic_data/geometry.py
@@ -25,8 +25,7 @@ def from_chrom_sizes(cls, chrom_sizes: ChromosomeSize):
         ----------
         chrom_sizes : ChromosomeSize
         """
-        return cls({chrom_size.name.to_string():chrom_size.size
-                    for chrom_size in chrom_sizes})
+        return cls({str(chrom_size.name): chrom_size.size for chrom_size in chrom_sizes})
 
     def chrom_size(self, chromsome: str) -> int:
         """Return the size of the given chromosome
diff --git a/bionumpy/genomic_data/global_offset.py b/bionumpy/genomic_data/global_offset.py
index af8c59d58..16843cb99 100644
--- a/bionumpy/genomic_data/global_offset.py
+++ b/bionumpy/genomic_data/global_offset.py
@@ -3,6 +3,7 @@
 from ..encodings.string_encodings import StringEncoding
 from ..encoded_array import EncodedArray, as_encoded_array
 from ..bnpdataclass import replace
+
 global_encoding = StringEncoding(["global"])
 
 
@@ -35,20 +36,32 @@ def get_size(self, seq_name):
         return self._sizes[seq_name.raw()]
 
     def from_local_coordinates(self, sequence_name, local_offset):
+        mask = local_offset >= self.get_size(sequence_name)
+        if np.any(np.atleast_1d(mask)):
+            raise Exception('Coordinate outside of reference:', local_offset[mask], self.get_size(sequence_name)[mask])
         return self.get_offset(sequence_name) + local_offset
 
     def to_local_interval(self, global_interval):
-        chromosome_idxs = np.searchsorted(self._offset, global_interval.start, side="right")-1
-        start = global_interval.start-self._offset[chromosome_idxs]
-        stop = global_interval.stop-self._offset[chromosome_idxs]
+        chromosome_idxs = np.searchsorted(self._offset, global_interval.start, side="right") - 1
+        start = global_interval.start - self._offset[chromosome_idxs]
+        stop = global_interval.stop - self._offset[chromosome_idxs]
         assert np.all(stop <= self._sizes[chromosome_idxs])
         chromosome = EncodedArray(chromosome_idxs, self._old_encoding)
-        return dataclasses.replace(global_interval,
-                                   chromosome=chromosome,
-                                   start=start,
-                                   stop=stop)
+        return replace(global_interval,
+                       chromosome=chromosome,
+                       start=start,
+                       stop=stop)
 
     def from_local_interval(self, interval, do_clip=False):
+        start_offsets, stop_offsets = self.start_ends_from_intervals(interval, do_clip)
+        zeros = EncodedArray(np.broadcast_to(np.array(0, dtype=np.uint8), (len(interval),)), global_encoding)
+        return replace(
+            interval,
+            chromosome=zeros,
+            start=start_offsets,
+            stop=stop_offsets)
+
+    def start_ends_from_intervals(self, interval, do_clip=False):
         chromosome = as_encoded_array(interval.chromosome, target_encoding=self._old_encoding)
         offsets = self.get_offset(chromosome)
         sizes = self.get_size(chromosome)
@@ -60,9 +73,6 @@ def from_local_interval(self, interval, do_clip=False):
             stop = np.minimum(stop, sizes)
         else:
             assert np.all(stop <= sizes)
-        return replace(
-            interval,
-            chromosome=EncodedArray(
-                np.broadcast_to(np.array(0,dtype=np.uint8), (len(interval), )), global_encoding),
-            start=interval.start+offsets,
-            stop=stop+offsets)
+        start_offsets = interval.start + offsets
+        stop_offsets = stop + offsets
+        return start_offsets, stop_offsets
diff --git a/bionumpy/io/__init__.py b/bionumpy/io/__init__.py
index 5b94730d9..75c9aae43 100644
--- a/bionumpy/io/__init__.py
+++ b/bionumpy/io/__init__.py
@@ -2,8 +2,8 @@
 from .indexed_files import open_indexed
 from .multiline_buffer import MultiLineFastaBuffer
 from .file_buffers import (TwoLineFastaBuffer, FastQBuffer)
-from .delimited_buffers import (BedBuffer, Bed6Buffer, VCFBuffer, PhasedVCFMatrixBuffer, VCFMatrixBuffer,
-                                GfaSequenceBuffer, get_bufferclass_for_datatype, NarrowPeakBuffer)
+from .delimited_buffers import (BedBuffer, Bed6Buffer, GfaSequenceBuffer, get_bufferclass_for_datatype, NarrowPeakBuffer)
+from .vcf_buffers import VCFBuffer, VCFMatrixBuffer, PhasedVCFMatrixBuffer
 from .indexed_fasta import IndexedFasta
 from .matrix_dump import read_matrix
 from .motifs import read_motif
diff --git a/bionumpy/io/delimited_buffers.py b/bionumpy/io/delimited_buffers.py
index 66a174134..4856fe878 100644
--- a/bionumpy/io/delimited_buffers.py
+++ b/bionumpy/io/delimited_buffers.py
@@ -2,25 +2,22 @@
 import logging
 import dataclasses
 from typing import List
-from npstructures import RaggedArray, ragged_slice
+from npstructures import RaggedArray, ragged_slice, RaggedShape
 from ..bnpdataclass import bnpdataclass, BNPDataClass
-from ..datatypes import (Interval, VCFGenotypeEntry,
-                         SequenceEntry, VCFEntry, Bed12, Bed6, BedGraph,
-                         GTFEntry, GFFEntry, SAMEntry, ChromosomeSize, NarrowPeak, PhasedVCFGenotypeEntry,
-                         GfaPath)
+from ..bnpdataclass.lazybnpdataclass import LazyBNPDataClass
+from ..datatypes import (Interval, SequenceEntry, Bed12, Bed6, BedGraph,
+                         GTFEntry, GFFEntry, SAMEntry, ChromosomeSize, NarrowPeak, GfaPath)
 from ..encoded_array import EncodedArray, EncodedRaggedArray
 from ..encoded_array import as_encoded_array
 from ..encodings import Encoding
 from ..encodings.exceptions import EncodingError
-from ..encodings.vcf_encoding import GenotypeRowEncoding, PhasedGenotypeRowEncoding
-from ..encodings.alphabet_encoding import get_alphabet_encodings, DigitEncoding
-from ..encoded_array import get_base_encodings, BaseEncoding
+from ..encodings.alphabet_encoding import DigitEncoding
+from ..encoded_array import BaseEncoding
 from ..util import is_subclass_or_instance
-from .file_buffers import FileBuffer, NEWLINE
+from .file_buffers import FileBuffer, NEWLINE, TextBufferExtractor
 from .strops import (
-    ints_to_strings, split, str_to_int, str_to_float,
-    int_lists_to_strings, float_to_strings)
-from .dump_csv import dump_csv
+    split, str_to_int, str_to_float)
+from .dump_csv import dump_csv, join_columns
 from .exceptions import FormatException
 import numpy as np
 
@@ -38,7 +35,7 @@ class DelimitedBuffer(FileBuffer):
     HAS_UNCOMMENTED_HEADER_LINE = False
     n_lines_per_entry = 1
 
-    def __init__(self, data: EncodedArray, new_lines: np.ndarray, delimiters: np.ndarray = None, header_data=None):
+    def __init__(self, data: EncodedArray, new_lines: np.ndarray = None, delimiters: np.ndarray = None, header_data=None):
         super().__init__(data, new_lines)
         if delimiters is None:
             delimiters = np.concatenate(
@@ -47,6 +44,35 @@ def __init__(self, data: EncodedArray, new_lines: np.ndarray, delimiters: np.nda
             delimiters.sort(kind="mergesort")
         self._delimiters = delimiters
         self._header_data = header_data
+        self.__buffer_extractor = None
+
+    @property
+    def _buffer_extractor(self):
+        if self.__buffer_extractor is None:
+            self.__buffer_extractor = self._get_buffer_extractor()
+        return self.__buffer_extractor
+
+    def _get_buffer_extractor(self) -> TextBufferExtractor:
+        self.validate_if_not()
+        starts = self._delimiters[:-1].reshape(-1, self._n_cols)+1
+        ends = self._delimiters[1:].reshape(-1, self._n_cols)
+        return TextBufferExtractor(self._data, starts, ends)
+
+    def __getitem__(self, idx):
+        self.validate_if_not()
+        cell_lens = np.diff(self._delimiters).reshape(-1, self._n_cols)[idx]
+        entries = self.entries[idx].ravel()
+        delimiters = np.insert(np.cumsum(cell_lens)-1, 0, -1)
+        new_lines = delimiters[self._n_cols:: self._n_cols]
+        return self.__class__(entries, new_lines, delimiters)
+
+    @property
+    def entries(self):
+        if not hasattr(self, "_entries"):
+            lengths = np.diff(self._new_lines)
+            lengths = np.insert(lengths, 0, self._new_lines[0]+1)
+            self._entries = EncodedRaggedArray(self._data, RaggedShape(lengths))
+        return self._entries
 
     @classmethod
     def from_raw_buffer(cls, chunk: np.ndarray, header_data=None) -> "DelimitedBuffer":
@@ -151,12 +177,19 @@ def get_text(self, col, fixed_length=True, keep_sep=False):
         else:
             return self._move_intervals_to_ragged_array(starts, ends)
 
+    @classmethod
+    def join_fields(cls, fields_list: List[EncodedRaggedArray]):
+        return join_columns(fields_list, cls.DELIMITER).ravel()
+
     def _col_starts(self, col):
         return self._delimiters[:-1].reshape(-1, self._n_cols)[:, col] + 1
 
     def _col_ends(self, col):
         return self._delimiters[1:].reshape(-1, self._n_cols)[:, col]
 
+    def get_field_range_as_text(self, *args, **kwargs):
+        return self.get_column_range_as_text(*args, **kwargs)
+
     def get_column_range_as_text(self, col_start, col_end, keep_sep=False):
         """Get multiple columns as text
 
@@ -170,7 +203,7 @@ def get_column_range_as_text(self, col_start, col_end, keep_sep=False):
             keep seperator at end
         """
         self.validate_if_not()
-        assert col_start < col_end
+        assert col_start < col_end, (col_start, col_end)
         assert col_start < self._n_cols, self._n_cols
         assert col_end <= self._n_cols
 
@@ -250,6 +283,8 @@ def from_data(cls, data: BNPDataClass) -> "DelimitedBuffer":
         data : bnpdataclass
             Data
         """
+        if isinstance(data, LazyBNPDataClass):
+            return cls.from_data(data.get_data_object())
         data_dict = [(field.type, getattr(data, field.name)) for field in dataclasses.fields(data)]
         return dump_csv(data_dict, cls.DELIMITER)
 
@@ -273,22 +308,7 @@ def get_data(self) -> BNPDataClass:
         columns = {}
         fields = dataclasses.fields(self.dataclass)
         for col_number, field in enumerate(fields):
-            if field.type is None:
-                col = None
-            elif field.type == str or is_subclass_or_instance(field.type, Encoding):
-                col = self.get_text(col_number, fixed_length=False)
-            elif field.type == int:
-                col = self.get_integers(col_number).ravel()
-            elif field.type == float:
-                col = self.get_floats(col_number).ravel()
-            elif field.type == -1:
-                col = self.get_integers(col_number).ravel() - 1
-            elif field.type == List[int]:
-                col = self.get_split_ints(col_number)
-            elif field.type == List[bool]:
-                col = self.get_split_ints(col_number, sep="").astype(bool)
-            else:
-                assert False, field
+            col = self._get_field_by_number(col_number, field.type)
             columns[field.name] = col
         n_entries = len(next(col for col in columns if col is not None))
         columns = {c: value if c is not None else np.empty((n_entries, 0))
@@ -297,6 +317,34 @@ def get_data(self) -> BNPDataClass:
         data.set_context("header", self._header_data)
         return data
 
+    def _get_field_by_number(self, col_number, field_type):
+        if field_type is None:
+            col = None
+        elif field_type == str or is_subclass_or_instance(field_type, Encoding):
+            col = self.get_text(col_number, fixed_length=False)
+        elif field_type == int:
+            col = self.get_integers(col_number).ravel()
+        elif field_type == float:
+            col = self.get_floats(col_number).ravel()
+        elif field_type == bool:
+            col = self.get_integers(col_number).ravel().astype(bool)
+        elif field_type == -1:
+            col = self.get_integers(col_number).ravel() - 1
+        elif field_type == List[int]:
+            col = self.get_split_ints(col_number)
+        elif field_type == List[bool]:
+            col = self.get_split_ints(col_number, sep="").astype(bool)
+        else:
+            assert False, field_type
+        return col
+
+    def get_field_by_number(self, field_nr: int, field_type: type=object):
+        self.validate_if_not()
+        if field_type is None:
+            field_type = dataclasses.fields(self.dataclass)[field_nr]
+        return self._get_field_by_number(
+            field_nr, field_type)
+
     def get_split_ints(self, col: int, sep: str = ",") -> RaggedArray:
         """Split a column of separated integers into a raggedarray
 
@@ -331,12 +379,15 @@ def count_entries(self) -> int:
 
 class GfaSequenceBuffer(DelimitedBuffer):
     dataclass = SequenceEntry
-
+    SKIP_LAZY = True
     def get_data(self):
         ids = self.get_text(1, fixed_length=False)
         sequences = self.get_text(col=2, fixed_length=False)
         return SequenceEntry(ids, sequences)
 
+    def get_field_by_number(self, field_nr: int, field_type: type=object):
+        return super().get_field_by_number(field_nr+1, field_type)
+
     @classmethod
     def from_data(cls, data: SequenceEntry) -> EncodedArray:
         return dump_csv([(str, as_encoded_array(["S"]*len(data))),
@@ -426,6 +477,13 @@ def set_fields_from_header(self, columns: List[str]):
             self.fields = [next(field for field in fields if field.name == col) for col in columns]
             assert np.array_equal(columns, [field.name for field in self.fields])
 
+        def get_field_by_number(self, field_nr: int, field_type: type=object):
+            if self.fields is None:
+                return super().get_field_by_number(field_nr, field_type)
+            col_id, t = next((i, field.type) for i, field in enumerate(dataclasses.fields(self.dataclass)) if field.name == self.fields[field_nr].name)
+            return super().get_field_by_number(col_id, t)
+            # fields = self.fields if self.fields is not None else dataclasses.fields(self.dataclass)
+
         def get_data(self) -> _dataclass:
             """Parse the data in the buffer according to the fields in _dataclass
 
@@ -447,6 +505,8 @@ def get_data(self) -> _dataclass:
                     col = self.get_integers(col_number).ravel()
                 elif field.type == float:
                     col = self.get_floats(col_number).ravel()
+                elif field.type == bool:
+                    col = self.get_integers(col_number).ravel().astype(bool)
                 elif field.type == List[int]:
                     col = self.get_split_ints(col_number)
                 elif field.type == List[bool]:
@@ -501,58 +561,11 @@ class Bed6Buffer(BedBuffer):
 class Bed12Buffer(Bed6Buffer):
     dataclass = Bed12
 
+
 class BdgBuffer(BedBuffer):
     dataclass = BedGraph
 
 
-class VCFBuffer(DelimitedBuffer):
-    dataclass = VCFEntry
-
-    def get_data(self):
-        data = super().get_data()
-        data.position -= 1
-        return data
-
-    @classmethod
-    def from_data(cls, data: BNPDataClass) -> "DelimitedBuffer":
-        data = dataclasses.replace(data, position=data.position+1)
-        return super().from_data(data)
-
-
-class VCFMatrixBuffer(VCFBuffer):
-    dataclass = VCFEntry
-    genotype_dataclass = VCFGenotypeEntry
-    genotype_encoding = GenotypeRowEncoding
-
-    @classmethod
-    def __read_header(cls, file_object):
-        # NOT IMPLEMENTED version of read_header()
-        prev_line = None
-        for line in file_object:
-            line = line.decode()
-            if line[0] != cls.COMMENT:
-                file_object.seek(-len(line), 1)
-                break
-            prev_line = line
-        if prev_line is None:
-            return []
-
-        sample_names = prev_line.split("\t")[9:]
-        return sample_names
-
-    def get_data(self):
-        data = super().get_data()
-        genotypes = self.get_column_range_as_text(9, self._n_cols, keep_sep=True)
-        genotypes = EncodedArray(self.genotype_encoding.encode(genotypes), self.genotype_encoding)
-        return self.genotype_dataclass(*data.shallow_tuple(), genotypes)
-
-
-class PhasedVCFMatrixBuffer(VCFMatrixBuffer):
-    dataclass = VCFEntry
-    genotype_dataclass = PhasedVCFGenotypeEntry
-    genotype_encoding = PhasedGenotypeRowEncoding
-
-
 class NarrowPeakBuffer(DelimitedBuffer):
     dataclass = NarrowPeak
 
@@ -585,15 +598,14 @@ def _calculate_col_starts_and_ends(self, data, delimiters):
         return start_delimiters+1, end_delimiters
 
     def _col_ends(self, col):
-        return self.__col_ends[:, col]
+        return self._wig_col_ends[:, col]
 
     def _col_starts(self, col):
-        return self.__col_starts[:, col]
+        return self._wig_col_starts[:, col]
 
     def _validate(self):
         self._validated = True
 
-
     def __init__(self, data: EncodedArray, new_lines: np.ndarray, delimiters: np.ndarray = None, header_data=None):
         delimiters_mask = data == self.DELIMITER
         delimiters_mask[new_lines] = True
@@ -602,8 +614,11 @@ def __init__(self, data: EncodedArray, new_lines: np.ndarray, delimiters: np.nda
         starts, ends = self._calculate_col_starts_and_ends(data, delimiters)
         n_fields = next(i for i, d in enumerate(ends) if data[d] == '\n') + 1
         self._n_cols = n_fields
-        self.__col_starts = starts.reshape(-1, n_fields)
-        self.__col_ends = ends.reshape(-1, n_fields)
+        self._wig_col_starts = starts.reshape(-1, n_fields)
+        self._wig_col_ends = ends.reshape(-1, n_fields)
+
+    def _get_buffer_extractor(self) -> TextBufferExtractor:
+        return TextBufferExtractor(self._data, self._wig_col_starts, self._wig_col_ends)
 
     @classmethod
     def from_raw_buffer(cls, chunk: np.ndarray, header_data=None) -> "DelimitedBuffer":
diff --git a/bionumpy/io/dump_csv.py b/bionumpy/io/dump_csv.py
index ee5df1936..227317a12 100644
--- a/bionumpy/io/dump_csv.py
+++ b/bionumpy/io/dump_csv.py
@@ -15,22 +15,16 @@ def str_func(column):
         return column.encoding.decode(column)
     assert False
 
-def dump_csv(data_dict: List[Tuple], sep: str = "\t") -> EncodedArray:
-    """Put each field of the dataclass into a column in a buffer.
-
-    Parameters
-    data : bnpdataclass
-        Data
-    """
-
-    funcs = {int: ints_to_strings,
-             str: str_func,  #lambda x: x,
-             List[int]: int_lists_to_strings,
-             float: float_to_strings,
-             List[bool]: lambda x: int_lists_to_strings(x.astype(int), sep="")
-             }
 
+def get_column(values, field_type) -> EncodedRaggedArray:
     def get_func_for_datatype(datatype):
+        funcs = {int: ints_to_strings,
+                 str: str_func,  # lambda x: x,
+                 List[int]: int_lists_to_strings,
+                 float: float_to_strings,
+                 List[bool]: lambda x: int_lists_to_strings(x.astype(int), sep=""),
+                 bool: lambda x: ints_to_strings(x.astype(int))
+                 }
         if is_subclass_or_instance(datatype, Encoding):
             encoding = datatype
 
@@ -38,18 +32,35 @@ def dynamic(x):
                 if isinstance(x, EncodedRaggedArray):
                     return EncodedRaggedArray(EncodedArray(encoding.decode(x.ravel()), BaseEncoding), x.shape)
                 return EncodedArray(encoding.decode(x), BaseEncoding)
+
             return dynamic
         else:
             return funcs[datatype]
 
-    columns = [get_func_for_datatype(key)(value) for key, value in data_dict]
+    return get_func_for_datatype(field_type)(values)
+
+
+def dump_csv(data_dict: List[Tuple], sep: str = "\t") -> EncodedArray:
+    """Put each field of the dataclass into a column in a buffer.
+
+    Parameters
+    data : bnpdataclass
+        Data
+    """
+
+    columns = [get_column(value, key) for key, value in data_dict]
+    lines = join_columns(columns, sep)
+    return lines.ravel()
+
+
+def join_columns(columns: List[EncodedRaggedArray], sep: str) -> EncodedRaggedArray:
     lengths = np.concatenate([((column.lengths if
                                 isinstance(column, RaggedArray)
                                 else np.array([
-                                            column.shape[-1] if len(column.shape) == 2 else 1
-                                              ]*len(column))) + 1
+                                                  column.shape[-1] if len(column.shape) == 2 else 1
+                                              ] * len(column))) + 1
                                )[:, np.newaxis]
-                              
+
                               for column in columns], axis=-1).ravel()
     lines = EncodedRaggedArray(EncodedArray(np.empty(lengths.sum(), dtype=np.uint8), BaseEncoding),
                                lengths)
@@ -60,4 +71,4 @@ def dynamic(x):
         lines[i::n_columns, :-1] = column
     lines[:, -1] = sep
     lines[(n_columns - 1)::n_columns, -1] = "\n"
-    return lines.ravel()
+    return lines
diff --git a/bionumpy/io/file_buffers.py b/bionumpy/io/file_buffers.py
index 4c864a419..115924f56 100644
--- a/bionumpy/io/file_buffers.py
+++ b/bionumpy/io/file_buffers.py
@@ -1,11 +1,13 @@
 import sys
+from typing import List
+
 import numpy as np
 from io import FileIO
 from npstructures import RaggedArray, RaggedView, RaggedShape
 from .exceptions import FormatException
 from ..bnpdataclass import bnpdataclass
 from ..datatypes import SequenceEntry, SequenceEntryWithQuality
-from ..encoded_array import EncodedArray, EncodedRaggedArray
+from ..encoded_array import EncodedArray, EncodedRaggedArray, as_encoded_array
 from ..encodings import QualityEncoding, BaseEncoding
 
 NEWLINE = "\n"
@@ -33,12 +35,28 @@ def __init__(self, data: EncodedArray, new_lines: np.ndarray):
         self._data = data
         self._new_lines = np.asanyarray(new_lines)
         self._is_validated = False
-        self.size = self._data.size
+
+    @property
+    def size(self):
+        return self.data.size
+
+    @property
+    def data(self):
+        return self._data
+
+    def __getitem__(self, idx):
+        return NotImplemented
 
     def raise_if(condition, *args, **kwargs):
         if condition:
             raise FormatException(*args, **kwargs)
 
+    @property
+    def header_data(self):
+        if hasattr(self, "_header_data"):
+            return self._header_data
+        return None
+
     @property
     def n_lines(self):
         return len(self._new_lines)
@@ -160,15 +178,81 @@ def _move_2d_array_to_intervals(self, array, starts, ends):
 
     @classmethod
     def contains_complete_entry(cls, chunks):
-        n_new_lines = sum(np.sum(EncodedArray(chunk, BaseEncoding) == NEWLINE) for chunk in chunks)
+        n_new_lines = sum(np.count_nonzero(EncodedArray(chunk, BaseEncoding) == NEWLINE) for chunk in chunks)
         return n_new_lines >= cls.n_lines_per_entry
 
+    #def get_field_by_number(self, field_nr: int, field_type: type=object):
+    #    raise NotImplementedError
+
+
+class IncompleteEntryException(Exception):
+    pass
+
+
+class TextBufferExtractor:
+    def __init__(self, data: EncodedArray, field_starts: np.ndarray, field_ends: np.ndarray):
+        '''
+        field_starts: n_entries x n_fields
+        field_ends: n_entries x n_fields
+        '''
+        self._data = data
+        self._field_starts = field_starts
+        self._field_ends = field_ends
+        self._field_lens = field_ends-field_starts
+        self._n_fields = field_starts.shape[1]
+
+    @property
+    def data(self):
+        return self._data
+
+    def __len__(self):
+        return len(self._field_starts)
+
+    def __getitem__(self, idx):
+        return self.__class__(self._data, field_starts=self._field_starts[idx], field_ends=self._field_ends[idx])
+
+    def get_field_by_number(self, field_nr: int):
+        e = EncodedRaggedArray(self._data, RaggedView(self._field_starts.ravel(), self._field_lens.ravel()))
+        return e[field_nr::self._n_fields]
+
 
 class OneLineBuffer(FileBuffer):
     """ Base class for file formats where data fields are contained in lines."""
 
     n_lines_per_entry = 2
     _buffer_divisor = 32
+    _line_offsets = (1, 0)
+    _empty_lines = []
+
+    def __init__(self, buffer_extractor: TextBufferExtractor):
+        # super().__init__(*args, **kwargs)
+        self._is_validated=False
+        self._buffer_extractor = buffer_extractor
+        self._is_validated = True
+
+    @property
+    def n_lines(self):
+        return len(self._buffer_extractor)*self.n_lines_per_entry
+
+    @property
+    def data(self):
+        return self._buffer_extractor.data
+
+    @classmethod
+    def _get_buffer_extractor(cls, data, new_lines):
+        tmp = np.insert(new_lines, 0, -1)
+        ends = tmp[1:].reshape(-1, cls.n_lines_per_entry)
+        starts = tmp[:-1].reshape(-1, cls.n_lines_per_entry)+(np.array(cls._line_offsets)+1)
+        return TextBufferExtractor(data, starts, ends)
+
+    @classmethod
+    def contains_complete_entry(cls, chunks):
+        if len(chunks) == 1:
+            try:
+                return True, cls.from_raw_buffer(chunks[0])
+            except IncompleteEntryException:
+                return False
+        return super().contains_complete_entry(chunks)
 
     @classmethod
     def from_raw_buffer(cls, chunk, header_data=None) -> "OneLineBuffer":
@@ -195,24 +279,50 @@ def from_raw_buffer(cls, chunk, header_data=None) -> "OneLineBuffer":
         chunk = EncodedArray(chunk, BaseEncoding)
         new_lines = np.flatnonzero(chunk == NEWLINE)
         n_lines = new_lines.size
-        assert n_lines >= cls.n_lines_per_entry, "No complete entry in buffer. Try increasing chunk_size."
+        if n_lines < cls.n_lines_per_entry:
+            raise IncompleteEntryException("No complete entry in buffer. Try increasing chunk_size.")
         new_lines = new_lines[: n_lines - (n_lines % cls.n_lines_per_entry)]
         chunk = chunk[: new_lines[-1] + 1]
-        return cls(chunk[: new_lines[-1] + 1], new_lines)
+        data = chunk[: new_lines[-1] + 1]
+        cls._validate(data, new_lines)
+        return cls(cls._get_buffer_extractor(data, new_lines))
+
+    @property
+    def lines(self):
+        if not hasattr(self, "__lines"):
+            starts = np.insert(self._new_lines, 0, -1)
+            lengths = np.diff(starts)
+            self.__lines = EncodedRaggedArray(self._data, RaggedShape(lengths))
+        return self.__lines
+
+    @property
+    def entries(self):
+        if not hasattr(self, "__entries"):
+            lengths = np.diff(self._new_lines[self.n_lines_per_entry-1::self.n_lines_per_entry])
+            lengths = np.insert(lengths, 0, self._new_lines[self.n_lines_per_entry-1]+1)
+            self.__entries = EncodedRaggedArray(self._data, RaggedShape(lengths))
+        return self.__entries
 
     def get_data(self) -> bnpdataclass:
         """Get and parse fields from each line"""
-        self.validate_if_not()
-        starts = np.insert(self._new_lines, 0, -1)
-        lengths = np.diff(starts)
-        self.lines = EncodedRaggedArray(self._data, RaggedShape(lengths))
-        sequences = self.lines[1::self.n_lines_per_entry, :-1]
-        headers = self.lines[:: self.n_lines_per_entry, 1:-1]
+        headers, sequences = [self._buffer_extractor.get_field_by_number(i) for i in (0, 1)]
         return SequenceEntry(headers, sequences)
 
+    def get_field_by_number(self, i: int, t: type=object):
+        """ Get a field indexed by number"""
+        return self._buffer_extractor.get_field_by_number(i)
+
+    def __getitem__(self, idx):
+        return self.__class__(self._buffer_extractor[idx])
+
     def count_entries(self) -> int:
         """Count number of entries in file"""
-        return len(self._new_lines)//self.n_lines_per_entry
+        return len(self._buffer_extractor)
+
+    def get_field_range_as_text(self, start, end):
+        """Get a range of fields as text"""
+        assert end == start+1
+        return self.get_text_field_by_number(start)
 
     @classmethod
     def from_data(cls, entries: bnpdataclass) -> "OneLineBuffer":
@@ -232,37 +342,64 @@ def from_data(cls, entries: bnpdataclass) -> "OneLineBuffer":
             A ASCII encoded buffer
         """
 
-        name_lengths = entries.name.lengths
-        sequence_lengths = entries.sequence.lengths
+        names = entries.name
+        sequences = entries.sequence
+        name_lengths = names.lengths
+        sequence_lengths = sequences.lengths
         line_lengths = np.hstack(
             (name_lengths[:, None] + 2, sequence_lengths[:, None] + 1)
         ).ravel()
         buf = EncodedArray(np.empty(line_lengths.sum(), dtype=np.uint8), BaseEncoding)
         lines = EncodedRaggedArray(buf, line_lengths)
         step = cls.n_lines_per_entry
-        lines[0::step, 1:-1] = entries.name
+        lines[0::step, 1:-1] = names
         lines[1::step, :-1] = EncodedRaggedArray(
-            EncodedArray(entries.sequence.encoding.decode(entries.sequence.ravel()),  entries.sequence.encoding),entries.sequence.shape)
+            EncodedArray(sequences.encoding.decode(sequences.ravel()), sequences.encoding), sequences.shape)
 
         lines[0::step, 0] = ">"
 
         lines[:, -1] = "\n"
         return buf
 
-    def _validate(self):
-        n_lines = self._new_lines.size
-        assert n_lines % self.n_lines_per_entry == 0, "Wrong number of lines in buffer"
+    @classmethod
+    def join_fields(cls, fields: List[EncodedRaggedArray]):
+        field_lengths = np.hstack([field.shape[1][:, None] for field in fields])
+        line_lengths = field_lengths+1
+        for i in range(len(fields)):
+            line_lengths[:, i] += cls._line_offsets[i]
+        entry_lengths = line_lengths.sum(axis=-1)
+        buffer_size = entry_lengths.sum()
+        buf = EncodedArray(np.empty(buffer_size, dtype=np.uint8), BaseEncoding)
+        lines = EncodedRaggedArray(buf, line_lengths.ravel())
+        step = cls.n_lines_per_entry
+        for i, field in enumerate(fields):
+            lines[i::step, cls._line_offsets[i]:-1] = field
+        lines[0::step, 0] = cls.HEADER
+        lines[:, -1] = "\n"
+        return buf
+
+    @classmethod
+    def _validate(cls, data, new_lines):
+        header = cls.HEADER
+        if data.size == 0 and new_lines.size==0:
+            return
+        n_lines = new_lines.size
+        n_lines_per_entry = cls.n_lines_per_entry
+        assert n_lines % n_lines_per_entry == 0, "Wrong number of lines in buffer"
         header_idxs = (
-            self._new_lines[self.n_lines_per_entry - 1 : -1 : self.n_lines_per_entry]
-            + 1
+                new_lines[n_lines_per_entry - 1: -1: n_lines_per_entry]
+                + 1
         )
-        if np.any(self._data[header_idxs] != self.HEADER) or self._data[0] != self.HEADER:
-            if self._data[0] != self.HEADER:
+
+        if np.any(data[header_idxs] != header) or data[0] != header:
+            if data[0] != header:
                 line_number = 0
             else:
-                line_number = (np.flatnonzero(self._data[header_idxs] != self.HEADER)[0]+1)*self.n_lines_per_entry
-            raise FormatException(f"Expected header line to start with {self.HEADER}" % self._data, line_number=line_number)
-        self._is_validated = True
+                line_number = (np.flatnonzero(data[header_idxs] != header)[0] + 1) * n_lines_per_entry
+            raise FormatException(f"Expected header line to start with {header}" % data, line_number=line_number)
+
+    def get_text_field_by_number(self, i):
+        return self.get_field_by_number(i)
 
 
 class TwoLineFastaBuffer(OneLineBuffer):
@@ -275,14 +412,25 @@ class FastQBuffer(OneLineBuffer):
     HEADER = "@"
     n_lines_per_entry = 4
     dataclass = SequenceEntryWithQuality
+    _line_offsets = (1, 0, 0, 0)
+    _empty_lines = [2]
+
+    def get_text_field_by_number(self, i: int) -> EncodedRaggedArray:
+        if i == 2:
+            return self._buffer_extractor.get_field_by_number(3)
+        return super().get_text_field_by_number(i)
+
+    def get_field_by_number(self, i: int, t: type=object):
+        if i == 2:
+            return QualityEncoding.encode(self.get_text_field_by_number(i))
+        else:
+            return super().get_field_by_number(i, t)
 
     def get_data(self):
         seq_entry = super().get_data()
-        quality = self.lines[3 :: self.n_lines_per_entry, :-1]
+        quality = self.get_field_by_number(2, QualityEncoding)
         return SequenceEntryWithQuality(
-            seq_entry.name, seq_entry.sequence,
-            #quality)
-            QualityEncoding.encode(quality))
+            seq_entry.name, seq_entry.sequence,quality)
 
     @classmethod
     def _get_line_lens(cls, entries):
@@ -300,13 +448,20 @@ def _get_line_lens(cls, entries):
             + 1
         )
 
-    def _validate(self):
-        super()._validate()
-        if np.any(self._data[self._new_lines[1::self.n_lines_per_entry] + 1] != "+"):
-            entry_number = np.flatnonzero(self._data[self._new_lines[1::self.n_lines_per_entry] + 1] != "+")[0]
-            line_number = 2+entry_number*self.n_lines_per_entry
-            raise FormatException(f"Expected '+' at third line of entry in {self._data}", line_number=line_number)
-        
+    @classmethod
+    def _validate(cls, data, new_lines):
+        super()._validate(data, new_lines)
+        n_lines_per_entry = cls.n_lines_per_entry
+        if np.any(data[new_lines[1::n_lines_per_entry] + 1] != "+"):
+            entry_number = np.flatnonzero(data[new_lines[1::n_lines_per_entry] + 1] != "+")[0]
+            line_number = 2 + entry_number * n_lines_per_entry
+            raise FormatException(f"Expected '+' at third line of entry in {data}", line_number=line_number)
+
+    @classmethod
+    def join_fields(cls, fields: List[EncodedRaggedArray]):
+        plus_line = as_encoded_array(['+']*len(fields[0]))
+        return super().join_fields(fields[:2]+[plus_line]+fields[2:])
+
     @classmethod
     def from_data(cls, entries):
         line_lengths = cls._get_line_lens(entries)
diff --git a/bionumpy/io/files.py b/bionumpy/io/files.py
index 1a9f96ba1..579ad9f8c 100644
--- a/bionumpy/io/files.py
+++ b/bionumpy/io/files.py
@@ -1,12 +1,14 @@
 from pathlib import PurePath
-import gzip
+from .gzip_reading import gzip
+# import gzip
 import dataclasses
 from .file_buffers import FastQBuffer, FileBuffer
 from .multiline_buffer import MultiLineFastaBuffer
 from .bam import BamBuffer, BamIntervalBuffer
-from .delimited_buffers import (VCFBuffer, BedBuffer, GfaSequenceBuffer,
+from .delimited_buffers import (BedBuffer, GfaSequenceBuffer,
                                 GFFBuffer, SAMBuffer, ChromosomeSizeBuffer,
                                 NarrowPeakBuffer, BdgBuffer, GTFBuffer)
+from .vcf_buffers import VCFBuffer
 from .wig import WigBuffer
 from .parser import NumpyFileReader, NpBufferedWriter
 from .exceptions import FormatException
diff --git a/bionumpy/io/gzip_reading.py b/bionumpy/io/gzip_reading.py
new file mode 100644
index 000000000..c70d1d4fe
--- /dev/null
+++ b/bionumpy/io/gzip_reading.py
@@ -0,0 +1,4 @@
+try:
+    from isal import igzip as gzip
+except ImportError:
+    import gzip
diff --git a/bionumpy/io/multiline_buffer.py b/bionumpy/io/multiline_buffer.py
index 9c6de232e..39e2f90fe 100644
--- a/bionumpy/io/multiline_buffer.py
+++ b/bionumpy/io/multiline_buffer.py
@@ -9,7 +9,8 @@
 
 
 class MultiLineBuffer(FileBuffer):
-    pass
+    SKIP_LAZY = True
+
 
 class MultiLineFastaBuffer(MultiLineBuffer):
     _new_entry_marker = ">"
diff --git a/bionumpy/io/npdataclassreader.py b/bionumpy/io/npdataclassreader.py
index 38f7ca7f0..ba3adf528 100644
--- a/bionumpy/io/npdataclassreader.py
+++ b/bionumpy/io/npdataclassreader.py
@@ -1,8 +1,12 @@
 from itertools import takewhile, repeat
+
+from .delimited_buffers import DelimitedBufferWithInernalComments
 from .parser import NumpyFileReader
 from ..bnpdataclass import BNPDataClass
 from .exceptions import FormatException
+from ..datatypes import GTFEntry, VCFGenotypeEntry
 from ..streams import NpDataclassStream
+from ..bnpdataclass.lazybnpdataclass import create_lazy_class, ItemGetter
 
 
 class NpDataclassReader:
@@ -15,6 +19,7 @@ class NpDataclassReader:
     """
     def __init__(self, numpyfilereader: NumpyFileReader):
         self._reader = numpyfilereader
+        self.__lazy_class = None
 
     def __enter__(self):
         return self
@@ -40,7 +45,26 @@ def read(self) -> BNPDataClass:
         4
 
         """
-        return self._reader.read().get_data()
+        chunk = self._reader.read()
+        should_be_lazy = self._should_be_lazy(chunk)
+        if should_be_lazy:
+            return self._get_lazy_class(chunk.dataclass)(ItemGetter(chunk, chunk.dataclass))
+        return chunk.get_data()
+
+    def _get_lazy_class(self, dataclass, header=None):
+        if self.__lazy_class is None:
+            self.__lazy_class = create_lazy_class(dataclass, header=header)
+        return self.__lazy_class
+
+    def _should_be_lazy(self, chunk):
+        should_be_lazy = False
+        if hasattr(chunk, 'get_field_by_number') and hasattr(chunk, 'dataclass'):
+            if not issubclass(chunk.dataclass, (GTFEntry)):
+                if not hasattr(chunk, 'genotype_dataclass') and not (
+                        hasattr(chunk, 'HAS_UNCOMMENTED_HEADER_LINE') and chunk.HAS_UNCOMMENTED_HEADER_LINE):
+                    if not isinstance(chunk, DelimitedBufferWithInernalComments):
+                        should_be_lazy = True
+        return should_be_lazy
 
     def read_chunk(self, min_chunk_size: int = 5000000, max_chunk_size: int = None) -> BNPDataClass:
         """Read a single chunk into memory
@@ -70,6 +94,9 @@ def read_chunk(self, min_chunk_size: int = 5000000, max_chunk_size: int = None)
         if chunk is None:
             return self._reader._buffer_type.dataclass.empty()
         try:
+            if self._should_be_lazy(chunk):
+                return self._get_lazy_class(chunk.dataclass, header=chunk.header_data)(ItemGetter(chunk, chunk.dataclass, n_lines_read))
+
             data = chunk.get_data()
         except FormatException as e:
             e.line_number += n_lines_read
diff --git a/bionumpy/io/parser.py b/bionumpy/io/parser.py
index 92632970e..b6bc283cf 100644
--- a/bionumpy/io/parser.py
+++ b/bionumpy/io/parser.py
@@ -2,6 +2,8 @@
 import numpy as np
 from typing.io import IO
 from npstructures import npdataclass
+
+from .exceptions import FormatException
 from ..streams import BnpStream
 from ..encoded_array import EncodedArray
 from ..streams.grouped import grouped_stream
@@ -108,22 +110,38 @@ def read_chunk(self, min_chunk_size: int = 5000000, max_chunk_size: int = None)
         complete_entry_found = False
         temp_chunks = []
         local_bytes_read = 0
-        local_lines_read = 0
         if len(self._prepend):
             temp_chunks.append(self._prepend)
+        made_buffer = None
         while not complete_entry_found:
             chunk = self._get_buffer(min_chunk_size, max_chunk_size)
             if chunk is None:
                 return None
             temp_chunks.append(chunk)
-            if max_chunk_size is not None and sum(chunk.size for chunk in chunks) > max_chunk_size:
+            if max_chunk_size is not None and sum(chunk.size for chunk in temp_chunks) > max_chunk_size:
                 raise Exception("No complete entry found")
             local_bytes_read += chunk.size
-            #local_lines_read += chunk.n_lines
-            complete_entry_found = self._buffer_type.contains_complete_entry(temp_chunks)
-            
-        chunk = np.concatenate(temp_chunks)
-        buff = self._buffer_type.from_raw_buffer(chunk, header_data=self._header_data)
+            try:
+                complete_entry_found = self._buffer_type.contains_complete_entry(temp_chunks)
+            except FormatException as e:
+                e.line_number += self.n_lines_read
+                raise e
+            if isinstance(complete_entry_found, tuple):
+                complete_entry_found, made_buffer = complete_entry_found
+
+        if made_buffer is not None:
+            buff = made_buffer
+        else:
+            if len(temp_chunks) == 1:
+                chunk = temp_chunks[0]
+            else:
+                chunk = np.concatenate(temp_chunks)
+            try:
+                buff = self._buffer_type.from_raw_buffer(chunk, header_data=self._header_data)
+            except FormatException as e:
+                e.line_number += self.n_lines_read
+                raise e
+
         self._prepend = []
         if not self._is_finished:
             if not self._do_prepend:
@@ -231,9 +249,19 @@ def write(self, data: npdataclass):
                 #getattr(self._buffer_type, 'HAS_UNCOMMENTED_HEADER_LINE', False):
             header_array = self._buffer_type.make_header(data)
             self._file_obj.write(header_array)
+        if len(data) == 0:
+            return
 
+        if hasattr(data, 'get_data_object'):
+            bytes_array = data.get_buffer(buffer_class=self._buffer_type)
+            # if not hasattr(self._buffer_type, 'get_column_range_as_text'):
+            #     data = data.get_data_object()
+            #     bytes_array = self._buffer_type.from_data(data)
+            # else:
+            #
+        else:
+            bytes_array = self._buffer_type.from_data(data)
 
-        bytes_array = self._buffer_type.from_data(data)
         if isinstance(bytes_array, EncodedArray):
             bytes_array = bytes_array.raw()
         self._file_obj.write(bytes(bytes_array))  # .tofile(self._file_obj)
diff --git a/bionumpy/io/vcf_buffers.py b/bionumpy/io/vcf_buffers.py
new file mode 100644
index 000000000..4aea5dca8
--- /dev/null
+++ b/bionumpy/io/vcf_buffers.py
@@ -0,0 +1,86 @@
+import dataclasses
+
+import numpy as np
+
+from ..encoded_array import EncodedArray
+from ..bnpdataclass import BNPDataClass
+from ..datatypes import VCFEntry, VCFGenotypeEntry, PhasedVCFGenotypeEntry, PhasedVCFHaplotypeEntry
+from ..encodings.vcf_encoding import GenotypeRowEncoding, PhasedGenotypeRowEncoding, PhasedHaplotypeRowEncoding
+from .delimited_buffers import DelimitedBuffer
+
+
+class VCFBuffer(DelimitedBuffer):
+    dataclass = VCFEntry
+
+    def get_data(self):
+        data = super().get_data()
+        data.position -= 1
+        return data
+
+    def get_field_by_number(self, field_nr: int, field_type: type=object):
+        val = super().get_field_by_number(field_nr, field_type)
+        if field_nr == 1:
+            val -= 1
+        return val
+
+    @classmethod
+    def from_data(cls, data: BNPDataClass) -> "DelimitedBuffer":
+        data = dataclasses.replace(data, position=data.position+1)
+        return super().from_data(data)
+
+
+class VCFMatrixBuffer(VCFBuffer):
+    dataclass = VCFEntry
+    genotype_dataclass = VCFGenotypeEntry
+    genotype_encoding = GenotypeRowEncoding
+
+    @classmethod
+    def __read_header(cls, file_object):
+        # NOT IMPLEMENTED version of read_header()
+        prev_line = None
+        for line in file_object:
+            line = line.decode()
+            if line[0] != cls.COMMENT:
+                file_object.seek(-len(line), 1)
+                break
+            prev_line = line
+        if prev_line is None:
+            return []
+
+        sample_names = prev_line.split("\t")[9:]
+        return sample_names
+
+    def get_data(self):
+        data = super().get_data()
+        genotypes = self.get_column_range_as_text(9, self._n_cols, keep_sep=True)
+        genotypes = EncodedArray(self.genotype_encoding.encode(genotypes), self.genotype_encoding)
+        return self.genotype_dataclass(*data.shallow_tuple(), genotypes)
+
+    def get_field_by_number(self, field_nr: int, field_type: type=object):
+        if field_nr != 9:
+            return super().get_field_by_number(field_nr, field_type)
+        else:
+            genotypes = self.get_column_range_as_text(9, self._n_cols, keep_sep=True)
+            genotypes = EncodedArray(self.genotype_encoding.encode(genotypes), self.genotype_encoding)
+            return genotypes
+
+
+#@bnpdataclass
+class Genotypes:
+    genotype: np.ndarray
+    phased: bool
+
+
+
+
+class PhasedVCFMatrixBuffer(VCFMatrixBuffer):
+    dataclass = VCFEntry
+    genotype_dataclass = PhasedVCFGenotypeEntry
+    genotype_encoding = PhasedGenotypeRowEncoding
+
+
+class PhasedHaplotypeVCFMatrixBuffer(VCFMatrixBuffer):
+    """Encodes genotype info using one column per haplotype (not genotype)"""
+    dataclass = VCFEntry
+    genotype_dataclass = PhasedVCFHaplotypeEntry
+    genotype_encoding = PhasedHaplotypeRowEncoding
diff --git a/bionumpy/io/vcf_header.py b/bionumpy/io/vcf_header.py
new file mode 100644
index 000000000..8add6429f
--- /dev/null
+++ b/bionumpy/io/vcf_header.py
@@ -0,0 +1,222 @@
+from collections import OrderedDict
+from dataclasses import dataclass
+from typing import Any, Callable, Dict, Set, Tuple, List, Mapping, Optional
+
+import re
+import warnings
+@dataclass
+class VCFHeader:
+    """VCFHeader"""
+    fileformat: str
+    source: str
+    fileDate: str
+    reference: str
+    FILTER: Mapping[str, Any]
+    FORMAT: Mapping[str, Any]
+    INFO: Mapping[str, Any]
+    contig: Mapping[str, Any]
+    optional: Mapping[str, List[Any]]
+
+def preprocess_number(x: str) -> Optional[int]:
+    """Preprocess the number in field such as Number. Return None if x
+    cannot be formated as integer.
+
+    Parameters
+    ----------
+    x: str
+        the value for the field. e.g. Number={$x}
+
+    Returns
+    -------
+    Optional[int]
+        the integer representation of x. None if x cannot be formated 
+        as integer.
+    """
+    regex = re.match('(\d)', x)
+    if regex:
+        return int(regex.group(1))
+    else:
+        return None
+
+def preprocess_type(x: str) -> type:
+    """Preprocess the type in field such as Type. Reuturn the class 
+    type.
+
+    Parameters
+    ----------
+    x: str
+        the value for the field. e.g. Type={$x}
+
+    Returns
+    -------
+    type
+        the corresponding class type.
+    """
+    mapping = {'Float': float, 'Integer': int, 'Flag': bool, 'String': str}
+    return mapping[x]
+
+
+def extract_identifier_and_content(line: str) -> Tuple[Optional[str]]:
+    """Separate the identifier and the content. Return (None, None) if
+    not in the correct VCF header format (##{$identifier}={$content}).
+
+    Parameters
+    ----------
+    line: str
+        the header line. e.g. ##{$identifier}={$content}
+
+    Returns
+    -------
+    Tuple[str]
+        the first element is the identifier, the second element is the 
+        content.  Return (None, None) if not in the correct VCF header 
+        format (##{$identifier}={$content}).
+
+    """
+    regex = re.search('^##(\S+?)=(.*)$', line)
+    if regex:
+        return regex.group(1), regex.group(2)
+    else:
+        message = "".join((
+            f"Header is not in the correct format. Expect ",
+            f'^##(\S+?)=(.*)$, get {line}. Ignore this line.'))
+        warnings.warn(message, RuntimeWarning)
+        return None, None
+
+def parse_mapping_header(
+        content: str,
+        preprocessors: Dict[str, Callable],
+        field_regex: str) -> Dict[str, Any]:
+    """Parse the mapping header.
+
+    Parameters
+    ----------
+    content : str
+        the content of the identifier. e.g. 
+        '##($identifier)=<$content>$'
+
+    preprocessors : Dict[str, Callable]
+        preprocessor for the value, the keys should be the field, the
+        values shoul be a function that takes the string content as 
+        inputs. e.g. {'Type':  preprocess_type}, where 
+        preprocess_type('Float')
+
+    identifier_regex : str
+        regular expression to extract the string content. e.g. 
+        '=(.+?)[,>]' would extract 'Float' from '..Type=Float,..'
+
+    Returns
+    -------
+    Dict[str, Any]
+        the mapping for the identifier.
+
+    """
+    results = dict()
+    for field, preprocess_fn in preprocessors.items():
+        regex = re.search(f'{field}{field_regex[field]}', content)
+        if regex:
+            results.setdefault(field, preprocess_fn(regex.group(1)))
+
+    return results
+
+class HeaderReflection:
+    """HeaderReflection that deals with the mapping of preprocessor, 
+    regular expression, string identifiers, mapping identifiers.
+
+    Attributes
+    ----------
+    string_identifiers: Set[str]
+        the set for identifiers that has a string content.
+
+    mapping_identifiers: Set[str]
+        the set for identifiers that has a mapping content.
+
+    preprocessor: Dict[str, Callable]
+        the preprocessor for different fields in the mapping headers.
+    
+    regex: Dict[str, str]
+        the regular expression used to extract the values of mapping 
+        headers.
+    """
+    string_identifiers: Set[str] = {'fileformat', 'fileDate', 'source', 'reference'}
+    mapping_identifiers: Set[str] = {'FILTER', 'FORMAT', 'INFO', 'contig'}
+
+    preprocessor: Dict[str, Callable] = {
+        'ID': lambda x: x,
+        'Number': preprocess_number,
+        'Type': preprocess_type,
+        'Description': lambda x: x
+    }
+    regex: Dict[str, str] = {
+        'ID': '=(.+?)[,>]',
+        'Number': '=(.+?)[,>]',
+        'Type': '=(.+?)[,>]',
+        'Description': '="(.+?)"'
+    }
+
+    @classmethod
+    def is_string_identifier(cls, field: str) -> bool:
+        return field in cls.string_identifiers
+
+    @classmethod
+    def is_mapping_identifier(cls, field: str) -> bool:
+        return field in cls.mapping_identifiers
+
+def parse_header(lines: str) -> VCFHeader:
+    """Parse the VCF header.
+
+    Parameters
+    ----------
+    lines: str
+        string of the VCF file.
+
+    Returns
+    -------
+    VCFHeader
+        The parsed VCFHeader
+    """
+    headers = dict()
+    for line in lines.split('\n'):
+        if not line.startswith('##'):
+            continue
+        identifier, content = extract_identifier_and_content(line)
+        if identifier is None:
+            continue
+
+        if HeaderReflection.is_string_identifier(identifier):
+            # for line '##FIELD=STRING
+            headers[identifier] = content
+        elif HeaderReflection.is_mapping_identifier(identifier):
+            # for line '##FIELD=<KEY1=VALUE1,KEY2=VALUE2...>'
+            if identifier not in headers:
+                headers[identifier] = OrderedDict()
+
+            mapping = parse_mapping_header(
+                    content,
+                    HeaderReflection.preprocessor,
+                    HeaderReflection.regex)
+
+            if mapping.get('ID', None):
+                headers[identifier][mapping.get('ID')] = mapping
+            else:
+                absent_id = 'Without ID'
+                if absent_id not in headers[identifier]:
+                    headers[identifier][absent_id] = list()
+                headers[identifier][absent_id].append(mapping)
+
+        else:
+            # for line ##FIELD=STRING where FIELD not in string identifier or mapping identifier
+            optional_identifier = 'optional'
+            if optional_identifier not in headers:
+                headers[optional_identifier] = dict()
+            if identifier not in headers[optional_identifier]:
+                headers[optional_identifier][identifier] = list()
+            headers[optional_identifier][identifier].append(content)
+    
+    optional_identifier = 'optional'
+    for identifier, contents in headers[optional_identifier].items():
+        if len(contents) > 1:
+            message = f"Duplicated values for {identifier} ({len(contents)} values)"
+            warnings.warn(message, RuntimeWarning)    
+
+    return VCFHeader(**headers)
\ No newline at end of file
diff --git a/bionumpy/plotting.py b/bionumpy/plotting.py
index effa53fe4..65084debd 100644
--- a/bionumpy/plotting.py
+++ b/bionumpy/plotting.py
@@ -1,28 +1,44 @@
 import numpy as np
 from npstructures import RunLengthArray
 from .genomic_data import GenomicArray, GenomicIntervals
-from .computation_graph import Node
+from .genomic_data.binned_genome import BinnedGenome
 from .io.matrix_dump import Matrix
 from .sequence.count_encoded import EncodedCounts
 from .encoded_array import EncodedRaggedArray
 import dataclasses
 
+
 @dataclasses.dataclass
 class Vector:
     data: np.ndarray
     names: list
 
+
 class Plotter:
-    def __init__(self, plt):
-        self.plt = plt
+    def __init__(self, plt=None):
+        self._plt = plt
         self._show = True
+        self._tried = False
+
+    @property
+    def plt(self):
+        if self._plt is None and not self._tried:
+            try:
+                import matplotlib.pyplot as _plt
+                plt.style.use('Solarize_Light2')
+                self._plt = _plt
+            except:
+                pass
+            self._tried = True
+        return self._plt
+
 
     def set_config(self, **kwargs):
         accepted_configs = {'show'}
         for key, value in kwargs.items():
             assert key in accepted_configs
             if key == 'show':
-                self._show =value
+                self._show = value
 
     def show(self, f=None):
         if not self._show:
@@ -55,7 +71,7 @@ def _plot_heatmap(self, matrix):
         data = np.asanyarray(matrix.data)
         n_rows, n_cols = data.shape
         im = ax.imshow(data)
-        
+
         # Show all ticks and label them with the respective list entries
 
         ax.set_xticks(np.arange(n_cols))
@@ -71,13 +87,13 @@ def _plot_heatmap(self, matrix):
                 ax.set_yticklabels(matrix.row_names.tolist())
             else:
                 ax.set_yticklabels(matrix.row_names)
-        
+
         # Rotate the tick labels and set their alignment.
         self.plt.setp(ax.get_xticklabels(), rotation=45, ha="right",
-                 rotation_mode="anchor")
-        
+                      rotation_mode="anchor")
+
         # Loop over data dimensions and create text annotations.
-        if n_rows*n_cols<100:
+        if n_rows * n_cols < 100:
             for i in range(n_rows):
                 for j in range(n_cols):
                     text = ax.text(j, i, data[i, j],
@@ -85,7 +101,6 @@ def _plot_heatmap(self, matrix):
         fig.tight_layout()
         self.show()
 
-
     def _plot_list(self, data_list):
         data_list = [self._conversion(data) for data in data_list]
         if isinstance(data_list[0], GenomicArray):
@@ -99,7 +114,7 @@ def _plot_list(self, data_list):
                 axes[i].title.set_text(chromosome)
             self.show(f)
         else:
-            raise NotImplementedError(f'{type(data)} not supported')        
+            raise NotImplementedError(f'{type(data)} not supported')
 
     def _plot_dict(self, data_dict):
         data_dict = {key: self._conversion(data) for key, data in data_dict.items()}
@@ -125,15 +140,15 @@ def _plot_dict(self, data_dict):
             f.legend()
             self.show()
         else:
-            raise NotImplementedError(f'{type(data)} not supported')        
-        pass 
+            raise NotImplementedError(f'{type(data)} not supported')
+        pass
 
     def _subplots_for_genome(self, gc):
         return self.plt.subplots(len(gc.chrom_sizes), 1, sharey=True, sharex=True)
-        
+
     def _plot_single(self, data, ax=None, label=None):
         data = self._conversion(data)
-        if isinstance(data, GenomicArray):
+        if isinstance(data, (GenomicArray, BinnedGenome)):
             f, axes = self._subplots_for_genome(data.genome_context)
             # f, axes = self.plt.subplots(1, len(data.genome_context.chrom_sizes), sharey=True, sharex=True)
             for i, chromosome in enumerate(data.genome_context.chrom_sizes.keys()):
@@ -142,7 +157,7 @@ def _plot_single(self, data, ax=None, label=None):
             self.show(f)
         elif isinstance(data, RunLengthArray):
             if ax is None:
-                f, ax = self.plt.subplots()# self._subplots_for_genome(data.genome_context)
+                f, ax = self.plt.subplots()  # self._subplots_for_genome(data.genome_context)
                 ax.plot(np.asarray(data), label=label)
                 self.show()
             else:
@@ -160,7 +175,6 @@ def _plot_single(self, data, ax=None, label=None):
 
     def __call__(self, data):
         return self.plot(data)
-    
 
     def plot(self, data):
         if isinstance(data, list):
@@ -183,26 +197,40 @@ def new_func(*args, **kwargs):
         args = [convert(a) for a in args]
         kwargs = {key: convert(val) for key, val in kwargs.items()}
         return func(*args, **kwargs)
+
     return new_func
 
 
 class PXWrapper:
-    def __init__(self, px):
+    def __init__(self, px=None):
         self._px = px
+        self._tried = False
+
+    @property
+    def px(self):
+        if self._px is None and not self._tried:
+            try:
+                import plotly.express as _px
+                self._px = _px
+            except:
+                pass
+            self._tried = True
+        return self._px
 
     def __getattr__(self, name):
-        func = getattr(self._px, name)
+        func = getattr(self.px, name)
         return convert_args_to_array(func)
 
-try:
-    import plotly.express as _px
-    px = PXWrapper(_px)
-except:
-    px = None
-    
-try:
-    import matplotlib.pyplot as plt
-    plt.style.use('Solarize_Light2')
-    plot = Plotter(plt)
-except:
-    plot = None
+# px = PXWrapper()
+# try:
+#     import plotly.express as _px
+#     px = PXWrapper(_px)
+# except:
+#     px = None
+
+# try:
+#     import matplotlib.pyplot as plt
+#     plt.style.use('Solarize_Light2')
+plot = Plotter()
+# except:
+#     plot = None
diff --git a/bionumpy/sequence/bloom_filter.py b/bionumpy/sequence/bloom_filter.py
index a843073ea..83eeda4da 100644
--- a/bionumpy/sequence/bloom_filter.py
+++ b/bionumpy/sequence/bloom_filter.py
@@ -1,27 +1,45 @@
+from functools import reduce
+
 import numpy as np
 
+'''
+P(X | M) = 
+'''
+'''
+8 bits per kmer
+cache optimized bloom filter
+
+'''
 
-def hash_function(offset, mod):
+def hash_function(offset):
     def f(kmer):
-        return (kmer ^ offset) % mod
+        return kmer ^ offset
+
     return f
 
 
 class BloomFilter:
-    def __init__(self, hash_functions, mask):
+    def __init__(self, mask_size, hash_functions):
         self._hash_functions = hash_functions
-        self._mask = mask
+        self._mask = np.zeros(mask_size, dtype=bool)
+
+    @classmethod
+    def from_m_and_k(cls, m, k, seed=12345):
+        offsets = np.random.RandomState(seed).randint(0, m, k)
+        return cls(m, [hash_function(offset, m) for offset in offsets])
 
     @classmethod
     def from_hash_functions_and_seqeuences(cls, hash_functions, sequence, mask_size):
-        mask = np.zeros(mask_size, dtype=bool)
-        for function in hash_functions:
-            mask[function(sequence)] = True
-        return cls(hash_functions, mask)
+        bloom_filter = cls(mask_size, hash_functions)
+        bloom_filter.insert(sequence)
+        return bloom_filter
+
+    def insert(self, sequences):
+        for function in self._hash_functions:
+            self._mask[function(sequences) % self._mask.size] = True
 
     def __getitem__(self, idx):
-        return np.all([self._mask[h(idx)] for h in self._hash_functions],
-                      axis=0)
+        return reduce(np.logical_and, (self._mask[h(idx) % self._mask.size] for h in self._hash_functions))
 
 
 class InterleavedBloomFilter:
@@ -41,4 +59,3 @@ def __getitem__(self, idx):
         kmer, seq_idx = idx
         return np.all([self._mask[h(kmer), seq_idx] for h in self._hash_functions],
                       axis=0)
-    
diff --git a/bionumpy/simulate/sequences.py b/bionumpy/simulate/sequences.py
index 368ae323b..64a5623c4 100644
--- a/bionumpy/simulate/sequences.py
+++ b/bionumpy/simulate/sequences.py
@@ -57,6 +57,9 @@ def simulate_reads_from_genome(genome: GenomicSequence, length: int = 150, n_rea
                                ignore_reads_with_n=False):
     """ Simulates reads on a genome. Yields chunks of SequenceEntryWithQuality objects
     """
+    if isinstance(rng, int):
+        rng = default_rng(rng)
+
     chromosomes = genome.genome_context.chrom_sizes
     genome_size = sum([size for size in chromosomes.values()])
 
diff --git a/bionumpy/streams/groupby_func.py b/bionumpy/streams/groupby_func.py
index 9efbe6981..048527db2 100644
--- a/bionumpy/streams/groupby_func.py
+++ b/bionumpy/streams/groupby_func.py
@@ -32,8 +32,12 @@ def get_ragged_changes(ragged_array):
 
 def join_groupbys(grouped_generator):
     double_grouped = itertools.groupby(itertools.chain.from_iterable(grouped_generator), lambda x: x[0])
-    return grouped_stream(((key, np.concatenate([g[1] for g in groups]))
-                           for key, groups in double_grouped),
+    def f(groups):
+        groups_ = [g[1] for g in groups]
+        return np.concatenate(groups_)
+
+    grouped_ = ((key, f(groups)) for key, groups in double_grouped)
+    return grouped_stream(grouped_,
                           grouped_generator.attribute_name if hasattr(grouped_generator, "attribute_name") else None)
 
 
diff --git a/bionumpy/util/ascii_hash.py b/bionumpy/util/ascii_hash.py
index b830834c7..a77317408 100644
--- a/bionumpy/util/ascii_hash.py
+++ b/bionumpy/util/ascii_hash.py
@@ -26,6 +26,8 @@ def _get_power_array(n, mod):
 
 
 def get_ascii_hash(encoded_array, mod):
+    if len(encoded_array) == 0:
+        return np.array([], dtype=int)
     powers = _get_power_array(np.max(encoded_array.shape[-1]), mod)
     if isinstance(encoded_array, EncodedRaggedArray):
         col_indices = column_index_array(encoded_array.shape)
diff --git a/docs_source/source/intervals.rst b/docs_source/source/intervals.rst
index 695fc07ea..25cba3002 100644
--- a/docs_source/source/intervals.rst
+++ b/docs_source/source/intervals.rst
@@ -13,10 +13,9 @@ Independent Operations
 Independent operations covers any thing we want to do with a set of intervals that can be done one interval at a time. This covers simple geometric operations such as shifting and resizing; filtering based on properties of intervals. Since both the start and stop (and strand) attribute are numpy arrays, we can implement alot of functionality simply by using numpy indexing and ufuncs:
 
     >>> import bionumpy as bnp
-    >>> from dataclasses import replace
     >>> intervals = bnp.open("example_data/small_interval.bed").read()
-    >>> extended_right = replace(intervals, stop=intervals.stop+10)
-    >>> shifted = replace(intervals, start=intervals.start+5, stop=intervals.stop+5)
+    >>> extended_right = bnp.replace(intervals, stop=intervals.stop+10)
+    >>> shifted = bnp.replace(intervals, start=intervals.start+5, stop=intervals.stop+5)
     >>> small = intervals[(intervals.stop-intervals.start)<50]
 
 Even though it is possible to write these kinds of functions using simple numpy functionality, it is often better to wrap them in functions with instructive names. In the `bnp.interval` module, a set of utility funcitons ffor working with intervals is provided.
diff --git a/scripts/gzip_debugging.py b/scripts/gzip_debugging.py
new file mode 100644
index 000000000..692120cb5
--- /dev/null
+++ b/scripts/gzip_debugging.py
@@ -0,0 +1,50 @@
+import gzip
+import time
+import bionumpy as bnp
+import mgzip
+import numpy as np
+
+from bionumpy import MultiLineFastaBuffer
+from bionumpy.io.npdataclassreader import NpDataclassReader
+from bionumpy.io.parser import NumpyFileReader
+
+
+def np_filereader(file):
+    with gzip.open(file, "rb") as f:
+        return np.frombuffer(f.read(), dtype=np.uint8)
+    
+    
+def np_filereader_read_chunks(file):
+    with gzip.open(file, "rb") as f:
+        return NumpyFileReader(f, buffer_type=MultiLineFastaBuffer).read_chunks()
+
+
+def read_python(file):
+    with gzip.open(file, "rb") as f:
+        return f.read()
+
+
+def read_bionumpy(file):
+    with bnp.open(file) as f:
+        return f.read()
+
+
+def read_mgzip(file):
+    with mgzip.open(file, "rt", thread=8) as f:
+        return f.read()
+
+
+funcs = [np_filereader_read_chunks, np_filereader, read_python, read_mgzip]
+
+
+def benchmark(func, file):
+    t0 = time.perf_counter()
+    func(file)
+    print("Time ", func.__name__, time.perf_counter() - t0)
+
+
+if __name__ == "__main__":
+    large_file = "../example_data/hg38.fa.gz"
+
+    for func in funcs:
+        benchmark(func, large_file)
\ No newline at end of file
diff --git a/scripts/jaccard_all_vs_all_example.py b/scripts/jaccard_all_vs_all_example.py
index 8156d82aa..8b50591c9 100644
--- a/scripts/jaccard_all_vs_all_example.py
+++ b/scripts/jaccard_all_vs_all_example.py
@@ -3,6 +3,8 @@
 Reads all bedfiles into memory.
 Files should already be sorted.
 """
+import itertools
+from typing import List
 
 import bionumpy as bnp
 from bionumpy.arithmetics.similarity_measures import jaccard
@@ -24,7 +26,7 @@ def jaccard_all_vs_all(chrom_sizes, bed_file_names):
 
     for file1, data1 in all_data.items():
         for file2, data2 in all_data.items():
-            print(file1, file2)
+            # print(file1, file2)
             if file1 == file2:
                 continue
 
@@ -38,13 +40,30 @@ def jaccard_all_vs_all(chrom_sizes, bed_file_names):
     return results
 
 
+def jaccard_func(mask_a, mask_b):
+    return (mask_a & mask_b).sum() / (mask_a | mask_b).sum()
+
+
+def jaccard(chrom_sizes_file: str, bed_files: List[str]):
+    genome = bnp.Genome.from_file(chrom_sizes_file)
+    masks = {filename: genome.read_intervals(filename).get_mask() for filename in bed_files}
+    results = {(file_a, file_b): jaccard_func(masks[file_a], masks[file_b]) for file_a, file_b in itertools.combinations(masks.keys(), r=2)}
+    return results
+
+
+def _test_profiling():
+    chrom, *inputs = '../example_data/hg38_unix_sorted.chrom.sizes ../benchmarks/results/intervals/ENCFF143HTO_mapped_reads_100k.bed ../benchmarks/results/intervals/ENCFF227NIG_mapped_reads_100k.bed'.split()
+    jaccard(chrom, inputs)
+
+
 if __name__ == "__main__":
     chrom_sizes = sys.argv[1]
     files = sys.argv[2:]
     # chrom_sizes = "example_data/hg38_unix_sorted.chrom.sizes"
     # files = "benchmarks/results/bed_files/ENCFF120PGJ.sorted.bed  benchmarks/results/bed_files/ENCFF405ZPR.sorted.bed  benchmarks/results/bed_files/ENCFF539MIO.sorted.bed  benchmarks/results/bed_files/ENCFF965PER.sorted.bed benchmarks/results/bed_files/ENCFF193LLN.sorted.bed  benchmarks/results/bed_files/ENCFF410SWS.sorted.bed".split()
-    result = jaccard_all_vs_all(chrom_sizes, files)
-    for r in result:
+    # result = jaccard_all_vs_all(chrom_sizes, files)
+    result = jaccard(chrom_sizes, files)
+    for r in result.items():
         print(r)
 
 # chrom_sizes = "example_data/hg38_unix_sorted.chrom.sizes"
diff --git a/scripts/mutation_signature_example.py b/scripts/mutation_signature_example.py
index 5b6910f29..725ebf3a7 100644
--- a/scripts/mutation_signature_example.py
+++ b/scripts/mutation_signature_example.py
@@ -2,7 +2,7 @@
 import numpy as np
 from bionumpy.streams.multistream import MultiStream
 from bionumpy.variants import count_mutation_types_genomic, count_mutation_types, MutationTypeEncoding
-from bionumpy.io.delimited_buffers import PhasedVCFMatrixBuffer
+from bionumpy.io import PhasedVCFMatrixBuffer
 from bionumpy.io.matrix_dump import matrix_to_csv, read_matrix
 from bionumpy import Genome
 import logging
diff --git a/scripts/plot_bins.py b/scripts/plot_bins.py
new file mode 100644
index 000000000..49d7cec14
--- /dev/null
+++ b/scripts/plot_bins.py
@@ -0,0 +1,6 @@
+import bionumpy as bnp
+import matplotlib
+
+genome = bnp.genomic_data.BinnedGenome.from_file('example_data/sacCer3.fa', bin_size=10000)
+genome.count_file('example_data/sacCer3_sample_variants.vcf.gz')
+bnp.plot(genome).show()
diff --git a/scripts/reverse_compliment_example.py b/scripts/reverse_compliment_example.py
new file mode 100644
index 000000000..0a64bd12d
--- /dev/null
+++ b/scripts/reverse_compliment_example.py
@@ -0,0 +1,15 @@
+import bionumpy as bnp
+input_filename = None
+output_filname = None
+
+
+def reverse_complement(input_filename, output_filename):
+    sequence_entries = bnp.open(input_filename, buffer_type=bnp.TwoLineFastaBuffer).read_chunks()
+    reversed_entries = bnp.sequence.get_reverse_complement(sequence_entries)
+    with bnp.open(output_filename, "w", buffer_type=bnp.TwoLineFastaBuffer) as outfile:
+        outfile.write(reversed_entries)
+
+def _test():
+    assert False
+    reverse_complement('benchmarks/results/dna_sequences/length150_nreads500000.fa', 'tmp.fa')
+
diff --git a/scripts/sequence_length_distribution_example.py b/scripts/sequence_length_distribution_example.py
index 2009eeda4..6cdadc04d 100644
--- a/scripts/sequence_length_distribution_example.py
+++ b/scripts/sequence_length_distribution_example.py
@@ -1,7 +1,11 @@
+import gzip
+
 import bionumpy as bnp
 import numpy as np
 import sys
 
+from bionumpy.io.parser import NumpyFileReader
+
 
 @bnp.streamable()
 def get_sequence_lengths(chunk):
@@ -20,6 +24,15 @@ def get_sorted_sequence_length_distribution(file_name):
 def test():
     get_sorted_sequence_length_distribution("example_data/big.fq.gz")
 
+def _test_2():
+    get_sorted_sequence_length_distribution("benchmarks/results/dna_sequences/ENCFF689IPX.fq.gz")
+
+def _test_3():
+    filename = "benchmarks/results/dna_sequences/ENCFF689IPX.fq.gz"
+    reader = NumpyFileReader(gzip.open(filename), buffer_type=bnp.FastQBuffer)
+    reader.set_prepend_mode()
+    for chunk in reader.read_chunks():
+        pass
 
 if __name__ == "__main__":
     file_name = sys.argv[1]
diff --git a/scripts/subsample_reads.py b/scripts/subsample_reads.py
new file mode 100644
index 000000000..17fae3823
--- /dev/null
+++ b/scripts/subsample_reads.py
@@ -0,0 +1,37 @@
+import numpy as np
+import bionumpy as bnp
+snakemake = None
+
+
+def subsample_reads(input_filename, output_filename):
+    n_entries = bnp.count_entries(input_filename, buffer_type=bnp.TwoLineFastaBuffer)
+    n_to_subsample = n_entries // 2
+    n_entries_left = n_entries
+    n_samples_left = n_to_subsample
+    out_file = bnp.open(output_filename, "w", buffer_type=bnp.TwoLineFastaBuffer)
+    for chunk in bnp.open(input_filename, buffer_type=bnp.TwoLineFastaBuffer).read_chunks():
+        if n_entries_left > len(chunk):
+            n_samples_in_this_chunk = np.random.hypergeometric(len(chunk), n_entries_left - len(chunk), n_samples_left)
+        else:
+            n_samples_in_this_chunk = n_samples_left
+        rows = np.random.choice(np.arange(len(chunk)), n_samples_in_this_chunk, replace=False)
+        n_samples_left -= n_samples_in_this_chunk
+        n_entries_left -= len(chunk)
+        out_file.write(chunk[rows])
+    assert n_samples_left == 0, n_samples_left
+    assert n_entries_left == 0, n_entries_left
+    out_file.close()
+
+def _test_count_entries():
+    assert False
+    filename = 'length150_nreads5000000.fa'
+    bnp.count_entries(f'benchmarks/results/dna_sequences/{filename}', buffer_type=bnp.TwoLineFastaBuffer)
+
+
+def _test():
+    assert False
+    # filename = 'ENCFF689IPX.fq.gz'
+    # filename = 'length150_nreads500000.fa'
+    filename = 'length150_nreads5000000.fa'
+    subsample_reads(f'benchmarks/results/dna_sequences/{filename}',
+                    'tmp.fa')
diff --git a/scripts/unique_intersect_example.py b/scripts/unique_intersect_example.py
new file mode 100644
index 000000000..91d00a063
--- /dev/null
+++ b/scripts/unique_intersect_example.py
@@ -0,0 +1,30 @@
+from functools import reduce
+
+import bionumpy as bnp
+
+
+def unique_intersect(filename_a: str, filename_b: str, chrom_sizes_file: str, out_filename: str):
+    genome = bnp.Genome.from_file(chrom_sizes_file, filter_function=None)
+    intervals = genome.read_intervals(filename_a, buffer_type=bnp.Bed6Buffer)
+    genome_mask = genome.read_intervals(filename_b).get_mask()
+    mask = genome_mask[intervals].any(axis=-1)
+    kept_intervals = intervals[mask]
+    with bnp.open(out_filename, 'w') as f:
+        f.write(kept_intervals.get_data())
+
+
+def _test_profiling():
+    name_a = 'ENCFF143HTO_mapped_reads_1m'
+    name_b = 'ENCFF491EEI'
+    unique_intersect(
+                     f'../benchmarks/results/intervals/{name_a}.bed',
+                     f'../benchmarks/results/bed_files/{name_b}.bed.gz',
+                     '../example_data/hg38.chrom.sizes',
+                     'tmp.bed')
+
+
+# test_profiling()
+
+if __name__ == '__main__':
+    import sys
+    unique_intersect(*sys.argv[1:])
diff --git a/scripts/vcf_allele_frequency_filtering_example.py b/scripts/vcf_allele_frequency_filtering_example.py
index 8bd29aaa4..c6d5bcefd 100644
--- a/scripts/vcf_allele_frequency_filtering_example.py
+++ b/scripts/vcf_allele_frequency_filtering_example.py
@@ -1,6 +1,8 @@
 import sys
 import numpy as np
 import bionumpy as bnp
+import bionumpy.io.vcf_buffers
+
 
 # Filter a vcf file by min allele count
 
@@ -14,11 +16,10 @@ def filter_on_allele_count(chunk, min_ac=10):
 
 
 def filter_file_on_allele_count(input_file, output_file, min_ac=10):
-    f = bnp.open(input_file, buffer_type=bnp.io.PhasedVCFMatrixBuffer)
-    output_file = bnp.open(output_file, "w", buffer_type=bnp.PhasedVCFMatrixBuffer)
+    f = bnp.open(input_file, buffer_type=bionumpy.io.vcf_buffers.PhasedVCFMatrixBuffer)
+    output_file = bnp.open(output_file, "w", buffer_type=bionumpy.io.vcf_buffers.PhasedVCFMatrixBuffer)
     chunks = f.read_chunks(min_chunk_size=200000000)  # [f.read()]
     for chunk in chunks:
-        print("chunk")
         filtered = filter_on_allele_count(chunk, min_ac=min_ac)
         output_file.write(filtered)
 
@@ -26,6 +27,10 @@ def filter_file_on_allele_count(input_file, output_file, min_ac=10):
 def test():
     filter_file_on_allele_count("example_data/variants_phased.vcf", "test.vcf.tmp", min_ac=1)
 
+def _test_profile():
+    filter_file_on_allele_count("../benchmarks/results/vcfs/big_phased.vcf.gz", "test.vcf.tmp", min_ac=10)
+
+
 
 if __name__ == "__main__":
     filter_file_on_allele_count(sys.argv[1], sys.argv[2], int(sys.argv[3]))
diff --git a/scripts/vcf_filtering_example.py b/scripts/vcf_filtering_example.py
index 8e9e0bae6..8543324fd 100644
--- a/scripts/vcf_filtering_example.py
+++ b/scripts/vcf_filtering_example.py
@@ -1,6 +1,8 @@
 import bionumpy as bnp
 import numpy as np
 
+import bionumpy.io.vcf_buffers
+
 
 def remove_variants_with_N(chunk):
     # find and remove variants that have N in ref or alt sequence
@@ -9,7 +11,7 @@ def remove_variants_with_N(chunk):
     new_variants = chunk[~mask]
 
     # write new variants to file:
-    with bnp.open("tmp_variants.vcf", "w", buffer_type=bnp.VCFMatrixBuffer) as f:
+    with bnp.open("tmp_variants.vcf", "w", buffer_type=bionumpy.io.vcf_buffers.VCFMatrixBuffer) as f:
         f.write(new_variants)
 
 
@@ -22,7 +24,7 @@ def test():
     analyses = [remove_variants_with_N, phsased_vcf_example]
     analyses = [remove_variants_with_N]
     for func in analyses:
-        f = bnp.open("example_data/variants_with_header.vcf", buffer_type=bnp.VCFMatrixBuffer)
+        f = bnp.open("example_data/variants_with_header.vcf", buffer_type=bionumpy.io.vcf_buffers.VCFMatrixBuffer)
         data = f.read()
         func(data)
 
diff --git a/scripts/vcf_reading.py b/scripts/vcf_reading.py
index 28672ab55..5170c9b88 100644
--- a/scripts/vcf_reading.py
+++ b/scripts/vcf_reading.py
@@ -2,10 +2,12 @@
 from bionumpy import IndexedFasta, bnp_open, chromosome_map
 import numpy as np
 
+
 @chromosome_map(reduction=sum)
 def func(variants, reference):
     return np.bincount(reference[variants.position], minlength=256)
 
+
 if __name__ == "__main__":
     vcf_filename = sys.argv[1]
     reference_filename = sys.argv[2]
@@ -13,11 +15,11 @@ def func(variants, reference):
     reference = bnp_open(reference_filename, remove_chr=True)
     print(func(vcf_file, reference))
 
-    #for chromosome, variants in vcf_file:
+    # for chromosome, variants in vcf_file:
     #     print(func(variants, reference[chromosome])) 
-        # 
-        # ref_seq = reference[chromosome]
-        # s = ref_seq[variants.position]
-        # print(chromosome, len(variants))
-        # print(s[:10])
-        # print(variants[:10])
+    #
+    # ref_seq = reference[chromosome]
+    # s = ref_seq[variants.position]
+    # print(chromosome, len(variants))
+    # print(s[:10])
+    # print(variants[:10])
diff --git a/setup.py b/setup.py
index ebc0aa861..871a10701 100644
--- a/setup.py
+++ b/setup.py
@@ -10,7 +10,7 @@
 with open('HISTORY.rst') as history_file:
     history = history_file.read()
 
-    requirements = ['numpy>=1.20,<1.24',
+    requirements = ['numpy>=1.20',
                     'npstructures>=0.2.9']
 # 'npstructures @ git+https://github.com/knutdrand/npstructures.git']
 
@@ -48,6 +48,7 @@
     url='https://github.com/bionumpy/bionumpy',
     version='0.2.21',
     zip_safe=False,
+    extras_require={'full': ['isal']}
 )
 
 # python -m build


Test patch:
diff --git a/tests/api_tests/immune_ml.py.pyi b/tests/api_tests/immune_ml.py.pyi
new file mode 100644
index 000000000..e69de29bb
diff --git a/tests/buffers.py b/tests/buffers.py
index 917bc5e23..7e7758fa7 100644
--- a/tests/buffers.py
+++ b/tests/buffers.py
@@ -1,7 +1,8 @@
 import pytest
 from bionumpy.io.file_buffers import FastQBuffer, TwoLineFastaBuffer
 from bionumpy.datatypes import SequenceEntry, SequenceEntryWithQuality, Interval, SNP, SAMEntry, VCFEntry, Bed12, Bed6, BedGraph
-from bionumpy.io.delimited_buffers import BedBuffer, VCFBuffer, GfaSequenceBuffer, Bed12Buffer, Bed6Buffer
+from bionumpy.io.delimited_buffers import BedBuffer, GfaSequenceBuffer, Bed12Buffer, Bed6Buffer
+from bionumpy.io import VCFBuffer
 from bionumpy.io.wig import WigBuffer
 from bionumpy.io.multiline_buffer import MultiLineFastaBuffer
 from bionumpy.encoded_array import EncodedArray
diff --git a/tests/genomic_fixtures.py b/tests/genomic_fixtures.py
index 0134070d1..0e2651a78 100644
--- a/tests/genomic_fixtures.py
+++ b/tests/genomic_fixtures.py
@@ -71,6 +71,12 @@ def windows(genome_context):
     return w
 
 
+@pytest.fixture
+def location_entries():
+    return LocationEntry(
+        ['chr1', 'chr2', 'chr2'],
+        [19, 20, 50])
+
 @pytest.fixture
 def genomic_locations(genome_context):
     return GenomicLocation.from_fields(
diff --git a/tests/test_binned_genome.py b/tests/test_binned_genome.py
new file mode 100644
index 000000000..54d2d741e
--- /dev/null
+++ b/tests/test_binned_genome.py
@@ -0,0 +1,14 @@
+import numpy as np
+
+from bionumpy.genomic_data.binned_genome import BinnedGenome
+from .genomic_fixtures import location_entries, genome_context, chrom_sizes
+
+
+def test_count(location_entries, genome_context):
+    bg = BinnedGenome(genome_context, bin_size=10)
+    bg.count(location_entries)
+    true = {
+        'chr1': np.array([0, 1, 0, 0, 0, 0, 0, 0, 0, 0]),
+        'chr2': np.array([0, 0, 1, 0, 0, 1, 0, 0, 0, 0])}
+    for key in set(true.keys()) | set(bg.count_dict.keys()):
+        np.testing.assert_equal(true[key], bg.count_dict[key])
diff --git a/tests/test_bloom_filter.py b/tests/test_bloom_filter.py
index adbaa0dfb..545e1b3df 100644
--- a/tests/test_bloom_filter.py
+++ b/tests/test_bloom_filter.py
@@ -18,15 +18,21 @@ def kmer_hashes_list():
 def hash_functions():
     return [hash_function(offset, 17) for offset in (13, 51)]
 
-
+@pytest.mark.skip
 def test_bloom_filter_precence(kmer_hashes, hash_functions):
     f = BloomFilter.from_hash_functions_and_seqeuences(hash_functions,
                                                        kmer_hashes, 17)
     assert np.all(f[kmer_hashes])
 
 
+@pytest.mark.skip
 def test_interleaved_bloom_filter_precence(kmer_hashes_list, hash_functions):
     f = InterleavedBloomFilter.from_hash_functions_and_seqeuences(hash_functions,
                                                                   kmer_hashes_list, 17)
     for i, seq in enumerate(kmer_hashes_list):
         assert np.all(f[seq, i])
+
+@pytest.mark.skip
+def test_needle_acceptance(sequences):
+    minimizers = get_minimizers(sequences, k=31, w=41)
+
diff --git a/tests/test_bnpdataclass.py b/tests/test_bnpdataclass.py
index 8279539b8..a442ea045 100644
--- a/tests/test_bnpdataclass.py
+++ b/tests/test_bnpdataclass.py
@@ -7,7 +7,7 @@
 from bionumpy.bnpdataclass.bnpdataclassfunction import bnpdataclassfunction
 from numpy.testing import assert_equal
 from bionumpy.util.testing import assert_bnpdataclass_equal
-import pandas as pd
+# import pandas as pd
 import bionumpy as bnp
 
 
@@ -68,7 +68,7 @@ def test_make_dataclass():
 def add(a, b):
     return a+b
 
-
+@pytest.mark.skip
 def test_from_pandas():
     persons = {'name': ['knut', 'marit'],
                'age': [35, 30]}
diff --git a/tests/test_encodedarrayfunctions.py b/tests/test_encodedarrayfunctions.py
index 5ff629583..4104b9723 100644
--- a/tests/test_encodedarrayfunctions.py
+++ b/tests/test_encodedarrayfunctions.py
@@ -2,6 +2,8 @@
 import dataclasses
 import bionumpy as bnp
 import numpy as np
+
+from bionumpy.bnpdataclass import replace
 from bionumpy.encodings.string_encodings import StringEncoding
 
 
@@ -36,7 +38,6 @@ def chrom_sizes():
 def test_interval_sort(intervals, chrom_sizes):
     truth = bnp.arithmetics.sort_intervals(intervals, sort_order=chrom_sizes.name.tolist())
     encoding = StringEncoding(chrom_sizes.name)
-    new_intervals = dataclasses.replace(intervals,
-                                        chromosome=encoding.encode(intervals.chromosome))
+    new_intervals = replace(intervals, chromosome=encoding.encode(intervals.chromosome))
     result = bnp.arithmetics.sort_intervals(new_intervals)
     np.testing.assert_array_equal(truth.start, result.start)
diff --git a/tests/test_encodings.py b/tests/test_encodings.py
index 3ddea99df..9900165cf 100644
--- a/tests/test_encodings.py
+++ b/tests/test_encodings.py
@@ -130,9 +130,10 @@ def test_encoding_sequence_entry():
 
     assert type(s.quality) == RaggedArray
     assert np.all(s.quality.ravel() == correct)
-
-    data = FastQBuffer.dataclass.stack_with_ragged(s)
-    assert np.all(data.quality == correct)
+    data = s
+    print(data.quality, data.quality.shape)
+    # data = FastQBuffer.dataclass.stack_with_ragged(s)
+    assert_raggedarray_equal(data.quality, [correct])
 
 
 @pytest.mark.parametrize("data", ["!!@-^", ["!!@@", "!+"]])
diff --git a/tests/test_forbes_acceptance.py b/tests/test_forbes_acceptance.py
new file mode 100644
index 000000000..2f4513f33
--- /dev/null
+++ b/tests/test_forbes_acceptance.py
@@ -0,0 +1,22 @@
+import pytest
+
+import bionumpy as bnp
+from bionumpy.arithmetics import forbes, sort_intervals
+
+# Read the bed files
+
+@pytest.mark.skip
+def test():
+    a = bnp.open("example_data/ctcf.bed.gz").read()
+    b = bnp.open("example_data/znf263.bed.gz").read()
+
+    # We also need chromosomes to do sorting
+    chromosome_sizes = bnp.open("example_data/hg38.chrom.sizes").read()
+    print(chromosome_sizes)
+
+    # sort the bed files
+    a_sorted = sort_intervals(a, sort_order=chromosome_sizes.name.tolist())
+    b_sorted = sort_intervals(b, sort_order=chromosome_sizes.name.tolist())
+
+    similarity = forbes(chromosome_sizes, a_sorted, b_sorted)
+    print(similarity)
diff --git a/tests/test_groupby.py b/tests/test_groupby.py
index 3419179d5..2d569e2fe 100644
--- a/tests/test_groupby.py
+++ b/tests/test_groupby.py
@@ -62,3 +62,12 @@ def test_join_groupbys(ragged_array, grouped, do_encode):
             true = EncodedRaggedArray(EncodedArray(true.ravel(), bnp.encodings.BaseEncoding), true.shape)
 
         assert_raggedarray_equal(g, true)
+
+
+#@pytest.mark.xfail
+def test_groupby_many_chunks():
+    file = "example_data/variants_with_header.vcf"
+    chunks = bnp.open(file).read_chunks(100)
+    for chromosome, variants in bnp.groupby(chunks, "chromosome"):
+        print(chromosome)
+        print(variants)
diff --git a/tests/test_io.py b/tests/test_io.py
index a8420103f..32e3017b1 100644
--- a/tests/test_io.py
+++ b/tests/test_io.py
@@ -5,27 +5,33 @@
 import numpy as np
 from io import BytesIO
 import bionumpy as bnp
+from bionumpy.bnpdataclass import bnpdataclass
 from bionumpy.io.files import NumpyFileReader, NpDataclassReader, NpBufferedWriter
-from .buffers import buffer_texts, combos, big_fastq_text, SequenceEntryWithQuality
+from .buffers import buffer_texts, combos, big_fastq_text, SequenceEntryWithQuality, VCFEntry
 from bionumpy.io.matrix_dump import matrix_to_csv, parse_matrix
-from bionumpy.util.testing import assert_bnpdataclass_equal, assert_encoded_array_equal, assert_encoded_raggedarray_equal
+from bionumpy.util.testing import assert_bnpdataclass_equal, assert_encoded_array_equal, \
+    assert_encoded_raggedarray_equal
 from numpy.testing import assert_equal
 
 
 @pytest.mark.parametrize("file_format", combos.keys())
-def test_read_write_roundtrip(file_format):
+@pytest.mark.parametrize("chunked", [False, True])
+def test_read_write_roundtrip(file_format, chunked):
     if file_format in ("multiline_fasta", "bed12", 'wig'):
         return
     _, _, buf_type = combos[file_format]
-    buffer_text = buffer_texts[file_format]*100
+    buffer_text = buffer_texts[file_format] * 100
     input_buffer = bytes(buffer_text, encoding="utf8")
     in_obj = BytesIO(input_buffer)
     out_buffer = bytes()
     out_obj = BytesIO(out_buffer)
     reader = NpDataclassReader(NumpyFileReader(in_obj, buffer_type=buf_type))
     writer = NpBufferedWriter(out_obj, buf_type)
-    for chunk in reader.read_chunks(200):
-        writer.write(chunk)
+    if chunked:
+        for chunk in reader.read_chunks(200):
+            writer.write(chunk)
+    else:
+        writer.write(reader.read())
     print(out_obj.getvalue())
     print(input_buffer)
     assert out_obj.getvalue() == input_buffer
@@ -49,13 +55,14 @@ def integers():
 @pytest.fixture
 def matrix_text(header, integers):
     sep = '\t'
-    return sep.join(header) + "\n" + "\n".join(sep.join(str(i) for i in row) for row in integers)+"\n"
+    return sep.join(header) + "\n" + "\n".join(sep.join(str(i) for i in row) for row in integers) + "\n"
 
 
 @pytest.fixture
 def matrix_text2(header, integers, row_names):
     sep = '\t'
-    return sep.join(['type'] + header) + "\n" + "\n".join(sep.join([row_name] + [str(i) for i in row]) for row_name, row in zip(row_names,integers))+"\n"
+    return sep.join(['type'] + header) + "\n" + "\n".join(
+        sep.join([row_name] + [str(i) for i in row]) for row_name, row in zip(row_names, integers)) + "\n"
 
 
 def test_matrix_to_csv(header, integers, row_names):
@@ -63,11 +70,12 @@ def test_matrix_to_csv(header, integers, row_names):
     integers = np.arange(20).reshape(4, 5)
     text = matrix_to_csv(integers, header=header)
 
-    assert text.to_string() == ",".join(header) + "\n" + "\n".join(",".join(str(i) for i in row) for row in integers)+"\n"
+    assert text.to_string() == ",".join(header) + "\n" + "\n".join(
+        ",".join(str(i) for i in row) for row in integers) + "\n"
 
 
 def test_read_matrix(header, integers, matrix_text):
-    matrix = parse_matrix(matrix_text, rowname_type=None,field_type=int)
+    matrix = parse_matrix(matrix_text, rowname_type=None, field_type=int)
     assert_encoded_raggedarray_equal(matrix.col_names, header)
     assert_equal(integers, matrix.data)
 
@@ -79,15 +87,14 @@ def test_read_matrix_with_row_names(header, integers, matrix_text2, row_names):
     assert_encoded_raggedarray_equal(matrix.row_names, row_names)
 
 
-
-@pytest.mark.parametrize("buffer_name", ["bed", "vcf2", "vcf", "fastq", "fasta", "gfa_sequence", "multiline_fasta", 'wig'])
+@pytest.mark.parametrize("buffer_name",
+                         ["bed", "vcf2", "vcf", "fastq", "fasta", "gfa_sequence", "multiline_fasta", 'wig'])
 def test_buffer_read(buffer_name):
     _, true_data, buf_type = combos[buffer_name]
     text = buffer_texts[buffer_name]
     io_obj = BytesIO(bytes(text, encoding="utf8"))
     data = NpDataclassReader(NumpyFileReader(io_obj, buf_type)).read()
     for line, true_line in zip(data, true_data):
-        # print("#####", line, true_line)
         assert_bnpdataclass_equal(line, true_line)
 
 
@@ -103,7 +110,7 @@ def test_buffer_read_chunks(buffer_name, min_chunk_size):
 
 
 def test_read_big_fastq():
-    io_obj = BytesIO(bytes(big_fastq_text*20, encoding="utf8"))
+    io_obj = BytesIO(bytes(big_fastq_text * 20, encoding="utf8"))
     for b in NpDataclassReader(NumpyFileReader(io_obj, combos["fastq"][2])).read_chunks(min_chunk_size=1000):
         print(b)
 
@@ -147,9 +154,35 @@ def test_append_to_file(buffer_name):
 
 
 def test_write_dna_fastq():
-    _, data, buf_type= combos["fastq"]
+    _, data, buf_type = combos["fastq"]
     entry = SequenceEntryWithQuality(["name"], ["ACGT"], ["!!!!"])
     entry.sequence = bnp.as_encoded_array(entry.sequence, bnp.DNAEncoding)
     result = buf_type.from_raw_buffer(buf_type.from_data(entry)).get_data()
     print(result)
     assert np.all(entry.sequence == result.sequence)
+
+
+def test_write_empty():
+    entry = VCFEntry([], [], [], [],
+                     [], [], [], [])
+    with bnp.open('tmp.vcf', 'w') as f:
+        f.write(entry)
+
+
+def test_read_write_bool():
+    @bnpdataclass
+    class BoolDC:
+        param: bool
+
+    obj = BoolDC(param=[True, False, True])
+    buf_type = bnp.io.delimited_buffers.get_bufferclass_for_datatype(BoolDC, delimiter='\t', has_header=True)
+
+    with bnp.open('tmp_bool.tsv', 'w', buf_type) as f:
+        f.write(obj)
+
+    with bnp.open('tmp_bool.tsv', 'r', buf_type) as f:
+        obj2 = f.read()
+
+    assert np.array_equal(obj.param, obj2.param)
+
+    os.remove('tmp_bool.tsv')
diff --git a/tests/test_io_exceptions.py b/tests/test_io_exceptions.py
index ab7a7f0bf..3b247714a 100644
--- a/tests/test_io_exceptions.py
+++ b/tests/test_io_exceptions.py
@@ -112,7 +112,7 @@ def test_npdataclass_raises_format_exception_bed(data):
     npfilereader = NumpyFileReader(fobj, buffer_type=BedBuffer)
     reader = NpDataclassReader(npfilereader)
     with pytest.raises(FormatException) as e:
-        for _ in reader.read_chunks(200):
-            pass
+        for chunk in reader.read_chunks(200):
+            chunk.stop
     assert e.value.line_number == 4*100+line_number
 
diff --git a/tests/test_lazybnpdataclass.py b/tests/test_lazybnpdataclass.py
new file mode 100644
index 000000000..8e057a2b3
--- /dev/null
+++ b/tests/test_lazybnpdataclass.py
@@ -0,0 +1,169 @@
+import dataclasses
+
+import numpy as np
+from numpy.testing import assert_array_equal
+
+import bionumpy as bnp
+import pytest
+
+from bionumpy.bnpdataclass import bnpdataclass
+from bionumpy.bnpdataclass.lazybnpdataclass import create_lazy_class, ItemGetter
+from bionumpy.util.testing import assert_bnpdataclass_equal, assert_encoded_raggedarray_equal
+from .buffers import combos
+
+
+@dataclasses.dataclass
+class DummyClass:
+    a: int
+    b: str
+
+
+@bnpdataclass
+class DummyBNP:
+    a: int
+    b: str
+
+
+@pytest.mark.skip
+def test_buffer_backed_descriptor():
+    buffer = [0, 1, 2]
+
+    class NewClass:
+        attr_a: BufferBackedDescriptor(buffer, 0, str)
+        attr_b: BufferBackedDescriptor(buffer, 1, int)
+        attr_c: BufferBackedDescriptor(buffer, 2, float)
+
+
+def itemgetter(dataclass):
+    str_getter = lambda i: f'str({i})'
+    int_getter = lambda i: i * 2
+    return {str: str_getter, int: int_getter}
+
+
+def silly_property_getter():
+    i = 0
+
+    def my_prop():
+        nonlocal i
+        i += 1
+        return i
+
+    return my_prop()
+
+
+def simple_backend(var_name):
+    if var_name == 'a':
+        return 10
+    elif var_name == 'b':
+        return 'hei'
+    return None
+
+
+class BufferMock:
+
+    def validate_if_not(self):
+        pass
+
+    def get_field_by_number(self, i: int, t: type):
+        if i == 0:
+            return 10
+        elif i == 1:
+            return 'hei'
+        return None
+
+
+class BetterBufferMock:
+    def __init__(self, *data):
+        self.data = data
+
+    def get_field_by_number(self, i, t):
+        return self.data[i]
+
+    def validate_if_not(self):
+        pass
+
+@pytest.mark.parametrize('backend', [simple_backend, ItemGetter(BufferMock(), DummyClass)])
+def test_create_lazy_class(backend):
+    lazy_obj = create_lazy_class(DummyClass)(backend)
+    assert lazy_obj.a == 10
+    assert lazy_obj.b == 'hei'
+    lazy_obj.b = 20
+    assert lazy_obj.b == 20
+
+
+@pytest.fixture
+def lazy_bnp():
+    b = bnp.as_encoded_array(['hei'])
+    a = np.array([10])
+    mock = BetterBufferMock(a, b)
+    return create_lazy_class(DummyBNP)(ItemGetter(mock, DummyBNP))
+
+
+def test_bnp_lazy(lazy_bnp):
+    assert_array_equal(lazy_bnp.a, [10])
+    assert_encoded_raggedarray_equal(lazy_bnp.b, ['hei'])
+    lazy_bnp.a = [20]
+    assert_array_equal(lazy_bnp.a, [20])
+
+
+def test_lazy_with_fasta_buffer():
+    text, data, bt = combos['fasta']
+    b = bt.from_raw_buffer(bnp.as_encoded_array(text))
+    lazy = create_lazy_class(bnp.datatypes.SequenceEntry)(ItemGetter(b, bnp.datatypes.SequenceEntry))
+    assert_encoded_raggedarray_equal(lazy.name, [data.name for data in data])
+    lazy.name = bnp.as_encoded_array(['deg'])
+    assert_encoded_raggedarray_equal(lazy.name, ['deg'])
+
+
+def get_pairwise_data(buffer_type):
+    text, data, bt = combos[buffer_type]
+    b = bt.from_raw_buffer(bnp.as_encoded_array(text))
+    cls = bt.dataclass
+    lazy = create_lazy_class(bt.dataclass)(ItemGetter(b, bt.dataclass))
+    tuples = [tuple(getattr(row, field.name) for field in dataclasses.fields(row)) for row in data]
+    return lazy, cls.from_entry_tuples(tuples)
+
+
+def test_lazy_with_bed_buffer():
+    text, data, bt = combos['bed']
+    b = bt.from_raw_buffer(bnp.as_encoded_array(text))
+    lazy = create_lazy_class(bnp.datatypes.Bed6)(ItemGetter(b, bnp.datatypes.Bed6))
+    assert_encoded_raggedarray_equal(lazy.chromosome, [data.chromosome for data in data])
+    lazy.chromosome = bnp.as_encoded_array(['deg'])
+    assert_encoded_raggedarray_equal(lazy.chromosome, ['deg'])
+
+
+@pytest.fixture
+def pair():
+    return get_pairwise_data('bed')
+
+
+@pytest.fixture
+def fastq_pair():
+    return get_pairwise_data('bed')
+
+
+@pytest.mark.parametrize('pair', [get_pairwise_data('fasta'), get_pairwise_data('bed')])
+@pytest.mark.parametrize('idx', [[0], [1], [0, 1]])
+def test_indexing(pair, idx):
+    lazy, dataclass = pair
+    assert_bnpdataclass_equal(lazy.get_data_object(), dataclass)
+    data_object = lazy[idx].get_data_object()
+    assert_bnpdataclass_equal(data_object, dataclass[idx])
+
+
+@pytest.mark.parametrize('idx', [[0], [1], [0, 1]])
+def test_mutated_indexing(idx):
+    lazy, dataclass = get_pairwise_data('bed')
+    new_starts = np.full(len(dataclass), 100000)
+    lazy.start = new_starts
+    dataclass.start = new_starts
+    assert_bnpdataclass_equal(lazy.get_data_object(), dataclass)
+    data_object = lazy[idx].get_data_object()
+    assert_bnpdataclass_equal(data_object, dataclass[idx])
+
+@pytest.mark.parametrize('pair', [get_pairwise_data('fasta'), get_pairwise_data('bed')])
+def test_str(pair):
+    lazy, dataclass = pair
+    assert str(lazy) == str(dataclass)
+    assert repr(lazy) == str(dataclass)
diff --git a/tests/test_parsers.py b/tests/test_parsers.py
index e13e71ff4..13b935704 100644
--- a/tests/test_parsers.py
+++ b/tests/test_parsers.py
@@ -4,7 +4,8 @@
 import numpy as np
 from bionumpy.io.file_buffers import FastQBuffer, TwoLineFastaBuffer
 from bionumpy.datatypes import Interval, SNP, SequenceEntry, VCFEntry
-from bionumpy.io.delimited_buffers import BedBuffer, VCFBuffer, GfaSequenceBuffer, get_bufferclass_for_datatype
+from bionumpy.io.delimited_buffers import BedBuffer, GfaSequenceBuffer, get_bufferclass_for_datatype
+from bionumpy.io import VCFBuffer
 from bionumpy.io.files import bnp_open
 from bionumpy.util.testing import assert_bnpdataclass_equal, assert_encoded_raggedarray_equal
 from .buffers import fastq_buffer, twoline_fasta_buffer, bed_buffer, vcf_buffer, vcf_buffer2, gfa_sequence_buffer, combos, data
@@ -183,7 +184,8 @@ def test_read_chunk_after_read_chunks_returns_empty_dataclass():
     file = bnp.open("example_data/reads.fq")
     chunks = list(file.read_chunks())
     new_chunk = file.read_chunk()
-    assert isinstance(new_chunk, type(chunks[0]))
+    assert isinstance(chunks[0],
+                      type(new_chunk))
 
 
 def test_read_gtf():
diff --git a/tests/test_vcf.py b/tests/test_vcf.py
index 5040e3e47..9c41742c8 100644
--- a/tests/test_vcf.py
+++ b/tests/test_vcf.py
@@ -1,4 +1,4 @@
-# from bionumpy.vcf import parse_header, parse_info_line
+from bionumpy.io.vcf_header import parse_header
 import pytest
 
 
@@ -49,7 +49,7 @@
 ##bcftools_concatVersion=1.4.1+htslib-1.4.1
 ##bcftools_concatCommand=concat -o all_variants_with_genotypes.vcf.gz -O z variants_with_genotypes_chr10_0.001.vcf.gz variants_with_genotypes_chr1_0.001.vcf.gz variants_with_genotypes_chr11_0.001.vcf.gz variants_with_genotypes_chr12_0.001.vcf.gz variants_with_genotypes_chr13_0.001.vcf.gz variants_with_genotypes_chr14_0.001.vcf.gz variants_with_genotypes_chr15_0.001.vcf.gz variants_with_genotypes_chr16_0.001.vcf.gz variants_with_genotypes_chr17_0.001.vcf.gz variants_with_genotypes_chr18_0.001.vcf.gz variants_with_genotypes_chr19_0.001.vcf.gz variants_with_genotypes_chr20_0.001.vcf.gz variants_with_genotypes_chr2_0.001.vcf.gz variants_with_genotypes_chr21_0.001.vcf.gz variants_with_genotypes_chr22_0.001.vcf.gz variants_with_genotypes_chr3_0.001.vcf.gz variants_with_genotypes_chr4_0.001.vcf.gz variants_with_genotypes_chr5_0.001.vcf.gz variants_with_genotypes_chr6_0.001.vcf.gz variants_with_genotypes_chr7_0.001.vcf.gz variants_with_genotypes_chr8_0.001.vcf.gz variants_with_genotypes_chr9_0.001.vcf.gz variants_with_genotypes_chrX_0.001.vcf.gz; Date=Fri Jan 22 16:54:29 2021
 ##bcftools_viewCommand=view -O z --regions 1:1-5000000 data/variants.vcf.gz; Date=Mon May 30 20:18:21 2022\
-""".split("\n")
+"""
 
 info_line = '''##INFO=<ID=AF,Number=A,Type=Float,Description="Estimated allele frequency in the range (0,1)">'''
 
@@ -58,7 +58,28 @@
 def test_parse_info():
     assert parse_info_line(info_line) == {"ID": "AF", "Number": "A", "Type": "Float"}
 
-@pytest.mark.skip("unimplemented")
+# @pytest.mark.skip("unimplemented")
 def test_parse_header():
-    assert parse_header(lines) == {"ID": "AF", "Number": "A", "Type": "Float"}
+    # INFO is not tested
+    header = parse_header(lines)
+    assert header.fileformat == 'VCFv4.3'
+    assert header.source == 'IGSRpipeline'
+    assert header.fileDate == '27022019_15h52m43s'
+    assert header.reference == 'ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.fa'
     
+    header_ans = dict()
+    header_ans['FILTER'] = {'PASS': {'ID': 'PASS', 'Description': 'All filters passed'}}
+    header_ans['FORMAT'] = {'GT': {'ID': 'GT', 'Number': 1, 'Type': str, 'Description': 'Phased Genotype'}}
+    header_ans['contig'] = {}
+    for c in [str(x) for x in range(1, 23)] + ['X']:
+        header_ans['contig'][c] = {'ID': c}
+
+    header_ans['optional'] = {}
+    header_ans['optional']['bcftools_viewVersion'] = ['1.4.1+htslib-1.4.1']
+    header_ans['optional']['bcftools_viewCommand'] = ['view -q 0.001 -O z ALL.chr10.shapeit2_integrated_snvindels_v2a_27022019.GRCh38.phased.vcf.gz; Date=Fri Jan 22 15:41:31 2021', 'view -O z --regions 1:1-5000000 data/variants.vcf.gz; Date=Mon May 30 20:18:21 2022']
+    header_ans['optional']['bcftools_concatVersion'] = ['1.4.1+htslib-1.4.1']
+    header_ans['optional']['bcftools_concatCommand'] = ['concat -o all_variants_with_genotypes.vcf.gz -O z variants_with_genotypes_chr10_0.001.vcf.gz variants_with_genotypes_chr1_0.001.vcf.gz variants_with_genotypes_chr11_0.001.vcf.gz variants_with_genotypes_chr12_0.001.vcf.gz variants_with_genotypes_chr13_0.001.vcf.gz variants_with_genotypes_chr14_0.001.vcf.gz variants_with_genotypes_chr15_0.001.vcf.gz variants_with_genotypes_chr16_0.001.vcf.gz variants_with_genotypes_chr17_0.001.vcf.gz variants_with_genotypes_chr18_0.001.vcf.gz variants_with_genotypes_chr19_0.001.vcf.gz variants_with_genotypes_chr20_0.001.vcf.gz variants_with_genotypes_chr2_0.001.vcf.gz variants_with_genotypes_chr21_0.001.vcf.gz variants_with_genotypes_chr22_0.001.vcf.gz variants_with_genotypes_chr3_0.001.vcf.gz variants_with_genotypes_chr4_0.001.vcf.gz variants_with_genotypes_chr5_0.001.vcf.gz variants_with_genotypes_chr6_0.001.vcf.gz variants_with_genotypes_chr7_0.001.vcf.gz variants_with_genotypes_chr8_0.001.vcf.gz variants_with_genotypes_chr9_0.001.vcf.gz variants_with_genotypes_chrX_0.001.vcf.gz; Date=Fri Jan 22 16:54:29 2021']
+    for field in ['FILTER', 'FORMAT', 'contig', 'optional']:
+      target = getattr(header, field)
+      for key in header_ans[field].keys():
+        assert header_ans[field][key] == target[key]
diff --git a/tests/test_vcf_encoding.py b/tests/test_vcf_encoding.py
index d9d4efc94..13a623808 100644
--- a/tests/test_vcf_encoding.py
+++ b/tests/test_vcf_encoding.py
@@ -2,12 +2,13 @@
 
 import bionumpy as bnp
 import bionumpy.encoded_array
-from bionumpy.encodings.vcf_encoding import PhasedGenotypeRowEncoding, GenotypeRowEncoding
+import bionumpy.io.vcf_buffers
+from bionumpy.encodings.vcf_encoding import PhasedGenotypeRowEncoding, GenotypeRowEncoding, PhasedHaplotypeRowEncoding
 
 
 def test_vcf_matrix_buffer():
     f = bnp.open("example_data/variants_with_header.vcf",
-                 buffer_type=bnp.io.delimited_buffers.PhasedVCFMatrixBuffer)
+                 buffer_type=bionumpy.io.vcf_buffers.PhasedVCFMatrixBuffer)
 
     out = bnp.open("test1.vcf", mode="w")
 
@@ -22,7 +23,7 @@ def test_vcf_matrix_buffer():
 
 def test_vcf_matrix_buffer_stream():
     f = bnp.open("example_data/variants_with_header.vcf",
-                 buffer_type=bnp.io.delimited_buffers.PhasedVCFMatrixBuffer)
+                 buffer_type=bionumpy.io.vcf_buffers.PhasedVCFMatrixBuffer)
 
     out = bnp.open("test1.vcf", mode="w")
     out.write(f.read_chunks())
@@ -44,6 +45,21 @@ def test_phased_genotype_encoding():
     assert np.all(encoded == correct)
 
 
+def test_phased_haplotype_encoding():
+    data = bionumpy.encoded_array.as_encoded_array([
+        "0|0\t0|1\t2|0\t",
+        "0|0\t4|1\t1|1\t"
+    ])
+    encoded = PhasedHaplotypeRowEncoding.encode(data)
+    correct = np.array(
+        [
+            [0, 0, 0, 1, 2, 0],
+            [0, 0, 4, 1, 1, 1]
+        ]
+    )
+    assert np.all(encoded == correct)
+
+
 def test_genotype_encoding():
     raw_data = [
         "0|0\t0/1\t./.\t",
@@ -62,11 +78,11 @@ def test_genotype_encoding():
 
 def test_parse_unphased_vcf():
     # example_data/variants.vcf has messy unphased and missing genotypes
-    f = bnp.open("example_data/variants.vcf", buffer_type=bnp.VCFMatrixBuffer)
+    f = bnp.open("example_data/variants.vcf", buffer_type=bionumpy.io.vcf_buffers.VCFMatrixBuffer)
     data = f.read()
     decoded = GenotypeRowEncoding.decode(data.genotypes)
     print(repr(decoded))
-    out_file = bnp.open("test.tmp", "w", buffer_type=bnp.VCFMatrixBuffer)
+    out_file = bnp.open("test.tmp", "w", buffer_type=bionumpy.io.vcf_buffers.VCFMatrixBuffer)
     out_file.write(data)
 
     written = open("test.tmp").read().split("\n")
@@ -75,7 +91,7 @@ def test_parse_unphased_vcf():
 
 
 def test_parse_phased_vcf():
-    f = bnp.open("example_data/variants_phased.vcf", buffer_type=bnp.PhasedVCFMatrixBuffer)
+    f = bnp.open("example_data/variants_phased.vcf", buffer_type=bionumpy.io.vcf_buffers.PhasedVCFMatrixBuffer)
     data = f.read()
     data = data.genotypes.raw()
     print(data)