Replace custom merge/chimera Python code with C-level dada2 implementations

dada2_denoise.py: remove 170 lines of custom merge_pairs() and is_bimera(),
use dada2py.merge_pairs() (C-level NW alignment) and
dada2py.remove_bimera_denovo() (C-level consensus chimera detection).

remove_chimeras.py: remove 130 lines of custom sparse chimera detection,
use dada2py.remove_bimera_denovo() with vectorized matrix construction.

Validated at 98.3-100% Jaccard with R dada2 across 10-75 samples.

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>

Author: rec3141

nextflow/bin/py/dada2_denoise.py

@@ -5,9 +5,9 @@
 # Python mirror of dada2_denoise.R. For each sample on a plate:
 #   1. Dereplicate forward and reverse filtered reads via dada2py
 #   2. Apply the DADA2 denoising algorithm using pre-learned error models
-#   3. Merge denoised forward/reverse pairs (overlap alignment)
+#   3. Merge denoised forward/reverse pairs via C-level NW alignment
 #   4. Combine all merged samples into a sequence table (pandas DataFrame)
-#   5. Per-plate chimera removal (abundance-ratio bimera check)
+#   5. Per-plate chimera removal via C-level consensus bimera detection
 #
 # Usage:
 #   dada2_denoise.py <plate_id> <errF.pkl> <errR.pkl> <min_overlap> <cpus>
@@ -82,173 +82,9 @@
 
 
 # ---------------------------------------------------------------------------
-# Merge pairs implementation
+# Merge pairs and chimera removal use dada2's C-level implementations
+# via the dada2py library (py/paired.py and py/chimera.py).
 # ---------------------------------------------------------------------------
-COMPLEMENT = str.maketrans("ACGTacgt", "TGCAtgca")
-
-
-def reverse_complement(seq):
-    """Return the reverse complement of a DNA sequence."""
-    return seq.translate(COMPLEMENT)[::-1]
-
-
-def find_best_overlap(fwd_seq, rev_rc, min_ov):
-    """Find the best overlap between the 3' end of fwd and 5' end of rev_rc.
-
-    Slides rev_rc along fwd to find the position with the most matches.
-    Returns (overlap_length, n_mismatches, merged_seq) or None if no
-    valid overlap is found.
-    """
-    fwd_len = len(fwd_seq)
-    rev_len = len(rev_rc)
-    best_overlap = 0
-    best_mismatches = float('inf')
-    best_pos = -1
-
-    # The overlap region: the end of fwd overlaps with the start of rev_rc.
-    # overlap_len ranges from min_ov to min(fwd_len, rev_len).
-    max_overlap = min(fwd_len, rev_len)
-
-    for ov in range(min_ov, max_overlap + 1):
-        fwd_region = fwd_seq[fwd_len - ov:]
-        rev_region = rev_rc[:ov]
-
-        mismatches = sum(1 for a, b in zip(fwd_region, rev_region) if a != b)
-        match_ratio = 1.0 - mismatches / ov
-
-        if match_ratio > 0.9 and (ov > best_overlap or
-                                   (ov == best_overlap and
-                                    mismatches < best_mismatches)):
-            best_overlap = ov
-            best_mismatches = mismatches
-            best_pos = ov
-
-    if best_pos < 0:
-        return None
-
-    # Build merged sequence: fwd non-overlap + overlap (use fwd) + rev non-overlap
-    ov = best_overlap
-    merged = fwd_seq + rev_rc[ov:]
-    return (ov, best_mismatches, merged)
-
-
-def merge_pairs(fwd_denoised, rev_denoised, min_ov=10):
-    """Merge denoised forward and reverse reads by overlap alignment.
-
-    Args:
-        fwd_denoised: dict {sequence: abundance} from dada() on forward reads
-        rev_denoised: dict {sequence: abundance} from dada() on reverse reads
-        min_ov: minimum overlap length
-
-    Returns:
-        dict {merged_sequence: abundance}
-    """
-    merged = {}
-
-    fwd_seqs = sorted(fwd_denoised.keys(),
-                       key=lambda s: fwd_denoised[s], reverse=True)
-    rev_seqs = sorted(rev_denoised.keys(),
-                       key=lambda s: rev_denoised[s], reverse=True)
-
-    for fwd_seq in fwd_seqs:
-        fwd_abund = fwd_denoised[fwd_seq]
-        for rev_seq in rev_seqs:
-            rev_abund = rev_denoised[rev_seq]
-
-            # Use the minimum abundance of the pair
-            pair_abund = min(fwd_abund, rev_abund)
-            if pair_abund <= 0:
-                continue
-
-            rev_rc = reverse_complement(rev_seq)
-            result = find_best_overlap(fwd_seq, rev_rc, min_ov)
-
-            if result is not None:
-                ov, mm, merged_seq = result
-                merged[merged_seq] = merged.get(merged_seq, 0) + pair_abund
-
-    return merged
-
-
-# ---------------------------------------------------------------------------
-# Per-plate chimera removal (bimera detection)
-#
-# A simple abundance-ratio approach: for each ASV, check if it can be
-# reconstructed as a chimera of two more-abundant parent ASVs. A sequence
-# is flagged as a bimera if its left half closely matches one parent and
-# its right half closely matches another.
-# ---------------------------------------------------------------------------
-def is_bimera(query, parents, min_fold=1.5, min_match_frac=0.9):
-    """Check if query is a bimera (chimera of two parents).
-
-    Splits the query at every possible breakpoint and checks if the left
-    half matches one parent and the right half matches another parent.
-    """
-    qlen = len(query)
-    if qlen < 20:
-        return False
-
-    for bp in range(10, qlen - 10):
-        left = query[:bp]
-        right = query[bp:]
-        left_len = len(left)
-        right_len = len(right)
-
-        for i, p1 in enumerate(parents):
-            if len(p1) < qlen:
-                continue
-            # Check left half against p1
-            p1_left = p1[:left_len]
-            left_matches = sum(1 for a, b in zip(left, p1_left) if a == b)
-            if left_matches / left_len < min_match_frac:
-                continue
-
-            for j, p2 in enumerate(parents):
-                if i == j:
-                    continue
-                if len(p2) < qlen:
-                    continue
-                # Check right half against p2 (aligned from the right end)
-                p2_right = p2[len(p2) - right_len:]
-                right_matches = sum(1 for a, b in zip(right, p2_right)
-                                    if a == b)
-                if right_matches / right_len >= min_match_frac:
-                    return True
-
-    return False
-
-
-def remove_bimeras_from_table(seqtab_dict, min_fold=1.5):
-    """Remove bimeric sequences from a sequence table.
-
-    Args:
-        seqtab_dict: dict {sequence: total_abundance}
-
-    Returns:
-        set of chimeric sequences
-    """
-    # Sort by abundance descending
-    sorted_seqs = sorted(seqtab_dict.keys(),
-                          key=lambda s: seqtab_dict[s], reverse=True)
-
-    chimeras = set()
-
-    for idx, query in enumerate(sorted_seqs):
-        query_abund = seqtab_dict[query]
-        # Parents must be more abundant (by min_fold)
-        parents = [s for s in sorted_seqs[:idx]
-                   if seqtab_dict[s] >= min_fold * query_abund]
-
-        if len(parents) < 2:
-            continue
-
-        # Limit to top 20 parents to keep runtime reasonable
-        parents = parents[:20]
-
-        if is_bimera(query, parents, min_fold=min_fold):
-            chimeras.add(query)
-
-    return chimeras
 
 
 # ---------------------------------------------------------------------------
@@ -271,18 +107,24 @@ def remove_bimeras_from_table(seqtab_dict, min_fold=1.5):
         print(f"[WARNING] {sample_name}: dada() failed: {e}")
         continue
 
-    # Extract denoised sequences (may be missing if sample is too small)
-    fwd_denoised = ddF.get("denoised", {})
-    rev_denoised = ddR.get("denoised", {})
-
-    if not fwd_denoised or not rev_denoised:
+    # Check for denoised output
+    if not ddF.get("cluster_seqs") or not ddR.get("cluster_seqs"):
         print(f"[WARNING] {sample_name}: no denoised sequences, skipping")
         continue
 
-    merged = merge_pairs(fwd_denoised, rev_denoised, min_ov=min_overlap)
+    # Merge pairs using C-level NW alignment (matches R's mergePairs)
+    merged_list = dada2py.merge_pairs(ddF, derepF, ddR, derepR,
+                                       min_overlap=min_overlap,
+                                       trim_overhang=True, verbose=False)
+
+    # Collect accepted merged reads into {sequence: abundance}
+    accepted = {}
+    for m in merged_list:
+        if m["accept"]:
+            accepted[m["sequence"]] = accepted.get(m["sequence"], 0) + m["abundance"]
 
-    if len(merged) > 0:
-        all_merged[sample_name] = merged
+    if accepted:
+        all_merged[sample_name] = accepted
 
 # Drop samples that produced no merged reads
 if len(all_merged) == 0:
@@ -311,15 +153,25 @@ def remove_bimeras_from_table(seqtab_dict, min_fold=1.5):
       f"{n_asvs} unique ASVs, {total_reads} total reads")
 
 # ---------------------------------------------------------------------------
-# Per-plate chimera removal
+# Per-plate chimera removal using dada2's C-level consensus method
 # ---------------------------------------------------------------------------
-seq_totals = dt.groupby("sequence")["count"].sum().to_dict()
+all_seqs = sorted(set(s for d in all_merged.values() for s in d))
+sample_list = list(all_merged.keys())
+seq_to_idx = {s: i for i, s in enumerate(all_seqs)}
+
+mat = np.zeros((len(sample_list), len(all_seqs)), dtype=np.int32)
+for si, sam in enumerate(sample_list):
+    for seq, ab in all_merged[sam].items():
+        mat[si, seq_to_idx[seq]] = ab
 
-chimeric_seqs = remove_bimeras_from_table(seq_totals)
+seqtab = {"table": mat, "seqs": all_seqs}
+chim_result = dada2py.remove_bimera_denovo(seqtab, verbose=True)
 
+is_chimera = np.array(chim_result["is_chimera"])
+chimeric_seqs = set(s for s, c in zip(all_seqs, is_chimera) if c)
 dt_nochim = dt[~dt["sequence"].isin(chimeric_seqs)].copy()
 
-n_chimeras = len(chimeric_seqs)
+n_chimeras = int(is_chimera.sum())
 n_asvs_after = dt_nochim["sequence"].nunique()
 reads_after = dt_nochim["count"].sum()
 pct = round(reads_after / max(total_reads, 1) * 100, 1)