I recently ran NCBI's AMRFinderPlus to detect antimicrobial resistance on two different datasets of clinical isolates of Mycobacterium tuberculosis.
Since the input file requires is a fasta file, I had fastq input files which I then aligned to a reference genome using BWA-MEM, used bcftools to create vcf files (variant calling), compressed and indexed them to then run bcftools consensus to generate a consensus genome assembly fasta file for each single isolate to use as input for AMRFinder. However, the two different datasets had very similar results (which doesn't make much biological sense) and I was wondering if I should follow a different approach to generate the fasta input files from the fastq isolate files. I saw SPAdes is common in this aspect but I read it's best for de novo assembly, and since I'm interested in SNPs I have to do a reference-based assembly.
Here's one of the output tsv files from AMRFinder:
Protein id Contig id Start Stop Strand Element symbol Element name Scope Type Subtype Class Subclass Method Target length Reference sequence length % Coverage of reference % Identity to reference Alignment length Closest reference accession Closest reference name HMM accession HMM description
NA NC_000962.3 314347 314889 - aac(2')-Ic aminoglycoside N-acetyltransferase AAC(2')-Ic core AMR AMR AMINOGLYCOSIDE GENTAMICIN/TOBRAMYCIN EXACTX 181 181 100.00 100.00 181 WP_003899880.1 aminoglycoside N-acetyltransferase AAC(2')-Ic NA NA
NA NC_000962.3 2231620 2232156 + erm(37) 23S rRNA (adenine(2058)-N(6))-methyltransferase Erm(37) core AMR AMR LINCOSAMIDE/MACROLIDE AZITHROMYCIN/CLARITHROMYCIN/CLINDAMYCIN/ERYTHROMYCIN/TELITHROMYCIN EXACTX 179 179 100.00 100.00 179 WP_003900446.1 23S rRNA (adenine(2058)-N(6))-methyltransferase Erm(37) NA NA
NA NC_000962.3 2325827 2326747 - blaC class A beta-lactamase BlaC core AMR AMR BETA-LACTAM BETA-LACTAM EXACTX 307 307 100.00 100.00 307 WP_003410677.1 class A beta-lactamase BlaC NA NA
the only differences I saw among the two datasets were ONLY in the start and stop columns, the rest - subclass, closest reference accession, closest reference name... - were exactly the same. One of the datasets is multiresistant and the other is monoresistant, I've already verified with Mykrobe, TBProfiler and ResFinder.
Has anyone here used AMRFinder with clinical isolate fastq files before? If so, how did you obtain the fasta files to use as input?
2 answers
Those three hits are themselves the explanation. aac(2')-Ic, erm(37) and blaC are intrinsic to M. tuberculosis - they sit in H37Rv and in every clinical isolate - so getting the identical set out of both datasets is the correct result rather than a bug. Start/stop shifting slightly is just indels moving coordinates in your consensus.
The deeper issue is that TB resistance is almost entirely point mutations in chromosomal genes: rpoB, katG, inhA, gyrA, pncA. A default AMRFinderPlus run only screens for acquired resistance genes, so if you didn't pass --organism it never looked for any of those at all. Worth running --list_organisms on your version to see whether Mycobacterium tuberculosis is in there.
Honestly though, TBProfiler and Mykrobe are built around the WHO mutation catalogue and are the right tools here, and they already separated your two datasets. I'd take their answer as the real one rather than trying to make AMRFinderPlus reproduce it.
Your fasta-generation pipeline (BWA-MEM -> bcftools call -> bcftools consensus) is a reasonably standard approach and probably isn't your core issue. I think the real problem is upstream of that, in how AMRFinderPlus is being run.
aac(2')-Ic, erm(37), and blaC are intrinsic, species-wide genes present in essentially every M. tuberculosis genome — including the H37Rv reference itself (which is exactly what NC_000962.3 is). AMRFinderPlus's default mode screens for acquired/core AMR genes via BLAST-based presence/absence, which is the right model for organisms like Enterobacteriaceae where resistance genes are horizontally acquired and vary between strains. But it's the wrong lens for MTB: essentially all clinically relevant MTB resistance comes from point mutations in specific resistance-determining regions (rpoB S450L-type mutations for rifampicin, katG/inhA for isoniazid, gyrA/gyrB for fluoroquinolones, embB, pncA, etc.), not from presence/absence of acquired genes. Since those three genes are always present and unchanged in essentially every MTB isolate, AMRFinderPlus's core screen will report near-identical hits regardless of the isolate's actual resistance profile — which is exactly what you're seeing.
Mykrobe, TBProfiler, and ResFinder all differ from this by querying curated MTB-specific point-mutation catalogs (e.g., the WHO mutation catalogue) rather than doing generic gene presence/absence — that's why they correctly separate your multi- vs. mono-resistant datasets and AMRFinderPlus doesn't.
The fix: rerun AMRFinderPlus with the --organism Mycobacterium_tuberculosis flag. This activates AMRFinderPlus's organism-specific point-mutation detection module, which screens the actual resistance-determining positions rather than relying solely on the generic acquired-gene database. If you ran it without that flag, you were only ever going to see the intrinsic core genes, not the mutation-driven differences between your isolates.
Separately, worth double-checking coverage depth at the key resistance loci in your bcftools consensus step — low-coverage positions can silently default to the reference allele rather than the true isolate variant unless you're explicitly masking low-coverage regions, which could mask real resistance SNPs even with the --organism flag enabled. But I'd fix the --organism flag first and see if that alone resolves the discordance before troubleshooting the assembly step further.
Log in to answer this question.
Please read through this thread --> Generating consensus sequence from bam file
You may want to do
samtools consensusinstead of the procedure above and see if that helps. Michael's answer explains what happens when one follows the procedure you used above.