Repository: taxprofiler/taxpasta

Issue:
[Feature]  Add a taxonomy table to the BIOM output
### Checklist

- [X] There are [no similar issues or pull requests](https://github.com/taxprofiler/taxpasta/issues) for this yet.

### Problem

I would like to create a phyloseq object directly from a BIOM format file. This fails at the moment because of missing taxonomy information.

### Solution

Add the taxonomy table to the file.

### Alternatives

_No response_

### Anything else?

_No response_

Discussion:

Patch:
diff --git a/src/taxpasta/domain/service/taxonomy_service.py b/src/taxpasta/domain/service/taxonomy_service.py
index 982cec3..9f26450 100644
--- a/src/taxpasta/domain/service/taxonomy_service.py
+++ b/src/taxpasta/domain/service/taxonomy_service.py
@@ -20,7 +20,7 @@
 
 
 from abc import ABC, abstractmethod
-from typing import List, Optional, TypeVar
+from typing import Dict, List, Optional, TypeVar
 
 from pandera.typing import DataFrame
 
@@ -81,6 +81,12 @@ class TaxonomyService(ABC):
     def add_rank_lineage(self, table: DataFrame[ResultTable]) -> DataFrame[ResultTable]:
         """Add a column for the taxon lineage as ranks to the given table."""
 
+    @abstractmethod
+    def format_biom_taxonomy(
+        self, table: DataFrame[ResultTable]
+    ) -> List[Dict[str, List[str]]]:
+        """Format the taxonomy as BIOM observation metadata."""
+
     @abstractmethod
     def summarise_at(
         self, profile: DataFrame[StandardProfile], rank: str
diff --git a/src/taxpasta/infrastructure/application/wide_observation_table_writer/biom_wide_observation_table_writer.py b/src/taxpasta/infrastructure/application/wide_observation_table_writer/biom_wide_observation_table_writer.py
index 1c572dc..c4bac0e 100644
--- a/src/taxpasta/infrastructure/application/wide_observation_table_writer/biom_wide_observation_table_writer.py
+++ b/src/taxpasta/infrastructure/application/wide_observation_table_writer/biom_wide_observation_table_writer.py
@@ -47,15 +47,16 @@ class BIOMWideObservationTableWriter(WideObservationTableWriter):
         # Drop unclassified reads.
         matrix = matrix.loc[matrix.iloc[:, 0] != 0].copy()
         if taxonomy is not None:
-            observation_meta = None
-            # FIXME (Moritz): All lists need to be of equal length, i.e., per key in
-            #  the dict across all observations (aka samples).
+            observation_meta = taxonomy.format_biom_taxonomy(matrix)
+            tmp = taxonomy.add_rank_lineage(matrix)
+            for ranks, meta in zip(tmp.rank_lineage, observation_meta):
+                meta["rank_lineage"] = ranks
         else:
             observation_meta = None
         result = Table(
             data=matrix.iloc[:, 1:].values,
             observation_ids=matrix.iloc[:, 0].astype(str),
-            sample_ids=matrix.columns[1:],
+            sample_ids=matrix.columns[1:].astype(str),
             observation_metadata=observation_meta,
             create_date=datetime.utcnow().isoformat(timespec="microseconds"),
         )
diff --git a/src/taxpasta/infrastructure/domain/service/taxopy_taxonomy_service.py b/src/taxpasta/infrastructure/domain/service/taxopy_taxonomy_service.py
index 4141e9f..450cc17 100644
--- a/src/taxpasta/infrastructure/domain/service/taxopy_taxonomy_service.py
+++ b/src/taxpasta/infrastructure/domain/service/taxopy_taxonomy_service.py
@@ -24,7 +24,7 @@ from __future__ import annotations
 import logging
 from collections import defaultdict
 from pathlib import Path
-from typing import List, Optional
+from typing import Dict, List, Optional
 
 import pandas as pd
 import taxopy
@@ -68,20 +68,30 @@ class TaxopyTaxonomyService(TaxonomyService):
         return self._tax_db.taxid2rank.get(taxonomy_id)
 
     def get_taxon_name_lineage(self, taxonomy_id: int) -> Optional[List[str]]:
-        """Return the lineage of a given taxonomy identifier as names."""
+        """
+        Return the lineage of a given taxonomy identifier as names.
+
+        Only names with associated ranks are included.
+
+        """
         try:
             taxon = taxopy.Taxon(taxid=taxonomy_id, taxdb=self._tax_db)
         except TaxidError:
             return None
-        return list(reversed(taxon.name_lineage))
+        return list(reversed(taxon.rank_name_dictionary.values()))
 
     def get_taxon_identifier_lineage(self, taxonomy_id: int) -> Optional[List[int]]:
-        """Return the lineage of a given taxonomy identifier as identifiers."""
+        """
+        Return the lineage of a given taxonomy identifier as identifiers.
+
+        Only identifiers with associated ranks are included.
+
+        """
         try:
             taxon = taxopy.Taxon(taxid=taxonomy_id, taxdb=self._tax_db)
         except TaxidError:
             return None
-        return list(reversed(taxon.taxid_lineage))
+        return list(reversed(taxon.rank_taxid_dictionary.values()))
 
     def get_taxon_rank_lineage(self, taxonomy_id: int) -> Optional[List[str]]:
         """Return the lineage of a given taxonomy identifier as ranks."""
@@ -164,6 +174,29 @@ class TaxopyTaxonomyService(TaxonomyService):
         else:
             return None
 
+    def format_biom_taxonomy(
+        self, table: DataFrame[ResultTable]
+    ) -> List[Dict[str, List[str]]]:
+        """Format the taxonomy as BIOM observation metadata."""
+        lineages = [self._get_rank_name(tax_id) for tax_id in table.taxonomy_id]
+        result = [{} if lineage is None else lineage for lineage in lineages]
+        longest_lineage = max(result, key=len)
+        max_ranks = list(reversed(longest_lineage.keys()))
+        return [
+            {"taxonomy": self._pad_lineage(lineage, max_ranks)} for lineage in result
+        ]
+
+    def _get_rank_name(self, taxonomy_id: int) -> Optional[Dict[str, str]]:
+        try:
+            taxon = taxopy.Taxon(taxid=taxonomy_id, taxdb=self._tax_db)
+        except TaxidError:
+            return None
+        return taxon.rank_name_dictionary
+
+    def _pad_lineage(self, lineage: Dict[str, str], max_ranks: List[str]) -> List[str]:
+        """Pad a lineage to match the length of the longest lineage."""
+        return [lineage.get(rank, "") for rank in max_ranks]
+
     def summarise_at(
         self, profile: DataFrame[StandardProfile], rank: str
     ) -> DataFrame[StandardProfile]:


Test patch:
diff --git a/tests/system/cli/test_taxpasta_merge.py b/tests/system/cli/test_taxpasta_merge.py
index dab3c0e..dc4b918 100644
--- a/tests/system/cli/test_taxpasta_merge.py
+++ b/tests/system/cli/test_taxpasta_merge.py
@@ -67,7 +67,7 @@ def test_merge_profiles_wide(
     runner: CliRunner,
     profiler: SupportedProfiler,
     profiles: List[str],
-    wide_observation_table_format,
+    wide_observation_table_format: WideObservationTableFileFormat,
     tmp_path: Path,
     monkeypatch: pytest.MonkeyPatch,
 ):
diff --git a/tests/unit/infrastructure/domain/service/test_taxopy_taxonomy_service.py b/tests/unit/infrastructure/domain/service/test_taxopy_taxonomy_service.py
index 1c3e2ab..a56755a 100644
--- a/tests/unit/infrastructure/domain/service/test_taxopy_taxonomy_service.py
+++ b/tests/unit/infrastructure/domain/service/test_taxopy_taxonomy_service.py
@@ -21,7 +21,7 @@
 
 from collections import OrderedDict
 from pathlib import Path
-from typing import List
+from typing import Dict, List
 
 import pandas as pd
 import pytest
@@ -77,19 +77,18 @@ def test_get_taxon_rank(tax_service: TaxopyTaxonomyService, tax_id: int, expecte
 @pytest.mark.parametrize(
     ("tax_id", "expected"),
     [
-        (1, ["root"]),
+        (1, []),
         (42, None),
         (
             86398254,
             [
-                "root",
                 "Bacteria",
                 "Proteobacteria",
                 "Gammaproteobacteria",
                 "Pseudomonadales",
             ],
         ),
-        (1199096325, ["root", "Eukaryota", "Ascomycota", "Saccharomycetes"]),
+        (1199096325, ["Eukaryota", "Ascomycota", "Saccharomycetes"]),
     ],
 )
 def test_get_taxon_name_lineage(
@@ -102,10 +101,10 @@ def test_get_taxon_name_lineage(
 @pytest.mark.parametrize(
     ("tax_id", "expected"),
     [
-        (1, [1]),
+        (1, []),
         (42, None),
-        (86398254, [1, 609216830, 1641076285, 329474883, 86398254]),
-        (1199096325, [1, 476817098, 432158898, 1199096325]),
+        (86398254, [609216830, 1641076285, 329474883, 86398254]),
+        (1199096325, [476817098, 432158898, 1199096325]),
     ],
 )
 def test_get_taxon_identifier_lineage(
@@ -143,11 +142,11 @@ def test_get_taxon_rank_lineage(
                         (
                             "lineage",
                             [
-                                "root",
                                 None,
-                                "root;Bacteria;Proteobacteria;Gammaproteobacteria;"
+                                None,
+                                "Bacteria;Proteobacteria;Gammaproteobacteria;"
                                 "Pseudomonadales",
-                                "root;Eukaryota;Ascomycota;Saccharomycetes",
+                                "Eukaryota;Ascomycota;Saccharomycetes",
                             ],
                         ),
                     ]
@@ -175,10 +174,10 @@ def test_add_name_lineage(
                         (
                             "id_lineage",
                             [
-                                "1",
                                 None,
-                                "1;609216830;1641076285;329474883;86398254",
-                                "1;476817098;432158898;1199096325",
+                                None,
+                                "609216830;1641076285;329474883;86398254",
+                                "476817098;432158898;1199096325",
                             ],
                         ),
                     ]
@@ -223,3 +222,61 @@ def test_add_rank_lineage(
 ):
     """Expect that we can add rank lineages to a result table."""
     assert_frame_equal(tax_service.add_rank_lineage(result), expected)
+
+
+@pytest.mark.parametrize(
+    ("result", "expected"),
+    [
+        (
+            pd.DataFrame(
+                OrderedDict(
+                    [("taxonomy_id", [1, 42, 1887621118, 1945799576, 492356122])]
+                )
+            ),
+            [
+                {"taxonomy": [""] * 7},
+                {"taxonomy": [""] * 7},
+                {
+                    "taxonomy": [
+                        "Bacteria",
+                        "Proteobacteria",
+                        "Gammaproteobacteria",
+                        "Pseudomonadales",
+                        "Pseudomonadaceae",
+                        "Pseudomonas",
+                        "Pseudomonas putida",
+                    ]
+                },
+                {
+                    "taxonomy": [
+                        "Bacteria",
+                        "Proteobacteria",
+                        "Gammaproteobacteria",
+                        "Enterobacterales",
+                        "Enterobacteriaceae",
+                        "Escherichia",
+                        "Escherichia coli",
+                    ]
+                },
+                {
+                    "taxonomy": [
+                        "Eukaryota",
+                        "Ascomycota",
+                        "Saccharomycetes",
+                        "Saccharomycetales",
+                        "Saccharomycetaceae",
+                        "Saccharomyces",
+                        "Saccharomyces cerevisiae",
+                    ]
+                },
+            ],
+        ),
+    ],
+)
+def test_format_biom_taxonomy(
+    tax_service: TaxopyTaxonomyService,
+    result: pd.DataFrame,
+    expected: List[Dict[str, List[str]]],
+):
+    """Expect that we can add rank lineages to a result table."""
+    assert tax_service.format_biom_taxonomy(result) == expected