Issue with reporting resistance OXA-48 #149

talnor · 2021-11-26T13:18:33Z

Describe the bug
Perfect matches to OXA-48 are no longer reporter in microSALT, instead microSALT reports OXA-566 resistance of lower %identities.

Message from KS in Ticket 722237:

...angående att OXA-48 har börjat rapporteras som OXA-566 (identity 99,88%). Kör jag fastq-filerna från er i Resfinder så rapporterar den OXA-48 (100% match). Jag skulle då meddela ungefär när denna förändring inträffat och vad jag kan se så rapporteras prov 20ET200231 (Ticket 112544, v44 2020) som OXA-48 medan prov 20ET500236 (Ticket 935485, v47 2020) rapporteras som OXA-566. Är det något som ändrats i assembly under denna period?

To Reproduce
Steps to reproduce the behavior:

Run sample 20ET200231 or 20ET500236 with microSALT > 3.1.3 and microSALT 3.1.0

Expected behavior
OXA-48 should be reported if this resistances has higher %identity.

Additional context
Is OXA-48 in the raw blast results?

Adding some of the commits made during the period below. Has any affected OXA-48 reporting? Scrape_blast changes likely.

Scrape_blast: Attempt to make resistance searches more compatible
Scrape_blast: Locilengths now ignores reference name
Scrape_blast: Rolledback expec locilength generous lookup
Scrape_blast: Added padder safety
Scrape_blast: from Clinical-Genomics/resistance-scrape
Referencer: Updated and fixed bugs in pubmlst and ncbi downloads
Referencer: Clarified logic in referencer
Referencer: from Clinical-Genomics/fix_referencer

talnor · 2021-11-26T13:59:16Z

Currently blocked by installation issues

moahaegglund · 2022-10-19T10:17:34Z

No longer blocked, installation issues solved in #142.

pbiology · 2023-06-29T07:53:14Z

Will be resolved in Clinical-Genomics/CG-JASEN#5

samuell · 2024-06-19T14:47:18Z

Re-opening this, as it continues to cause problems, and we might need to prioritize it if Jasen is not completed soon.

samuell · 2024-06-19T14:58:17Z

I have been running Spades (4.0.0) assembly + Blasting against the resfinder DB semi-manually, using the same commands as used in microSALT, both with and without the --isolate flag to Spades, and the OXA-48 gene is reported correctly with both methods.

This thus strengthens the suspicion that the error might be in the scraping logic.

(Unless there are some differences between Spades 3.13.1 used in microSALT right now, and 4.0.0.)

samuell · 2024-06-25T11:56:47Z

(Unless there are some differences between Spades 3.13.1 used in microSALT right now, and 4.0.0.)

Tested this also with Spades 3.13.1 with the same result, so that is apparently not the issue.

samuell · 2024-06-25T14:27:59Z

I have now confirmed that an existing OXA-48 is removed in this part of the code (for cleaning up overlapping hits):

microSALT/microSALT/utils/scraper.py

Lines 344 to 410 in 2f8b19e

    
           # Cleanup of overlapping hits 
        
           if type == "seq_type": 
        
               identifier = "loci" 
        
           elif type == "resistance" or type == "expec": 
        
               identifier = "gene" 
        
           ind = 0 
        
           while ind < len(hypo) - 1: 
        
               targ = ind + 1 
        
               while targ < len(hypo): 
        
                   ignore = False 
        
                   if ( 
        
                       hypo[ind]["contig_name"] == hypo[targ]["contig_name"] 
        
                       or hypo[ind][identifier] == hypo[targ][identifier] 
        
                   ): 
        
                       # Overlapping or shared gene 
        
                       if ( 
        
                           ( 
        
                               hypo[ind].get("contig_start") 
        
                               >= hypo[targ].get("contig_start") 
        
                               and hypo[ind].get("contig_start") 
        
                               <= hypo[targ].get("contig_end") 
        
                           ) 
        
                           or ( 
        
                               hypo[ind].get("contig_end") 
        
                               >= hypo[targ].get("contig_start") 
        
                               and hypo[ind].get("contig_end") 
        
                               <= hypo[targ].get("contig_end") 
        
                           ) 
        
                           or (hypo[ind].get(identifier) == hypo[targ].get(identifier)) 
        
                       ): 
        
                           # Rightmost is worse 
        
                           if float(hypo[ind].get("identity")) * ( 
        
                               1 - abs(1 - hypo[ind].get("span")) 
        
                           ) > float(hypo[targ].get("identity")) * ( 
        
                               1 - abs(1 - hypo[targ].get("span")) 
        
                           ): 
        
                               del hypo[targ] 
        
                               ignore = True 
        
                           # Leftmost is worse 
        
                           elif float(hypo[ind].get("identity")) * ( 
        
                               1 - abs(1 - hypo[ind].get("span")) 
        
                           ) < float(hypo[targ].get("identity")) * ( 
        
                               1 - abs(1 - hypo[targ].get("span")) 
        
                           ): 
        
                               del hypo[ind] 
        
                               targ = ind + 1 
        
                               ignore = True 
        
                           # Identical identity and span, seperating based on contig coverage 
        
                           else: 
        
                               # Rightmost is worse 
        
                               if float(hypo[ind].get("contig_coverage")) >= float( 
        
                                   hypo[targ].get("contig_coverage") 
        
                               ): 
        
                                   del hypo[targ] 
        
                                   ignore = True 
        
                               # Leftmost is worse 
        
                               elif float(hypo[ind].get("contig_coverage")) < float( 
        
                                   hypo[targ].get("contig_coverage") 
        
                               ): 
        
                                   del hypo[ind] 
        
                                   targ = ind + 1 
        
                                   ignore = True 
        
                   if not ignore: 
        
                       targ += 1 
        
                   else: 
        
                       pass 
        
               ind += 1

Perhaps there is an overlap with OXA-566 ... we'll see.

samuell · 2024-06-25T14:53:20Z

And in one of these commits:

deed81b Removed newline char in log
ab280db Update pull_request_template.md
7e38873 Clarified logic in referencer
9757f6b Updated and fixed bugs in pubmlst and ncbi downloads
27e4685 Merge pull request #121 from Clinical-Genomics/resistance-scrape
170b013 Version bump
b4e1231 Added padder safety
5749721 Quast backwards fix
2e626ef Rolledback expec locilength generous lookup
395a218 Locilengths now ignores reference name
a05ef89 Update init.py
3f1488c Attempt to make resistance searches more compatible

samuell · 2024-06-26T07:18:06Z

It is introduced with this one:

395a218 Locilengths now ignores reference name

(Tested by cherry-picking the commits one by one into a test branch, using a real-world resistance blast output containing both OXA-566 and OXA-48 instead of the test one, and looking for occurrences of OXA-48 in caplog.text).

This needs to be follow up with additional checks The same change might need to be done for other matching as well.

… prefix

samuell · 2024-07-11T20:43:38Z

I think I'm starting to have a reasonable test now, that now fails (before implementing the fix), saying:

>     assert "blaOXA-48" in genes
E     assert 'blaOXA-48' in ["aph(3')-III", 'ant(6)-Ia', 'blaOXA-566']

(See https://github.com/Clinical-Genomics/microSALT/actions/runs/9898158769/job/27344305202 )

This is triggered by using some real-world example data where both OXA-48 and OXA-566 are included. But re-assuring to see that we reproduce the inclusion of blaOXA-566 there, as the issue described.

And, with the fix implemented, this test passes.

talnor added the bug label Nov 26, 2021

talnor added the KS requests label Nov 26, 2021

pbiology closed this as completed Jun 29, 2023

samuell reopened this Jun 19, 2024

samuell added a commit that referenced this issue Jun 26, 2024

Fix #149: Avoid promiscuous matchin of gene name

81fb2ff

This needs to be follow up with additional checks The same change might need to be done for other matching as well.

samuell added a commit that referenced this issue Jun 26, 2024

Fix #149: Avoid promiscuous matching of gene name

ca8a56f

This needs to be follow up with additional checks The same change might need to be done for other matching as well.

samuell linked a pull request Jun 26, 2024 that will close this issue

Fix #149: Missing or mis-reported resistances #173

Draft

9 tasks

samuell added a commit that referenced this issue Jul 5, 2024

Fix #149: Avoid promiscuous matching of gene name

264cb2a

This needs to be follow up with additional checks The same change might need to be done for other matching as well.

samuell added a commit that referenced this issue Jul 11, 2024

Mark where #149 can be fixed

e00b092

This needs to be follow up with additional checks The same change might need to be done for other matching as well.

samuell added a commit that referenced this issue Jul 11, 2024

Fix #149: Make regex more specific to not match other genes with same…

c9ee900

… prefix

samuell added a commit that referenced this issue Jul 11, 2024

Fix #149: Make regex more specific to not match other genes with same…

036d5c9

… prefix

samuell added a commit that referenced this issue Jul 11, 2024

Fix #149 by making regex specific and not matching multiple genes

0610b60

samuell added a commit that referenced this issue Jul 11, 2024

Fix #149 by making regex specific and not matching multiple genes

c32b23f

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Issue with reporting resistance OXA-48 #149

Issue with reporting resistance OXA-48 #149

talnor commented Nov 26, 2021 •

edited

Loading

talnor commented Nov 26, 2021

moahaegglund commented Oct 19, 2022

pbiology commented Jun 29, 2023

samuell commented Jun 19, 2024 •

edited

Loading

samuell commented Jun 19, 2024

samuell commented Jun 25, 2024

samuell commented Jun 25, 2024 •

edited

Loading

samuell commented Jun 25, 2024

samuell commented Jun 26, 2024

samuell commented Jul 11, 2024 •

edited

Loading

Issue with reporting resistance OXA-48 #149

Issue with reporting resistance OXA-48 #149

Comments

talnor commented Nov 26, 2021 • edited Loading

talnor commented Nov 26, 2021

moahaegglund commented Oct 19, 2022

pbiology commented Jun 29, 2023

samuell commented Jun 19, 2024 • edited Loading

samuell commented Jun 19, 2024

samuell commented Jun 25, 2024

samuell commented Jun 25, 2024 • edited Loading

samuell commented Jun 25, 2024

samuell commented Jun 26, 2024

samuell commented Jul 11, 2024 • edited Loading

talnor commented Nov 26, 2021 •

edited

Loading

samuell commented Jun 19, 2024 •

edited

Loading

samuell commented Jun 25, 2024 •

edited

Loading

samuell commented Jul 11, 2024 •

edited

Loading