Bug in ARG + taxonomy merging #491

Activity

1. jasmezz

   added

   bugSomething isn't workingSomething isn't working

   on Aug 14, 2025
2. jasmezz commented on Aug 14, 2025

   

   jasmezz

   on Aug 14, 2025

   CollaboratorAuthor

   More actions

   This might be easy to solve @Darcy220606 ? Looking at the merging script, line 269 seems responsible. Maybe just have to change str(x) in tool_contig_id into str(x) == tool_contig_id?
3. jfy133

   added

   requires-upstream-fix

   on Nov 7, 2025
4. cjfields commented on Sep 16, 2026

   

   cjfields

   on Sep 16, 2026

   More actions

   Hi, I wanted to confirm that we ran into a pretty extreme version of this bug.
   The one-line fix proposed in the thread, on it's own, unfortunately isn't
   sufficient since it misses another bug that results in mixing contigs with
   the same name from different samples.

   Full disclosure: I want to mention this was sorted out with some help from
   Claude Code, but the results were verified by me. We have a fix for it but
   the code is specific for AMR, and I think it's worth a discussion first
   before any changes are addressed.

   Sample background and funcscan results

   We have results from 263 assemblies in a recent run (r4.0.0), submitted across
   three batches, with one batch consisting of MEGAHIT assemblies(mixed cultures
   of enterics) and the other samples as Unicycler/spaDES assemblies
   (pure cultures). All sample names are generic: 'Sample1', 'Sample2', all the
   way to 'Sample263'.

   In the worst batch, all Unicycler assemblies, 17,437 real AMR hits came back
   as 6,327,544 rows (363x) with the taxonomy columns wrong. Nothing was
   lost — every hit is still present, checked against the per-sample
   hAMRonization reports — but the taxonomy isn't usable in all three batches.
   Assemblers that number contigs 1..N(Unicycler) are hit hardest though MEGAHIT's
   k141_* names still inflated ~11x. The only reason we noticed this is due to
   the huge increase in the number of records leading to a larger-than expected
   file size for the one worst-case batch mentioned above.

   Narrowing it down, in bin/merge_taxonomy.py::hamronization_taxa there are
   two different issues:

   samples_taxa = taxa_df['sample_id'].unique()
   for sampleID in samples_taxa:                                            # 252
       subset_tool = tool_df[tool_df['sample_id'].str.contains(sampleID)]   # 254
       subset_taxa = taxa_df[taxa_df['sample_id'].str.contains(sampleID)]   # 256
       ...
       matches = subset_taxa[
           subset_taxa['contig_id_taxa'].apply(lambda x: str(x) in tool_contig_id)  # 269
       ]

   1. In Line 269 — the contig test is a substring test, so contig 70 picks up
      taxonomy for contigs 7 and 70.
   2. Lines 254/256 — the sample scoping is also a substring test. The
      iteration for sampleID="Sample1" pulls in every sample whose name contains
      Sample1 (88 of them in our worst-case batch) and cross-joins within that
      block, and is what puts another sample's lineage on a hit.

   Running this manually reproduced our issue exactly: 88 for any Sample1* hit,
   7 for Sample21*, 1 for Sample200. The outer loop also re-emits each tool
   row once per matching sampleID, which is the remaining duplication.

   Testing the == fix mentioned above on a small set (samples Sample1,
   Sample13, Sample133, Sample200; contigs 1/7/70; one hit each on
   contig 70):

   variant	rows for 4 hits	rows for the Sample133 hit	distinct taxa on it
   current	30	12	6 (truth: 1)
   == at line 269	15	6	3 (truth: 1)
   correct	4	1	1

   This halves the row count, however this leaves cross-sample contamination:
   contig 70 exact-matches contig 70 in Sample1 and Sample13.

   So, is the logic above meant to deal with differences in the output
   extensions? If so I'd recommend it might be better to have a dedicated lookup
   table for each analysis with a mapping to the extension to strip, followed by
   exact string matches, but I wanted to make sure about this.

   I also notices that deeparg results report ORF locus tags, so it's not
   matching anything. Fixing deeparg means resolving locus tag -> contig
   mappings from the annotation before the merge. funcscan should have this
   GFF3 in hand at that point, and our in-house fix does add this as an option.

   I should add: it's possible AMP and BGC are also affected, they also seem to
   use str.contains, but I haven't confirmed, our use case was strictly AMR.
5. cjfields commented on Sep 16, 2026

   

   cjfields

   on Sep 16, 2026

   More actions

   I should add, a very similar issue, paleobiotechnology/AMPcombi#105, has been reported upstream in 2025 and which has not been addressed.