This tutorial covers the bam_coverage_evaluation pipeline, which takes a set of pre-aligned BAM files and evaluates the resistome across a grid of alignment filter thresholds. Instead of committing to a single cutoff, the evaluation shows you how the number of detected genes, groups, mechanisms and classes drops off as you tighten each threshold, so you can choose sensible cutoffs and see how sensitive each sample is to them. It produces combined result matrices and a set of drop-off plots.
Before you start: This pipeline runs on BAM files you have already produced with AMR++ (for example, the
*_alignment_sorted.bamfiles underAlignment/BAM_files/Standard/). It also requires a few extra tools that are not bundled in the default AMR++ conda environment. See Additional tools required below.
- Additional tools required
- What the evaluation does
- Running the evaluation
- The three filters
- Choosing what to evaluation
- Outputs
- Interpreting the plots
- Applying the selected filter thresholds
The evaluation uses a Python analysis step and an R plotting step whose dependencies are not included in the AMR++ conda environment, because adding them, especially the R plotting stack, would substantially increase an already large environment. Install them separately before running the evaluation.
Python (for coverage_threshold_evaluation.py and combine_evaluation_results.py):
pysamfor reading BAM alignmentsduckdbfor fast aggregation of the per-alignment table
# into your active AMR++ environment, or any python3 env you point the pipeline at
pip install pysam duckdbR (for plot_evaluation_dropoff.R):
data.tableggplot2scalespatchworkfor the two-panel heat maps
# from an R session, or Rscript -e '...'
install.packages(c("data.table", "ggplot2", "scales", "patchwork"))If any of these are missing, the step that needs them will fail. The Python step is required to produce any output. The R step only affects the figures, so without R you still get the combined CSV matrices. If
patchworkalone is missing, the plots are still written, but the count and percentage heat maps come out as separate files rather than paired panels.
Make sure Rscript is on your PATH. The pipeline calls it through the $RSCRIPT environment variable, so set RSCRIPT in your config or profile if your Rscript lives somewhere non-standard.
For each BAM, the pipeline:
- Streams every primary alignment once and records, per alignment, how much of the read aligned, how well it matched, and which reference positions it covered.
- Re-counts the resistome at every combination in the threshold grid, using the same fragment-level, group-aware, edge-aware counting as the main pipeline.
- Writes per-sample result CSVs.
Then across all samples:
- Combines the per-sample CSVs into single matrices with
combine_evaluation_results.py. - Produces drop-off plots and summary tables with
plot_evaluation_dropoff.R.
Because the BAM is parsed only once and the whole grid is evaluated from the resulting intermediate table, the full grid costs roughly the same as a single counting run.
A sample that ends up with no usable alignments, for example one where every alignment was to a gene requiring SNP confirmation and evaluation_exclude_snp removed it, is handled gracefully. The evaluation writes header-only outputs for that sample, prints a warning naming the filters responsible, and the run continues.
nextflow run main_AMR++.nf \
-profile local \
--pipeline bam_coverage_evaluation \
--bam_files "/path/to/Alignment/BAM_files/Standard/*_alignment_sorted.bam" \
--output Filtering_sensitivity_analysis--bam_filesis a glob matching the BAMs to evaluation. Quote it so the shell does not expand it.--outputis where results are written.- On a cluster, swap
-profile localfor your SLURM profile. See Running with SLURM.
Example with the SLURM profile:
nextflow run main_AMR++.nf \
--pipeline bam_coverage_evaluation \
--bam_files "AMR++_output/Alignment/BAM_files/Standard/*" \
--output Filtering_sensitivity_analysis \
-profile local_slurmAMR++ applies three independent alignment filters. The evaluation explores the same three, so a threshold you pick from a evaluation can be applied directly to a counting run.
All values are proportions from 0 to 1. A value of 0 turns a filter off.
| Filter | What it asks |
|---|---|
| gene fraction | What proportion of the gene was covered by retained alignments? |
| query coverage | What proportion of the read aligned? |
| identity | What proportion of the aligned portion matched the reference? |
These are independent rather than substitutes for one another. A short fragment that aligns perfectly scores high on identity and low on query coverage. A full-length read carrying many mismatches scores the reverse.
Query coverage has two definitions, and you choose between them rather than applying both. The threshold value means the same thing either way. Only the metric it is compared against changes.
| Parameter | Calculation | Behavior |
|---|---|---|
evaluation_query_coverage |
aligned_length / read_length |
Counts mismatches inside the aligned region as covered |
evaluation_match_qcov |
(aligned_length - NM) / read_length |
Counts only genuine matches |
Because the edit distance is never negative, the match-based value can never exceed the standard one, so a given threshold is always at least as strict under evaluation_match_qcov.
In the counting step the same choice is made with the match_qcov parameter, which is "Y" or "N" rather than a list of values.
evaluation_match_qcov and evaluation_identity are both computed from the NM tag, which records each alignment's edit distance. NM is an optional SAM field, so it is not always present.
With the filter off, alignments lacking NM still pass. At any threshold above 0 they are excluded, because their identity cannot be verified. BWA-MEM emits NM by default, so BAMs produced by the standard AMR++ alignment step are fine, but BAMs rewritten by other tools may not carry it.
If you need to add it:
samtools calmd -b input.bam reference.fasta > output.nm.bam
samtools index output.nm.bamThe evaluation always tests a two-dimensional grid: gene fraction against one read-level filter.
Axis 1 is always gene fraction. evaluation_gene_fraction is evaluated in every run and appears in every figure.
Axis 2 is one of the three read-level filters. Set your choice to a list of values and leave the others at 0.
The relevant params.config block:
/* AXIS 1: always evaluated */
evaluation_gene_fraction = "0,0.1,0.25,0.5,0.8"
/* AXIS 2: evaluation ONE of these three; leave the others at "0",
* or give one a single value to apply it as a fixed filter */
evaluation_query_coverage = "0,0.5,0.6,0.7,0.8,0.9,0.95"
evaluation_match_qcov = "0"
evaluation_identity = "0"
/* Applied at every point in the grid regardless of what is evaluated */
evaluation_edge_aware_qcov = "Y"
evaluation_exclude_snp = "N"evaluation_query_coverage and evaluation_match_qcov measure the same property, so evaluating both is meaningless. Pick whichever definition of query coverage you want to filter on. evaluation_identity measures something different and is a genuine alternative axis.
A single non-zero value is a filter applied at every point in the grid rather than a evaluated axis. This is how you hold one filter constant while exploring another:
// evaluation gene fraction against query coverage, with identity pinned at 0.9
evaluation_gene_fraction = "0,0.1,0.25,0.5,0.8"
evaluation_query_coverage = "0,0.5,0.6,0.7,0.8,0.9,0.95"
evaluation_identity = "0.9"The figures plot gene fraction against query coverage, and the plot subtitle records that identity was fixed at 0.9 so the constraint stays visible.
nextflow run main_AMR++.nf --pipeline bam_coverage_evaluation \
--bam_files "AMR++_output/Alignment/BAM_files/Standard/*" \
--evaluation_gene_fraction '0,0.5,0.9' \
--evaluation_query_coverage '0,0.6,0.7,0.8,0.85,0.9,0.95' \
--output Filtering_sensitivity_analysisQuote the comma-separated values. Most shells handle them unquoted, but quoting removes any ambiguity, particularly inside SLURM submission scripts.
The plotting script reads the gene-fraction grid from the results themselves, so changing the grid needs no corresponding edit anywhere else.
Every combination is evaluated, so the grid grows multiplicatively. Rows written to combined_results.csv per sample:
| Configuration | Grid | Rows per sample |
|---|---|---|
| Gene fraction against query coverage (default) | 5 x 7 | 35 |
| Gene fraction against identity | 5 x 5 | 25 |
| A finer grid, 8 gene fraction by 10 query coverage | 8 x 10 | 80 |
combined_gene_detail.csv multiplies these by the number of detected genes, so it grows quickly. Runtime scales far better than file size, since the BAM is parsed once, but keep an eye on disk usage for large sample sets.
Under --output, in CoverageEvaluation/:
PerSample/contains one set of CSVs per BAM:<prefix>_results.csv,_gene_detail.csv,_redundancy.csvand_length_quantiles.csv.Combined/contains the merged matrices across all samples:combined_results.csv,combined_gene_detail.csvandcombined_length_quantiles.csv, each with asample_idcolumn.Summary/contains the drop-off plots as PNG files and the summary tables as CSV.
Drop-off curves (dropoff_<level>_count_vs_primary.png and _pct_vs_primary.png) show, for each annotation level, how the count falls as the read-level threshold tightens. The bold line is the mean across samples and each thin grey line is one sample. One line is drawn per gene-fraction value, so both axes appear together.
Heat maps (dropoff_<level>_heatmap.png and reads_retained_combined_heatmap.png) show the full two-dimensional grid as a pair of panels: the absolute feature count on the left and the percentage of baseline retained on the right. The percentage panel uses a fixed 0 to 100 colour scale, so colour carries the same meaning across every level and every figure.
Reads retained and lost figures show what fraction of fragment hits survive each threshold relative to the no-filter baseline.
Expect the annotation levels to behave differently. Gene accession counts usually fall fastest, because accession-level detection is the most sensitive to database redundancy, while group, mechanism and class counts hold up longer. A threshold that removes many accessions but few groups is generally preferable to the reverse, since aggregation to the group level or above is what we recommend for statistical comparison in any case.
The evaluation shows what a threshold costs on your data. It cannot tell you which threshold is correct, because that depends on the question you are asking.
Once you have chosen thresholds, re-run counting on the same BAM files with bam_resistome. The parameter names match the evaluation parameters without the evaluation_ prefix, and use the same 0 to 1 proportion scale.
| Evaluation parameter | Counting parameter |
|---|---|
evaluation_gene_fraction |
--min_gene_fraction |
evaluation_query_coverage |
--min_query_coverage |
evaluation_identity |
--min_identity |
evaluation_match_qcov |
--match_qcov Y together with --min_query_coverage |
A run using a gene fraction of 0.5 and a query coverage of 0.8:
nextflow run main_AMR++.nf \
--pipeline bam_resistome \
--bam_files "AMR++_output/Alignment/BAM_files/Standard/*bam" \
--min_gene_fraction 0.5 \
--min_query_coverage 0.8 \
--output AMR++_gf50_qcov80 \
--prefix AMR_gf50_qcov80 \
-profile local_slurmIf your evaluation used evaluation_match_qcov rather than evaluation_query_coverage, set match_qcov when applying so the counting step computes query coverage the same way the evaluation did:
nextflow run main_AMR++.nf \
--pipeline bam_resistome \
--bam_files "AMR++_output/Alignment/BAM_files/Standard/*bam" \
--min_gene_fraction 0.5 \
--min_query_coverage 0.8 \
--match_qcov Y \
--output AMR++_gf50_matchqcov80 \
--prefix AMR_gf50_matchqcov80 \
-profile local_slurmThis produces a new count matrix under Results/ reflecting the selected filters, ready for aggregation and statistical analysis.
The count matrix filename does not encode the filter settings. Every run writes
<prefix>_analytic_matrix.csv, which by default isAMR_analytic_matrix.csvno matter which thresholds were applied. Two runs with different filters produce two files with identical names, and if they share an output directory the second silently overwrites the first. Nothing inside the file records how it was generated.Three habits prevent this from becoming a problem.
Give each run its own
--outputdirectory. The simplest safeguard, and it also keeps the per-sample intermediate files separated.Encode the thresholds in
--prefix. The prefix becomes part of the matrix filename, so--prefix AMR_gf50_qcov80yieldsAMR_gf50_qcov80_analytic_matrix.csv, which stays self-describing even after the file is moved or shared. Both examples above do this.Save the command. Keep the exact invocation alongside the results, or move the settings into a params file you can archive:
nextflow run main_AMR++.nf -params-file my_thresholds.yaml --pipeline bam_resistome ...