ENH: Added action normalize that can normalize a FeatureTable[Frequency] by gene length, library size and composition.

README.md

@@ -15,7 +15,7 @@ To install _q2-amr_, follow the steps described below.
 mamba create -yn q2-amr \
   -c https://packages.qiime2.org/qiime2/2024.2/shotgun/released/  \
   -c qiime2 -c conda-forge -c bioconda -c defaults \
-  qiime2 q2cli q2templates q2-types rgi
+  qiime2 q2cli q2templates q2-types rgi rnanorm
 
 conda activate q2-amr
 
@@ -40,7 +40,7 @@ qiime info
 CONDA_SUBDIR=osx-64 mamba create -yn q2-amr \
   -c https://packages.qiime2.org/qiime2/2024.2/shotgun/released/ \
   -c qiime2 -c conda-forge -c bioconda -c defaults \
-  qiime2 q2cli q2templates q2-types rgi
+  qiime2 q2cli q2templates q2-types rgi rnanorm
 
 conda activate q2-amr
 conda config --env --set subdir osx-64
@@ -62,15 +62,17 @@ qiime info
 ## Functionality
 This QIIME 2 plugin contains actions used to annotate short single/paired-end
 sequencing reads and MAGs with antimicrobial resistance genes. Currently, the [CARD](https://card.mcmaster.ca) database is supported  (for details on
-the implementation and usage, please refer to the [rgi](https://github.com/arpcard/rgi) documentation). Below you will
+the implementation and usage, please refer to the [RGI](https://github.com/arpcard/rgi) documentation). Below you will
 find an overview of actions available in the plugin.
 
-| Action                     | Description                                                                          | Underlying tool                       | Used function                        |
-|----------------------------|--------------------------------------------------------------------------------------|---------------------------------------|--------------------------------------|
-| fetch-card-db              | Download and preprocess CARD and WildCARD data.                                      | [rgi](https://github.com/arpcard/rgi) | card_annotation, wildcard_annotation |
-| annotate-mags-card         | Annotate MAGs with antimicrobial resistance gene information from CARD.              | [rgi](https://github.com/arpcard/rgi) | main, load                           |
-| annotate-reads-card        | Annotate metagenomic reads with antimicrobial resistance gene information from CARD. | [rgi](https://github.com/arpcard/rgi) | bwt, load                            |
-| heatmap                    | Create a heatmap from annotate-mags-card output files.                               | [rgi](https://github.com/arpcard/rgi) | heatmap                              |
+| Action              | Description                                                                          | Underlying tool                           | Used function                        |
+|---------------------|--------------------------------------------------------------------------------------|-------------------------------------------|--------------------------------------|
+| fetch-card-db       | Download and preprocess CARD and WildCARD data.                                      | [RGI](https://github.com/arpcard/rgi)     | card_annotation, wildcard_annotation |
+| annotate-mags-card  | Annotate MAGs with antimicrobial resistance gene information from CARD.              | [RGI](https://github.com/arpcard/rgi)     | main, load                           |
+| annotate-reads-card | Annotate metagenomic reads with antimicrobial resistance gene information from CARD. | [RGI](https://github.com/arpcard/rgi)     | bwt, load                            |
+| heatmap             | Create a heatmap from annotate-mags-card output files.                               | [RGI](https://github.com/arpcard/rgi)     | heatmap                              |
+| normalize           | Normalize feature table with [FPKM](https://www.nature.com/articles/nmeth.1226), [TPM](https://link.springer.com/article/10.1007/s12064-012-0162-3), [TMM](https://genomebiology.biomedcentral.com/articles/10.1186/gb-2010-11-3-r25), [UQ](https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-11-94), [CUF](https://genomebiology.biomedcentral.com/articles/10.1186/s13059-021-02568-9/), or [CTF](https://genomebiology.biomedcentral.com/articles/10.1186/s13059-021-02568-9/) method.                 | [rnaNORM](https://github.com/genialis/RNAnorm) | fpkm, tpm, tmm, uq, cuf, ctf         |
+
 
 ## Dev environment
 This repository follows the _black_ code style. To make the development slightly easier

ci/recipe/meta.yaml

@@ -23,6 +23,7 @@ requirements:
     - q2templates {{ qiime2_epoch }}.*
     - q2cli {{ qiime2_epoch }}.*
     - rgi
+    - rnanorm
     - tqdm
 
 test:

q2_amr/card/heatmap.py

@@ -8,7 +8,7 @@
 import pkg_resources
 import q2templates
 
-from q2_amr.card.utils import run_command
+from q2_amr.card.utils import InvalidParameterCombinationError, run_command
 from q2_amr.types import CARDAnnotationDirectoryFormat
 
 
@@ -50,12 +50,6 @@ def heatmap(
     q2templates.render(templates, output_dir, context=context)
 
 
-class InvalidParameterCombinationError(Exception):
-    def __init__(self, message="Invalid parameter combination"):
-        self.message = message
-        super().__init__(self.message)
-
-
 def run_rgi_heatmap(tmp, json_files_dir, clus, cat, display, frequency):
     cmd = [
         "rgi",

q2_amr/card/normalization.py

@@ -0,0 +1,90 @@
+import os
+
+import pandas as pd
+from q2_types.feature_data import SequenceCharacteristicsDirectoryFormat
+from rnanorm import CPM, CTF, CUF, FPKM, TMM, TPM, UQ
+
+
+def normalize(
+    table: pd.DataFrame,
+    method: str,
+    m_trim: float = None,
+    a_trim: float = None,
+    gene_length: SequenceCharacteristicsDirectoryFormat = None,
+) -> pd.DataFrame:
+    # Validate parameter combinations and set trim parameters
+    m_trim, a_trim = _validate_parameters(method, m_trim, a_trim, gene_length)
+
+    # Process gene_lengths input and define methods that need gene_lengths input
+    if method in ["tpm", "fpkm"]:
+        lengths = _convert_lengths(table, gene_length)
+
+        methods = {
+            "tpm": TPM(gene_lengths=lengths),
+            "fpkm": FPKM(gene_lengths=lengths),
+        }
+
+    # Define remaining methods that don't need gene_lengths input
+    else:
+        methods = {
+            "tmm": TMM(m_trim=m_trim, a_trim=a_trim),
+            "ctf": CTF(m_trim=m_trim, a_trim=a_trim),
+            "uq": UQ(),
+            "cuf": CUF(),
+            "cpm": CPM(),
+        }
+
+    # Run normalization method on frequency table
+    normalized = methods[method].set_output(transform="pandas").fit_transform(table)
+    normalized.index.name = "sample_id"
+
+    return normalized
+
+
+def _validate_parameters(method, m_trim, a_trim, gene_length):
+    # Raise Error if gene-length is missing when using methods TPM or FPKM
+    if method in ["tpm", "fpkm"] and not gene_length:
+        raise ValueError("gene-length input is missing.")
+
+    # Raise Error if gene-length is given when using methods TMM, UQ, CUF, CPM or CTF
+    if method in ["tmm", "uq", "cuf", "ctf", "cpm"] and gene_length:
+        raise ValueError(
+            "gene-length input can only be used with FPKM and TPM methods."
+        )
+
+    # Raise Error if m_trim or a_trim are given when not using methods TMM or CTF
+    if (method not in ["tmm", "ctf"]) and (m_trim is not None or a_trim is not None):
+        raise ValueError(
+            "Parameters m-trim and a-trim can only be used with methods TMM and CTF."
+        )
+
+    # Set m_trim and a_trim to their default values for methods TMM and CTF
+    if method in ["tmm", "ctf"]:
+        m_trim = 0.3 if m_trim is None else m_trim
+        a_trim = 0.05 if a_trim is None else a_trim
+
+    return m_trim, a_trim
+
+
+def _convert_lengths(table, gene_length):
+    # Read in table from sequence_characteristics.tsv as a pd.Series
+    lengths = pd.read_csv(
+        os.path.join(gene_length.path, "sequence_characteristics.tsv"),
+        sep="\t",
+        header=None,
+        names=["index", "values"],
+        index_col="index",
+        squeeze=True,
+        skiprows=1,
+    )
+
+    # Check if all gene IDs that are present in the table are also present in
+    # the lengths
+    if not set(table.columns).issubset(set(lengths.index)):
+        only_in_counts = set(table.columns) - set(lengths.index)
+        raise ValueError(
+            f"There are genes present in the FeatureTable that are not present "
+            f"in the gene-length input. Missing lengths for genes: "
+            f"{only_in_counts}"
+        )
+    return lengths

q2_amr/card/tests/data/feature-table.tsv

@@ -0,0 +1,3 @@
+ID	ARO:3000026|ID:377|Name:mepA|NCBI:AY661734.1	ARO:3000027|ID:1757|Name:emrA|NCBI:AP009048.1
+sample1	2.0	0.0
+sample2	2.0	0.0