diff --git a/bdpy/bdata/bdata.py b/bdpy/bdata/bdata.py
index 31e21adb..0dc1ed27 100644
--- a/bdpy/bdata/bdata.py
+++ b/bdpy/bdata/bdata.py
@@ -8,19 +8,66 @@
 
 
 import datetime
+import functools
 import inspect
 import os
 import re
 import time
 import warnings
-from typing import Optional, Tuple, Union
+from typing import (
+    TYPE_CHECKING,
+    Any,
+    Callable,
+    Dict,
+    Iterable,
+    KeysView,
+    List,
+    Optional,
+    Sequence,
+    Tuple,
+    TypeVar,
+    Union,
+    cast,
+    overload,
+)
 
 import h5py
 import numpy as np
 import scipy.io as sio
+from typing_extensions import Literal
 
 from .featureselector import FeatureSelector
-from .metadata import MetaData
+from .metadata import MetaData, MetaDataValue
+
+# Misc -----------------------------------------------------------------
+
+_F = TypeVar("_F", bound=Callable[..., Any])
+_T = TypeVar("_T")
+
+def _obsoleted_method(alternative: str) -> Callable[[_F], _F]:
+    """Return a decorator that warns about an obsolete method."""
+    def decorator(func: _F) -> _F:
+        @functools.wraps(func)
+        def wrapper(*args: Any, **kwargs: Any) -> Any:  # noqa: ANN401
+            funcname = func.__name__
+            warnings.warn(
+                f"'{funcname}' is obsoleted and kept for compatibility. Use '{alternative}' instead.",
+                UserWarning,
+                stacklevel=2
+            )
+            return func(*args, **kwargs)
+        return cast(_F, wrapper)
+    return decorator
+
+if TYPE_CHECKING:
+    # `numpy.typing` requires numpy >= 1.20 (`NDArray` requires >= 1.21),
+    # which is not guaranteed at runtime. Keep these imports type-check only.
+    from numpy.typing import NDArray
+
+    ApplyFuncIndex = Union[Sequence[int], NDArray[np.integer]]
+    ApplyFuncResult = Union[np.ndarray, Tuple[np.ndarray, ApplyFuncIndex]]
+
+SelectionOperand = Union[np.ndarray, float]
 
 # BData class ##########################################################
 
@@ -53,8 +100,8 @@ def __init__(self, file_name: Optional[str] = None, file_type: Optional[str] = N
         """Initialize BData instance."""
         self.__dataset: np.ndarray = np.ndarray((0, 0), dtype=float)
         self.__metadata = MetaData()
-        self.__header: dict = {}
-        self.__vmap: dict = {}
+        self.__header: Dict[str, Any] = {}
+        self.__vmap: Dict[str, Dict[float, str]] = {}
 
         if file_name is not None:
             self.load(file_name, file_type)
@@ -63,59 +110,76 @@ def __init__(self, file_name: Optional[str] = None, file_type: Optional[str] = N
 
     # dataset
     @property
-    def dataset(self):
+    def dataset(self) -> np.ndarray:
+        """Data matrix stored in the BData instance.
+
+        Rows correspond to samples, and columns correspond to data variables or
+        features. Column groups and attributes are described by `metadata`.
+        """
         return self.__dataset
 
     @dataset.setter
-    def dataset(self, value):
+    def dataset(self, value: np.ndarray) -> None:
         self.__dataset = value
 
     # metadata
     @property
-    def metadata(self):
+    def metadata(self) -> MetaData:
+        """Meta-data describing columns of `dataset`.
+
+        The metadata stores keys, descriptions, and per-column values. These
+        values are aligned with the columns of `dataset` and are used by
+        methods such as `select`, `get`, and `get_metadata`.
+        """
         return self.__metadata
 
     @metadata.setter
-    def metadata(self, value):
+    def metadata(self, value: MetaData) -> None:
         self.__metadata = value
 
     # header
     @property
-    def header(self):
+    def header(self) -> Dict[str, Any]:
+        """Header information associated with the BData instance.
+
+        The header stores auxiliary information such as creation time,
+        call stack, and values loaded from BData files.
+        Header keys are strings, while values are implementation-defined and
+        may include strings, numbers, lists, or values loaded from
+        external files.
+        """
         return self.__header
 
     # dataSet
     @property
-    def dataSet(self):
+    def dataSet(self) -> np.ndarray:  # noqa: N802
+        """Alias for `dataset` kept for backward compatibility."""
+        warnings.warn(
+            "'dataSet' is obsoleted and kept for compatibility. Use 'dataset' instead.",
+            UserWarning,
+            stacklevel=2
+        )
         return self.__dataset
 
     @dataSet.setter
-    def dataSet(self, value):
+    def dataSet(self, value: np.ndarray) -> None:  # noqa: N802
         self.__dataset = value
 
     # metaData
     @property
-    def metaData(self):
+    def metaData(self) -> MetaData:  # noqa: N802
+        """Alias for `metadata` kept for backward compatibility."""
+        warnings.warn(
+            "'metaData' is obsoleted and kept for compatibility. Use 'metadata' instead.",
+            UserWarning,
+            stacklevel=2
+        )
         return self.__metadata
 
     @metaData.setter
-    def metaData(self, value):
+    def metaData(self, value: MetaData) -> None:  # noqa: N802
         self.__metadata = value
 
-    # Misc -------------------------------------------------------------
-
-    def __obsoleted_method(alternative):
-        """Decorator for obsoleted functions."""
-        def __obsoleted_method_in(func):
-            import functools
-            @functools.wraps(func)
-            def wrapper(*args, **kwargs):
-                funcname = func.__name__
-                warnings.warn("'%s' is obsoleted and kept for compatibility. Use '%s' instead." % (funcname, alternative), UserWarning, stacklevel=2)
-                return func(*args, **kwargs)
-            return wrapper
-        return __obsoleted_method_in
-
 
     # Data modification ------------------------------------------------
 
@@ -153,7 +217,7 @@ def add(self, x: np.ndarray, name: str) -> None:
         self.metadata.set(name, column_value, column_description,
                           lambda x, y: np.hstack((y[:colnum_has], x[-colnum_add:])))
 
-    @__obsoleted_method('add')
+    @_obsoleted_method('add')
     def add_dataset(self, x: np.ndarray, attribute_key: str) -> None:
         """Add `x` to dataset with attribute meta-data key `attribute_key`.
 
@@ -184,17 +248,28 @@ def update(self, key: str, dat: np.ndarray) -> None:
         -------
         None
         """
-        mdind = [a == 1 for a in self.get_metadata(key)]
+        md = self.get_metadata(key)
+        if md is None:
+            raise ValueError(f"Meta-data key '{key}' not found.")
+
+        mdind = [a == 1 for a in md]
         self.dataset[:, np.array(mdind)] = dat
 
-    def add_metadata(self, key: str, value: np.ndarray, description: str = '', where: Optional[str] = None, attribute: Optional[str] = None) -> None:
+    def add_metadata(
+        self,
+        key: str,
+        value: MetaDataValue,
+        description: str = '',
+        where: Optional[str] = None,
+        attribute: Optional[str] = None,
+    ) -> None:
         """Add meta-data with `key`, `description`, and `value` to metadata.
 
         Parameters
         ----------
         key : str
             Meta-data key.
-        value : numpy.ndarray
+        value : array_like
             Meta-data array.
         description : str, optional
             Meta-data description.
@@ -219,6 +294,7 @@ def add_metadata(self, key: str, value: np.ndarray, description: str = '', where
             else:
                 where = attribute
 
+        add_value: MetaDataValue
         if where is not None:
             attr_ind = self.metadata.get(where, 'value') == 1
             add_value = np.array([np.nan for _ in range(self.metadata.get_value_len())])
@@ -228,7 +304,14 @@ def add_metadata(self, key: str, value: np.ndarray, description: str = '', where
 
         self.metadata.set(key, add_value, description)
 
-    def merge_metadata(self, key: str, sources, description: str = '', where: Optional[str] = None, method: str = 'logical_or') -> None:
+    def merge_metadata(
+        self,
+        key: str,
+        sources: Iterable[str],
+        description: str = '',
+        where: Optional[str] = None,
+        method: str = 'logical_or',
+    ) -> None:
         """Merage metadata rows."""
         if not method == 'logical_or':
             raise NotImplementedError('Only `logical_or` is implemented')
@@ -273,7 +356,7 @@ def set_metadatadescription(self, key: str, description: str) -> None:
         """
         self.metadata.set(key, None, description, lambda x, y: y)
 
-    @__obsoleted_method('set_metadatadescription')
+    @_obsoleted_method('set_metadatadescription')
     def edit_metadatadescription(self, metakey: str, description: str) -> None:
         """Set description of metadata specified by `key`.
 
@@ -290,11 +373,16 @@ def edit_metadatadescription(self, metakey: str, description: str) -> None:
         """
         self.set_metadatadescription(metakey, description)
 
-    def update_header(self, header) -> None:
+    def update_header(self, header: Dict[str, Any]) -> None:
         """Update header."""
         self.__header.update(header)
 
-    def applyfunc(self, func, where=None, **kargs):
+    def applyfunc(
+        self,
+        func: "Callable[..., ApplyFuncResult]",
+        where: Optional[Union[str, List[str]]] = None,
+        **kargs: Any,  # noqa: ANN401
+    ) -> "BData":
         """Apply `func` to the dataset."""
         if where is None:
             # FIXME
@@ -327,7 +415,7 @@ def applyfunc(self, func, where=None, **kargs):
                 #import pdb; pdb.set_trace()
 
                 ds[:, index] = fout[0]
-                ds[:, ~index] = self.dataset[np.ix_(ind_map, ~index)]
+                ds[:, ~index] = self.dataset[np.ix_(ind_map, ~index)]  # type: ignore[arg-type]
 
                 self.dataset = ds
             else:
@@ -338,7 +426,39 @@ def applyfunc(self, func, where=None, **kargs):
 
     # Data access ------------------------------------------------------
 
-    def select(self, condition: str, return_index: bool = False, verbose: bool = True) -> Union[np.ndarray, Tuple[np.ndarray, list]]:
+    @overload
+    def select(
+        self,
+        condition: str,
+        return_index: Literal[False] = False,
+        verbose: bool = True,
+    ) -> np.ndarray:
+        ...
+
+    @overload
+    def select(
+        self,
+        condition: str,
+        return_index: Literal[True],
+        verbose: bool = True,
+    ) -> Tuple[np.ndarray, np.ndarray]:
+        ...
+
+    @overload
+    def select(
+        self,
+        condition: str,
+        return_index: bool,
+        verbose: bool = True,
+    ) -> Union[np.ndarray, Tuple[np.ndarray, np.ndarray]]:
+        ...
+
+    def select(
+        self,
+        condition: str,
+        return_index: bool = False,
+        verbose: bool = True
+    ) -> Union[np.ndarray, Tuple[np.ndarray, np.ndarray]]:
         """Select data (columns) from dataset.
 
         Parameters
@@ -366,99 +486,105 @@ def select(self, condition: str, return_index: bool = False, verbose: bool = Tru
         - = (equal)
         - @ (conditional)
         """
-        expr_rpn = FeatureSelector(condition).rpn
+        rpn_tokens = FeatureSelector(condition).rpn
 
-        stack: list = []
-        buf_sel = []
+        stack: List[SelectionOperand] = []
+        buf_sel: List[int] = []
 
-        for i in expr_rpn:
-            if i == '=':
-                r = stack.pop()
-                l = stack.pop()
+        for token in rpn_tokens:
+            if token == '=':
+                right = stack.pop()
+                left = cast(np.ndarray, stack.pop())
 
-                stack.append(np.array([n == r for n in l], dtype=bool))
+                stack.append(np.array([n == right for n in left], dtype=bool))
 
-            elif i == 'top':
+            elif token == 'top':
                 # Dirty solution
 
                 # Need fix on handling 'None'
 
-                n = int(stack.pop()) # Num of elements to be selected
-                v = stack.pop()
+                num_selected = int(cast(float, stack.pop())) # Num of elements to be selected
+                values_to_rank = cast(np.ndarray, stack.pop())
 
-                order = self.__get_order(v)
+                order = self.__get_order(values_to_rank)
 
                 stack.append(order)
-                buf_sel.append(n)
+                buf_sel.append(num_selected)
 
-            elif i in ['|', '&', '-']:
-                r = stack.pop()
-                l = stack.pop()
+            elif token in ['|', '&', '-']:
+                right = cast(np.ndarray, stack.pop())
+                left = cast(np.ndarray, stack.pop())
 
-                if r.dtype != 'bool':
-                    # 'r' should be an order vector
+                if right.dtype != 'bool':
+                    # 'right' should be an order vector
                     num_sel = buf_sel.pop()
-                    r = self.__get_top_elm_from_order(r, num_sel)
+                    right = self.__get_top_elm_from_order(right, num_sel)
                     #r = np.array([ n < num_sel for n in r ], dtype = bool)
 
-                if l.dtype != 'bool':
-                    # 'l' should be an order vector
+                if left.dtype != 'bool':
+                    # 'left' should be an order vector
                     num_sel = buf_sel.pop()
-                    l = self.__get_top_elm_from_order(l, num_sel)
-                    #l = np.array([ n < num_sel for n in l ], dtype = bool)
+                    left = self.__get_top_elm_from_order(left, num_sel)
+                    #left = np.array([ n < num_sel for n in left ], dtype = bool)
 
-                if i == '|':
-                    result = np.logical_or(l, r)
-                elif i == '&':
-                    result = np.logical_and(l, r)
-                elif i == '-':
-                    result = np.logical_and(l, np.logical_not(r))
+                if token == '|':
+                    result = np.logical_or(left, right)
+                elif token == '&':
+                    result = np.logical_and(left, right)
+                elif token == '-':
+                    result = np.logical_and(left, np.logical_not(right))
 
                 stack.append(result)
 
-            elif i == '@':
+            elif token == '@':
                 # FIXME
                 # In the current version, the right term of '@' is assumed to
                 # be a boolean, and the left is to be an order vector.
 
-                r = stack.pop() # Boolean
-                l = stack.pop() # Float
+                right = cast(np.ndarray, stack.pop()) # Boolean
+                left = cast(np.ndarray, stack.pop()) # Float
 
-                l[~r] = np.inf
+                left[~right] = np.inf
 
-                selind = self.__get_top_elm_from_order(l, buf_sel.pop())
+                conditional_mask = self.__get_top_elm_from_order(left, buf_sel.pop())
 
-                stack.append(np.array(selind))
+                stack.append(np.array(conditional_mask))
 
             else:
-                if isinstance(i, str):
-                    if i.isdigit():
-                        # 'i' should be a criteria value
-                        i = float(i)
-                    else:
-                        # 'i' should be a meta-data key
-                        i = self.__metadata_key_to_bool_vector(i)
+                if isinstance(token, str) and token.isdigit():
+                    # 'token' should be a criteria value
+                    operand: SelectionOperand = float(token)
+                elif isinstance(token, str):
+                    # 'token' should be a meta-data key
+                    operand = self.__metadata_key_to_bool_vector(token)
+                else:
+                    operand = cast(SelectionOperand, token)
 
-                stack.append(i)
+                stack.append(operand)
 
-        selected_index = stack.pop()
+        selected_mask = stack.pop()
 
         # If buf_sel still has an element, `select_index` should be an order vector.
         # Select N elements based on the order vector.
         if buf_sel:
             num_sel = buf_sel.pop()
-            selected_index = [n < num_sel for n in selected_index]
+            selected_mask = np.array([n < num_sel for n in cast(np.ndarray, selected_mask)])
 
         # Very dirty solution
-        selected_index = np.array(selected_index) == True  # Should use "==" instead of "is" here.
+        selected_mask = np.array(selected_mask) == True  # Should use "==" instead of "is" here.
 
         if return_index:
-            return self.dataset[:, np.array(selected_index)], selected_index
+            return self.dataset[:, np.array(selected_mask)], selected_mask
         else:
-            return self.dataset[:, np.array(selected_index)]
-
-    @__obsoleted_method('select')
-    def select_dataset(self, condition: str, return_index: bool = False, verbose: bool = True) -> Union[np.ndarray, Tuple[np.ndarray, list]]:
+            return self.dataset[:, np.array(selected_mask)]
+
+    @_obsoleted_method('select')
+    def select_dataset(
+        self,
+        condition: str,
+        return_index: bool = False,
+        verbose: bool = True
+    ) -> Union[np.ndarray, Tuple[np.ndarray, np.ndarray]]:
         """Select data (columns) from dataset.
 
         Parameters
@@ -488,8 +614,13 @@ def select_dataset(self, condition: str, return_index: bool = False, verbose: bo
         """
         return self.select(condition, return_index, verbose)
 
-    @__obsoleted_method('select')
-    def select_feature(self, condition: str, return_index: bool = False, verbose: bool = True) -> Union[np.ndarray, Tuple[np.ndarray, list]]:
+    @_obsoleted_method('select')
+    def select_feature(
+        self,
+        condition: str,
+        return_index: bool = False,
+        verbose: bool = True
+    ) -> Union[np.ndarray, Tuple[np.ndarray, np.ndarray]]:
         """Select data (columns) from dataset.
 
         Parameters
@@ -531,7 +662,7 @@ def get(self, key: Optional[str] = None) -> np.ndarray:
             query = '%s = 1' % key
             return self.select(query, return_index=False, verbose=False)
 
-    @__obsoleted_method('get')
+    @_obsoleted_method('get')
     def get_dataset(self, key: Optional[str] = None) -> np.ndarray:
         """Get dataset.
 
@@ -540,7 +671,15 @@ def get_dataset(self, key: Optional[str] = None) -> np.ndarray:
         """
         return self.get(key)
 
-    def get_metadata(self, key: str, where: Optional[str] = None) -> np.ndarray:
+    @overload
+    def get_metadata(self, key: str, where: str) -> np.ndarray:
+        ...
+
+    @overload
+    def get_metadata(self, key: str, where: None = None) -> Optional[np.ndarray]:
+        ...
+
+    def get_metadata(self, key: str, where: Optional[str] = None) -> Optional[np.ndarray]:
         """Get value of meta-data specified by `key`.
 
         Parameters
@@ -553,16 +692,26 @@ def get_metadata(self, key: str, where: Optional[str] = None) -> np.ndarray:
 
         Returns
         -------
-        numpy.ndarray
+        numpy.ndarray or None
+            Meta-data value. If `key` is not found and `where` is not
+            given, None is returned.
+
+        Raises
+        ------
+        ValueError
+            If `key` is not found and `where` is given.
         """
         md = self.metadata.get(key, 'value')
 
-        if where is not None:
-            # Mask the metadata array with columns specified with `where`
-            ind = self.metadata.get(where, 'value') == 1
-            md = md[ind]
+        if where is None:
+            return md
+
+        if md is None:
+            raise ValueError(f"Meta-data key '{key}' not found.")
 
-        return md
+        # Mask the metadata array with columns specified with `where`
+        ind = self.metadata.get(where, 'value') == 1
+        return cast(np.ndarray, md[ind])
 
     def show_metadata(self) -> None:
         """Show all the key and description in metadata."""
@@ -580,7 +729,7 @@ def show_metadata(self) -> None:
 
     # Value-label map --------------------------------------------------------
 
-    def get_labels(self, key: str) -> list:
+    def get_labels(self, key: str) -> List[str]:
         """Get `key` as labels."""
         if key not in self.__vmap:
             raise ValueError('Key not found in vmap: %s' % key)
@@ -598,17 +747,18 @@ def get_label(self, key: str) -> list:
         """Get `key` as labels."""
         return self.get_labels(key)
 
-    def get_vmap(self, key: str) -> dict:
-        """Returns vmap of `key`."""
+    def get_vmap(self, key: str) -> Dict[float, str]:
+        """Return the value-label map for a metadata key."""
         if key in self.__vmap:
             return self.__vmap[key]
         else:
             raise ValueError('%s not found in vmap' % key)
 
-    def get_vmap_keys(self):
+    def get_vmap_keys(self) -> KeysView[str]:
+        """Return keys of value-label maps."""
         return self.__vmap.keys()
 
-    def add_vmap(self, key: str, vmap: dict) -> None:
+    def add_vmap(self, key: str, vmap: Dict[float, str]) -> None:
         """Add vmap."""
         if key not in self.__metadata.key:
             raise ValueError('%s not found in metadata.' % key)
@@ -633,7 +783,7 @@ def add_vmap(self, key: str, vmap: dict) -> None:
 
         return None
 
-    def __get_act_vmap(self, key: str, vmap: dict) -> dict:
+    def __get_act_vmap(self, key: str, vmap: Dict[float, str]) -> Dict[float, str]:
         values = np.unique(self.get(key))
         try:
             vmap_add = {}
@@ -645,7 +795,7 @@ def __get_act_vmap(self, key: str, vmap: dict) -> dict:
             raise ValueError('Invalid vmap: label for %f not found.' % val) from err
         return vmap_add
 
-    def __check_vmap_consistency(self, vmap_new: dict, vmap_old: dict) -> bool:
+    def __check_vmap_consistency(self, vmap_new: Dict[float, str], vmap_old: Dict[float, str]) -> bool:
         for key in vmap_new.keys():
             if key not in vmap_old:
                 continue
@@ -680,6 +830,8 @@ def save(self, file_name: str, file_type: Optional[str] = None) -> None:
         callstack = []
         callstack_code = []
         f = inspect.currentframe()
+        if f is None:
+            raise RuntimeError('Failed to get the current frame for call stack information.')
         while True:
             f = f.f_back
             if f is None:
@@ -717,14 +869,18 @@ def __metadata_key_to_bool_vector(self, key: str) -> np.ndarray:
         keys = [k for k in self.metadata.key if re.match(key_esc, k)]
         if len(keys) == 0:
             raise RuntimeError('Meta-data %s not found' % key)
+        # `keys` only contains existing meta-data keys, so `get_metadata`
+        # never returns None here.
         vals = np.vstack([
-            self.get_metadata(k)
+            cast(np.ndarray, self.get_metadata(k))
             for k in keys])
         vals = (vals == 1)
         vec = np.sum(vals, axis=0).astype(bool)
-        return vec
+        return cast(np.ndarray, vec)
 
     def __get_order(self, v: np.ndarray, sort_order: str = 'descend') -> np.ndarray:
+        if sort_order != "descend":
+            raise NotImplementedError('Only "descend" is implemented for `sort_order`.')
 
         # 'np.nan' comes to the last of an acending series, and thus the top of a decending series.
         # To avoid that, convert 'np.nan' to -Inf.
@@ -747,7 +903,7 @@ def __get_top_elm_from_order(self, order: np.ndarray, n: int) -> np.ndarray:
 
         return index
 
-    def __save_h5(self, file_name: str, header: Optional[dict] = None) -> None:
+    def __save_h5(self, file_name: str, header: Optional[Dict[str, Any]] = None) -> None:
         """Save data in HDF5 format (*.h5)."""
         with h5py.File(file_name, 'w') as h5file:
             # dataset
@@ -766,18 +922,27 @@ def __save_h5(self, file_name: str, header: Optional[dict] = None) -> None:
             # header
             if header is not None:
                 h5file.create_group('/header')
-                for k, v in header.items():
-                    if isinstance(v, list):
-                        h5file.create_dataset('/header/' + k, data=[self.__to_bytes(x) for x in v])
+                for header_key, header_value in header.items():
+                    if isinstance(header_value, list):
+                        h5file.create_dataset(
+                            '/header/' + header_key,
+                            data=[self.__to_bytes(x) for x in header_value]
+                        )
                     else:
-                        h5file.create_dataset('/header/' + k, data=self.__to_bytes(v)) # FIXME: save unicode str as is
+                        h5file.create_dataset(
+                            '/header/' + header_key,
+                            data=self.__to_bytes(header_value)
+                        ) # FIXME: save unicode str as is
 
             # vmap
             h5file.create_group('/vmap')
-            for mk, vm in self.__vmap.items():
-                h5file.create_group('/vmap/' + mk)
-                for k, v in vm.items():
-                    h5file.create_dataset('/vmap/' + mk + '/' + str(k), data=self.__to_bytes(v)) # FIXME: save unicode str as is
+            for metadata_key, value_map in self.__vmap.items():
+                h5file.create_group('/vmap/' + metadata_key)
+                for value_key, label in value_map.items():
+                    h5file.create_dataset(
+                        '/vmap/' + metadata_key + '/' + str(value_key),
+                        data=self.__to_bytes(label)
+                    ) # FIXME: save unicode str as is
 
     def __load_mat(self, load_filename: str) -> None:
         """Load dataset and metadata from Matlab file."""
@@ -842,27 +1007,49 @@ def __load_h5(self, load_filename: str) -> None:
         self.__metadata.value = md_values
         self.__metadata.description = md_descs
 
-    def __to_unicode(self, s):
-        """Convert s (bytes) to unicode str."""
+    @overload
+    def __to_unicode(self, s: bytes) -> str:
+        ...
+
+    @overload
+    def __to_unicode(self, s: _T) -> _T:
+        ...
+
+    def __to_unicode(self, s: Any) -> Any:
+        """Convert `s` to a unicode str only when it is bytes.
+
+        Values of any other type (e.g., numeric header values) are returned
+        unchanged.
+        """
         if isinstance(s, bytes):
             return s.decode('utf-8')
         return s
 
-    def __to_bytes(self, s):
-        """Convert s (unicode str) to bytes."""
+    @overload
+    def __to_bytes(self, s: str) -> bytes:
+        ...
+
+    @overload
+    def __to_bytes(self, s: _T) -> _T:
+        ...
+
+    def __to_bytes(self, s: Any) -> Any:
+        """Convert `s` to bytes only when it is a unicode str.
+
+        Values of any other type (e.g., numeric header values) are returned
+        unchanged.
+        """
         if isinstance(s, str):
             return s.encode('utf-8')
         return s
 
-    def __get_filetype(self, file_name: str):
+    def __get_filetype(self, file_name: str) -> Literal['Matlab', 'HDF5']:
         """Return the type of `file_name` based on the file extension."""
         _, ext = os.path.splitext(file_name)
 
         if ext == ".mat":
-            file_type = "Matlab"
+            return "Matlab"
         elif ext == ".h5":
-            file_type = "HDF5"
+            return "HDF5"
         else:
             raise ValueError("Unknown file extension: %s" % (ext))
-
-        return file_type
diff --git a/bdpy/bdata/featureselector.py b/bdpy/bdata/featureselector.py
index f3fc5937..9819fd09 100644
--- a/bdpy/bdata/featureselector.py
+++ b/bdpy/bdata/featureselector.py
@@ -1,13 +1,15 @@
 """
-Feature selector class
+Feature selector class.
 
 This file is a part of BdPy
 """
 
 
+from typing import ClassVar, Dict, List, Optional, Tuple
+
+
 class FeatureSelector(object):
-    """
-    Feature selector class
+    """Feature selector class.
 
     Parameters
     ----------
@@ -25,31 +27,33 @@ class FeatureSelector(object):
     """
 
     # Class variables ##################
-    signs = ('(', ')')
-    operators = ('=', '|', '&', '+', '-', '@')
-
-    __op_order = {'=': 10,
-                  '|': 5,
-                  '&': 5,
-                  '+': 5,
-                  '-': 5,
-                  '@': 3,
-                  '(': -1,
-                  ')': -1}
+    signs: ClassVar[Tuple[str, ...]] = ('(', ')')
+    operators: ClassVar[Tuple[str, ...]] = ('=', '|', '&', '+', '-', '@')
+
+    __op_order: ClassVar[Dict[str, int]] = {
+        '=': 10,
+        '|': 5,
+        '&': 5,
+        '+': 5,
+        '-': 5,
+        '@': 3,
+        '(': -1,
+        ')': -1,
+    }
 
     # Methods ##########################
 
-    def __init__(self, expression):
+    def __init__(self, expression: str) -> None:
         self.expression = expression
         self.token = self.lexical_analysis(self.expression)
         self.rpn = self.parse(self.token)
 
-        self.index = None
+        self.index: Optional[int] = None
 
-    def lexical_analysis(self, expression):
-        """Lexical analyser"""
+    def lexical_analysis(self, expression: str) -> Tuple[str, ...]:
+        """Tokenize selection command."""
         str_buf = ''
-        output_buf = []
+        output_buf: List[str] = []
 
         i = 0
         while i < len(expression):
@@ -91,10 +95,10 @@ def lexical_analysis(self, expression):
 
         return tuple(output_buf)
 
-    def parse(self, token_list):
-        """Parser for selection command"""
-        out_que = []
-        op_stack = []
+    def parse(self, token_list: Tuple[str, ...]) -> Tuple[str, ...]:
+        """Parse selection command."""
+        out_que: List[str] = []
+        op_stack: List[str] = []
 
         for token in token_list:
 
diff --git a/bdpy/bdata/metadata.py b/bdpy/bdata/metadata.py
index 2f927c57..1705663b 100644
--- a/bdpy/bdata/metadata.py
+++ b/bdpy/bdata/metadata.py
@@ -1,22 +1,31 @@
-"""
-MetaData class
+"""MetaData class.
 
 This file is a part of BdPy
 """
 
 
+from typing import Callable, List, Optional, Sequence, Union, overload
+
 import numpy as np
+from typing_extensions import Literal
+
+MetaDataValue = Union[np.ndarray, Sequence[float]]
+MetaDataUpdater = Callable[[np.ndarray, np.ndarray], MetaDataValue]
 
 
 class MetaData(object):
-    """
-    MetaData class
+    """MetaData class.
 
     'MetaData' is a list of dictionaries. Each element has three keys: 'key',
     'value', and 'description'.
     """
 
-    def __init__(self, key=None, value=None, description=None):
+    def __init__(
+        self,
+        key: Optional[List[str]] = None,
+        value: Optional[np.ndarray] = None,
+        description: Optional[List[str]] = None,
+    ) -> None:
         if key is None:
             key = []
         if value is None:
@@ -29,59 +38,73 @@ def __init__(self, key=None, value=None, description=None):
         self.__description = description
 
     @property
-    def key(self):
+    def key(self) -> List[str]:
+        """Meta-data keys."""
         return self.__key
 
     @key.setter
-    def key(self, x):
+    def key(self, x: List[str]) -> None:
         self.__key = x
 
     @property
-    def value(self):
+    def value(self) -> np.ndarray:
+        """Meta-data values."""
         return self.__value
 
     @value.setter
-    def value(self, x):
+    def value(self, x: np.ndarray) -> None:
         self.__value = x
 
     @property
-    def description(self):
+    def description(self) -> List[str]:
+        """Meta-data descriptions."""
         return self.__description
 
     @description.setter
-    def description(self, x):
+    def description(self, x: List[str]) -> None:
         self.__description = x
 
-    def set(self, key, value, description, updater=None):
-        """
-        Set meta-data with `key`, `description`, and `value`
+    def set(
+        self,
+        key: str,
+        value: Optional[MetaDataValue],
+        description: str,
+        updater: Optional[MetaDataUpdater] = None,
+    ) -> None:
+        """Set meta-data with `key`, `description`, and `value`.
 
         Parameters
         ----------
         key : str
             Meta-data key
-        value : array_like
-            Meta-data value
+        value : array_like or None
+            Meta-data value. If None, only the description of existing
+            meta-data is updated.
         description : str
             Meta-data description
         updater : function
             Function applied to meta-data value when meta-data named `key` already exists.
             It should take two args: new and old meta-data values.
+
+        Raises
+        ------
+        ValueError
+            If `value` is None and meta-data named `key` does not exist.
         """
         # If `value` is None, `set` does not update the value.
         is_novalue = True if value is None else False
 
-        value = np.array(value)
+        value_array = np.array(value)
 
         if key in self.__key:
             # Update existing metadata
 
-            ind = [i for i, k in enumerate(self.__key) if k == key]
+            indices = [i for i, k in enumerate(self.__key) if k == key]
 
-            if len(ind) > 1:
+            if len(indices) > 1:
                 raise ValueError('Multiple meta-data with the same key is not supported')
 
-            ind = ind[0]
+            ind = indices[0]
 
             self.__description[ind] = description
 
@@ -89,33 +112,47 @@ def set(self, key, value, description, updater=None):
             if is_novalue:
                 return None
 
-            if value.shape[0] > self.get_value_len():
-                cols = np.empty((self.__value.shape[0], value.shape[0] - self.get_value_len()))
+            if value_array.shape[0] > self.get_value_len():
+                cols = np.empty((self.__value.shape[0], value_array.shape[0] - self.get_value_len()))
                 cols[:] = np.nan
 
                 self.__value = np.hstack([self.__value, cols])
 
             if updater is None:
-                self.__value[ind, :] = value
+                self.__value[ind, :] = value_array
             else:
-                self.__value[ind, :] = np.array(updater(value, self.__value[ind, :]), dtype=float)
+                self.__value[ind, :] = np.array(updater(value_array, self.__value[ind, :]), dtype=float)
         else:
             # Add new metadata
+            if is_novalue:
+                raise ValueError(f"Cannot add new meta-data '{key}' without a value.")
+
             self.__key.append(key)
             self.__description.append(description)
 
-            if value.shape[0] > self.get_value_len():
-                cols = np.empty((self.__value.shape[0], value.shape[0] - self.get_value_len()))
+            if value_array.shape[0] > self.get_value_len():
+                cols = np.empty((self.__value.shape[0], value_array.shape[0] - self.get_value_len()))
                 cols[:] = np.nan
 
                 self.__value = np.hstack([self.__value, cols])
 
-            self.__value = np.vstack([self.__value, value])
+            self.__value = np.vstack([self.__value, value_array])
 
 
-    def get(self, key, field):
-        """
-        Returns meta-data specified by `key`
+    @overload
+    def get(self, key: str, field: Literal['value']) -> Optional[np.ndarray]:
+        ...
+
+    @overload
+    def get(self, key: str, field: Literal['description']) -> Optional[str]:
+        ...
+
+    @overload
+    def get(self, key: str, field: str) -> Union[np.ndarray, str, None]:
+        ...
+
+    def get(self, key: str, field: str) -> Union[np.ndarray, str, None]:
+        """Return meta-data specified by `key`.
 
         Parameters
         ----------
@@ -141,12 +178,14 @@ def get(self, key, field):
         if field == 'description':
             return self.__description[ind]
 
+        return None
+
 
-    def get_value_len(self):
-        """Returns length of meta-data value"""
+    def get_value_len(self) -> int:
+        """Return length of meta-data value."""
         return self.__value.shape[1]
 
 
-    def keylist(self):
-        """Returns a list of keys"""
+    def keylist(self) -> List[str]:
+        """Return a list of keys."""
         return self.__key
diff --git a/bdpy/bdata/utils.py b/bdpy/bdata/utils.py
index 1aa70d39..ea7c027f 100644
--- a/bdpy/bdata/utils.py
+++ b/bdpy/bdata/utils.py
@@ -2,14 +2,20 @@
 
 
 import copy
-from typing import List
+from typing import Dict, List, Optional, Sequence, cast
 
 import numpy as np
+from typing_extensions import Literal
 
 from .bdata import BData
 
 
-def vstack(bdata_list, successive=[], metadata_merge='strict', ignore_metadata_description=False):
+def vstack(
+    bdata_list: Sequence[BData],
+    successive: Optional[Sequence[str]] = None,
+    metadata_merge: Literal['strict', 'minimal'] = 'strict',
+    ignore_metadata_description: bool = False,
+) -> BData:
     """Concatenate datasets vertically.
 
     Currently, `concat_dataset` does not validate the consistency of meta-data
@@ -39,6 +45,9 @@ def vstack(bdata_list, successive=[], metadata_merge='strict', ignore_metadata_d
 
         data = vstack([data0, data1, data2], successive=['Session', 'Run', 'Block'])
     """
+    if successive is None:
+        successive = []
+
     suc_cols = {s : 0 for s in successive}
 
     dat = BData()  # Concatenated BData
@@ -80,10 +89,10 @@ def vstack(bdata_list, successive=[], metadata_merge='strict', ignore_metadata_d
                         raise ValueError('Inconsistent meta-data description (%s)' % mkey)
                     try:
                         np.testing.assert_equal(d0_value, d1_value)
-                    except AssertionError:
-                        raise ValueError('Inconsistent meta-data value (%s)' % mkey)
-                    shared_mdesc.append(d0_desc)
-                    shared_mvalue_lst.append(d0_value)
+                    except AssertionError as err:
+                        raise ValueError('Inconsistent meta-data value (%s)' % mkey) from err
+                    shared_mdesc.append(cast(str, d0_desc))
+                    shared_mvalue_lst.append(cast(np.ndarray, d0_value))
                 shared_mvalue = np.vstack(shared_mvalue_lst)
 
                 dat.metadata.key = shared_mkeys
@@ -108,7 +117,7 @@ def vstack(bdata_list, successive=[], metadata_merge='strict', ignore_metadata_d
     return dat
 
 
-def resolve_vmap(bdata_list):
+def resolve_vmap(bdata_list: List[BData]) -> List[BData]:
     """Replace the conflicting vmaps for multiple bdata with non-conflicting vmaps.
 
     Parameters
@@ -126,7 +135,7 @@ def resolve_vmap(bdata_list):
 
     # Check each vmap key.
     for vmap_key in vmap_keys:
-        new_vmap = {}
+        new_vmap: Dict[float, str] = {}
         # Check each bdata vmap.
         for ds in bdata_list:
             vmap = ds.get_vmap(vmap_key)
@@ -166,13 +175,13 @@ def resolve_vmap(bdata_list):
             vmap = ds.get_vmap(vmap_key)
             if not np.array_equal(sorted(list(vmap.keys())), sorted(list(new_vmap.keys()))):
                 # If the present vmap is different from new_vmap, update it.
-                ds._BData__vmap[vmap_key] = new_vmap # BDataクラスにvmapのsetterがあると良い
+                ds._BData__vmap[vmap_key] = new_vmap # type: ignore[attr-defined] # BDataクラスにvmapのsetterがあると良い
 
     return bdata_list
 
 
-def concat_dataset(data_list, successive=[]):
-    """Concatenate datasets
+def concat_dataset(data_list: Sequence[BData], successive: Optional[Sequence[str]] = None) -> BData:
+    """Concatenate datasets.
 
     Currently, `concat_dataset` does not validate the consistency of meta-data
     among data.
@@ -198,7 +207,7 @@ def concat_dataset(data_list, successive=[]):
     return vstack(data_list, successive=successive)
 
 
-def metadata_equal(d0, d1, strict=False):
+def metadata_equal(d0: BData, d1: BData, strict: bool = False) -> bool:
     """Check whether `d0` and `d1` share the same meta-data.
 
     Parameters
@@ -235,8 +244,10 @@ def metadata_equal(d0, d1, strict=False):
         return False
 
     for mkey in d0_mkeys:
-        d0_mdesc, d1_mdesc = d0.metadata.get(mkey, 'description'), d1.metadata.get(mkey, 'description')
-        d0_mval, d1_mval = d0.metadata.get(mkey, 'value'), d1.metadata.get(mkey, 'value')
+        d0_mdesc = cast(str, d0.metadata.get(mkey, 'description'))
+        d1_mdesc = cast(str, d1.metadata.get(mkey, 'description'))
+        d0_mval = cast(np.ndarray, d0.metadata.get(mkey, 'value'))
+        d1_mval = cast(np.ndarray, d1.metadata.get(mkey, 'value'))
 
         if not d0_mdesc == d1_mdesc:
             return False
diff --git a/tests/bdata/test_bdata.py b/tests/bdata/test_bdata.py
index f0c52fb9..fe20a65c 100644
--- a/tests/bdata/test_bdata.py
+++ b/tests/bdata/test_bdata.py
@@ -1,7 +1,10 @@
 '''Tests for bdpy.bdata.bdata.'''
 
 
+import os
+import tempfile
 import unittest
+import warnings
 
 import numpy as np
 from numpy.testing import assert_array_equal
@@ -42,6 +45,18 @@ def test_add_get(self):
         assert_array_equal(b.get_metadata('Data_Y'), np.array([np.nan] * 10 + [1] * 8 + [np.nan] * 20))
         assert_array_equal(b.get_metadata('Data_Z'), np.array([np.nan] * 10 + [np.nan] * 8 + [1] * 20))
 
+    def test_add_dataset_is_obsoleted(self):
+        '''Test that BData.add_dataset warns once and delegates to add.'''
+        b = BData()
+
+        with warnings.catch_warnings(record=True) as recorded_warnings:
+            warnings.simplefilter('always')
+            b.add_dataset(np.ones((2, 3)), 'Data')
+
+        self.assertEqual(len(recorded_warnings), 1)
+        self.assertIn("'add_dataset' is obsoleted", str(recorded_warnings[0].message))
+        assert_array_equal(b.get('Data'), np.ones((2, 3)))
+
     def test_metadata_add_get(self):
         '''Test for add/get_metadata.'''
 
@@ -88,6 +103,29 @@ def test_metadata_add_get_where(self):
         assert_array_equal(b.get_metadata('Metadata_A', where='Data_X'), metadata_a)
         assert_array_equal(b.get_metadata('Metadata_B', where='Data_Y'), metadata_b)
 
+    def test_get_metadata_notfound(self):
+        '''Test for BData.get_metadata with a missing key.'''
+        b = BData()
+        b.add(np.ones((2, 3)), 'Data')
+
+        self.assertIsNone(b.get_metadata('Metadata_NotFound'))
+
+    def test_get_metadata_notfound_where(self):
+        '''Test for BData.get_metadata with a missing key and where option.'''
+        b = BData()
+        b.add(np.ones((2, 3)), 'Data')
+
+        with self.assertRaises(ValueError):
+            b.get_metadata('Metadata_NotFound', where='Data')
+
+    def test_update_notfound(self):
+        '''Test for BData.update with a missing key.'''
+        b = BData()
+        b.add(np.ones((2, 3)), 'Data')
+
+        with self.assertRaises(ValueError):
+            b.update('Metadata_NotFound', np.zeros((2, 3)))
+
     def test_set_metadatadescription_1(self):
         '''Test for set_metadatadescription.'''
 
@@ -138,6 +176,92 @@ def test_select(self):
         assert_array_equal(b.select('ROI_0:5 + ROI_3:8'), data_x[:, 0:8])
         assert_array_equal(b.select('ROI_0:5 - ROI_3:8'), data_x[:, 0:3])
 
+    def test_select_return_index(self):
+        '''Test for BData.select with return_index=True.'''
+        data_x = np.arange(20, dtype=float).reshape(4, 5)
+        data_y = np.arange(8, dtype=float).reshape(4, 2)
+
+        b = BData()
+
+        # BData.add stacks arrays column-wise: Data_X columns come first,
+        # followed by Data_Y columns.
+        b.add(data_x, 'Data_X')
+        b.add(data_y, 'Data_Y')
+
+        # add_metadata(..., where='Data_X') defines the ROI only inside the
+        # Data_X column group, not against the whole dataset.
+        b.add_metadata('ROI_1:4', [0, 1, 1, 1, 0], where='Data_X')
+
+        selected_data, selected_index = b.select('ROI_1:4', return_index=True)
+
+        expected_data_x_index = np.array([False, True, True, True, False])
+        expected_data_y_index = np.array([False, False])
+        expected_dataset_index = np.hstack([expected_data_x_index, expected_data_y_index])
+
+        assert_array_equal(selected_data, data_x[:, 1:4])
+        self.assertEqual(len(selected_index), b.dataset.shape[1])
+        assert_array_equal(selected_index, expected_dataset_index)
+        self.assertEqual(selected_index.dtype, np.dtype(bool))
+
+    def test_applyfunc_with_selected_columns(self):
+        '''Test for BData.applyfunc with a selected column group.'''
+        data_x = np.arange(6, dtype=float).reshape(3, 2)
+        data_y = np.arange(3, dtype=float).reshape(3, 1)
+
+        b = BData()
+        b.add(data_x, 'Data_X')
+        b.add(data_y, 'Data_Y')
+
+        b.applyfunc(lambda x: x + 10, where='Data_X')
+
+        assert_array_equal(b.get('Data_X'), data_x + 10)
+        assert_array_equal(b.get('Data_Y'), data_y)
+
+    def test_applyfunc_tuple_result_reindexes_all_columns(self):
+        '''Test that BData.applyfunc reindexes all columns when func returns an index map.'''
+        data_x = np.arange(6, dtype=float).reshape(3, 2)
+        data_y = np.arange(3, dtype=float).reshape(3, 1)
+        row_index = np.array([2, 0])
+
+        b = BData()
+        b.add(data_x, 'Data_X')
+        b.add(data_y, 'Data_Y')
+
+        # A tuple result means that the selected column group is replaced by
+        # the first element, and every non-selected column follows the returned
+        # row index map so that rows remain aligned across the whole dataset.
+        b.applyfunc(lambda x: (x[row_index], row_index), where='Data_X')
+
+        assert_array_equal(b.get('Data_X'), data_x[row_index])
+        assert_array_equal(b.get('Data_Y'), data_y[row_index])
+
+    def test_save_load_hdf5_header_and_vmap(self):
+        '''Test for HDF5 roundtrip of header values and vmap.'''
+        data = np.arange(6, dtype=float).reshape(3, 2)
+        label = np.array([1, 2, 1], dtype=float).reshape(3, 1)
+
+        bdata = BData()
+        bdata.add(data, 'Data')
+        bdata.add(label, 'Label')
+        bdata.add_vmap('Label', {1: 'label-1', 2: 'label-2'})
+        bdata.update_header({
+            'source': 'manual',
+            'indices': [1, 2],
+            'scale': 1.5,
+        })
+
+        with tempfile.TemporaryDirectory() as temp_dir:
+            h5_path = os.path.join(temp_dir, 'test_bdata.h5')
+            bdata.save(h5_path, 'HDF5')
+            loaded_bdata = BData(h5_path, 'HDF5')
+
+        assert_array_equal(loaded_bdata.get('Data'), data)
+        assert_array_equal(loaded_bdata.get('Label'), label)
+        self.assertEqual(loaded_bdata.get_vmap('Label'), {1.0: 'label-1', 2.0: 'label-2'})
+        self.assertEqual(loaded_bdata.header['source'], 'manual')
+        self.assertEqual(loaded_bdata.header['indices'], [1, 2])
+        self.assertEqual(loaded_bdata.header['scale'], 1.5)
+
     # Tests for vmap
     def test_vmap_add_get(self):
         bdata = BData()
@@ -152,6 +276,7 @@ def test_vmap_add_get(self):
 
         bdata.add_vmap('Label', label_map)
         assert bdata.get_vmap('Label') == label_map
+        self.assertEqual(set(bdata.get_vmap_keys()), {'Label'})
 
         # Get labels
         np.testing.assert_array_equal(bdata.get_label('Label'), label)
@@ -250,4 +375,4 @@ def test_vmap_add_invalid_name_vmap(self):
 
 
 if __name__ == "__main__":
-    unittest.main()
\ No newline at end of file
+    unittest.main()
diff --git a/tests/bdata/test_metadata.py b/tests/bdata/test_metadata.py
index 7d3d6b65..2cc75149 100644
--- a/tests/bdata/test_metadata.py
+++ b/tests/bdata/test_metadata.py
@@ -79,6 +79,25 @@ def test_set_get_update(self):
         assert_array_equal(md.get('MetaData_B', 'value'), [0] * 10 + [1] * 5)
         assert_array_equal(md.get('MetaData_B', 'description'), 'Test metadata B')
 
+    def test_set_description_only(self):
+        '''Test for MetaData.set(); updating description without value.'''
+        md = metadata.MetaData()
+        md.set('MetaData_A', [1] * 10 + [0] * 5, 'Test metadata A')
+
+        md.set('MetaData_A', None, 'Updated metadata A')
+
+        assert_array_equal(md.get('MetaData_A', 'value'), [1] * 10 + [0] * 5)
+        assert_array_equal(md.get('MetaData_A', 'description'), 'Updated metadata A')
+
+    def test_set_novalue_new_key(self):
+        '''Test for MetaData.set(); value=None is not allowed for a new key.'''
+        md = metadata.MetaData()
+
+        with self.assertRaises(ValueError):
+            md.set('MetaData_A', None, 'Test metadata A')
+
+        self.assertEqual(md.keylist(), [])
+
     def test_get_notfound(self):
         '''Test for MetaData.get(); key not found case.'''
         md = metadata.MetaData()
@@ -88,6 +107,13 @@ def test_get_notfound(self):
         assert_array_equal(md.get('MetaData_NotFound', 'value'), None)
         assert_array_equal(md.get('MetaData_NotFound', 'description'), None)
 
+    def test_get_unknown_field(self):
+        '''Test for MetaData.get(); unknown field case.'''
+        md = metadata.MetaData()
+        md.set('MetaData_A', [1] * 10 + [0] * 5, 'Test metadata A')
+
+        assert_array_equal(md.get('MetaData_A', 'unknown'), None)
+
     def test_get_value_len(self):
         '''Test for get_value_len().'''
         md = metadata.MetaData()
diff --git a/tests/bdata/test_utils.py b/tests/bdata/test_utils.py
index 81999bae..d3dd4994 100644
--- a/tests/bdata/test_utils.py
+++ b/tests/bdata/test_utils.py
@@ -62,6 +62,21 @@ def test_vstack_successive(self):
                                       np.vstack([x0_run,
                                                  x1_run + len(x0_run)]))
 
+    def test_vstack_successive_none(self):
+        x0_data = np.random.rand(10, 20)
+        x1_data = np.random.rand(10, 20)
+
+        bdata0 = BData()
+        bdata0.add(x0_data, 'Data')
+
+        bdata1 = BData()
+        bdata1.add(x1_data, 'Data')
+
+        bdata_merged = vstack([bdata0, bdata1], successive=None)
+
+        np.testing.assert_array_equal(bdata_merged.select('Data'),
+                                      np.vstack([x0_data, x1_data]))
+
     def test_vstack_minimal(self):
         x0_data = np.random.rand(5, 10)
         x0_label = np.random.rand(5, 1)
@@ -93,6 +108,19 @@ def test_vstack_minimal(self):
         self.assertFalse('key only in 0' in bdata_merged.metadata.key)
         self.assertFalse('key only in 1' in bdata_merged.metadata.key)
 
+    def test_vstack_unknown_metadata_merge(self):
+        x0_data = np.random.rand(5, 10)
+        x1_data = np.random.rand(5, 10)
+
+        bdata0 = BData()
+        bdata0.add(x0_data, 'Data')
+
+        bdata1 = BData()
+        bdata1.add(x1_data, 'Data')
+
+        with self.assertRaises(ValueError):
+            vstack([bdata0, bdata1], metadata_merge='unknown')
+
     def test_vstack_vmap(self):
         x0_data = np.random.rand(10, 20)
         x0_label = np.random.permutation(np.arange(10)).reshape(10, 1) + 1
@@ -271,4 +299,4 @@ def test_metadata_equal_loose_notequal_value(self):
 
 
 if __name__ == "__main__":
-    unittest.main()
\ No newline at end of file
+    unittest.main()