diff --git a/CHANGES.md b/CHANGES.md
index d679dc27..b45a34bd 100644
--- a/CHANGES.md
+++ b/CHANGES.md
@@ -1,6 +1,8 @@
 
 # 4.0.0 (TBD)
 
+* add `apply_expression` method in `rio_tiler.models.ImageData` class
+
 **breaking changes**
 
 * remove python 3.7 support
@@ -8,6 +10,99 @@
 * `rio_tiler.readers.read()`, `rio_tiler.readers.part()`, `rio_tiler.readers.preview()` now return a ImageData object
 * remove `minzoom` and `maxzoom` attribute in `rio_tiler.io.SpatialMixin` base class
 * remove `minzoom` and `maxzoom` attribute in `rio_tiler.io.COGReader` (now defined as properties).
+* use `b` prefix for band names in `rio_tiler.models.ImageData` class (and in rio-tiler's readers)
+    ```python
+    # before
+    with COGReader("cog.tif") as cog:
+        img = cog.read()
+        print(cog.band_names)
+        >>> ["1", "2", "3"]
+
+        print(cog.info().band_metadata)
+        >>> [("1", {}), ("2", {}), ("3", {})]
+
+        print(cog.info().band_descriptions)
+        >>> [("1", ""), ("2", ""), ("3", "")]
+
+        print(list(cog.statistics()))
+        >>> ["1", "2", "3"]
+
+    # now
+    with COGReader("cog.tif") as cog:
+        img = cog.read()
+        print(img.band_names)
+        >>> ["b1", "b2", "b3"]
+
+        print(cog.info().band_metadata)
+        >>> [("b1", {}), ("b2", {}), ("b3", {})]
+
+        print(cog.info().band_descriptions)
+        >>> [("b1", ""), ("b2", ""), ("b3", "")]
+
+        print(list(cog.statistics()))
+        >>> ["b1", "b2", "b3"]
+
+    with STACReader("stac.json") as stac:
+        print(stac.tile(701, 102, 8, assets=("green", "red")).band_names)
+        >>> ["green_b1", "red_b1"]
+    ```
+
+* depreciate `asset_expression` in MultiBaseReader. Use of expression is now possible
+* `expression` for MultiBaseReader must be in form of `{asset}_b{index}`
+
+    ```python
+    # before
+    with STACReader("stac.json") as stac:
+        stac.tile(701, 102, 8, expression="green/red")
+
+    # now
+    with STACReader("stac.json") as stac:
+        stac.tile(701, 102, 8, expression="green_b1/red_b1")
+    ```
+
+* `rio_tiler.reader.point()` (and all Reader's point methods) now return a **Tuple** of values and band names
+
+    ```python
+    # before
+    with rasterio.open("cog.tif") as src::
+        v = rio_tiler.reader.point(10.20, -42.0)
+        print(v)
+        >>> [0, 0, 0]
+
+    with COGReader("cog.tif") as cog:
+        print(cog.point(10.20, -42.0))
+        >>> [0, 0, 0]
+
+    # now
+    with rasterio.open("cog.tif") as src::
+        v, band_names = rio_tiler.reader.point(10.20, -42)
+        print(v)
+        >>> [0, 0, 0]
+        print(band_names)
+        >>> ["b1", "b2", "b3"]
+
+    with COGReader("cog.tif") as cog:
+        print(cog.point(10.20, -42.0))
+        >>> ([0, 0, 0], ["b1", "b2", "b3"])
+    ```
+
+* `MultiBaseReader.point()` method now returns data as flat (merged) list (instead of a list of list)
+
+    ```python
+    # before
+    with STACReader("stac.json") as stac:
+        pt = stac.point(10.20, -42, assets=("green", "red"))
+        print(pt)
+        >>> [[0], [0]]
+
+    # now
+    with STACReader("stac.json") as stac:
+        pt, names = stac.point(10.20, -42, assets=("green", "red"))
+        print(pt)
+        >>> [0, 0]
+        print(names)
+        >>> ["green_b1", "red_b1"]
+    ```
 
 # 3.1.6 (2022-07-22)
 
diff --git a/rio_tiler/expression.py b/rio_tiler/expression.py
index 7b33c868..d733be08 100644
--- a/rio_tiler/expression.py
+++ b/rio_tiler/expression.py
@@ -2,7 +2,7 @@
 
 import re
 import warnings
-from typing import List, Sequence, Tuple, Union
+from typing import List, Sequence, Tuple
 
 import numexpr
 import numpy
@@ -59,7 +59,7 @@ def get_expression_blocks(expression: str) -> List[str]:
 
 def apply_expression(
     blocks: Sequence[str],
-    bands: Sequence[Union[str, int]],
+    bands: Sequence[str],
     data: numpy.ndarray,
 ) -> numpy.ndarray:
     """Apply rio-tiler expression.
@@ -74,6 +74,11 @@ def apply_expression(
         numpy.array: output data.
 
     """
+    if len(bands) != data.shape[0]:
+        raise ValueError(
+            f"Incompatible number of bands ({bands}) and data shape {data.shape}"
+        )
+
     return numpy.array(
         [
             numpy.nan_to_num(
diff --git a/rio_tiler/io/base.py b/rio_tiler/io/base.py
index 75986be1..5cb5720a 100644
--- a/rio_tiler/io/base.py
+++ b/rio_tiler/io/base.py
@@ -1,6 +1,7 @@
 """rio_tiler.io.base: ABC class for rio-tiler readers."""
 
 import abc
+import itertools
 import re
 import warnings
 from typing import Any, Dict, List, Optional, Sequence, Tuple, Type, Union
@@ -25,6 +26,13 @@
 from ..utils import get_array_statistics
 
 
+def _AssetExpressionWarning():
+    warnings.warn(
+        "asset_expression is deprecated and will be removed in 4.0. Use pure Expression",
+        DeprecationWarning,
+    )
+
+
 @attr.s
 class SpatialMixin:
     """Spatial Info Mixin.
@@ -194,7 +202,7 @@ def preview(self, **kwargs: Any) -> ImageData:
         ...
 
     @abc.abstractmethod
-    def point(self, lon: float, lat: float, **kwargs: Any) -> List:
+    def point(self, lon: float, lat: float, **kwargs: Any) -> Tuple[List, List[str]]:
         """Read a value from a Dataset.
 
         Args:
@@ -203,6 +211,7 @@ def point(self, lon: float, lat: float, **kwargs: Any) -> List:
 
         Returns:
             list: Pixel value per bands/assets.
+            list: band names
 
         """
         ...
@@ -262,8 +271,8 @@ def _get_asset_url(self, asset: str) -> str:
 
     def parse_expression(self, expression: str) -> Tuple:
         """Parse rio-tiler band math expression."""
-        assets = "|".join([rf"\b{asset}\b" for asset in self.assets])
-        _re = re.compile(assets.replace("\\\\", "\\"))
+        assets = "|".join([asset for asset in self.assets])
+        _re = re.compile(rf"\b({assets})_b\d+")
         return tuple(set(re.findall(_re, expression)))
 
     def info(  # type: ignore
@@ -360,7 +369,7 @@ def merged_statistics(  # type: ignore
             assets (sequence of str or str): assets to fetch info from.
             expression (str, optional): rio-tiler expression for the asset list (e.g. asset1/asset2+asset3).
             asset_indexes (dict, optional): Band indexes for each asset (e.g {"asset1": 1, "asset2": (1, 2,)}).
-            asset_expression (dict, optional): rio-tiler expression for each asset (e.g. {"asset1": "b1/b2+b3", "asset2": ...}).
+            asset_expression (dict, optional): rio-tiler expression for each asset (e.g. {"asset1": "b1/b2+b3", "asset2": ...}). **Deprecated**
             categorical (bool): treat input data as categorical data. Defaults to False.
             categories (list of numbers, optional): list of categories to return value for.
             percentiles (list of numbers, optional): list of percentile values to calculate. Defaults to `[2, 98]`.
@@ -373,6 +382,9 @@ def merged_statistics(  # type: ignore
             Dict[str, rio_tiler.models.BandStatistics]: bands statistics.
 
         """
+        if asset_expression:
+            _AssetExpressionWarning()
+
         if not expression:
             if not assets:
                 warnings.warn(
@@ -385,7 +397,6 @@ def merged_statistics(  # type: ignore
             assets=assets,
             expression=expression,
             asset_indexes=asset_indexes,
-            asset_expression=asset_expression,
             max_size=max_size,
             **kwargs,
         )
@@ -425,13 +436,16 @@ def tile(
             assets (sequence of str or str, optional): assets to fetch info from.
             expression (str, optional): rio-tiler expression for the asset list (e.g. asset1/asset2+asset3).
             asset_indexes (dict, optional): Band indexes for each asset (e.g {"asset1": 1, "asset2": (1, 2,)}).
-            asset_expression (dict, optional): rio-tiler expression for each asset (e.g. {"asset1": "b1/b2+b3", "asset2": ...}).
+            asset_expression (dict, optional): rio-tiler expression for each asset (e.g. {"asset1": "b1/b2+b3", "asset2": ...}). **Deprecated**
             kwargs (optional): Options to forward to the `self.reader.tile` method.
 
         Returns:
             rio_tiler.models.ImageData: ImageData instance with data, mask and tile spatial info.
 
         """
+        if asset_expression:
+            _AssetExpressionWarning()
+
         if not self.tile_exists(tile_x, tile_y, tile_z):
             raise TileOutsideBounds(
                 f"Tile {tile_z}/{tile_x}/{tile_y} is outside image bounds"
@@ -455,35 +469,20 @@ def tile(
             )
 
         asset_indexes = asset_indexes or {}
-        asset_expression = asset_expression or {}
 
         def _reader(asset: str, *args: Any, **kwargs: Any) -> ImageData:
             url = self._get_asset_url(asset)
+            idx = asset_indexes.get(asset) or kwargs.pop("indexes", None)  # type: ignore
             with self.reader(url, tms=self.tms, **self.reader_options) as cog:  # type: ignore
-                data = cog.tile(
-                    *args,
-                    indexes=asset_indexes.get(asset, kwargs.pop("indexes", None)),  # type: ignore
-                    expression=asset_expression.get(asset),  # type: ignore
-                    **kwargs,
-                )
+                data = cog.tile(*args, indexes=idx, **kwargs)
                 data.band_names = [f"{asset}_{n}" for n in data.band_names]
                 return data
 
-        output = multi_arrays(
-            assets,
-            _reader,
-            tile_x,
-            tile_y,
-            tile_z,
-            **kwargs,
-        )
-
+        img = multi_arrays(assets, _reader, tile_x, tile_y, tile_z, **kwargs)
         if expression:
-            blocks = get_expression_blocks(expression)
-            output.data = apply_expression(blocks, assets, output.data)
-            output.band_names = blocks
+            return img.apply_expression(expression)
 
-        return output
+        return img
 
     def part(
         self,
@@ -501,13 +500,16 @@ def part(
             assets (sequence of str or str, optional): assets to fetch info from.
             expression (str, optional): rio-tiler expression for the asset list (e.g. asset1/asset2+asset3).
             asset_indexes (dict, optional): Band indexes for each asset (e.g {"asset1": 1, "asset2": (1, 2,)}).
-            asset_expression (dict, optional): rio-tiler expression for each asset (e.g. {"asset1": "b1/b2+b3", "asset2": ...}).
+            asset_expression (dict, optional): rio-tiler expression for each asset (e.g. {"asset1": "b1/b2+b3", "asset2": ...}).  **Deprecated**
             kwargs (optional): Options to forward to the `self.reader.part` method.
 
         Returns:
             rio_tiler.models.ImageData: ImageData instance with data, mask and tile spatial info.
 
         """
+        if asset_expression:
+            _AssetExpressionWarning()
+
         if isinstance(assets, str):
             assets = (assets,)
 
@@ -526,28 +528,20 @@ def part(
             )
 
         asset_indexes = asset_indexes or {}
-        asset_expression = asset_expression or {}
 
         def _reader(asset: str, *args: Any, **kwargs: Any) -> ImageData:
             url = self._get_asset_url(asset)
+            idx = asset_indexes.get(asset) or kwargs.pop("indexes", None)  # type: ignore
             with self.reader(url, tms=self.tms, **self.reader_options) as cog:  # type: ignore
-                data = cog.part(
-                    *args,
-                    indexes=asset_indexes.get(asset, kwargs.pop("indexes", None)),  # type: ignore
-                    expression=asset_expression.get(asset),  # type: ignore
-                    **kwargs,
-                )
+                data = cog.part(*args, indexes=idx, **kwargs)
                 data.band_names = [f"{asset}_{n}" for n in data.band_names]
                 return data
 
-        output = multi_arrays(assets, _reader, bbox, **kwargs)
-
+        img = multi_arrays(assets, _reader, bbox, **kwargs)
         if expression:
-            blocks = get_expression_blocks(expression)
-            output.data = apply_expression(blocks, assets, output.data)
-            output.band_names = blocks
+            return img.apply_expression(expression)
 
-        return output
+        return img
 
     def preview(
         self,
@@ -563,13 +557,16 @@ def preview(
             assets (sequence of str or str, optional): assets to fetch info from.
             expression (str, optional): rio-tiler expression for the asset list (e.g. asset1/asset2+asset3).
             asset_indexes (dict, optional): Band indexes for each asset (e.g {"asset1": 1, "asset2": (1, 2,)}).
-            asset_expression (dict, optional): rio-tiler expression for each asset (e.g. {"asset1": "b1/b2+b3", "asset2": ...}).
+            asset_expression (dict, optional): rio-tiler expression for each asset (e.g. {"asset1": "b1/b2+b3", "asset2": ...}). **Deprecated**
             kwargs (optional): Options to forward to the `self.reader.preview` method.
 
         Returns:
             rio_tiler.models.ImageData: ImageData instance with data, mask and tile spatial info.
 
         """
+        if asset_expression:
+            _AssetExpressionWarning()
+
         if isinstance(assets, str):
             assets = (assets,)
 
@@ -588,27 +585,20 @@ def preview(
             )
 
         asset_indexes = asset_indexes or {}
-        asset_expression = asset_expression or {}
 
         def _reader(asset: str, **kwargs: Any) -> ImageData:
             url = self._get_asset_url(asset)
+            idx = asset_indexes.get(asset) or kwargs.pop("indexes", None)  # type: ignore
             with self.reader(url, tms=self.tms, **self.reader_options) as cog:  # type: ignore
-                data = cog.preview(
-                    indexes=asset_indexes.get(asset, kwargs.pop("indexes", None)),  # type: ignore
-                    expression=asset_expression.get(asset),  # type: ignore
-                    **kwargs,
-                )
+                data = cog.preview(indexes=idx, **kwargs)
                 data.band_names = [f"{asset}_{n}" for n in data.band_names]
                 return data
 
-        output = multi_arrays(assets, _reader, **kwargs)
-
+        img = multi_arrays(assets, _reader, **kwargs)
         if expression:
-            blocks = get_expression_blocks(expression)
-            output.data = apply_expression(blocks, assets, output.data)
-            output.band_names = blocks
+            return img.apply_expression(expression)
 
-        return output
+        return img
 
     def point(
         self,
@@ -619,7 +609,7 @@ def point(
         asset_indexes: Optional[Dict[str, Indexes]] = None,  # Indexes for each asset
         asset_expression: Optional[Dict[str, str]] = None,  # Expression for each asset
         **kwargs: Any,
-    ) -> List:
+    ) -> Tuple[List, List[str]]:
         """Read pixel value from multiple assets.
 
         Args:
@@ -628,13 +618,17 @@ def point(
             assets (sequence of str or str, optional): assets to fetch info from.
             expression (str, optional): rio-tiler expression for the asset list (e.g. asset1/asset2+asset3).
             asset_indexes (dict, optional): Band indexes for each asset (e.g {"asset1": 1, "asset2": (1, 2,)}).
-            asset_expression (dict, optional): rio-tiler expression for each asset (e.g. {"asset1": "b1/b2+b3", "asset2": ...}).
+            asset_expression (dict, optional): rio-tiler expression for each asset (e.g. {"asset1": "b1/b2+b3", "asset2": ...}). **Deprecated**
             kwargs (optional): Options to forward to the `self.reader.point` method.
 
         Returns:
             list: Pixel values per assets.
+            list: band names
 
         """
+        if asset_expression:
+            _AssetExpressionWarning()
+
         if isinstance(assets, str):
             assets = (assets,)
 
@@ -653,26 +647,34 @@ def point(
             )
 
         asset_indexes = asset_indexes or {}
-        asset_expression = asset_expression or {}
 
-        def _reader(asset: str, *args, **kwargs: Any) -> Dict:
+        def _reader(asset: str, *args, **kwargs: Any) -> Tuple[List, List[str]]:
             url = self._get_asset_url(asset)
+            idx = asset_indexes.get(asset) or kwargs.pop("indexes", None)  # type: ignore
             with self.reader(url, tms=self.tms, **self.reader_options) as cog:  # type: ignore
-                return cog.point(
-                    *args,
-                    indexes=asset_indexes.get(asset, kwargs.pop("indexes", None)),  # type: ignore
-                    expression=asset_expression.get(asset),  # type: ignore
-                    **kwargs,
+                points, band_names = cog.point(*args, indexes=idx, **kwargs)
+                return (
+                    points,
+                    [f"{asset}_{n}" for n in band_names],
                 )
 
         data = multi_values(assets, _reader, lon, lat, **kwargs)
+        values, band_names = zip(*[(v, b) for _, (v, b) in data.items()])
+
+        # Create an array with all the point values
+        values = numpy.array(list(itertools.chain.from_iterable(values)))
+
+        # Create a list of all point's band name
+        band_names = list(itertools.chain.from_iterable(band_names))
 
-        values = [numpy.array(d) for _, d in data.items()]
         if expression:
             blocks = get_expression_blocks(expression)
-            values = apply_expression(blocks, assets, values)
+            return (
+                apply_expression(blocks, band_names, values).tolist(),
+                blocks,
+            )
 
-        return [v.tolist() for v in values]
+        return values.tolist(), band_names
 
     def feature(
         self,
@@ -690,13 +692,16 @@ def feature(
             assets (sequence of str or str, optional): assets to fetch info from.
             expression (str, optional): rio-tiler expression for the asset list (e.g. asset1/asset2+asset3).
             asset_indexes (dict, optional): Band indexes for each asset (e.g {"asset1": 1, "asset2": (1, 2,)}).
-            asset_expression (dict, optional): rio-tiler expression for each asset (e.g. {"asset1": "b1/b2+b3", "asset2": ...}).
+            asset_expression (dict, optional): rio-tiler expression for each asset (e.g. {"asset1": "b1/b2+b3", "asset2": ...}). **Deprecated**
             kwargs (optional): Options to forward to the `self.reader.feature` method.
 
         Returns:
             rio_tiler.models.ImageData: ImageData instance with data, mask and tile spatial info.
 
         """
+        if asset_expression:
+            _AssetExpressionWarning()
+
         if isinstance(assets, str):
             assets = (assets,)
 
@@ -715,28 +720,20 @@ def feature(
             )
 
         asset_indexes = asset_indexes or {}
-        asset_expression = asset_expression or {}
 
         def _reader(asset: str, *args: Any, **kwargs: Any) -> ImageData:
             url = self._get_asset_url(asset)
+            idx = asset_indexes.get(asset) or kwargs.pop("indexes", None)  # type: ignore
             with self.reader(url, tms=self.tms, **self.reader_options) as cog:  # type: ignore
-                data = cog.feature(
-                    *args,
-                    indexes=asset_indexes.get(asset, kwargs.pop("indexes", None)),  # type: ignore
-                    expression=asset_expression.get(asset),  # type: ignore
-                    **kwargs,
-                )
+                data = cog.feature(*args, indexes=idx, **kwargs)
                 data.band_names = [f"{asset}_{n}" for n in data.band_names]
                 return data
 
-        output = multi_arrays(assets, _reader, shape, **kwargs)
-
+        img = multi_arrays(assets, _reader, shape, **kwargs)
         if expression:
-            blocks = get_expression_blocks(expression)
-            output.data = apply_expression(blocks, assets, output.data)
-            output.band_names = blocks
+            return img.apply_expression(expression)
 
-        return output
+        return img
 
 
 @attr.s
@@ -942,17 +939,15 @@ def _reader(band: str, *args: Any, **kwargs: Any) -> ImageData:
             url = self._get_band_url(band)
             with self.reader(url, tms=self.tms, **self.reader_options) as cog:  # type: ignore
                 data = cog.tile(*args, **kwargs)
-                data.band_names = [band]
+                data.band_names = [band]  # use `band` as name instead of band index
                 return data
 
-        output = multi_arrays(bands, _reader, tile_x, tile_y, tile_z, **kwargs)
+        img = multi_arrays(bands, _reader, tile_x, tile_y, tile_z, **kwargs)
 
         if expression:
-            blocks = get_expression_blocks(expression)
-            output.data = apply_expression(blocks, bands, output.data)
-            output.band_names = blocks
+            return img.apply_expression(expression)
 
-        return output
+        return img
 
     def part(
         self,
@@ -994,17 +989,15 @@ def _reader(band: str, *args: Any, **kwargs: Any) -> ImageData:
             url = self._get_band_url(band)
             with self.reader(url, tms=self.tms, **self.reader_options) as cog:  # type: ignore
                 data = cog.part(*args, **kwargs)
-                data.band_names = [band]
+                data.band_names = [band]  # use `band` as name instead of band index
                 return data
 
-        output = multi_arrays(bands, _reader, bbox, **kwargs)
+        img = multi_arrays(bands, _reader, bbox, **kwargs)
 
         if expression:
-            blocks = get_expression_blocks(expression)
-            output.data = apply_expression(blocks, bands, output.data)
-            output.band_names = blocks
+            return img.apply_expression(expression)
 
-        return output
+        return img
 
     def preview(
         self,
@@ -1044,17 +1037,15 @@ def _reader(band: str, **kwargs: Any) -> ImageData:
             url = self._get_band_url(band)
             with self.reader(url, tms=self.tms, **self.reader_options) as cog:  # type: ignore
                 data = cog.preview(**kwargs)
-                data.band_names = [band]
+                data.band_names = [band]  # use `band` as name instead of band index
                 return data
 
-        output = multi_arrays(bands, _reader, **kwargs)
+        img = multi_arrays(bands, _reader, **kwargs)
 
         if expression:
-            blocks = get_expression_blocks(expression)
-            output.data = apply_expression(blocks, bands, output.data)
-            output.band_names = blocks
+            return img.apply_expression(expression)
 
-        return output
+        return img
 
     def point(
         self,
@@ -1063,7 +1054,7 @@ def point(
         bands: Union[Sequence[str], str] = None,
         expression: Optional[str] = None,
         **kwargs: Any,
-    ) -> List:
+    ) -> Tuple[List, List[str]]:
         """Read a pixel values from multiple bands.
 
         Args:
@@ -1075,6 +1066,7 @@ def point(
 
         Returns:
             list: Pixel value per bands.
+            list: band names
 
         """
         if isinstance(bands, str):
@@ -1094,19 +1086,27 @@ def point(
                 "bands must be passed either via expression or bands options."
             )
 
-        def _reader(band: str, *args, **kwargs: Any) -> Dict:
+        def _reader(band: str, *args, **kwargs: Any) -> Tuple[numpy.array, List[str]]:
             url = self._get_band_url(band)
             with self.reader(url, tms=self.tms, **self.reader_options) as cog:  # type: ignore
-                return cog.point(*args, **kwargs)[0]  # We only return the first value
+                values, _ = cog.point(*args, **kwargs)
+                return (
+                    values[0],  # We only return the first value
+                    [band],
+                )
 
         data = multi_values(bands, _reader, lon, lat, **kwargs)
+        values, band_names = zip(*[(v, b) for b, (v, _) in data.items()])
+        values = numpy.array(values)
 
-        values = [numpy.array(d) for _, d in data.items()]
         if expression:
             blocks = get_expression_blocks(expression)
-            values = apply_expression(blocks, bands, values)
+            return (
+                apply_expression(blocks, band_names, values).tolist(),
+                blocks,
+            )
 
-        return [v.tolist() for v in values]
+        return [v.tolist() for v in values], list(band_names)
 
     def feature(
         self,
@@ -1148,14 +1148,12 @@ def _reader(band: str, *args: Any, **kwargs: Any) -> ImageData:
             url = self._get_band_url(band)
             with self.reader(url, tms=self.tms, **self.reader_options) as cog:  # type: ignore
                 data = cog.feature(*args, **kwargs)
-                data.band_names = [band]
+                data.band_names = [band]  # use `band` as name instead of band index
                 return data
 
-        output = multi_arrays(bands, _reader, shape, **kwargs)
+        img = multi_arrays(bands, _reader, shape, **kwargs)
 
         if expression:
-            blocks = get_expression_blocks(expression)
-            output.data = apply_expression(blocks, bands, output.data)
-            output.band_names = blocks
+            return img.apply_expression(expression)
 
-        return output
+        return img
diff --git a/rio_tiler/io/cogeo.py b/rio_tiler/io/cogeo.py
index 8674d315..6dbedbfb 100644
--- a/rio_tiler/io/cogeo.py
+++ b/rio_tiler/io/cogeo.py
@@ -2,7 +2,7 @@
 
 import contextlib
 import warnings
-from typing import Any, Callable, Dict, List, Optional, Sequence, Union
+from typing import Any, Callable, Dict, List, Optional, Sequence, Tuple, Union
 
 import attr
 import numpy
@@ -23,13 +23,7 @@
 from ..expression import apply_expression, get_expression_blocks, parse_expression
 from ..models import BandStatistics, ImageData, Info
 from ..types import BBox, DataMaskType, Indexes, NoData, NumType
-from ..utils import (
-    create_cutline,
-    get_array_statistics,
-    get_bands_names,
-    has_alpha_band,
-    has_mask_band,
-)
+from ..utils import create_cutline, get_array_statistics, has_alpha_band, has_mask_band
 from .base import BaseReader
 
 
@@ -259,10 +253,10 @@ def _get_descr(ix):
             "minzoom": self.minzoom,
             "maxzoom": self.maxzoom,
             "band_metadata": [
-                (f"{ix}", self.dataset.tags(ix)) for ix in self.dataset.indexes
+                (f"b{ix}", self.dataset.tags(ix)) for ix in self.dataset.indexes
             ],
             "band_descriptions": [
-                (f"{ix}", _get_descr(ix)) for ix in self.dataset.indexes
+                (f"b{ix}", _get_descr(ix)) for ix in self.dataset.indexes
             ],
             "dtype": self.dataset.meta["dtype"],
             "colorinterp": [
@@ -446,22 +440,12 @@ def part(
             buffer=buffer,
             **kwargs,
         )
+        img.assets = [self.input]
 
-        if expression and indexes:
-            blocks = get_expression_blocks(expression)
-            bands = [f"b{bidx}" for bidx in indexes]
-            img.data = apply_expression(blocks, bands, img.data)
-
-        return ImageData(
-            img.data,
-            img.mask,
-            bounds=img.bounds,
-            crs=img.crs,
-            assets=[self.input],
-            band_names=get_bands_names(
-                indexes=indexes, expression=expression, count=img.count
-            ),
-        )
+        if expression:
+            return img.apply_expression(expression)
+
+        return img
 
     def preview(
         self,
@@ -508,22 +492,12 @@ def preview(
             height=height,
             **kwargs,
         )
+        img.assets = [self.input]
 
-        if expression and indexes:
-            blocks = get_expression_blocks(expression)
-            bands = [f"b{bidx}" for bidx in indexes]
-            img.data = apply_expression(blocks, bands, img.data)
-
-        return ImageData(
-            img.data,
-            img.mask,
-            bounds=img.bounds,
-            crs=img.crs,
-            assets=[self.input],
-            band_names=get_bands_names(
-                indexes=indexes, expression=expression, count=img.count
-            ),
-        )
+        if expression:
+            return img.apply_expression(expression)
+
+        return img
 
     def point(
         self,
@@ -533,7 +507,7 @@ def point(
         indexes: Optional[Indexes] = None,
         expression: Optional[str] = None,
         **kwargs: Any,
-    ) -> List:
+    ) -> Tuple[List, List[str]]:
         """Read a pixel value from a COG.
 
         Args:
@@ -546,6 +520,7 @@ def point(
 
         Returns:
             list: Pixel value per band indexes.
+            list: band names
 
         """
         kwargs = {**self._kwargs, **kwargs}
@@ -562,16 +537,18 @@ def point(
         if expression:
             indexes = parse_expression(expression)
 
-        point = reader.point(
+        values, band_names = reader.point(
             self.dataset, (lon, lat), indexes=indexes, coord_crs=coord_crs, **kwargs
         )
 
-        if expression and indexes:
+        if expression:
             blocks = get_expression_blocks(expression)
-            bands = [f"b{bidx}" for bidx in indexes]
-            point = apply_expression(blocks, bands, numpy.array(point)).tolist()
+            return (
+                apply_expression(blocks, band_names, numpy.array(values)).tolist(),
+                blocks,
+            )
 
-        return point
+        return values, band_names
 
     def feature(
         self,
@@ -660,22 +637,12 @@ def read(
             indexes = parse_expression(expression)
 
         img = reader.read(self.dataset, indexes=indexes, **kwargs)
+        img.assets = [self.input]
 
-        if expression and indexes:
-            blocks = get_expression_blocks(expression)
-            bands = [f"b{bidx}" for bidx in indexes]
-            img.data = apply_expression(blocks, bands, img.data)
-
-        return ImageData(
-            img.data,
-            img.mask,
-            bounds=img.bounds,
-            crs=img.crs,
-            assets=[self.input],
-            band_names=get_bands_names(
-                indexes=indexes, expression=expression, count=img.count
-            ),
-        )
+        if expression:
+            return img.apply_expression(expression)
+
+        return img
 
 
 @attr.s
diff --git a/rio_tiler/models.py b/rio_tiler/models.py
index e67dee30..e9f28d63 100644
--- a/rio_tiler/models.py
+++ b/rio_tiler/models.py
@@ -18,6 +18,7 @@
 from rio_color.utils import scale_dtype, to_math_type
 
 from .errors import InvalidDatatypeWarning
+from .expression import apply_expression, get_expression_blocks
 from .types import ColorMapType, GDALColorMapType, IntervalTuple, NumType
 from .utils import linear_rescale, render, resize_array
 
@@ -163,7 +164,7 @@ def _validate_data(self, attribute, value):
 
     @band_names.default
     def _default_names(self):
-        return [f"{ix + 1}" for ix in range(self.count)]
+        return [f"b{ix + 1}" for ix in range(self.count)]
 
     @mask.default
     def _default_mask(self):
@@ -274,6 +275,19 @@ def apply_color_formula(self, color_formula: Optional[str]):
         for ops in parse_operations(color_formula):
             self.data = scale_dtype(ops(to_math_type(self.data)), numpy.uint8)
 
+    def apply_expression(self, expression: str) -> "ImageData":
+        """Apply expression to the image data."""
+        blocks = get_expression_blocks(expression)
+        return ImageData(
+            apply_expression(blocks, self.band_names, self.data),
+            self.mask,
+            assets=self.assets,
+            crs=self.crs,
+            bounds=self.bounds,
+            band_names=blocks,
+            metadata=self.metadata,
+        )
+
     def post_process(
         self,
         in_range: Optional[Sequence[IntervalTuple]] = None,
diff --git a/rio_tiler/reader.py b/rio_tiler/reader.py
index ee14d452..4193849b 100644
--- a/rio_tiler/reader.py
+++ b/rio_tiler/reader.py
@@ -71,6 +71,8 @@ def read(
                 "Alpha band was removed from the output data array", AlphaBandWarning
             )
 
+    band_names = [f"b{idx}" for idx in indexes]
+
     vrt_params = dict(add_alpha=True, resampling=Resampling[resampling_method])
     nodata = nodata if nodata is not None else src_dst.nodata
     if nodata is not None:
@@ -127,6 +129,7 @@ def read(
             mask,
             bounds=windows.bounds(window, vrt.transform),
             crs=vrt.crs,
+            band_names=band_names,
         )
 
 
@@ -324,7 +327,7 @@ def point(
     post_process: Optional[
         Callable[[numpy.ndarray, numpy.ndarray], DataMaskType]
     ] = None,
-) -> List:
+) -> Tuple[List, List[str]]:
     """Read a pixel value for a point.
 
     Args:
@@ -341,6 +344,7 @@ def point(
 
     Returns:
         list: Pixel value per band indexes.
+        list: band names
 
     """
     if isinstance(indexes, int):
@@ -357,6 +361,8 @@ def point(
 
     indexes = indexes if indexes is not None else src_dst.indexes
 
+    band_names = [f"b{idx}" for idx in indexes]
+
     vrt_params = dict(add_alpha=True, resampling=Resampling[resampling_method])
     nodata = nodata if nodata is not None else src_dst.nodata
     if nodata is not None:
@@ -383,4 +389,4 @@ def point(
     if post_process:
         point_values, _ = post_process(point_values, mask)
 
-    return point_values.tolist()
+    return point_values.tolist(), band_names
diff --git a/rio_tiler/utils.py b/rio_tiler/utils.py
index 0a4ae31b..56f9bc99 100644
--- a/rio_tiler/utils.py
+++ b/rio_tiler/utils.py
@@ -67,6 +67,10 @@ def get_bands_names(
     count: Optional[int] = None,
 ) -> List[str]:
     """Define bands names based on expression, indexes or band count."""
+    warnings.warn(
+        "`get_bands_names` is deprecated, and will be removed in rio-tiler 4.0`.",
+        DeprecationWarning,
+    )
     if expression:
         return get_expression_blocks(expression)
 
diff --git a/tests/fixtures/scene_b1.tif b/tests/fixtures/scene_band1.tif
similarity index 100%
rename from tests/fixtures/scene_b1.tif
rename to tests/fixtures/scene_band1.tif
diff --git a/tests/fixtures/scene_b2.tif b/tests/fixtures/scene_band2.tif
similarity index 100%
rename from tests/fixtures/scene_b2.tif
rename to tests/fixtures/scene_band2.tif
diff --git a/tests/test_io_MultiBand.py b/tests/test_io_MultiBand.py
index 749562ca..43be36e0 100644
--- a/tests/test_io_MultiBand.py
+++ b/tests/test_io_MultiBand.py
@@ -59,86 +59,98 @@ def _get_band_url(self, band: str) -> str:
 def test_MultiBandReader():
     """Should work as expected."""
     with BandFileReader(PREFIX) as cog:
-        assert cog.bands == ["b1", "b2"]
+        assert cog.bands == ["band1", "band2"]
         assert cog.minzoom is not None
         assert cog.maxzoom is not None
         assert cog.bounds
         assert cog.bounds
         assert cog.crs
 
-        assert sorted(cog.parse_expression("b1/b2")) == ["b1", "b2"]
+        assert sorted(cog.parse_expression("band1/band2")) == ["band1", "band2"]
 
         with pytest.warns(UserWarning):
             meta = cog.info()
-        assert meta.band_descriptions == [("b1", ""), ("b2", "")]
+        assert meta.band_descriptions == [("band1", ""), ("band2", "")]
 
-        meta = cog.info(bands="b1")
-        assert meta.band_descriptions == [("b1", "")]
+        meta = cog.info(bands="band1")
+        assert meta.band_descriptions == [("band1", "")]
 
-        meta = cog.info(bands=("b1", "b2"))
-        assert meta.band_descriptions == [("b1", ""), ("b2", "")]
+        meta = cog.info(bands=("band1", "band2"))
+        assert meta.band_descriptions == [("band1", ""), ("band2", "")]
 
         with pytest.warns(UserWarning):
             stats = cog.statistics()
-            assert stats["b1"]
-            assert stats["b2"]
+            assert stats["band1"]
+            assert stats["band2"]
 
-        stats = cog.statistics(bands="b1")
-        assert "b1" in stats
-        assert isinstance(stats["b1"], BandStatistics)
+        stats = cog.statistics(bands="band1")
+        assert "band1" in stats
+        assert isinstance(stats["band1"], BandStatistics)
 
-        stats = cog.statistics(bands=("b1", "b2"))
-        assert stats["b1"]
-        assert stats["b2"]
+        stats = cog.statistics(bands=("band1", "band2"))
+        assert stats["band1"]
+        assert stats["band2"]
 
-        stats = cog.statistics(expression="b1;b1+b2;b1-100")
-        assert stats["b1"]
-        assert stats["b1+b2"]
-        assert stats["b1-100"]
+        stats = cog.statistics(expression="band1;band1+band2;band1-100")
+        assert stats["band1"]
+        assert stats["band1+band2"]
+        assert stats["band1-100"]
 
         with pytest.raises(MissingBands):
             cog.tile(238, 218, 9)
 
-        tile = cog.tile(238, 218, 9, bands="b1")
+        tile = cog.tile(238, 218, 9, bands="band1")
         assert tile.data.shape == (1, 256, 256)
-        assert tile.band_names == ["b1"]
+        assert tile.band_names == ["band1"]
 
         with pytest.warns(ExpressionMixingWarning):
-            tile = cog.tile(238, 218, 9, bands="b1", expression="b1*2")
+            tile = cog.tile(238, 218, 9, bands="band1", expression="band1*2")
         assert tile.data.shape == (1, 256, 256)
-        assert tile.band_names == ["b1*2"]
+        assert tile.band_names == ["band1*2"]
 
         with pytest.raises(MissingBands):
             cog.part((-11.5, 24.5, -11.0, 25.0))
 
-        tile = cog.part((-11.5, 24.5, -11.0, 25.0), bands="b1")
+        tile = cog.part((-11.5, 24.5, -11.0, 25.0), bands="band1")
         assert tile.data.any()
-        assert tile.band_names == ["b1"]
+        assert tile.band_names == ["band1"]
 
         with pytest.warns(ExpressionMixingWarning):
-            tile = cog.part((-11.5, 24.5, -11.0, 25.0), bands="b1", expression="b1*2")
+            tile = cog.part(
+                (-11.5, 24.5, -11.0, 25.0), bands="band1", expression="band1*2"
+            )
         assert tile.data.any()
-        assert tile.band_names == ["b1*2"]
+        assert tile.band_names == ["band1*2"]
 
         with pytest.raises(MissingBands):
             cog.preview()
 
-        tile = cog.preview(bands="b1")
+        tile = cog.preview(bands="band1")
         assert tile.data.any()
-        assert tile.band_names == ["b1"]
+        assert tile.band_names == ["band1"]
 
         with pytest.warns(ExpressionMixingWarning):
-            tile = cog.preview(bands="b1", expression="b1*2")
+            tile = cog.preview(bands="band1", expression="band1*2")
         assert tile.data.any()
-        assert tile.band_names == ["b1*2"]
+        assert tile.band_names == ["band1*2"]
 
         with pytest.raises(MissingBands):
             cog.point(-11.5, 24.5)
 
-        assert cog.point(-11.5, 24.5, bands="b1")
+        pts, names = cog.point(-11.5, 24.5, bands="band1")
+        assert len(pts) == 1
+        assert names == ["band1"]
+
+        pts, names = cog.point(-11.5, 24.5, bands=("band1", "band2"))
+        assert len(pts) == 2
+        assert names == ["band1", "band2"]
+
+        pts, names = cog.point(-11.5, 24.5, expression="band1/band2")
+        assert len(pts) == 1
+        assert names == ["band1/band2"]
 
         with pytest.warns(ExpressionMixingWarning):
-            assert cog.point(-11.5, 24.5, bands="b1", expression="b1*2")
+            assert cog.point(-11.5, 24.5, bands="band1", expression="band1*2")
 
         feat = {
             "type": "Feature",
@@ -162,11 +174,11 @@ def test_MultiBandReader():
         with pytest.raises(MissingBands):
             cog.feature(feat)
 
-        img = cog.feature(feat, bands="b1")
+        img = cog.feature(feat, bands="band1")
         assert img.data.any()
         assert not img.mask.all()
-        assert img.band_names == ["b1"]
+        assert img.band_names == ["band1"]
 
         with pytest.warns(ExpressionMixingWarning):
-            img = cog.feature(feat, bands="b1", expression="b1*2")
-            assert img.band_names == ["b1*2"]
+            img = cog.feature(feat, bands="band1", expression="band1*2")
+            assert img.band_names == ["band1*2"]
diff --git a/tests/test_io_cogeo.py b/tests/test_io_cogeo.py
index 8bf1868f..303984e6 100644
--- a/tests/test_io_cogeo.py
+++ b/tests/test_io_cogeo.py
@@ -115,7 +115,7 @@ def test_info_valid():
         assert meta.offset
         assert meta.band_metadata
         band_meta = meta.band_metadata[0]
-        assert band_meta[0] == "1"
+        assert band_meta[0] == "b1"
         assert "STATISTICS_MAXIMUM" in band_meta[1]
 
     with COGReader(COG_ALPHA) as cog:
@@ -142,7 +142,7 @@ def test_tile_valid_default():
         img = cog.tile(43, 24, 7)
         assert img.data.shape == (1, 256, 256)
         assert img.mask.all()
-        assert img.band_names == ["1"]
+        assert img.band_names == ["b1"]
 
         # Validate that Tile and Part gives the same result
         tile_bounds = WEB_MERCATOR_TMS.xy_bounds(43, 24, 7)
@@ -183,7 +183,7 @@ def test_tile_valid_default():
             ),
         )
         assert img.data.shape == (2, 256, 256)
-        assert img.band_names == ["1", "1"]
+        assert img.band_names == ["b1", "b1"]
 
     # We are using a file that is aligned with the grid so no resampling should be involved
     with COGReader(COG_WEB) as cog:
@@ -232,20 +232,24 @@ def test_point_valid():
     lon = -56.624124590533825
     lat = 73.52687881825946
     with COGReader(COG_NODATA) as cog:
-        pts = cog.point(lon, lat)
+        pts, names = cog.point(lon, lat)
         assert len(pts) == 1
+        assert names == ["b1"]
 
-        pts = cog.point(lon, lat, expression="b1*2;b1-100")
+        pts, names = cog.point(lon, lat, expression="b1*2;b1-100")
         assert len(pts) == 2
+        assert names == ["b1*2", "b1-100"]
 
         with pytest.warns(ExpressionMixingWarning):
-            pts = cog.point(lon, lat, indexes=(1, 2, 3), expression="b1*2")
+            pts, names = cog.point(lon, lat, indexes=(1, 2, 3), expression="b1*2")
             assert len(pts) == 1
+            assert names == ["b1*2"]
 
-        pts = cog.point(lon, lat, indexes=1)
+        pts, names = cog.point(lon, lat, indexes=1)
         assert len(pts) == 1
+        assert names == ["b1"]
 
-        pts = cog.point(
+        pts, names = cog.point(
             lon,
             lat,
             indexes=(
@@ -254,6 +258,7 @@ def test_point_valid():
             ),
         )
         assert len(pts) == 2
+        assert names == ["b1", "b1"]
 
 
 def test_area_valid():
@@ -267,7 +272,7 @@ def test_area_valid():
     with COGReader(COG_NODATA) as cog:
         img = cog.part(bbox)
         assert img.data.shape == (1, 11, 40)
-        assert img.band_names == ["1"]
+        assert img.band_names == ["b1"]
 
         data, mask = cog.part(bbox, dst_crs=cog.dataset.crs)
         assert data.shape == (1, 28, 30)
@@ -295,7 +300,7 @@ def test_area_valid():
             ),
         )
         assert img.data.shape == (2, 11, 40)
-        assert img.band_names == ["1", "1"]
+        assert img.band_names == ["b1", "b1"]
 
 
 def test_preview_valid():
@@ -303,7 +308,7 @@ def test_preview_valid():
     with COGReader(COGEO) as cog:
         img = cog.preview(max_size=128)
         assert img.data.shape == (1, 128, 128)
-        assert img.band_names == ["1"]
+        assert img.band_names == ["b1"]
 
         data, mask = cog.preview()
         assert data.shape == (1, 1024, 1021)
@@ -328,7 +333,7 @@ def test_preview_valid():
             ),
         )
         assert img.data.shape == (2, 128, 128)
-        assert img.band_names == ["1", "1"]
+        assert img.band_names == ["b1", "b1"]
 
 
 def test_statistics():
@@ -336,21 +341,21 @@ def test_statistics():
     with COGReader(COGEO) as cog:
         stats = cog.statistics()
         assert len(stats) == 1
-        assert isinstance(stats["1"], BandStatistics)
-        assert stats["1"].percentile_2
-        assert stats["1"].percentile_98
+        assert isinstance(stats["b1"], BandStatistics)
+        assert stats["b1"].percentile_2
+        assert stats["b1"].percentile_98
 
     with COGReader(COGEO) as cog:
         stats = cog.statistics(percentiles=[3])
-        assert stats["1"].percentile_3
+        assert stats["b1"].percentile_3
 
     with COGReader(COGEO) as cog:
         stats = cog.statistics(percentiles=[3])
-        assert stats["1"].percentile_3
+        assert stats["b1"].percentile_3
 
     with COGReader(COG_CMAP) as cog:
         stats = cog.statistics(categorical=True)
-        assert stats["1"].histogram[1] == [
+        assert stats["b1"].histogram[1] == [
             1.0,
             3.0,
             4.0,
@@ -365,7 +370,7 @@ def test_statistics():
         ]
 
         stats = cog.statistics(categorical=True, categories=[1, 3])
-        assert stats["1"].histogram[1] == [
+        assert stats["b1"].histogram[1] == [
             1.0,
             3.0,
         ]
@@ -373,7 +378,7 @@ def test_statistics():
     # make sure kwargs are passed to `preview`
     with COGReader(COGEO) as cog:
         stats = cog.statistics(width=100, height=100, max_size=None)
-        assert stats["1"].count == 10000.0
+        assert stats["b1"].count == 10000.0
 
     # Check results for expression
     with COGReader(COGEO) as cog:
@@ -406,11 +411,11 @@ def test_COGReader_Options():
         assert not numpy.array_equal(data_default, data)
 
     with COGReader(COG_SCALE, unscale=True) as cog:
-        p = cog.point(310000, 4100000, coord_crs=cog.dataset.crs)
+        p, _ = cog.point(310000, 4100000, coord_crs=cog.dataset.crs)
         assert round(p[0], 3) == 1000.892
 
         # passing unscale in method should overwrite the defaults
-        p = cog.point(310000, 4100000, coord_crs=cog.dataset.crs, unscale=False)
+        p, _ = cog.point(310000, 4100000, coord_crs=cog.dataset.crs, unscale=False)
         assert p[0] == 8917
 
     cutline = "POLYGON ((13 1685, 1010 6, 2650 967, 1630 2655, 13 1685))"
@@ -432,10 +437,10 @@ def callback(data, mask):
     lon = -56.624124590533825
     lat = 73.52687881825946
     with COGReader(COG_NODATA, post_process=callback) as cog:
-        pts = cog.point(lon, lat)
+        pts, _ = cog.point(lon, lat)
 
     with COGReader(COG_NODATA) as cog:
-        pts_init = cog.point(lon, lat)
+        pts_init, _ = cog.point(lon, lat)
     assert pts[0] == pts_init[0] * 2
 
 
@@ -451,7 +456,7 @@ def test_cog_with_internal_gcps():
 
         metadata = cog.info()
         assert len(metadata.band_metadata) == 1
-        assert metadata.band_descriptions == [("1", "")]
+        assert metadata.band_descriptions == [("b1", "")]
 
         tile_z = 8
         tile_x = 183
@@ -479,7 +484,7 @@ def test_cog_with_internal_gcps():
 
                 metadata = cog.info()
                 assert len(metadata.band_metadata) == 1
-                assert metadata.band_descriptions == [("1", "")]
+                assert metadata.band_descriptions == [("b1", "")]
 
                 tile_z = 8
                 tile_x = 183
@@ -611,7 +616,7 @@ def test_feature_valid():
     with COGReader(COG_NODATA) as cog:
         img = cog.feature(feature, max_size=1024)
         assert img.data.shape == (1, 348, 1024)
-        assert img.band_names == ["1"]
+        assert img.band_names == ["b1"]
 
         img = cog.feature(feature, dst_crs=cog.dataset.crs, max_size=1024)
         assert img.data.shape == (1, 1024, 869)
@@ -642,7 +647,7 @@ def test_feature_valid():
             max_size=1024,
         )
         assert img.data.shape == (2, 348, 1024)
-        assert img.band_names == ["1", "1"]
+        assert img.band_names == ["b1", "b1"]
 
         # feature overlaping on mask area
         mask_feat = {
diff --git a/tests/test_io_stac.py b/tests/test_io_stac.py
index daf93a52..8f098b35 100644
--- a/tests/test_io_stac.py
+++ b/tests/test_io_stac.py
@@ -138,22 +138,29 @@ def test_tile_valid(rio):
         img = stac.tile(71, 102, 8, assets="green")
         assert img.data.shape == (1, 256, 256)
         assert img.mask.shape == (256, 256)
-        assert img.band_names == ["green_1"]
+        assert img.band_names == ["green_b1"]
 
-        data, mask = stac.tile(71, 102, 8, assets=("green",))
-        assert data.shape == (1, 256, 256)
-        assert mask.shape == (256, 256)
+        img = stac.tile(71, 102, 8, assets=("green",))
+        assert img.data.shape == (1, 256, 256)
+        assert img.mask.shape == (256, 256)
+        assert img.band_names == ["green_b1"]
 
-        img = stac.tile(71, 102, 8, expression="green/red")
+        img = stac.tile(71, 102, 8, assets=("green", "red"))
+        assert img.data.shape == (2, 256, 256)
+        assert img.mask.shape == (256, 256)
+        assert img.band_names == ["green_b1", "red_b1"]
+
+        img = stac.tile(71, 102, 8, expression="green_b1/red_b1")
         assert img.data.shape == (1, 256, 256)
         assert img.mask.shape == (256, 256)
-        # Note: Here we loose the information about the band
-        assert img.band_names == ["green/red"]
+        assert img.band_names == ["green_b1/red_b1"]
 
         with pytest.warns(ExpressionMixingWarning):
-            img = stac.tile(71, 102, 8, assets=("green", "red"), expression="green/red")
+            img = stac.tile(
+                71, 102, 8, assets=("green", "red"), expression="green_b1/red_b1"
+            )
             assert img.data.shape == (1, 256, 256)
-            assert img.band_names == ["green/red"]
+            assert img.band_names == ["green_b1/red_b1"]
 
         img = stac.tile(
             71,
@@ -170,7 +177,7 @@ def test_tile_valid(rio):
         )
         assert img.data.shape == (3, 256, 256)
         assert img.mask.shape == (256, 256)
-        assert img.band_names == ["green_1", "green_1", "red_1"]
+        assert img.band_names == ["green_b1", "green_b1", "red_b1"]
 
         # check backward compatibility for `indexes`
         img = stac.tile(
@@ -178,23 +185,27 @@ def test_tile_valid(rio):
             102,
             8,
             assets=("green", "red"),
-            indexes=1,
+            indexes=(1, 1),
         )
-        assert img.data.shape == (2, 256, 256)
+        assert img.data.shape == (4, 256, 256)
         assert img.mask.shape == (256, 256)
-        assert img.band_names == ["green_1", "red_1"]
+        assert img.band_names == ["green_b1", "green_b1", "red_b1", "red_b1"]
 
-        img = stac.tile(
-            71,
-            102,
-            8,
-            assets=("green", "red"),
-            asset_expression={"green": "b1*2;b1", "red": "b1*2"},
-        )
+        img = stac.tile(71, 102, 8, expression="green_b1*2;green_b1;red_b1*2")
         assert img.data.shape == (3, 256, 256)
         assert img.mask.shape == (256, 256)
         assert img.band_names == ["green_b1*2", "green_b1", "red_b1*2"]
 
+        # Should raise KeyError because of missing band 2
+        with pytest.raises(KeyError):
+            img = stac.tile(
+                71,
+                102,
+                8,
+                expression="green_b1/red_b2",
+                asset_indexes={"green": 1, "red": 1},
+            )
+
 
 @patch("rio_tiler.io.cogeo.rasterio")
 def test_part_valid(rio):
@@ -214,25 +225,25 @@ def test_part_valid(rio):
         img = stac.part(bbox, assets="green")
         assert img.data.shape == (1, 73, 83)
         assert img.mask.shape == (73, 83)
-        assert img.band_names == ["green_1"]
+        assert img.band_names == ["green_b1"]
 
-        data, mask = stac.part(bbox, assets=("green",))
-        assert data.shape == (1, 73, 83)
-        assert mask.shape == (73, 83)
+        img = stac.part(bbox, assets=("green",))
+        assert img.data.shape == (1, 73, 83)
+        assert img.mask.shape == (73, 83)
 
-        img = stac.part(bbox, expression="green/red")
+        img = stac.part(bbox, expression="green_b1/red_b1")
         assert img.data.shape == (1, 73, 83)
         assert img.mask.shape == (73, 83)
-        assert img.band_names == ["green/red"]
+        assert img.band_names == ["green_b1/red_b1"]
 
-        data, mask = stac.part(bbox, assets="green", max_size=30)
-        assert data.shape == (1, 27, 30)
-        assert mask.shape == (27, 30)
+        img = stac.part(bbox, assets="green", max_size=30)
+        assert img.data.shape == (1, 27, 30)
+        assert img.mask.shape == (27, 30)
 
         with pytest.warns(ExpressionMixingWarning):
-            img = stac.part(bbox, assets=("green", "red"), expression="green/red")
+            img = stac.part(bbox, assets=("green", "red"), expression="green_b1/red_b1")
             assert img.data.shape == (1, 73, 83)
-            assert img.band_names == ["green/red"]
+            assert img.band_names == ["green_b1/red_b1"]
 
         img = stac.part(
             bbox,
@@ -247,18 +258,14 @@ def test_part_valid(rio):
         )
         assert img.data.shape == (3, 73, 83)
         assert img.mask.shape == (73, 83)
-        assert img.band_names == ["green_1", "green_1", "red_1"]
+        assert img.band_names == ["green_b1", "green_b1", "red_b1"]
 
         img = stac.part(bbox, assets=("green", "red"), indexes=1)
         assert img.data.shape == (2, 73, 83)
         assert img.mask.shape == (73, 83)
-        assert img.band_names == ["green_1", "red_1"]
+        assert img.band_names == ["green_b1", "red_b1"]
 
-        img = stac.part(
-            bbox,
-            assets=("green", "red"),
-            asset_expression={"green": "b1*2;b1", "red": "b1*2"},
-        )
+        img = stac.part(bbox, expression="green_b1*2;green_b1;red_b1*2")
         assert img.data.shape == (3, 73, 83)
         assert img.mask.shape == (73, 83)
         assert img.band_names == ["green_b1*2", "green_b1", "red_b1*2"]
@@ -280,21 +287,21 @@ def test_preview_valid(rio):
         img = stac.preview(assets="green")
         assert img.data.shape == (1, 259, 255)
         assert img.mask.shape == (259, 255)
-        assert img.band_names == ["green_1"]
+        assert img.band_names == ["green_b1"]
 
-        data, mask = stac.preview(assets=("green",))
-        assert data.shape == (1, 259, 255)
-        assert mask.shape == (259, 255)
+        img = stac.preview(assets=("green",))
+        assert img.data.shape == (1, 259, 255)
+        assert img.mask.shape == (259, 255)
 
-        img = stac.preview(expression="green/red")
+        img = stac.preview(expression="green_b1/red_b1")
         assert img.data.shape == (1, 259, 255)
         assert img.mask.shape == (259, 255)
-        assert img.band_names == ["green/red"]
+        assert img.band_names == ["green_b1/red_b1"]
 
         with pytest.warns(ExpressionMixingWarning):
-            img = stac.preview(assets=("green", "red"), expression="green/red")
+            img = stac.preview(assets=("green", "red"), expression="green_b1/red_b1")
             assert img.data.shape == (1, 259, 255)
-            assert img.band_names == ["green/red"]
+            assert img.band_names == ["green_b1/red_b1"]
 
         img = stac.preview(
             assets=("green", "red"),
@@ -308,17 +315,14 @@ def test_preview_valid(rio):
         )
         assert img.data.shape == (3, 259, 255)
         assert img.mask.shape == (259, 255)
-        assert img.band_names == ["green_1", "green_1", "red_1"]
+        assert img.band_names == ["green_b1", "green_b1", "red_b1"]
 
         img = stac.preview(assets=("green", "red"), indexes=1)
         assert img.data.shape == (2, 259, 255)
         assert img.mask.shape == (259, 255)
-        assert img.band_names == ["green_1", "red_1"]
+        assert img.band_names == ["green_b1", "red_b1"]
 
-        img = stac.preview(
-            assets=("green", "red"),
-            asset_expression={"green": "b1*2;b1", "red": "b1*2"},
-        )
+        img = stac.preview(expression="green_b1*2;green_b1;red_b1*2")
         assert img.data.shape == (3, 259, 255)
         assert img.mask.shape == (259, 255)
         assert img.band_names == ["green_b1*2", "green_b1", "red_b1*2"]
@@ -337,50 +341,51 @@ def test_point_valid(rio):
         with pytest.raises(MissingAssets):
             stac.point(-80.477, 33.4453)
 
-        data = stac.point(-80.477, 33.4453, assets="green")
+        data, names = stac.point(-80.477, 33.4453, assets="green")
         assert len(data) == 1
+        assert names == ["green_b1"]
 
-        data = stac.point(-80.477, 33.4453, assets=("green",))
+        data, names = stac.point(-80.477, 33.4453, assets=("green",))
         assert len(data) == 1
+        assert names == ["green_b1"]
 
-        data = stac.point(-80.477, 33.4453, expression="green/red")
+        data, names = stac.point(-80.477, 33.4453, assets=("green", "red"))
+        assert len(data) == 2
+        assert data == [7994, 7003]
+        assert names == ["green_b1", "red_b1"]
+
+        data, names = stac.point(-80.477, 33.4453, expression="green_b1/red_b1")
         assert len(data) == 1
+        assert data == [7994 / 7003]
+        assert names == ["green_b1/red_b1"]
 
         with pytest.warns(ExpressionMixingWarning):
-            data = stac.point(
-                -80.477, 33.4453, assets=("green", "red"), expression="green/red"
+            data, names = stac.point(
+                -80.477, 33.4453, assets=("green", "red"), expression="green_b1/red_b1"
             )
             assert len(data) == 1
+            assert names == ["green_b1/red_b1"]
 
-        data = stac.point(
+        data, names = stac.point(
             -80.477,
             33.4453,
             assets=("green", "red"),
             asset_indexes={"green": (1, 1), "red": 1},
         )
-        assert len(data) == 2
-        assert len(data[0]) == 2
-        assert len(data[1]) == 1
+        assert len(data) == 3
+        assert data == [7994, 7994, 7003]
+        assert names == ["green_b1", "green_b1", "red_b1"]
 
-        data = stac.point(
-            -80.477,
-            33.4453,
-            assets=("green", "red"),
-            indexes=1,
-        )
+        data, names = stac.point(-80.477, 33.4453, assets=("green", "red"), indexes=1)
         assert len(data) == 2
-        assert len(data[0]) == 1
-        assert len(data[1]) == 1
+        assert data == [7994, 7003]
+        assert names == ["green_b1", "red_b1"]
 
-        data = stac.point(
-            -80.477,
-            33.4453,
-            assets=("green", "red"),
-            asset_expression={"green": "b1*2;b1", "red": "b1*2"},
+        data, names = stac.point(
+            -80.477, 33.4453, expression="green_b1*2;green_b1;red_b1*2"
         )
-        assert len(data) == 2
-        assert len(data[0]) == 2
-        assert len(data[1]) == 1
+        assert len(data) == 3
+        assert names == ["green_b1*2", "green_b1", "red_b1*2"]
 
 
 @patch("rio_tiler.io.cogeo.rasterio")
@@ -400,11 +405,11 @@ def test_statistics_valid(rio):
 
         stats = stac.statistics(assets="green")
         assert stats["green"]
-        assert isinstance(stats["green"]["1"], BandStatistics)
+        assert isinstance(stats["green"]["b1"], BandStatistics)
 
         stats = stac.statistics(assets=("green", "red"), hist_options={"bins": 20})
         assert len(stats) == 2
-        assert len(stats["green"]["1"]["histogram"][0]) == 20
+        assert len(stats["green"]["b1"]["histogram"][0]) == 20
 
         # Check that asset_expression is passed
         stats = stac.statistics(
@@ -419,14 +424,14 @@ def test_statistics_valid(rio):
             assets=("green", "red"), asset_indexes={"green": 1, "red": 1}
         )
         assert stats["green"]
-        assert isinstance(stats["green"]["1"], BandStatistics)
-        assert isinstance(stats["red"]["1"], BandStatistics)
+        assert isinstance(stats["green"]["b1"], BandStatistics)
+        assert isinstance(stats["red"]["b1"], BandStatistics)
 
         # Check that asset_indexes is passed
         stats = stac.statistics(assets=("green", "red"), indexes=1)
         assert stats["green"]
-        assert isinstance(stats["green"]["1"], BandStatistics)
-        assert isinstance(stats["red"]["1"], BandStatistics)
+        assert isinstance(stats["green"]["b1"], BandStatistics)
+        assert isinstance(stats["red"]["b1"], BandStatistics)
 
 
 @patch("rio_tiler.io.cogeo.rasterio")
@@ -438,28 +443,25 @@ def test_merged_statistics_valid(rio):
         with pytest.warns(UserWarning):
             stats = stac.merged_statistics()
             assert len(stats) == 3
-            assert isinstance(stats["red_1"], BandStatistics)
-            assert stats["red_1"]
-            assert stats["green_1"]
-            assert stats["blue_1"]
+            assert isinstance(stats["red_b1"], BandStatistics)
+            assert stats["red_b1"]
+            assert stats["green_b1"]
+            assert stats["blue_b1"]
 
         with pytest.raises(InvalidAssetName):
             stac.merged_statistics(assets="vert")
 
         stats = stac.merged_statistics(assets="green")
-        assert isinstance(stats["green_1"], BandStatistics)
+        assert isinstance(stats["green_b1"], BandStatistics)
 
         stats = stac.merged_statistics(
             assets=("green", "red"), hist_options={"bins": 20}
         )
         assert len(stats) == 2
-        assert len(stats["green_1"]["histogram"][0]) == 20
-        assert len(stats["red_1"]["histogram"][0]) == 20
+        assert len(stats["green_b1"]["histogram"][0]) == 20
+        assert len(stats["red_b1"]["histogram"][0]) == 20
 
-        # Check that asset_expression is passed
-        stats = stac.merged_statistics(
-            assets=("green", "red"), asset_expression={"green": "b1*2", "red": "b1+100"}
-        )
+        stats = stac.merged_statistics(expression="green_b1*2;green_b1;red_b1+100")
         assert isinstance(stats["green_b1*2"], BandStatistics)
         assert isinstance(stats["red_b1+100"], BandStatistics)
 
@@ -467,19 +469,8 @@ def test_merged_statistics_valid(rio):
         stats = stac.merged_statistics(
             assets=("green", "red"), asset_indexes={"green": 1, "red": 1}
         )
-        assert isinstance(stats["green_1"], BandStatistics)
-        assert isinstance(stats["red_1"], BandStatistics)
-
-        # Check Expression
-        stats = stac.merged_statistics(expression="green/red")
-        assert isinstance(stats["green/red"], BandStatistics)
-
-        # Check that we can use expression and asset_expression
-        stats = stac.merged_statistics(
-            expression="green/red",
-            asset_expression={"green": "b1*2", "red": "b1+100"},
-        )
-        assert isinstance(stats["green/red"], BandStatistics)
+        assert isinstance(stats["green_b1"], BandStatistics)
+        assert isinstance(stats["red_b1"], BandStatistics)
 
 
 @patch("rio_tiler.io.cogeo.rasterio")
@@ -506,14 +497,22 @@ def test_info_valid(rio):
 
 
 def test_parse_expression():
-    """."""
+    """Parse assets expressions."""
     with STACReader(STAC_PATH) as stac:
-        assert sorted(stac.parse_expression("green*red+red/blue+2.0")) == [
+        assert sorted(
+            stac.parse_expression("green_b1*red_b1+red_b1/blue_b1+2.0;red_b1")
+        ) == [
             "blue",
             "green",
             "red",
         ]
 
+    # make sure we match full word only
+    with STACReader(STAC_PATH) as stac:
+        assert sorted(
+            stac.parse_expression("greenish_b1*red_b1+red_b1/blue_b1+2.0;red_b1")
+        ) == ["blue", "red"]
+
 
 @patch("rio_tiler.io.cogeo.rasterio")
 def test_feature_valid(rio):
@@ -553,43 +552,41 @@ def test_feature_valid(rio):
         img = stac.feature(feat, assets="green")
         assert img.data.shape == (1, 118, 96)
         assert img.mask.shape == (118, 96)
-        assert img.band_names == ["green_1"]
+        assert img.band_names == ["green_b1"]
 
-        data, mask = stac.feature(feat, assets=("green",))
-        assert data.shape == (1, 118, 96)
-        assert mask.shape == (118, 96)
+        img = stac.feature(feat, assets=("green",))
+        assert img.data.shape == (1, 118, 96)
+        assert img.mask.shape == (118, 96)
 
-        img = stac.feature(feat, expression="green/red")
+        img = stac.feature(feat, expression="green_b1/red_b1")
         assert img.data.shape == (1, 118, 96)
         assert img.mask.shape == (118, 96)
-        assert img.band_names == ["green/red"]
+        assert img.band_names == ["green_b1/red_b1"]
 
-        data, mask = stac.feature(feat, assets="green", max_size=30)
-        assert data.shape == (1, 30, 25)
-        assert mask.shape == (30, 25)
+        img = stac.feature(feat, assets="green", max_size=30)
+        assert img.data.shape == (1, 30, 25)
+        assert img.mask.shape == (30, 25)
 
         with pytest.warns(ExpressionMixingWarning):
-            img = stac.feature(feat, assets=("green", "red"), expression="green/red")
+            img = stac.feature(
+                feat, assets=("green", "red"), expression="green_b1/red_b1"
+            )
             assert img.data.shape == (1, 118, 96)
-            assert img.band_names == ["green/red"]
+            assert img.band_names == ["green_b1/red_b1"]
 
         img = stac.feature(
             feat, assets=("green", "red"), asset_indexes={"green": (1, 1), "red": 1}
         )
         assert img.data.shape == (3, 118, 96)
         assert img.mask.shape == (118, 96)
-        assert img.band_names == ["green_1", "green_1", "red_1"]
+        assert img.band_names == ["green_b1", "green_b1", "red_b1"]
 
         img = stac.feature(feat, assets=("green", "red"), indexes=1)
         assert img.data.shape == (2, 118, 96)
         assert img.mask.shape == (118, 96)
-        assert img.band_names == ["green_1", "red_1"]
+        assert img.band_names == ["green_b1", "red_b1"]
 
-        img = stac.feature(
-            feat,
-            assets=("green", "red"),
-            asset_expression={"green": "b1*2;b1", "red": "b1*2"},
-        )
+        img = stac.feature(feat, expression="green_b1*2;green_b1;red_b1*2")
         assert img.data.shape == (3, 118, 96)
         assert img.mask.shape == (118, 96)
         assert img.band_names == ["green_b1*2", "green_b1", "red_b1*2"]
diff --git a/tests/test_models.py b/tests/test_models.py
index eea50727..3ba00ae0 100644
--- a/tests/test_models.py
+++ b/tests/test_models.py
@@ -121,3 +121,17 @@ def test_merge_with_diffsize():
         img2 = ImageData(numpy.zeros((1, 256, 256)))
         img = ImageData.create_from_list([img1, img2])
     assert len(w) == 0
+
+
+def test_apply_expression():
+    """Apply expression"""
+    img = ImageData(numpy.zeros((2, 256, 256)))
+    img2 = img.apply_expression("b1+b2")
+    assert img.count == 2
+    assert img.width == 256
+    assert img.height == 256
+    assert img.band_names == ["b1", "b2"]
+    assert img2.count == 1
+    assert img2.width == 256
+    assert img2.height == 256
+    assert img2.band_names == ["b1+b2"]
diff --git a/tests/test_mosaic.py b/tests/test_mosaic.py
index d2d76f39..da147095 100644
--- a/tests/test_mosaic.py
+++ b/tests/test_mosaic.py
@@ -69,10 +69,10 @@ def test_mosaic_tiler():
     assert t.dtype == m.dtype
 
     img, _ = mosaic.mosaic_reader(assets, _read_tile, x, y, z)
-    assert img.band_names == ["1", "2", "3"]
+    assert img.band_names == ["b1", "b2", "b3"]
 
     img, _ = mosaic.mosaic_reader(assets, _read_tile, x, y, z, indexes=[1])
-    assert img.band_names == ["1"]
+    assert img.band_names == ["b1"]
 
     img, _ = mosaic.mosaic_reader(assets, _read_tile, x, y, z, expression="b1*3")
     assert img.band_names == ["b1*3"]
@@ -281,7 +281,7 @@ def _reader(src_path: str, *args, **kwargs) -> ImageData:
         assets="green",
         threads=0,
     )
-    assert img.band_names == ["green_1"]
+    assert img.band_names == ["green_b1"]
 
     img, _ = mosaic.mosaic_reader(
         [stac_asset],
@@ -289,23 +289,11 @@ def _reader(src_path: str, *args, **kwargs) -> ImageData:
         71,
         102,
         8,
-        assets=["green"],
-        asset_expression={"green": "b1*2"},
+        expression="green_b1*2",
         threads=0,
     )
     assert img.band_names == ["green_b1*2"]
 
-    img, _ = mosaic.mosaic_reader(
-        [stac_asset],
-        _reader,
-        71,
-        102,
-        8,
-        expression="green*2",
-        threads=0,
-    )
-    assert img.band_names == ["green*2"]
-
 
 def test_mosaic_tiler_Stdev():
     """Test Stdev mosaic methods."""
diff --git a/tests/test_reader.py b/tests/test_reader.py
index 2c1ffd4c..90d4e30a 100644
--- a/tests/test_reader.py
+++ b/tests/test_reader.py
@@ -369,24 +369,32 @@ def test_tile_read_vrt_option():
 def test_point():
     """Read point values"""
     with rasterio.open(COG_SCALE) as src_dst:
-        p = reader.point(src_dst, [310000, 4100000], coord_crs=src_dst.crs, indexes=1)
+        p, name = reader.point(
+            src_dst, [310000, 4100000], coord_crs=src_dst.crs, indexes=1
+        )
         assert p == [8917]
+        assert name == ["b1"]
 
-        p = reader.point(src_dst, [310000, 4100000], coord_crs=src_dst.crs)
+        p, name = reader.point(src_dst, [310000, 4100000], coord_crs=src_dst.crs)
         assert p == [8917]
+        assert name == ["b1"]
 
         with pytest.raises(PointOutsideBounds):
             reader.point(src_dst, [810000, 4100000], coord_crs=src_dst.crs)
 
     with rasterio.open(S3_MASK_PATH) as src_dst:
         # Test with COG + internal mask
-        assert not reader.point(src_dst, [-104.7753105, 38.953548])[0]
-        assert reader.point(src_dst, [-104.7753105415, 38.953548], masked=False)[0] == 0
+        assert not reader.point(src_dst, [-104.7753105, 38.953548])[0][0]
+        assert (
+            reader.point(src_dst, [-104.7753105415, 38.953548], masked=False)[0][0] == 0
+        )
 
     with rasterio.open(S3_ALPHA_PATH) as src_dst:
         # Test with COG + Alpha Band
-        assert not reader.point(src_dst, [-104.77519499, 38.95367054])[0]
-        assert reader.point(src_dst, [-104.77519499, 38.95367054], masked=False)[0] == 0
+        assert not reader.point(src_dst, [-104.77519499, 38.95367054])[0][0]
+        assert (
+            reader.point(src_dst, [-104.77519499, 38.95367054], masked=False)[0][0] == 0
+        )
 
 
 def test_part_with_buffer():